my_star/StRoot/StFemtoDstMaker/StO97FemtoDstMaker.cxx

#include "StO97FemtoDstMaker.h"
#include "TBranch.h"


ClassImp(StO97FemtoDstMaker)

//
// Magnetic field -- values obtainted from StFemtoDstPythia6Maker
//
const Float_t BFIELD = 4.510600000;
const Float_t RBFIELD = -4.510600000;


//
// Set maximum file size to 1.9 GB (Root has a 2GB limit)
//
#define MAXFILESIZE 1900000000


//_________________
StO97FemtoDstMaker::StO97FemtoDstMaker(const Char_t *iFileName,
                                       const Char_t *oFileName)
{
    std::cout << "[StO97FemtoDstMaker] Creating an instance... ";
    mOutFileName = oFileName;
    mInFileName  = iFileName;
    mTree = 0;
    mOutFile = 0;
    mInFile = 0;
    mCompression = 9;
    mNEvents = 0;
    mStop = false;
    mDb       = TDatabasePDG::Instance();
    mRandy    = new TRandom3(0);

    //
    // RefMult counters
    //
    mRefMult     = 0;
    mRefMult2    = 0;
    mRefMultPos  = 0;
    mRefMult2Pos = 0;

    //
    // Setup sphericity matrix
    //
    mMatrix = new TMatrixTSym<double>(2);
    mPtCut = 0.1;
    //
    // Setup default cut's
    //
    mTrackMomentum[0] = 0.05;
    mTrackMomentum[1] = 2.1;
    mTrackDca[0] = 0.;
    mTrackDca[1] = 3.;
    mTrackEta[0] = -1.1;
    mTrackEta[1] = 1.1;

    std::cout << "[OK]" << std::endl;
}



StO97FemtoDstMaker::~StO97FemtoDstMaker()
{
    delete mMatrix;
}



Int_t StO97FemtoDstMaker::Init()
{
    std::cout << "[StO97FemtoDstMaker] ---=== Initializing maker ===---" << std::endl;
    mEvent = new StFemtoEvent();
    //
    // Input file
    //
    std::cout << "[StO97FemtoDstMaker] Openning an input file... ";
    if (mInFile)
        fclose(mInFile);
    mInFile = fopen(mInFileName, "r");
    if (mInFile)
        std::cout << "[OK]" << std::endl;
    else
    {
        std::cout << "[FAIL]" << std::endl;
        return kStFatal;
    }
    //
    // Output file
    //
    std::cout << "[StO97FemtoDstMaker] Creating output file... ";
    if (mOutFile)
    {
        mOutFile->Close();
        mOutFile->Open(mOutFileName, "RECREATE");
    }
    else
        mOutFile = new TFile(mOutFileName, "RECREATE");
    if (mOutFile)
        std::cout << "[OK]" << std::endl;
    else
    {
        std::cout << "[FAIL]" << std::endl;
        return kStFatal;
    }
    SetCompressionlevel(mCompression); // default value from constructor
    //
    // TTree setup
    //
    std::cout << "[StO97FemtoDstMaker] Creating TTree... ";
    if (mTree)
        delete mTree;
    mTree = new TTree("StFemtoDst", "StFemtoDst");
    if (!mTree) {
        std::cout << "[FAIL]" << std::endl;
        return kStFatal;
    }
    mTree->SetMaxTreeSize(MAXFILESIZE);
    mTree->Branch("StFemtoEvent", "StFemtoEvent", &mEvent);
    mNBytes = 0;
    std::cout << "[OK]" << std::endl;
    //
    // Finalization of initialization :)
    //
    std::cout << "[StO97FemtoDstMaker] ---=== Initialization complite ===---" << std::endl;

    return StMaker::Init();
}



void StO97FemtoDstMaker::SetCompressionlevel(Int_t comp)
{
    mCompression = comp;
    std::cout << "[StO97FemtoDstMaker] Compression level = " << mCompression << std::endl;
    mOutFile->SetCompressionLevel(mCompression);
}



Int_t StO97FemtoDstMaker::Make()
{
    //
    // Read event
    //
    mNEvents++;
    mStop = ReadEvent();
    if (mStop)
        return kStOk;
    //
    // Read tracks
    //
    mStop = ReadTracks();
    if (!mStop)
        mNBytes += mTree->Fill();

    return kStOk;
}



