diff --git a/Generators/include/Generators/TPCLoopers.h b/Generators/include/Generators/TPCLoopers.h index 3e8685257b829..d979f0ff569e6 100644 --- a/Generators/include/Generators/TPCLoopers.h +++ b/Generators/include/Generators/TPCLoopers.h @@ -112,6 +112,15 @@ class GenTPCLoopers // check if a vertex lies in the TPC volume where ionisation can be recorded bool isInTPCActiveVolume(double vx, double vy) const; + // false if every collision of this simulation is owned by an earlier source of the collision context + static bool ownsContextCollisions(); + + // true if the loopers of a collision are to be generated by the given source of the collision context + static bool ownsCollision(const std::vector& parts, int source); + + // work out the share of the timeframe each event of this simulation has to cover with loopers + void setupContextCoverage(); + unsigned int getNLoopers() const { return (mNLoopersPairs + mNLoopersCompton); } // loopers dropped by the geometrical protection since the last reset @@ -143,6 +152,8 @@ class GenTPCLoopers double mIntTimeRecMean = 1.0; // Average interaction time record used for the reference double mTimeLimit = 0.0; // Time limit for the current event double mTimeEnd = 0.0; // Time limit for the last event + int mSourceID = -1; // ID of this simulation in the collision context (-1 if unknown) + std::vector mEventLooperTime; // per event, the time of the timeframe it has to cover with loopers (empty if there is no context) float mLoopsFractionPairs = 0.08; // Fraction of loopers from Pairs int mInteractionRate = 50000; // Interaction rate in Hz bool mGeomProtection = true; // Skip loopers generated outside the TPC active volume diff --git a/Generators/src/Generator.cxx b/Generators/src/Generator.cxx index 63b0fa36d5016..02f3c1a1e2970 100644 --- a/Generators/src/Generator.cxx +++ b/Generators/src/Generator.cxx @@ -164,6 +164,14 @@ bool Generator::initTPCLoopersGen() return kFALSE; } } + // Loopers are a detector background continuous in time, and a single simulation covers the whole + // timeframe: when several of them feed the same collision context (embedding), only the source + // owning a collision generates its loopers, otherwise each source will generate them for the same time intervals. + if (!o2::eventgen::GenTPCLoopers::ownsContextCollisions()) { + LOG(info) << "All collisions of this simulation are also taken by a source which comes earlier in the collision context and generates their loopers."; + LOG(info) << "Fast loopers generator will remain OFF."; + return kFALSE; + } std::array multiplier = {loopersParam.multiplier[0], loopersParam.multiplier[1]}; unsigned int nLoopersPairs = loopersParam.fixedNLoopers[0]; unsigned int nLoopersCompton = loopersParam.fixedNLoopers[1]; diff --git a/Generators/src/TPCLoopers.cxx b/Generators/src/TPCLoopers.cxx index 8f1125bf94467..f7e2e8d53c6df 100644 --- a/Generators/src/TPCLoopers.cxx +++ b/Generators/src/TPCLoopers.cxx @@ -15,8 +15,12 @@ #include "CCDB/CCDBTimeStampUtils.h" #include "CCDB/CcdbApi.h" #include "DetectorsRaw/HBFUtils.h" +#include #include "TF1.h" #include +#include +#include +#include #include #include "SimulationDataFormat/MCGenProperties.h" #include @@ -160,6 +164,105 @@ void GenTPCLoopers::setGeomProtection(bool protect) } } +bool GenTPCLoopers::ownsCollision(const std::vector& parts, int source) +{ + // A collision is owned by the source with the lowest ID among the event parts it holds, so that + // exactly one of the simulations feeding a collision context (i.e. bkg, sgn) generates its loopers. + int lowest = std::numeric_limits::max(); + for (const auto& part : parts) { + lowest = std::min(lowest, part.sourceID); + } + return lowest == source; +} + +bool GenTPCLoopers::ownsContextCollisions() +{ + // Loopers are a detector background which is continuous in time: event by event this generator + // covers the time interval between consecutive collisions of the collision context, so a single + // simulation tiles the whole timeframe exactly once. When several transport simulations feed the + // same context (embedding), each of them would tile that same time axis again and the loopers + // would end up in the digits as many times as there are simulations. + // Every collision is therefore owned by exactly one source, the one with the lowest source ID + // among the event parts present, and only the owner generates loopers for it. Embedding without + // replacement ("@" patterns) leaves all collisions to the background, since its event part is + // present everywhere; when the signal replaces the background ("r" patterns) the background part + // is taken out of the collisions it loses, so each source covers exactly its own share. + const auto& simConfig = o2::conf::SimConfig::Instance(); + const auto contextAndPrefix = simConfig.getCollContextFilenameAndEventPrefix(); + const std::string contextName = contextAndPrefix.first.empty() ? "collisioncontext.root" : contextAndPrefix.first; + const std::string& simPrefix = contextAndPrefix.second; + if (!std::filesystem::exists(contextName)) { + return true; // no context to share with anybody + } + std::unique_ptr context(o2::steer::DigitizationContext::loadFromFile(contextName)); + if (!context) { + LOG(warn) << "Could not read the collision context from " << contextName + << "; loopers will be generated without checking for other sources"; + return true; + } + const auto source = context->findSimPrefix(simPrefix); + if (source < 0) { + LOG(warn) << "Simulation prefix " << simPrefix << " not found in the collision context " << contextName + << "; loopers will be generated without checking for other sources"; + return true; + } + // count the collisions this source takes part in and those it owns + unsigned int nPresent = 0, nOwned = 0; + for (const auto& parts : context->getEventParts()) { + if (ownsCollision(parts, source)) { + nOwned++; + } + if (std::find_if(parts.begin(), parts.end(), [source](const auto& part) { return part.sourceID == source; }) != parts.end()) { + nPresent++; + } + } + LOG(info) << "Loopers: source " << simPrefix << " owns " << nOwned << " of the " << nPresent + << " collisions it takes part in, out of " << context->getEventParts().size() << " in " << contextName; + return nOwned > 0; +} + +void GenTPCLoopers::setupContextCoverage() +{ + // Work out how much of the time axis every event of this simulation has to cover with loopers. + // A collision stands for the time from itself to the next one (the last one to the end of the + // frame), and an event covers the time of the collisions it is used at: when it is used at + // several of them, as the background of an embedding production is, the loopers it carries are + // digitised once per use, so it covers their average. Collisions owned by another source count + // as zero. The loopers reaching the digits then add up to the reference number per orbit times + // the time this source owns, however its events are spread over the collision context. + mEventLooperTime.clear(); + mSourceID = mCollisionContext ? mCollisionContext->findSimPrefix(o2::conf::SimConfig::Instance().getCollContextFilenameAndEventPrefix().second) : -1; + if (mSourceID < 0) { + return; // no context information: every event keeps the interval of the collision of same index + } + const auto& parts = mCollisionContext->getEventParts(); + std::vector nUses; + double ownedTime = 0.; + for (size_t c = 0; c < parts.size() && c < mInteractionTimeRecords.size(); ++c) { + const double gap = (c + 1 < mInteractionTimeRecords.size()) + ? mInteractionTimeRecords[c + 1].bc2ns() - mInteractionTimeRecords[c].bc2ns() + : mTimeEnd - mInteractionTimeRecords[c].bc2ns(); + const double owned = ownsCollision(parts[c], mSourceID) ? gap : 0.; + ownedTime += owned; + for (const auto& part : parts[c]) { + if (part.sourceID != mSourceID) { + continue; + } + if (part.entryID >= (int)mEventLooperTime.size()) { // entryIDs are dense, see SimConfig + mEventLooperTime.resize(part.entryID + 1, 0.); + nUses.resize(part.entryID + 1, 0); + } + mEventLooperTime[part.entryID] += owned; + nUses[part.entryID]++; + } + } + for (size_t e = 0; e < mEventLooperTime.size(); ++e) { + mEventLooperTime[e] /= std::max(nUses[e], 1); + } + LOG(info) << "Loopers: the " << mEventLooperTime.size() << " events of this simulation share the " + << ownedTime << " ns of the timeframe it owns"; +} + GenTPCLoopers::GenTPCLoopers(std::string model_pairs, std::string model_compton, std::string poisson, std::string gauss, std::string scaler_pair, std::string scaler_compton) @@ -233,9 +336,12 @@ Bool_t GenTPCLoopers::generateEvent() if (mFlatGas) { unsigned int nLoopers, nLoopersPairs, nLoopersCompton; LOG(debug) << "mCurrentEvent is " << mCurrentEvent; - 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"; - LOG(debug) << "Current time offset wrt BC: " << mInteractionTimeRecords[mCurrentEvent].getTimeOffsetWrtBC() << " ns"; - mTimeLimit = (mCurrentEvent < mInteractionTimeRecords.size() - 1) ? mInteractionTimeRecords[mCurrentEvent + 1].bc2ns() - mInteractionTimeRecords[mCurrentEvent].bc2ns() : mTimeEnd - mInteractionTimeRecords[mCurrentEvent].bc2ns(); + // the share of the timeframe this event has to cover, or, with no collision context to look + // at, the interval between the collision of same index and the next one + mTimeLimit = (mCurrentEvent < (int)mEventLooperTime.size()) + ? mEventLooperTime[mCurrentEvent] + : ((mCurrentEvent < mInteractionTimeRecords.size() - 1) ? mInteractionTimeRecords[mCurrentEvent + 1].bc2ns() - mInteractionTimeRecords[mCurrentEvent].bc2ns() : mTimeEnd - mInteractionTimeRecords[mCurrentEvent].bc2ns()); + LOG(debug) << "Current event time: " << mTimeLimit << " ns"; // With flat gas the number of loopers are adapted based on time interval widths // The denominator is either the LHC orbit (if mFlatGasOrbit is true) or the mean interaction time record interval 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) o2::InteractionRecord finalOrbitIR(0, lastIR.orbit + 2); // Final orbit, BC = 0 mTimeEnd = finalOrbitIR.bc2ns(); LOG(debug) << "Final orbit start time: " << mTimeEnd << " ns while last interaction record time is " << mInteractionTimeRecords.back().bc2ns() << " ns"; + setupContextCoverage(); } } } else {