123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362 |
- #include "mfPicoDstReader.h"
- #include "mfNamespace.h"
- #include "mfProcessEvent.h"
- #include "mfInputReader.h"
- #include "mfSelector.h"
- #include "mfRecentering.h"
- #include "mfFlattening.h"
- #include "mfAcceptanceFilter.h"
- #include <TStopwatch.h>
- #include <TProfile.h>
- #include <TClonesArray.h>
- #include <TFile.h>
- #include <iostream>
- #include <TMath.h>
- #include <PicoDstMCEvent.h>
- #include <PicoDstRecoEvent.h>
- #include <PicoDstMCTrack.h>
- #include <PicoDstRecoTrack.h>
- #include <PicoDstFHCal.h>
- using mfdca::Ndim;
- using mfEp::EpNames;
- using mfEp::NetaGaps;
- using mfEp::Nflattening;
- using mfEp::NsubEvents;
- using mfEp::TpcEtaGap;
- using mfEp::TpcLEp;
- using mfEp::TpcREp;
- using mfEp::SpTpcLEp;
- using mfEp::SpTpcREp;
- using mfRes::NresBin;
- using TMath::ATan2;
- using TMath::Cos;
- using TMath::Pi;
- using TMath::Sin;
- int main(int argc, char **argv)
- {
- Bool_t isAcceptanceFilterOn = false;
- TString iFileName, oFileName, dcaFileName, recenteringFileName = "", flatteningFileName = "";
- if (argc < 7)
- {
- std::cerr << "./get-res -i INPUTFILE/LIST -o OUTPUT.root -dca DCA_FIT.root [OPTIONAL: -recentering RECENTERING.root -flattening FLATTENING.root --bad-acceptance]" << std::endl;
- return 1;
- }
- for (int i = 1; i < argc; i++)
- {
- if (std::string(argv[i]) != "-i" &&
- std::string(argv[i]) != "-o" &&
- std::string(argv[i]) != "-dca" &&
- std::string(argv[i]) != "-recentering" &&
- std::string(argv[i]) != "-flattening" &&
- std::string(argv[i]) != "--bad-acceptance")
- {
- 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]) == "-dca" && i != argc - 1)
- {
- dcaFileName = argv[++i];
- continue;
- }
- if (std::string(argv[i]) == "-dca" && i == argc - 1)
- {
- std::cerr << "\n[ERROR]: Input file name was not specified " << std::endl;
- return 1;
- }
- if (std::string(argv[i]) == "-recentering" && i != argc - 1)
- {
- recenteringFileName = argv[++i];
- continue;
- }
- if (std::string(argv[i]) == "-recentering" && i == argc - 1)
- {
- std::cerr << "\n[ERROR]: Input file name was not specified " << std::endl;
- return 1;
- }
- if (std::string(argv[i]) == "-flattening" && i != argc - 1)
- {
- flatteningFileName = argv[++i];
- continue;
- }
- if (std::string(argv[i]) == "-flattening" && i == argc - 1)
- {
- std::cerr << "\n[ERROR]: Input file name was not specified " << std::endl;
- return 1;
- }
- if (std::string(argv[i]) == "--bad-acceptance")
- {
- std::cout << "Message: Acceptance filter is ON." << std::endl;
- isAcceptanceFilterOn = true;
- continue;
- }
- }
- }
- TStopwatch timer;
- timer.Start();
- mfPicoDstReader *input = new mfPicoDstReader();
- TChain *inChain = (TChain *)input->Open(iFileName);
- if (!inChain)
- {
- std::cerr << "[ERROR]: nothing to read. Maybe the input TTree(s) is(are) empty?" << std::endl;
- return 10;
- }
- // Read recentering
- mfRecentering *recentering = new mfRecentering();
- recentering->ReadRecenteringFile(recenteringFileName);
- // Read flattening
- mfFlattening *flattening = new mfFlattening();
- flattening->ReadFlatteningFile(flatteningFileName);
- // Init output profiles
- TProfile *pResRecoSqrTpc[NetaGaps];
- TProfile *pResRecoSqrSpTpc[NetaGaps];
- TProfile *pResRpTpc[NetaGaps];
- TProfile *pResRpTpcR[NetaGaps];
- TProfile *pResRpTpcL[NetaGaps];
- TProfile *pResRecoSqrFHCal;
- TProfile *pResV1RpFHCal;
- TProfile *pResV1RpFHCalR;
- TProfile *pResV1RpFHCalL;
- TProfile *pResV2RpFHCal;
- TProfile *pResV2RpFHCalR;
- TProfile *pResV2RpFHCalL;
- TProfile *pResRecoSqrMcTpc[NetaGaps];
- TProfile *pResRpMcTpc[NetaGaps];
- TProfile *pResRpMcTpcR[NetaGaps];
- TProfile *pResRpMcTpcL[NetaGaps];
- for (int i = 0; i < NetaGaps; i++)
- {
- pResRecoSqrTpc[i] = new TProfile(Form("pResRecoSqrtTpc_gap%i", i),
- Form("pResRecoSqrtTpc_gap%i", i),
- NresBin, 0., 100.);
- pResRecoSqrSpTpc[i] = new TProfile(Form("pResRecoSqrtSpTpc_gap%i", i),
- Form("pResRecoSqrtSpTpc_gap%i", i),
- NresBin, 0., 100.);
- pResRpTpc[i] = new TProfile(Form("pResRpTpc_gap%i", i),
- Form("pResRpTpc_gap%i", i),
- NresBin, 0, 100);
- pResRpTpcR[i] = new TProfile(Form("pResRpTpcR_gap%i", i),
- Form("pResRpTpcR_gap%i", i),
- NresBin, 0, 100);
- pResRpTpcL[i] = new TProfile(Form("pResRpTpcL_gap%i", i),
- Form("pResRpTpcL_gap%i", i),
- NresBin, 0, 100);
- pResRecoSqrMcTpc[i] = new TProfile(Form("pResRecoSqrtMcTpc_gap%i", i),
- Form("pResRecoSqrtMcTpc_gap%i", i),
- NresBin, 0., 100.);
- pResRpMcTpc[i] = new TProfile(Form("pResRpMcTpc_gap%i", i),
- Form("pResRpMcTpc_gap%i", i),
- NresBin, 0, 100);
- pResRpMcTpcR[i] = new TProfile(Form("pResRpMcTpcR_gap%i", i),
- Form("pResRpMcTpcR_gap%i", i),
- NresBin, 0, 100);
- pResRpMcTpcL[i] = new TProfile(Form("pResRpMcTpcL_gap%i", i),
- Form("pResRpMcTpcL_gap%i", i),
- NresBin, 0, 100);
- }
- pResRecoSqrFHCal = new TProfile(Form("pResRecoSqrFHCal"),
- Form("pResRecoSqrFHCal"),
- NresBin, 0., 100.);
- pResV1RpFHCal = new TProfile(Form("pResV1RpFHCal"),
- Form("pResV1RpFHCal"),
- NresBin, 0., 100.);
- pResV2RpFHCal = new TProfile(Form("pResV2RpFHCal"),
- Form("pResV2RpFHCal"),
- NresBin, 0., 100.);
- pResV1RpFHCalR = new TProfile(Form("pResV1RpFHCalR"),
- Form("pResV1RpFHCalR"),
- NresBin, 0., 100.);
- pResV1RpFHCalL = new TProfile(Form("pResV1RpFHCalL"),
- Form("pResV1RpFHCalL"),
- NresBin, 0., 100.);
- pResV2RpFHCalR = new TProfile(Form("pResV2RpFHCalR"),
- Form("pResV2RpFHCalR"),
- NresBin, 0., 100.);
- pResV2RpFHCalL = new TProfile(Form("pResV2RpFHCalL"),
- Form("pResV2RpFHCalL"),
- NresBin, 0., 100.);
- PicoDstMCEvent *mcEvent = nullptr;
- // PicoDstRecoEvent *recoEvent = nullptr;
- TClonesArray *recoTracks = nullptr;
- TClonesArray *mcTracks = nullptr;
- TClonesArray *fhcalmodules = nullptr;
- inChain->SetBranchAddress("mcevent.", &mcEvent);
- // inChain->SetBranchAddress("recoevent.", &recoEvent);
- inChain->SetBranchAddress("mctracks", &mcTracks);
- inChain->SetBranchAddress("recotracks", &recoTracks);
- inChain->SetBranchAddress("FHCalModules", &fhcalmodules);
- // Start event loop
- int n_entries = inChain->GetEntries();
- mfSelector *selector = new mfSelector();
- mfInputReader *inputReader = new mfInputReader();
- inputReader->SetDcaFitFile(dcaFileName);
- TVector2 vQ_SP_L, vQ_SP_R;
- for (int iEv = 0; iEv < n_entries; iEv++)
- {
- if (iEv % 1000 == 0)
- std::cout << "Event [" << iEv << "/" << n_entries << "]" << std::endl;
- inChain->GetEntry(iEv);
- mfProcessEvent *procEvent = new mfProcessEvent();
- float cent = procEvent->GetCentB(mcEvent);
- if (cent == -1)
- continue;
- procEvent->SetInputReader(inputReader);
- mfAcceptanceFilter *AccFilter = new mfAcceptanceFilter();
- if (isAcceptanceFilterOn)
- AccFilter->Activate();
- procEvent->SetAcceptaceFilter(AccFilter);
- // std::cout << "\tNtracks = " << recoTracks->GetEntriesFast() << std::endl;
- mfEventInfo evInfoTpcRaw = procEvent->ProcessRecoTpcEP(recoTracks, 2.);
- if (evInfoTpcRaw.GetQvectors().size() == 0)
- continue;
- mfEventInfo evInfoTpcRec = recentering->ProcessRecentering(evInfoTpcRaw, cent);
- mfEventInfo evInfoTpcFlat = flattening->ProcessFlattening(evInfoTpcRec, cent);
- for (int iGap = 0; iGap < NetaGaps; iGap++)
- {
- // Ntracks >= 4 for TpcR and TpcL EP
- if (!selector->isGoodTpcEp(evInfoTpcRec.GetEpMult(TpcREp[iGap]), evInfoTpcRec.GetEpMult(TpcLEp[iGap])))
- continue;
- double PsiR = evInfoTpcFlat.GetQvector(TpcREp[iGap]).GetPhi();
- double PsiL = evInfoTpcFlat.GetQvector(TpcLEp[iGap]).GetPhi();
- double resReco = Cos(2. * (PsiR - PsiL));
- vQ_SP_L.Set(evInfoTpcFlat.GetQvector(SpTpcLEp[iGap]).X(), evInfoTpcFlat.GetQvector(SpTpcLEp[iGap]).Y());
- vQ_SP_R.Set(evInfoTpcFlat.GetQvector(SpTpcREp[iGap]).X(), evInfoTpcFlat.GetQvector(SpTpcREp[iGap]).Y());
- double resRecoSP = (double)(vQ_SP_R * vQ_SP_L);
- double resRpR = Cos(2. * (PsiR - mcEvent->GetPhiRP()));
- double resRpL = Cos(2. * (PsiL - mcEvent->GetPhiRP()));
- pResRecoSqrTpc[iGap]->Fill(cent, resReco);
- pResRecoSqrSpTpc[iGap]->Fill(cent, resRecoSP);
- pResRpTpcR[iGap]->Fill(cent, resRpR);
- pResRpTpcL[iGap]->Fill(cent, resRpL);
- pResRpTpc[iGap]->Fill(cent, resRpR);
- pResRpTpc[iGap]->Fill(cent, resRpL);
- }
- mfEventInfo evInfoFHCalRaw = procEvent->ProcessRecoFHCalEP(fhcalmodules, 1.);
- // FHCal energy > 0 in FHCal_L and FHCal_R
- if (!selector->isGoodFHCalEp(evInfoFHCalRaw.GetEpWeight(EpNames::kFHCal_R), evInfoFHCalRaw.GetEpWeight(EpNames::kFHCal_R)))
- continue;
- mfEventInfo evInfoFHCalRec = recentering->ProcessRecentering(evInfoFHCalRaw, cent);
- mfEventInfo evInfoFHCalFlat = flattening->ProcessFlattening(evInfoFHCalRec, cent);
- double PsiFHCalR = evInfoFHCalFlat.GetQvector(EpNames::kFHCal_R).GetPhi();
- double PsiFHCalL = evInfoFHCalFlat.GetQvector(EpNames::kFHCal_L).GetPhi();
- double PsiFHCal = evInfoFHCalFlat.GetQvector(EpNames::kFHCal).GetPhi();
- double resEpFHCal = Cos(1. * (PsiFHCalR - PsiFHCalL));
- double resV1RpFHCal = Cos(1. * (PsiFHCal - mcEvent->GetPhiRP()));
- double resV1RpFHCalR = Cos(1. * (PsiFHCalR - mcEvent->GetPhiRP()));
- double resV1RpFHCalL = Cos(1. * (PsiFHCalL - mcEvent->GetPhiRP()));
- double resV2RpFHCal = Cos(2. * (PsiFHCal - mcEvent->GetPhiRP()));
- double resV2RpFHCalR = Cos(2. * (PsiFHCalR - mcEvent->GetPhiRP()));
- double resV2RpFHCalL = Cos(2. * (PsiFHCalL - mcEvent->GetPhiRP()));
- pResRecoSqrFHCal->Fill(cent, resEpFHCal);
- pResV1RpFHCal->Fill(cent, resV1RpFHCal);
- pResV1RpFHCalR->Fill(cent, resV1RpFHCalR);
- pResV1RpFHCalL->Fill(cent, resV1RpFHCalL);
- pResV2RpFHCal->Fill(cent, resV2RpFHCal);
- pResV2RpFHCalR->Fill(cent, resV2RpFHCalR);
- pResV2RpFHCalL->Fill(cent, resV2RpFHCalL);
- mfEventInfo evInfoMcTpc = procEvent->ProcessMcTpcEP(mcTracks, 2.);
- if (evInfoMcTpc.GetQvectors().size() == 0)
- continue;
- for (int iGap = 0; iGap < NetaGaps; iGap++)
- {
- // Ntracks >= 4 for TpcR and TpcL EP
- if (!selector->isGoodTpcEp(evInfoMcTpc.GetEpMult(TpcREp[iGap]), evInfoMcTpc.GetEpMult(TpcLEp[iGap])))
- continue;
- double PsiR = evInfoMcTpc.GetQvector(TpcREp[iGap]).GetPhi();
- double PsiL = evInfoMcTpc.GetQvector(TpcLEp[iGap]).GetPhi();
- double resReco = Cos(2. * (PsiR - PsiL));
- double resRpR = Cos(2. * (PsiR - mcEvent->GetPhiRP()));
- double resRpL = Cos(2. * (PsiL - mcEvent->GetPhiRP()));
- pResRecoSqrMcTpc[iGap]->Fill(cent, resReco);
- pResRpMcTpcR[iGap]->Fill(cent, resRpR);
- pResRpMcTpcL[iGap]->Fill(cent, resRpL);
- pResRpMcTpc[iGap]->Fill(cent, resRpR);
- pResRpMcTpc[iGap]->Fill(cent, resRpL);
- }
- }
- // Init output file
- TFile *fo = new TFile(oFileName.Data(), "recreate");
- fo->cd();
- for (int iGap = 0; iGap < NetaGaps; iGap++)
- {
- pResRecoSqrTpc[iGap]->Write();
- pResRecoSqrSpTpc[iGap]->Write();
- pResRpTpcR[iGap]->Write();
- pResRpTpcL[iGap]->Write();
- pResRpTpc[iGap]->Write();
- }
- pResRecoSqrFHCal->Write();
- pResV1RpFHCal->Write();
- pResV1RpFHCalR->Write();
- pResV1RpFHCalL->Write();
- pResV2RpFHCal->Write();
- pResV2RpFHCalR->Write();
- pResV2RpFHCalL->Write();
- for (int iGap = 0; iGap < NetaGaps; iGap++)
- {
- pResRecoSqrMcTpc[iGap]->Write();
- pResRpMcTpcR[iGap]->Write();
- pResRpMcTpcL[iGap]->Write();
- pResRpMcTpc[iGap]->Write();
- }
- fo->Close();
- timer.Stop();
- timer.Print();
- return 0;
- }
|