diff --git a/PWGCF/FemtoDream/Core/femtoDreamCollisionSelection.h b/PWGCF/FemtoDream/Core/femtoDreamCollisionSelection.h index 4e67c43b40c..398847e7ada 100644 --- a/PWGCF/FemtoDream/Core/femtoDreamCollisionSelection.h +++ b/PWGCF/FemtoDream/Core/femtoDreamCollisionSelection.h @@ -229,25 +229,26 @@ class FemtoDreamCollisionSelection /// Initializes histograms for the flow calculation /// \param registry Histogram registry to be passed - void initFlow(HistogramRegistry* registry, bool doQnSeparation, int mumQnBins = 10, int binPt = 100, int binEta = 32) + void initFlow(HistogramRegistry* registry, bool doQnSeparation, int mumQnBins = 10, int centBins = 10) { if (!mCutsSet) { LOGF(error, "Event selection not set - quitting!"); } - mReQthisEvt = new TH2D("ReQthisEvt", "", binPt, 0., 5., binEta, -0.8, 0.8); - mImQthisEvt = new TH2D("ImQthisEvt", "", binPt, 0., 5., binEta, -0.8, 0.8); - mReQ2thisEvt = new TH2D("ReQ2thisEvt", "", binPt, 0., 5., binEta, -0.8, 0.8); - mImQ2thisEvt = new TH2D("ImQ2thisEvt", "", binPt, 0., 5., binEta, -0.8, 0.8); - mMQthisEvt = new TH2D("MQthisEvt", "", binPt, 0., 5., binEta, -0.8, 0.8); - mMQWeightthisEvt = new TH2D("MQWeightthisEvt", "", binPt, 0., 5., binEta, -0.8, 0.8); mHistogramQn = registry; - mHistogramQn->add("Event/profileC22", "; cent; c22", kTProfile, {{10, 0, 100}}, "s"); - mHistogramQn->add("Event/profileC24", "; cent; c24", kTProfile, {{10, 0, 100}}, "s"); + mHistogramQn->add("Event/hN2allQn", ";centrality; #sum Re(Q_{2,A} Q_{2,B}^{*})", kTH1F, {{centBins, 0, 100}}); + mHistogramQn->add("Event/hD2allQn", ";centrality; #sum (W_{A} W_{B})", kTH1F, {{centBins, 0, 100}}); + mHistogramQn->get(HIST("Event/hN2allQn"))->Sumw2(); + mHistogramQn->get(HIST("Event/hD2allQn"))->Sumw2(); 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}}, "s")); + hN2.push_back(mHistogramQn->add(("Qn/hN2_" + std::to_string(iqn)).c_str(), ";centrality; #sum Re(Q_{2,A} Q_{2,B}^{*})", kTH1F, {{centBins, 0, 100}})); + hD2.push_back(mHistogramQn->add(("Qn/hD2_" + std::to_string(iqn)).c_str(), ";centrality; #sum (W_{A} W_{B})", kTH1F, {{centBins, 0, 100}})); + } + for (int iqn(0); iqn < mumQnBins; ++iqn) { + std::get>(hN2[iqn])->Sumw2(); + std::get>(hD2[iqn])->Sumw2(); } } return; @@ -421,45 +422,6 @@ class FemtoDreamCollisionSelection mHistogramQn->fill(HIST("Event/psiEP"), psiEP); } - /// \todo to be implemented! - /// Fill cumulants histo for flow calculation - /// Reset hists event-by-event - /// \tparam T1 type of the collision - /// \tparam T2 type of the tracks - /// \param tracks All tracks - template - bool fillCumulants(T1 const& col, T2 const& tracks, float fHarmonic = 2.f) - { - int numOfTracks = col.numContrib(); - if (numOfTracks < 3) - return false; - - mReQthisEvt->Reset(); - mImQthisEvt->Reset(); - mReQ2thisEvt->Reset(); - mImQ2thisEvt->Reset(); - mMQthisEvt->Reset(); - mMQWeightthisEvt->Reset(); - - for (auto const& track : tracks) { - double weight = 1; // Will implement NUA&NUE correction - double phi = track.phi(); - double pt = track.pt(); - double eta = track.eta(); - double cosnphi = weight * TMath::Cos(fHarmonic * phi); - double sinnphi = weight * TMath::Sin(fHarmonic * phi); - double cos2nphi = weight * TMath::Cos(2 * fHarmonic * phi); - double sin2nphi = weight * TMath::Sin(2 * fHarmonic * phi); - mReQthisEvt->Fill(pt, eta, cosnphi); - mImQthisEvt->Fill(pt, eta, sinnphi); - mReQ2thisEvt->Fill(pt, eta, cos2nphi); - mImQ2thisEvt->Fill(pt, eta, sin2nphi); - mMQthisEvt->Fill(pt, eta); - mMQWeightthisEvt->Fill(pt, eta, weight); - } - return true; - } - /// \todo to be implemented! /// Do cumulants for flow calculation /// \tparam T type of the collision @@ -467,38 +429,59 @@ class FemtoDreamCollisionSelection /// \param qnBin should be in 0-9 /// \param fEtaGap eta gap for flow cumulant template - void doCumulants(T1 const& col, T2 const& tracks, float centrality, bool doQnSeparation = false, int numQnBins = 10, float fEtaGap = 0.3f, int binPt = 100, int binEta = 32) + void doCumulants(T1 const& col, T2 const& tracks, float centrality, bool doQnSeparation = false, int numQnBins = 10, float fEtaGap = 0.5f, float ptMin = 0.2f, float ptMax = 5.0f, float harmonic = 2.0f) { - if (!fillCumulants(col, tracks)) + int numOfTracks = col.numContrib(); + if (numOfTracks < 3) return; - if (mMQthisEvt->Integral(1, binPt, 1, binEta) < 2) - return; + double Q2A_re = 0., Q2A_im = 0., WA = 0.; + double Q2B_re = 0., Q2B_im = 0., WB = 0.; - double allReQ = mReQthisEvt->Integral(1, binPt, 1, binEta); - double allImQ = mImQthisEvt->Integral(1, binPt, 1, binEta); - TComplex Q(allReQ, allImQ); - TComplex QStar = TComplex::Conjugate(Q); + int nA = 0, nB = 0; - double posEtaRe = mReQthisEvt->Integral(1, binPt, mReQthisEvt->GetYaxis()->FindBin(fEtaGap + 1e-6), binEta); - double posEtaIm = mImQthisEvt->Integral(1, binPt, mImQthisEvt->GetYaxis()->FindBin(fEtaGap + 1e-6), binEta); - if (mMQthisEvt->Integral(1, binPt, mMQthisEvt->GetYaxis()->FindBin(fEtaGap + 1e-6), binEta) < 2) - return; - float posEtaMQ = mMQWeightthisEvt->Integral(1, binPt, mMQthisEvt->GetYaxis()->FindBin(fEtaGap + 1e-6), binEta); - TComplex posEtaQ = TComplex(posEtaRe, posEtaIm); - TComplex posEtaQStar = TComplex::Conjugate(posEtaQ); + for (auto const& trk : tracks) { + const double pt = trk.pt(); + const double eta = trk.eta(); + if (pt < ptMin || pt > ptMax) { + continue; + } + + const double w = 1.0; // TODO: NUA/NUE weight if you want + const double phi = trk.phi(); + const double c = w * TMath::Cos(harmonic * phi); + const double s = w * TMath::Sin(harmonic * phi); + + if (eta > fEtaGap) { + Q2A_re += c; + Q2A_im += s; + WA += w; + nA++; + } else if (eta < -1 * fEtaGap) { + Q2B_re += c; + Q2B_im += s; + WB += w; + nB++; + } + } - double negEtaRe = mReQthisEvt->Integral(1, binPt, 1, mReQthisEvt->GetYaxis()->FindBin(-1 * fEtaGap - 1e-6)); - double negEtaIm = mImQthisEvt->Integral(1, binPt, 1, mImQthisEvt->GetYaxis()->FindBin(-1 * fEtaGap - 1e-6)); - if (mMQthisEvt->Integral(1, binPt, 1, mMQthisEvt->GetYaxis()->FindBin(-1 * fEtaGap - 1e-6)) < 2) + // need at least 1 track on each side to form pairs; for stability, require >=2 + if (nA < 2 || nB < 2) { return; - float negEtaMQ = mMQWeightthisEvt->Integral(1, binPt, 1, mMQthisEvt->GetYaxis()->FindBin(-1 * fEtaGap - 1e-6)); - TComplex negEtaQ = TComplex(negEtaRe, negEtaIm); - TComplex negEtaQStar = TComplex::Conjugate(negEtaQ); + } + const double D2_evt = WA * WB; + if (D2_evt <= 0.) { + return; + } + + // N2_evt = Re(Q2A * conj(Q2B)) = Q2A_re*Q2B_re + Q2A_im*Q2B_im + const double N2_evt = Q2A_re * Q2B_re + Q2A_im * Q2B_im; - mHistogramQn->get(HIST("Event/profileC22"))->Fill(centrality, (negEtaQ * posEtaQStar).Re() / (negEtaMQ * posEtaMQ), (negEtaMQ * posEtaMQ)); + mHistogramQn->fill(HIST("Event/hN2allQn"), centrality, N2_evt); + mHistogramQn->fill(HIST("Event/hD2allQn"), centrality, D2_evt); if (doQnSeparation && mQnBin >= 0 && mQnBin < numQnBins) { - std::get>(profilesC22[mQnBin])->Fill(centrality, (negEtaQ * posEtaQStar).Re() / (negEtaMQ * posEtaMQ), (negEtaMQ * posEtaMQ)); + std::get>(hN2[mQnBin])->Fill(centrality, N2_evt); + std::get>(hD2[mQnBin])->Fill(centrality, D2_evt); } return; } @@ -516,13 +499,8 @@ class FemtoDreamCollisionSelection float mSphericityPtmin = 0.f; int mQnBin = -999; HistogramRegistry* mHistogramQn = nullptr; ///< For flow cumulant output - 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 - TH2D* mReQ2thisEvt = nullptr; ///< For flow cumulant in an event - TH2D* mImQ2thisEvt = nullptr; ///< For flow cumulant in an event - TH2D* mMQthisEvt = nullptr; ///< For flow cumulant in an event - TH2D* mMQWeightthisEvt = nullptr; ///< For flow cumulant in an event + std::vector hN2; ///< Histograms of c22 per Qn bin + std::vector hD2; ///< Histograms of c22 per Qn bin }; } // namespace o2::analysis::femtoDream