Skip to content

[PWGCF] Fix flow code in femtoDream #14708

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wenyaCern wants to merge 3 commits into
base: master
Choose a base branch
from
+57 −79
Open

[PWGCF] Fix flow code in femtoDream #14708

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
1dadae1
change flow code in femtoDream
wenyaCern Jan 30, 2026
c61d2f2
Please consider the following formatting changes
alibuild Jan 30, 2026
ce1cacb
Merge pull request #14 from alibuild/alibot-cleanup-14708
wenyaCern Jan 30, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
136 changes: 57 additions & 79 deletions PWGCF/FemtoDream/Core/femtoDreamCollisionSelection.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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<TProfile>("Event/profileC22", "; cent; c22", kTProfile, {{10, 0, 100}}, "s");
mHistogramQn->add<TProfile>("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<TH1>(HIST("Event/hN2allQn"))->Sumw2();
mHistogramQn->get<TH1>(HIST("Event/hD2allQn"))->Sumw2();

if (doQnSeparation) {
for (int iqn(0); iqn < mumQnBins; ++iqn) {
profilesC22.push_back(mHistogramQn->add<TProfile>(("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<std::shared_ptr<TH1>>(hN2[iqn])->Sumw2();
std::get<std::shared_ptr<TH1>>(hD2[iqn])->Sumw2();
}
}
return;
Expand Down Expand Up @@ -421,84 +422,66 @@ 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 <typename T1, typename T2>
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
/// \param doQnSeparation to fill flow in divied qn bins
/// \param qnBin should be <int> in 0-9
/// \param fEtaGap eta gap for flow cumulant
template <typename T1, typename T2>
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<TProfile>(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<std::shared_ptr<TProfile>>(profilesC22[mQnBin])->Fill(centrality, (negEtaQ * posEtaQStar).Re() / (negEtaMQ * posEtaMQ), (negEtaMQ * posEtaMQ));
std::get<std::shared_ptr<TH1>>(hN2[mQnBin])->Fill(centrality, N2_evt);
std::get<std::shared_ptr<TH1>>(hD2[mQnBin])->Fill(centrality, D2_evt);
}
return;
}
Expand All @@ -516,13 +499,8 @@ class FemtoDreamCollisionSelection
float mSphericityPtmin = 0.f;
int mQnBin = -999;
HistogramRegistry* mHistogramQn = nullptr; ///< For flow cumulant output
std::vector<HistPtr> 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<HistPtr> hN2; ///< Histograms of c22 per Qn bin
std::vector<HistPtr> hD2; ///< Histograms of c22 per Qn bin
};
} // namespace o2::analysis::femtoDream

Expand Down
Loading