my_star/StRoot/StFemtoDstMaker/StFemtoDstMaker_StarGen.cxx

#include "StFemtoDstMaker_StarGen.h"
#include "StBTofHeader.h"
#include "TBranch.h"
#include "StPhysicalHelixD.hh"
#include "StDcaGeometry.h"

ClassImp(StFemtoDstMaker_StarGen)

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



StFemtoDstMaker_StarGen::StFemtoDstMaker_StarGen(const char *dirName, const char *fileName,
                                                 const char *oFileName,
                                                 const char *filter, int maxFiles)
{
    std::cout << "StFemtoDstMaker_StarGen::StFemtoDstMaker_StarGen - Creating an instance...";

    mOutFileName = oFileName;
    mCompression = 9;

    // initialize single-particle cut variables
    mTrackMomentum[0] = 0.1;
    mTrackMomentum[1] = 2.;
    mTrackEta[0] = -1.1;
    mTrackEta[1] = 1.1;

    mDir      = string(dirName);
    mFileName = string(fileName);
    mFilter   = filter;
    mMaxFiles = maxFiles;
    mDb       = TDatabasePDG::Instance();
    mRandy    = new TRandom3(0);
    mMatrix   = new TMatrixTSym<double>(2);

    mTotalTracks = 0.;

    mTree       = 0;
    mTChain     = 0;
    mSGEvent    = 0;
    mEventIndex = 0;

    InitRead(mDir, mFileName, mFilter, mMaxFiles);

    mNEvents = (unsigned int)mTChain->GetEntries();
    std::cout << "\t[DONE]\n";
}




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



int StFemtoDstMaker_StarGen::FillChain(TChain *chain, char *dir,
                                       const char *filter, int maxFiles)
{
    TSystem *gSystem = new TSystem();
    if (!gSystem)
    {
        cout << "StFemtoDstMaker_StarGen[ERROR]: cannot allocate memory for TSystem\n";
        return 0;
    }
    void *lDir = gSystem->OpenDirectory(dir);
    if (!lDir)
    {
        cout
            << "StFemtoDstMaker_StarGen[ERROR]: can't open directory "
            << dir << endl;
        delete gSystem;
        return 0;
    }

    int mCount = 0;
    char *lFileName;

    while ((lFileName = (char *)gSystem->GetDirEntry(dir)))
    {
        if (strcmp(lFileName, ".") == 0 || strcmp(lFileName, "..") == 0)
            continue;
        if (strstr(lFileName, filter))
        {
            char *lFullName = gSystem->ConcatFileName(dir, lFileName);
            cout
                << "StFemtoDstMaker_StarGen[INFO]: Adding "
                << lFullName << " to the chain\n";
            chain->Add(lFullName);
            delete lFullName;
            mCount++;
            if ((maxFiles > 0) && (mCount > maxFiles)) break;
        }
    }

    cout
        << "StFemtoDstMaker_StarGen[INFO]: Added " << mCount
        << " files to the chain\n";
    delete gSystem;
    return mCount;
}



int StFemtoDstMaker_StarGen::FillChain(TChain *chain, const char *fileName, int maxFiles)
{
    ifstream lInStream(fileName);

    if (!lInStream.is_open())
    {
        cout
            << "StFemtoDstMaker_StarGen[ERROR]: can't open file list: "
            << fileName << endl;
        return 0;
    }
    cout
        << "StFemtoDstMaker_StarGen[INFO]: using file list: "
        << fileName << endl;

    int mCount = 0;
    char buf[0xFF];
    for ( ; lInStream.good(); )
    {
        lInStream.getline(buf, 0xFF);
        TString lFileName(buf);
        if (lFileName.Contains("root"))
        {
            chain->Add(buf);
            mCount++;
            if ((maxFiles > 0) && (mCount > maxFiles)) break;
        }
    }

    cout
        << "StFemtoDstMaker_StarGen[INFO]: Added " << mCount
        << " files to the chain\n";
    return mCount;
}



