Skip to content

[PWGCF] Added option of poissonian sampling of efficiency to corrected multiplicity estimate in filterCorrelations.cxx #18672

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
EmilGorm wants to merge 1 commit into
base: master
Choose a base branch
from
+36 −6
Open

[PWGCF] Added option of poissonian sampling of efficiency to corrected multiplicity estimate in filterCorrelations.cxx #18672

Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
42 changes: 36 additions & 6 deletions PWGCF/TableProducer/filterCorrelations.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -55,6 +55,7 @@
#include <experimental/type_traits> // required for is_detected
#include <limits>
#include <memory>
#include <random>
#include <string>
#include <type_traits>
#include <vector>
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -160,6 +163,7 @@ struct FilterCF {
// Own local histograms independently of their input file. CCDB owns its objects.
std::unique_ptr<THn> localMultiplicityEfficiency;
THn* mEfficiency = nullptr;
std::mt19937_64 mMultiplicityGenerator;
static constexpr int MultiplicityEfficiencyDimensions = 4;

// persistent caches
Expand Down Expand Up @@ -389,11 +393,12 @@ struct FilterCF {
}

template <typename TCollision, typename TTracks>
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)) {
Expand All @@ -403,14 +408,39 @@ struct FilterCF {
// Keep the original estimator as the map coordinate, including for centrality.
const std::array<double, MultiplicityEfficiencyDimensions> 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<double>(std::numeric_limits<uint64_t>::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<uint64_t>(collision.globalIndex());
std::seed_seq seed{cfgMultiplicityPoissonSeed.value, static_cast<uint32_t>(runNumber),
static_cast<uint32_t>(globalBC), static_cast<uint32_t>(globalBC >> 32),
static_cast<uint32_t>(collisionIndex), static_cast<uint32_t>(collisionIndex >> 32)};
// Speed up the seeding by creating one 64bit seed from the event information
std::array<uint32_t, 2> seedWords{};
seed.generate(seedWords.begin(), seedWords.end());
const uint64_t eventSeed = (static_cast<uint64_t>(seedWords[0]) << std::numeric_limits<uint32_t>::digits) | seedWords[1];
mMultiplicityGenerator.seed(eventSeed);
correctedMultiplicity += std::poisson_distribution<uint64_t>(totalMeanMissing)(mMultiplicityGenerator);
}
}
if (cfgVerbosity > 0 && skippedTracks > 0) {
LOGF(warning, "Skipped %zu tracks with invalid efficiency while correcting collision %lld", skippedTracks, static_cast<int64_t>(collision.globalIndex()));
LOGF(warning, "Skipped %zu tracks with invalid or unsupported efficiency while correcting collision %lld", skippedTracks, static_cast<int64_t>(collision.globalIndex()));
}
if (!std::isfinite(correctedMultiplicity) || correctedMultiplicity > std::numeric_limits<float>::max()) {
LOGF(fatal, "Corrected multiplicity cannot be represented as a float: %g", correctedMultiplicity);
Expand Down Expand Up @@ -454,7 +484,7 @@ struct FilterCF {
auto bc = collision.template bc_as<aod::BCsWithTimestamps>();
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<HasMultTables, C1>::value) {
Expand Down Expand Up @@ -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();
Expand Down
Loading