Skip to content

Commit b195d68

Browse files
committed
Switch to seeding per dataframe
1 parent 03bf8c5 commit b195d68

1 file changed

Lines changed: 39 additions & 15 deletions

File tree

‎PWGCF/TableProducer/filterCorrelations.cxx‎

Lines changed: 39 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -121,8 +121,8 @@ struct FilterCF {
121121
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");
122122
O2_DEFINE_CONFIGURABLE(cfgLocalEfficiency, int, 0, "0 = CCDB efficiency, 1 = local ROOT efficiency");
123123
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)");
124-
O2_DEFINE_CONFIGURABLE(cfgMultiplicityPoisson, bool, false, "Use (reproducibly) seeded 1 + Poisson(sum(1/eff - 1)) for multiplicity correction; skips eff > 1 tracks");
125-
O2_DEFINE_CONFIGURABLE(cfgMultiplicityPoissonSeed, uint32_t, 1, "Base seed for reproducible per-collision Poisson multiplicity sampling");
124+
O2_DEFINE_CONFIGURABLE(cfgMultiplicityPoisson, bool, false, "Use reproducibly seeded Nvalid + Poisson(sum(1/eff - 1)) for multiplicity correction; skips eff > 1 tracks");
125+
O2_DEFINE_CONFIGURABLE(cfgMultiplicityPoissonSeed, uint32_t, 1, "Base seed for reproducible per-dataframe Poisson multiplicity sampling");
126126

127127
// Filters and input definitions
128128
Filter collisionZVtxFilter = nabs(aod::collision::posZ) < cfgCutVertex;
@@ -172,6 +172,9 @@ struct FilterCF {
172172

173173
void init(InitContext& initContext)
174174
{
175+
if (cfgMultiplicityPoisson && !doprocessDataframeSeed) {
176+
LOGF(fatal, "Poisson multiplicity sampling requires processDataframeSeed");
177+
}
175178
if (!cfgEfficiencyMultiplicity.value.empty()) {
176179
bool processTracksEnabled = false;
177180
if (!o2::common::core::getTaskOptionValue(initContext, "multiplicity-selector", "processTracks", processTracksEnabled, false)) {
@@ -393,7 +396,7 @@ struct FilterCF {
393396
}
394397

395398
template <typename TCollision, typename TTracks>
396-
float getCorrectedMultiplicity(const TCollision& collision, const TTracks& tracks, uint64_t timestamp, uint64_t globalBC, int runNumber)
399+
float getCorrectedMultiplicity(const TCollision& collision, const TTracks& tracks, uint64_t timestamp)
397400
{
398401
auto* efficiency = loadMultiplicityEfficiency(timestamp);
399402
double correctedMultiplicity = 0.;
@@ -426,16 +429,6 @@ struct FilterCF {
426429
}
427430
// Sample from sum of missing count (Poisson sum rule)
428431
if (totalMeanMissing > 0.) {
429-
// Reseed per collision
430-
const auto collisionIndex = static_cast<uint64_t>(collision.globalIndex());
431-
std::seed_seq seed{cfgMultiplicityPoissonSeed.value, static_cast<uint32_t>(runNumber),
432-
static_cast<uint32_t>(globalBC), static_cast<uint32_t>(globalBC >> 32),
433-
static_cast<uint32_t>(collisionIndex), static_cast<uint32_t>(collisionIndex >> 32)};
434-
// Speed up the seeding by creating one 64bit seed from the event information
435-
std::array<uint32_t, 2> seedWords{};
436-
seed.generate(seedWords.begin(), seedWords.end());
437-
const uint64_t eventSeed = (static_cast<uint64_t>(seedWords[0]) << std::numeric_limits<uint32_t>::digits) | seedWords[1];
438-
mMultiplicityGenerator.seed(eventSeed);
439432
correctedMultiplicity += std::poisson_distribution<uint64_t>(totalMeanMissing)(mMultiplicityGenerator);
440433
}
441434
}
@@ -464,6 +457,37 @@ struct FilterCF {
464457
return true;
465458
}
466459

460+
// This process must run before the collision processes below.
461+
void processDataframeSeed(aod::Origins const& origins, aod::BCs const& bcs)
462+
{
463+
if (!cfgMultiplicityPoisson) {
464+
return;
465+
}
466+
if (origins.size() == 0) {
467+
LOGF(fatal, "Cannot seed Poisson multiplicity sampling without a dataframe ID");
468+
return;
469+
}
470+
471+
const uint32_t runNumber = bcs.size() > 0 ? static_cast<uint32_t>(bcs.iteratorAt(0).runNumber()) : 0u;
472+
std::vector<uint32_t> seedInputs;
473+
seedInputs.reserve(2 + 2 * origins.size());
474+
seedInputs.push_back(cfgMultiplicityPoissonSeed.value);
475+
seedInputs.push_back(runNumber);
476+
for (const auto& origin : origins) {
477+
const uint64_t dataframeId = origin.dataframeID();
478+
seedInputs.push_back(static_cast<uint32_t>(dataframeId));
479+
seedInputs.push_back(static_cast<uint32_t>(dataframeId >> std::numeric_limits<uint32_t>::digits));
480+
}
481+
482+
std::seed_seq seed(seedInputs.begin(), seedInputs.end());
483+
std::array<uint32_t, 2> seedWords{};
484+
seed.generate(seedWords.begin(), seedWords.end());
485+
// Speed up seeding by creating one 64-bit seed from the dataframe information.
486+
const uint64_t dataframeSeed = (static_cast<uint64_t>(seedWords[0]) << std::numeric_limits<uint32_t>::digits) | seedWords[1];
487+
mMultiplicityGenerator.seed(dataframeSeed);
488+
}
489+
PROCESS_SWITCH(FilterCF, processDataframeSeed, "Seed Poisson multiplicity sampling for each dataframe", false);
490+
467491
template <class T>
468492
using HasMultTables = decltype(std::declval<T&>().multNTracksPV());
469493

@@ -484,7 +508,7 @@ struct FilterCF {
484508
auto bc = collision.template bc_as<aod::BCsWithTimestamps>();
485509
outputCollisions(bc.runNumber(), collision.posZ(), collision.multiplicity(), bc.timestamp());
486510
if (!cfgEfficiencyMultiplicity.value.empty()) {
487-
outputCollisionsExtra(getCorrectedMultiplicity(collision, tracks, bc.timestamp(), bc.globalBC(), bc.runNumber()));
511+
outputCollisionsExtra(getCorrectedMultiplicity(collision, tracks, bc.timestamp()));
488512
}
489513

490514
if constexpr (std::experimental::is_detected<HasMultTables, C1>::value) {
@@ -682,7 +706,7 @@ struct FilterCF {
682706
// NOTE works only when we store all MC collisions (as we do here)
683707
outputCollisions(bc.runNumber(), collision.posZ(), collision.multiplicity(), bc.timestamp());
684708
if (!cfgEfficiencyMultiplicity.value.empty()) {
685-
outputCollisionsExtra(getCorrectedMultiplicity(collision, groupedTracks, bc.timestamp(), bc.globalBC(), bc.runNumber()));
709+
outputCollisionsExtra(getCorrectedMultiplicity(collision, groupedTracks, bc.timestamp()));
686710
}
687711

688712
const auto mcCollisionId = collision.mcCollisionId();

0 commit comments

Comments
 (0)