int StFemtoDstMaker_StarGen::InitRead(string aDir, string fileName,
                                      string aFilter, int maxFiles)
{
    mEventIndex = 0;
    mTChain = new TChain("genevents", "genevents");

    cout << "StFemtoDstMaker_StarGen[INFO]: Call InitRead\n";

    int lNumFiles = 0;
    if (!fileName.empty())
    {
        if (strstr(fileName.c_str(), ".lis") ||
            strstr(fileName.c_str(), ".list"))
        {
            lNumFiles = FillChain( mTChain,
                                  (aDir + fileName).c_str(),
                                  maxFiles );
        }
        else
        {
            mTChain->Add((aDir + fileName).c_str());
            lNumFiles++;
        }
    } else
    {
        lNumFiles = FillChain(mTChain, (char *)aDir.c_str(), // not cool
                              aFilter.c_str(), maxFiles);
    }
    mTChain->SetBranchAddress("primaryEvent", &mSGEvent);
    mNEvents = (unsigned int)mTChain->GetEntries();
    cout
        << "StFemtoDstMaker_StarGen::InitRead [INFO]"
        << " NEvents to read: " << mNEvents << endl;
    return lNumFiles;
}



void StFemtoDstMaker_StarGen::UninitRead()
{
    if (mTChain)   delete mTChain;
    mSGEvent = 0;
    mTChain  = 0;
}



int StFemtoDstMaker_StarGen::GetNEvents()
{
    if (mTChain)   return mNEvents;
    else           return -1;
}



Int_t StFemtoDstMaker_StarGen::Init()
{
    std::cout << "StFemtoDstMaker_StarGen::Init - Initializing the maker\n"
              << "Creating output file: " << mOutFileName;

    mFemtoEvent = new StFemtoEvent();
    mOutFile = new TFile(mOutFileName, "RECREATE");
    mOutFile->cd();
    mOutFile->SetCompressionLevel(mCompression);
    std::cout << "\t[DONE]" << std::endl;

    mTree = new TTree("StFemtoDst","StFemtoDst");
    mTree->SetMaxTreeSize(MAXFILESIZE);
    mTree->Branch("StFemtoEvent","StFemtoEvent", &mFemtoEvent);

    mNBytes = 0;
    std::cout << "StFemtoDstMaker_StarGen::Init - Initialization has been finished\n";
    return StMaker::Init();
}



void StFemtoDstMaker_StarGen::Clear(Option_t *option)
{
    StMaker::Clear();
}



