123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347 |
- #include <map>
- #include <iostream>
- #include <TString.h>
- #include <TMath.h>
- #include "mfProcessEvent.h"
- #include "mfQvector.h"
- #include "mfSelector.h"
- using mfcent::CentB;
- using mfdca::NsigmaDcaCut;
- using mfEp::BadFHCalModules;
- using mfEp::EpNames;
- using mfEp::GetFHCalPhi;
- using mfEp::NetaGaps;
- using mfEp::Nmodules;
- using mfEp::NsubEvents;
- using mfEp::NsubEventsTpc;
- using mfEp::SpTpcLEp;
- using mfEp::SpTpcREp;
- using mfEp::TpcEtaGap;
- using mfEp::TpcLEp;
- using mfEp::TpcREp;
- using TMath::Abs;
- using TMath::Cos;
- using TMath::Sin;
- ClassImp(mfProcessEvent);
- mfProcessEvent::mfProcessEvent(/* args */)
- {
- }
- mfProcessEvent::~mfProcessEvent()
- {
- }
- Float_t mfProcessEvent::GetCentB(PicoDstMCEvent *const &mcEvent)
- {
- return CentB(mcEvent->GetB());
- }
- mfEventInfo mfProcessEvent::ProcessRecoTpcEP(TClonesArray *const &recoTracks, Double_t harmonic)
- {
- int mult = 0;
- std::vector<mfQvector> Qvectors;
- std::vector<double> Qx, Qy, Wt;
- std::vector<int> EpMult;
- mfSelector *selector = new mfSelector();
- mfEventInfo fEventInfo;
- if (!recoTracks)
- return fEventInfo;
- if (recoTracks->GetEntriesFast() == 0)
- return fEventInfo;
- if (!fInputReader)
- return fEventInfo;
- if (!fAccFilter)
- {
- fAccFilter = new mfAcceptanceFilter();
- }
- for (int i = 0; i < NsubEvents; i++)
- {
- Qx.push_back(0.);
- Qy.push_back(0.);
- Wt.push_back(0.);
- EpMult.push_back(0);
- Qvectors.push_back(mfQvector(harmonic, 0., 0.));
- }
- for (int iTrack = 0; iTrack < recoTracks->GetEntriesFast(); iTrack++)
- {
- auto recoTrack = (PicoDstRecoTrack *)recoTracks->UncheckedAt(iTrack);
- if (!recoTrack)
- continue;
- if (!fAccFilter->isGoodTpc(recoTrack->GetPhi()))
- continue;
- // Get refMult
- if (selector->isRefMult(recoTrack))
- {
- mult++;
- }
- // Primary track selection
- if (!selector->isRecoPrimary(recoTrack->GetDCAx(), fInputReader->GetDcaSigmaX(recoTrack)))
- continue;
- if (!selector->isRecoPrimary(recoTrack->GetDCAy(), fInputReader->GetDcaSigmaY(recoTrack)))
- continue;
- if (!selector->isRecoPrimary(recoTrack->GetDCAz(), fInputReader->GetDcaSigmaZ(recoTrack)))
- continue;
- // Good reco track selection
- if (!selector->isGoodRecoTrack(recoTrack))
- continue;
- for (int iGap = 0; iGap < NetaGaps; iGap++)
- {
- // Tpc left
- if (recoTrack->GetEta() < -1. * TpcEtaGap[iGap])
- {
- Qx.at(TpcLEp[iGap]) += recoTrack->GetPt() * Cos(harmonic * recoTrack->GetPhi());
- Qy.at(TpcLEp[iGap]) += recoTrack->GetPt() * Sin(harmonic * recoTrack->GetPhi());
- Wt.at(TpcLEp[iGap]) += recoTrack->GetPt();
- Qx.at(SpTpcLEp[iGap]) += 1. * Cos(harmonic * recoTrack->GetPhi());
- Qy.at(SpTpcLEp[iGap]) += 1. * Sin(harmonic * recoTrack->GetPhi());
- Wt.at(SpTpcLEp[iGap]) += 1.;
- EpMult.at(TpcLEp[iGap]) += 1;
- EpMult.at(SpTpcLEp[iGap]) += 1;
- }
- // Tpc right
- if (recoTrack->GetEta() > TpcEtaGap[iGap])
- {
- Qx.at(TpcREp[iGap]) += recoTrack->GetPt() * Cos(harmonic * recoTrack->GetPhi());
- Qy.at(TpcREp[iGap]) += recoTrack->GetPt() * Sin(harmonic * recoTrack->GetPhi());
- Wt.at(TpcREp[iGap]) += recoTrack->GetPt();
- Qx.at(SpTpcREp[iGap]) += 1. * Cos(harmonic * recoTrack->GetPhi());
- Qy.at(SpTpcREp[iGap]) += 1. * Sin(harmonic * recoTrack->GetPhi());
- Wt.at(SpTpcREp[iGap]) += 1.;
- EpMult.at(TpcREp[iGap]) += 1;
- EpMult.at(SpTpcREp[iGap]) += 1;
- }
- }
- }
- for (int i = EpNames::kTpc_R_EtaGap1; i <= EpNames::kTpc_L_EtaGap3; i++)
- {
- if (Wt.at(i) == 0 || EpMult.at(i) == 0)
- {
- return fEventInfo;
- }
- else
- {
- Qx.at(i) /= Wt.at(i);
- Qy.at(i) /= Wt.at(i);
- Qvectors.at(i).Set(Qx.at(i), Qy.at(i));
- }
- }
- for (int i = EpNames::kSpTpc_R_EtaGap1; i <= EpNames::kSpTpc_L_EtaGap3; i++)
- {
- if (Wt.at(i) == 0 || EpMult.at(i) == 0)
- {
- return fEventInfo;
- }
- else
- {
- Qvectors.at(i).Set(Qx.at(i), Qy.at(i));
- }
- }
- fEventInfo.SetMultiplicity(mult);
- fEventInfo.SetQv(Qvectors);
- fEventInfo.SetEpMult(EpMult);
- fEventInfo.SetEpWeight(Wt);
- return fEventInfo;
- }
- mfEventInfo mfProcessEvent::ProcessRecoFHCalEP(TClonesArray *const &fhcalModules, Double_t harmonic)
- {
- std::vector<mfQvector> Qvectors;
- std::vector<double> Qx, Qy, Wt;
- std::vector<int> EpMult;
- mfEventInfo fEventInfo;
- if (!fhcalModules)
- return fEventInfo;
- if (fhcalModules->GetEntriesFast() == 0)
- return fEventInfo;
- if (!fInputReader)
- return fEventInfo;
- for (int i = 0; i < NsubEvents; i++)
- {
- Qx.push_back(0.);
- Qy.push_back(0.);
- Wt.push_back(0.);
- EpMult.push_back(0);
- Qvectors.push_back(mfQvector(harmonic, 0., 0.));
- }
- for (int iMod = 0; iMod < fhcalModules->GetEntriesFast(); iMod++)
- {
- auto fhcalModule = (PicoDstFHCal *)fhcalModules->UncheckedAt(iMod);
- if (!fhcalModule)
- continue;
- if ((iMod == BadFHCalModules.at(0)) || (iMod == BadFHCalModules.at(1)))
- continue;
- if (!fAccFilter->isGoodFHCal(iMod))
- continue;
- double phi = GetFHCalPhi(iMod);
- int w_sign = (iMod < Nmodules) ? 1 : -1;
- // for left FHCal modules
- if (iMod < Nmodules)
- {
- Qx.at(EpNames::kFHCal_L) += w_sign * fhcalModule->GetEnergy() * Cos(harmonic * phi);
- Qy.at(EpNames::kFHCal_L) += w_sign * fhcalModule->GetEnergy() * Sin(harmonic * phi);
- Wt.at(EpNames::kFHCal_L) += fhcalModule->GetEnergy();
- }
- // for right FHCal modules
- if (iMod >= Nmodules)
- {
- Qx.at(EpNames::kFHCal_R) += w_sign * fhcalModule->GetEnergy() * Cos(harmonic * phi);
- Qy.at(EpNames::kFHCal_R) += w_sign * fhcalModule->GetEnergy() * Sin(harmonic * phi);
- Wt.at(EpNames::kFHCal_R) += fhcalModule->GetEnergy();
- }
- // Total FHCal modules
- Qx.at(EpNames::kFHCal) += w_sign * fhcalModule->GetEnergy() * Cos(harmonic * phi);
- Qy.at(EpNames::kFHCal) += w_sign * fhcalModule->GetEnergy() * Sin(harmonic * phi);
- Wt.at(EpNames::kFHCal) += fhcalModule->GetEnergy();
- }
- for (int i = EpNames::kFHCal_R; i <= EpNames::kFHCal; i++)
- {
- if (Wt.at(i) == 0)
- {
- // std::cerr << "No tracks for EP" << std::endl;
- return fEventInfo;
- }
- else
- {
- Qx.at(i) /= Wt.at(i);
- Qy.at(i) /= Wt.at(i);
- Qvectors.at(i).Set(Qx.at(i), Qy.at(i));
- }
- }
- fEventInfo.SetQv(Qvectors);
- fEventInfo.SetEpWeight(Wt);
- return fEventInfo;
- }
- mfEventInfo mfProcessEvent::ProcessMcTpcEP(TClonesArray *const &mcTracks, Double_t harmonic)
- {
- // int mult = 0;
- std::vector<mfQvector> Qvectors;
- std::vector<double> Qx, Qy, Wt;
- std::vector<int> EpMult;
- mfSelector *selector = new mfSelector();
- mfEventInfo fEventInfo;
- if (!mcTracks)
- return fEventInfo;
- if (mcTracks->GetEntriesFast() == 0)
- return fEventInfo;
- if (!fInputReader)
- return fEventInfo;
- for (int i = 0; i < NsubEvents; i++)
- {
- Qx.push_back(0.);
- Qy.push_back(0.);
- Wt.push_back(0.);
- EpMult.push_back(0);
- Qvectors.push_back(mfQvector(harmonic, 0., 0.));
- }
- for (int iTrack = 0; iTrack < mcTracks->GetEntriesFast(); iTrack++)
- {
- auto mcTrack = (PicoDstMCTrack *)mcTracks->UncheckedAt(iTrack);
- if (!mcTrack)
- continue;
- // Get refMult
- // if (selector->isRefMult(mcTrack))
- // {
- // mult++;
- // }
- // Primary track selection
- if (!selector->isMcPrimary(mcTrack->GetMotherId()))
- continue;
- // Good mc track selection
- if (!selector->isGoodMcTrack(mcTrack))
- continue;
- for (int iGap = 0; iGap < NetaGaps; iGap++)
- {
- // Tpc left
- if (mcTrack->GetEta() < -1. * TpcEtaGap[iGap])
- {
- Qx.at(TpcLEp[iGap]) += mcTrack->GetPt() * Cos(harmonic * mcTrack->GetPhi());
- Qy.at(TpcLEp[iGap]) += mcTrack->GetPt() * Sin(harmonic * mcTrack->GetPhi());
- Wt.at(TpcLEp[iGap]) += mcTrack->GetPt();
- Qx.at(SpTpcLEp[iGap]) += 1. * Cos(harmonic * mcTrack->GetPhi());
- Qy.at(SpTpcLEp[iGap]) += 1. * Sin(harmonic * mcTrack->GetPhi());
- Wt.at(TpcLEp[iGap]) += 1.;
- Wt.at(SpTpcLEp[iGap]) += 1.;
- EpMult.at(TpcLEp[iGap]) += 1;
- EpMult.at(SpTpcLEp[iGap]) += 1;
- }
- // Tpc right
- if (mcTrack->GetEta() > TpcEtaGap[iGap])
- {
- Qx.at(TpcREp[iGap]) += mcTrack->GetPt() * Cos(harmonic * mcTrack->GetPhi());
- Qy.at(TpcREp[iGap]) += mcTrack->GetPt() * Sin(harmonic * mcTrack->GetPhi());
- Wt.at(TpcREp[iGap]) += mcTrack->GetPt();
- Qx.at(SpTpcREp[iGap]) += 1. * Cos(harmonic * mcTrack->GetPhi());
- Qy.at(SpTpcREp[iGap]) += 1. * Sin(harmonic * mcTrack->GetPhi());
- Wt.at(SpTpcREp[iGap]) += 1.;
- Wt.at(SpTpcREp[iGap]) += 1.;
- EpMult.at(TpcREp[iGap]) += 1;
- EpMult.at(SpTpcREp[iGap]) += 1;
- }
- }
- }
- for (int i = EpNames::kTpc_R_EtaGap1; i <= EpNames::kTpc_L_EtaGap3; i++)
- {
- if (Wt.at(i) == 0 || EpMult.at(i) == 0)
- {
- return fEventInfo;
- }
- else
- {
- Qx.at(i) /= Wt.at(i);
- Qy.at(i) /= Wt.at(i);
- Qvectors.at(i).Set(Qx.at(i), Qy.at(i));
- }
- }
- for (int i = EpNames::kSpTpc_R_EtaGap1; i <= EpNames::kSpTpc_L_EtaGap3; i++)
- {
- if (Wt.at(i) == 0 || EpMult.at(i) == 0)
- {
- return fEventInfo;
- }
- else
- {
- Qvectors.at(i).Set(Qx.at(i), Qy.at(i));
- }
- }
- // fEventInfo.SetMultiplicity(mult);
- fEventInfo.SetQv(Qvectors);
- fEventInfo.SetEpMult(EpMult);
- fEventInfo.SetEpWeight(Wt);
- return fEventInfo;
- }
|