From 0e8194a9ef33b0cd40cf31b651d5163356ccefc4 Mon Sep 17 00:00:00 2001 From: Emil Gorm Nielsen Date: Fri, 28 Aug 2026 15:23:50 +0200 Subject: [PATCH 1/4] Generalise user axis and use for event classifier --- .../Tasks/twoParticleCorrelationsMpi.cxx | 176 ++++++++++++++++-- 1 file changed, 157 insertions(+), 19 deletions(-) diff --git a/PWGCF/TwoParticleCorrelations/Tasks/twoParticleCorrelationsMpi.cxx b/PWGCF/TwoParticleCorrelations/Tasks/twoParticleCorrelationsMpi.cxx index 7406b6521bb..67bebeba4e9 100644 --- a/PWGCF/TwoParticleCorrelations/Tasks/twoParticleCorrelationsMpi.cxx +++ b/PWGCF/TwoParticleCorrelations/Tasks/twoParticleCorrelationsMpi.cxx @@ -124,7 +124,7 @@ struct TwoParticleCorrelationsMpi { Configurable cfgDecayParticleMask{"cfgDecayParticleMask", 0, "Selection bitmask for the decay particles: 0 = no selection"}; Configurable cfgV0RapidityMax{"cfgV0RapidityMax", 0.8, "Maximum rapidity for the decay particles (0 = no selection)"}; - Configurable cfgMassAxis{"cfgMassAxis", 0, "Use invariant mass axis (0 = OFF, 1 = ON)"}; + Configurable cfgUserAxis{"cfgUserAxis", 0, "Additional user axis: 0 = OFF, 1 = invariant mass, 2 = event seed"}; Configurable> cfgMcTriggerPDGs{"cfgMcTriggerPDGs", {}, "MC PDG codes to use exclusively as trigger particles and exclude from associated particles. Empty = no selection."}; ConfigurableAxis axisVertex{"axisVertex", {7, -7, 7}, "vertex axis for histograms"}; @@ -138,7 +138,7 @@ struct TwoParticleCorrelationsMpi { ConfigurableAxis axisEtaEfficiency{"axisEtaEfficiency", {20, -1.0, 1.0}, "eta axis for efficiency histograms"}; ConfigurableAxis axisPtEfficiency{"axisPtEfficiency", {VARIABLE_WIDTH, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.25, 1.5, 1.75, 2.0, 2.25, 2.5, 2.75, 3.0, 3.25, 3.5, 3.75, 4.0, 4.5, 5.0, 6.0, 7.0, 8.0}, "pt axis for efficiency histograms"}; - ConfigurableAxis axisInvMass{"axisInvMass", {VARIABLE_WIDTH, 1.7, 1.75, 1.8, 1.85, 1.9, 1.95, 2.0, 5.0}, "invariant mass axis for histograms"}; + ConfigurableAxis axisUser{"axisUser", {VARIABLE_WIDTH, 1.7, 1.75, 1.8, 1.85, 1.9, 1.95, 2.0, 5.0}, "additional user axis for histograms"}; ConfigurableAxis axisMultCorrCent{"axisMultCorrCent", {100, 0, 100}, "multiplicity correlation axis for centralities"}; ConfigurableAxis axisMultCorrV0{"axisMultCorrV0", {1000, 0, 100000}, "multiplicity correlation axis for V0 multiplicities"}; @@ -237,6 +237,23 @@ struct TwoParticleCorrelationsMpi { } }; + struct PendingSeedTriggerFill { + float pt; + float multiplicity; + float posZ; + float weight; + }; + + struct PendingSeedPairFill { + float deltaEta; + float assocPt; + float triggerPt; + float multiplicity; + float deltaPhi; + float posZ; + float weight; + }; + std::vector yieldTemplates; std::vector> pairAcceptanceMaps; std::vector> pairAcceptanceEtaVertexMaps; @@ -249,6 +266,11 @@ struct TwoParticleCorrelationsMpi { Dd, Ddbar }; + enum UserAxisMode { + NoUserAxis = 0, + InvariantMassAxis, + EventSeedAxis + }; HistogramRegistry registry{"registry"}; PairCuts mPairCuts; @@ -263,6 +285,9 @@ struct TwoParticleCorrelationsMpi { void init(o2::framework::InitContext&) { + if (cfgUserAxis < NoUserAxis || cfgUserAxis > EventSeedAxis) { + LOGF(fatal, "Unsupported cfgUserAxis=%d; use 0 (off), 1 (invariant mass), or 2 (event seed)", cfgUserAxis.value); + } if (doprocessMCSameDerived && (doprocessSameDerived || doprocessSameDerivedMultSet)) { LOGF(fatal, "processMCSameDerived is mutually exclusive with the reconstructed derived same-event processes because it also fills those outputs"); } @@ -279,6 +304,9 @@ struct TwoParticleCorrelationsMpi { LOGF(fatal, "MAP-EM configuration requires non-negative prior exposure, at least one iteration, and positive tolerance"); } eventSeedEstimatorEnabled = !cfgNuncSeedsTemplateFile.value.empty() || !cfgNuncSeedsTemplate.value.empty(); + if (cfgUserAxis == EventSeedAxis && !eventSeedEstimatorEnabled) { + LOGF(fatal, "cfgUserAxis=2 requires an event-seed template configured through cfgNuncSeedsTemplateFile or cfgNuncSeedsTemplate"); + } LOGF(info, "Event seed estimator histogram booking: %s (local template='%s', CCDB template='%s')", eventSeedEstimatorEnabled ? "enabled" : "disabled", cfgNuncSeedsTemplateFile.value.c_str(), cfgNuncSeedsTemplate.value.c_str()); @@ -443,9 +471,12 @@ struct TwoParticleCorrelationsMpi { std::vector userAxis; std::vector userMixingAxis; - if (cfgMassAxis != 0) { - userAxis.emplace_back(axisInvMass, "m (GeV/c^2)"); - userMixingAxis.emplace_back(axisInvMass, "m (GeV/c^2)"); + if (cfgUserAxis == InvariantMassAxis) { + userAxis.emplace_back(axisUser, "m (GeV/c^2)"); + userMixingAxis.emplace_back(axisUser, "m (GeV/c^2)"); + } else if (cfgUserAxis == EventSeedAxis) { + userAxis.emplace_back(axisUser, "N_{seed}"); + userMixingAxis.emplace_back(axisUser, "N_{seed}"); } same.setObject(new CorrelationContainer("sameEvent", "sameEvent", corrAxis, effAxis, userAxis)); mixed.setObject(new CorrelationContainer("mixedEvent", "mixedEvent", corrAxis, effAxis, userMixingAxis)); @@ -929,7 +960,7 @@ struct TwoParticleCorrelationsMpi { } void addPairProbabilities(EventSeedEstimate& estimate, const YieldTemplate& yieldTemplate, - double multiplicity, double deltaPhi, double deltaEta, double posZ) + double multiplicity, double deltaPhi, double deltaEta, double posZ, bool fillAcceptanceQA = true) { const auto& parameters = yieldTemplate.parameters; const double near = std::max(0.0, evaluateGaussian(deltaPhi, parameters.data()) + evaluateGaussian(deltaPhi, parameters.data() + 3)); @@ -979,7 +1010,9 @@ struct TwoParticleCorrelationsMpi { } estimate.sumAcceptanceWeights += acceptanceWeight; ++estimate.nAcceptanceCorrectedPairs; - registry.fill(HIST("eventSeedAcceptanceWeight"), acceptanceWeight); + if (fillAcceptanceQA) { + registry.fill(HIST("eventSeedAcceptanceWeight"), acceptanceWeight); + } } void finalizeEventSeedEstimate(EventSeedEstimate& estimate) @@ -1196,8 +1229,32 @@ struct TwoParticleCorrelationsMpi { } } + template + void flushSeedAxisFills(TTarget target, const std::vector& triggerFills, const std::vector& pairFills, double eventSeed) + { + for (const auto& fill : triggerFills) { + target->getTriggerHist()->Fill(step, fill.pt, fill.multiplicity, fill.posZ, eventSeed, fill.weight); + } + for (const auto& fill : pairFills) { + target->getPairHist()->Fill(step, fill.deltaEta, fill.assocPt, fill.triggerPt, fill.multiplicity, fill.deltaPhi, fill.posZ, eventSeed, fill.weight); + } + } + + template + double estimateEventSeedWithoutFilling(TTarget target, TTracks& tracks, float multiplicity, float posZ, int magField) + { + EventSeedEstimate estimate; + std::vector discardedTriggerFills; + std::vector discardedPairFills; + discardedTriggerFills.reserve(tracks.size()); + discardedPairFills.reserve(tracks.size() * tracks.size()); + fillCorrelations(target, tracks, tracks, multiplicity, posZ, magField, 1.0f, &estimate, &discardedTriggerFills, &discardedPairFills, -1.f, false, false); + finalizeEventSeedEstimate(estimate); + return estimate.nuncSeeds(); + } + template - void fillCorrelations(TTarget target, TTracks1& tracks1, TTracks2& tracks2, float multiplicity, float posZ, int magField, float eventWeight, EventSeedEstimate* seedEstimate = nullptr) + void fillCorrelations(TTarget target, TTracks1& tracks1, TTracks2& tracks2, float multiplicity, float posZ, int magField, float eventWeight, EventSeedEstimate* seedEstimate = nullptr, std::vector* pendingTriggerFills = nullptr, std::vector* pendingPairFills = nullptr, double eventSeed = -1.0, bool fillEstimatorAcceptanceQA = true, bool fillLoopQA = true) { const float containerMultiplicity = getCorrelationContainerMultiplicity(multiplicity); if (containerMultiplicity < 0.f) { @@ -1257,7 +1314,9 @@ struct TwoParticleCorrelationsMpi { if (t && std::abs(y) > cfgV0RapidityMax) { continue; // V0s are not allowed to be outside the rapidity range } - registry.fill(HIST("yvspt"), y, track1.pt()); + if (fillLoopQA) { + registry.fill(HIST("yvspt"), y, track1.pt()); + } } } @@ -1274,7 +1333,13 @@ struct TwoParticleCorrelationsMpi { } } - if (cfgMassAxis) { + if (cfgUserAxis == EventSeedAxis) { + if (pendingTriggerFills) { + pendingTriggerFills->push_back({track1.pt(), containerMultiplicity, posZ, triggerWeight}); + } else { + target->getTriggerHist()->Fill(step, track1.pt(), containerMultiplicity, posZ, eventSeed, triggerWeight); + } + } else if (cfgUserAxis == InvariantMassAxis) { if constexpr (std::experimental::is_detected::value) { target->getTriggerHist()->Fill(step, track1.pt(), containerMultiplicity, posZ, track1.invMass(), triggerWeight); } else if constexpr (std::experimental::is_detected::value) { @@ -1282,7 +1347,7 @@ struct TwoParticleCorrelationsMpi { // target->getTriggerHist()->Fill(step, track1.pt(), multiplicity, posZ, p->Mass(), triggerWeight); target->getTriggerHist()->Fill(step, track1.pt(), containerMultiplicity, posZ, 1.8, triggerWeight); } else { - LOGF(fatal, "Can not fill mass axis without invMass column. Disable cfgMassAxis."); + LOGF(fatal, "Can not fill invariant-mass user axis without invMass column. Disable cfgUserAxis or select another mode."); } } else { target->getTriggerHist()->Fill(step, track1.pt(), containerMultiplicity, posZ, triggerWeight); @@ -1391,20 +1456,26 @@ struct TwoParticleCorrelationsMpi { if (triggerHasTemplate) { ++seedEstimate->nCandidatePairs; if (const auto* yieldTemplate = findYieldTemplate(multiplicity, track1.pt(), track2.pt())) { - addPairProbabilities(*seedEstimate, *yieldTemplate, multiplicity, deltaPhi, deltaEta, posZ); + addPairProbabilities(*seedEstimate, *yieldTemplate, multiplicity, deltaPhi, deltaEta, posZ, fillEstimatorAcceptanceQA); } else { ++seedEstimate->nPairsWithoutTemplate; } } // last param is the weight - if (cfgMassAxis) { + if (cfgUserAxis == EventSeedAxis) { + if (pendingPairFills) { + pendingPairFills->push_back({deltaEta, track2.pt(), track1.pt(), containerMultiplicity, deltaPhi, posZ, associatedWeight}); + } else { + target->getPairHist()->Fill(step, deltaEta, track2.pt(), track1.pt(), containerMultiplicity, deltaPhi, posZ, eventSeed, associatedWeight); + } + } else if (cfgUserAxis == InvariantMassAxis) { if constexpr (std::experimental::is_detected::value) { target->getPairHist()->Fill(step, deltaEta, track2.pt(), track1.pt(), containerMultiplicity, deltaPhi, posZ, track1.invMass(), associatedWeight); } else if constexpr (std::experimental::is_detected::value) { target->getPairHist()->Fill(step, deltaEta, track2.pt(), track1.pt(), containerMultiplicity, deltaPhi, posZ, 1.8, associatedWeight); // p->Mass() } else { - LOGF(fatal, "Can not fill mass axis without invMass column. Disable cfgMassAxis."); + LOGF(fatal, "Can not fill invariant-mass user axis without invMass column. Disable cfgUserAxis or select another mode."); } } else { target->getPairHist()->Fill(step, deltaEta, track2.pt(), track1.pt(), containerMultiplicity, deltaPhi, posZ, associatedWeight); @@ -1482,8 +1553,19 @@ struct TwoParticleCorrelationsMpi { registry.fill(HIST("eventcount_same"), -2); fillQA(collision, multiplicity, tracks); EventSeedEstimate seedEstimate; - fillCorrelations(same, tracks, tracks, multiplicity, collision.posZ(), getMagneticField(bc.timestamp()), 1.0f, &seedEstimate); + std::vector pendingTriggerFills; + std::vector pendingPairFills; + if (cfgUserAxis == EventSeedAxis) { + pendingTriggerFills.reserve(tracks.size()); + pendingPairFills.reserve(tracks.size() * tracks.size()); + } + fillCorrelations(same, tracks, tracks, multiplicity, collision.posZ(), getMagneticField(bc.timestamp()), 1.0f, &seedEstimate, + cfgUserAxis == EventSeedAxis ? &pendingTriggerFills : nullptr, + cfgUserAxis == EventSeedAxis ? &pendingPairFills : nullptr); finalizeEventSeedEstimate(seedEstimate); + if (cfgUserAxis == EventSeedAxis) { + flushSeedAxisFills(same, pendingTriggerFills, pendingPairFills, seedEstimate.nuncSeeds()); + } fillEventSeedEstimatorQA(multiplicity, seedEstimate); if (trueNMPI) { fillMCValidation(multiplicity, seedEstimate, *trueNMPI); @@ -1511,8 +1593,19 @@ struct TwoParticleCorrelationsMpi { const auto generatedMultiplicity = mcCollision.multiplicity(); fillContainerEvent(same, generatedMultiplicity, CorrelationContainer::kCFStepAll); EventSeedEstimate seedEstimate; - fillCorrelations(same, mcParticles, mcParticles, generatedMultiplicity, mcCollision.posZ(), 0, 1.0f, &seedEstimate); + std::vector pendingTriggerFills; + std::vector pendingPairFills; + if (cfgUserAxis == EventSeedAxis) { + pendingTriggerFills.reserve(mcParticles.size()); + pendingPairFills.reserve(mcParticles.size() * mcParticles.size()); + } + fillCorrelations(same, mcParticles, mcParticles, generatedMultiplicity, mcCollision.posZ(), 0, 1.0f, &seedEstimate, + cfgUserAxis == EventSeedAxis ? &pendingTriggerFills : nullptr, + cfgUserAxis == EventSeedAxis ? &pendingPairFills : nullptr); finalizeEventSeedEstimate(seedEstimate); + if (cfgUserAxis == EventSeedAxis) { + flushSeedAxisFills(same, pendingTriggerFills, pendingPairFills, seedEstimate.nuncSeeds()); + } fillGeneratedMCValidation(generatedMultiplicity, seedEstimate, mcCollision.nMPI()); } PROCESS_SWITCH(TwoParticleCorrelationsMpi, processSameGenMC, "Process generated MC events from derived data and validate against the stored HepMC N MPI", false); @@ -1548,6 +1641,35 @@ struct TwoParticleCorrelationsMpi { const bool fillReco = !(cfgDropStepRECO && hasEfficiency); EventSeedEstimate seedEstimate; + if (cfgUserAxis == EventSeedAxis) { + std::vector pendingTriggerFills; + std::vector pendingPairFills; + pendingTriggerFills.reserve(tracks1.size()); + pendingPairFills.reserve(tracks1.size() * tracks2.size()); + + if (fillReco) { + fillContainerEvent(same, multiplicity, CorrelationContainer::kCFStepReconstructed); + fillCorrelations(same, tracks1, tracks2, multiplicity, collision.posZ(), field, 1.0f, &seedEstimate, &pendingTriggerFills, &pendingPairFills); + finalizeEventSeedEstimate(seedEstimate); + flushSeedAxisFills(same, pendingTriggerFills, pendingPairFills, seedEstimate.nuncSeeds()); + } else if (hasEfficiency) { + fillContainerEvent(same, multiplicity, CorrelationContainer::kCFStepCorrected); + fillCorrelations(same, tracks1, tracks2, multiplicity, collision.posZ(), field, 1.0f, &seedEstimate, &pendingTriggerFills, &pendingPairFills); + finalizeEventSeedEstimate(seedEstimate); + flushSeedAxisFills(same, pendingTriggerFills, pendingPairFills, seedEstimate.nuncSeeds()); + } + + if (fillReco && hasEfficiency) { + fillContainerEvent(same, multiplicity, CorrelationContainer::kCFStepCorrected); + fillCorrelations(same, tracks1, tracks2, multiplicity, collision.posZ(), field, 1.0f, nullptr, nullptr, nullptr, seedEstimate.nuncSeeds()); + } + fillEventSeedEstimatorQA(multiplicity, seedEstimate); + if (trueNMPI) { + fillMCValidation(multiplicity, seedEstimate, *trueNMPI); + } + return; + } + if (fillReco) { fillContainerEvent(same, multiplicity, CorrelationContainer::kCFStepReconstructed); fillCorrelations(same, tracks1, tracks2, multiplicity, collision.posZ(), field, 1.0f, &seedEstimate); @@ -1589,6 +1711,7 @@ struct TwoParticleCorrelationsMpi { SameKindPair pairs{configurableBinning, cfgNumMixedEvents, -1, collisions, tracksTuple, &cache}; // -1 is the number of the bin to skip int skipID = -1; + double triggerEventSeed = -1.0; for (auto it = pairs.begin(); it != pairs.end(); it++) { auto& [collision1, tracks1, collision2, tracks2] = *it; int bin = configurableBinning.getBin({collision1.posZ(), collision1.centRun2V0M()}); @@ -1605,6 +1728,11 @@ struct TwoParticleCorrelationsMpi { skipID = collision1.globalIndex(); continue; } + if (cfgUserAxis == EventSeedAxis) { + auto bc = collision1.bc_as(); + loadCcdbYieldTemplates(bc.timestamp()); + triggerEventSeed = estimateEventSeedWithoutFilling(mixed, tracks1, collision1.centRun2V0M(), collision1.posZ(), getMagneticField(bc.timestamp())); + } } if (!collision2.alias_bit(kINT7) || !collision2.sel7()) { continue; @@ -1616,7 +1744,7 @@ struct TwoParticleCorrelationsMpi { // LOGF(info, "Tracks: %d and %d entries", tracks1.size(), tracks2.size()); - fillCorrelations(mixed, tracks1, tracks2, collision1.centRun2V0M(), collision1.posZ(), getMagneticField(bc.timestamp()), 1.0f / it.currentWindowNeighbours()); + fillCorrelations(mixed, tracks1, tracks2, collision1.centRun2V0M(), collision1.posZ(), getMagneticField(bc.timestamp()), 1.0f / it.currentWindowNeighbours(), nullptr, nullptr, nullptr, triggerEventSeed); } } PROCESS_SWITCH(TwoParticleCorrelationsMpi, processMixedAOD, "Process mixed events on AOD", false); @@ -1644,6 +1772,7 @@ struct TwoParticleCorrelationsMpi { using TB = std::tuple_element - 1, decltype(tracksTuple)>::type; Pair pairs{configurableBinningDerived, cfgNumMixedEvents, -1, collisions, tracksTuple, &cache}; // -1 is the number of the bin to skip + double triggerEventSeed = -1.0; for (auto it = pairs.begin(); it != pairs.end(); it++) { auto& [collision1, tracks1, collision2, tracks2] = *it; float multiplicity = getMultiplicity(collision1); @@ -1666,6 +1795,15 @@ struct TwoParticleCorrelationsMpi { hasEfficiencyMixed = (cfg.mEfficiencyAssociated != nullptr || cfg.mEfficiencyTrigger != nullptr); fillRecoMixed = !(cfgDropStepRECO && hasEfficiencyMixed); + if (cfgUserAxis == EventSeedAxis) { + loadCcdbYieldTemplates(collision1.timestamp()); + if constexpr (std::is_same_v, std::remove_cvref_t>) { + triggerEventSeed = estimateEventSeedWithoutFilling(mixed, tracks1, collision1.multiplicity(), collision1.posZ(), field); + } else { + LOGF(fatal, "Event-seed user axis for mixed events requires the same trigger and associated track table so the trigger event can be estimated independently"); + } + } + if (fillRecoMixed) { fillContainerEvent(mixed, collision1.multiplicity(), CorrelationContainer::kCFStepReconstructed); } @@ -1676,14 +1814,14 @@ struct TwoParticleCorrelationsMpi { registry.fill(HIST("eventcount_mixed"), bin); registry.fill(HIST("trackcount_mixed"), bin, tracks1.size(), tracks2.size()); if (fillRecoMixed) { - fillCorrelations(mixed, tracks1, tracks2, collision1.multiplicity(), collision1.posZ(), field, eventWeight); + fillCorrelations(mixed, tracks1, tracks2, collision1.multiplicity(), collision1.posZ(), field, eventWeight, nullptr, nullptr, nullptr, triggerEventSeed); } if (hasEfficiencyMixed) { if (it.isNewWindow()) { fillContainerEvent(mixed, collision1.multiplicity(), CorrelationContainer::kCFStepCorrected); } - fillCorrelations(mixed, tracks1, tracks2, collision1.multiplicity(), collision1.posZ(), field, eventWeight); + fillCorrelations(mixed, tracks1, tracks2, collision1.multiplicity(), collision1.posZ(), field, eventWeight, nullptr, nullptr, nullptr, triggerEventSeed); } } } From 5e1c145b73549dd3bb72ec8c5a61633183932ba1 Mon Sep 17 00:00:00 2001 From: Emil Gorm Nielsen Date: Mon, 31 Aug 2026 11:25:19 +0200 Subject: [PATCH 2/4] make user axis takes event seed quantiles --- .../Tasks/twoParticleCorrelationsMpi.cxx | 70 +++++++++++++++++-- 1 file changed, 63 insertions(+), 7 deletions(-) diff --git a/PWGCF/TwoParticleCorrelations/Tasks/twoParticleCorrelationsMpi.cxx b/PWGCF/TwoParticleCorrelations/Tasks/twoParticleCorrelationsMpi.cxx index 67bebeba4e9..ead110e4f59 100644 --- a/PWGCF/TwoParticleCorrelations/Tasks/twoParticleCorrelationsMpi.cxx +++ b/PWGCF/TwoParticleCorrelations/Tasks/twoParticleCorrelationsMpi.cxx @@ -117,6 +117,8 @@ struct TwoParticleCorrelationsMpi { Configurable cfgEventSeedPriorExposure{"cfgEventSeedPriorExposure", 1.f, "Gamma-prior strength in equivalent event exposures for MAP-EM"}; Configurable cfgEventSeedEMMaxIterations{"cfgEventSeedEMMaxIterations", 10, "Maximum number of weighted MAP-EM iterations per event"}; Configurable cfgEventSeedEMTolerance{"cfgEventSeedEMTolerance", 1.e-5f, "Relative component-yield convergence tolerance for MAP-EM"}; + Configurable> cfgEventSeedPercentileLower{"cfgEventSeedPercentileLower", {}, "Lower event-seed percentile boundary in each axisMultiplicity bin; empty keeps the raw Seed user axis"}; + Configurable> cfgEventSeedPercentileUpper{"cfgEventSeedPercentileUpper", {}, "Upper event-seed percentile boundary in each axisMultiplicity bin; empty keeps the raw Seed user axis"}; Configurable cfgNumMixedEvents{"cfgNumMixedEvents", 5, "Number of mixed events per event"}; @@ -260,6 +262,7 @@ struct TwoParticleCorrelationsMpi { int pairAcceptanceSchemaVersion = 0; const TList* loadedCcdbYieldTemplateObject = nullptr; bool eventSeedEstimatorEnabled = false; + bool eventSeedPercentileAxisEnabled = false; enum CorrelationMethod { All = 0, @@ -307,6 +310,36 @@ struct TwoParticleCorrelationsMpi { if (cfgUserAxis == EventSeedAxis && !eventSeedEstimatorEnabled) { LOGF(fatal, "cfgUserAxis=2 requires an event-seed template configured through cfgNuncSeedsTemplateFile or cfgNuncSeedsTemplate"); } + const bool hasLowerPercentileBoundaries = !cfgEventSeedPercentileLower->empty(); + const bool hasUpperPercentileBoundaries = !cfgEventSeedPercentileUpper->empty(); + if (hasLowerPercentileBoundaries != hasUpperPercentileBoundaries) { + LOGF(fatal, "Configure both cfgEventSeedPercentileLower and cfgEventSeedPercentileUpper, or leave both empty"); + } + eventSeedPercentileAxisEnabled = hasLowerPercentileBoundaries; + if (eventSeedPercentileAxisEnabled) { + const int nMultiplicityBins = AxisSpec(axisMultiplicity).getNbins(); + if (static_cast(cfgEventSeedPercentileLower->size()) != nMultiplicityBins || + static_cast(cfgEventSeedPercentileUpper->size()) != nMultiplicityBins) { + LOGF(fatal, "Event-seed percentile boundary vectors must each contain exactly %d values, one per axisMultiplicity bin", nMultiplicityBins); + } + for (int multBin = 0; multBin < nMultiplicityBins; ++multBin) { + if (!std::isfinite(cfgEventSeedPercentileLower->at(multBin)) || + !std::isfinite(cfgEventSeedPercentileUpper->at(multBin)) || + cfgEventSeedPercentileLower->at(multBin) >= cfgEventSeedPercentileUpper->at(multBin)) { + LOGF(fatal, "Invalid event-seed percentile boundaries in multiplicity bin %d: lower=%g upper=%g", + multBin, cfgEventSeedPercentileLower->at(multBin), cfgEventSeedPercentileUpper->at(multBin)); + } + } + if (cfgUserAxis == EventSeedAxis) { + const std::vector expectedEdges{-1.5, -0.5, 0.5, 1.5, 2.5}; + const auto& configuredEdges = AxisSpec(axisUser).binEdges; + if (configuredEdges.size() != expectedEdges.size() || + !std::equal(configuredEdges.begin(), configuredEdges.end(), expectedEdges.begin(), + [](double lhs, double rhs) { return std::abs(lhs - rhs) < 1.e-6; })) { + LOGF(fatal, "Percentile-class Seed axis requires axisUser edges {-1.5,-0.5,0.5,1.5,2.5}"); + } + } + } LOGF(info, "Event seed estimator histogram booking: %s (local template='%s', CCDB template='%s')", eventSeedEstimatorEnabled ? "enabled" : "disabled", cfgNuncSeedsTemplateFile.value.c_str(), cfgNuncSeedsTemplate.value.c_str()); @@ -342,6 +375,7 @@ struct TwoParticleCorrelationsMpi { registry.add("eventSeedEstimatorStatus", "event-estimator status", {HistType::kTH1F, {{5, -0.5, 4.5, "status"}}}); registry.add("eventSeedTemplateCoverage", "template coverage vs multiplicity", {HistType::kTH2F, {{100, 0, 100, "multiplicity"}, {102, -0.01, 1.01, "matched/candidate pairs"}}}); registry.add("eventSeedEstimateVsMultiplicity", "event-level seed estimate vs multiplicity", {HistType::kTH2F, {{100, 0, 100, "multiplicity"}, {200, 0, 100, "N_{uncorrelated seeds}"}}}); + registry.add("eventSeedFixedProbabilityVsMultiplicity", "fixed-probability event Seed versus multiplicity;multiplicity;N_{seed}^{probability sum}", {HistType::kTH2F, {axisMultiplicity, {1000, 0, 100}}}); registry.add("eventNearYieldVsMultiplicity", "event-level near yield vs multiplicity", {HistType::kTH2F, {{100, 0, 100, "multiplicity"}, {200, 0, 20, "Y_{near}"}}}); registry.add("eventAwayYieldVsMultiplicity", "event-level away yield vs multiplicity", {HistType::kTH2F, {{100, 0, 100, "multiplicity"}, {200, 0, 20, "Y_{away}"}}}); registry.add("profileEventNTriggers", "mean event trigger count", {HistType::kTProfile, {axisMultiplicity}}); @@ -365,6 +399,7 @@ struct TwoParticleCorrelationsMpi { registry.add("profileEventEMAwayPrior", "mean MAP-EM away-component prior mode", {HistType::kTProfile, {axisMultiplicity}}); registry.add("profileEventEMBaselinePrior", "mean MAP-EM baseline-component prior mode", {HistType::kTProfile, {axisMultiplicity}}); registry.add("eventSeedMAPVsProbability", "MAP-EM versus fixed-probability seed estimate;N_{seed}^{probability sum};N_{seed}^{MAP-EM}", {HistType::kTH2F, {{200, 0, 100}, {200, 0, 100}}}); + registry.add("eventSeedMAPEMVsMultiplicity", "MAP-EM event Seed versus multiplicity;multiplicity;N_{seed}^{MAP-EM}", {HistType::kTH2F, {axisMultiplicity, {1000, 0, 100}}}); registry.add("profileEventProbabilityNuncSeeds", "mean fixed-probability seed estimate", {HistType::kTProfile, {axisMultiplicity}}); } auto* estimatorStatus = registry.get(HIST("eventSeedEstimatorStatus")).get(); @@ -475,8 +510,9 @@ struct TwoParticleCorrelationsMpi { userAxis.emplace_back(axisUser, "m (GeV/c^2)"); userMixingAxis.emplace_back(axisUser, "m (GeV/c^2)"); } else if (cfgUserAxis == EventSeedAxis) { - userAxis.emplace_back(axisUser, "N_{seed}"); - userMixingAxis.emplace_back(axisUser, "N_{seed}"); + const char* title = eventSeedPercentileAxisEnabled ? "N_{seed} percentile class" : "N_{seed}"; + userAxis.emplace_back(axisUser, title); + userMixingAxis.emplace_back(axisUser, title); } same.setObject(new CorrelationContainer("sameEvent", "sameEvent", corrAxis, effAxis, userAxis)); mixed.setObject(new CorrelationContainer("mixedEvent", "mixedEvent", corrAxis, effAxis, userMixingAxis)); @@ -929,6 +965,24 @@ struct TwoParticleCorrelationsMpi { return static_cast(std::distance(edges.begin(), upper)) - 1; } + double getEventSeedUserAxisValue(double multiplicity, double eventSeed) const + { + if (!eventSeedPercentileAxisEnabled || eventSeed < 0.0) { + return eventSeed; + } + const int multBin = findPairAcceptanceMultiplicityBin(multiplicity); + if (multBin < 0 || multBin >= static_cast(cfgEventSeedPercentileLower->size())) { + return -1.0; + } + if (eventSeed < cfgEventSeedPercentileLower->at(multBin)) { + return 0.0; + } + if (eventSeed >= cfgEventSeedPercentileUpper->at(multBin)) { + return 2.0; + } + return 1.0; + } + double getPairAcceptance(double multiplicity, double deltaPhi, double deltaEta, double posZ) const { const int multBin = findPairAcceptanceMultiplicityBin(multiplicity); @@ -1125,6 +1179,7 @@ struct TwoParticleCorrelationsMpi { registry.fill(HIST("eventSeedEstimator"), multiplicity, estimate.nTriggers, estimate.nearYield(), estimate.awayYield(), estimate.nuncSeeds()); registry.fill(HIST("eventSeedPairProbabilities"), estimate.probabilityBaselinePairs, estimate.probabilityNearPairs, estimate.probabilityAwayPairs); registry.fill(HIST("eventSeedEstimateVsMultiplicity"), multiplicity, estimate.nuncSeeds()); + registry.fill(HIST("eventSeedFixedProbabilityVsMultiplicity"), multiplicity, estimate.probabilityNuncSeeds()); registry.fill(HIST("eventNearYieldVsMultiplicity"), multiplicity, estimate.nearYield()); registry.fill(HIST("eventAwayYieldVsMultiplicity"), multiplicity, estimate.awayYield()); registry.fill(HIST("profileEventNearYield"), multiplicity, estimate.nearYield()); @@ -1138,6 +1193,7 @@ struct TwoParticleCorrelationsMpi { registry.fill(HIST("profileEventEMAwayPrior"), multiplicity, estimate.priorComponentCounts[1]); registry.fill(HIST("profileEventEMBaselinePrior"), multiplicity, estimate.priorComponentCounts[2]); registry.fill(HIST("eventSeedMAPVsProbability"), estimate.probabilityNuncSeeds(), estimate.nuncSeeds()); + registry.fill(HIST("eventSeedMAPEMVsMultiplicity"), multiplicity, estimate.nuncSeeds()); } } @@ -1233,10 +1289,10 @@ struct TwoParticleCorrelationsMpi { void flushSeedAxisFills(TTarget target, const std::vector& triggerFills, const std::vector& pairFills, double eventSeed) { for (const auto& fill : triggerFills) { - target->getTriggerHist()->Fill(step, fill.pt, fill.multiplicity, fill.posZ, eventSeed, fill.weight); + target->getTriggerHist()->Fill(step, fill.pt, fill.multiplicity, fill.posZ, getEventSeedUserAxisValue(fill.multiplicity, eventSeed), fill.weight); } for (const auto& fill : pairFills) { - target->getPairHist()->Fill(step, fill.deltaEta, fill.assocPt, fill.triggerPt, fill.multiplicity, fill.deltaPhi, fill.posZ, eventSeed, fill.weight); + target->getPairHist()->Fill(step, fill.deltaEta, fill.assocPt, fill.triggerPt, fill.multiplicity, fill.deltaPhi, fill.posZ, getEventSeedUserAxisValue(fill.multiplicity, eventSeed), fill.weight); } } @@ -1661,7 +1717,7 @@ struct TwoParticleCorrelationsMpi { if (fillReco && hasEfficiency) { fillContainerEvent(same, multiplicity, CorrelationContainer::kCFStepCorrected); - fillCorrelations(same, tracks1, tracks2, multiplicity, collision.posZ(), field, 1.0f, nullptr, nullptr, nullptr, seedEstimate.nuncSeeds()); + fillCorrelations(same, tracks1, tracks2, multiplicity, collision.posZ(), field, 1.0f, nullptr, nullptr, nullptr, getEventSeedUserAxisValue(multiplicity, seedEstimate.nuncSeeds())); } fillEventSeedEstimatorQA(multiplicity, seedEstimate); if (trueNMPI) { @@ -1731,7 +1787,7 @@ struct TwoParticleCorrelationsMpi { if (cfgUserAxis == EventSeedAxis) { auto bc = collision1.bc_as(); loadCcdbYieldTemplates(bc.timestamp()); - triggerEventSeed = estimateEventSeedWithoutFilling(mixed, tracks1, collision1.centRun2V0M(), collision1.posZ(), getMagneticField(bc.timestamp())); + triggerEventSeed = getEventSeedUserAxisValue(collision1.centRun2V0M(), estimateEventSeedWithoutFilling(mixed, tracks1, collision1.centRun2V0M(), collision1.posZ(), getMagneticField(bc.timestamp()))); } } if (!collision2.alias_bit(kINT7) || !collision2.sel7()) { @@ -1798,7 +1854,7 @@ struct TwoParticleCorrelationsMpi { if (cfgUserAxis == EventSeedAxis) { loadCcdbYieldTemplates(collision1.timestamp()); if constexpr (std::is_same_v, std::remove_cvref_t>) { - triggerEventSeed = estimateEventSeedWithoutFilling(mixed, tracks1, collision1.multiplicity(), collision1.posZ(), field); + triggerEventSeed = getEventSeedUserAxisValue(collision1.multiplicity(), estimateEventSeedWithoutFilling(mixed, tracks1, collision1.multiplicity(), collision1.posZ(), field)); } else { LOGF(fatal, "Event-seed user axis for mixed events requires the same trigger and associated track table so the trigger event can be estimated independently"); } From cc520b54000cf81d2c0715ed5ba0238ae301b384 Mon Sep 17 00:00:00 2001 From: Emil Gorm Nielsen Date: Mon, 31 Aug 2026 11:31:40 +0200 Subject: [PATCH 3/4] linter --- .../Tasks/twoParticleCorrelationsMpi.cxx | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/PWGCF/TwoParticleCorrelations/Tasks/twoParticleCorrelationsMpi.cxx b/PWGCF/TwoParticleCorrelations/Tasks/twoParticleCorrelationsMpi.cxx index ead110e4f59..ce25dcd8a12 100644 --- a/PWGCF/TwoParticleCorrelations/Tasks/twoParticleCorrelationsMpi.cxx +++ b/PWGCF/TwoParticleCorrelations/Tasks/twoParticleCorrelationsMpi.cxx @@ -333,9 +333,10 @@ struct TwoParticleCorrelationsMpi { if (cfgUserAxis == EventSeedAxis) { const std::vector expectedEdges{-1.5, -0.5, 0.5, 1.5, 2.5}; const auto& configuredEdges = AxisSpec(axisUser).binEdges; + const float threshold = 1.e-6; if (configuredEdges.size() != expectedEdges.size() || !std::equal(configuredEdges.begin(), configuredEdges.end(), expectedEdges.begin(), - [](double lhs, double rhs) { return std::abs(lhs - rhs) < 1.e-6; })) { + [](double lhs, double rhs) { return std::abs(lhs - rhs) < threshold; })) { LOGF(fatal, "Percentile-class Seed axis requires axisUser edges {-1.5,-0.5,0.5,1.5,2.5}"); } } From 95174f2c3e7cd7255bdfba8570c41829658c26ec Mon Sep 17 00:00:00 2001 From: Emil Gorm Nielsen Date: Mon, 31 Aug 2026 11:47:47 +0200 Subject: [PATCH 4/4] fix lambda capture --- .../Tasks/twoParticleCorrelationsMpi.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PWGCF/TwoParticleCorrelations/Tasks/twoParticleCorrelationsMpi.cxx b/PWGCF/TwoParticleCorrelations/Tasks/twoParticleCorrelationsMpi.cxx index ce25dcd8a12..908ce848d7d 100644 --- a/PWGCF/TwoParticleCorrelations/Tasks/twoParticleCorrelationsMpi.cxx +++ b/PWGCF/TwoParticleCorrelations/Tasks/twoParticleCorrelationsMpi.cxx @@ -336,7 +336,7 @@ struct TwoParticleCorrelationsMpi { const float threshold = 1.e-6; if (configuredEdges.size() != expectedEdges.size() || !std::equal(configuredEdges.begin(), configuredEdges.end(), expectedEdges.begin(), - [](double lhs, double rhs) { return std::abs(lhs - rhs) < threshold; })) { + [&threshold](double lhs, double rhs) { return std::abs(lhs - rhs) < threshold; })) { LOGF(fatal, "Percentile-class Seed axis requires axisUser edges {-1.5,-0.5,0.5,1.5,2.5}"); } }