From db966bf6d0764842de77703b05e93ab83696c32c Mon Sep 17 00:00:00 2001 From: jmcarcell Date: Sat, 22 Jun 2024 10:42:52 +0200 Subject: [PATCH] Update with some comments and name changes --- k4Reco/Overlay/components/OverlayTiming.cpp | 44 ++++++++------- k4Reco/Overlay/components/OverlayTiming.h | 62 +++++++++++++-------- 2 files changed, 64 insertions(+), 42 deletions(-) diff --git a/k4Reco/Overlay/components/OverlayTiming.cpp b/k4Reco/Overlay/components/OverlayTiming.cpp index 8db6e0a..24c8bac 100644 --- a/k4Reco/Overlay/components/OverlayTiming.cpp +++ b/k4Reco/Overlay/components/OverlayTiming.cpp @@ -55,18 +55,25 @@ StatusCode OverlayTiming::initialize() { } } + 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); + return StatusCode::FAILURE; + } + if (m_Noverlay.empty()) { info() << "Using the default number of overlay events (1) for each group, since none was specified with " "NumberBackground " << endmsg; - m_Noverlay.value() = std::vector(m_bkgEvents->size(), 1); + m_Noverlay = std::vector(m_bkgEvents->size(), 1); } if (m_Poisson.empty()) { info() << "Using the default overlay mode (no Poission distribution) for each group, since none was specified with " "Poisson_random_Noverlay" << endmsg; - m_Poisson.value() = std::vector(m_bkgEvents->size(), false); + m_Poisson = std::vector(m_bkgEvents->size(), false); } return StatusCode::SUCCESS; @@ -87,7 +94,6 @@ retType OverlayTiming::operator()(const edm4hep::EventHeaderCollection& for (size_t i = 0; i < simCaloHits.size(); ++i) { ocaloHitContribs.emplace_back(edm4hep::CaloHitContributionCollection()); } - // Copy MCParticles for physics event into a new collection for (auto&& part : particles) { @@ -154,17 +160,17 @@ retType OverlayTiming::operator()(const edm4hep::EventHeaderCollection& // Iterate over each group of files and parameters for (size_t groupIndex = 0; groupIndex < m_bkgEvents->size(); groupIndex++) { - if (_randomBX) { - _BX_phys = std::uniform_int_distribution(0, _nBunchTrain - 1)(m_engine); - debug() << "Physics Event was placed in the " << _BX_phys << " bunch crossing!" << endmsg; + if (m_randomBX) { + m_physBX = std::uniform_int_distribution(0, _nBunchTrain - 1)(m_engine); + debug() << "Physics Event was placed in the " << m_physBX << " bunch crossing!" << endmsg; } // define a permutation for the events to overlay -- the physics event is per definition at position 0 std::vector permutation; // Permutation has negative values and the last one is 0 - // if (!_randomBX) then _BX_phys (default = 1) - for (int i = -(_BX_phys - 1); i < _nBunchTrain - (_BX_phys - 1); ++i) { + // if (!m_randomBX) then m_physBX (default = 1) + for (int i = -(m_physBX - 1); i < _nBunchTrain - (m_physBX - 1); ++i) { permutation.push_back(i); } std::shuffle(permutation.begin(), permutation.end(), m_engine); @@ -190,19 +196,19 @@ retType OverlayTiming::operator()(const edm4hep::EventHeaderCollection& NOverlay_to_this_BX = m_Noverlay[groupIndex]; } - debug() << "Will overlay " << NOverlay_to_this_BX << " events to BX number " << BX_number_in_train + _BX_phys + debug() << "Will overlay " << NOverlay_to_this_BX << " events to BX number " << BX_number_in_train + m_physBX << endmsg; for (int k = 0; k < NOverlay_to_this_BX; ++k) { - info() << "Overlaying event " << m_bkgEvents->m_nextEntry << " to BX number " << BX_number_in_train + _BX_phys + info() << "Overlaying event " << m_bkgEvents->m_nextEntry << " to BX number " << BX_number_in_train + m_physBX << endmsg; - auto backgroundEvent = - m_bkgEvents->m_rootFileReaders[groupIndex].readEvent(m_bkgEvents->m_nextEntry); - // m_bkgEvents->m_nextEntry++; - // if(m_bkgEvents->m_nextEntry >= m_bkgEvents->m_totalNumberOfEvents[groupIndex] && !m_allowReusingBackgroundFiles) { - // throw GaudiException("No more events in background file", name(), StatusCode::FAILURE); - // } - // m_bkgEvents->m_nextEntry %= m_bkgEvents->m_totalNumberOfEvents[groupIndex]; + auto backgroundEvent = m_bkgEvents->m_rootFileReaders[groupIndex].readEvent(m_bkgEvents->m_nextEntry); + m_bkgEvents->m_nextEntry++; + if (m_bkgEvents->m_nextEntry >= m_bkgEvents->m_totalNumberOfEvents[groupIndex] && + !m_allowReusingBackgroundFiles) { + throw GaudiException("No more events in background file", name(), StatusCode::FAILURE); + } + m_bkgEvents->m_nextEntry %= m_bkgEvents->m_totalNumberOfEvents[groupIndex]; auto availableCollections = backgroundEvent.getAvailableCollections(); // Either 0 or negative @@ -262,7 +268,7 @@ retType OverlayTiming::operator()(const edm4hep::EventHeaderCollection& } 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 - _BX_phys) * _T_diff) { + if (this_stop <= (BX_number_in_train - m_physBX) * _T_diff) { info() << "Skipping collection " << name << " as it is not in the integration window" << endmsg; continue; } @@ -291,7 +297,7 @@ retType OverlayTiming::operator()(const edm4hep::EventHeaderCollection& } 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 - _BX_phys) * _T_diff) { + if (this_stop <= (BX_number_in_train - m_physBX) * _T_diff) { info() << "Skipping collection " << name << " as it is not in the integration window" << endmsg; continue; } diff --git a/k4Reco/Overlay/components/OverlayTiming.h b/k4Reco/Overlay/components/OverlayTiming.h index e1bc503..9a7342e 100644 --- a/k4Reco/Overlay/components/OverlayTiming.h +++ b/k4Reco/Overlay/components/OverlayTiming.h @@ -1,3 +1,18 @@ +/** Background overlay algorithm + + This algorithm overlays background events on top of the signal events. The + background events are read from a set of input files, and the signal events + are the input in the main event loop. + + The MCParticleCollection in signal are background are overlaid into one + collection. The SimTrackerHit collections are cropped and overlayed if they + are in the time window. The SimCalorimeterHit collections are overlayed + based on the cellID. If a signal hit has the same cellID as a background + hit, they are combined into a single hit. Only hits that have + CaloHitContributions in the time range are considered. + +**/ + #include "podio/Frame.h" #include "podio/ROOTReader.h" #include "podio/Reader.h" @@ -18,10 +33,10 @@ #include struct EventHolder { - std::vector> m_fileNames; - std::vector m_rootFileReaders; - std::vector m_totalNumberOfEvents; - std::map m_events; + std::vector> m_fileNames; + std::vector m_rootFileReaders; + std::vector m_totalNumberOfEvents; + std::map m_events; int m_nextEntry = 0; @@ -33,26 +48,28 @@ struct EventHolder { } EventHolder() = default; + // TODO: Cache functionality + // podio::Frame& read + size_t size() const { return m_fileNames.size(); } }; -using retType = std::tuple, - std::vector, - std::vector>; +using retType = + std::tuple, + std::vector, std::vector>; struct OverlayTiming : public k4FWCore::MultiTransformer&, const std::vector& // const std::map& - )> { + )> { OverlayTiming(const std::string& name, ISvcLocator* svcLoc) : MultiTransformer( name, svcLoc, {KeyValues("EventHeader", {"EventHeader"}), KeyValues("MCParticles", {"DefaultMCParticles"}), KeyValues("SimTrackerHits", {"DefaultSimTrackerHits"}), - KeyValues("SimCalorimeterHits", {"DefaultSimCalorimeterHits"}) - }, + KeyValues("SimCalorimeterHits", {"DefaultSimCalorimeterHits"})}, {KeyValues("OutputMCParticles", {"NewMCParticles"}), KeyValues("OutputSimTrackerHits", {"MCParticles1"}), KeyValues("OutputSimCalorimeterHits", {"MCParticles2"}), KeyValues("OutputCaloHitContributions", {"OverlayCaloHitContributions"})}) {} @@ -63,27 +80,23 @@ struct OverlayTiming : public k4FWCore::MultiTransformer& simTrackerHits, - const std::vector& simCalorimeterHits - ) const final; + const std::vector& simTrackerHits, + const std::vector& simCalorimeterHits) const final; std::pair define_time_windows(const std::string& Collection_name) const; private: - // These correspond to the index position in the argument list constexpr static int TRACKERHIT_INDEX_POSITION = 2; constexpr static int SIMCALOHIT_INDEX_POSITION = 3; - Gaudi::Property _randomBX{this, "RandomBx", false, + Gaudi::Property m_randomBX{this, "RandomBx", false, "Place the physics event at an random position in the train: overrides PhysicsBX"}; - mutable Gaudi::Property _BX_phys{this, "PhysicsBX", 1, "Number of the Bunch crossing of the physics event"}; - Gaudi::Property _tpcVdrift_mm_ns{this, "TPCDriftvelocity", float(5.0e-2), - "Drift velocity in TPC [mm/ns] - default 5.0e-2 (5cm/us)"}; + mutable Gaudi::Property m_physBX{this, "PhysicsBX", 1, "Number of the Bunch crossing of the physics event"}; Gaudi::Property _nBunchTrain{this, "NBunchtrain", 1, "Number of bunches in a bunch train"}; - Gaudi::Property m_startWithBackgroundFile{this, "StartBackgroundFileIndex", -1, - "Which background file to startWith"}; - Gaudi::Property m_startWithBackgroundEvent{this, "StartBackgroundEventIndex", -1, + // Gaudi::Property m_startWithBackgroundFile{this, "StartBackgroundFileIndex", -1, + // "Which background file to startWith"}; + Gaudi::Property m_startWithBackgroundEvent{this, "StartBackgroundEventIndex", -1, "Which background event to startWith"}; Gaudi::Property>> m_inputFileNames{ @@ -116,7 +129,10 @@ struct OverlayTiming : public k4FWCore::MultiTransformer m_allowReusingBackgroundFiles{ this, "AllowReusingBackgroundFiles", false, "If true the same background file can be used for the same event"}; + // Gaudi::Property m_maxCachedFrames{ + // this, "MaxCachedFrames", 0, "Maximum number of frames cached from background files"}; + private: - mutable std::mt19937 m_engine; - SmartIF m_uidSvc; + inline static thread_local std::mt19937 m_engine; + SmartIF m_uidSvc; };