Skip to content

Commit 0c4311d

Browse files
committed
Fix loopers double counting when using embedding
1 parent f98960a commit 0c4311d

3 files changed

Lines changed: 129 additions & 3 deletions

File tree

‎Generators/include/Generators/TPCLoopers.h‎

Lines changed: 11 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -112,6 +112,15 @@ class GenTPCLoopers
112112
// check if a vertex lies in the TPC volume where ionisation can be recorded
113113
bool isInTPCActiveVolume(double vx, double vy) const;
114114

115+
// false if every collision of this simulation is owned by an earlier source of the collision context
116+
static bool ownsContextCollisions();
117+
118+
// true if the loopers of a collision are to be generated by the given source of the collision context
119+
static bool ownsCollision(const std::vector<o2::steer::EventPart>& parts, int source);
120+
121+
// work out the share of the timeframe each event of this simulation has to cover with loopers
122+
void setupContextCoverage();
123+
115124
unsigned int getNLoopers() const { return (mNLoopersPairs + mNLoopersCompton); }
116125

117126
// loopers dropped by the geometrical protection since the last reset
@@ -143,6 +152,8 @@ class GenTPCLoopers
143152
double mIntTimeRecMean = 1.0; // Average interaction time record used for the reference
144153
double mTimeLimit = 0.0; // Time limit for the current event
145154
double mTimeEnd = 0.0; // Time limit for the last event
155+
int mSourceID = -1; // ID of this simulation in the collision context (-1 if unknown)
156+
std::vector<double> mEventLooperTime; // per event, the time of the timeframe it has to cover with loopers (empty if there is no context)
146157
float mLoopsFractionPairs = 0.08; // Fraction of loopers from Pairs
147158
int mInteractionRate = 50000; // Interaction rate in Hz
148159
bool mGeomProtection = true; // Skip loopers generated outside the TPC active volume

‎Generators/src/Generator.cxx‎

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -164,6 +164,14 @@ bool Generator::initTPCLoopersGen()
164164
return kFALSE;
165165
}
166166
}
167+
// Loopers are a detector background continuous in time, and a single simulation covers the whole
168+
// timeframe: when several of them feed the same collision context (embedding), only the source
169+
// owning a collision generates its loopers, otherwise each source will generate them for the same time intervals.
170+
if (!o2::eventgen::GenTPCLoopers::ownsContextCollisions()) {
171+
LOG(info) << "All collisions of this simulation are also taken by a source which comes earlier in the collision context and generates their loopers.";
172+
LOG(info) << "Fast loopers generator will remain OFF.";
173+
return kFALSE;
174+
}
167175
std::array<float, 2> multiplier = {loopersParam.multiplier[0], loopersParam.multiplier[1]};
168176
unsigned int nLoopersPairs = loopersParam.fixedNLoopers[0];
169177
unsigned int nLoopersCompton = loopersParam.fixedNLoopers[1];

‎Generators/src/TPCLoopers.cxx‎

Lines changed: 110 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -15,8 +15,12 @@
1515
#include "CCDB/CCDBTimeStampUtils.h"
1616
#include "CCDB/CcdbApi.h"
1717
#include "DetectorsRaw/HBFUtils.h"
18+
#include <SimConfig/SimConfig.h>
1819
#include "TF1.h"
1920
#include <filesystem>
21+
#include <algorithm>
22+
#include <limits>
23+
#include <memory>
2024
#include <SimulationDataFormat/ParticleStatus.h>
2125
#include "SimulationDataFormat/MCGenProperties.h"
2226
#include <iostream>
@@ -160,6 +164,105 @@ void GenTPCLoopers::setGeomProtection(bool protect)
160164
}
161165
}
162166