Int_t StO97FemtoDstMaker::Finish()
{
    std::cout << "[StO97FemtoDstMaker] ---=== Finalization of maker ===---"
        << "[StO97FemtoDstMaker] total number of events = "
        << mNEvents << std::endl
        << "[StO97FemtoDstMaker] total bytes = "
        << mNBytes
        << std::endl;
    if (mInFile)
    {
        fclose(mInFile);
        mInFile = 0;
    }
    if (mTree)
    {
        if (mOutFile)
            mOutFile->Write();
        delete mTree;
        mTree = 0;
    }
    if (mOutFile)
    {
        mOutFile->Close();
        delete mOutFile;
        mOutFile = 0;
    }

    return kStOk;
}



Bool_t StO97FemtoDstMaker::ReadEvent()
{
    Int_t eventNumber;
    Float_t impPar;
    Float_t rot;
    Char_t buf[128];

    while (1)
    {
        Int_t ret = fscanf(mInFile, " %i %i %f %f",
                           &eventNumber, &mNTracks, &impPar, &rot);
        if (ret <= 0)               // if can't read a line
            fgets(buf, 127, mInFile); // skip it
        else if (ret == 4)
        { // 4 - total number of event parameters
            mEvent->Clear();
            mEvent->SetEventId(eventNumber);
            mEvent->SetRunId(0xDEADBEEF);
            mEvent->SetCollisionType(false);
            mEvent->SetCent16(0);
            mEvent->SetMagneticField(BFIELD);
            mEvent->SetPrimaryVertexPosition(impPar, rot, 0.); // why not?
            mEvent->SetVpdVz(0.);
            mEvent->SetTriggerId(0);
            mEvent->SetPrimaryVertexRanking(1e4);
            return false;
        }
        if (feof(mInFile)) // if EOF then return stop flag
            return true;
    }
}



Bool_t StO97FemtoDstMaker::ReadTracks()
{
    Int_t id;
    Int_t pdgId;
    Float_t px, py, pz;
    Float_t energy, mass;
    Float_t xFr, yFr, zFr, tFr;

    mSphericity    = 0;
    float pt_full  = 0;
    mRefMult2Pos   = 0;
    mRefMultPos    = 0;
    mRefMult2      = 0;
    mRefMult       = 0;
    int goodTracks = 0;

    mMatrix->Zero();

    for (Int_t iTrack = 0; iTrack < mNTracks; iTrack++)
    {
        Int_t ret = fscanf(mInFile, " %i %i %f %f %f %f %f %f %f %f %f",
                           &id, &pdgId, &px, &py, &pz,
                           &energy, &mass, &xFr, &yFr, &zFr, &tFr);
        if (ret != 11)
            return true;
        else
        {
            TParticlePDG *partPdg = mDb->GetParticle(pdgId);
            if (partPdg == 0x0)   continue;
            int   charge  = TDatabasePDG::Instance()->GetParticle(pdgId)->Charge();
            float ptot    = sqrt(px*px + py*py + pz*pz);
            float pt      = sqrt(px*px + py*py);
            float eta     = ptot > 0. && pz > 0. ? 0.5*log((ptot + pz)/(ptot - pz)) : -999.;
            float dca     = sqrt(xFr*xFr + yFr*yFr)*1.e-13; // DCA of track to 0 in [cm]
            float beta    = ptot/energy;

            if (charge == 0)   continue;
            charge /= abs(charge);

            if (fabs(eta) < 1.0 && dca < 3.0)
            {
                mRefMult2++;
                if (charge > 0)   mRefMult2Pos++;
                if (fabs(eta) < 0.5)
                {
                    mRefMult++;
                    if (charge > 0)   mRefMultPos++;
                }
            }

            //
            // Calculate sphericity
            //
            if (pt > mPtCut)
            {
                if (fabs(eta) < 1.0)
                {
                    (*mMatrix)(0, 0) += px*px/pt;
                    (*mMatrix)(1, 1) += py*py/pt;
                    (*mMatrix)(0, 1) += px*py/pt;
                    (*mMatrix)(1, 0) += py*px/pt;

                    pt_full += pt;
                }
            }

            float nSigmaElectron = -999.;
            float nSigmaPion     = -999.;
            float nSigmaKaon     = -999.;
            float nSigmaProton   = -999.;

            switch (abs(pdgId))
            {
                case 211:
                    nSigmaPion = 0.0;
                    break;
                case 321:
                    nSigmaKaon = 0.0;
                    break;
                case 2212:
                    nSigmaProton = 0.0;
                    break;
                default:
                    continue;
            }

            if  (ptot > mTrackMomentum[0] && ptot < mTrackMomentum[1] &&
                 eta > mTrackEta[0] && eta < mTrackEta[1] &&
                 charge != 0 &&
                 dca > mTrackDca[0] && dca < mTrackDca[1])
            {
                unsigned int tmap0 = 0; // This is
                unsigned int tmap1 = 0; //  topology map
                float dedx = dedxMean(mass, ptot);
                TopologyMap(&tmap0, &tmap1, ptot);

                StFemtoTrack *track = mEvent->AddFemtoTrack();
                track->SetId(iTrack);
                track->SetNHits(45*charge);
                track->SetNHitsPoss(45);
                track->SetNSigmaElectron(nSigmaElectron);
                track->SetNSigmaPion(nSigmaPion);
                track->SetNSigmaKaon(nSigmaKaon);
                track->SetNSigmaProton(nSigmaProton);
                track->SetDedx(dedx);
                track->SetMap0(tmap0);
                track->SetMap1(tmap1);
                track->SetP(px, py, pz);
                track->SetDCAGlobal(xFr*1.e-13, yFr*1.e-13, zFr*1.e-13);
                track->SetGlobalP(px, py, pz);
                track->SetBeta(beta);
                goodTracks++;
            }
        }
    }

    //
    // Calculate eigenvalues
    //
    *mMatrix *= 1./pt_full;
    TMatrixDSymEigen death_machine(*mMatrix);
    TVectorD eigen = death_machine.GetEigenValues();
    float mSphericity = 2.*eigen.Min()/(eigen[0] + eigen[1]);

    mEvent->SetSphericity(mSphericity);
    mEvent->SetRefMult(mRefMult);
    mEvent->SetRefMultCorr(mRefMult);
    mEvent->SetRefMultCorrWeight(1.0);
    mEvent->SetRefMultPos(mRefMultPos);
    mEvent->SetNumberOfBTofHit(mRefMult);
    mEvent->SetNumberOfPrimaryTracks(goodTracks);
    mEvent->SetNumberOfGlobalTracks(goodTracks);
    mEvent->SetRefMult2(mRefMult2);
    mEvent->SetRefMult2Pos(mRefMult2Pos);

    return false;
}



