my_star/StRoot/StFemtoDstMaker/StO97DstQAMaker.cxx

#include "StO97DstQAMaker.h"

ClassImp(StO97DstQAMaker)

//_________________
StO97DstQAMaker::StO97DstQAMaker(const char *dirName,
                                 const char *fileName,
                                 const char *filter,
                                 int        maxFiles) {
    mDir = string(dirName);
    mFileName = string(fileName);
    mFilter = string(filter);

    mTree     = 0;
    mChain    = 0;
    mO97Event = 0;
    mEventId  = 0;
}


//_________________
Int_t StO97DstQAMaker::Init() {
    std::cout << "[StO97DstQAMaker]: Initialization... ";
    InitRead(mDir, mFileName, mFilter, mMaxFiles);
    mNEvents = (unsigned int)mChain->GetEntries();
    mOutFile = new TFile(mOutFileName, "RECREATE");
    if (!mOutFile) {
        std::cout << "FAIL" << std::endl;
        return kStErr;
    }
    else
        std::cout << "OK" << std::endl;
    //
    // Load PDG database
    //
    pdgDb = new TDatabasePDG;
    pdgDb->ReadPDGTable("StRoot/StFemtoDstMaker/pdg_table.txt");
    //
    // Creating event histograms
    //
		int refMultBins = 500;
		float refMultLo = 0.;
		float refMultHi = 500;
    hRefMultP90  = new TH1F("hRefMultP90", "reference multiplicity (greg) p_{#perp} > 0.09 GeV/c;RefMult;#",  refMultBins, refMultLo, refMultHi);
    hRefMultP100 = new TH1F("hRefMultP100", "reference multiplicity (greg) p_{#perp} > 0.1 GeV/c;RefMult;#",  refMultBins, refMultLo, refMultHi);
    hRefMultP120 = new TH1F("hRefMultP120", "reference multiplicity (greg) p_{#perp} > 0.12 GeV/c;RefMult;#", refMultBins, refMultLo, refMultHi);
    hRefMultP150 = new TH1F("hRefMultP150", "reference multiplicity (greg) p_{#perp} > 0.15 GeV/c;RefMult;#", refMultBins, refMultLo, refMultHi);
    hImpPar      = new TH1F("hImpPar", "impact parameter;imp (fm/c);#", 100, 0., 20.);
    hSph         = new TH1F("hSph",
                            "Transverse sphericity (|#eta| < 0.5, p_{#perp} > 0.15 GeV/c);S_{#perp},#frac{dN}{dS_{#perp}}",
                            100, 0., 1.);
    hSph         = new TH1F("hSph",
                            "Transverse sphericity (|#eta| < 0.5, p_{#perp} > 0.15 GeV/c);S_{#perp},#frac{dN}{dS_{#perp}}",
                            100, 0., 1.);
    hSph2        = new TH1F("hSph2",
                            "Transverse sphericity (|#eta| < 1.0, p_{#perp} > 0.15 GeV/c);S_{#perp},#frac{dN}{dS_{#perp}}",
                            100, 0., 1.);
		hImpVsRef    = new TH2F("hImpVsRef", "p_{#perp} > 0.15 [GeV/c] && |#eta| < 0.5;b [fm];RefMult", 100, 0., 20., refMultBins, refMultLo, refMultHi);

    //
    // Creating track histograms
    //
    hVx = new TH1F("hVx", "x freeze-out position;V_{x} (fm);#", 500, -50., 50.);
    hVy = new TH1F("hVy", "y freeze-out position;V_{y} (fm);#", 500, -50., 50.);
    hVz = new TH1F("hVz", "z freeze-out position;V_{z} (fm);#", 500, -50., 50.);
    hPx = new TH1F("hPx", "single particle p_{x};p_{x} (GeV/c);#", 200, -2., 2.);
    hPy = new TH1F("hPy", "single particle p_{y};p_{y} (GeV/c);#", 200, -2., 2.);
    hPz = new TH1F("hPz", "single particle p_{z};p_{z} (GeV/c);#", 200, -2., 2.);
    hP  = new TH1F("hP", "single particle p;p (GeV/c);#", 400, 0., 2.);
    hPt = new TH1F("hPt", "single particle p_{T};p_{T} (GeV/c);#", 400, 0., 2.);  
    hEta = new TH1F("hEta", ";#eta;#", 200, -10., 10.);
    hEtaInt = new TH1F("hEtaInt", ";#eta;#", 200, -10., 10.);  
    hMSqrVsPt = new TH2F("hMSqrVsPt",
                         "m^{2} vs p_{T};p_{T}*charge (GeV/c);m^{2} (GeV^{2})",
                         200, -2., 2., // Pt
                         100, 0., 2.); // msqr
    hMomElectron = new TH1F("hMomElectron", "Momentum of e^{#pm};p (GeV/c);#frac{dN}{dp}",
                            250, 0., 2.5);
    hMomPion     = new TH1F("hMomPion", "Momentum of #pi^{#pm};p (GeV/c);#frac{dN}{dp}",
                            250, 0., 2.5);
    hMomKaon     = new TH1F("hMomKaon", "Momentum of K^{#pm};p (GeV/c);#frac{dN}{dp}",
                            250, 0., 2.5);
    hMomProton   = new TH1F("hMomProton", "Momentum of p^{#pm};p (GeV/c);#frac{dN}{dp}",
                            250, 0., 2.5);

    hRfrVsTfr_pion = new TH2F("hRfrVsTfr_pion", ";R^{#pi^{ch}}_{freeze out} [fm];#tau^{#pi^{ch}}_{freeze out} [fm/c]",
                              300, 0., 15.,
                              300, 0., 15.);
    hRfrVsTfr_kaon = new TH2F("hRfrVsTfr_kaon", ";R^{K^{ch}}_{freeze out} [fm];#tau^{K^{ch}}_{freeze out} [fm/c]",
                              300, 0., 15.,
                              300, 0., 15.);
    hRfrVsTfr_proton = new TH2F("hRfrVsTfr_proton", ";R^{p^{#pm}}_{freeze out} [fm];#tau^{p^{#pm}}_{freeze out} [fm/c]",
                                300, 0., 15.,
                                300, 0., 15.);

    hdNdTfr_pion = new TH1F("hdNdTfr_pion", ";#frac{dN}{d#tau^{#pi^{ch}_{freezeout}}};#tau^{#pi^{ch}}_{freezeout} [fm/c]",
                            1000, 0., 1000.);
    hdNdTfr_kaon = new TH1F("hdNdTfr_kaon", ";#frac{dN}{d#tau^{K^{ch}_{freezeout}}};#tau^{K^{ch}}_{freezeout} [fm/c]",
                            1000, 0., 1000.);
    hdNdTfr_proton = new TH1F("hdNdTfr_proton", ";#frac{dN}{d#tau^{p^{ch}_{freezeout}}};#tau^{p^{ch}}_{freezeout} [fm/c]",
                              1000, 0., 1000.);

    return kStOk;
}