167+
bool GenTPCLoopers::ownsCollision(const std::vector<o2::steer::EventPart>& parts, int source)
168+
{
169+
// A collision is owned by the source with the lowest ID among the event parts it holds, so that
170+
// exactly one of the simulations feeding a collision context (i.e. bkg, sgn) generates its loopers.
171+
int lowest = std::numeric_limits<int>::max();
172+
for (const auto& part : parts) {
173+
lowest = std::min(lowest, part.sourceID);
174+
}
175+
return lowest == source;
176+
}
177+
178+
bool GenTPCLoopers::ownsContextCollisions()
179+
{
180+
// Loopers are a detector background which is continuous in time: event by event this generator
181+
// covers the time interval between consecutive collisions of the collision context, so a single
182+
// simulation tiles the whole timeframe exactly once. When several transport simulations feed the
183+
// same context (embedding), each of them would tile that same time axis again and the loopers
184+
// would end up in the digits as many times as there are simulations.
185+
// Every collision is therefore owned by exactly one source, the one with the lowest source ID
186+
// among the event parts present, and only the owner generates loopers for it. Embedding without
187+
// replacement ("@" patterns) leaves all collisions to the background, since its event part is
188+
// present everywhere; when the signal replaces the background ("r" patterns) the background part
189+
// is taken out of the collisions it loses, so each source covers exactly its own share.
190+
const auto& simConfig = o2::conf::SimConfig::Instance();
191+
const auto contextAndPrefix = simConfig.getCollContextFilenameAndEventPrefix();
192+
const std::string contextName = contextAndPrefix.first.empty() ? "collisioncontext.root" : contextAndPrefix.first;
193+
const std::string& simPrefix = contextAndPrefix.second;
194+
if (!std::filesystem::exists(contextName)) {
195+
return true; // no context to share with anybody
196+
}
197+
std::unique_ptr<o2::steer::DigitizationContext> context(o2::steer::DigitizationContext::loadFromFile(contextName));
198+
if (!context) {
199+
LOG(warn) << "Could not read the collision context from " << contextName
200+
<< "; loopers will be generated without checking for other sources";
201+
return true;
202+
}
203+
const auto source = context->findSimPrefix(simPrefix);
204+
if (source < 0) {
205+
LOG(warn) << "Simulation prefix " << simPrefix << " not found in the collision context " << contextName
206+
<< "; loopers will be generated without checking for other sources";
207+
return true;
208+
}
209+
// count the collisions this source takes part in and those it owns
210+
unsigned int nPresent = 0, nOwned = 0;
211+
for (const auto& parts : context->getEventParts()) {
212+
if (ownsCollision(parts, source)) {
213+
nOwned++;
214+
}
215+
if (std::find_if(parts.begin(), parts.end(), [source](const auto& part) { return part.sourceID == source; }) != parts.end()) {
216+
nPresent++;
217+
}
218+
}
219+
LOG(info) << "Loopers: source " << simPrefix << " owns " << nOwned << " of the " << nPresent
220+
<< " collisions it takes part in, out of " << context->getEventParts().size() << " in " << contextName;
221+
return nOwned > 0;
222+
}
223+
224+
void GenTPCLoopers::setupContextCoverage()
225+
{
226+
// Work out how much of the time axis every event of this simulation has to cover with loopers.
227+
// A collision stands for the time from itself to the next one (the last one to the end of the
228+
// frame), and an event covers the time of the collisions it is used at: when it is used at
229+
// several of them, as the background of an embedding production is, the loopers it carries are
230+
// digitised once per use, so it covers their average. Collisions owned by another source count
231+
// as zero. The loopers reaching the digits then add up to the reference number per orbit times
232+
// the time this source owns, however its events are spread over the collision context.
233+
mEventLooperTime.clear();
234+
mSourceID = mCollisionContext ? mCollisionContext->findSimPrefix(o2::conf::SimConfig::Instance().getCollContextFilenameAndEventPrefix().second) : -1;
235+
if (mSourceID < 0) {
236+
return; // no context information: every event keeps the interval of the collision of same index
237+
}
238+
const auto& parts = mCollisionContext->getEventParts();
239+
std::vector<int> nUses;
240+
double ownedTime = 0.;
241+
for (size_t c = 0; c < parts.size() && c < mInteractionTimeRecords.size(); ++c) {
242+
const double gap = (c + 1 < mInteractionTimeRecords.size())
243+
? mInteractionTimeRecords[c + 1].bc2ns() - mInteractionTimeRecords[c].bc2ns()
244+
: mTimeEnd - mInteractionTimeRecords[c].bc2ns();
245+
const double owned = ownsCollision(parts[c], mSourceID) ? gap : 0.;
246+
ownedTime += owned;
247+
for (const auto& part : parts[c]) {
248+
if (part.sourceID != mSourceID) {
249+
continue;
250+
}
251+
if (part.entryID >= (int)mEventLooperTime.size()) { // entryIDs are dense, see SimConfig
252+
mEventLooperTime.resize(part.entryID + 1, 0.);
253+
nUses.resize(part.entryID + 1, 0);
254+
}
255+
mEventLooperTime[part.entryID] += owned;
256+
nUses[part.entryID]++;
257+
}
258+
}
259+
for (size_t e = 0; e < mEventLooperTime.size(); ++e) {
260+
mEventLooperTime[e] /= std::max(nUses[e], 1);
261+
}
262+
LOG(info) << "Loopers: the " << mEventLooperTime.size() << " events of this simulation share the "
263+
<< ownedTime << " ns of the timeframe it owns";
264+
}
265+
163266
GenTPCLoopers::GenTPCLoopers(std::string model_pairs, std::string model_compton,
164267
std::string poisson, std::string gauss, std::string scaler_pair,
165268
std::string scaler_compton)
@@ -233,9 +336,12 @@ Bool_t GenTPCLoopers::generateEvent()
233336
if (mFlatGas) {
234337
unsigned int nLoopers, nLoopersPairs, nLoopersCompton;
235338
LOG(debug) << "mCurrentEvent is " << mCurrentEvent;
236-
LOG(debug) << "Current event time: " << ((mCurrentEvent < mInteractionTimeRecords.size() - 1) ? std::to_string(mInteractionTimeRecords[mCurrentEvent + 1].bc2ns() - mInteractionTimeRecords[mCurrentEvent].bc2ns()) : std::to_string(mTimeEnd - mInteractionTimeRecords[mCurrentEvent].bc2ns())) << " ns";
237-
LOG(debug) << "Current time offset wrt BC: " << mInteractionTimeRecords[mCurrentEvent].getTimeOffsetWrtBC() << " ns";
238-
mTimeLimit = (mCurrentEvent < mInteractionTimeRecords.size() - 1) ? mInteractionTimeRecords[mCurrentEvent + 1].bc2ns() - mInteractionTimeRecords[mCurrentEvent].bc2ns() : mTimeEnd - mInteractionTimeRecords[mCurrentEvent].bc2ns();
339+
// the share of the timeframe this event has to cover, or, with no collision context to look
340+
// at, the interval between the collision of same index and the next one
341+
mTimeLimit = (mCurrentEvent < (int)mEventLooperTime.size())
342+
? mEventLooperTime[mCurrentEvent]
343+
: ((mCurrentEvent < mInteractionTimeRecords.size() - 1) ? mInteractionTimeRecords[mCurrentEvent + 1].bc2ns() - mInteractionTimeRecords[mCurrentEvent].bc2ns() : mTimeEnd - mInteractionTimeRecords[mCurrentEvent].bc2ns());
344+
LOG(debug) << "Current event time: " << mTimeLimit << " ns";
239345
// With flat gas the number of loopers are adapted based on time interval widths
240346
// The denominator is either the LHC orbit (if mFlatGasOrbit is true) or the mean interaction time record interval
241347
nLoopers = mFlatGasOrbit ? (mFlatGasNumber * (mTimeLimit / o2::constants::lhc::LHCOrbitNS)) : (mFlatGasNumber * (mTimeLimit / mIntTimeRecMean));
@@ -472,6 +578,7 @@ void GenTPCLoopers::setFlatGas(Bool_t flat, Int_t number, Int_t nloopers_orbit)
472578
o2::InteractionRecord finalOrbitIR(0, lastIR.orbit + 2); // Final orbit, BC = 0
473579
mTimeEnd = finalOrbitIR.bc2ns();
474580
LOG(debug) << "Final orbit start time: " << mTimeEnd << " ns while last interaction record time is " << mInteractionTimeRecords.back().bc2ns() << " ns";
581+
setupContextCoverage();
475582
}
476583
}
477584
} else {

0 commit comments

Comments
 (0)