Int_t StFemtoDstMaker_StarGen::Make()
{
    if (!mNEvents)
    {
        cout << "StHbtStarGenReader[INFO]: no events to read\n";
        return 0;
    }

    int lBytes = mTChain->GetEntry(mEventIndex++);

    if (mNEvents < mEventIndex)
    {
        cout << "StHbtStarGenReader[INFO]: EOF\n";
        return 0;
    }

    if (!lBytes)
    {
        cout << "StHbtStarGenReader[INFO]: no event\n";
        return 0;
    }

    if (!mSGEvent)   return 0;

    const float magRFF = -4.98;   // [10^3 G], ref: p+Au 200 GeV y2015
    unsigned int nTracks = mSGEvent->GetNumberOfParticles(); // number of tracks in an event

    mFemtoEvent->SetEventId(mSGEvent->GetEventNumber());
    mFemtoEvent->SetRunId(mSGEvent->GetRunNumber());
    mFemtoEvent->SetCollisionType(false);
    mFemtoEvent->SetCent16(0);
    mFemtoEvent->SetMagneticField(magRFF);
    mFemtoEvent->SetPrimaryVertexPosition(0.0, 0.0, 0.0);
    mFemtoEvent->SetVpdVz(0.0);
    mFemtoEvent->SetTriggerId(0);
    mFemtoEvent->SetPrimaryVertexRanking(666);

    unsigned int refMult2    = 0;
    unsigned int refMult2Pos = 0;
    unsigned int refMult     = 0;
    unsigned int refMultPos  = 0;
    unsigned int nGoodTracks = 0;

    StFemtoTrack *mFTrack = NULL;

    //
    // For sphericity
    //
    float sph     = -999.;
    float pt_full = 0.;
    mMatrix->Zero();

    //
    // Loop over primary tracks
    //
    for (unsigned int iTrk = 0; iTrk < nTracks; iTrk++)
    {
        StarGenParticle *sParticle = (*mSGEvent)[iTrk];
        //
        // Get charge
        //
        int                   pdg     = sParticle->GetId();
        TParticlePDG         *partPdg = mDb->GetParticle(pdg);
        if (partPdg == 0x0)   continue;
        char charge = partPdg->Charge();
        if (charge == 0)      continue;
        charge /= abs(charge);

        //
        // Get momentum, beta and DCA
        //
        const TLorentzVector &p = sParticle->momentum();
        float px     = p.Px();
        float py     = p.Py();
        float pz     = p.Pz();
        float pt     = p.Pt();
        if (pt < 1.0e-5)   continue;
        float ptot   = p.P();
        float eta    = p.PseudoRapidity();
        float energy = p.Energy();
        float beta   = ptot/energy;
        float xDCA   = sParticle->GetVx();
        float yDCA   = sParticle->GetVy();
        float zDCA   = sParticle->GetVz();
        float dca    = sqrt(xDCA*xDCA + yDCA*yDCA + zDCA*zDCA);

        //
        // Calculate sphericity
        //
        // find lead particle and calculate matrix
        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;
            }
        }

        if (fabs(eta) < 1.0 && dca < 3.0)
        {
            refMult2++;
            if (charge > 0)   refMult2Pos++;
            if (fabs(eta) < 0.5)
            {
                refMult++;
                if (charge > 0)   refMultPos++;
            }
        }

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

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

        if (!AcceptTrack(p))   continue;
        else                   nGoodTracks++;

        //
        // Calculate topology map
        //
        unsigned int tmap1 = 0;
        unsigned int tmap2 = 0;
        TopologyMap(&tmap1, &tmap2, ptot);

        // add femto track
        mFTrack = mFemtoEvent->AddFemtoTrack();

        mFTrack->SetId(iTrk);
        mFTrack->SetNHits(charge*45);
        mFTrack->SetNHitsPoss(45);
        mFTrack->SetNSigmaElectron(nSigmaElectron);
        mFTrack->SetNSigmaPion(nSigmaPion);
        mFTrack->SetNSigmaKaon(nSigmaKaon);
        mFTrack->SetNSigmaProton(nSigmaProton);
        mFTrack->SetDedx(dedxMean(sParticle->GetMass(), ptot));
        mFTrack->SetMap0(tmap1);
        mFTrack->SetMap1(tmap2);
        mFTrack->SetP(px, py, pz);
        mFTrack->SetDCAxGlobal(sParticle->GetVx());
        mFTrack->SetDCAyGlobal(sParticle->GetVy());
        mFTrack->SetDCAzGlobal(sParticle->GetVz());
        mFTrack->SetGlobalP(px, py, pz);
        mFTrack->SetBeta(beta);
    }

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

    //
    // Continue to fill femto event information
    //
    mFemtoEvent->SetSphericity(sph);
    mFemtoEvent->SetRefMult(refMult);
    mFemtoEvent->SetRefMultCorr(refMult);
    mFemtoEvent->SetRefMultCorrWeight(1.0);
    mFemtoEvent->SetRefMultPos(refMultPos);
    mFemtoEvent->SetNumberOfBTofHit(refMult);
    mFemtoEvent->SetNumberOfPrimaryTracks(nGoodTracks);
    mFemtoEvent->SetNumberOfGlobalTracks(nGoodTracks);
    mFemtoEvent->SetRefMult2(refMult2);
    mFemtoEvent->SetRefMult2Pos(refMult2Pos);

    mNBytes += mTree->Fill();
    mFemtoEvent->Clear();

    return kStOk;
}



bool StFemtoDstMaker_StarGen::AcceptTrack(const TLorentzVector &p)
{
    return p.P() >= mTrackMomentum[0] && p.P() <= mTrackMomentum[1] &&
        p.PseudoRapidity() >= mTrackEta[0] && p.PseudoRapidity() <= mTrackEta[1];
}



Int_t StFemtoDstMaker_StarGen::Finish()
{
    mOutFile->cd();
    mOutFile->Write();
    mOutFile->Close();

    std::cout
        << "***\n"
        << "*\n"
        << "*  StFemtoDstMaker_StarGen has been finished\n"
        << "*    bytes has been written: " << mNBytes << "\n"
        << "***\n";

    return kStOk;
}



void StFemtoDstMaker_StarGen::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
        }
    }
}



float StFemtoDstMaker_StarGen::dedxMean(float mass, float momentum)
{
    double dedxMean;
    double tpcDedxGain   = 0.174325e-06;
    double tpcDedxOffset = -2.71889;
    double tpcDedxRise   = 776.626;

    double gamma = TMath::Sqrt(momentum*momentum/(mass*mass) + 1.0);
    double beta  = TMath::Sqrt(1.0 - 1.0/(gamma*gamma));
    double rise  = tpcDedxRise* beta*beta*gamma*gamma;

    if (beta > 0.0)
        dedxMean = tpcDedxGain/(beta*beta)*(0.5*TMath::Log(rise) - beta*beta
                                            - tpcDedxOffset);
    else
        dedxMean = 1000.0;

    return dedxMean;

}