//_________________
Int_t StO97DstQAMaker::Finish() {
    UninitRead();

    mOutFile->Write();
    mOutFile->Close();

    delete pdgDb;

    return kStOk;
}


//_________________
int StO97DstQAMaker::FillChain(TChain *chain, char *dir,
                               const char *filter, int maxFiles) {
    TSystem *gSystem = new TSystem();
    if (!gSystem) {
        std::cout << "[StO97DstQAMaker]: can't allocate memory for TSystem" << std::endl;;
        return 0;
    }

    void *gDir = gSystem->OpenDirectory(dir);
    if (!gDir) {
        std::cout << "[StO97DstQAMaker]: can't open directory " 
            << dir << std::endl;
        delete gSystem;
        return 0;
    }

    int count = 0;
    char *file;

    while ((file = (char *)gSystem->GetDirEntry(dir))) {
        if (strcmp(file, ".") == 0 || strcmp(file, "..") == 0)
            continue;
        if (strstr(file, filter)) {
            char *gFullName = gSystem->ConcatFileName(dir, file);
            std::cout << "[StO97DstQAMaker]: Adding "
                << gFullName << " to the chain" << std::endl;
            chain->Add(gFullName);
            delete gFullName;
            count++;
            if ((maxFiles > 0) && (count > maxFiles)) break;
        }
    }

    std::cout << "[StO97DstQAMaker]: Added " << count 
        << " files to the chain" << std::endl;
    delete gSystem;
    return count;
}


//_________________
int StO97DstQAMaker::FillChain(TChain *chain, const char *fileName,
                               int maxFiles) {
    std::ifstream inStream(fileName);


    if (!inStream.is_open()) {
        std::cout << "[StFemtoDstQAMaker]: can't open file list: "
            << fileName << std::endl;
        return 0;
    }
    std::cout << "[StFemtoDstQAMaker]: using file list: "
        << fileName << std::endl;

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

    std::cout << "[StO97DstQAMaker]: Added " << count
        << " files to the chain" << std::endl;
    return count;
}



//_________________
int StO97DstQAMaker::InitRead(string dir, string fileName,
                              string filter, int maxFiles) {
    mEventId = 0;

    mChain = new TChain("StO97Dst", "StO97Dst");

    std::cout << "[StO97DstQAMaker]: Call InitRead" << std::endl;

    int numFiles = 0;
    if (!fileName.empty()) {
        if (strstr(fileName.c_str(), ".lis") ||
            strstr(fileName.c_str(), ".list")) {
            numFiles = FillChain(mChain, 
                                 (dir + fileName).c_str(), 
                                 maxFiles);
        }
        else {
            mChain->Add((dir + fileName).c_str());
            numFiles++;
        }
    } else {
        numFiles = FillChain(mChain, (char *)dir.c_str(), // not cool
                             filter.c_str(), maxFiles);
    }
    mChain->SetBranchAddress("StO97Event", &mO97Event);
    mNEvents = (unsigned int)mChain->GetEntries();
    return numFiles;
}


