From 5c76efae21cd25031f57e344c291447426ea1d81 Mon Sep 17 00:00:00 2001 From: Federico Meloni Date: Mon, 6 Jul 2026 12:57:24 +0200 Subject: [PATCH 01/15] merging randomMix functionalities into overlaytiming --- doc/OverlayTiming.md | 21 +- doc/ReleaseNotes.md | 5 + k4FWCore/components/OverlayTiming.cpp | 192 +++++++++++------- k4FWCore/components/OverlayTiming.h | 154 +++++++++++++- test/k4FWCoreTest/CMakeLists.txt | 2 + .../options/TestOverlayTimingMT.py | 77 +++++++ .../options/TestOverlayTimingRandomMix.py | 62 ++++++ 7 files changed, 434 insertions(+), 79 deletions(-) create mode 100644 test/k4FWCoreTest/options/TestOverlayTimingMT.py create mode 100644 test/k4FWCoreTest/options/TestOverlayTimingRandomMix.py diff --git a/doc/OverlayTiming.md b/doc/OverlayTiming.md index 320edccb..f388dbab 100644 --- a/doc/OverlayTiming.md +++ b/doc/OverlayTiming.md @@ -33,7 +33,9 @@ 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 event source and pick a random file for every overlaid 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 | | `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 +95,20 @@ ApplicationMgr( OutputLevel=INFO, ) ``` + +## Random background mixing + +For setups where the background is split across a large number of files, each +containing a single 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 file 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/"]] +``` + +All background ROOT I/O is serialized on a dedicated worker thread, so the +algorithm may be run with intra-event multithreading. diff --git a/doc/ReleaseNotes.md b/doc/ReleaseNotes.md index 1de9fc8e..4848001f 100644 --- a/doc/ReleaseNotes.md +++ b/doc/ReleaseNotes.md @@ -1,5 +1,10 @@ # v01-06 +* 2026-07-06 Federico Meloni + - Extend `OverlayTiming` with random background-file mixing: the new `RandomMixBackgroundFiles` option treats each file in a background group as an independent event source and picks a random file for every overlaid event, and `BackgroundFileNames` entries may now be directories (their `.root` files are used). + - Add the `MergeMCParticles` option to `OverlayTiming` (default `true`); when `false`, background MCParticles are not stored, tracker hits keep the momentum of their originating particle and calorimeter contributions get an empty particle. + - Serialize all background ROOT I/O on a dedicated worker thread so `OverlayTiming` is safe to run with intra-event multithreading. + * 2026-04-16 Juan Miguel Carceller ([PR#397](https://github.com/key4hep/k4FWCore/pull/397)) - Bump the version of k4FWCore for versioning the usage of new features diff --git a/k4FWCore/components/OverlayTiming.cpp b/k4FWCore/components/OverlayTiming.cpp index c688cdd7..82f8e09c 100644 --- a/k4FWCore/components/OverlayTiming.cpp +++ b/k4FWCore/components/OverlayTiming.cpp @@ -32,10 +32,27 @@ #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 +75,46 @@ StatusCode OverlayTiming::initialize() { error() << "Unable to get UniqueIDGenSvc" << endmsg; } + // 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()) { @@ -202,14 +232,21 @@ 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; } } @@ -229,17 +266,17 @@ 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; + } 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]; + // In random-mix mode pick a random file of the group; in sequential mode + // the file index is ignored. All bounds checking, reading and cursor + // advancing happen thread-safely inside getFrame (serialized ROOT I/O). + const int fileIndex = m_randomMix ? fileIndices[k % fileIndices.size()] : 0; + debug() << "Overlaying background event from group " << groupIndex << " to BX " << bxInTrain << endmsg; + const auto backgroundEvent = m_bkgEvents->getFrame(groupIndex, fileIndex); const auto availableCollections = backgroundEvent.getAvailableCollections(); // Either 0 or negative @@ -254,39 +291,41 @@ retType OverlayTiming::operator()(const edm4hep::EventHeaderCollection& headers, 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; + 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); } - oparticles[index].addToParents(oparticles[oldToNewMap[parent]]); + oldToNewMap[i] = j; + j++; } - for (const auto& daughter : daughters) { - if (parentDaughterMap.find(oldToNewMap[daughter]) == parentDaughterMap.end()) { - // warning() << "Parent " << daughter << " not found in background event" << endmsg; - continue; + 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]]); } - // info() << "Adding (daughter) " << daughter << " to " << index << endmsg; - oparticles[index].addToDaughters(oparticles[oldToNewMap[daughter]]); } } @@ -313,7 +352,16 @@ retType OverlayTiming::operator()(const edm4hep::EventHeaderCollection& headers, auto nhit = simTrackerHit.clone(false); nhit.setOverlay(true); nhit.setTime(simTrackerHit.getTime() + timeOffset); - nhit.setParticle(oparticles[oldToNewMap[simTrackerHit.getParticle().getObjectID().index]]); + 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); } } @@ -344,7 +392,11 @@ retType OverlayTiming::operator()(const edm4hep::EventHeaderCollection& headers, add = true; // TODO: Make sure a contribution is not added twice auto newContrib = contrib.clone(false); - newContrib.setParticle(oparticles[oldToNewMap[contrib.getParticle().getObjectID().index]]); + 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); @@ -363,7 +415,11 @@ retType OverlayTiming::operator()(const edm4hep::EventHeaderCollection& headers, 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]]); + 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); diff --git a/k4FWCore/components/OverlayTiming.h b/k4FWCore/components/OverlayTiming.h index d507f73b..842be20a 100644 --- a/k4FWCore/components/OverlayTiming.h +++ b/k4FWCore/components/OverlayTiming.h @@ -44,35 +44,157 @@ #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 +#include +#include #include +#include #include +// Holds the background events and serializes all ROOT I/O on a single worker +// thread. ROOT TFile access is not thread-safe, so routing every read through +// one worker thread makes a shared EventHolder safe when the (const) overlay +// operator() is scheduled concurrently under intra-event multi-threading. +// +// 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. 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. + // A total of 0 means "not yet determined" (filled lazily in random-mix mode). + std::vector> m_totalNumberOfEvents; + std::vector> m_nextEntry; + + // Single worker thread serializing all ROOT I/O. + struct Request { + int groupIndex; + int fileIndex; + std::promise prom; + }; + std::queue m_requests; + std::mutex m_queueMutex; + std::condition_variable m_queueCV; + std::thread m_worker; + bool m_stop{false}; - 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")); + 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); + + m_worker = std::thread([this]() { + while (true) { + Request req; + { + std::unique_lock lock(m_queueMutex); + m_queueCV.wait(lock, [this]() { return m_stop || !m_requests.empty(); }); + if (m_stop && m_requests.empty()) { + return; + } + req = std::move(m_requests.front()); + m_requests.pop(); + } + try { + req.prom.set_value(read(req.groupIndex, req.fileIndex)); + } catch (...) { + req.prom.set_exception(std::current_exception()); + } + } + }); } EventHolder() = default; - // TODO: Cache functionality - // podio::Frame& read + ~EventHolder() { + { + std::lock_guard lock(m_queueMutex); + m_stop = true; + } + m_queueCV.notify_all(); + if (m_worker.joinable()) { + m_worker.join(); + } + } + + // Thread-safe: enqueue a read request and block until the worker fulfills it. + // In sequential mode fileIndex is ignored (always slot 0). + podio::Frame getFrame(int groupIndex, int fileIndex) { + Request req{groupIndex, m_randomMix ? fileIndex : 0, std::promise()}; + auto fut = req.prom.get_future(); + { + std::lock_guard lock(m_queueMutex); + m_requests.push(std::move(req)); + } + m_queueCV.notify_one(); + return fut.get(); + } size_t size() const { return m_fileNames.size(); } + +private: + // Runs exclusively on the worker thread, so cursor bookkeeping needs no extra + // locking. + podio::Frame read(int group, int file) { + size_t& total = m_totalNumberOfEvents[group][file]; + size_t& entry = m_nextEntry[group][file]; + + podio::Reader* reader = nullptr; + std::optional ondemand; + if (m_randomMix) { + ondemand.emplace(podio::makeReader(m_fileNames[group][file])); + reader = &ondemand.value(); + if (total == 0) { + total = reader->getEntries("events"); + } + } else { + reader = &m_rootFileReaders[group]; + } + + if (total == 0) { + throw GaudiException("No events found in background file " + m_fileNames[group][file], m_algName, + StatusCode::FAILURE); + } + if (entry >= total && !m_allowReuse) { + throw GaudiException("No more events in background file", m_algName, StatusCode::FAILURE); + } + podio::Frame frame = reader->readEvent(entry % total); + entry = (entry + 1) % total; + return frame; + } }; using retType = @@ -149,6 +271,18 @@ 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 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_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 8fc9c6e5..7ab35e25 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..7e252138 --- /dev/null +++ b/test/k4FWCoreTest/options/TestOverlayTimingRandomMix.py @@ -0,0 +1,62 @@ +# +# 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, and MCParticles are not merged (MergeMCParticles = False). + +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 + +ApplicationMgr( + TopAlg=[header, overlay], + EvtSel="NONE", + EvtMax=3, + ExtSvc=[EventDataSvc("EventDataSvc"), uid_svc], + OutputLevel=INFO, +) From 54afad5618ee5d3b2dadabc4f77b338d8b26f0b5 Mon Sep 17 00:00:00 2001 From: Federico Meloni Date: Mon, 6 Jul 2026 13:04:50 +0200 Subject: [PATCH 02/15] docs --- doc/OverlayTiming.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/OverlayTiming.md b/doc/OverlayTiming.md index f388dbab..65aa8fb2 100644 --- a/doc/OverlayTiming.md +++ b/doc/OverlayTiming.md @@ -34,7 +34,7 @@ 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. 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 event source and pick a random file for every overlaid event (one-event-per-file mixing) | +| `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 | | `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` | From 4c170509d757049c01e113ce322845ab4313ef8c Mon Sep 17 00:00:00 2001 From: Federico Meloni Date: Mon, 6 Jul 2026 13:06:31 +0200 Subject: [PATCH 03/15] more docs --- doc/OverlayTiming.md | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/doc/OverlayTiming.md b/doc/OverlayTiming.md index 65aa8fb2..5d624b06 100644 --- a/doc/OverlayTiming.md +++ b/doc/OverlayTiming.md @@ -99,9 +99,10 @@ ApplicationMgr( ## Random background mixing For setups where the background is split across a large number of files, each -containing a single event (e.g. Muon Collider beam-induced background), set +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 file is chosen for every overlaid event. +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: From 9ef106394e06bebb950bc94fc6577448246c647c Mon Sep 17 00:00:00 2001 From: Federico Meloni Date: Mon, 6 Jul 2026 13:08:25 +0200 Subject: [PATCH 04/15] undo over-eager release notes --- doc/ReleaseNotes.md | 5 ----- 1 file changed, 5 deletions(-) diff --git a/doc/ReleaseNotes.md b/doc/ReleaseNotes.md index 4848001f..1de9fc8e 100644 --- a/doc/ReleaseNotes.md +++ b/doc/ReleaseNotes.md @@ -1,10 +1,5 @@ # v01-06 -* 2026-07-06 Federico Meloni - - Extend `OverlayTiming` with random background-file mixing: the new `RandomMixBackgroundFiles` option treats each file in a background group as an independent event source and picks a random file for every overlaid event, and `BackgroundFileNames` entries may now be directories (their `.root` files are used). - - Add the `MergeMCParticles` option to `OverlayTiming` (default `true`); when `false`, background MCParticles are not stored, tracker hits keep the momentum of their originating particle and calorimeter contributions get an empty particle. - - Serialize all background ROOT I/O on a dedicated worker thread so `OverlayTiming` is safe to run with intra-event multithreading. - * 2026-04-16 Juan Miguel Carceller ([PR#397](https://github.com/key4hep/k4FWCore/pull/397)) - Bump the version of k4FWCore for versioning the usage of new features From fdfde2743c7ff364c00df045fe81bd0d47b96149 Mon Sep 17 00:00:00 2001 From: Federico Meloni Date: Mon, 6 Jul 2026 13:13:02 +0200 Subject: [PATCH 05/15] phrasing of pseudo-events --- k4FWCore/components/OverlayTiming.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/k4FWCore/components/OverlayTiming.h b/k4FWCore/components/OverlayTiming.h index 842be20a..be29acfe 100644 --- a/k4FWCore/components/OverlayTiming.h +++ b/k4FWCore/components/OverlayTiming.h @@ -273,7 +273,7 @@ struct OverlayTiming : public k4FWCore::MultiTransformer m_randomMix{ this, "RandomMixBackgroundFiles", false, - "Treat each file in a background group as an independent event source and pick a random file for every " + "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."}; From a315f27cbd9a23f266328eb10f9d27bdaaaae8d7 Mon Sep 17 00:00:00 2001 From: Federico Meloni Date: Mon, 6 Jul 2026 14:34:51 +0200 Subject: [PATCH 06/15] use TBB --- CMakeLists.txt | 1 + k4FWCore/CMakeLists.txt | 2 +- k4FWCore/components/OverlayTiming.cpp | 391 ++++++++++++++++---------- k4FWCore/components/OverlayTiming.h | 150 ++++------ 4 files changed, 300 insertions(+), 244 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 29bfe291..0266dcf1 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/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 82f8e09c..cca24bb4 100644 --- a/k4FWCore/components/OverlayTiming.cpp +++ b/k4FWCore/components/OverlayTiming.cpp @@ -31,9 +31,15 @@ #include "k4FWCore/MetadataUtils.h" #include +#include + +#include #include +#include +#include #include +#include #include #include #include @@ -75,6 +81,11 @@ 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. @@ -134,6 +145,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, @@ -215,6 +381,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) { @@ -270,166 +450,67 @@ retType OverlayTiming::operator()(const edm4hep::EventHeaderCollection& headers, warning() << "No background files available for group " << groupIndex << ", skipping overlay" << endmsg; continue; } + const float timeOffset = BX_number_in_train * m_deltaT; for (int k = 0; k < NOverlay_to_this_BX; ++k) { // In random-mix mode pick a random file of the group; in sequential mode - // the file index is ignored. All bounds checking, reading and cursor - // advancing happen thread-safely inside getFrame (serialized ROOT I/O). - const int fileIndex = m_randomMix ? fileIndices[k % fileIndices.size()] : 0; - debug() << "Overlaying background event from group " << groupIndex << " to BX " << bxInTrain << endmsg; - const auto backgroundEvent = m_bkgEvents->getFrame(groupIndex, fileIndex); - 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; - - 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()) { - // 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); - 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; - } + // the file index is ignored. reserve() advances the per-file cursor now + // (serially); the actual ROOT read happens in phase 2, possibly on + // several threads. + const int fileIndex = m_randomMix ? fileIndices[k % fileIndices.size()] : 0; + const size_t entry = m_bkgEvents->reserve(groupIndex, fileIndex); + reads.push_back({static_cast(groupIndex), fileIndex, entry, timeOffset, BX_number_in_train, physBX}); + } + } + } - 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; + // 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 size_t 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 {}; } - } 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()); + return std::make_shared(Item{&reads[idx], {}}); + }) & + tbb::make_filter, std::shared_ptr>( + tbb::filter_mode::parallel, + [&](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); } - newContrib.setTime(contrib.getTime() + timeOffset); - calhit.addToContributions(newContrib); - calHitContribs.push_back(newContrib); - } - } - } - } - } - } - } + return item; + }) & + tbb::make_filter, void>( + tbb::filter_mode::serial_in_order, [&](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 diff --git a/k4FWCore/components/OverlayTiming.h b/k4FWCore/components/OverlayTiming.h index be29acfe..57a51de9 100644 --- a/k4FWCore/components/OverlayTiming.h +++ b/k4FWCore/components/OverlayTiming.h @@ -50,20 +50,17 @@ #include "Gaudi/Parsers/Factory.h" #include "Gaudi/Property.h" -#include -#include #include #include -#include -#include #include -#include #include -// Holds the background events and serializes all ROOT I/O on a single worker -// thread. ROOT TFile access is not thread-safe, so routing every read through -// one worker thread makes a shared EventHolder safe when the (const) overlay -// operator() is scheduled concurrently under intra-event multi-threading. +// 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 @@ -72,7 +69,8 @@ // 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. +// 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}; @@ -83,21 +81,11 @@ struct EventHolder { std::vector m_rootFileReaders; // [group][file]. In sequential mode the inner vector has a single element. - // A total of 0 means "not yet determined" (filled lazily in random-mix mode). std::vector> m_totalNumberOfEvents; std::vector> m_nextEntry; - // Single worker thread serializing all ROOT I/O. - struct Request { - int groupIndex; - int fileIndex; - std::promise prom; - }; - std::queue m_requests; - std::mutex m_queueMutex; - std::condition_variable m_queueCV; - std::thread m_worker; - bool m_stop{false}; + // 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) @@ -116,84 +104,51 @@ struct EventHolder { m_nextEntry[group].push_back(0); } } - - m_worker = std::thread([this]() { - while (true) { - Request req; - { - std::unique_lock lock(m_queueMutex); - m_queueCV.wait(lock, [this]() { return m_stop || !m_requests.empty(); }); - if (m_stop && m_requests.empty()) { - return; - } - req = std::move(m_requests.front()); - m_requests.pop(); - } - try { - req.prom.set_value(read(req.groupIndex, req.fileIndex)); - } catch (...) { - req.prom.set_exception(std::current_exception()); - } - } - }); } EventHolder() = default; - ~EventHolder() { - { - std::lock_guard lock(m_queueMutex); - m_stop = true; - } - m_queueCV.notify_all(); - if (m_worker.joinable()) { - m_worker.join(); - } - } + size_t size() const { return m_fileNames.size(); } - // Thread-safe: enqueue a read request and block until the worker fulfills it. - // In sequential mode fileIndex is ignored (always slot 0). - podio::Frame getFrame(int groupIndex, int fileIndex) { - Request req{groupIndex, m_randomMix ? fileIndex : 0, std::promise()}; - auto fut = req.prom.get_future(); - { - std::lock_guard lock(m_queueMutex); - m_requests.push(std::move(req)); - } - m_queueCV.notify_one(); - return fut.get(); + // 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; } - size_t size() const { return m_fileNames.size(); } - -private: - // Runs exclusively on the worker thread, so cursor bookkeeping needs no extra - // locking. - podio::Frame read(int group, int file) { - size_t& total = m_totalNumberOfEvents[group][file]; - size_t& entry = m_nextEntry[group][file]; - - podio::Reader* reader = nullptr; - std::optional ondemand; + // 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) { - ondemand.emplace(podio::makeReader(m_fileNames[group][file])); - reader = &ondemand.value(); + podio::Reader reader = podio::makeReader(m_fileNames[group][file]); + const size_t total = reader.getEntries("events"); if (total == 0) { - total = reader->getEntries("events"); + throw GaudiException("No events found in background file " + m_fileNames[group][file], m_algName, + StatusCode::FAILURE); } - } else { - reader = &m_rootFileReaders[group]; - } - - 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); } - if (entry >= total && !m_allowReuse) { + 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); } - podio::Frame frame = reader->readEvent(entry % total); - entry = (entry + 1) % total; - return frame; + 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)); } }; @@ -228,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; @@ -283,6 +250,13 @@ struct OverlayTiming : public k4FWCore::MultiTransformer 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"}; From 24b2107291b710bebd1bb72fd86731be54d37a93 Mon Sep 17 00:00:00 2001 From: Federico Meloni Date: Mon, 6 Jul 2026 15:18:55 +0200 Subject: [PATCH 07/15] clang format --- doc/OverlayTiming.md | 2 +- k4FWCore/components/OverlayTiming.cpp | 99 +++++++++++++-------------- k4FWCore/components/OverlayTiming.h | 24 +++---- 3 files changed, 62 insertions(+), 63 deletions(-) diff --git a/doc/OverlayTiming.md b/doc/OverlayTiming.md index 5d624b06..d7a9953e 100644 --- a/doc/OverlayTiming.md +++ b/doc/OverlayTiming.md @@ -101,7 +101,7 @@ ApplicationMgr( 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) +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: diff --git a/k4FWCore/components/OverlayTiming.cpp b/k4FWCore/components/OverlayTiming.cpp index cca24bb4..79092baf 100644 --- a/k4FWCore/components/OverlayTiming.cpp +++ b/k4FWCore/components/OverlayTiming.cpp @@ -103,7 +103,8 @@ StatusCode OverlayTiming::initialize() { inputFiles.push_back(std::move(expanded)); } - m_bkgEvents = make_unique(inputFiles, m_randomMix.value(), m_allowReusingBackgroundFiles.value(), name()); + 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 @@ -120,8 +121,9 @@ StatusCode OverlayTiming::initialize() { return StatusCode::FAILURE; } } - if (std::any_of(counts.begin(), counts.end(), - [this](const size_t& val) { return this->m_startWithBackgroundEvent >= static_cast(val); })) { + 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); } @@ -147,12 +149,11 @@ StatusCode OverlayTiming::initialize() { 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, + const std::vector& simTrackerHits, + const std::vector& simCaloHits, + edm4hep::MCParticleCollection& oparticles, std::vector& osimTrackerHits, std::map>& cellIDsMap, - std::vector& ocaloHitContribs) const { + std::vector& ocaloHitContribs) const { const auto availableCollections = backgroundEvent.getAvailableCollections(); if (std::find(availableCollections.begin(), availableCollections.end(), m_MCParticleCollectionName) == @@ -161,12 +162,12 @@ void OverlayTiming::mergeBackgroundFrame( } // To fix the relations we will need to have a map from old to new particle index - std::map oldToNewMap; + std::map oldToNewMap; std::map, std::vector>> parentDaughterMap; if (m_mergeMCParticles) { const auto& bgParticles = backgroundEvent.get(m_MCParticleCollectionName); - int j = oparticles.size(); + int j = oparticles.size(); for (size_t i = 0; i < bgParticles.size(); ++i) { auto npart = bgParticles[i].clone(false); @@ -250,13 +251,13 @@ void OverlayTiming::mergeBackgroundFrame( continue; } - auto& calHitMap = cellIDsMap[i]; + 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; + bool add = false; for (const auto& contrib : simCaloHit.getContributions()) { if ((contrib.getTime() + timeOffset > this_start) && (contrib.getTime() + timeOffset < this_stop)) { add = true; @@ -386,12 +387,12 @@ retType OverlayTiming::operator()(const edm4hep::EventHeaderCollection& headers, // 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; + int group; + int fileIndex; size_t entry; - float timeOffset; - int bxNumber; - int physBX; + float timeOffset; + int bxNumber; + int physBX; }; std::vector reads; @@ -456,8 +457,8 @@ retType OverlayTiming::operator()(const edm4hep::EventHeaderCollection& headers, // the file index is ignored. reserve() advances the per-file cursor now // (serially); the actual ROOT read happens in phase 2, possibly on // several threads. - const int fileIndex = m_randomMix ? fileIndices[k % fileIndices.size()] : 0; - const size_t entry = m_bkgEvents->reserve(groupIndex, fileIndex); + const int fileIndex = m_randomMix ? fileIndices[k % fileIndices.size()] : 0; + const size_t entry = m_bkgEvents->reserve(groupIndex, fileIndex); reads.push_back({static_cast(groupIndex), fileIndex, entry, timeOffset, BX_number_in_train, physBX}); } } @@ -477,40 +478,38 @@ retType OverlayTiming::operator()(const edm4hep::EventHeaderCollection& headers, // serial in-order merge. struct Item { const BkgRead* read = nullptr; - podio::Frame frame; + podio::Frame frame; }; std::atomic nextRead{0}; - const size_t ntokens = static_cast(m_overlayThreads.value()); + const size_t 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, - [&](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, [&](std::shared_ptr item) { - const auto& r = *item->read; - mergeBackgroundFrame(item->frame, r.timeOffset, r.bxNumber, r.physBX, simTrackerHits, simCaloHits, - oparticles, osimTrackerHits, cellIDsMap, ocaloHitContribs); - })); + 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, + [&](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, [&](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 diff --git a/k4FWCore/components/OverlayTiming.h b/k4FWCore/components/OverlayTiming.h index 57a51de9..f1bb5645 100644 --- a/k4FWCore/components/OverlayTiming.h +++ b/k4FWCore/components/OverlayTiming.h @@ -73,9 +73,9 @@ // 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; + bool m_randomMix{false}; + bool m_allowReuse{false}; + std::string m_algName; // Sequential mode only: one persistent reader per group. std::vector m_rootFileReaders; @@ -114,10 +114,10 @@ struct EventHolder { // 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 + 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; } @@ -127,7 +127,7 @@ struct EventHolder { 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"); + 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); @@ -138,7 +138,7 @@ struct EventHolder { return reader.readEvent(rawEntry % total); } std::lock_guard lock(m_ioMutex); - const size_t total = m_totalNumberOfEvents[group][0]; + 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); } @@ -188,10 +188,10 @@ struct OverlayTiming : public k4FWCore::MultiTransformer& simTrackerHits, + const std::vector& simTrackerHits, const std::vector& simCaloHits, - edm4hep::MCParticleCollection& oparticles, - std::vector& osimTrackerHits, + edm4hep::MCParticleCollection& oparticles, + std::vector& osimTrackerHits, std::map>& cellIDsMap, std::vector& ocaloHitContribs) const; From b517240769823b7bd26ec8f4a12f0faaa4734bd8 Mon Sep 17 00:00:00 2001 From: Federico Meloni Date: Mon, 6 Jul 2026 15:39:27 +0200 Subject: [PATCH 08/15] update MT docs --- doc/OverlayTiming.md | 25 +++++++++++++++++-- .../options/TestOverlayTimingRandomMix.py | 5 +++- 2 files changed, 27 insertions(+), 3 deletions(-) diff --git a/doc/OverlayTiming.md b/doc/OverlayTiming.md index d7a9953e..0fa29f7c 100644 --- a/doc/OverlayTiming.md +++ b/doc/OverlayTiming.md @@ -36,6 +36,7 @@ It uses [`UniqueIDGenSvc`](uniqueIDGen.md) to seed the internal random number ge | `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 | @@ -111,5 +112,25 @@ overlay.RandomMixBackgroundFiles = True overlay.BackgroundFileNames = [["/path/to/bib_files/"]] ``` -All background ROOT I/O is serialized on a dedicated worker thread, so the -algorithm may be run with intra-event multithreading. +## 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/test/k4FWCoreTest/options/TestOverlayTimingRandomMix.py b/test/k4FWCoreTest/options/TestOverlayTimingRandomMix.py index 7e252138..be807602 100644 --- a/test/k4FWCoreTest/options/TestOverlayTimingRandomMix.py +++ b/test/k4FWCoreTest/options/TestOverlayTimingRandomMix.py @@ -20,7 +20,8 @@ # 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, and MCParticles are not merged (MergeMCParticles = False). +# 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 @@ -52,6 +53,8 @@ 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], From a983a902f31f1fba3940f93950d219fcb1e1bccd Mon Sep 17 00:00:00 2001 From: Federico Meloni Date: Fri, 31 Jul 2026 15:40:58 +0200 Subject: [PATCH 09/15] implement suggestion from ArinaPon --- k4FWCore/components/OverlayTiming.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/k4FWCore/components/OverlayTiming.cpp b/k4FWCore/components/OverlayTiming.cpp index 79092baf..24bf5b14 100644 --- a/k4FWCore/components/OverlayTiming.cpp +++ b/k4FWCore/components/OverlayTiming.cpp @@ -433,6 +433,7 @@ retType OverlayTiming::operator()(const edm4hep::EventHeaderCollection& headers, } // Overlay the background events to each bunchcrossing in the bunch train + size_t fileCursor = 0; for (int bxInTrain = 0; bxInTrain < m_NBunchTrain; ++bxInTrain) { const int BX_number_in_train = permutation.at(bxInTrain); @@ -457,7 +458,7 @@ retType OverlayTiming::operator()(const edm4hep::EventHeaderCollection& headers, // the file index is ignored. reserve() advances the per-file cursor now // (serially); the actual ROOT read happens in phase 2, possibly on // several threads. - const int fileIndex = m_randomMix ? fileIndices[k % fileIndices.size()] : 0; + const int fileIndex = m_randomMix ? fileIndices[fileCursor++ % fileIndices.size()] : 0; const size_t entry = m_bkgEvents->reserve(groupIndex, fileIndex); reads.push_back({static_cast(groupIndex), fileIndex, entry, timeOffset, BX_number_in_train, physBX}); } From 9aab53b3cf60f3e8522b28e83a4cccd38dc7c38d Mon Sep 17 00:00:00 2001 From: Federico Meloni Date: Fri, 31 Jul 2026 22:56:10 +0200 Subject: [PATCH 10/15] Fix clang-format in OverlayTiming.cpp --- k4FWCore/components/OverlayTiming.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/k4FWCore/components/OverlayTiming.cpp b/k4FWCore/components/OverlayTiming.cpp index 24bf5b14..e7fa1b74 100644 --- a/k4FWCore/components/OverlayTiming.cpp +++ b/k4FWCore/components/OverlayTiming.cpp @@ -458,7 +458,7 @@ retType OverlayTiming::operator()(const edm4hep::EventHeaderCollection& headers, // the file index is ignored. reserve() advances the per-file cursor now // (serially); the actual ROOT read happens in phase 2, possibly on // several threads. - const int fileIndex = m_randomMix ? fileIndices[fileCursor++ % fileIndices.size()] : 0; + const int fileIndex = m_randomMix ? fileIndices[fileCursor++ % fileIndices.size()] : 0; const size_t entry = m_bkgEvents->reserve(groupIndex, fileIndex); reads.push_back({static_cast(groupIndex), fileIndex, entry, timeOffset, BX_number_in_train, physBX}); } From 349e5894e0bd3a5780062878a88488c03ac2f297 Mon Sep 17 00:00:00 2001 From: Federico Meloni Date: Fri, 31 Jul 2026 23:19:42 +0200 Subject: [PATCH 11/15] updated comments and added shuffle on wrap --- k4FWCore/components/OverlayTiming.cpp | 24 ++++++++++++++++++------ 1 file changed, 18 insertions(+), 6 deletions(-) diff --git a/k4FWCore/components/OverlayTiming.cpp b/k4FWCore/components/OverlayTiming.cpp index e7fa1b74..963f8e8c 100644 --- a/k4FWCore/components/OverlayTiming.cpp +++ b/k4FWCore/components/OverlayTiming.cpp @@ -432,7 +432,11 @@ retType OverlayTiming::operator()(const edm4hep::EventHeaderCollection& headers, } } - // 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); @@ -454,11 +458,19 @@ retType OverlayTiming::operator()(const edm4hep::EventHeaderCollection& headers, } const float timeOffset = BX_number_in_train * m_deltaT; for (int k = 0; k < NOverlay_to_this_BX; ++k) { - // In random-mix mode pick a random file of the group; in sequential mode - // the file index is ignored. reserve() advances the per-file cursor now - // (serially); the actual ROOT read happens in phase 2, possibly on - // several threads. - const int fileIndex = m_randomMix ? fileIndices[fileCursor++ % fileIndices.size()] : 0; + // 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}); } From b8b4f57cf37504b16a6c6f0d54729d9c7cd6d995 Mon Sep 17 00:00:00 2001 From: Federico Meloni Date: Fri, 31 Jul 2026 23:26:24 +0200 Subject: [PATCH 12/15] make pre-commit happy again --- k4FWCore/components/OverlayTiming.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/k4FWCore/components/OverlayTiming.cpp b/k4FWCore/components/OverlayTiming.cpp index 963f8e8c..78a4a506 100644 --- a/k4FWCore/components/OverlayTiming.cpp +++ b/k4FWCore/components/OverlayTiming.cpp @@ -459,7 +459,7 @@ retType OverlayTiming::operator()(const edm4hep::EventHeaderCollection& headers, const float timeOffset = BX_number_in_train * m_deltaT; for (int k = 0; k < NOverlay_to_this_BX; ++k) { // In random-mix mode walk the shuffled permutation so that consecutive - // overlaid events draw a distinct set of files. + // 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. From 6e6b20ceec71dd2c1d3084876f7866a5f33c3625 Mon Sep 17 00:00:00 2001 From: Federico Meloni Date: Mon, 10 Aug 2026 13:53:43 +0200 Subject: [PATCH 13/15] preempt clang tidy --- k4FWCore/components/OverlayTiming.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/k4FWCore/components/OverlayTiming.cpp b/k4FWCore/components/OverlayTiming.cpp index 07ddd903..a9975267 100644 --- a/k4FWCore/components/OverlayTiming.cpp +++ b/k4FWCore/components/OverlayTiming.cpp @@ -460,6 +460,7 @@ retType OverlayTiming::operator()(const edm4hep::EventHeaderCollection& headers, 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) { // In random-mix mode walk the shuffled permutation so that consecutive // overlaid events draw a distinct set of files. @@ -510,7 +511,7 @@ retType OverlayTiming::operator()(const edm4hep::EventHeaderCollection& headers, }) & tbb::make_filter, std::shared_ptr>( tbb::filter_mode::parallel, - [&](std::shared_ptr item) -> std::shared_ptr { + [&](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 @@ -521,7 +522,7 @@ retType OverlayTiming::operator()(const edm4hep::EventHeaderCollection& headers, return item; }) & tbb::make_filter, void>( - tbb::filter_mode::serial_in_order, [&](std::shared_ptr item) { + 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); From 3b80ce75d171f3817da523c5494f08340ffb7029 Mon Sep 17 00:00:00 2001 From: Federico Meloni Date: Mon, 10 Aug 2026 14:09:39 +0200 Subject: [PATCH 14/15] fix clang-tidy --- k4FWCore/components/OverlayTiming.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/k4FWCore/components/OverlayTiming.cpp b/k4FWCore/components/OverlayTiming.cpp index a9975267..349b6e59 100644 --- a/k4FWCore/components/OverlayTiming.cpp +++ b/k4FWCore/components/OverlayTiming.cpp @@ -498,7 +498,7 @@ retType OverlayTiming::operator()(const edm4hep::EventHeaderCollection& headers, podio::Frame frame; }; std::atomic nextRead{0}; - const size_t ntokens = static_cast(m_overlayThreads.value()); + 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 { From caa0afb1bf25a3486139bf6ede81cedf6efc8832 Mon Sep 17 00:00:00 2001 From: Federico Meloni Date: Tue, 1 Sep 2026 12:48:48 +0200 Subject: [PATCH 15/15] add isolate --- doc/OverlayTiming.md | 27 ++++++++++- k4FWCore/components/OverlayTiming.cpp | 66 +++++++++++++++------------ 2 files changed, 63 insertions(+), 30 deletions(-) diff --git a/doc/OverlayTiming.md b/doc/OverlayTiming.md index 0fa29f7c..62518248 100644 --- a/doc/OverlayTiming.md +++ b/doc/OverlayTiming.md @@ -128,9 +128,32 @@ 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. +algorithm also remains safe to run under Gaudi's intra-event multithreading. 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. + +### Interplay with the Gaudi scheduler + +`OverlayThreads` and the scheduler's `ThreadPoolSize` are not independent: they +draw from the same pool of threads. `ThreadPoolSvc` sets a process-wide TBB +`global_control` of `ThreadPoolSize + maxParallelismExtra + 1` and creates the +task arena into which `AvalancheSchedulerSvc` enqueues every algorithm, and the +parallel reading runs inside that same arena. This means the machine is never +oversubscribed, but it also means: + +- Gaudi has no way of knowing about the extra parallelism. The scheduler counts + algorithms in flight, not threads, so it keeps dispatching other algorithms + while `OverlayTiming` is fanned out. +- `OverlayThreads` bounds the number of background reads in flight, not the + number of threads actually available. With `ThreadPoolSize = 1` the arena has + two threads, and a large `OverlayThreads` will not buy more than that. + +There is currently no way for a functional algorithm to declare its internal +parallelism to the scheduler (`Asynchronous` is the Boost.Fiber path for +offloaded work, not this). In practice `OverlayThreads > 1` pays off when +`ThreadPoolSize` is small -- branching out inside a single event rather than +running many events concurrently -- and raising both just repartitions the same +threads. If the extra threads should come on top of the scheduler's pool, +`AvalancheSchedulerSvc.maxParallelismExtra` raises the TBB limit accordingly. diff --git a/k4FWCore/components/OverlayTiming.cpp b/k4FWCore/components/OverlayTiming.cpp index 349b6e59..3b78583a 100644 --- a/k4FWCore/components/OverlayTiming.cpp +++ b/k4FWCore/components/OverlayTiming.cpp @@ -34,6 +34,7 @@ #include #include +#include #include #include @@ -499,34 +500,43 @@ retType OverlayTiming::operator()(const edm4hep::EventHeaderCollection& headers, }; 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); - })); + // Run the pipeline in an isolated region. This is defensive: while this thread + // waits here TBB is otherwise free to steal another algorithm's task from the + // scheduler's arena, and Gaudi's AlgTask overwrites the thread-local + // EventContext and whiteboard partition without restoring them, so in principle + // we could resume on the wrong event slot. It has not been observed in practice, + // but isolating costs nothing. + tbb::this_task_arena::isolate([&] { + 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