123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481 |
- #include <iostream>
- #include <TString.h>
- #include <TH1F.h>
- #include <TH1I.h>
- #include <TH2F.h>
- #include <TProfile.h>
- #include <TFile.h>
- #include <TStopwatch.h>
- #include <qaParticle.h>
- #include <qaEvent.h>
- #include <qaReader_manager.h>
- #include <qaReader_smash_root.h>
- #include <qaReader_mcpico.h>
- #include <Utility.h>
- #ifdef _MCINI_
- #include <qaReader_mcini.h>
- #endif
- int main(int argc, char **argv)
- {
- TString iFileName, oFileName, configFileName = "";
- if (argc < 7)
- {
- #ifndef _MCINI_
- std::cerr << "./qaProcess -i input.list -o qa_output.root -format [mcpico/particles] [OPTIONAL: -config qa.cfg]" << std::endl;
- #endif
- #ifdef _MCINI_
- std::cerr << "./qaProcess -i input.list -o qa_output.root -format [mcpico/mcini/particles] [OPTIONAL: -config qa.cfg]" << std::endl;
- #endif
- return 1;
- }
- for (int i = 1; i < argc; i++)
- {
- if (std::string(argv[i]) != "-i" &&
- std::string(argv[i]) != "-o" &&
- std::string(argv[i]) != "-format" &&
- std::string(argv[i]) != "-config")
- {
- std::cerr << "\n[ERROR]: Unknown parameter " << i << ": " << argv[i] << std::endl;
- return 2;
- }
- else
- {
- if (std::string(argv[i]) == "-i" && i != argc - 1)
- {
- iFileName = argv[++i];
- continue;
- }
- if (std::string(argv[i]) == "-i" && i == argc - 1)
- {
- std::cerr << "\n[ERROR]: Input file name was not specified " << std::endl;
- return 1;
- }
- if (std::string(argv[i]) == "-o" && i != argc - 1)
- {
- oFileName = argv[++i];
- continue;
- }
- if (std::string(argv[i]) == "-o" && i == argc - 1)
- {
- std::cerr << "\n[ERROR]: Output file name was not specified " << std::endl;
- return 1;
- }
- if (std::string(argv[i]) == "-format" && i != argc - 1)
- {
- qaUtility::GetInstance()->format = argv[++i];
- continue;
- }
- if (std::string(argv[i]) == "-format" && i == argc - 1)
- {
- std::cerr << "\n[ERROR]: Output file name was not specified " << std::endl;
- return 1;
- }
- if (std::string(argv[i]) == "-config" && i != argc - 1)
- {
- configFileName = argv[++i];
- continue;
- }
- if (std::string(argv[i]) == "-config" && i == argc - 1)
- {
- std::cerr << "\n[ERROR]: Output file name was not specified " << std::endl;
- return 1;
- }
- }
- }
- TStopwatch timer;
- timer.Start();
- if (configFileName.Length() > 0)
- {
- qaUtility::GetInstance()->ReadConfig(configFileName);
- }
- if (qaUtility::GetInstance()->debug)
- {
- std::cout << "Format = " << qaUtility::GetInstance()->format << std::endl;
- std::cout << "Configuration:" << std::endl;
- std::cout << "debug = " << qaUtility::GetInstance()->debug << std::endl;
- std::cout << "Nevents = " << qaUtility::GetInstance()->Nevents << std::endl;
- std::cout << "Cut_minbias_Event_bmin = " << qaUtility::GetInstance()->Cut_minbias_Event_bmin << std::endl;
- std::cout << "Cut_minbias_Event_bmax = " << qaUtility::GetInstance()->Cut_minbias_Event_bmax << std::endl;
- std::cout << "Cut_minbias_Particle_ptmin = " << qaUtility::GetInstance()->Cut_minbias_Particle_ptmin << std::endl;
- std::cout << "Cut_minbias_Particle_ptmax = " << qaUtility::GetInstance()->Cut_minbias_Particle_ptmax << std::endl;
- std::cout << "Cut_minbias_Particle_etamin = " << qaUtility::GetInstance()->Cut_minbias_Particle_etamin << std::endl;
- std::cout << "Cut_minbias_Particle_etamax = " << qaUtility::GetInstance()->Cut_minbias_Particle_etamax << std::endl;
- std::cout << "Cut_refmult_Event_bmin = " << qaUtility::GetInstance()->Cut_refmult_Event_bmin << std::endl;
- std::cout << "Cut_refmult_Event_bmax = " << qaUtility::GetInstance()->Cut_refmult_Event_bmax << std::endl;
- std::cout << "Cut_refmult_Particle_ptmin = " << qaUtility::GetInstance()->Cut_refmult_Particle_ptmin << std::endl;
- std::cout << "Cut_refmult_Particle_ptmax = " << qaUtility::GetInstance()->Cut_refmult_Particle_ptmax << std::endl;
- std::cout << "Cut_refmult_Particle_etamin = " << qaUtility::GetInstance()->Cut_refmult_Particle_etamin << std::endl;
- std::cout << "Cut_refmult_Particle_etamax = " << qaUtility::GetInstance()->Cut_refmult_Particle_etamax << std::endl;
- std::cout << std::endl;
- }
- TFile *fo = new TFile(oFileName, "recreate");
- TH1F *h_minbias_Event_b = new TH1F("h_minbias_Event_b", "dN/db minbias;b, fm; dN/db", 200, 0., 20.);
- TH1I *h_minbias_Event_Mult = new TH1I("h_minbias_Event_Mult", "dN/dN_{particles} minbias;N_{particles}; dN/dN_{particles}", 2500, 0, 2500);
- TH2F *h2_minbias_Event_b_Mult = new TH2F("h2_minbias_Event_b_Mult", "b vs N_{particles} minbias;N_{particles};b, fm;", 2500, 0., 2500., 200, 0., 20.);
- TH1F *h_minbias_Particle_pt = new TH1F("h_minbias_Particle_pt", "dN/dp_{T} minbias;p_{T}, GeV/c; dN/dp_{T}", 800, 0., 8.);
- TH1F *h_minbias_Particle_eta = new TH1F("h_minbias_Particle_eta", "dN/d#eta minbias;#eta; dN/d#eta", 2000, -10., 10.);
- TH1F *h_minbias_Particle_E = new TH1F("h_minbias_Particle_E", "dN/dE minbias;E, GeV; dN/dE", 800, 0., 8.);
- TH1F *h_minbias_Particle_px = new TH1F("h_minbias_Particle_px", "dN/dp_{x} minbias;p_{x}, GeV/c; dN/dp_{x}", 800, -8., 8.);
- TH1F *h_minbias_Particle_py = new TH1F("h_minbias_Particle_py", "dN/dp_{y} minbias;p_{y}, GeV/c; dN/dp_{y}", 800, -8., 8.);
- TH1F *h_minbias_Particle_pz = new TH1F("h_minbias_Particle_pz", "dN/dp_{z} minbias;p_{z}, GeV/c; dN/dp_{z}", 800, -8., 8.);
- TH1F *h_minbias_Particle_t = new TH1F("h_minbias_Particle_t", "dN/dt minbias;t, fm/c; dN/dt", 800, 0., 800.);
- TH1F *h_minbias_Particle_x = new TH1F("h_minbias_Particle_x", "dN/dx minbias;x, fm; dN/dx", 800, -8., 8.);
- TH1F *h_minbias_Particle_y = new TH1F("h_minbias_Particle_y", "dN/dy minbias;y, fm; dN/dy", 800, -8., 8.);
- TH1F *h_minbias_Particle_z = new TH1F("h_minbias_Particle_z", "dN/dz minbias;z, fm; dN/dz", 800, -8., 8.);
- TH1I *h_minbias_Particle_pdg = new TH1I("h_minbias_Particle_pdg", "PDG minbias;PDG;N_{entries}", 7000, -3500, 3500);
- TH1F *h_minbias_Particle_PID_pt[qaUtility::GetInstance()->npid];
- TH1F *h_minbias_Particle_PID_eta[qaUtility::GetInstance()->npid];
- TH1F *h_minbias_Particle_PID_E[qaUtility::GetInstance()->npid];
- TH1F *h_minbias_Particle_PID_px[qaUtility::GetInstance()->npid];
- TH1F *h_minbias_Particle_PID_py[qaUtility::GetInstance()->npid];
- TH1F *h_minbias_Particle_PID_pz[qaUtility::GetInstance()->npid];
- TH1F *h_minbias_Particle_PID_t[qaUtility::GetInstance()->npid];
- TH1F *h_minbias_Particle_PID_x[qaUtility::GetInstance()->npid];
- TH1F *h_minbias_Particle_PID_y[qaUtility::GetInstance()->npid];
- TH1F *h_minbias_Particle_PID_z[qaUtility::GetInstance()->npid];
- TH1F *h_refmult_Event_b = new TH1F("h_refmult_Event_b", "dN/db refmult;b, fm; dN/db", 200, 0., 20.);
- TH1I *h_refmult_Event_Mult = new TH1I("h_refmult_Event_Mult", "dN/dN_{particles} refmult;N_{particles}; dN/dN_{particles}", 2500, 0, 2500);
- TH2F *h2_refmult_Event_b_Mult = new TH2F("h2_refmult_Event_b_Mult", "b vs N_{particles} refmult;N_{particles};b, fm;", 2500, 0., 2500., 200, 0., 20.);
- TH1F *h_refmult_Particle_pt = new TH1F("h_refmult_Particle_pt", "dN/dp_{T} refmult;p_{T}, GeV/c; dN/dp_{T}", 800, 0., 8.);
- TH1F *h_refmult_Particle_eta = new TH1F("h_refmult_Particle_eta", "dN/d#eta refmult;#eta; dN/d#eta", 2000, -10., 10.);
- TH1F *h_refmult_Particle_E = new TH1F("h_refmult_Particle_E", "dN/dE refmult;E, GeV; dN/dE", 800, 0., 8.);
- TH1F *h_refmult_Particle_px = new TH1F("h_refmult_Particle_px", "dN/dp_{x} refmult;p_{x}, GeV/c; dN/dp_{x}", 800, -8., 8.);
- TH1F *h_refmult_Particle_py = new TH1F("h_refmult_Particle_py", "dN/dp_{y} refmult;p_{y}, GeV/c; dN/dp_{y}", 800, -8., 8.);
- TH1F *h_refmult_Particle_pz = new TH1F("h_refmult_Particle_pz", "dN/dp_{z} refmult;p_{z}, GeV/c; dN/dp_{z}", 800, -8., 8.);
- TH1F *h_refmult_Particle_t = new TH1F("h_refmult_Particle_t", "dN/dt refmult;t, fm/c; dN/dt", 800, 0., 800.);
- TH1F *h_refmult_Particle_x = new TH1F("h_refmult_Particle_x", "dN/dx refmult;x, fm; dN/dx", 800, -8., 8.);
- TH1F *h_refmult_Particle_y = new TH1F("h_refmult_Particle_y", "dN/dy refmult;y, fm; dN/dy", 800, -8., 8.);
- TH1F *h_refmult_Particle_z = new TH1F("h_refmult_Particle_z", "dN/dz refmult;z, fm; dN/dz", 800, -8., 8.);
- TH1I *h_refmult_Particle_pdg = new TH1I("h_refmult_Particle_pdg", "PDG refmult;PDG;N_{entries}", 7000, -3500, 3500);
- TH1F *h_refmult_Particle_PID_pt[qaUtility::GetInstance()->npid];
- TH1F *h_refmult_Particle_PID_eta[qaUtility::GetInstance()->npid];
- TH1F *h_refmult_Particle_PID_E[qaUtility::GetInstance()->npid];
- TH1F *h_refmult_Particle_PID_px[qaUtility::GetInstance()->npid];
- TH1F *h_refmult_Particle_PID_py[qaUtility::GetInstance()->npid];
- TH1F *h_refmult_Particle_PID_pz[qaUtility::GetInstance()->npid];
- TH1F *h_refmult_Particle_PID_t[qaUtility::GetInstance()->npid];
- TH1F *h_refmult_Particle_PID_x[qaUtility::GetInstance()->npid];
- TH1F *h_refmult_Particle_PID_y[qaUtility::GetInstance()->npid];
- TH1F *h_refmult_Particle_PID_z[qaUtility::GetInstance()->npid];
- for (int i=0; i<qaUtility::GetInstance()->npid; i++)
- {
- h_minbias_Particle_PID_pt[i] = new TH1F(Form("h_minbias_Particle_PID_pt_%i",i),Form("h_minbias_Particle_PID_pt_%i;p_{T}, GeV/c; dN/dp_{T}",i), 800, 0., 8.);
- h_minbias_Particle_PID_eta[i] = new TH1F(Form("h_minbias_Particle_PID_eta_%i",i),Form("h_minbias_Particle_PID_eta_%i;#eta; dN/d#eta",i), 2000, -10., 10.);
- h_minbias_Particle_PID_E[i] = new TH1F(Form("h_minbias_Particle_PID_E_%i",i),Form("h_minbias_Particle_PID_E_%i;E, GeV; dN/dE",i), 800, 0., 8.);
- h_minbias_Particle_PID_px[i] = new TH1F(Form("h_minbias_Particle_PID_px_%i",i),Form("h_minbias_Particle_PID_px_%i;p_{x}, GeV/c; dN/dp_{x}",i), 800, -8., 8.);
- h_minbias_Particle_PID_py[i] = new TH1F(Form("h_minbias_Particle_PID_py_%i",i),Form("h_minbias_Particle_PID_py_%i;p_{y}, GeV/c; dN/dp_{y}",i), 800, -8., 8.);
- h_minbias_Particle_PID_pz[i] = new TH1F(Form("h_minbias_Particle_PID_pz_%i",i),Form("h_minbias_Particle_PID_pz_%i;p_{z}, GeV/c; dN/dp_{z}",i), 800, -8., 8.);
- h_minbias_Particle_PID_t[i] = new TH1F(Form("h_minbias_Particle_PID_t_%i",i),Form("h_minbias_Particle_PID_t_%i;t, fm/c; dN/dt",i), 800, 0., 8.);
- h_minbias_Particle_PID_x[i] = new TH1F(Form("h_minbias_Particle_PID_x_%i",i),Form("h_minbias_Particle_PID_x_%i;x, fm; dN/dx",i), 800, -8., 8.);
- h_minbias_Particle_PID_y[i] = new TH1F(Form("h_minbias_Particle_PID_y_%i",i),Form("h_minbias_Particle_PID_y_%i;y, fm; dN/dy",i), 800, -8., 8.);
- h_minbias_Particle_PID_z[i] = new TH1F(Form("h_minbias_Particle_PID_z_%i",i),Form("h_minbias_Particle_PID_z_%i;z, fm; dN/dz",i), 800, -8., 8.);
- h_refmult_Particle_PID_pt[i] = new TH1F(Form("h_refmult_Particle_PID_pt_%i",i),Form("h_refmult_Particle_PID_pt_%i;p_{T}, GeV/c; dN/dp_{T}",i), 800, 0., 8.);
- h_refmult_Particle_PID_eta[i] = new TH1F(Form("h_refmult_Particle_PID_eta_%i",i),Form("h_refmult_Particle_PID_eta_%i;#eta; dN/d#eta",i), 2000, -10., 10.);
- h_refmult_Particle_PID_E[i] = new TH1F(Form("h_refmult_Particle_PID_E_%i",i),Form("h_refmult_Particle_PID_E_%i;E, GeV; dN/dE",i), 800, 0., 8.);
- h_refmult_Particle_PID_px[i] = new TH1F(Form("h_refmult_Particle_PID_px_%i",i),Form("h_refmult_Particle_PID_px_%i;p_{x}, GeV/c; dN/dp_{x}",i), 800, -8., 8.);
- h_refmult_Particle_PID_py[i] = new TH1F(Form("h_refmult_Particle_PID_py_%i",i),Form("h_refmult_Particle_PID_py_%i;p_{y}, GeV/c; dN/dp_{y}",i), 800, -8., 8.);
- h_refmult_Particle_PID_pz[i] = new TH1F(Form("h_refmult_Particle_PID_pz_%i",i),Form("h_refmult_Particle_PID_pz_%i;p_{z}, GeV/c; dN/dp_{z}",i), 800, -8., 8.);
- h_refmult_Particle_PID_t[i] = new TH1F(Form("h_refmult_Particle_PID_t_%i",i),Form("h_refmult_Particle_PID_t_%i;t, fm/c; dN/dt",i), 800, 0., 8.);
- h_refmult_Particle_PID_x[i] = new TH1F(Form("h_refmult_Particle_PID_x_%i",i),Form("h_refmult_Particle_PID_x_%i;x, fm; dN/dx",i), 800, -8., 8.);
- h_refmult_Particle_PID_y[i] = new TH1F(Form("h_refmult_Particle_PID_y_%i",i),Form("h_refmult_Particle_PID_y_%i;y, fm; dN/dy",i), 800, -8., 8.);
- h_refmult_Particle_PID_z[i] = new TH1F(Form("h_refmult_Particle_PID_z_%i",i),Form("h_refmult_Particle_PID_z_%i;z, fm; dN/dz",i), 800, -8., 8.);
- }
- qaReader_manager *readerManager;
- if (qaUtility::GetInstance()->format == "mcpico")
- {
- readerManager = new qaReader_mcpico();
- }
- #ifdef _MCINI_
- if (qaUtility::GetInstance()->format == "mcini")
- {
- readerManager = new qaReader_mcini();
- }
- #endif
- if (qaUtility::GetInstance()->format == "particles")
- {
- readerManager = new qaReader_smash_root();
- }
- if (!readerManager)
- {
- std::cerr << "This input format is not found!" << std::endl;
- return 20;
- }
- readerManager->SetChain(iFileName.Data());
- std::vector<Long64_t> vRejectedEvents;
- Long64_t Nentries_chain = readerManager->GetEntries();
- Long64_t Nentries = (qaUtility::GetInstance()->Nevents > Nentries_chain) ? Nentries_chain : qaUtility::GetInstance()->Nevents;
- if (qaUtility::GetInstance()->Nevents == -1) Nentries = Nentries_chain;
- Int_t Nparticles;
- Int_t Ncounter_minbias, Ncounter_refmult;
- qaEvent *event = nullptr;
- qaParticle *particle = nullptr;
- Long64_t Absolute_counter = 0;
- Long64_t Minbias_counter = 0;
- Int_t ipid;
- while(Minbias_counter <= Nentries)
- {
- if (Minbias_counter % 1000 == 0)
- std::cout << "Event [" << Minbias_counter << "/" << Nentries << "]" << std::endl;
- event = (qaEvent *)readerManager->ReadEvent(Absolute_counter);
- Absolute_counter++;
- if (Absolute_counter > Nentries_chain)
- break;
- if (!event)
- continue;
- Ncounter_minbias = 0;
- Ncounter_refmult = 0;
- Nparticles = event->GetNparticles();
- for (int iparticle = 0; iparticle < Nparticles; iparticle++)
- {
- particle = readerManager->ReadParticle(iparticle);
- if (!particle)
- continue;
- if (qaUtility::GetInstance()->Cut_Event_minbias(event) &&
- qaUtility::GetInstance()->Cut_Particle_minbias(particle))
- {
- Ncounter_minbias++;
- h_minbias_Particle_pt ->Fill(particle->GetPt());
- h_minbias_Particle_eta->Fill(particle->GetEta());
- h_minbias_Particle_E ->Fill(particle->GetEnergy());
- h_minbias_Particle_px ->Fill(particle->GetPx());
- h_minbias_Particle_py ->Fill(particle->GetPy());
- h_minbias_Particle_pz ->Fill(particle->GetPz());
- h_minbias_Particle_t ->Fill(particle->GetT());
- h_minbias_Particle_x ->Fill(particle->GetX());
- h_minbias_Particle_y ->Fill(particle->GetY());
- h_minbias_Particle_z ->Fill(particle->GetZ());
- h_minbias_Particle_pdg->Fill(particle->GetPdg());
- if (qaUtility::GetInstance()->GetCharge(particle->GetPdg()) > 0)
- {
- ipid = 0;
- }
- if (qaUtility::GetInstance()->GetCharge(particle->GetPdg()) < 0)
- {
- ipid = 4;
- }
- if (ipid == 0 || ipid == 4)
- {
- h_minbias_Particle_PID_pt[ipid] ->Fill(particle->GetPt());
- h_minbias_Particle_PID_eta[ipid]->Fill(particle->GetEta());
- h_minbias_Particle_PID_E[ipid] ->Fill(particle->GetEnergy());
- h_minbias_Particle_PID_px[ipid] ->Fill(particle->GetPx());
- h_minbias_Particle_PID_py[ipid] ->Fill(particle->GetPy());
- h_minbias_Particle_PID_pz[ipid] ->Fill(particle->GetPz());
- h_minbias_Particle_PID_t[ipid] ->Fill(particle->GetT());
- h_minbias_Particle_PID_x[ipid] ->Fill(particle->GetX());
- h_minbias_Particle_PID_y[ipid] ->Fill(particle->GetY());
- h_minbias_Particle_PID_z[ipid] ->Fill(particle->GetZ());
- }
- ipid = qaUtility::GetInstance()->GetPdgId(particle->GetPdg());
- if (ipid != -1)
- {
- h_minbias_Particle_PID_pt[ipid] ->Fill(particle->GetPt());
- h_minbias_Particle_PID_eta[ipid]->Fill(particle->GetEta());
- h_minbias_Particle_PID_E[ipid] ->Fill(particle->GetEnergy());
- h_minbias_Particle_PID_px[ipid] ->Fill(particle->GetPx());
- h_minbias_Particle_PID_py[ipid] ->Fill(particle->GetPy());
- h_minbias_Particle_PID_pz[ipid] ->Fill(particle->GetPz());
- h_minbias_Particle_PID_t[ipid] ->Fill(particle->GetT());
- h_minbias_Particle_PID_x[ipid] ->Fill(particle->GetX());
- h_minbias_Particle_PID_y[ipid] ->Fill(particle->GetY());
- h_minbias_Particle_PID_z[ipid] ->Fill(particle->GetZ());
- }
- }
- if (qaUtility::GetInstance()->Cut_Event_refmult(event) &&
- qaUtility::GetInstance()->Cut_Particle_refmult(particle))
- {
- Ncounter_refmult++;
- h_refmult_Particle_pt->Fill(particle->GetPt());
- h_refmult_Particle_eta->Fill(particle->GetEta());
- h_refmult_Particle_E->Fill(particle->GetEnergy());
- h_refmult_Particle_px->Fill(particle->GetPx());
- h_refmult_Particle_py->Fill(particle->GetPy());
- h_refmult_Particle_pz->Fill(particle->GetPz());
- h_refmult_Particle_t->Fill(particle->GetT());
- h_refmult_Particle_x->Fill(particle->GetX());
- h_refmult_Particle_y->Fill(particle->GetY());
- h_refmult_Particle_z->Fill(particle->GetZ());
- h_refmult_Particle_pdg->Fill(particle->GetPdg());
- if (qaUtility::GetInstance()->GetCharge(particle->GetPdg()) > 0)
- {
- ipid = 0;
- }
- if (qaUtility::GetInstance()->GetCharge(particle->GetPdg()) < 0)
- {
- ipid = 4;
- }
- if (ipid == 0 || ipid == 4)
- {
- h_refmult_Particle_PID_pt[ipid] ->Fill(particle->GetPt());
- h_refmult_Particle_PID_eta[ipid]->Fill(particle->GetEta());
- h_refmult_Particle_PID_E[ipid] ->Fill(particle->GetEnergy());
- h_refmult_Particle_PID_px[ipid] ->Fill(particle->GetPx());
- h_refmult_Particle_PID_py[ipid] ->Fill(particle->GetPy());
- h_refmult_Particle_PID_pz[ipid] ->Fill(particle->GetPz());
- h_refmult_Particle_PID_t[ipid] ->Fill(particle->GetT());
- h_refmult_Particle_PID_x[ipid] ->Fill(particle->GetX());
- h_refmult_Particle_PID_y[ipid] ->Fill(particle->GetY());
- h_refmult_Particle_PID_z[ipid] ->Fill(particle->GetZ());
- }
- ipid = qaUtility::GetInstance()->GetPdgId(particle->GetPdg());
- if (ipid >= 0)
- {
- h_refmult_Particle_PID_pt[ipid] ->Fill(particle->GetPt());
- h_refmult_Particle_PID_eta[ipid]->Fill(particle->GetEta());
- h_refmult_Particle_PID_E[ipid] ->Fill(particle->GetEnergy());
- h_refmult_Particle_PID_px[ipid] ->Fill(particle->GetPx());
- h_refmult_Particle_PID_py[ipid] ->Fill(particle->GetPy());
- h_refmult_Particle_PID_pz[ipid] ->Fill(particle->GetPz());
- h_refmult_Particle_PID_t[ipid] ->Fill(particle->GetT());
- h_refmult_Particle_PID_x[ipid] ->Fill(particle->GetX());
- h_refmult_Particle_PID_y[ipid] ->Fill(particle->GetY());
- h_refmult_Particle_PID_z[ipid] ->Fill(particle->GetZ());
- }
- }
- delete particle;
- }
- if (qaUtility::GetInstance()->Cut_Event_minbias(event))
- {
- Minbias_counter++;
- h_minbias_Event_b->Fill(event->GetB());
- h_minbias_Event_Mult->Fill(Ncounter_minbias);
- h2_minbias_Event_b_Mult->Fill(Ncounter_minbias, event->GetB());
- }
- else
- {
- vRejectedEvents.push_back(Absolute_counter-1);
- }
- if (qaUtility::GetInstance()->Cut_Event_refmult(event))
- {
- h_refmult_Event_b->Fill(event->GetB());
- h_refmult_Event_Mult->Fill(Ncounter_refmult);
- h2_refmult_Event_b_Mult->Fill(Ncounter_refmult, event->GetB());
- }
- delete event;
- }
- std::cout << "Loop is closed, " << Minbias_counter << " minbias events were counted (" << Absolute_counter-1 << " events in total)." << std::endl;
- if (qaUtility::GetInstance()->debug)
- {
- if (vRejectedEvents.size() > 0)
- {
- std::cout << vRejectedEvents.size() << " rejected events:" << std::endl;
- for (Long64_t ire=0; ire < (Long64_t)vRejectedEvents.size(); ire++)
- {
- std::cout << "Event " << vRejectedEvents.at(ire) << std::endl;
- }
- }
- else
- {
- std::cout << "No rejected events." << std::endl;
- }
- }
- fo->mkdir("minbias");
- fo->cd("minbias");
- h_minbias_Event_b->Write();
- h_minbias_Event_Mult->Write();
- h2_minbias_Event_b_Mult->Write();
- h_minbias_Particle_pt->Write();
- h_minbias_Particle_eta->Write();
- h_minbias_Particle_E->Write();
- h_minbias_Particle_px->Write();
- h_minbias_Particle_py->Write();
- h_minbias_Particle_pz->Write();
- h_minbias_Particle_t->Write();
- h_minbias_Particle_x->Write();
- h_minbias_Particle_y->Write();
- h_minbias_Particle_z->Write();
- h_minbias_Particle_pdg->Write();
- for (int i=0; i<qaUtility::GetInstance()->npid; i++)
- {
- h_minbias_Particle_PID_pt[i] ->Write();
- h_minbias_Particle_PID_eta[i]->Write();
- h_minbias_Particle_PID_E[i] ->Write();
- h_minbias_Particle_PID_px[i] ->Write();
- h_minbias_Particle_PID_py[i] ->Write();
- h_minbias_Particle_PID_pz[i] ->Write();
- h_minbias_Particle_PID_t[i] ->Write();
- h_minbias_Particle_PID_x[i] ->Write();
- h_minbias_Particle_PID_y[i] ->Write();
- h_minbias_Particle_PID_z[i] ->Write();
- }
- fo->mkdir("refmult");
- fo->cd("refmult");
- h_refmult_Event_b->Write();
- h_refmult_Event_Mult->Write();
- h2_refmult_Event_b_Mult->Write();
- h_refmult_Particle_pt->Write();
- h_refmult_Particle_eta->Write();
- h_refmult_Particle_E->Write();
- h_refmult_Particle_px->Write();
- h_refmult_Particle_py->Write();
- h_refmult_Particle_pz->Write();
- h_refmult_Particle_t->Write();
- h_refmult_Particle_x->Write();
- h_refmult_Particle_y->Write();
- h_refmult_Particle_z->Write();
- h_refmult_Particle_pdg->Write();
- for (int i=0; i<qaUtility::GetInstance()->npid; i++)
- {
- h_refmult_Particle_PID_pt[i] ->Write();
- h_refmult_Particle_PID_eta[i]->Write();
- h_refmult_Particle_PID_E[i] ->Write();
- h_refmult_Particle_PID_px[i] ->Write();
- h_refmult_Particle_PID_py[i] ->Write();
- h_refmult_Particle_PID_pz[i] ->Write();
- h_refmult_Particle_PID_t[i] ->Write();
- h_refmult_Particle_PID_x[i] ->Write();
- h_refmult_Particle_PID_y[i] ->Write();
- h_refmult_Particle_PID_z[i] ->Write();
- }
- fo->Close();
- timer.Stop();
- timer.Print();
- return 0;
- }
|