diff --git a/PWGCF/DataModel/FemtoDerived.h b/PWGCF/DataModel/FemtoDerived.h index be637c38d8c..13e3e65b3bb 100644 --- a/PWGCF/DataModel/FemtoDerived.h +++ b/PWGCF/DataModel/FemtoDerived.h @@ -37,7 +37,8 @@ enum CollisionBinning { kMult, //! Bin collision in number of charged tracks for mixing kMultPercentile, //! Bin collision in multiplicity percentile for mixing kMultMultPercentile, //! Bin collision in number of charged tracks and multiplicity percentile for mixing - kMultPercentileQn, //! Bin collision in multiplicity percentile an qn value for mixing + kMultPercentileQn, //! Bin collision in multiplicity percentile and qn value for mixing + kMultPercentileEP, //! Bin collision in multiplicity percentile and event plane (deg) for mixing kNCollisionBinning }; @@ -54,8 +55,9 @@ DECLARE_SOA_COLUMN(BitMaskTrackThree, bitmaskTrackThree, BitMaskType); //! Bit f DECLARE_SOA_COLUMN(Downsample, downsample, bool); //! Flag for downsampling -DECLARE_SOA_COLUMN(QnVal, qnVal, int); //! qn values for dividing events -DECLARE_SOA_COLUMN(Occupancy, occupancy, int); //! Occupancy of the event +DECLARE_SOA_COLUMN(QnVal, qnVal, double); //! qn values for dividing events +DECLARE_SOA_COLUMN(Occupancy, occupancy, int); //! Occupancy of the event +DECLARE_SOA_COLUMN(EventPlane, eventPlane, double); //! Event-plane of the event (deg) } // namespace femtodreamcollision DECLARE_SOA_TABLE_STAGED(FDCollisions, "FDCOLLISION", @@ -71,6 +73,9 @@ DECLARE_SOA_TABLE(FDExtQnCollisions, "AOD", "FDEXTQNCOLLISION", femtodreamcollision::QnVal, femtodreamcollision::Occupancy); +DECLARE_SOA_TABLE(FDExtEPCollisions, "AOD", "FDEXTEPCOLLISION", + femtodreamcollision::EventPlane); + DECLARE_SOA_TABLE(FDColMasks, "AOD", "FDCOLMASK", femtodreamcollision::BitMaskTrackOne, femtodreamcollision::BitMaskTrackTwo, @@ -228,47 +233,47 @@ enum CharmHadronMassHypo { lcToPiKP = 2, dplusToPiKPi = 4 }; -DECLARE_SOA_COLUMN(GIndexCol, gIndexCol, int); //! Global index for the collision -DECLARE_SOA_COLUMN(TimeStamp, timeStamp, int64_t); //! Timestamp for the collision -DECLARE_SOA_COLUMN(VertexZ, vertexZ, float); //! VertexZ for the collision -DECLARE_SOA_COLUMN(TrackId, trackId, int); //! track id to match associate particle with charm hadron prongs -DECLARE_SOA_COLUMN(Charge, charge, int8_t); //! Charge of charm hadron -DECLARE_SOA_COLUMN(Prong0Id, prong0Id, int); //! Track id of charm hadron prong0 -DECLARE_SOA_COLUMN(Prong1Id, prong1Id, int); //! Track id of charm hadron prong1 -DECLARE_SOA_COLUMN(Prong2Id, prong2Id, int); //! Track id of charm hadron prong2 -DECLARE_SOA_COLUMN(Prong0Pt, prong0Pt, float); //! Track pT of charm hadron prong0 -DECLARE_SOA_COLUMN(Prong1Pt, prong1Pt, float); //! Track pT of charm hadron prong1 -DECLARE_SOA_COLUMN(Prong2Pt, prong2Pt, float); //! Track pT of charm hadron prong2 -DECLARE_SOA_COLUMN(Prong0Eta, prong0Eta, float); //! Track eta of charm hadron prong0 -DECLARE_SOA_COLUMN(Prong1Eta, prong1Eta, float); //! Track eta of charm hadron prong1 -DECLARE_SOA_COLUMN(Prong2Eta, prong2Eta, float); //! Track eta of charm hadron prong2 -DECLARE_SOA_COLUMN(Prong0Phi, prong0Phi, float); //! Track phi of charm hadron prong0 -DECLARE_SOA_COLUMN(Prong1Phi, prong1Phi, float); //! Track phi of charm hadron prong1 -DECLARE_SOA_COLUMN(Prong2Phi, prong2Phi, float); //! Track phi of charm hadron prong2 -DECLARE_SOA_COLUMN(CandidateSelFlag, candidateSelFlag, int); //! Selection of mass hypothesis for charm hadron (1 for Lc -> pkpi, 2 for Lc -> pikp, 4 for D+ -> pikpi) -DECLARE_SOA_COLUMN(BDTBkg, bdtBkg, float); //! Background score using Boosted Decision Tree for charm hadron -DECLARE_SOA_COLUMN(BDTPrompt, bdtPrompt, float); //! Prompt signal score using Boosted Decision Tree for charm hadron -DECLARE_SOA_COLUMN(BDTFD, bdtFD, float); //! Feed-down score using Boosted Decision Tree for charm hadron -DECLARE_SOA_COLUMN(FlagMc, flagMc, int); //! To select MC particle among charm hadrons, { DplusToPiKPi = 1, LcToPKPi = 17, DsToKKPi = 6, XicToPKPi = 21, N3ProngD = 2ecays }; -DECLARE_SOA_COLUMN(OriginMcRec, originMcRec, int); //! flag for reconstruction level matching (1 for prompt, 2 for non-prompt) -DECLARE_SOA_COLUMN(OriginMcGen, originMcGen, int); //! flag for generator level matching (1 for prompt, 2 for non-prompt) -DECLARE_SOA_COLUMN(IsCandidateSwapped, isCandidateSwapped, int); //! swapping of the prongs order (0 for Lc -> pkpi, 1 for Lc -> pikp) -DECLARE_SOA_COLUMN(TrkPt, trkPt, float); //! Transverse momentum of associate femto particle -DECLARE_SOA_COLUMN(TrkEta, trkEta, float); //! Eta of associate femto particle -DECLARE_SOA_COLUMN(TrkPhi, trkPhi, float); //! Phi of associate femto particle -DECLARE_SOA_COLUMN(Kstar, kstar, float); //! Relative momentum in particles pair frame -DECLARE_SOA_COLUMN(KT, kT, float); //! kT distribution of particle pairs -DECLARE_SOA_COLUMN(MT, mT, float); //! Transverse mass distribution -DECLARE_SOA_COLUMN(CharmM, charmM, float); //! Charm hadron mass -DECLARE_SOA_COLUMN(CharmTrkM, charmtrkM, float); //! Charm hadron track mass -DECLARE_SOA_COLUMN(CharmPt, charmPt, float); //! Transverse momentum of charm hadron for result task -DECLARE_SOA_COLUMN(CharmEta, charmEta, float); //! Eta of charm hadron for result task -DECLARE_SOA_COLUMN(CharmPhi, charmPhi, float); //! Phi of charm hadron for result task -DECLARE_SOA_COLUMN(Mult, mult, int); //! Charge particle multiplicity -DECLARE_SOA_COLUMN(MultPercentile, multPercentile, float); //! Multiplicity precentile -DECLARE_SOA_COLUMN(PairSign, pairSign, int8_t); //! Selection between like sign (1) and unlike sign pair (2) -DECLARE_SOA_COLUMN(ProcessType, processType, int64_t); //! Selection between same-event (1), and mixed-event (2) -DECLARE_SOA_DYNAMIC_COLUMN(M, m, //! +DECLARE_SOA_COLUMN(GIndexCol, gIndexCol, int); //! Global index for the collision +DECLARE_SOA_COLUMN(TimeStamp, timeStamp, int64_t); //! Timestamp for the collision +DECLARE_SOA_COLUMN(VertexZ, vertexZ, float); //! VertexZ for the collision +DECLARE_SOA_COLUMN(TrackId, trackId, int); //! track id to match associate particle with charm hadron prongs +DECLARE_SOA_COLUMN(Charge, charge, int8_t); //! Charge of charm hadron +DECLARE_SOA_COLUMN(Prong0Id, prong0Id, int); //! Track id of charm hadron prong0 +DECLARE_SOA_COLUMN(Prong1Id, prong1Id, int); //! Track id of charm hadron prong1 +DECLARE_SOA_COLUMN(Prong2Id, prong2Id, int); //! Track id of charm hadron prong2 +DECLARE_SOA_COLUMN(Prong0Pt, prong0Pt, float); //! Track pT of charm hadron prong0 +DECLARE_SOA_COLUMN(Prong1Pt, prong1Pt, float); //! Track pT of charm hadron prong1 +DECLARE_SOA_COLUMN(Prong2Pt, prong2Pt, float); //! Track pT of charm hadron prong2 +DECLARE_SOA_COLUMN(Prong0Eta, prong0Eta, float); //! Track eta of charm hadron prong0 +DECLARE_SOA_COLUMN(Prong1Eta, prong1Eta, float); //! Track eta of charm hadron prong1 +DECLARE_SOA_COLUMN(Prong2Eta, prong2Eta, float); //! Track eta of charm hadron prong2 +DECLARE_SOA_COLUMN(Prong0Phi, prong0Phi, float); //! Track phi of charm hadron prong0 +DECLARE_SOA_COLUMN(Prong1Phi, prong1Phi, float); //! Track phi of charm hadron prong1 +DECLARE_SOA_COLUMN(Prong2Phi, prong2Phi, float); //! Track phi of charm hadron prong2 +DECLARE_SOA_COLUMN(CandidateSelFlag, candidateSelFlag, int); //! Selection of mass hypothesis for charm hadron (1 for Lc -> pkpi, 2 for Lc -> pikp, 4 for D+ -> pikpi) +DECLARE_SOA_COLUMN(BDTBkg, bdtBkg, float); //! Background score using Boosted Decision Tree for charm hadron +DECLARE_SOA_COLUMN(BDTPrompt, bdtPrompt, float); //! Prompt signal score using Boosted Decision Tree for charm hadron +DECLARE_SOA_COLUMN(BDTFD, bdtFD, float); //! Feed-down score using Boosted Decision Tree for charm hadron +DECLARE_SOA_COLUMN(FlagMc, flagMc, int); //! To select MC particle among charm hadrons, { DplusToPiKPi = 1, LcToPKPi = 17, DsToKKPi = 6, XicToPKPi = 21, N3ProngD = 2ecays }; +DECLARE_SOA_COLUMN(OriginMcRec, originMcRec, int); //! flag for reconstruction level matching (1 for prompt, 2 for non-prompt) +DECLARE_SOA_COLUMN(OriginMcGen, originMcGen, int); //! flag for generator level matching (1 for prompt, 2 for non-prompt) +DECLARE_SOA_COLUMN(IsCandidateSwapped, isCandidateSwapped, int); //! swapping of the prongs order (0 for Lc -> pkpi, 1 for Lc -> pikp) +DECLARE_SOA_COLUMN(TrkPt, trkPt, float); //! Transverse momentum of associate femto particle +DECLARE_SOA_COLUMN(TrkEta, trkEta, float); //! Eta of associate femto particle +DECLARE_SOA_COLUMN(TrkPhi, trkPhi, float); //! Phi of associate femto particle +DECLARE_SOA_COLUMN(Kstar, kstar, float); //! Relative momentum in particles pair frame +DECLARE_SOA_COLUMN(KT, kT, float); //! kT distribution of particle pairs +DECLARE_SOA_COLUMN(MT, mT, float); //! Transverse mass distribution +DECLARE_SOA_COLUMN(CharmM, charmM, float); //! Charm hadron mass +DECLARE_SOA_COLUMN(CharmTrkM, charmtrkM, float); //! Charm hadron track mass +DECLARE_SOA_COLUMN(CharmPt, charmPt, float); //! Transverse momentum of charm hadron for result task +DECLARE_SOA_COLUMN(CharmEta, charmEta, float); //! Eta of charm hadron for result task +DECLARE_SOA_COLUMN(CharmPhi, charmPhi, float); //! Phi of charm hadron for result task +DECLARE_SOA_COLUMN(Mult, mult, int); //! Charge particle multiplicity +DECLARE_SOA_COLUMN(MultPercentile, multPercentile, float); //! Multiplicity precentile +DECLARE_SOA_COLUMN(PairSign, pairSign, int8_t); //! Selection between like sign (1) and unlike sign pair (2) +DECLARE_SOA_COLUMN(ProcessType, processType, int64_t); //! Selection between same-event (1), and mixed-event (2) +DECLARE_SOA_DYNAMIC_COLUMN(M, m, //! [](float pt0, float phi0, float eta0, float pt1, float phi1, float eta1, float pt2, float phi2, float eta2, const std::array& m) -> float { return RecoDecay::m(std::array{ RecoDecayPtEtaPhi::pVector(pt0, eta0, phi0), RecoDecayPtEtaPhi::pVector(pt1, eta1, phi1), diff --git a/PWGCF/FemtoDream/Core/femtoDreamCollisionSelection.h b/PWGCF/FemtoDream/Core/femtoDreamCollisionSelection.h index c8a3947b141..4e67c43b40c 100644 --- a/PWGCF/FemtoDream/Core/femtoDreamCollisionSelection.h +++ b/PWGCF/FemtoDream/Core/femtoDreamCollisionSelection.h @@ -24,6 +24,8 @@ #include "Framework/HistogramRegistry.h" #include "Framework/Logger.h" +#include "TMath.h" + #include #include #include @@ -212,18 +214,16 @@ class FemtoDreamCollisionSelection /// Initializes histograms for qn bin /// \param registry Histogram registry to be passed - void initQn(HistogramRegistry* registry, int mumQnBins = 10) + void initEPQA(HistogramRegistry* registry) { mHistogramQn = registry; - mHistogramQn->add("Event/centFT0CBefore", "; cent", kTH1F, {{10, 0, 100}}); - mHistogramQn->add("Event/centFT0CAfter", "; cent", kTH1F, {{10, 0, 100}}); + mHistogramQn->add("Event/centFT0CBeforeQn", "; cent", kTH1F, {{10, 0, 100}}); + mHistogramQn->add("Event/centFT0CAfterQn", "; cent", kTH1F, {{10, 0, 100}}); mHistogramQn->add("Event/centVsqn", "; cent; qn", kTH2F, {{10, 0, 100}, {100, 0, 1000}}); mHistogramQn->add("Event/centVsqnVsSpher", "; cent; qn; Sphericity", kTH3F, {{10, 0, 100}, {100, 0, 1000}, {100, 0, 1}}); mHistogramQn->add("Event/qnBin", "; qnBin; entries", kTH1F, {{20, 0, 20}}); + mHistogramQn->add("Event/psiEP", "; #Psi_{EP} (deg); entries", kTH1F, {{100, 0, 180}}); - for (int iqn(0); iqn < mumQnBins; ++iqn) { - qnMults.push_back(mHistogramQn->add(("Qn/mult_" + std::to_string(iqn)).c_str(), "; cent; c22", kTH1F, {{100, 0, 3500}})); - } return; } @@ -242,12 +242,12 @@ class FemtoDreamCollisionSelection mMQWeightthisEvt = new TH2D("MQWeightthisEvt", "", binPt, 0., 5., binEta, -0.8, 0.8); mHistogramQn = registry; - mHistogramQn->add("Event/profileC22", "; cent; c22", kTProfile, {{10, 0, 100}}); - mHistogramQn->add("Event/profileC24", "; cent; c24", kTProfile, {{10, 0, 100}}); + mHistogramQn->add("Event/profileC22", "; cent; c22", kTProfile, {{10, 0, 100}}, "s"); + mHistogramQn->add("Event/profileC24", "; cent; c24", kTProfile, {{10, 0, 100}}, "s"); if (doQnSeparation) { for (int iqn(0); iqn < mumQnBins; ++iqn) { - profilesC22.push_back(mHistogramQn->add(("Qn/profileC22_" + std::to_string(iqn)).c_str(), "; cent; c22", kTProfile, {{10, 0, 100}})); + profilesC22.push_back(mHistogramQn->add(("Qn/profileC22_" + std::to_string(iqn)).c_str(), "; cent; c22", kTProfile, {{10, 0, 100}}, "s")); } } return; @@ -324,7 +324,6 @@ class FemtoDreamCollisionSelection return spt; } - /// \todo to be implemented! /// Compute the qn-vector(FT0C) of an event /// \tparam T type of the collision /// \param col Collision @@ -336,7 +335,21 @@ class FemtoDreamCollisionSelection return qn; } - /// \todo to be implemented! + /// Compute the event plane of an event + /// \tparam T type of the collision + /// \param col Collision + /// \param nmode EP in which harmonic(default 2nd harmonic) + /// \return angle of the event plane (rad) of FT0C of the event + template + float computeEP(T const& col, int nmode) + { + double EP = ((1. / nmode) * (TMath::ATan2(col.qvecFT0CImVec()[0], col.qvecFT0CReVec()[0]))); + if (EP < 0) + EP += TMath::Pi(); + // atan2 return in rad -pi/2-pi/2, then make it 0-pi + return EP; + } + /// \return the 1-d qn-vector separator to 2-d std::vector> getQnBinSeparator2D(std::vector flat, const int numQnBins = 10) { @@ -359,11 +372,11 @@ class FemtoDreamCollisionSelection return res; } - /// \todo to be implemented! /// Get the bin number of qn-vector(FT0C) of an event /// \param centBinWidth centrality bin width, example: per 1%, per 10% ... /// \return bin number of qn-vector of the event - int myqnBin(float centrality, float centMax, std::vector qnBinSeparator, bool doFillHisto, float fSpher, float qn, const int numQnBins, float mult, float centBinWidth = 1.f) + // add a param : bool doFillHisto ? + int myqnBin(float centrality, float centMax, bool doFillCent, std::vector qnBinSeparator, float qn, const int numQnBins, float centBinWidth = 1.f) { auto twoDSeparator = getQnBinSeparator2D(qnBinSeparator, numQnBins); if (twoDSeparator.empty() || twoDSeparator[0][0] == -999.) { @@ -371,9 +384,9 @@ class FemtoDreamCollisionSelection return -999; // safe fallback } - if (doFillHisto) - mHistogramQn->fill(HIST("Event/centFT0CBefore"), centrality); - + // if (doFillHisto) + // mHistogramQn->fill(HIST("Event/centFT0CBefore"), centrality); + // add a param : bool doFillHisto ? int qnBin = -999; int mycentBin = static_cast(centrality / centBinWidth); if (mycentBin >= static_cast(centMax / centBinWidth)) @@ -382,6 +395,9 @@ class FemtoDreamCollisionSelection if (mycentBin > static_cast(twoDSeparator.size()) - 1) return qnBin; + if (doFillCent) + mHistogramQn->fill(HIST("Event/centFT0CAfterQn"), centrality); + for (int iqn(0); iqn < static_cast(twoDSeparator[mycentBin].size()) - 1; ++iqn) { if (qn > twoDSeparator[mycentBin][iqn] && qn <= twoDSeparator[mycentBin][iqn + 1]) { qnBin = iqn; @@ -392,20 +408,19 @@ class FemtoDreamCollisionSelection } mQnBin = qnBin; - - if (doFillHisto) { - mHistogramQn->fill(HIST("Event/centFT0CAfter"), centrality); - mHistogramQn->fill(HIST("Event/centVsqn"), centrality, qn); - mHistogramQn->fill(HIST("Event/centVsqnVsSpher"), centrality, qn, fSpher); - mHistogramQn->fill(HIST("Event/qnBin"), qnBin); - if (qnBin >= 0 && qnBin < numQnBins) { - std::get>(qnMults[qnBin])->Fill(mult); - } - } - return qnBin; } + /// \fill event-wise informations + void fillEPQA(float centrality, float fSpher, float qn, float psiEP) + { + mHistogramQn->fill(HIST("Event/centFT0CBeforeQn"), centrality); + mHistogramQn->fill(HIST("Event/centVsqn"), centrality, qn); + mHistogramQn->fill(HIST("Event/centVsqnVsSpher"), centrality, qn, fSpher); + mHistogramQn->fill(HIST("Event/qnBin"), mQnBin + 0.f); + mHistogramQn->fill(HIST("Event/psiEP"), psiEP); + } + /// \todo to be implemented! /// Fill cumulants histo for flow calculation /// Reset hists event-by-event @@ -501,7 +516,6 @@ class FemtoDreamCollisionSelection float mSphericityPtmin = 0.f; int mQnBin = -999; HistogramRegistry* mHistogramQn = nullptr; ///< For flow cumulant output - std::vector qnMults; /// Histograms of multiplicity (TH1F) per Qn bin. Stored as HistPtr (variant of shared_ptr) from HistogramManager. std::vector profilesC22; /// Pofile Histograms of c22 per Qn bin TH2D* mReQthisEvt = nullptr; ///< For flow cumulant in an event TH2D* mImQthisEvt = nullptr; ///< For flow cumulant in an event diff --git a/PWGCF/FemtoDream/Core/femtoDreamContainer.h b/PWGCF/FemtoDream/Core/femtoDreamContainer.h index 2b58f0e9c5c..ca8d99492bf 100644 --- a/PWGCF/FemtoDream/Core/femtoDreamContainer.h +++ b/PWGCF/FemtoDream/Core/femtoDreamContainer.h @@ -183,18 +183,18 @@ class FemtoDreamContainer } } - /// Initialize the histograms for pairs in divided qn bins + /// Initialize the histograms for pairs in divided qn or phi-psi bins template - void init_base_qn(std::string folderName, std::string femtoObs, - T& femtoObsAxis, T& mTAxi4D, T& multPercentileAxis4D, T& qnAxis4D) + void init_base_EP(std::string folderName, std::string femtoObs, + T& femtoObsAxis, T& mTAxi4D, T& multPercentileAxis4D, T& epObsAxis, std::string epObs) { - mHistogramRegistry->add((folderName + "/relPairkstarmTMultMultPercentileQn").c_str(), ("; " + femtoObs + "; #it{m}_{T} (GeV/#it{c}); Centrality; qn").c_str(), kTHnSparseF, {femtoObsAxis, mTAxi4D, multPercentileAxis4D, qnAxis4D}); + mHistogramRegistry->add((folderName + "/relPairkstarmTMultMultPercentileQn").c_str(), ("; " + femtoObs + "; #it{m}_{T} (GeV/#it{c}); Centrality;" + epObs).c_str(), kTHnSparseF, {femtoObsAxis, mTAxi4D, multPercentileAxis4D, epObsAxis}); } template - void init_qn(HistogramRegistry* registry, - T& kstarBins4D, T& mTBins4D, T& multPercentileBins4D, - bool isMC, float highkstarCut, ConfigurableAxis qnBins4D = {"qnBins4D", {10, 0, 10}, "qn binning"}) + void init_EP(HistogramRegistry* registry, + T& kstarBins4D, T& mTBins4D, T& multPercentileBins4D, T& epObsBins, std::string epObs, + bool isMC, float highkstarCut) { mHistogramRegistry = registry; std::string femtoObs; @@ -206,17 +206,56 @@ class FemtoDreamContainer framework::AxisSpec kstarAxis4D = {kstarBins4D, femtoObs}; framework::AxisSpec mTAxis4D = {mTBins4D, "#it{m}_{T} (GeV/#it{c})"}; framework::AxisSpec multPercentileAxis4D = {multPercentileBins4D, "Centralty(%)"}; - framework::AxisSpec qnAxis4D = {qnBins4D, "qn"}; + framework::AxisSpec epObsAxis = {epObsBins, epObs}; - std::string folderName = static_cast(mFolderSuffix[mEventType]) + static_cast(o2::aod::femtodreamMCparticle::MCTypeName[o2::aod::femtodreamMCparticle::MCType::kRecon]) + static_cast("_qn"); + std::string folderName = static_cast(mFolderSuffix[mEventType]) + static_cast(o2::aod::femtodreamMCparticle::MCTypeName[o2::aod::femtodreamMCparticle::MCType::kRecon]) + static_cast("_EP"); - init_base_qn(folderName, femtoObs, - kstarAxis4D, mTAxis4D, multPercentileAxis4D, qnAxis4D); + init_base_EP(folderName, femtoObs, + kstarAxis4D, mTAxis4D, multPercentileAxis4D, epObsAxis, epObs); if (isMC) { - folderName = static_cast(mFolderSuffix[mEventType]) + static_cast(o2::aod::femtodreamMCparticle::MCTypeName[o2::aod::femtodreamMCparticle::MCType::kTruth]) + static_cast("_qn"); - init_base_qn(folderName, femtoObs, - kstarAxis4D, mTAxis4D, multPercentileAxis4D, qnAxis4D); + folderName = static_cast(mFolderSuffix[mEventType]) + static_cast(o2::aod::femtodreamMCparticle::MCTypeName[o2::aod::femtodreamMCparticle::MCType::kTruth]) + static_cast("_EP"); + init_base_EP(folderName, femtoObs, + kstarAxis4D, mTAxis4D, multPercentileAxis4D, epObsAxis, epObs); + } + } + + /// Initialize the histograms for pairs with 3D component in divided qn bins + template + void init_base_3Dqn(std::string folderName, std::string femtoDKout, std::string femtoDKside, std::string femtoDKlong, + T& femtoDKoutAxis, T& femtoDKsideAxis, T& femtoDKlongAxis, T& mTAxi4D, T& multPercentileAxis4D, T& qnAxis, T& pairPhiAxis) + { + mHistogramRegistry->add((folderName + "/relPair3dRmTMultPercentileQnPairphi").c_str(), ("; " + femtoDKout + femtoDKside + femtoDKlong + "; #it{m}_{T} (GeV/#it{c}); Centrality; qn; #varphi_{pair} - #Psi_{EP}").c_str(), kTHnSparseF, {femtoDKoutAxis, femtoDKsideAxis, femtoDKlongAxis, mTAxi4D, multPercentileAxis4D, qnAxis, pairPhiAxis}); + } + + template + void init_3Dqn(HistogramRegistry* registry, + T& DKoutBins, T& DKsideBins, T& DKlongBins, T& mTBins4D, T& multPercentileBins4D, + bool isMC, ConfigurableAxis qnBins = {"qnBins", {10, 0, 10}, "qn binning"}, ConfigurableAxis pairPhiBins = {"phiBins", {10, 0 - 0.05, TMath::Pi() + 0.05}, "pair phi binning"}) + { + mHistogramRegistry = registry; + + std::string femtoObsDKout = "DK_{out} (GeV/#it{c})"; + std::string femtoObsDKside = "DK_{side} (GeV/#it{c})"; + std::string femtoObsDKlong = "DK_{long} (GeV/#it{c})"; + + framework::AxisSpec DKoutAxis = {DKoutBins, femtoObsDKout}; + framework::AxisSpec DKsideAxis = {DKsideBins, femtoObsDKside}; + framework::AxisSpec DKlongAxis = {DKlongBins, femtoObsDKlong}; + framework::AxisSpec mTAxis4D = {mTBins4D, "#it{m}_{T} (GeV/#it{c})"}; + framework::AxisSpec multPercentileAxis4D = {multPercentileBins4D, "Centralty(%)"}; + framework::AxisSpec qnAxis = {qnBins, "qn"}; + framework::AxisSpec pairPhiAxis = {pairPhiBins, "#varphi_{pair} - #Psi_{EP} (rad)"}; + + std::string folderName = static_cast(mFolderSuffix[mEventType]) + static_cast(o2::aod::femtodreamMCparticle::MCTypeName[o2::aod::femtodreamMCparticle::MCType::kRecon]) + static_cast("_3Dqn"); + + init_base_3Dqn(folderName, femtoObsDKout, femtoObsDKside, femtoObsDKlong, + DKoutAxis, DKsideAxis, DKlongAxis, mTAxis4D, multPercentileAxis4D, qnAxis, pairPhiAxis); + + if (isMC) { + folderName = static_cast(mFolderSuffix[mEventType]) + static_cast(o2::aod::femtodreamMCparticle::MCTypeName[o2::aod::femtodreamMCparticle::MCType::kTruth]) + static_cast("_3Dqn"); + init_base_3Dqn(folderName, femtoObsDKout, femtoObsDKside, femtoObsDKlong, + DKoutAxis, DKsideAxis, DKlongAxis, mTAxis4D, multPercentileAxis4D, qnAxis, pairPhiAxis); } } @@ -332,7 +371,7 @@ class FemtoDreamContainer } const float mTMC = FemtoDreamMath::getmT(part1.fdMCParticle(), mMassOne, part2, mMassTwo); - if (abs(part1.fdMCParticle().pdgMCTruth()) == mPDGOne) { // Note: all pair-histogramms are filled with MC truth information ONLY in case of non-fake candidates + if (std::abs(part1.fdMCParticle().pdgMCTruth()) == mPDGOne) { // Note: all pair-histogramms are filled with MC truth information ONLY in case of non-fake candidates setPair_base(femtoObsMC, mTMC, part1.fdMCParticle(), part2, mult, multPercentile, use4dplots, extendedplots); setPair_MC(femtoObsMC, femtoObs, mT, mult, part1.fdMCParticle().partOriginMCTruth(), part2.flagMc(), smearingByOrigin); } else { @@ -346,7 +385,7 @@ class FemtoDreamContainer } const float mTMC = FemtoDreamMath::getmT(part1.fdMCParticle(), mMassOne, part2.fdMCParticle(), mMassTwo); - if (abs(part1.fdMCParticle().pdgMCTruth()) == mPDGOne && abs(part2.fdMCParticle().pdgMCTruth()) == mPDGTwo) { // Note: all pair-histogramms are filled with MC truth information ONLY in case of non-fake candidates + if (std::abs(part1.fdMCParticle().pdgMCTruth()) == mPDGOne && std::abs(part2.fdMCParticle().pdgMCTruth()) == mPDGTwo) { // Note: all pair-histogramms are filled with MC truth information ONLY in case of non-fake candidates setPair_base(femtoObsMC, mTMC, part1.fdMCParticle(), part2.fdMCParticle(), mult, multPercentile, use4dplots, extendedplots); setPair_MC(femtoObsMC, femtoObs, mT, mult, part1.fdMCParticle().partOriginMCTruth(), part2.fdMCParticle().partOriginMCTruth(), smearingByOrigin); } else { @@ -360,19 +399,59 @@ class FemtoDreamContainer } } - /// Pass a pair to the container and compute all the relevant observables in divided qn bins + /// Pass a pair to the container and compute all the relevant observables in divided qn&phi-psi bins template - void setPair_qn_base(const float femtoObs, const float mT, const float multPercentile, const int myQnBin, const int numQnBins = 10) + void setPair_EP_base(const float femtoObs, const float mT, const float multPercentile, const float myEPObs) { - if (myQnBin >= 0 && myQnBin < numQnBins) { - mHistogramRegistry->fill(HIST(mFolderSuffix[mEventType]) + HIST(o2::aod::femtodreamMCparticle::MCTypeName[mc]) + HIST("_qn") + HIST("/relPairkstarmTMultMultPercentileQn"), femtoObs, mT, multPercentile, myQnBin); - } else { - return; + mHistogramRegistry->fill(HIST(mFolderSuffix[mEventType]) + HIST(o2::aod::femtodreamMCparticle::MCTypeName[mc]) + HIST("_EP") + HIST("/relPairkstarmTMultMultPercentileQn"), femtoObs, mT, multPercentile, myEPObs); + } + + template + void setPair_EP(T1 const& part1, T2 const& part2, const float multPercentile, const bool doQnSeparation, float myEPObs) + { + float femtoObs, femtoObsMC; + // Calculate femto observable and the mT with reconstructed information + if constexpr (mFemtoObs == femtoDreamContainer::Observable::kstar) { + femtoObs = FemtoDreamMath::getkstar(part1, mMassOne, part2, mMassTwo); + } + if (mHighkstarCut > 0) { + if (femtoObs > mHighkstarCut) { + return; + } + } + const float mT = FemtoDreamMath::getmT(part1, mMassOne, part2, mMassTwo); + + if (!doQnSeparation) { + myEPObs = FemtoDreamMath::getPairPhiEP(part1, mMassOne, part2, mMassTwo, myEPObs); + } + + if (mHistogramRegistry) { + setPair_EP_base(femtoObs, mT, multPercentile, myEPObs); + + if constexpr (isMC) { + if (part1.has_fdMCParticle() && part2.has_fdMCParticle()) { + // calculate the femto observable and the mT with MC truth information + if constexpr (mFemtoObs == femtoDreamContainer::Observable::kstar) { + femtoObsMC = FemtoDreamMath::getkstar(part1.fdMCParticle(), mMassOne, part2.fdMCParticle(), mMassTwo); + } + const float mTMC = FemtoDreamMath::getmT(part1.fdMCParticle(), mMassOne, part2.fdMCParticle(), mMassTwo); + + if (std::abs(part1.fdMCParticle().pdgMCTruth()) == mPDGOne && std::abs(part2.fdMCParticle().pdgMCTruth()) == mPDGTwo) { // Note: all pair-histogramms are filled with MC truth information ONLY in case of non-fake candidates + setPair_EP_base(femtoObsMC, mTMC, multPercentile, myEPObs); + } else { + mHistogramRegistry->fill(HIST(mFolderSuffix[mEventType]) + HIST(o2::aod::femtodreamMCparticle::MCTypeName[o2::aod::femtodreamMCparticle::MCType::kTruth]) + HIST("/hFakePairsCounter"), 0); + } + + } else { + mHistogramRegistry->fill(HIST(mFolderSuffix[mEventType]) + HIST(o2::aod::femtodreamMCparticle::MCTypeName[o2::aod::femtodreamMCparticle::MCType::kTruth]) + HIST("/hNoMCtruthPairsCounter"), 0); + } + } } } + // while doing mixing for EP, we have to compute the phi angular of a pair according to plane-calibarated second particle template - void setPair_qn(T1 const& part1, T2 const& part2, const float multPercentile, const int myQnBin, const int numQnBins = 10) + void setPair_EP(T1 const& part1, T2 const& part2, const float multPercentile, const bool doQnSeparation, float EP1, float EP2) { float femtoObs, femtoObsMC; // Calculate femto observable and the mT with reconstructed information @@ -386,8 +465,12 @@ class FemtoDreamContainer } const float mT = FemtoDreamMath::getmT(part1, mMassOne, part2, mMassTwo); + if (doQnSeparation) + LOG(fatal) << "Wrong usage of setPair_EP_mix, this is only for phi-psi mixing"; + EP1 = FemtoDreamMath::getPairPhiEP(part1, mMassOne, part2, mMassTwo, EP1, EP2); + if (mHistogramRegistry) { - setPair_qn_base(femtoObs, mT, multPercentile, myQnBin, numQnBins); + setPair_EP_base(femtoObs, mT, multPercentile, EP1); if constexpr (isMC) { if (part1.has_fdMCParticle() && part2.has_fdMCParticle()) { @@ -397,8 +480,87 @@ class FemtoDreamContainer } const float mTMC = FemtoDreamMath::getmT(part1.fdMCParticle(), mMassOne, part2.fdMCParticle(), mMassTwo); - if (abs(part1.fdMCParticle().pdgMCTruth()) == mPDGOne && abs(part2.fdMCParticle().pdgMCTruth()) == mPDGTwo) { // Note: all pair-histogramms are filled with MC truth information ONLY in case of non-fake candidates - setPair_qn_base(femtoObsMC, mTMC, multPercentile, myQnBin, numQnBins); + if (std::abs(part1.fdMCParticle().pdgMCTruth()) == mPDGOne && std::abs(part2.fdMCParticle().pdgMCTruth()) == mPDGTwo) { // Note: all pair-histogramms are filled with MC truth information ONLY in case of non-fake candidates + setPair_EP_base(femtoObsMC, mTMC, multPercentile, EP1); + } else { + mHistogramRegistry->fill(HIST(mFolderSuffix[mEventType]) + HIST(o2::aod::femtodreamMCparticle::MCTypeName[o2::aod::femtodreamMCparticle::MCType::kTruth]) + HIST("/hFakePairsCounter"), 0); + } + + } else { + mHistogramRegistry->fill(HIST(mFolderSuffix[mEventType]) + HIST(o2::aod::femtodreamMCparticle::MCTypeName[o2::aod::femtodreamMCparticle::MCType::kTruth]) + HIST("/hNoMCtruthPairsCounter"), 0); + } + } + } + } + + /// Pass a pair to the container and compute all the relevant observables in divided qn bins + template + void setPair_3Dqn_base(const float femtoDKout, const float femtoDKside, const float femtoDKlong, const float mT, const float multPercentile, const float myQnBin, const float pairPhiEP) + { + mHistogramRegistry->fill(HIST(mFolderSuffix[mEventType]) + HIST(o2::aod::femtodreamMCparticle::MCTypeName[mc]) + HIST("_3Dqn") + HIST("/relPair3dRmTMultPercentileQnPairphi"), femtoDKout, femtoDKside, femtoDKlong, mT, multPercentile, myQnBin, pairPhiEP); + } + + template + void setPair_3Dqn(T1 const& part1, T2 const& part2, const float multPercentile, bool IsSameSpecies, const float myQnBin, const float eventPlane) + { + + std::vector k3d = FemtoDreamMath::newpairfunc(part1, mMassOne, part2, mMassTwo, IsSameSpecies); + float DKout = k3d[1]; + float DKside = k3d[2]; + float DKlong = k3d[3]; + + const float mT = FemtoDreamMath::getmT(part1, mMassOne, part2, mMassTwo); + + const float pairPhiEP = FemtoDreamMath::getPairPhiEP(part1, mMassOne, part2, mMassTwo, eventPlane); + + if (mHistogramRegistry) { + setPair_3Dqn_base(DKout, DKside, DKlong, mT, multPercentile, myQnBin, pairPhiEP); + + if constexpr (isMC) { + if (part1.has_fdMCParticle() && part2.has_fdMCParticle()) { + + std::vector k3dMC = FemtoDreamMath::newpairfunc(part1.fdMCParticle(), mMassOne, part2.fdMCParticle(), mMassTwo, IsSameSpecies); + const float mTMC = FemtoDreamMath::getmT(part1.fdMCParticle(), mMassOne, part2.fdMCParticle(), mMassTwo); + const float pairPhiEPMC = FemtoDreamMath::getPairPhiEP(part1.fdMCParticle(), mMassOne, part2.fdMCParticle(), mMassTwo, eventPlane); + + if (std::abs(part1.fdMCParticle().pdgMCTruth()) == mPDGOne && std::abs(part2.fdMCParticle().pdgMCTruth()) == mPDGTwo) { // Note: all pair-histogramms are filled with MC truth information ONLY in case of non-fake candidates + setPair_3Dqn_base(k3dMC[1], k3dMC[2], k3dMC[3], mTMC, multPercentile, myQnBin, pairPhiEPMC); + } else { + mHistogramRegistry->fill(HIST(mFolderSuffix[mEventType]) + HIST(o2::aod::femtodreamMCparticle::MCTypeName[o2::aod::femtodreamMCparticle::MCType::kTruth]) + HIST("/hFakePairsCounter"), 0); + } + + } else { + mHistogramRegistry->fill(HIST(mFolderSuffix[mEventType]) + HIST(o2::aod::femtodreamMCparticle::MCTypeName[o2::aod::femtodreamMCparticle::MCType::kTruth]) + HIST("/hNoMCtruthPairsCounter"), 0); + } + } + } + } + + template + void setPair_3Dqn(T1 const& part1, T2 const& part2, const float multPercentile, bool IsSameSpecies, const float myQnBin, const float EP1, const float EP2) + { + + std::vector k3d = FemtoDreamMath::newpairfunc(part1, mMassOne, part2, mMassTwo, IsSameSpecies); + float DKout = k3d[1]; + float DKside = k3d[2]; + float DKlong = k3d[3]; + + const float mT = FemtoDreamMath::getmT(part1, mMassOne, part2, mMassTwo); + + const float pairPhiEP = FemtoDreamMath::getPairPhiEP(part1, mMassOne, part2, mMassTwo, EP1, EP2); + + if (mHistogramRegistry) { + setPair_3Dqn_base(DKout, DKside, DKlong, mT, multPercentile, myQnBin, pairPhiEP); + + if constexpr (isMC) { + if (part1.has_fdMCParticle() && part2.has_fdMCParticle()) { + + std::vector k3dMC = FemtoDreamMath::newpairfunc(part1.fdMCParticle(), mMassOne, part2.fdMCParticle(), mMassTwo, IsSameSpecies); + const float mTMC = FemtoDreamMath::getmT(part1.fdMCParticle(), mMassOne, part2.fdMCParticle(), mMassTwo); + const float pairPhiEPMC = FemtoDreamMath::getPairPhiEP(part1.fdMCParticle(), mMassOne, part2.fdMCParticle(), mMassTwo, EP1, EP2); + + if (std::abs(part1.fdMCParticle().pdgMCTruth()) == mPDGOne && std::abs(part2.fdMCParticle().pdgMCTruth()) == mPDGTwo) { // Note: all pair-histogramms are filled with MC truth information ONLY in case of non-fake candidates + setPair_3Dqn_base(k3dMC[1], k3dMC[2], k3dMC[3], mTMC, multPercentile, myQnBin, pairPhiEPMC); } else { mHistogramRegistry->fill(HIST(mFolderSuffix[mEventType]) + HIST(o2::aod::femtodreamMCparticle::MCTypeName[o2::aod::femtodreamMCparticle::MCType::kTruth]) + HIST("/hFakePairsCounter"), 0); } diff --git a/PWGCF/FemtoDream/Core/femtoDreamMath.h b/PWGCF/FemtoDream/Core/femtoDreamMath.h index aba07f96c9c..02d8d715a34 100644 --- a/PWGCF/FemtoDream/Core/femtoDreamMath.h +++ b/PWGCF/FemtoDream/Core/femtoDreamMath.h @@ -16,12 +16,15 @@ #ifndef PWGCF_FEMTODREAM_CORE_FEMTODREAMMATH_H_ #define PWGCF_FEMTODREAM_CORE_FEMTODREAMMATH_H_ -#include - -#include "Math/Vector4D.h" #include "Math/Boost.h" +#include "Math/Vector4D.h" #include "TLorentzVector.h" #include "TMath.h" +#include "TVector2.h" +#include + +#include +#include namespace o2::analysis::femtoDream { @@ -149,6 +152,124 @@ class FemtoDreamMath return casc.M(); } + + /// Compute the 3d components of the pair momentum in LCMS and PRF + /// Copy from femto universe + /// \tparam T type of tracks + /// \param part1 Particle 1 + /// \param mass1 Mass of particle 1 + /// \param part2 Particle 2 + /// \param mass2 Mass of particle 2 + /// \param isiden Identical or non-identical particle pair + template + static std::vector newpairfunc(const T& part1, const float mass1, const T& part2, const float mass2, bool isiden) + { + const double e1 = std::sqrt(std::pow(part1.px(), 2) + std::pow(part1.py(), 2) + std::pow(part1.pz(), 2) + std::pow(mass1, 2)); + const double e2 = std::sqrt(std::pow(part2.px(), 2) + std::pow(part2.py(), 2) + std::pow(part2.pz(), 2) + std::pow(mass2, 2)); + + const ROOT::Math::PxPyPzEVector vecpart1(part1.px(), part1.py(), part1.pz(), e1); + const ROOT::Math::PxPyPzEVector vecpart2(part2.px(), part2.py(), part2.pz(), e2); + const ROOT::Math::PxPyPzEVector trackSum = vecpart1 + vecpart2; + + std::vector vect; + + const double tPx = trackSum.px(); + const double tPy = trackSum.py(); + const double tPz = trackSum.pz(); + const double tE = trackSum.E(); + + const double tPtSq = (tPx * tPx + tPy * tPy); + const double tMtSq = (tE * tE - tPz * tPz); + const double tM = std::sqrt(tMtSq - tPtSq); + const double tMt = std::sqrt(tMtSq); + const double tPt = std::sqrt(tPtSq); + + // Boost to LCMS + + const double beta_LCMS = tPz / tE; + const double gamma_LCMS = tE / tMt; + + const double fDKOut = (part1.px() * tPx + part1.py() * tPy) / tPt; + const double fDKSide = (-part1.px() * tPy + part1.py() * tPx) / tPt; + const double fDKLong = gamma_LCMS * (part1.pz() - beta_LCMS * e1); + const double fDE = gamma_LCMS * (e1 - beta_LCMS * part1.pz()); + + const double px1LCMS = fDKOut; + const double py1LCMS = fDKSide; + const double pz1LCMS = fDKLong; + const double pE1LCMS = fDE; + + const double px2LCMS = (part2.px() * tPx + part2.py() * tPy) / tPt; + const double py2LCMS = (part2.py() * tPx - part2.px() * tPy) / tPt; + const double pz2LCMS = gamma_LCMS * (part2.pz() - beta_LCMS * e2); + const double pE2LCMS = gamma_LCMS * (e2 - beta_LCMS * part2.pz()); + + const double fDKOutLCMS = px1LCMS - px2LCMS; + const double fDKSideLCMS = py1LCMS - py2LCMS; + const double fDKLongLCMS = pz1LCMS - pz2LCMS; + + // Boost to PRF + + const double betaOut = tPt / tMt; + const double gammaOut = tMt / tM; + + const double fDKOutPRF = gammaOut * (fDKOutLCMS - betaOut * (pE1LCMS - pE2LCMS)); + const double fDKSidePRF = fDKSideLCMS; + const double fDKLongPRF = fDKLongLCMS; + const double fKOut = gammaOut * (fDKOut - betaOut * fDE); + + const double qlcms = std::sqrt(fDKOutLCMS * fDKOutLCMS + fDKSideLCMS * fDKSideLCMS + fDKLongLCMS * fDKLongLCMS); + const double qinv = std::sqrt(fDKOutPRF * fDKOutPRF + fDKSidePRF * fDKSidePRF + fDKLongPRF * fDKLongPRF); + const double kstar = std::sqrt(fKOut * fKOut + fDKSide * fDKSide + fDKLong * fDKLong); + + if (isiden) { + vect.push_back(qinv); + vect.push_back(fDKOutLCMS); + vect.push_back(fDKSideLCMS); + vect.push_back(fDKLongLCMS); + vect.push_back(qlcms); + } else { + vect.push_back(kstar); + vect.push_back(fDKOut); + vect.push_back(fDKSide); + vect.push_back(fDKLong); + } + return vect; + } + + /// Compute the phi angular of a pair with respect to the event plane + /// \tparam T type of tracks + /// \param part1 Particle 1 + /// \param mass1 Mass of particle 1 + /// \param part2 Particle 2 + /// \param mass2 Mass of particle 2 + template + static float getPairPhiEP(const T1& part1, const float mass1, const T2& part2, const float mass2, const float Psi_ep) + { + const ROOT::Math::PtEtaPhiMVector vecpart1(part1.pt(), part1.eta(), part1.phi(), mass1); + const ROOT::Math::PtEtaPhiMVector vecpart2(part2.pt(), part2.eta(), part2.phi(), mass2); + const ROOT::Math::PtEtaPhiMVector trackSum = vecpart1 + vecpart2; + float phi_pair_onPsi = TVector2::Phi_mpi_pi(trackSum.Phi() - Psi_ep); + phi_pair_onPsi = TMath::Abs(phi_pair_onPsi); + return phi_pair_onPsi; + } + + /// Compute the phi angular of a pair according to plane-calibarated second particle + template + static float getPairPhiEP(const T1& part1, const float mass1, const T2& part2, const float mass2, const float Psi_ep1, const float Psi_ep2) + { + const ROOT::Math::PtEtaPhiMVector vecpart1(part1.pt(), part1.eta(), part1.phi(), mass1); + const ROOT::Math::PtEtaPhiMVector vecpart2(part2.pt(), part2.eta(), part2.phi(), mass2); + const float psidiff = Psi_ep2 - Psi_ep1; + // rotate phi of part2 + const float newPhi2 = TVector2::Phi_mpi_pi(vecpart2.Phi() - psidiff); + const ROOT::Math::PtEtaPhiMVector vecpart2_calibd(vecpart2.Pt(), vecpart2.Eta(), newPhi2, vecpart2.M()); + const ROOT::Math::PtEtaPhiMVector trackSum = vecpart1 + vecpart2_calibd; + // reCalibrate part2 phi with respect to part1 event plane + float phi_pair_onPsi = TVector2::Phi_mpi_pi(trackSum.Phi() - Psi_ep1); + phi_pair_onPsi = TMath::Abs(phi_pair_onPsi); + return phi_pair_onPsi; + } }; } // namespace o2::analysis::femtoDream diff --git a/PWGCF/FemtoDream/TableProducer/femtoDreamProducerTask.cxx b/PWGCF/FemtoDream/TableProducer/femtoDreamProducerTask.cxx index a0b362851ba..5ad7f46fc3b 100644 --- a/PWGCF/FemtoDream/TableProducer/femtoDreamProducerTask.cxx +++ b/PWGCF/FemtoDream/TableProducer/femtoDreamProducerTask.cxx @@ -95,6 +95,7 @@ struct femtoDreamProducerTask { Produces outputCollision; Produces outputExtQnCollision; + Produces outputExtEPCollision; Produces outputMCCollision; Produces outputCollsMCLabels; Produces outputParts; @@ -264,16 +265,16 @@ struct femtoDreamProducerTask { } rctCut; struct : o2::framework::ConfigurableGroup { - std::string prefix = std::string("qnCal"); - Configurable ConfFlowCalculate{"ConfFlowCalculate", false, "Evt sel: Cumulant of flow"}; // To do - Configurable ConfQnSeparation{"ConfQnSeparation", false, "Evt sel: Qn of event"}; - Configurable> ConfQnBinSeparator{"ConfQnBinSeparator", std::vector{-999.f, -999.f, -999.f}, "Qn bin separator"}; - Configurable ConfdoFillHisto{"ConfdoFillHisto", false, "Fill histos for Qn and sphericity and mult "}; - Configurable ConfCentralityMax{"ConfCentralityMax", 80.f, "Evt sel: Maximum Centrality cut"}; + std::string prefix = std::string("epCal"); + Configurable ConfFillFlowQA{"ConfFillFlowQA", false, "Evt sel: fill flow/event-plane related observables"}; + Configurable ConfQnSeparation{"ConfQnSeparation", false, "Evt sel: do qn separation of events"}; + Configurable ConfHarmonicOrder{"ConfHarmonicOrder", 2, "harmonic order for event plane calculation"}; + Configurable> ConfQnBinSeparator{"ConfQnBinSeparator", std::vector{-999.f, -999.f, -999.f}, "qn bin separator"}; + Configurable ConfDoCumlant{"ConfDoCumlant", false, "do cumulant for flow calculation"}; + Configurable ConfCentralityMax{"ConfCentralityMax", 100.f, "Evt sel: Maximum Centrality cut"}; Configurable ConfCentBinWidth{"ConfCentBinWidth", 1.f, "Centrality bin length for qn separator"}; - Configurable ConfQnBinMin{"ConfQnBinMin", 0, "Minimum qn bin"}; Configurable ConfNumQnBins{"ConfNumQnBins", 10, "Number of qn bins"}; - } qnCal; + } epCal; struct : o2::framework::ConfigurableGroup { std::string prefix = std::string("OptionEvtSpecialSelections"); @@ -290,7 +291,7 @@ struct femtoDreamProducerTask { HistogramRegistry V0Registry{"V0", {}, OutputObjHandlingPolicy::AnalysisObject}; HistogramRegistry ResoRegistry{"Reso", {}, OutputObjHandlingPolicy::AnalysisObject}; HistogramRegistry CascadeRegistry{"Cascade", {}, OutputObjHandlingPolicy::AnalysisObject}; - HistogramRegistry FlowRegistry{"Qn", {}, OutputObjHandlingPolicy::AnalysisObject}; + HistogramRegistry FlowRegistry{"Flow", {}, OutputObjHandlingPolicy::AnalysisObject}; int mRunNumber; float mMagField; @@ -299,10 +300,10 @@ struct femtoDreamProducerTask { void init(InitContext&) { - if (doprocessData == false && doprocessData_noCentrality == false && doprocessData_CentPbPb == false && doprocessData_CentPbPb_qvec == false && doprocessMC == false && doprocessMC_noCentrality == false && doprocessMC_CentPbPb == false) { + if (doprocessData == false && doprocessData_noCentrality == false && doprocessData_CentPbPb == false && doprocessData_CentPbPb_EP == false && doprocessMC == false && doprocessMC_noCentrality == false && doprocessMC_CentPbPb == false) { LOGF(fatal, "Neither processData nor processMC enabled. Please choose one."); } - if ((doprocessData == true && doprocessMC == true) || (doprocessData == true && doprocessMC_noCentrality == true) || (doprocessMC == true && doprocessMC_noCentrality == true) || (doprocessData_noCentrality == true && doprocessData == true) || (doprocessData_noCentrality == true && doprocessMC == true) || (doprocessData_noCentrality == true && doprocessMC_noCentrality == true) || (doprocessData_CentPbPb == true && doprocessData == true) || (doprocessData_CentPbPb == true && doprocessData_noCentrality == true) || (doprocessData_CentPbPb == true && doprocessMC == true) || (doprocessData_CentPbPb == true && doprocessMC_noCentrality == true) || (doprocessData_CentPbPb == true && doprocessMC_CentPbPb == true) || (doprocessData_CentPbPb_qvec == true && doprocessData == true) || (doprocessData_CentPbPb_qvec == true && doprocessData_noCentrality == true) || (doprocessData_CentPbPb_qvec == true && doprocessMC == true) || (doprocessData_CentPbPb_qvec == true && doprocessMC_noCentrality == true) || (doprocessData_CentPbPb_qvec == true && doprocessMC_CentPbPb == true) || (doprocessData_CentPbPb_qvec == true && doprocessData_CentPbPb == true)) { + if ((doprocessData == true && doprocessMC == true) || (doprocessData == true && doprocessMC_noCentrality == true) || (doprocessMC == true && doprocessMC_noCentrality == true) || (doprocessData_noCentrality == true && doprocessData == true) || (doprocessData_noCentrality == true && doprocessMC == true) || (doprocessData_noCentrality == true && doprocessMC_noCentrality == true) || (doprocessData_CentPbPb == true && doprocessData == true) || (doprocessData_CentPbPb == true && doprocessData_noCentrality == true) || (doprocessData_CentPbPb == true && doprocessMC == true) || (doprocessData_CentPbPb == true && doprocessMC_noCentrality == true) || (doprocessData_CentPbPb == true && doprocessMC_CentPbPb == true) || (doprocessData_CentPbPb_EP == true && doprocessData == true) || (doprocessData_CentPbPb_EP == true && doprocessData_noCentrality == true) || (doprocessData_CentPbPb_EP == true && doprocessMC == true) || (doprocessData_CentPbPb_EP == true && doprocessMC_noCentrality == true) || (doprocessData_CentPbPb_EP == true && doprocessMC_CentPbPb == true) || (doprocessData_CentPbPb_EP == true && doprocessData_CentPbPb == true)) { LOGF(fatal, "Cannot enable more than one process switch at the same time. " "Please choose one."); @@ -445,10 +446,9 @@ struct femtoDreamProducerTask { } } - if (qnCal.ConfFlowCalculate) { - colCuts.initFlow(&FlowRegistry, qnCal.ConfQnSeparation); - if (qnCal.ConfQnSeparation) - colCuts.initQn(&FlowRegistry, qnCal.ConfNumQnBins); + if (epCal.ConfFillFlowQA) { + colCuts.initFlow(&FlowRegistry, epCal.ConfQnSeparation); + colCuts.initEPQA(&FlowRegistry); } mRunNumber = 0; @@ -763,7 +763,7 @@ struct femtoDreamProducerTask { } if constexpr (doFlow) { - fillCollisionsFlow(col, tracks, mult, spher, multNtr); + fillCollisionsFlow(col, tracks, mult, spher, epCal.ConfHarmonicOrder); } std::vector childIDs = {0, 0}; // these IDs are necessary to keep track of the children @@ -1115,21 +1115,27 @@ struct femtoDreamProducerTask { } template - void fillCollisionsFlow(CollisionType const& col, TrackType const& tracks, float mult, float spher, float multNtr) + void fillCollisionsFlow(CollisionType const& col, TrackType const& tracks, float mult, float spher, int EPHarmonic) { - float myqn = -999.; - // Calculate and fill qn values - if (qnCal.ConfQnSeparation) { - myqn = colCuts.computeqnVec(col); + float myqn = colCuts.computeqnVec(col); + float myEP = TMath::RadToDeg() * colCuts.computeEP(col, EPHarmonic); + // psi from rad(0-pi) to deg, in table psi would be in deg,from 0-180 + + if ((myqn >= 0 && myqn < 1e6) || (myEP >= 0 && myEP < 180)) { outputExtQnCollision(myqn, col.trackOccupancyInTimeRange()); + outputExtEPCollision(myEP); } + // Calculate flow via cumulant - if (qnCal.ConfFlowCalculate) { - int qnBin = colCuts.myqnBin(mult, qnCal.ConfCentralityMax, qnCal.ConfQnBinSeparator, qnCal.ConfdoFillHisto, spher, myqn, qnCal.ConfNumQnBins, multNtr, qnCal.ConfCentBinWidth); - if (qnBin < qnCal.ConfQnBinMin || qnBin > qnCal.ConfNumQnBins) { - qnBin = -999; + + if (epCal.ConfQnSeparation) { + colCuts.myqnBin(mult, epCal.ConfCentralityMax, epCal.ConfFillFlowQA, epCal.ConfQnBinSeparator, myqn, epCal.ConfNumQnBins, epCal.ConfCentBinWidth); + } + if (epCal.ConfFillFlowQA) { + colCuts.fillEPQA(mult, spher, myqn, myEP); + if (epCal.ConfDoCumlant) { + colCuts.doCumulants(col, tracks, mult, epCal.ConfQnSeparation); } - colCuts.doCumulants(col, tracks, mult, qnCal.ConfQnSeparation); } } @@ -1195,11 +1201,11 @@ struct femtoDreamProducerTask { PROCESS_SWITCH(femtoDreamProducerTask, processData_CentPbPb, "Provide experimental data with centrality information for PbPb collisions", false); - void processData_CentPbPb_qvec(aod::FemtoFullCollision_CentPbPb_qvec const& col, - aod::BCsWithTimestamps const&, - aod::FemtoFullTracks const& tracks, - o2::aod::V0Datas const& fullV0s, - o2::aod::CascDatas const& fullCascades) + void processData_CentPbPb_EP(aod::FemtoFullCollision_CentPbPb_qvec const& col, + aod::BCsWithTimestamps const&, + aod::FemtoFullTracks const& tracks, + o2::aod::V0Datas const& fullV0s, + o2::aod::CascDatas const& fullCascades) { // get magnetic field for run initCCDB_Mag_Trig(col.bc_as()); @@ -1212,7 +1218,7 @@ struct femtoDreamProducerTask { fillCollisionsAndTracksAndV0AndCascade(col, tracks, tracks, fullV0s, fullCascades); } } - PROCESS_SWITCH(femtoDreamProducerTask, processData_CentPbPb_qvec, + PROCESS_SWITCH(femtoDreamProducerTask, processData_CentPbPb_EP, "Provide experimental data with centrality and q-vector table for PbPb collisions", false); void processMC(aod::FemtoFullCollisionMC const& col, diff --git a/PWGCF/FemtoDream/Tasks/femtoDreamPairTaskTrackTrack.cxx b/PWGCF/FemtoDream/Tasks/femtoDreamPairTaskTrackTrack.cxx index f75d924683b..f2100431efd 100644 --- a/PWGCF/FemtoDream/Tasks/femtoDreamPairTaskTrackTrack.cxx +++ b/PWGCF/FemtoDream/Tasks/femtoDreamPairTaskTrackTrack.cxx @@ -85,25 +85,33 @@ struct femtoDreamPairTaskTrackTrack { Filter EventMultiplicity = aod::femtodreamcollision::multNtr >= EventSel.MultMin && aod::femtodreamcollision::multNtr <= EventSel.MultMax && aod::femtodreamcollision::sphericity >= EventSel.SphericityMin; Filter EventMultiplicityPercentile = aod::femtodreamcollision::multV0M >= EventSel.MultPercentileMin && aod::femtodreamcollision::multV0M <= EventSel.MultPercentileMax; - /// qn-separator - FemtoDreamCollisionSelection qnBinCalculator; + /// qn&event_plane separator + FemtoDreamCollisionSelection epCalculator; struct : ConfigurableGroup { - std::string prefix = std::string("qnCal"); + std::string prefix = std::string("EPCal"); + Configurable do1DFemto{"do1DFemto", false, "Do 1D femtoscopy"}; + Configurable do3DFemto{"do3DFemto", false, "Do 3D femtoscopy"}; + Configurable fillFlowQA{"fillFlowQA", false, "Fill QA histos for flow/event-plane related observables"}; + Configurable storeEvtTrkInfo{"storeEvtTrkInfo", false, "Fill info of track1 and track2 while pariing in divided qn bins"}; Configurable doQnSeparation{"doQnSeparation", false, "Do qn separation"}; + Configurable doEPReClibForMixing{"doEPReClibForMixing", false, "While mixing, using respective event plane for participating particles azimuthal angle caulculation"}; Configurable> qnBinSeparator{"qnBinSeparator", std::vector{-999.f, -999.f, -999.f}, "Qn bin separator"}; - Configurable doFillHisto{"doFillHisto", false, "Fill histos for Qn and sphericity and mult "}; - Configurable storeEvtTrkInfo{"storeEvtTrkInfo", false, "Fill info of track1 and track2 while pariing in divided qn bins"}; Configurable numQnBins{"numQnBins", 10, "Number of qn bins"}; Configurable qnBinMin{"qnBinMin", 0, "Number of qn bins"}; - Configurable centMax{"centMax", 80.f, "Evt sel: Maximum Centrality cut"}; + Configurable centMax{"centMax", 100.f, "Evt sel: Maximum Centrality cut"}; Configurable centBinWidth{"centBinWidth", 1.f, "Centrality bin length for qn separator"}; - } qnCal; + ConfigurableAxis DKout{"DKout", {500, -2., 2.}, "binning DKout for the 3-D femtoscopy plot: R_out(LCMS) vs mT vs multiplicity percentile vs qnBin vs pait phi wrt EP (set <> to true)"}; + ConfigurableAxis DKside{"DKside", {500, -2., 2.}, "binning DKside for the 3-D femtoscopy plot: R_side(LCMS) vs mT vs multiplicity percentile vs qnBin vs pait phi wrt EP (set <> to true)"}; + ConfigurableAxis DKlong{"DKlong", {500, -2., 2.}, "binning DKlong for the 3-D femtoscopy plot: R_long(LCMS) vs mT vs multiplicity percentile vs qnBin vs pait phi wrt EP (set <> to true)"}; + ConfigurableAxis qnBins{"qnBins", {10, 0, 10}, "binning of qn interval"}; + ConfigurableAxis pairPhiBins{"pairPhiBins", {12, 0., TMath::Pi()}, "binning of pair phi"}; + } EPCal; using FilteredCollisions = soa::Filtered; using FilteredCollision = FilteredCollisions::iterator; using FilteredMCCollisions = soa::Filtered>; using FilteredMCCollision = FilteredMCCollisions::iterator; - using FilteredQnCollisions = soa::Filtered>; + using FilteredQnCollisions = soa::Filtered>; using FilteredQnCollision = FilteredQnCollisions::iterator; using FilteredMaskedCollisions = soa::Filtered>; @@ -233,14 +241,16 @@ struct femtoDreamPairTaskTrackTrack { ConfigurableAxis MultPercentileMixBins{"MultPercentileMixBins", {VARIABLE_WIDTH, 0.0f, 10.0f, 20.0f, 30.0f, 40.0f, 50.0f, 60.0f, 70.0f, 80.0f, 90.0f, 100.0f}, "Mixing bins - multiplicity percentile"}; ConfigurableAxis VztxMixBins{"VztxMixBins", {VARIABLE_WIDTH, -10.0f, -8.f, -6.f, -4.f, -2.f, 0.f, 2.f, 4.f, 6.f, 8.f, 10.f}, "Mixing bins - z-vertex"}; ConfigurableAxis QnMixBins{"QnMixBins", {VARIABLE_WIDTH, 0.50f, 68.50f, 100.50f, 126.50f, 151.50f, 176.50f, 203.50f, 232.50f, 269.50f, 322.50f, 833.50f}, "Mixing bins - qn-value"}; + ConfigurableAxis EPMixBins{"EPMixBins", {VARIABLE_WIDTH, 0., 10., 20., 30., 40., 50., 60., 70., 80., 90., 100., 110., 120., 130., 140., 150., 160., 170., 180.}, "Mixing bins - event plane (deg)"}; Configurable Depth{"Depth", 5, "Number of events for mixing"}; - Configurable Policy{"Policy", 0, "Binning policy for mixing - 0: multiplicity, 1: multipliciy percentile, 2: both, 3: multipliciy percentile and qn value"}; + Configurable Policy{"Policy", 0, "Binning policy for mixing - 0: multiplicity, 1: multipliciy percentile, 2: both, 3: multipliciy percentile and qn value, 4: multipliciy percentile and event plane"}; } Mixing; ColumnBinningPolicy colBinningMult{{Mixing.VztxMixBins, Mixing.MultMixBins}, true}; ColumnBinningPolicy colBinningMultPercentile{{Mixing.VztxMixBins, Mixing.MultPercentileMixBins}, true}; ColumnBinningPolicy colBinningMultMultPercentile{{Mixing.VztxMixBins, Mixing.MultMixBins, Mixing.MultPercentileMixBins}, true}; ColumnBinningPolicy colBinningMultPercentileqn{{Mixing.VztxMixBins, Mixing.MultPercentileMixBins, Mixing.QnMixBins}, true}; + ColumnBinningPolicy colBinningMultPercentileEP{{Mixing.VztxMixBins, Mixing.MultPercentileMixBins, Mixing.EPMixBins}, true}; FemtoDreamContainer sameEventCont; FemtoDreamContainer mixedEventCont; @@ -250,6 +260,7 @@ struct femtoDreamPairTaskTrackTrack { // Container for correlation functions in devided qn bins FemtoDreamContainer sameEventQnCont; + FemtoDreamContainer mixedEventQnCont; /// Histogram output HistogramRegistry Registry{"Output", {}, OutputObjHandlingPolicy::AnalysisObject}; @@ -264,6 +275,7 @@ struct femtoDreamPairTaskTrackTrack { colBinningMultPercentile = {{Mixing.VztxMixBins, Mixing.MultPercentileMixBins}, true}; colBinningMultMultPercentile = {{Mixing.VztxMixBins, Mixing.MultMixBins, Mixing.MultPercentileMixBins}, true}; colBinningMultPercentileqn = {{Mixing.VztxMixBins, Mixing.MultPercentileMixBins, Mixing.QnMixBins}, true}; + colBinningMultPercentileEP = {{Mixing.VztxMixBins, Mixing.MultPercentileMixBins, Mixing.EPMixBins}, true}; if (Option.RandomizePair.value) { random = new TRandom3(0); @@ -294,13 +306,29 @@ struct femtoDreamPairTaskTrackTrack { pairCloseRejectionME.init(&Registry, &Registry, Option.CPRdeltaPhiMax.value, Option.CPRdeltaEtaMax.value, Option.CPRPlotPerRadii.value, 2, Option.CPROld.value); } - if (qnCal.doQnSeparation) { - sameEventQnCont.init_qn(&Registry, - Binning4D.kstar, Binning4D.mT, Binning4D.multPercentile, + if (EPCal.do1DFemto) { + sameEventQnCont.init_EP(&Registry, + Binning4D.kstar, Binning4D.mT, Binning4D.multPercentile, EPCal.doQnSeparation ? EPCal.qnBins : EPCal.pairPhiBins, EPCal.doQnSeparation ? "qnBin" : "#phi_{pair}-#Psi_{EP} (rad)", Option.IsMC, Option.HighkstarCut); + mixedEventQnCont.init_EP(&Registry, + Binning4D.kstar, Binning4D.mT, Binning4D.multPercentile, EPCal.doQnSeparation ? EPCal.qnBins : EPCal.pairPhiBins, EPCal.doQnSeparation ? "qnBin" : "#phi_{pair}-#Psi_{EP} (rad)", + Option.IsMC, Option.HighkstarCut); + sameEventQnCont.setPDGCodes(Track1.PDGCode, Track2.PDGCode); + mixedEventQnCont.setPDGCodes(Track1.PDGCode, Track2.PDGCode); + if (EPCal.fillFlowQA) { + epCalculator.initEPQA(&Registry); + } + } + + if (EPCal.do3DFemto) { + sameEventQnCont.init_3Dqn(&Registry, EPCal.DKout, EPCal.DKside, EPCal.DKlong, + Binning4D.mT, Binning4D.multPercentile, Option.IsMC, EPCal.qnBins, EPCal.pairPhiBins); + mixedEventQnCont.init_3Dqn(&Registry, EPCal.DKout, EPCal.DKside, EPCal.DKlong, + Binning4D.mT, Binning4D.multPercentile, Option.IsMC, EPCal.qnBins, EPCal.pairPhiBins); sameEventQnCont.setPDGCodes(Track1.PDGCode, Track2.PDGCode); - if (qnCal.doFillHisto) { - qnBinCalculator.initQn(&Registry, qnCal.numQnBins); + mixedEventQnCont.setPDGCodes(Track1.PDGCode, Track2.PDGCode); + if (EPCal.fillFlowQA) { + epCalculator.initEPQA(&Registry); } } @@ -343,11 +371,11 @@ struct femtoDreamPairTaskTrackTrack { } } if ((doprocessSameEvent && doprocessSameEventMasked) || - (doprocessSameEvent && doprocessSameEventQn) || - (doprocessSameEventMasked && doprocessSameEventQn) || + (doprocessSameEvent && doprocessSameEventEP) || + (doprocessSameEventMasked && doprocessSameEventEP) || (doprocessMixedEvent && doprocessMixedEventMasked) || - (doprocessMixedEvent && doprocessMixedEventQn) || - (doprocessMixedEventMasked && doprocessMixedEventQn) || + (doprocessMixedEvent && doprocessMixedEventEP) || + (doprocessMixedEventMasked && doprocessMixedEventEP) || (doprocessSameEventMC && doprocessSameEventMCMasked) || (doprocessMixedEventMC && doprocessMixedEventMCMasked)) { LOG(fatal) << "Normal and masked processing cannot be activated simultaneously!"; @@ -662,9 +690,9 @@ struct femtoDreamPairTaskTrackTrack { /// This function processes the same event in divided qn bins /// col.multV0M() get the event centrality from ft0c for PbPb data template - void doSameEventQn(PartitionType SliceTrk1, PartitionType SliceTrk2, PartType parts, Collision col) + void doSameEventEP(PartitionType SliceTrk1, PartitionType SliceTrk2, PartType parts, Collision col) { - if (qnCal.storeEvtTrkInfo) { + if (EPCal.storeEvtTrkInfo) { for (auto& part : SliceTrk1) { trackHistoPartOne.fillQA(part, aod::femtodreamparticle::kPt, col.multNtr(), col.multV0M()); } @@ -676,9 +704,17 @@ struct femtoDreamPairTaskTrackTrack { } } - auto myqnBin = qnBinCalculator.myqnBin(col.multV0M(), qnCal.centMax, qnCal.qnBinSeparator, qnCal.doFillHisto, col.sphericity(), col.qnVal(), qnCal.numQnBins, col.multNtr(), qnCal.centBinWidth); - if (myqnBin < qnCal.qnBinMin || myqnBin > qnCal.numQnBins) { - myqnBin = -999; + auto myEP = TMath::DegToRad() * col.eventPlane(); + int myqnBin = -999; + if (EPCal.doQnSeparation || EPCal.do3DFemto) { + myqnBin = epCalculator.myqnBin(col.multV0M(), EPCal.centMax, EPCal.fillFlowQA, EPCal.qnBinSeparator, col.qnVal(), EPCal.numQnBins, EPCal.centBinWidth); + if (myqnBin < EPCal.qnBinMin || myqnBin > EPCal.numQnBins) { + myqnBin = -999; + } + } + + if (EPCal.fillFlowQA) { + epCalculator.fillEPQA(col.multV0M(), col.sphericity(), col.qnVal(), col.eventPlane()); } /// Now build the combinations @@ -698,9 +734,19 @@ struct femtoDreamPairTaskTrackTrack { rand = random->Rndm(); } if (rand <= 0.5) { - sameEventQnCont.setPair_qn(p1, p2, col.multV0M(), myqnBin, qnCal.numQnBins); + if (EPCal.do1DFemto) { + sameEventQnCont.setPair_EP(p1, p2, col.multV0M(), EPCal.doQnSeparation, EPCal.doQnSeparation ? myqnBin + 0.f : myEP); + } + if (EPCal.do3DFemto) { + sameEventQnCont.setPair_3Dqn(p1, p2, col.multV0M(), Option.SameSpecies.value, myqnBin + 0.f, myEP); + } } else { - sameEventQnCont.setPair_qn(p2, p1, col.multV0M(), myqnBin, qnCal.numQnBins); + if (EPCal.do1DFemto) { + sameEventQnCont.setPair_EP(p2, p1, col.multV0M(), EPCal.doQnSeparation, EPCal.doQnSeparation ? myqnBin + 0.f : myEP); + } + if (EPCal.do3DFemto) { + sameEventQnCont.setPair_3Dqn(p2, p1, col.multV0M(), Option.SameSpecies.value, myqnBin + 0.f, myEP); + } } } } else { @@ -714,17 +760,22 @@ struct femtoDreamPairTaskTrackTrack { if (!pairCleaner.isCleanPair(p1, p2, parts)) { continue; } - sameEventQnCont.setPair_qn(p1, p2, col.multV0M(), myqnBin, qnCal.numQnBins); + if (EPCal.do1DFemto) { + sameEventQnCont.setPair_EP(p1, p2, col.multV0M(), EPCal.doQnSeparation, EPCal.doQnSeparation ? myqnBin + 0.f : myEP); + } + if (EPCal.do3DFemto) { + sameEventQnCont.setPair_3Dqn(p1, p2, col.multV0M(), Option.SameSpecies.value, myEP, myqnBin); + } } } } - /// process function for to call doSameEventQn with Data + /// process function for to call doSameEventEP with Data /// \param col subscribe to the collision table (Data) /// \param parts subscribe to the femtoDreamParticleTable - void processSameEventQn(FilteredQnCollision& col, o2::aod::FDParticles& parts) + void processSameEventEP(FilteredQnCollision& col, o2::aod::FDParticles& parts) { - if (qnCal.storeEvtTrkInfo) { + if (EPCal.storeEvtTrkInfo) { fillCollision(col); } auto SliceTrk1 = PartitionTrk1->sliceByCached(aod::femtodreamparticle::fdCollisionId, col.globalIndex(), cache); @@ -732,35 +783,78 @@ struct femtoDreamPairTaskTrackTrack { if (SliceTrk1.size() == 0 && SliceTrk2.size() == 0) { return; } - if (qnCal.doQnSeparation) { - doSameEventQn(SliceTrk1, SliceTrk2, parts, col); + if (EPCal.do1DFemto || EPCal.do3DFemto) { + doSameEventEP(SliceTrk1, SliceTrk2, parts, col); } } - PROCESS_SWITCH(femtoDreamPairTaskTrackTrack, processSameEventQn, "Enable processing same event in divided qn bins", false); + PROCESS_SWITCH(femtoDreamPairTaskTrackTrack, processSameEventEP, "Enable processing same event wrt azimuthal angle and event-plane ", false); + + template + void doMixedEvent_NotMaskedEP(CollisionType& cols, PartType& parts, PartitionType& part1, PartitionType& part2, BinningType policy) + { + for (auto const& [collision1, collision2] : soa::selfCombinations(policy, Mixing.Depth.value, -1, cols, cols)) { + auto SliceTrk1 = part1->sliceByCached(aod::femtodreamparticle::fdCollisionId, collision1.globalIndex(), cache); + auto SliceTrk2 = part2->sliceByCached(aod::femtodreamparticle::fdCollisionId, collision2.globalIndex(), cache); + if (SliceTrk1.size() == 0 || SliceTrk2.size() == 0) { + continue; + } + + auto myEP_event1 = TMath::DegToRad() * collision1.eventPlane(); + auto myEP_event2 = TMath::DegToRad() * collision2.eventPlane(); + + for (auto& [p1, p2] : combinations(CombinationsFullIndexPolicy(SliceTrk1, SliceTrk2))) { + if (Option.CPROn.value) { + if (pairCloseRejectionME.isClosePair(p1, p2, parts, collision1.magField())) { + continue; + } + } + if (EPCal.doEPReClibForMixing) { + if (EPCal.do1DFemto) { + if (EPCal.doQnSeparation) + mixedEventQnCont.setPair_EP(p1, p2, collision1.multV0M(), EPCal.doQnSeparation, 0.f); + else + mixedEventQnCont.setPair_EP(p1, p2, collision1.multV0M(), EPCal.doQnSeparation, myEP_event1, myEP_event2); + } + if (EPCal.do3DFemto) { + mixedEventQnCont.setPair_3Dqn(p1, p2, collision1.multV0M(), Option.SameSpecies.value, 0.f, myEP_event1, myEP_event2); + } + } else { + if (EPCal.do1DFemto) + mixedEventQnCont.setPair_EP(p1, p2, collision1.multV0M(), EPCal.doQnSeparation, EPCal.doQnSeparation ? 0.f : myEP_event1); + if (EPCal.do3DFemto) { + mixedEventQnCont.setPair_3Dqn(p1, p2, collision1.multV0M(), Option.SameSpecies.value, 0.f, myEP_event1); + } + } + } + } + }; /// process function for to call doMixedEvent with Data /// @param cols subscribe to the collisions table (Data) /// @param parts subscribe to the femtoDreamParticleTable - void processMixedEventQn(FilteredQnCollisions& cols, o2::aod::FDParticles& parts) + void processMixedEventEP(FilteredQnCollisions& cols, o2::aod::FDParticles& parts) { switch (Mixing.Policy.value) { case femtodreamcollision::kMult: - doMixedEvent_NotMasked(cols, parts, PartitionTrk1, PartitionTrk2, colBinningMult); + doMixedEvent_NotMaskedEP(cols, parts, PartitionTrk1, PartitionTrk2, colBinningMult); break; case femtodreamcollision::kMultPercentile: - doMixedEvent_NotMasked(cols, parts, PartitionTrk1, PartitionTrk2, colBinningMultPercentile); + doMixedEvent_NotMaskedEP(cols, parts, PartitionTrk1, PartitionTrk2, colBinningMultPercentile); break; case femtodreamcollision::kMultMultPercentile: - doMixedEvent_NotMasked(cols, parts, PartitionTrk1, PartitionTrk2, colBinningMultMultPercentile); + doMixedEvent_NotMaskedEP(cols, parts, PartitionTrk1, PartitionTrk2, colBinningMultMultPercentile); break; case femtodreamcollision::kMultPercentileQn: - doMixedEvent_NotMasked(cols, parts, PartitionTrk1, PartitionTrk2, colBinningMultPercentileqn); + doMixedEvent_NotMaskedEP(cols, parts, PartitionTrk1, PartitionTrk2, colBinningMultPercentileqn); + break; + case femtodreamcollision::kMultPercentileEP: + doMixedEvent_NotMaskedEP(cols, parts, PartitionTrk1, PartitionTrk2, colBinningMultPercentileEP); break; default: LOG(fatal) << "Invalid binning policiy specifed. Breaking..."; } } - PROCESS_SWITCH(femtoDreamPairTaskTrackTrack, processMixedEventQn, "Enable processing mixed events", false); + PROCESS_SWITCH(femtoDreamPairTaskTrackTrack, processMixedEventEP, "Enable processing mixed events wrt azimuthal angle and event-plane", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) {