//_________________
void StO97DstQAMaker::UninitRead() {
    if (mChain)     delete mChain;
    mO97Event = 0;
    mChain    = 0;
}


//_________________
int StO97DstQAMaker::GetNEvents() {
    if (mChain)
        return mNEvents;
    else
        return -1;
}


//_________________
Int_t StO97DstQAMaker::Make() {
    if (!mNEvents) {
        std::cout << "[StO97DstQAMaker]: no events to read" << std::endl;
        return 0;
    }

    mO97Event->Clear();
    int bytes = mChain->GetEntry(mEventId++);

    if (mNEvents < mEventId) {
        std::cout << "[StO97DstQAMaker]: EOF" << std::endl;
        return 0;
    }

    if (!bytes) {
        std::cout << "[StO97DstQAMaker]: no event" << std::endl;
        return 0;
    }

    if (mO97Event) {
        int refMultP90 = 0;
        int refMultP100 = 0;
        int refMultP120 = 0;
        int refMultP150 = 0;
        float imp = mO97Event->GetImpactPar();

        TClonesArray *tracks = mO97Event->GetTracks();

        if (tracks) {
            int nTracks = tracks->GetEntries();
            for (int i = 0; i < nTracks; i++) {
                StO97Track *track = (StO97Track *)tracks->At(i);
                float px = track->GetPx();
                float py = track->GetPy();
                float pz = track->GetPz();
                float pt = sqrt(px*px + py*py);
                float p = sqrt(px*px + py*py + pz*pz);
                float massSqr = track->GetMass()*track->GetMass();
                float vx = track->GetXfr();
                float vy = track->GetYfr();
                float vz = track->GetZfr();
                float Rfr = TMath::Sqrt(vx*vx + vy*vy);
                float Tfr = track->GetTfr();
                float eta = pt > 0.05 /* i.e. pt > 0 */ ? 0.5*log((p + pz)/(p - pz)) : -999.;
                int   pdg = track->GetPdgId();
                int   charge = track->GetCharge(pdgDb);

                hEtaInt->Fill(eta);

                if (pt < 0.075 || fabs(eta) >= 0.5 /* STAR TPC acceptance */ ||
                    charge == 0. /* only charged particles */ ||
                    track->GetIsSpec()) {

                    continue;
                }

                if (pt >= 0.09)   refMultP90++;
                if (pt >= 0.10)   refMultP100++;
                if (pt >= 0.12)   refMultP120++;
                if (pt >= 0.15)   refMultP150++;

                if (pt >= 0.1 /* STAR TPC lower momentum limit from STAR TPC NIM */) {
                    hMSqrVsPt->Fill(pt*charge, massSqr);
                    hVx->Fill(vx);
                    hVy->Fill(vy);
                    hVz->Fill(vz);
                    hPx->Fill(px);
                    hPy->Fill(py);
                    hPz->Fill(pz);
                    hPt->Fill(pt);
                    hP->Fill(p);
                    hEta->Fill(eta);
                }

                switch (abs(pdg)) {
                    case 11: // Electron
                        hMomElectron->Fill(p);
                        break;
                    case 211: // Pion
                        hMomPion->Fill(p);
                        hdNdTfr_pion->Fill(Tfr);
                        hRfrVsTfr_pion->Fill(Rfr, Tfr);
                        break;
                    case 321: // Kaon
                        hMomKaon->Fill(p);
                        hdNdTfr_kaon->Fill(Tfr);
                        hRfrVsTfr_kaon->Fill(Rfr, Tfr);
                        break;
                    case 2212: // Proton
                        hMomProton->Fill(p);
                        hdNdTfr_proton->Fill(Tfr);
                        hRfrVsTfr_proton->Fill(Rfr, Tfr);
                } // switch (abs(pdg))
            } // for (int i = 0; i < nTracks; i++)
        } // if (tracks)
        else {
            std::cout << "[StO97DstQAMaker]: tracks == NULL" << std::endl;
        }

        if (refMultP90)
            hRefMultP90->Fill(refMultP90);
        if (refMultP100)
            hRefMultP100->Fill(refMultP100);
        if (refMultP120)
            hRefMultP120->Fill(refMultP120);
        if (refMultP150)
            hRefMultP150->Fill(refMultP150);

        hSph->Fill(mO97Event->GetTransverseSphericity());
        hSph2->Fill(mO97Event->GetTransverseSphericity2());

        hImpPar->Fill(imp);
				hImpVsRef->Fill(imp, refMultP150);
    } // if (mO97Event)

    return kStOk;
}