Double_t StO97FemtoDstMaker::dedxMean(Double_t mass, Double_t momentum)
{
    Double_t dedxMean;
    Double_t tpcDedxGain = 0.174325e-06;
    Double_t tpcDedxOffset = -2.71889; 
    Double_t tpcDedxRise = 776.626;

    Double_t gamma = TMath::Sqrt(momentum*momentum/(mass*mass)+1.);
    Double_t beta = TMath::Sqrt(1. - 1./(gamma*gamma));
    Double_t rise = tpcDedxRise* beta*beta*gamma*gamma;      
    if ( beta > 0)
        dedxMean = tpcDedxGain/(beta*beta) * (0.5*TMath::Log(rise) - beta*beta 
                                              - tpcDedxOffset);
    else
        dedxMean = 1000.;
    return dedxMean;
}



void StO97FemtoDstMaker::TopologyMap(unsigned int *tmap1, unsigned int *tmap2, float ptot)
{
    if (ptot < 0.15) // tracks only in the inner sector. R=1m
    {
        *tmap1 = 0;
        *tmap2 = 0;
        for (int iBit = 0; iBit < 32; iBit++)
        {
            if (mRandy->Rndm() < 0.95)
                *tmap1 |= 1 << iBit;
            else
                *tmap1 |= 0 << iBit;
        }
        *tmap1 |= 1 << 0; // default set to the first bit: primary-vertex-used
    }
    else
    {
        if (ptot < 0.3) // tracks may come to the end of the outer sector: R = 2m
        { 
            *tmap1 = 0;
            *tmap2 = 0;
            for (int iBit = 0; iBit < 32; iBit++)
            {
                if (mRandy->Rndm() < 0.95)
                    *tmap1 |= 1 << iBit;
                else
                    *tmap1 |= 0 << iBit;

                if (iBit < 27)
                {
                    if (mRandy->Rndm() < 0.85)
                        *tmap2 |= 1 << iBit;
                    else
                        *tmap2 |= 0 << iBit;
                }
            }
            *tmap1 |= 1 << 0; // default set to the first bit: primary-vertex-used
        }
        else
        {
            *tmap1 = 0;
            *tmap2 = 0;
            for (int iBit = 0; iBit < 32; iBit++)
            {
                if (mRandy->Rndm() < 0.95)
                    *tmap1 |= 1 << iBit;
                else
                    *tmap1 |= 0 << iBit;

                if(iBit < 27)
                {
                    if(mRandy->Rndm() < 0.95)
                        *tmap2 |= 1 << iBit;
                    else
                        *tmap2 |= 0 << iBit;
                }
            }
            *tmap1 |= 1 << 0; // default set to the first bit: primary-vertex-used
        }
    }
}