Skip to content

[PWGCF] Add MC gen process #16765

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 4 commits into
base: master
Choose a base branch
from
+142 −33
Open

[PWGCF] Add MC gen process #16765

Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
26fc012
Add cutoff of eta for q-vec and MC process
wenyaCern Jun 22, 2026
28c53e5
Add cutoff of eta for q-vec and MC process
wenyaCern Jun 22, 2026
56f4d8a
Please consider the following formatting changes
alibuild Jun 22, 2026
edb8a3c
Merge pull request #18 from alibuild/alibot-cleanup-16765
wenyaCern Jun 22, 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
175 changes: 142 additions & 33 deletions PWGCF/Flow/Tasks/flowFlucGfwPp.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -121,7 +121,7 @@

O2_DEFINE_CONFIGURABLE(cfgNbootstrap, int, 10, "Number of subsamples")
O2_DEFINE_CONFIGURABLE(cfgIsMC, bool, false, "Is MC event")
O2_DEFINE_CONFIGURABLE(cfgCentEstimator, int, 0, "0:FT0C; 1:FT0CVariant1; 2:FT0M; 3:FV0A, 4:NTPV, 5:NGlobals, 6:MFT")
O2_DEFINE_CONFIGURABLE(cfgCentEstimator, int, 0, "0:FT0C; 1:FT0CVariant1; 2:FT0M; 3:FV0A, 4:NTPV, 5:NGlobals")
O2_DEFINE_CONFIGURABLE(cfgUseNch, bool, false, "Do correlations as function of Nch")
O2_DEFINE_CONFIGURABLE(cfgQvecQA, bool, false, "Enable filling QA for q-Vec of TPC")
O2_DEFINE_CONFIGURABLE(cfgFillWeights, bool, false, "Fill NUA weights")
Expand Down Expand Up @@ -171,6 +171,7 @@
O2_DEFINE_CONFIGURABLE(cfgUseNegativeEtaHalfForq2, bool, true, "If true, use -eta half for qn selection; otherwise use +eta half");
O2_DEFINE_CONFIGURABLE(cfgQnSelectionHarmonic, int, 2, "Harmonic n used to build the reduced q_n vector for event shape selection, use 2 for q2 and 3 for q3");
O2_DEFINE_CONFIGURABLE(cfgQnHistMax, float, 6., "Upper range for q_n calibration histograms");
O2_DEFINE_CONFIGURABLE(cfgQnTrkAbsEtaMax, float, 0.5, "Upper range for abs eta of tracks contributing to q_n");
O2_DEFINE_CONFIGURABLE(cfgBypassQnSelection, bool, false, "Bypass q_n event shape selection and fill one integral q-bin");
O2_DEFINE_CONFIGURABLE(cfgMinPtOnTPC, float, 0.2, "minimum transverse momentum selection for TPC tracks participating in Q-vector reconstruction");
O2_DEFINE_CONFIGURABLE(cfgMaxPtOnTPC, float, 5., "maximum transverse momentum selection for TPC tracks participating in Q-vector reconstruction");
Expand Down Expand Up @@ -277,8 +278,7 @@
kCentFT0M,
kCentFV0A,
kCentNTPV,
kCentNGlobal,
kCentMFT
kCentNGlobal
};
enum EventSelFlags {
kFilteredEvent = 1,
Expand Down Expand Up @@ -337,11 +337,11 @@
o2::framework::expressions::Filter collisionFilter = nabs(aod::collision::posZ) < cfgVtxZ;
o2::framework::expressions::Filter trackFilter = nabs(aod::track::eta) < cfgEta && aod::track::pt > cfgPtmin&& aod::track::pt < cfgPtmax && (aod::track::itsChi2NCl < cfgChi2PrITSCls) && (aod::track::tpcChi2NCl < cfgChi2PrTPCCls) && nabs(aod::track::dcaZ) < cfgDCAz;

Preslice<aod::Tracks> perCollision = aod::track::collisionId;
o2::framework::expressions::Filter mcCollFilter = nabs(aod::mccollision::posZ) < cfgVtxZ;
o2::framework::expressions::Filter mcParticlesFilter = (aod::mcparticle::eta > o2::analysis::gfwflowflucpp::etalow && aod::mcparticle::eta < o2::analysis::gfwflowflucpp::etaup && aod::mcparticle::pt > o2::analysis::gfwflowflucpp::ptlow && aod::mcparticle::pt < o2::analysis::gfwflowflucpp::ptup);

using GFWTracks = soa::Filtered<soa::Join<aod::Tracks, aod::TracksExtra, aod::TrackSelection, aod::TrackSelectionExtension, aod::TracksDCA>>;
using GFWTracksMC = soa::Filtered<soa::Join<aod::Tracks, aod::TracksExtra, aod::TrackSelection, aod::TrackSelectionExtension, aod::TracksDCA, aod::McTrackLabels>>;

void init(InitContext const&)
{
Expand Down Expand Up @@ -415,8 +415,7 @@
{kCentFT0M, "FT0M"},
{kCentFV0A, "FV0A"},
{kCentNTPV, "NTPV"},
{kCentNGlobal, "NGlobals"},
{kCentMFT, "MFT"}};
{kCentNGlobal, "NGlobals"}};
sCentralityEstimator = centEstimatorMap.at(cfgCentEstimator);
sCentralityEstimator += " centrality (%)";
AxisSpec centAxis = {o2::analysis::gfwflowflucpp::centbinning, sCentralityEstimator.c_str()};
Expand Down Expand Up @@ -503,7 +502,7 @@
}
}

