From c1d81ff02b3089cd6e327b06cd157d2cb12960b6 Mon Sep 17 00:00:00 2001 From: Emil Gorm Nielsen Date: Sun, 8 Nov 2026 13:14:24 +0100 Subject: [PATCH 1/4] Add poissonian sampling of efficiency to corrected multiplicity estimate --- PWGCF/TableProducer/filterCorrelations.cxx | 42 ++++++++++++++++++---- 1 file changed, 36 insertions(+), 6 deletions(-) diff --git a/PWGCF/TableProducer/filterCorrelations.cxx b/PWGCF/TableProducer/filterCorrelations.cxx index d31f482fa9b..307a27d9e3d 100644 --- a/PWGCF/TableProducer/filterCorrelations.cxx +++ b/PWGCF/TableProducer/filterCorrelations.cxx @@ -55,6 +55,7 @@ #include // required for is_detected #include #include +#include #include #include #include @@ -120,6 +121,8 @@ struct FilterCF { O2_DEFINE_CONFIGURABLE(cfgEfficiencyMultiplicity, std::string, "", "Multiplicity efficiency (RecoAll / MC): CCDB path or local ROOT file with a 4D ccdb_object (eta, pT, multiplicity, z-vtx); empty disables CFCollisionsExtra output"); O2_DEFINE_CONFIGURABLE(cfgLocalEfficiency, int, 0, "0 = CCDB efficiency, 1 = local ROOT efficiency"); O2_DEFINE_CONFIGURABLE(cfgMultiplicityTrackBitMask, uint16_t, 0, "Required track-type bits for corrected multiplicity; match cfgTrackBitMask used to produce the efficiency (0 = all stored tracks)"); + O2_DEFINE_CONFIGURABLE(cfgMultiplicityPoisson, bool, false, "Use (reproducibly) seeded 1 + Poisson(sum(1/eff - 1)) for multiplicity correction; skips eff > 1 tracks"); + O2_DEFINE_CONFIGURABLE(cfgMultiplicityPoissonSeed, uint32_t, 1, "Base seed for reproducible per-collision Poisson multiplicity sampling"); // Filters and input definitions Filter collisionZVtxFilter = nabs(aod::collision::posZ) < cfgCutVertex; @@ -160,6 +163,7 @@ struct FilterCF { // Own local histograms independently of their input file. CCDB owns its objects. std::unique_ptr localMultiplicityEfficiency; THn* mEfficiency = nullptr; + std::mt19937_64 mMultiplicityGenerator; static constexpr int MultiplicityEfficiencyDimensions = 4; // persistent caches @@ -389,11 +393,12 @@ struct FilterCF { } template - float getCorrectedMultiplicity(const TCollision& collision, const TTracks& tracks, uint64_t timestamp) + float getCorrectedMultiplicity(const TCollision& collision, const TTracks& tracks, uint64_t timestamp, uint64_t globalBC, int runNumber) { auto* efficiency = loadMultiplicityEfficiency(timestamp); double correctedMultiplicity = 0.; size_t skippedTracks = 0; + double totalMeanMissing = 0.; for (const auto& track : tracks) { // Match the tracks written by the corresponding data/MC producer path. if (!isTrackSelected(track, true)) { @@ -403,14 +408,39 @@ struct FilterCF { // Keep the original estimator as the map coordinate, including for centrality. const std::array values{track.eta(), track.pt(), collision.multiplicity(), collision.posZ()}; const double eff = efficiency->GetBinContent(efficiency->GetBin(values.data())); - if (!std::isfinite(eff) || eff <= 0.) { + if (!std::isfinite(eff) || eff <= 0. || (cfgMultiplicityPoisson && eff > 1.)) { ++skippedTracks; continue; } - correctedMultiplicity += 1. / eff; + if (cfgMultiplicityPoisson) { + correctedMultiplicity += 1.; + totalMeanMissing += 1. / eff - 1.; + } else { + correctedMultiplicity += 1. / eff; + } + } + if (cfgMultiplicityPoisson) { + if (!std::isfinite(totalMeanMissing) || totalMeanMissing >= static_cast(std::numeric_limits::max())) { + LOGF(fatal, "Poisson multiplicity mean exceeds the sampling range: %g", totalMeanMissing); + return 0.f; + } + // Sample from sum of missing count (Poisson sum rule) + if (totalMeanMissing > 0.) { + // Reseed per collision + const auto collisionIndex = static_cast(collision.globalIndex()); + std::seed_seq seed{cfgMultiplicityPoissonSeed.value, static_cast(runNumber), + static_cast(globalBC), static_cast(globalBC >> 32), + static_cast(collisionIndex), static_cast(collisionIndex >> 32)}; + // Speed up the seeding by creating one 64bit seed from the event information + std::array seedWords{}; + seed.generate(seedWords.begin(), seedWords.end()); + const uint64_t eventSeed = (static_cast(seedWords[0]) << std::numeric_limits::digits) | seedWords[1]; + mMultiplicityGenerator.seed(eventSeed); + correctedMultiplicity += std::poisson_distribution(totalMeanMissing)(mMultiplicityGenerator); + } } if (cfgVerbosity > 0 && skippedTracks > 0) { - LOGF(warning, "Skipped %zu tracks with invalid efficiency while correcting collision %lld", skippedTracks, static_cast(collision.globalIndex())); + LOGF(warning, "Skipped %zu tracks with invalid or unsupported efficiency while correcting collision %lld", skippedTracks, static_cast(collision.globalIndex())); } if (!std::isfinite(correctedMultiplicity) || correctedMultiplicity > std::numeric_limits::max()) { LOGF(fatal, "Corrected multiplicity cannot be represented as a float: %g", correctedMultiplicity); @@ -454,7 +484,7 @@ struct FilterCF { auto bc = collision.template bc_as(); outputCollisions(bc.runNumber(), collision.posZ(), collision.multiplicity(), bc.timestamp()); if (!cfgEfficiencyMultiplicity.value.empty()) { - outputCollisionsExtra(getCorrectedMultiplicity(collision, tracks, bc.timestamp())); + outputCollisionsExtra(getCorrectedMultiplicity(collision, tracks, bc.timestamp(), bc.globalBC(), bc.runNumber())); } if constexpr (std::experimental::is_detected::value) { @@ -652,7 +682,7 @@ struct FilterCF { // NOTE works only when we store all MC collisions (as we do here) outputCollisions(bc.runNumber(), collision.posZ(), collision.multiplicity(), bc.timestamp()); if (!cfgEfficiencyMultiplicity.value.empty()) { - outputCollisionsExtra(getCorrectedMultiplicity(collision, groupedTracks, bc.timestamp())); + outputCollisionsExtra(getCorrectedMultiplicity(collision, groupedTracks, bc.timestamp(), bc.globalBC(), bc.runNumber())); } const auto mcCollisionId = collision.mcCollisionId(); From 03bf8c5caeec0421e5b548b8064dba5176628817 Mon Sep 17 00:00:00 2001 From: Emil Gorm Nielsen Date: Mon, 9 Nov 2026 23:32:01 +0100 Subject: [PATCH 2/4] codecheck --- PWGCF/TableProducer/filterCorrelations.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/PWGCF/TableProducer/filterCorrelations.cxx b/PWGCF/TableProducer/filterCorrelations.cxx index 307a27d9e3d..48097573838 100644 --- a/PWGCF/TableProducer/filterCorrelations.cxx +++ b/PWGCF/TableProducer/filterCorrelations.cxx @@ -130,7 +130,7 @@ struct FilterCF { // TODO how to have this in the second task? For now they are copied Filter trackFilter = (nabs(aod::track::eta) < cfgCutEta) && (aod::track::pt > cfgCutPt); - Filter trackSelection = (requireGlobalTrackInFilter()) || (aod::track::isGlobalTrackSDD == (uint8_t)true); + Filter trackSelection = (requireGlobalTrackInFilter()) || (aod::track::isGlobalTrackSDD == static_cast(true)); Filter mcCollisionFilter = nabs(aod::mccollision::posZ) < cfgCutVertex; @@ -798,7 +798,7 @@ struct MultiplicitySelector { O2_DEFINE_CONFIGURABLE(cfgCutEta, float, 0.8f, "Eta range for tracks"); Filter trackFilter = (nabs(aod::track::eta) < cfgCutEta) && (aod::track::pt > cfgCutPt); - Filter trackSelection = (requireGlobalTrackInFilter()) || (aod::track::isGlobalTrackSDD == (uint8_t)true); + Filter trackSelection = (requireGlobalTrackInFilter()) || (aod::track::isGlobalTrackSDD == static_cast(true)); void init(InitContext&) { From b195d68b0f3fa60bd680ef599c55151b4e6346bc Mon Sep 17 00:00:00 2001 From: Emil Gorm Nielsen Date: Tue, 10 Nov 2026 20:37:44 +0100 Subject: [PATCH 3/4] Switch to seeding per dataframe --- PWGCF/TableProducer/filterCorrelations.cxx | 54 ++++++++++++++++------ 1 file changed, 39 insertions(+), 15 deletions(-) diff --git a/PWGCF/TableProducer/filterCorrelations.cxx b/PWGCF/TableProducer/filterCorrelations.cxx index 48097573838..144b3a169e8 100644 --- a/PWGCF/TableProducer/filterCorrelations.cxx +++ b/PWGCF/TableProducer/filterCorrelations.cxx @@ -121,8 +121,8 @@ struct FilterCF { O2_DEFINE_CONFIGURABLE(cfgEfficiencyMultiplicity, std::string, "", "Multiplicity efficiency (RecoAll / MC): CCDB path or local ROOT file with a 4D ccdb_object (eta, pT, multiplicity, z-vtx); empty disables CFCollisionsExtra output"); O2_DEFINE_CONFIGURABLE(cfgLocalEfficiency, int, 0, "0 = CCDB efficiency, 1 = local ROOT efficiency"); O2_DEFINE_CONFIGURABLE(cfgMultiplicityTrackBitMask, uint16_t, 0, "Required track-type bits for corrected multiplicity; match cfgTrackBitMask used to produce the efficiency (0 = all stored tracks)"); - O2_DEFINE_CONFIGURABLE(cfgMultiplicityPoisson, bool, false, "Use (reproducibly) seeded 1 + Poisson(sum(1/eff - 1)) for multiplicity correction; skips eff > 1 tracks"); - O2_DEFINE_CONFIGURABLE(cfgMultiplicityPoissonSeed, uint32_t, 1, "Base seed for reproducible per-collision Poisson multiplicity sampling"); + O2_DEFINE_CONFIGURABLE(cfgMultiplicityPoisson, bool, false, "Use reproducibly seeded Nvalid + Poisson(sum(1/eff - 1)) for multiplicity correction; skips eff > 1 tracks"); + O2_DEFINE_CONFIGURABLE(cfgMultiplicityPoissonSeed, uint32_t, 1, "Base seed for reproducible per-dataframe Poisson multiplicity sampling"); // Filters and input definitions Filter collisionZVtxFilter = nabs(aod::collision::posZ) < cfgCutVertex; @@ -172,6 +172,9 @@ struct FilterCF { void init(InitContext& initContext) { + if (cfgMultiplicityPoisson && !doprocessDataframeSeed) { + LOGF(fatal, "Poisson multiplicity sampling requires processDataframeSeed"); + } if (!cfgEfficiencyMultiplicity.value.empty()) { bool processTracksEnabled = false; if (!o2::common::core::getTaskOptionValue(initContext, "multiplicity-selector", "processTracks", processTracksEnabled, false)) { @@ -393,7 +396,7 @@ struct FilterCF { } template - float getCorrectedMultiplicity(const TCollision& collision, const TTracks& tracks, uint64_t timestamp, uint64_t globalBC, int runNumber) + float getCorrectedMultiplicity(const TCollision& collision, const TTracks& tracks, uint64_t timestamp) { auto* efficiency = loadMultiplicityEfficiency(timestamp); double correctedMultiplicity = 0.; @@ -426,16 +429,6 @@ struct FilterCF { } // Sample from sum of missing count (Poisson sum rule) if (totalMeanMissing > 0.) { - // Reseed per collision - const auto collisionIndex = static_cast(collision.globalIndex()); - std::seed_seq seed{cfgMultiplicityPoissonSeed.value, static_cast(runNumber), - static_cast(globalBC), static_cast(globalBC >> 32), - static_cast(collisionIndex), static_cast(collisionIndex >> 32)}; - // Speed up the seeding by creating one 64bit seed from the event information - std::array seedWords{}; - seed.generate(seedWords.begin(), seedWords.end()); - const uint64_t eventSeed = (static_cast(seedWords[0]) << std::numeric_limits::digits) | seedWords[1]; - mMultiplicityGenerator.seed(eventSeed); correctedMultiplicity += std::poisson_distribution(totalMeanMissing)(mMultiplicityGenerator); } } @@ -464,6 +457,37 @@ struct FilterCF { return true; } + // This process must run before the collision processes below. + void processDataframeSeed(aod::Origins const& origins, aod::BCs const& bcs) + { + if (!cfgMultiplicityPoisson) { + return; + } + if (origins.size() == 0) { + LOGF(fatal, "Cannot seed Poisson multiplicity sampling without a dataframe ID"); + return; + } + + const uint32_t runNumber = bcs.size() > 0 ? static_cast(bcs.iteratorAt(0).runNumber()) : 0u; + std::vector seedInputs; + seedInputs.reserve(2 + 2 * origins.size()); + seedInputs.push_back(cfgMultiplicityPoissonSeed.value); + seedInputs.push_back(runNumber); + for (const auto& origin : origins) { + const uint64_t dataframeId = origin.dataframeID(); + seedInputs.push_back(static_cast(dataframeId)); + seedInputs.push_back(static_cast(dataframeId >> std::numeric_limits::digits)); + } + + std::seed_seq seed(seedInputs.begin(), seedInputs.end()); + std::array seedWords{}; + seed.generate(seedWords.begin(), seedWords.end()); + // Speed up seeding by creating one 64-bit seed from the dataframe information. + const uint64_t dataframeSeed = (static_cast(seedWords[0]) << std::numeric_limits::digits) | seedWords[1]; + mMultiplicityGenerator.seed(dataframeSeed); + } + PROCESS_SWITCH(FilterCF, processDataframeSeed, "Seed Poisson multiplicity sampling for each dataframe", false); + template using HasMultTables = decltype(std::declval().multNTracksPV()); @@ -484,7 +508,7 @@ struct FilterCF { auto bc = collision.template bc_as(); outputCollisions(bc.runNumber(), collision.posZ(), collision.multiplicity(), bc.timestamp()); if (!cfgEfficiencyMultiplicity.value.empty()) { - outputCollisionsExtra(getCorrectedMultiplicity(collision, tracks, bc.timestamp(), bc.globalBC(), bc.runNumber())); + outputCollisionsExtra(getCorrectedMultiplicity(collision, tracks, bc.timestamp())); } if constexpr (std::experimental::is_detected::value) { @@ -682,7 +706,7 @@ struct FilterCF { // NOTE works only when we store all MC collisions (as we do here) outputCollisions(bc.runNumber(), collision.posZ(), collision.multiplicity(), bc.timestamp()); if (!cfgEfficiencyMultiplicity.value.empty()) { - outputCollisionsExtra(getCorrectedMultiplicity(collision, groupedTracks, bc.timestamp(), bc.globalBC(), bc.runNumber())); + outputCollisionsExtra(getCorrectedMultiplicity(collision, groupedTracks, bc.timestamp())); } const auto mcCollisionId = collision.mcCollisionId(); From 21e51391ea83d5f9c60c0f67210e59316ac309ac Mon Sep 17 00:00:00 2001 From: Emil Gorm Nielsen Date: Tue, 10 Nov 2026 22:00:47 +0100 Subject: [PATCH 4/4] use only dataframe id for seeding --- PWGCF/TableProducer/filterCorrelations.cxx | 18 ++++-------------- 1 file changed, 4 insertions(+), 14 deletions(-) diff --git a/PWGCF/TableProducer/filterCorrelations.cxx b/PWGCF/TableProducer/filterCorrelations.cxx index 144b3a169e8..3271624b1b9 100644 --- a/PWGCF/TableProducer/filterCorrelations.cxx +++ b/PWGCF/TableProducer/filterCorrelations.cxx @@ -458,7 +458,7 @@ struct FilterCF { } // This process must run before the collision processes below. - void processDataframeSeed(aod::Origins const& origins, aod::BCs const& bcs) + void processDataframeSeed(aod::Origins const& origins) { if (!cfgMultiplicityPoisson) { return; @@ -467,19 +467,9 @@ struct FilterCF { LOGF(fatal, "Cannot seed Poisson multiplicity sampling without a dataframe ID"); return; } - - const uint32_t runNumber = bcs.size() > 0 ? static_cast(bcs.iteratorAt(0).runNumber()) : 0u; - std::vector seedInputs; - seedInputs.reserve(2 + 2 * origins.size()); - seedInputs.push_back(cfgMultiplicityPoissonSeed.value); - seedInputs.push_back(runNumber); - for (const auto& origin : origins) { - const uint64_t dataframeId = origin.dataframeID(); - seedInputs.push_back(static_cast(dataframeId)); - seedInputs.push_back(static_cast(dataframeId >> std::numeric_limits::digits)); - } - - std::seed_seq seed(seedInputs.begin(), seedInputs.end()); + const uint64_t dataframeId = origins.iteratorAt(0).dataframeID(); + std::seed_seq seed{cfgMultiplicityPoissonSeed.value, static_cast(dataframeId), + static_cast(dataframeId >> std::numeric_limits::digits)}; std::array seedWords{}; seed.generate(seedWords.begin(), seedWords.end()); // Speed up seeding by creating one 64-bit seed from the dataframe information.