diff --git a/CMakeLists.txt b/CMakeLists.txt index 3566b33b..b600f5dc 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -28,6 +28,7 @@ set(${PROJECT_NAME}_VERSION "${${PROJECT_NAME}_VERSION_MAJOR}.${${PROJECT_NAME}_ find_package(ROOT COMPONENTS RIO Tree REQUIRED) find_package(Gaudi REQUIRED) +find_package(TBB REQUIRED) find_package(podio 1.3 REQUIRED) find_package(EDM4HEP 1.0) if (NOT EDM4HEP_FOUND) diff --git a/doc/OverlayTiming.md b/doc/OverlayTiming.md index 320edccb..0fa29f7c 100644 --- a/doc/OverlayTiming.md +++ b/doc/OverlayTiming.md @@ -33,7 +33,10 @@ It uses [`UniqueIDGenSvc`](uniqueIDGen.md) to seed the internal random number ge | Property | Default | Description | |----------|---------|-------------| -| `BackgroundFileNames` | `[]` | List of groups of background input files, one group per overlay stream | +| `BackgroundFileNames` | `[]` | List of groups of background input files, one group per overlay stream. Entries may also be directories, in which case their `.root` files are used. | +| `RandomMixBackgroundFiles` | `false` | Treat each file in a background group as an independent pseudo-event source and pick a random file for every overlaid pseudo-event (one-event-per-file mixing) | +| `MergeMCParticles` | `true` | Merge background MCParticles into the output. If `false`, background particles are not stored: tracker hits keep the momentum of their originating particle instead of a particle link, and calorimeter contributions get an empty particle | +| `OverlayThreads` | `1` | Number of worker threads used to read and decompress background files within a single event (`1` = serial). Only the reading is parallelized; the merge stays serial and in-order, so the result is unchanged and deterministic. Most effective with `RandomMixBackgroundFiles` and many input files | | `NumberBackground` | `[]` | Number of background events to overlay per stream (fixed or Poisson mean) | | `Poisson_random_NOverlay` | `[]` | If true, draw the number of events from a Poisson distribution with mean `NumberBackground` | | `NBunchtrain` | `1` | Number of bunch crossings in the bunch train | @@ -93,3 +96,41 @@ ApplicationMgr( OutputLevel=INFO, ) ``` + +## Random background mixing + +For setups where the background is split across a large number of files, each +containing a single pseudo-event (e.g. Muon Collider beam-induced background), set +`RandomMixBackgroundFiles = True`. Each file in a group is then treated as an +independent event source, and a random number of files (set by NumberBackground) +is chosen for every overlaid event. +`BackgroundFileNames` entries may point at directories, whose `.root` files are +collected automatically: + +```python +overlay.RandomMixBackgroundFiles = True +overlay.BackgroundFileNames = [["/path/to/bib_files/"]] +``` + +## Parallel background reading + +With many large background files the algorithm is dominated by reading and +decompressing them. Set `OverlayThreads` to a value greater than 1 to read and +decompress the background files of a single event on several threads: + +```python +overlay.RandomMixBackgroundFiles = True +overlay.OverlayThreads = 4 +``` + +Only the reading is parallelized. The randomness (which files, how many, in +which bunch crossing) is drawn up front, and the merging of the background hits +into the output collections is always done serially and in the same order, so +the result is **identical and deterministic** regardless of `OverlayThreads`. +Because ROOT I/O is made thread-safe with `ROOT::EnableThreadSafety()`, the +algorithm also remains safe to run under Gaudi's intra-event multithreading; the +per-event parallelism composes with it via the shared task arena. + +The speed-up is largest when reading dominates (many large files, tight time +windows that keep the merge cheap); when the merge is the bottleneck the gain is +correspondingly smaller. diff --git a/k4FWCore/CMakeLists.txt b/k4FWCore/CMakeLists.txt index d19c08e5..13f29a46 100644 --- a/k4FWCore/CMakeLists.txt +++ b/k4FWCore/CMakeLists.txt @@ -39,7 +39,7 @@ gaudi_add_module(k4FWCorePlugins components/Reader.cpp components/UniqueIDGenSvc.cpp components/Writer.cpp - LINK Gaudi::GaudiKernel k4FWCore k4FWCore::k4Interface ROOT::Core ROOT::RIO ROOT::Tree EDM4HEP::edm4hep) + LINK Gaudi::GaudiKernel k4FWCore k4FWCore::k4Interface ROOT::Core ROOT::RIO ROOT::Tree EDM4HEP::edm4hep TBB::tbb) target_include_directories(k4FWCorePlugins PUBLIC $ diff --git a/k4FWCore/components/OverlayTiming.cpp b/k4FWCore/components/OverlayTiming.cpp index 30739929..349b6e59 100644 --- a/k4FWCore/components/OverlayTiming.cpp +++ b/k4FWCore/components/OverlayTiming.cpp @@ -31,11 +31,34 @@ #include "k4FWCore/MetadataUtils.h" #include +#include +#include + +#include +#include +#include +#include +#include +#include #include #include #include +namespace fs = std::filesystem; + +// Returns the .root files contained in a directory (non-recursive). +static std::vector filesInFolder(const std::string& folderPath) { + std::vector files; + for (const auto& entry : fs::directory_iterator(folderPath)) { + if (fs::is_regular_file(entry.path()) && entry.path().extension() == ".root") { + files.push_back(entry.path().string()); + } + } + std::sort(files.begin(), files.end()); + return files; +} + template inline float time_of_flight(const T& pos) { // Returns the time of flight to the radius in ns @@ -58,33 +81,53 @@ StatusCode OverlayTiming::initialize() { error() << "Unable to get UniqueIDGenSvc" << endmsg; } + // Make ROOT's global state safe for concurrent reads. IOSvc already does this, + // but call it defensively so the OverlayThreads read pipeline is safe even if + // OverlayTiming is driven differently. + ROOT::EnableThreadSafety(); + + // Expand any directory entries into their list of .root files. This is + // typically used together with RandomMixBackgroundFiles, where each file is + // an independent pseudo-event source. std::vector> inputFiles; - inputFiles = m_inputFileNames.value(); - // if (m_startWithBackgroundFile >= 0) { - // inputFiles = std::vector(m_inputFileNames.begin() + m_startWithBackgroundFile, - // m_inputFileNames.end()); - // } else { - // inputFiles = m_inputFileNames; - // } - // TODO:: shuffle input files - // std::shuffle(inputFiles.begin(), inputFiles.end(), rng_engine); - - m_bkgEvents = make_unique(inputFiles); - for (auto& val : m_bkgEvents->m_totalNumberOfEvents) { - if (val == 0) { - std::string err = "No events found in the background files"; - for (auto& file : m_inputFileNames.value()) { - err += " " + file[0]; + for (const auto& group : m_inputFileNames.value()) { + std::vector expanded; + for (const auto& entry : group) { + if (fs::is_directory(entry)) { + const auto found = filesInFolder(entry); + expanded.insert(expanded.end(), found.begin(), found.end()); + } else { + expanded.push_back(entry); } - error() << err << endmsg; - return StatusCode::FAILURE; } + inputFiles.push_back(std::move(expanded)); } - if (std::any_of(m_bkgEvents->m_totalNumberOfEvents.begin(), m_bkgEvents->m_totalNumberOfEvents.end(), - [this](const int& val) { return this->m_startWithBackgroundEvent >= val; })) { - throw GaudiException("StartWithBackgroundEvent is larger than the number of events in the background files", name(), - StatusCode::FAILURE); + m_bkgEvents = + make_unique(inputFiles, m_randomMix.value(), m_allowReusingBackgroundFiles.value(), name()); + + // In sequential mode the event counts are known upfront and can be validated + // here. In random-mix mode they are determined lazily on first read, so an + // empty file is reported at that point instead. + if (!m_randomMix) { + for (auto& counts : m_bkgEvents->m_totalNumberOfEvents) { + for (auto& val : counts) { + if (val == 0) { + std::string err = "No events found in the background files"; + for (auto& file : m_inputFileNames.value()) { + err += " " + file[0]; + } + error() << err << endmsg; + return StatusCode::FAILURE; + } + } + if (std::any_of(counts.begin(), counts.end(), [this](const size_t& val) { + return this->m_startWithBackgroundEvent >= static_cast(val); + })) { + throw GaudiException("StartWithBackgroundEvent is larger than the number of events in the background files", + name(), StatusCode::FAILURE); + } + } } if (m_Noverlay.empty()) { @@ -104,6 +147,161 @@ StatusCode OverlayTiming::initialize() { return StatusCode::SUCCESS; } +void OverlayTiming::mergeBackgroundFrame( + const podio::Frame& backgroundEvent, float timeOffset, int BX_number_in_train, int physBX, + const std::vector& simTrackerHits, + const std::vector& simCaloHits, + edm4hep::MCParticleCollection& oparticles, std::vector& osimTrackerHits, + std::map>& cellIDsMap, + std::vector& ocaloHitContribs) const { + const auto availableCollections = backgroundEvent.getAvailableCollections(); + + if (std::find(availableCollections.begin(), availableCollections.end(), m_MCParticleCollectionName) == + availableCollections.end()) { + warning() << "Collection " << m_MCParticleCollectionName << " not found in background event" << endmsg; + } + + // To fix the relations we will need to have a map from old to new particle index + std::map oldToNewMap; + std::map, std::vector>> parentDaughterMap; + + if (m_mergeMCParticles) { + const auto& bgParticles = backgroundEvent.get(m_MCParticleCollectionName); + int j = oparticles.size(); + for (size_t i = 0; i < bgParticles.size(); ++i) { + auto npart = bgParticles[i].clone(false); + + npart.setTime(bgParticles[i].getTime() + timeOffset); + npart.setOverlay(true); + oparticles.push_back(npart); + for (const auto& parent : bgParticles[i].getParents()) { + parentDaughterMap[j].first.push_back(parent.getObjectID().index); + } + for (const auto& daughter : bgParticles[i].getDaughters()) { + parentDaughterMap[j].second.push_back(daughter.getObjectID().index); + } + oldToNewMap[i] = j; + j++; + } + for (const auto& [index, parentsDaughters] : parentDaughterMap) { + const auto& [parents, daughters] = parentsDaughters; + for (const auto& parent : parents) { + if (parentDaughterMap.find(oldToNewMap[parent]) == parentDaughterMap.end()) { + continue; + } + oparticles[index].addToParents(oparticles[oldToNewMap[parent]]); + } + for (const auto& daughter : daughters) { + if (parentDaughterMap.find(oldToNewMap[daughter]) == parentDaughterMap.end()) { + continue; + } + oparticles[index].addToDaughters(oparticles[oldToNewMap[daughter]]); + } + } + } + + for (size_t i = 0; i < simTrackerHits.size(); ++i) { + const auto name = inputLocations(SIMTRACKERHIT_INDEX_POSITION)[i]; + debug() << "Processing collection " << name << endmsg; + if (std::find(availableCollections.begin(), availableCollections.end(), name) == availableCollections.end()) { + warning() << "Collection " << name << " not found in background event" << endmsg; + continue; + } + const auto [this_start, this_stop] = define_time_windows(name); + // There are only contributions to the readout if the hits are in the integration window + if (this_stop <= (BX_number_in_train - physBX) * m_deltaT) { + info() << "Skipping collection " << name << " as it is not in the integration window" << endmsg; + continue; + } + auto& ocoll = osimTrackerHits[i]; + for (const auto&& simTrackerHit : backgroundEvent.get(name)) { + const float tof = time_of_flight(simTrackerHit.getPosition()); + + if (!((simTrackerHit.getTime() + timeOffset > this_start + tof) && + (simTrackerHit.getTime() + timeOffset < this_stop + tof))) { + continue; + } + auto nhit = simTrackerHit.clone(false); + nhit.setOverlay(true); + nhit.setTime(simTrackerHit.getTime() + timeOffset); + if (m_mergeMCParticles) { + nhit.setParticle(oparticles[oldToNewMap[simTrackerHit.getParticle().getObjectID().index]]); + } else { + edm4hep::MCParticle mcp = simTrackerHit.getParticle(); + if (mcp.isAvailable()) { + // Preserve the momentum of the originating particle + edm4hep::Vector3d mom = mcp.getMomentum(); + nhit.setMomentum({(float)mom.x, (float)mom.y, (float)mom.z}); + } + } + ocoll.push_back(nhit); + } + } + + for (size_t i = 0; i < simCaloHits.size(); ++i) { + const auto name = inputLocations(SIMCALOHIT_INDEX_POSITION)[i]; + debug() << "Processing collection " << name << endmsg; + if (std::find(availableCollections.begin(), availableCollections.end(), name) == availableCollections.end()) { + warning() << "Collection " << name << " not found in background event" << endmsg; + continue; + } + const auto [this_start, this_stop] = define_time_windows(name); + // There are only contributions to the readout if the hits are in the integration window + if (this_stop <= (BX_number_in_train - physBX) * m_deltaT) { + info() << "Skipping collection " << name << " as it is not in the integration window" << endmsg; + continue; + } + + auto& calHitMap = cellIDsMap[i]; + auto& calHitContribs = ocaloHitContribs[i]; + for (const auto&& simCaloHit : backgroundEvent.get(name)) { + if (calHitMap.find(simCaloHit.getCellID()) == calHitMap.end()) { + // There is no hit at this position. The new hit can be added, if it is not outside the window + auto calhit = edm4hep::MutableSimCalorimeterHit(); + bool add = false; + for (const auto& contrib : simCaloHit.getContributions()) { + if ((contrib.getTime() + timeOffset > this_start) && (contrib.getTime() + timeOffset < this_stop)) { + add = true; + // TODO: Make sure a contribution is not added twice + auto newContrib = contrib.clone(false); + if (m_mergeMCParticles) { + newContrib.setParticle(oparticles[oldToNewMap[contrib.getParticle().getObjectID().index]]); + } else { + newContrib.setParticle(edm4hep::MCParticle()); + } + newContrib.setTime(contrib.getTime() + timeOffset); + calhit.addToContributions(newContrib); + calHitContribs.push_back(newContrib); + } + } + if (add) { + calhit.setCellID(simCaloHit.getCellID()); + calhit.setEnergy(simCaloHit.getEnergy()); + calhit.setPosition(simCaloHit.getPosition()); + calHitMap[calhit.getCellID()] = calhit; + } + } else { + // there is already a hit at this position + auto& calhit = calHitMap[simCaloHit.getCellID()]; + for (const auto& contrib : simCaloHit.getContributions()) { + if ((contrib.getTime() + timeOffset > this_start) && (contrib.getTime() + timeOffset < this_stop)) { + // TODO: Make sure a contribution is not added twice + auto newContrib = contrib.clone(false); + if (m_mergeMCParticles) { + newContrib.setParticle(oparticles[oldToNewMap[contrib.getParticle().getObjectID().index]]); + } else { + newContrib.setParticle(edm4hep::MCParticle()); + } + newContrib.setTime(contrib.getTime() + timeOffset); + calhit.addToContributions(newContrib); + calHitContribs.push_back(newContrib); + } + } + } + } + } +} + retType OverlayTiming::operator()(const edm4hep::EventHeaderCollection& headers, const edm4hep::MCParticleCollection& particles, const std::vector& simTrackerHits, @@ -187,6 +385,20 @@ retType OverlayTiming::operator()(const edm4hep::EventHeaderCollection& headers, } auto physBX = m_physBX.value(); + + // Phase 1: draw all randomness and build the ordered list of background reads. + // Building it up front keeps the RNG sequence -- and therefore the result -- + // identical whether the reads are later executed serially or in parallel. + struct BkgRead { + int group; + int fileIndex; + size_t entry; + float timeOffset; + int bxNumber; + int physBX; + }; + std::vector reads; + // Iterate over each group of files and parameters for (size_t groupIndex = 0; groupIndex < m_bkgEvents->size(); groupIndex++) { if (m_randomBX) { @@ -204,19 +416,31 @@ retType OverlayTiming::operator()(const edm4hep::EventHeaderCollection& headers, } std::shuffle(permutation.begin(), permutation.end(), rng_engine); + // In random-mix mode pick which files of the group to read from, at random. + // In sequential mode the file index is ignored, so this stays trivial. + std::vector fileIndices(m_bkgEvents->m_fileNames[groupIndex].size()); + std::iota(fileIndices.begin(), fileIndices.end(), 0); + if (m_randomMix) { + std::shuffle(fileIndices.begin(), fileIndices.end(), rng_engine); + } + // TODO: Check that there is anything to overlay - debug() << "Starting overlay at event: " << m_bkgEvents->m_nextEntry[groupIndex] << " for the background group " - << groupIndex << endmsg; + debug() << "Starting overlay for the background group " << groupIndex << endmsg; if (m_startWithBackgroundEvent >= 0) { info() << "Skipping to event: " << m_startWithBackgroundEvent << endmsg; - for (auto& entry : m_bkgEvents->m_nextEntry) { + for (auto& entry : m_bkgEvents->m_nextEntry[groupIndex]) { entry = m_startWithBackgroundEvent; } } - // Overlay the background events to each bunchcrossing in the bunch train + // Overlay the background events to each bunchcrossing in the bunch train. + // The file cursor is deliberately declared outside the BX loop: it has to + // keep advancing across bunch crossings, otherwise every BX would restart + // at the front of the permutation and reuse the same file for the whole + // train (which is what happens for NumberBackground = 1). + size_t fileCursor = 0; for (int bxInTrain = 0; bxInTrain < m_NBunchTrain; ++bxInTrain) { const int BX_number_in_train = permutation.at(bxInTrain); @@ -231,153 +455,79 @@ retType OverlayTiming::operator()(const edm4hep::EventHeaderCollection& headers, debug() << "Will overlay " << NOverlay_to_this_BX << " events to BX number " << BX_number_in_train + physBX << endmsg; + if (m_randomMix && fileIndices.empty()) { + warning() << "No background files available for group " << groupIndex << ", skipping overlay" << endmsg; + continue; + } + const float timeOffset = BX_number_in_train * m_deltaT; + reads.reserve(reads.size() + NOverlay_to_this_BX); for (int k = 0; k < NOverlay_to_this_BX; ++k) { - info() << "Overlaying background event " << m_bkgEvents->m_nextEntry[groupIndex] << " from group " << groupIndex - << " to BX " << bxInTrain << endmsg; - if (m_bkgEvents->m_nextEntry[groupIndex] >= m_bkgEvents->m_totalNumberOfEvents[groupIndex] && - !m_allowReusingBackgroundFiles) { - throw GaudiException("No more events in background file", name(), StatusCode::FAILURE); - } - const auto backgroundEvent = - m_bkgEvents->m_rootFileReaders[groupIndex].readEvent(m_bkgEvents->m_nextEntry[groupIndex]); - m_bkgEvents->m_nextEntry[groupIndex]++; - m_bkgEvents->m_nextEntry[groupIndex] %= m_bkgEvents->m_totalNumberOfEvents[groupIndex]; - const auto availableCollections = backgroundEvent.getAvailableCollections(); - - // Either 0 or negative - const auto timeOffset = BX_number_in_train * m_deltaT; - - if (std::find(availableCollections.begin(), availableCollections.end(), m_MCParticleCollectionName) == - availableCollections.end()) { - warning() << "Collection " << m_MCParticleCollectionName << " not found in background event" << endmsg; - } - - // To fix the relations we will need to have a map from old to new particle index - std::map oldToNewMap; - std::map, std::vector>> parentDaughterMap; - - const auto& bgParticles = backgroundEvent.get(m_MCParticleCollectionName); - int j = oparticles.size(); - for (size_t i = 0; i < bgParticles.size(); ++i) { - auto npart = bgParticles[i].clone(false); - - npart.setTime(bgParticles[i].getTime() + timeOffset); - npart.setOverlay(true); - oparticles.push_back(npart); - for (const auto& parent : bgParticles[i].getParents()) { - parentDaughterMap[j].first.push_back(parent.getObjectID().index); - } - for (const auto& daughter : bgParticles[i].getDaughters()) { - parentDaughterMap[j].second.push_back(daughter.getObjectID().index); - } - oldToNewMap[i] = j; - j++; - } - for (const auto& [index, parentsDaughters] : parentDaughterMap) { - const auto& [parents, daughters] = parentsDaughters; - for (const auto& parent : parents) { - if (parentDaughterMap.find(oldToNewMap[parent]) == parentDaughterMap.end()) { - // warning() << "Parent " << parent << " not found in background event" << endmsg; - continue; - } - oparticles[index].addToParents(oparticles[oldToNewMap[parent]]); - } - for (const auto& daughter : daughters) { - if (parentDaughterMap.find(oldToNewMap[daughter]) == parentDaughterMap.end()) { - // warning() << "Parent " << daughter << " not found in background event" << endmsg; - continue; - } - // info() << "Adding (daughter) " << daughter << " to " << index << endmsg; - oparticles[index].addToDaughters(oparticles[oldToNewMap[daughter]]); - } - } - - for (size_t i = 0; i < simTrackerHits.size(); ++i) { - const auto name = inputLocations(SIMTRACKERHIT_INDEX_POSITION)[i]; - debug() << "Processing collection " << name << endmsg; - if (std::find(availableCollections.begin(), availableCollections.end(), name) == availableCollections.end()) { - warning() << "Collection " << name << " not found in background event" << endmsg; - continue; - } - const auto [this_start, this_stop] = define_time_windows(name); - // There are only contributions to the readout if the hits are in the integration window - if (this_stop <= (BX_number_in_train - physBX) * m_deltaT) { - info() << "Skipping collection " << name << " as it is not in the integration window" << endmsg; - continue; - } - auto& ocoll = osimTrackerHits[i]; - for (const auto&& simTrackerHit : backgroundEvent.get(name)) { - const float tof = time_of_flight(simTrackerHit.getPosition()); - - if (!((simTrackerHit.getTime() + timeOffset > this_start + tof) && - (simTrackerHit.getTime() + timeOffset < this_stop + tof))) { - continue; - } - auto nhit = simTrackerHit.clone(false); - nhit.setOverlay(true); - nhit.setTime(simTrackerHit.getTime() + timeOffset); - nhit.setParticle(oparticles[oldToNewMap[simTrackerHit.getParticle().getObjectID().index]]); - ocoll.push_back(nhit); - } - } - - for (size_t i = 0; i < simCaloHits.size(); ++i) { - const auto name = inputLocations(SIMCALOHIT_INDEX_POSITION)[i]; - debug() << "Processing collection " << name << endmsg; - if (std::find(availableCollections.begin(), availableCollections.end(), name) == availableCollections.end()) { - warning() << "Collection " << name << " not found in background event" << endmsg; - continue; - } - const auto [this_start, this_stop] = define_time_windows(name); - // There are only contributions to the readout if the hits are in the integration window - if (this_stop <= (BX_number_in_train - physBX) * m_deltaT) { - info() << "Skipping collection " << name << " as it is not in the integration window" << endmsg; - continue; - } - - auto& calHitMap = cellIDsMap[i]; - auto& calHitContribs = ocaloHitContribs[i]; - for (const auto&& simCaloHit : backgroundEvent.get(name)) { - if (calHitMap.find(simCaloHit.getCellID()) == calHitMap.end()) { - // There is no hit at this position. The new hit can be added, if it is not outside the window - auto calhit = edm4hep::MutableSimCalorimeterHit(); - bool add = false; - for (const auto& contrib : simCaloHit.getContributions()) { - if ((contrib.getTime() + timeOffset > this_start) && (contrib.getTime() + timeOffset < this_stop)) { - add = true; - // TODO: Make sure a contribution is not added twice - auto newContrib = contrib.clone(false); - newContrib.setParticle(oparticles[oldToNewMap[contrib.getParticle().getObjectID().index]]); - newContrib.setTime(contrib.getTime() + timeOffset); - calhit.addToContributions(newContrib); - calHitContribs.push_back(newContrib); - } - } - if (add) { - calhit.setCellID(simCaloHit.getCellID()); - calhit.setEnergy(simCaloHit.getEnergy()); - calhit.setPosition(simCaloHit.getPosition()); - calHitMap[calhit.getCellID()] = calhit; - } - } else { - // there is already a hit at this position - auto& calhit = calHitMap[simCaloHit.getCellID()]; - for (const auto& contrib : simCaloHit.getContributions()) { - if ((contrib.getTime() + timeOffset > this_start) && (contrib.getTime() + timeOffset < this_stop)) { - // TODO: Make sure a contribution is not added twice - auto newContrib = contrib.clone(false); - newContrib.setParticle(oparticles[oldToNewMap[contrib.getParticle().getObjectID().index]]); - newContrib.setTime(contrib.getTime() + timeOffset); - calhit.addToContributions(newContrib); - calHitContribs.push_back(newContrib); - } - } - } + // In random-mix mode walk the shuffled permutation so that consecutive + // overlaid events draw a distinct set of files. + // Once the permutation is exhausted it is reshuffled, so every + // pass is an independent random set instead of a replay of the same + // order. In sequential mode the file index is ignored. + int fileIndex = 0; + if (m_randomMix) { + if (fileCursor == fileIndices.size()) { + std::shuffle(fileIndices.begin(), fileIndices.end(), rng_engine); + fileCursor = 0; } + fileIndex = fileIndices[fileCursor++]; } + const size_t entry = m_bkgEvents->reserve(groupIndex, fileIndex); + reads.push_back({static_cast(groupIndex), fileIndex, entry, timeOffset, BX_number_in_train, physBX}); } } } + + // Phase 2 + 3: read (and decompress) the background frames -- optionally on + // several threads -- and merge them into the outputs in list order. The merge + // is always serial and in-order, so the result is independent of OverlayThreads. + if (m_overlayThreads <= 1) { + for (const auto& r : reads) { + const auto backgroundEvent = m_bkgEvents->readAt(r.group, r.fileIndex, r.entry); + mergeBackgroundFrame(backgroundEvent, r.timeOffset, r.bxNumber, r.physBX, simTrackerHits, simCaloHits, oparticles, + osimTrackerHits, cellIDsMap, ocaloHitContribs); + } + } else { + // Parallel read + decompress (bounded by the number of in-flight tokens), + // serial in-order merge. + struct Item { + const BkgRead* read = nullptr; + podio::Frame frame; + }; + std::atomic nextRead{0}; + const auto ntokens = static_cast(m_overlayThreads.value()); + tbb::parallel_pipeline( + ntokens, tbb::make_filter>(tbb::filter_mode::serial_in_order, + [&](tbb::flow_control& fc) -> std::shared_ptr { + const size_t idx = nextRead++; + if (idx >= reads.size()) { + fc.stop(); + return {}; + } + return std::make_shared(Item{&reads[idx], {}}); + }) & + tbb::make_filter, std::shared_ptr>( + tbb::filter_mode::parallel, + [&](const std::shared_ptr& item) -> std::shared_ptr { + const auto& r = *item->read; + item->frame = m_bkgEvents->readAt(r.group, r.fileIndex, r.entry); + // Force decompression/materialization here (in parallel) so the + // serial merge only touches already in-memory data. + for (const auto& name : item->frame.getAvailableCollections()) { + item->frame.get(name); + } + return item; + }) & + tbb::make_filter, void>( + tbb::filter_mode::serial_in_order, [&](const std::shared_ptr& item) { + const auto& r = *item->read; + mergeBackgroundFrame(item->frame, r.timeOffset, r.bxNumber, r.physBX, simTrackerHits, + simCaloHits, oparticles, osimTrackerHits, cellIDsMap, ocaloHitContribs); + })); + } // Move the SimCalorimeterHitCollections to the output vector // So far they are stored in a map with the cellID as key // but they don't belong to any collection yet diff --git a/k4FWCore/components/OverlayTiming.h b/k4FWCore/components/OverlayTiming.h index d507f73b..f1bb5645 100644 --- a/k4FWCore/components/OverlayTiming.h +++ b/k4FWCore/components/OverlayTiming.h @@ -44,35 +44,112 @@ #include "k4FWCore/Transformer.h" #include "k4Interface/IUniqueIDGenSvc.h" +#include "GaudiKernel/GaudiException.h" + // Needed for some of the more complex properties #include "Gaudi/Parsers/Factory.h" #include "Gaudi/Property.h" #include +#include #include #include +// Holds the background events and provides thread-safe reads. ROOT TFile access +// is made safe process-wide by ROOT::EnableThreadSafety() (called in +// initialize()). In random-mix mode every read opens its own reader, so reads +// of different files proceed concurrently -- this is what lets OverlayThreads +// parallelize the (I/O-dominated) background reading. In sequential mode the +// shared per-group reader is used under a mutex. +// +// Two source strategies are supported, selected by m_randomMix: +// * sequential (default): each group is read as one logical stream through a +// persistent reader, advancing an internal cursor; +// * random mix: each file in a group is an independent event source, opened +// on demand, so the caller can pick a random file per overlay. +// +// Bookkeeping is stored per [group][file]. In sequential mode the file +// dimension has a single slot per group; in random-mix mode the event count is +// left at 0 ("unknown") and determined lazily at read time. struct EventHolder { std::vector> m_fileNames; + bool m_randomMix{false}; + bool m_allowReuse{false}; + std::string m_algName; + + // Sequential mode only: one persistent reader per group. std::vector m_rootFileReaders; - std::vector m_totalNumberOfEvents; - std::map m_events; - std::vector m_nextEntry; + // [group][file]. In sequential mode the inner vector has a single element. + std::vector> m_totalNumberOfEvents; + std::vector> m_nextEntry; - EventHolder(const std::vector>& fileNames) : m_fileNames(fileNames) { - for (auto& names : m_fileNames) { - m_rootFileReaders.emplace_back(podio::makeReader(names)); - m_totalNumberOfEvents.push_back(m_rootFileReaders.back().getEntries("events")); + // Guards the cursors and the shared sequential-mode readers. + std::mutex m_ioMutex; + + EventHolder(const std::vector>& fileNames, bool randomMix, bool allowReuse, + const std::string& algName) + : m_fileNames(fileNames), m_randomMix(randomMix), m_allowReuse(allowReuse), m_algName(algName) { + m_totalNumberOfEvents.resize(m_fileNames.size()); + m_nextEntry.resize(m_fileNames.size()); + for (size_t group = 0; group < m_fileNames.size(); ++group) { + if (m_randomMix) { + // One independent event source per file; counts are determined lazily. + m_totalNumberOfEvents[group].resize(m_fileNames[group].size(), 0); + m_nextEntry[group].resize(m_fileNames[group].size(), 0); + } else { + // The whole group is read as a single logical stream. + m_rootFileReaders.emplace_back(podio::makeReader(m_fileNames[group])); + m_totalNumberOfEvents[group].push_back(m_rootFileReaders.back().getEntries("events")); + m_nextEntry[group].push_back(0); + } } - m_nextEntry.resize(m_fileNames.size(), 0); } EventHolder() = default; - // TODO: Cache functionality - // podio::Frame& read - size_t size() const { return m_fileNames.size(); } + + // Advance the cursor for (group, file) and return the raw entry to read. + // Cheap and I/O-free, so calling it serially (during the work-list build) + // does not limit read parallelism. + size_t reserve(int group, int file) { + std::lock_guard lock(m_ioMutex); + size_t& entry = m_nextEntry[group][file]; + const size_t e = entry; + const size_t total = m_totalNumberOfEvents[group][file]; + entry = (total > 0) ? (entry + 1) % total : entry + 1; // wrap once the total is known + return e; + } + + // Read a specific (group, file, rawEntry). Thread-safe: in random-mix mode + // each call opens its own reader and can run concurrently; in sequential mode + // the shared per-group reader is used under the mutex. + podio::Frame readAt(int group, int file, size_t rawEntry) { + if (m_randomMix) { + podio::Reader reader = podio::makeReader(m_fileNames[group][file]); + const size_t total = reader.getEntries("events"); + if (total == 0) { + throw GaudiException("No events found in background file " + m_fileNames[group][file], m_algName, + StatusCode::FAILURE); + } + if (rawEntry >= total && !m_allowReuse) { + throw GaudiException("No more events in background file", m_algName, StatusCode::FAILURE); + } + return reader.readEvent(rawEntry % total); + } + std::lock_guard lock(m_ioMutex); + const size_t total = m_totalNumberOfEvents[group][0]; + if (rawEntry >= total && !m_allowReuse) { + throw GaudiException("No more events in background file", m_algName, StatusCode::FAILURE); + } + return m_rootFileReaders[group].readEvent(rawEntry % total); + } + + // Serial convenience: reserve + read in one call (used by the serial path). + podio::Frame getFrame(int group, int fileIndex) { + const int file = m_randomMix ? fileIndex : 0; + return readAt(group, file, reserve(group, file)); + } }; using retType = @@ -106,6 +183,18 @@ struct OverlayTiming : public k4FWCore::MultiTransformer define_time_windows(const std::string& Collection_name) const; + // Merge one already-read background frame into the output accumulators. This + // is the per-pseudo-event work; it is called serially and in-order (both by + // the serial path and by the in-order output stage of the parallel pipeline), + // so the result is independent of OverlayThreads. + void mergeBackgroundFrame(const podio::Frame& backgroundEvent, float timeOffset, int BX_number_in_train, int physBX, + const std::vector& simTrackerHits, + const std::vector& simCaloHits, + edm4hep::MCParticleCollection& oparticles, + std::vector& osimTrackerHits, + std::map>& cellIDsMap, + std::vector& ocaloHitContribs) const; + private: // These correspond to the index position in the argument list constexpr static int SIMTRACKERHIT_INDEX_POSITION = 2; @@ -149,6 +238,25 @@ struct OverlayTiming : public k4FWCore::MultiTransformer m_copyCellIDMetadata{this, "CopyCellIDMetadata", false, "Copy cell ID encoding metadata from input to output collections"}; + Gaudi::Property m_randomMix{ + this, "RandomMixBackgroundFiles", false, + "Treat each file in a background group as an independent (pseudo-)event source and pick a random file for every " + "overlaid event (one-event-per-file mixing). Entries of BackgroundFileNames may also be directories, " + "whose .root files are used."}; + + Gaudi::Property m_mergeMCParticles{ + this, "MergeMCParticles", true, + "Merge the background MCParticle collection into the output. If false, background particles are not " + "stored: tracker hits keep the momentum of their originating particle instead of a particle link, and " + "calorimeter contributions get an empty particle."}; + + Gaudi::Property m_overlayThreads{ + this, "OverlayThreads", 1, + "Number of worker threads used to read and decompress background files within a single event (1 = " + "serial, current behaviour). Only the (I/O-dominated) reading is parallelized; the merge into the " + "output collections stays serial and in-order, so results are unchanged and deterministic. Most " + "effective with RandomMixBackgroundFiles and many input files."}; + // Gaudi::Property m_maxCachedFrames{ // this, "MaxCachedFrames", 0, "Maximum number of frames cached from background files"}; diff --git a/test/k4FWCoreTest/CMakeLists.txt b/test/k4FWCoreTest/CMakeLists.txt index 2cd701ae..6fc90e54 100644 --- a/test/k4FWCoreTest/CMakeLists.txt +++ b/test/k4FWCoreTest/CMakeLists.txt @@ -272,6 +272,8 @@ add_test_fwcore(InvalidOutputCommandsNoCrash options/invalidOutputCommandsNoCras set_tests_properties(InvalidOutputCommandsNoCrash PROPERTIES FIXTURES_REQUIRED ProducerFile PASS_REGULAR_EXPRESSION "ERROR 'abc' is not a valid command for the KeepDropSwitch") add_test_fwcore(OverlayTiming options/TestOverlayTiming.py ADD_TO_CHECK_FILES PROPERTIES FIXTURES_REQUIRED ProducerMultipleFile) +add_test_fwcore(OverlayTimingRandomMix options/TestOverlayTimingRandomMix.py ADD_TO_CHECK_FILES PROPERTIES FIXTURES_REQUIRED ProducerMultipleFile) +add_test_fwcore(OverlayTimingMT options/TestOverlayTimingMT.py ADD_TO_CHECK_FILES PROPERTIES FIXTURES_REQUIRED ProducerMultipleFile) add_test(NAME check_broken_pipe COMMAND bash -c "[ $(${K4RUN} options/ExampleFunctionalProducer.py | head -n 2 | wc -l) = 2 ]" diff --git a/test/k4FWCoreTest/options/TestOverlayTimingMT.py b/test/k4FWCoreTest/options/TestOverlayTimingMT.py new file mode 100644 index 00000000..54c8f8b4 --- /dev/null +++ b/test/k4FWCoreTest/options/TestOverlayTimingMT.py @@ -0,0 +1,77 @@ +# +# Copyright (c) 2014-2024 Key4hep-Project. +# +# This file is part of Key4hep. +# See https://key4hep.github.io/key4hep-doc/ for further info. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# + +# Runs the OverlayTiming algorithm with intra-event multithreading (multiple +# event slots processed concurrently by the Avalanche scheduler). This exercises +# the const operator() from several worker threads sharing a single background +# EventHolder, checking that its worker-thread ROOT I/O serialization is safe. + +from Gaudi.Configuration import INFO, WARNING +from Configurables import EventDataSvc, EventHeaderCreator, OverlayTiming, UniqueIDGenSvc +from Configurables import HiveSlimEventLoopMgr, HiveWhiteBoard, AvalancheSchedulerSvc +from k4FWCore import ApplicationMgr, IOSvc + +evtslots = 6 +threads = 6 + +uid_svc = UniqueIDGenSvc("UniqueIDGenSvc") + +whiteboard = HiveWhiteBoard( + "EventDataSvc", + EventSlots=evtslots, + ForceLeaves=True, +) +slimeventloopmgr = HiveSlimEventLoopMgr( + "HiveSlimEventLoopMgr", + SchedulerName="AvalancheSchedulerSvc", + OutputLevel=WARNING, +) +scheduler = AvalancheSchedulerSvc(ThreadPoolSize=threads, OutputLevel=WARNING) + +iosvc = IOSvc("IOSvc") +iosvc.Input = "functional_producer_multiple.root" +iosvc.Output = "overlay_mt_output.root" + +header = EventHeaderCreator("EventHeaderCreator") + +overlay = OverlayTiming("OverlayTiming") +overlay.MCParticles = "MCParticles1" +overlay.SimTrackerHits = ["SimTrackerHits"] +overlay.SimCalorimeterHits = [] +overlay.OutputMCParticles = "OverlayMCParticles" +overlay.OutputSimTrackerHits = ["OverlaySimTrackerHits"] +overlay.OutputSimCalorimeterHits = [] +overlay.OutputCaloHitContributions = [] +overlay.BackgroundMCParticleCollectionName = "MCParticles1" +overlay.BackgroundFileNames = [["functional_producer_multiple.root"]] +overlay.NumberBackground = [1] +overlay.Poisson_random_NOverlay = [False] +overlay.NBunchtrain = 3 +overlay.TimeWindows = {"SimTrackerHits": [-10000, 10000]} +overlay.AllowReusingBackgroundFiles = True + +ApplicationMgr( + TopAlg=[header, overlay], + EvtSel="NONE", + EvtMax=-1, + ExtSvc=[whiteboard, uid_svc], + EventLoop=slimeventloopmgr, + MessageSvcType="InertMessageSvc", + OutputLevel=INFO, +) diff --git a/test/k4FWCoreTest/options/TestOverlayTimingRandomMix.py b/test/k4FWCoreTest/options/TestOverlayTimingRandomMix.py new file mode 100644 index 00000000..be807602 --- /dev/null +++ b/test/k4FWCoreTest/options/TestOverlayTimingRandomMix.py @@ -0,0 +1,65 @@ +# +# Copyright (c) 2014-2024 Key4hep-Project. +# +# This file is part of Key4hep. +# See https://key4hep.github.io/key4hep-doc/ for further info. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# + +# Tests the OverlayTiming algorithm in random-mix mode, using +# functional_producer_multiple.root as both signal and background. Each +# background file is treated as an independent event source and picked at +# random, MCParticles are not merged (MergeMCParticles = False), and the +# background files are read on 2 threads (OverlayThreads = 2). + +from Gaudi.Configuration import INFO +from Configurables import EventDataSvc, EventHeaderCreator, OverlayTiming, UniqueIDGenSvc +from k4FWCore import ApplicationMgr, IOSvc + +uid_svc = UniqueIDGenSvc("UniqueIDGenSvc") + +iosvc = IOSvc("IOSvc") +iosvc.Input = "functional_producer_multiple.root" +iosvc.Output = "overlay_randommix_output.root" + +header = EventHeaderCreator("EventHeaderCreator") + +overlay = OverlayTiming("OverlayTiming") +overlay.MCParticles = "MCParticles1" +overlay.SimTrackerHits = ["SimTrackerHits"] +overlay.SimCalorimeterHits = [] +overlay.OutputMCParticles = "OverlayMCParticles" +overlay.OutputSimTrackerHits = ["OverlaySimTrackerHits"] +overlay.OutputSimCalorimeterHits = [] +overlay.OutputCaloHitContributions = [] +overlay.BackgroundMCParticleCollectionName = "MCParticles1" +overlay.BackgroundFileNames = [["functional_producer_multiple.root"]] +overlay.NumberBackground = [1] +overlay.Poisson_random_NOverlay = [False] +overlay.NBunchtrain = 1 +overlay.TimeWindows = {"SimTrackerHits": [-10000, 10000]} +# Random-mix specific options +overlay.RandomMixBackgroundFiles = True +overlay.MergeMCParticles = False +overlay.AllowReusingBackgroundFiles = True +# Exercise the parallel background reading (result is independent of this) +overlay.OverlayThreads = 2 + +ApplicationMgr( + TopAlg=[header, overlay], + EvtSel="NONE", + EvtMax=3, + ExtSvc=[EventDataSvc("EventDataSvc"), uid_svc], + OutputLevel=INFO, +)