if (doprocessData) {
if (doprocessData || doprocessMC) {
registry.add("trackQA/before/phi_eta_vtxZ", "", {HistType::kTH3D, {phiAxis, etaAxis, vtxAxis}});
registry.add("trackQA/before/pt_dcaXY_dcaZ", "", {HistType::kTH3D, {ptAxis, dcaXYAXis, dcaZAXis}});
registry.add("trackQA/before/nch_pt", "#it{p}_{T} vs multiplicity; N_{ch}; #it{p}_{T}", {HistType::kTH2D, {nchAxis, ptAxis}});
Expand All @@ -530,24 +529,20 @@
registry.add("eventQA/before/globalTracks_multV0A", "", {HistType::kTH2D, {v0aAxis, nchAxis}});
registry.add("eventQA/before/multV0A_multT0A", "", {HistType::kTH2D, {t0aAxis, v0aAxis}});

if (doprocessData) {
registry.add("eventQA/before/centrality", "", {HistType::kTH1D, {centAxis}});
registry.add("eventQA/before/globalTracks_centT0C", "", {HistType::kTH2D, {centAxis, nchAxis}});
registry.add("eventQA/before/PVTracks_centT0C", "", {HistType::kTH2D, {centAxis, multpvAxis}});
registry.add("eventQA/before/multT0C_centT0C", "", {HistType::kTH2D, {centAxis, t0cAxis}});

registry.add("eventQA/before/centT0M_centT0C", "", {HistType::kTH2D, {centAxis, centAxis}});
registry.add("eventQA/before/centV0A_centT0C", "", {HistType::kTH2D, {centAxis, centAxis}});
registry.add("eventQA/before/centGlobal_centT0C", "", {HistType::kTH2D, {centAxis, centAxis}});
registry.add("eventQA/before/centNTPV_centT0C", "", {HistType::kTH2D, {centAxis, centAxis}});
registry.add("eventQA/before/centMFT_centT0C", "", {HistType::kTH2D, {centAxis, centAxis}});

if (cfgIsMC) {
registry.add("MCGen/trackQA/phi_eta_vtxZ", "", {HistType::kTH3D, {phiAxis, etaAxis, vtxAxis}});
registry.add("MCGen/trackQA/nch_pt", "#it{p}_{T} vs multiplicity; N_{ch}; #it{p}_{T}", {HistType::kTH2D, {nchAxis, ptAxis}});
registry.add("MCGen/trackQA/pt_ref", "", {HistType::kTH1D, {{100, o2::analysis::gfwflowflucpp::ptreflow, o2::analysis::gfwflowflucpp::ptrefup}}});
registry.add("MCGen/trackQA/pt_poi", "", {HistType::kTH1D, {{100, o2::analysis::gfwflowflucpp::ptpoilow, o2::analysis::gfwflowflucpp::ptpoiup}}});
}
registry.add("eventQA/before/centrality", "", {HistType::kTH1D, {centAxis}});
registry.add("eventQA/before/globalTracks_centT0C", "", {HistType::kTH2D, {centAxis, nchAxis}});
registry.add("eventQA/before/PVTracks_centT0C", "", {HistType::kTH2D, {centAxis, multpvAxis}});
registry.add("eventQA/before/multT0C_centT0C", "", {HistType::kTH2D, {centAxis, t0cAxis}});

registry.add("eventQA/before/centT0M_centT0C", "", {HistType::kTH2D, {centAxis, centAxis}});
registry.add("eventQA/before/centV0A_centT0C", "", {HistType::kTH2D, {centAxis, centAxis}});
registry.add("eventQA/before/centGlobal_centT0C", "", {HistType::kTH2D, {centAxis, centAxis}});
registry.add("eventQA/before/centNTPV_centT0C", "", {HistType::kTH2D, {centAxis, centAxis}});
if (cfgIsMC || doprocessMC) {
registry.add("MCGen/trackQA/phi_eta_vtxZ", "", {HistType::kTH3D, {phiAxis, etaAxis, vtxAxis}});
registry.add("MCGen/trackQA/nch_pt", "#it{p}_{T} vs multiplicity; N_{ch}; #it{p}_{T}", {HistType::kTH2D, {nchAxis, ptAxis}});
registry.add("MCGen/trackQA/pt_ref", "", {HistType::kTH1D, {{100, o2::analysis::gfwflowflucpp::ptreflow, o2::analysis::gfwflowflucpp::ptrefup}}});
registry.add("MCGen/trackQA/pt_poi", "", {HistType::kTH1D, {{100, o2::analysis::gfwflowflucpp::ptpoilow, o2::analysis::gfwflowflucpp::ptpoiup}}});
}

registry.addClone("eventQA/before/", "eventQA/after/");
Expand Down Expand Up @@ -1020,14 +1015,14 @@
registry.fill(HIST("qvecQA/ChTracks"), trk.pt(), trk.eta(), trk.phi());
}

if (trk.eta() > 0) {
if (trk.eta() > 0 && fabs(trk.eta()) < cfgQnTrkAbsEtaMax) {

Check failure on line 1018 in PWGCF/Flow/Tasks/flowFlucGfwPp.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[std-prefix] 

Use std:: prefix for names from the std namespace.
// In qVectorsTable this branch is additionally guarded by useDetector["QvectorTPCposs"] || useDetector["QvectorBPoss"].
// Here TPCpos is always computed because the downstream ESE selector can require it.
qvec.qVectTPCPos[0] += trk.pt() * std::cos(trk.phi() * nMode);
qvec.qVectTPCPos[1] += trk.pt() * std::sin(trk.phi() * nMode);
qvec.trkTPCPosLabel.push_back(trk.globalIndex());
qvec.nTrkTPCPos++;
} else if (trk.eta() < 0) {
} else if (trk.eta() < 0 && fabs(trk.eta()) < cfgQnTrkAbsEtaMax) {

Check failure on line 1025 in PWGCF/Flow/Tasks/flowFlucGfwPp.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[std-prefix] 

Use std:: prefix for names from the std namespace.
// In qVectorsTable this branch is additionally guarded by useDetector["QvectorTPCnegs"] || useDetector["QvectorBNegs"].
// Here TPCneg is always computed because the downstream ESE selector can require it.
qvec.qVectTPCNeg[0] += trk.pt() * std::cos(trk.phi() * nMode);
Expand Down Expand Up @@ -1172,6 +1167,48 @@
lRandom, qPtmp, run);
}

template <typename TCollision, typename TParticles>
void processGenCollision(TCollision collision, TParticles particles, const int& mcCollisionId, const XAxis& xaxis, const int& run, const int& qPtmp)
{
if (xaxis.multiplicity < cfgFixedMultMin || xaxis.multiplicity > cfgFixedMultMax)
return;

if (cfgFillQA && xaxis.centrality >= 0)
registry.fill(HIST("eventQA/after/centrality"), xaxis.centrality);
if (cfgFillQA)
registry.fill(HIST("eventQA/after/multiplicity"), xaxis.multiplicity);

fGFW->Clear();
float lRandom = fRndm->Rndm();
float vtxz = collision.posZ();

AcceptedTracks acceptedTracks{0, 0, 0, 0};
for (const auto& particle : particles) {
if (particle.mcCollisionId() != mcCollisionId)
continue;
processTrack(particle, vtxz, xaxis.multiplicity, run, acceptedTracks);
}

if (cfgConsistentEventFlag & kRequireBothEtaSides)
if (!acceptedTracks.nPos || !acceptedTracks.nNeg)
return;
if (cfgConsistentEventFlag & kRequireFullFourParticleTracks)
if (acceptedTracks.nFull < kMinTracksForFourParticleCorrelation)
return;
if (cfgConsistentEventFlag & kRequireTwoTracksInBothEtaSides)
if (acceptedTracks.nPos < kMinTracksPerEtaSideForGapCorrelation ||
acceptedTracks.nNeg < kMinTracksPerEtaSideForGapCorrelation)
return;
if (cfgConsistentEventFlag & kRequireTwoTracksInThreeEtaRegions)
if (acceptedTracks.nPos < kMinTracksPerEtaRegionForThreeSubevents ||
acceptedTracks.nMid < kMinTracksPerEtaRegionForThreeSubevents ||
acceptedTracks.nNeg < kMinTracksPerEtaRegionForThreeSubevents)
return;

fillOutputContainers<kGen>(cfgUseNch ? static_cast<float>(xaxis.multiplicity) : xaxis.centrality,
lRandom, qPtmp, run);
}

bool isStable(int pdg)
{
if (std::abs(pdg) == PDG_t::kPiPlus)
Expand Down Expand Up @@ -1334,8 +1371,6 @@
return collision.centNTPV();
case kCentNGlobal:
return collision.centNGlobal();
case kCentMFT:
return collision.centMFT();
default:
return collision.centFT0C();
}
Expand All @@ -1352,7 +1387,6 @@
registry.fill(HIST("eventQA/") + HIST(FillTimeName[ft]) + HIST("centV0A_centT0C"), collision.centFT0C(), collision.centFV0A());
registry.fill(HIST("eventQA/") + HIST(FillTimeName[ft]) + HIST("centGlobal_centT0C"), collision.centFT0C(), collision.centNGlobal());
registry.fill(HIST("eventQA/") + HIST(FillTimeName[ft]) + HIST("centNTPV_centT0C"), collision.centFT0C(), collision.centNTPV());
registry.fill(HIST("eventQA/") + HIST(FillTimeName[ft]) + HIST("centMFT_centT0C"), collision.centFT0C(), collision.centMFT());
}
registry.fill(HIST("eventQA/") + HIST(FillTimeName[ft]) + HIST("globalTracks_PVTracks"), collision.multNTracksPV(), xaxis.multiplicity);
registry.fill(HIST("eventQA/") + HIST(FillTimeName[ft]) + HIST("globalTracks_multT0A"), collision.multFT0A(), xaxis.multiplicity);
Expand Down Expand Up @@ -1406,8 +1440,7 @@

void processData(soa::Filtered<soa::Join<aod::Collisions, aod::EvSels, aod::Mults,
aod::CentFT0Cs, aod::CentFT0CVariant1s, aod::CentFT0Ms,
aod::CentFV0As, aod::CentNTPVs, aod::CentNGlobals,
aod::CentMFTs>>::iterator const& collision,
aod::CentFV0As, aod::CentNTPVs, aod::CentNGlobals>>::iterator const& collision,
aod::BCsWithTimestamps const&, GFWTracks const& tracks)
{
auto bc = collision.bc_as<aod::BCsWithTimestamps>();
Expand Down Expand Up @@ -1487,7 +1520,7 @@
}
PROCESS_SWITCH(FlowFlucGfwPp, processData, "Process analysis for non-derived data", false);

void processq2(soa::Filtered<soa::Join<aod::Collisions, aod::EvSels, aod::Mults, aod::CentFT0Cs, aod::CentFT0CVariant1s, aod::CentFT0Ms, aod::CentFV0As, aod::CentNTPVs, aod::CentNGlobals, aod::CentMFTs>>::iterator const& collision, aod::BCsWithTimestamps const&, GFWTracks const& tracks)
void processq2(soa::Filtered<soa::Join<aod::Collisions, aod::EvSels, aod::Mults, aod::CentFT0Cs, aod::CentFT0CVariant1s, aod::CentFT0Ms, aod::CentFV0As, aod::CentNTPVs, aod::CentNGlobals>>::iterator const& collision, aod::BCsWithTimestamps const&, GFWTracks const& tracks)
{
float count{0.5};
fillQnEventCounter(count++);
Expand Down Expand Up @@ -1521,6 +1554,82 @@
fillQnCalibrationHistograms(centr, multi, qvecPos, qvecNeg);
}
PROCESS_SWITCH(FlowFlucGfwPp, processq2, "Process analysis for filling q_n-vector calibration histograms", true);

void processMC(soa::Filtered<soa::Join<aod::Collisions, aod::EvSels, aod::Mults,
aod::CentFT0Cs, aod::CentFT0CVariant1s, aod::CentFT0Ms,
aod::CentFV0As, aod::CentNTPVs, aod::CentNGlobals,
aod::McCollisionLabels>>::iterator const& collision,
aod::BCsWithTimestamps const&, GFWTracksMC const& tracks, aod::McCollisions const&, aod::McParticles const& mcParticles)
{
auto bc = collision.bc_as<aod::BCsWithTimestamps>();
int run = bc.runNumber();
if (run != lastRun) {
lastRun = run;
LOGF(info, "run = %d", run);
if (cfgRunByRun) {
if (std::find(runNumbers.begin(), runNumbers.end(), run) == runNumbers.end()) {
LOGF(info, "Creating histograms for run %d", run);
createRunByRunHistograms(run);
runNumbers.push_back(run);
}
if (!cfgFillWeights)
loadCorrections(bc);
}
}
if (!cfgFillWeights && !cfgRunByRun)
loadCorrections(bc);

registry.fill(HIST("eventQA/eventSel"), kFilteredEvent);
if (cfgRunByRun)
th1sList[run][hEventSel]->Fill(kFilteredEvent);

if (!collision.sel8())
return;

registry.fill(HIST("eventQA/eventSel"), kSel8);
if (cfgRunByRun)
th1sList[run][hEventSel]->Fill(kSel8);

if (cfgDoOccupancySel) {
int occupancy = collision.trackOccupancyInTimeRange();
if (occupancy < 0 || occupancy > cfgOccupancySelection)
return;
}

registry.fill(HIST("eventQA/eventSel"), kOccupancy);
if (cfgRunByRun)
th1sList[run][hEventSel]->Fill(kOccupancy);

const XAxis xaxis{
getCentrality(collision),
collision.multNTracksPV(),
(cfgTimeDependent) ? getTimeSinceStartOfFill(bc.timestamp(), *firstRunOfCurrentFill) : -1.0};

if (cfgTimeDependent && run == *firstRunOfCurrentFill &&
firstRunOfCurrentFill != o2::analysis::gfwflowflucpp::firstRunsOfFill.end() - 1)
++firstRunOfCurrentFill;

if (cfgFillQA)
fillEventQA<kBefore>(collision, xaxis);

registry.fill(HIST("eventQA/before/centrality"), xaxis.centrality);
registry.fill(HIST("eventQA/before/multiplicity"), xaxis.multiplicity);

if (!eventSelected(collision, xaxis.multiplicity, xaxis.centrality, run))
return;

if (cfgFillQA)
fillEventQA<kAfter>(collision, xaxis);

processCollision<kReco>(collision, tracks, xaxis, run, 0);

if (!collision.has_mcCollision())
return;

const auto genCollision = collision.template mcCollision_as<aod::McCollisions>();
processGenCollision(genCollision, mcParticles, collision.mcCollisionId(), xaxis, run, 0);
}
PROCESS_SWITCH(FlowFlucGfwPp, processMC, "Process analysis for Monte-Carlo data", false);
};

WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
Expand Down
Loading