From a3baba052d43409a1e59796e16fa6d3a2b93c965 Mon Sep 17 00:00:00 2001 From: Wouter Deconinck Date: Mon, 11 Jul 2022 15:58:32 +0200 Subject: [PATCH 1/3] First integration of HEP-FCC/FWCore PileUp tools --- JugBase/JugBase/IEDMMergeTool.h | 23 +++ JugBase/JugBase/IPileUpTool.h | 36 ++++ JugBase/src/components/ConstPileUp.cpp | 24 +++ JugBase/src/components/ConstPileUp.h | 33 ++++ .../PileupDigiTrackHitMergeTool.cpp | 162 ++++++++++++++++++ .../components/PileupDigiTrackHitMergeTool.h | 99 +++++++++++ JugBase/src/components/PileupHitMergeTool.cpp | 119 +++++++++++++ JugBase/src/components/PileupHitMergeTool.h | 71 ++++++++ JugBase/src/components/PileupOverlayAlg.cpp | 97 +++++++++++ JugBase/src/components/PileupOverlayAlg.h | 73 ++++++++ .../components/PileupParticlesMergeTool.cpp | 106 ++++++++++++ .../src/components/PileupParticlesMergeTool.h | 51 ++++++ JugBase/src/components/PoissonPileUp.cpp | 46 +++++ JugBase/src/components/PoissonPileUp.h | 37 ++++ JugBase/src/components/RangePileUp.cpp | 28 +++ JugBase/src/components/RangePileUp.h | 33 ++++ 16 files changed, 1038 insertions(+) create mode 100644 JugBase/JugBase/IEDMMergeTool.h create mode 100644 JugBase/JugBase/IPileUpTool.h create mode 100644 JugBase/src/components/ConstPileUp.cpp create mode 100644 JugBase/src/components/ConstPileUp.h create mode 100644 JugBase/src/components/PileupDigiTrackHitMergeTool.cpp create mode 100644 JugBase/src/components/PileupDigiTrackHitMergeTool.h create mode 100644 JugBase/src/components/PileupHitMergeTool.cpp create mode 100644 JugBase/src/components/PileupHitMergeTool.h create mode 100644 JugBase/src/components/PileupOverlayAlg.cpp create mode 100644 JugBase/src/components/PileupOverlayAlg.h create mode 100644 JugBase/src/components/PileupParticlesMergeTool.cpp create mode 100644 JugBase/src/components/PileupParticlesMergeTool.h create mode 100644 JugBase/src/components/PoissonPileUp.cpp create mode 100644 JugBase/src/components/PoissonPileUp.h create mode 100644 JugBase/src/components/RangePileUp.cpp create mode 100644 JugBase/src/components/RangePileUp.h diff --git a/JugBase/JugBase/IEDMMergeTool.h b/JugBase/JugBase/IEDMMergeTool.h new file mode 100644 index 00000000..f529db98 --- /dev/null +++ b/JugBase/JugBase/IEDMMergeTool.h @@ -0,0 +1,23 @@ +#ifndef IEDMMergeTool_h +#define IEDMMergeTool_h + +#include "GaudiKernel/IAlgTool.h" + +/** @class IEDMMergeTool + * + * Interface for the tool used in the overlay algorithm. + * Must implement the correct collection and I/O and merging behavior, + * especially for the case when there are associations between parts of the EDM. + */ +class IEDMMergeTool : virtual public IAlgTool { +public: + /// read pileup collection from store and handle them internally + virtual StatusCode readPileupCollection(podio::EventStore& store) = 0; + /// merge internally stored pileup (and signal) collections and make them accessible (typically via event store) + virtual StatusCode mergeCollections() = 0; + + /// read signal collection from event store and + virtual StatusCode readSignal() = 0; +}; + +#endif // IEDMMergeTool_h diff --git a/JugBase/JugBase/IPileUpTool.h b/JugBase/JugBase/IPileUpTool.h new file mode 100644 index 00000000..34dbf9a1 --- /dev/null +++ b/JugBase/JugBase/IPileUpTool.h @@ -0,0 +1,36 @@ +#ifndef IPileUpTool_h +#define IPileUpTool_h + +#include "GaudiKernel/IAlgTool.h" + +/** @class IPileUpTool IPileUpTool.h "Generation/IPileUpTool.h" + * + * Abstract interface to pile up tools. Generates the number of pile-up + * interactions to generate for each event. + * + * @author Patrick Robbe + * @date 2005-08-17 + */ + +namespace Jug::Base { + +static const InterfaceID IID_IPileUpTool("IPileUpTool", 3, 0); + +class IPileUpTool : virtual public IAlgTool { +public: + /** Computes the number of pile-up interactions in the event. + * @param[out] currentLuminosity Luminosity of the current event. + * @return Number of pile-up interactions to generate. + */ + virtual unsigned int numberOfPileUp() = 0; + + virtual double getMeanPileUp() = 0; + + /// Print various counters at the end of the job + virtual void printPileUpCounters() = 0; + +}; + +} + +#endif // IPileUpTool_h diff --git a/JugBase/src/components/ConstPileUp.cpp b/JugBase/src/components/ConstPileUp.cpp new file mode 100644 index 00000000..b5140080 --- /dev/null +++ b/JugBase/src/components/ConstPileUp.cpp @@ -0,0 +1,24 @@ +#include "ConstPileUp.h" + +DECLARE_COMPONENT(ConstPileUp) + +ConstPileUp::ConstPileUp(const std::string& type, const std::string& name, const IInterface* parent) + : GaudiTool(type, name, parent) { + declareInterface(this); +} + +ConstPileUp::~ConstPileUp() { ; } + +StatusCode ConstPileUp::initialize() { + StatusCode sc = GaudiTool::initialize(); + printPileUpCounters(); + return sc; +} + +unsigned int ConstPileUp::numberOfPileUp() { return m_numPileUpEvents; } + +double ConstPileUp::getMeanPileUp() { return m_numPileUpEvents; } + +void ConstPileUp::printPileUpCounters() { + info() << "Current number of pileup events: " << m_numPileUpEvents << endmsg; +} diff --git a/JugBase/src/components/ConstPileUp.h b/JugBase/src/components/ConstPileUp.h new file mode 100644 index 00000000..c64a301f --- /dev/null +++ b/JugBase/src/components/ConstPileUp.h @@ -0,0 +1,33 @@ +#ifndef GENERATION_CONSTPILEUP_H +#define GENERATION_CONSTPILEUP_H + +#include "JugBase/IPileUpTool.h" +#include "GaudiAlg/GaudiTool.h" + +/** @class ConstPileUp + * + * Tool to generate number of pile-up events to be mixed with signal event. + * Concrete implementation of a IPileUpTool -- the most trivial one, actually, + * returning just a constant (that can be specified as a property). + * The interface is kept close to the originals in LHCb's Gauss(ino). + * + * @author Valentin Volkl + * @date 2015-12-16 + */ +class ConstPileUp : public GaudiTool, virtual public Jug::Base::IPileUpTool { +public: + ConstPileUp(const std::string& type, const std::string& name, const IInterface* parent); + virtual ~ConstPileUp(); + virtual StatusCode initialize(); + /// in the special case of constant pileup, + /// this method is identical to getMeanPileUp + virtual unsigned int numberOfPileUp(); + virtual double getMeanPileUp(); + virtual void printPileUpCounters(); + +private: + /// Number of min bias events to pile on signal event. + Gaudi::Property m_numPileUpEvents{this, "numPileUpEvents", 0, "number of pile-up events"}; +}; + +#endif // GENERATION_CONSTPILEUP_H diff --git a/JugBase/src/components/PileupDigiTrackHitMergeTool.cpp b/JugBase/src/components/PileupDigiTrackHitMergeTool.cpp new file mode 100644 index 00000000..4c88b226 --- /dev/null +++ b/JugBase/src/components/PileupDigiTrackHitMergeTool.cpp @@ -0,0 +1,162 @@ +#include "PileupDigiTrackHitMergeTool.h" + +#include "podio/EventStore.h" + +#include "edm4hep/DigiTrackHitAssociationCollection.h" +#include "edm4hep/MCParticle.h" +#include "edm4hep/MCParticleCollection.h" +#include "edm4hep/TrackHitCollection.h" + +DECLARE_COMPONENT(PileupDigiTrackHitMergeTool) +DECLARE_COMPONENT_WITH_ID(PileupDigiTrackHitMergeTool, "PileupDigiTrackHitMergeTool") + +PileupDigiTrackHitMergeTool::PileupDigiTrackHitMergeTool(const std::string& aType, const std::string& aName, + const IInterface* aParent) + : GaudiTool(aType, aName, aParent) { + declareProperty("signalPositionedHits", m_hitsSignal); + declareProperty("signalDigiHits", m_posHitsSignal); + declareProperty("signalParticles", m_particlesSignal); + declareProperty("mergedPositionedHits", m_hitsMerged); + declareProperty("mergedDigiHits", m_digiHitsMerged); + declareProperty("mergedParticles", m_particlesMerged); +} + +StatusCode PileupDigiTrackHitMergeTool::initialize() { return StatusCode::SUCCESS; } + +StatusCode PileupDigiTrackHitMergeTool::finalize() { return StatusCode::SUCCESS; } + +StatusCode PileupDigiTrackHitMergeTool::readPileupCollection(podio::EventStore& store) { + // local pointers, to be filled by the event store + const edm4hep::PositionedTrackHitCollection* hitCollection; + const edm4hep::DigiTrackHitAssociationCollection* digiHitCollection; + const edm4hep::MCParticleCollection* particleCollection; + + // get collection address and store it in container + bool hitCollectionPresent = store.get(m_pileupHitsBranchName, hitCollection); + if (hitCollectionPresent) { + m_hitCollections.push_back(hitCollection); + } else { + warning() << "No collection could be read from branch " << m_pileupHitsBranchName << endmsg; + return StatusCode::FAILURE; + } + + /// as above, for the positioned collection + bool digiHitCollectionPresent = store.get(m_pileupPosHitsBranchName, digiHitCollection); + if (digiHitCollectionPresent) { + m_digiTrackHitsCollections.push_back(digiHitCollection); + } else { + warning() << "No collection could be read from branch " << m_pileupPosHitsBranchName << endmsg; + return StatusCode::FAILURE; + } + + /// as above, for the particle collection (if wanted) + if (m_mergeParticles) { + bool particleCollectionPresent = store.get(m_pileupParticleBranchName, particleCollection); + if (particleCollectionPresent) { + m_particleCollections.push_back(particleCollection); + } else { + warning() << "No collection could be read from branch " << m_pileupParticleBranchName << endmsg; + return StatusCode::FAILURE; + } + } + + return StatusCode::SUCCESS; +} + +StatusCode PileupDigiTrackHitMergeTool::readSignal() { + // get collection from event store + auto collHitsSig = m_hitsSignal.get(); + auto collPosHitsSig = m_posHitsSignal.get(); + auto collPartSig = m_particlesSignal.get(); + // store them in internal container + m_hitCollections.push_back(collHitsSig); + m_digiTrackHitsCollections.push_back(collPosHitsSig); + // check if particles should be merged + if (m_mergeParticles) { + m_particleCollections.push_back(collPartSig); + } + + return StatusCode::SUCCESS; +} + +StatusCode PileupDigiTrackHitMergeTool::mergeCollections() { + + // ownership given to data service at end of execute + edm4hep::PositionedTrackHitCollection* collHitsMerged = new edm4hep::PositionedTrackHitCollection(); + edm4hep::DigiTrackHitAssociationCollection* collDigiHitsMerged = new edm4hep::DigiTrackHitAssociationCollection(); + edm4hep::MCParticleCollection* collParticlesMerged = new edm4hep::MCParticleCollection(); + + unsigned int collectionCounter = 0; + for (size_t j = 0; j < m_digiTrackHitsCollections.size(); j++) { + unsigned int trackIDCounter = 0; + std::unordered_map oldToNewTrackIDs; + + auto offset = collectionCounter * m_trackIDCollectionOffset; + auto digiHitColl = m_digiTrackHitsCollections.at(j); + auto hitColl = m_hitCollections.at(j); + for (int i = 0; i < digiHitColl->size(); i++) { + auto clon = hitColl->at(i).clone(); + // add pileup vertex counter with an offset + // i.e. for the signal event, 'bits' is just the trackID taken from geant + // for the n-th pileup event, 'bits' is the trackID + n * offset + // offset needs to be big enough to ensure uniqueness of trackID + if (isTrackerHit(hitColl->at(i).cellId())) { + if (trackIDCounter > m_trackIDCollectionOffset) { + error() << "Event contains too many tracks to guarantee a unique trackID"; + error() << " The offset width or trackID field size needs to be adjusted!" << endmsg; + return StatusCode::FAILURE; + } + auto search = oldToNewTrackIDs.find(clon.bits()); + if (search != oldToNewTrackIDs.end()) { + // was already present - use already created new trackID + unsigned newTrackID = (search->second + offset); + clon.bits(newTrackID); + } else { + // was not already present - create new trackID + oldToNewTrackIDs[clon.bits()] = trackIDCounter; + + unsigned newTrackID = (trackIDCounter + offset); + clon.bits(newTrackID); + + trackIDCounter++; + } + } + auto newDigiHit = digiHitColl->at(i).clone(); + newDigiHit.hit(clon); + collHitsMerged->push_back(clon); + collDigiHitsMerged->push_back(newDigiHit); + } + ///@todo find possibility in case not only tracker hits are produced + if (m_mergeParticles) { + auto partColl = m_particleCollections.at(j); + for (int i = 0; i < partColl->size(); i++) { + auto clonParticle = partColl->at(i).clone(); + auto search = oldToNewTrackIDs.find(clonParticle.bits()); + if (search != oldToNewTrackIDs.end()) { + // was not already present - create new trackID + clonParticle.bits(search->second + offset); + } else { + // was not already present - create new trackID + // this can happen if particles do not reach the tracker + oldToNewTrackIDs[clonParticle.bits()] = trackIDCounter; + + unsigned newTrackID = (trackIDCounter + offset); + clonParticle.bits(newTrackID); + + trackIDCounter++; + } + collParticlesMerged->push_back(clonParticle); + } + } + ++collectionCounter; + } + m_hitsMerged.put(collHitsMerged); + m_digiHitsMerged.put(collDigiHitsMerged); + m_particlesMerged.put(collParticlesMerged); + + m_hitCollections.clear(); + m_digiTrackHitsCollections.clear(); + m_particleCollections.clear(); + + return StatusCode::SUCCESS; +} diff --git a/JugBase/src/components/PileupDigiTrackHitMergeTool.h b/JugBase/src/components/PileupDigiTrackHitMergeTool.h new file mode 100644 index 00000000..f59612ff --- /dev/null +++ b/JugBase/src/components/PileupDigiTrackHitMergeTool.h @@ -0,0 +1,99 @@ +#ifndef FWCORE_PILEUPDIGITRACKERHITSMERGETOOL_H +#define FWCORE_PILEUPDIGITRACKERHITSMERGETOOL_H + +// Gaudi +#include "GaudiAlg/GaudiTool.h" + +// FCCSW +#include "JugBase/DataHandle.h" +#include "JugBase/IEDMMergeTool.h" + +namespace fcc { +class PositionedTrackHit; +class PositionedTrackHitCollection; +class DigiTrackHitAssociation; +class DigiTrackHitAssociationCollection; +class MCParticle; +class MCParticleCollection; +} + +/** @class PileupDigiTrackHitMergeTool + * + * Implemenation of the MergeTool for Digitization hits and, if requested, simulated hits. + * While merging, this algorithm tries to keep the trackIDs unique by adding the pileup event number with an offset. + * This should be transparent, but the trackIDs will be non-consecutive. + * + */ + +class PileupDigiTrackHitMergeTool : public GaudiTool, virtual public IEDMMergeTool { +public: + explicit PileupDigiTrackHitMergeTool(const std::string& aType, const std::string& aName, const IInterface* aParent); + virtual ~PileupDigiTrackHitMergeTool() = default; + + virtual StatusCode initialize() override final; + + virtual StatusCode finalize() override final; + + /// fill the member containers of collection pointers in this class with a collection from the event store + virtual StatusCode readPileupCollection(podio::EventStore& store) override final; + + /// merge the collections in the member container and put them on the event store + virtual StatusCode mergeCollections() override final; + + /// fill the member container of collection pointer with a collection from the event store + virtual StatusCode readSignal() override final; + +private: + /// Name of the branch from which to read pileup collection + Gaudi::Property m_pileupHitsBranchName{this, "pileupHitsBranch", "positionedHits", + "Name of the branch from which to read pileup collection"}; + /// Name of the branch from which to read pileup collection + Gaudi::Property m_pileupPosHitsBranchName{this, "pileupDigiHitsBranch", "digiHits", + "Name of the branch from which to read pileup collection"}; + + /// Name of the branch from which to read pileup collection + Gaudi::Property m_pileupParticleBranchName{this, "pileupParticleBranch", "simParticles", + "Name of the branch from which to read pileup collection"}; + + /// Flag indicating if particles are present and should be merged + Gaudi::Property m_mergeParticles{this, "mergeParticles", false, + "Flag indicating if particles are present and should be merged"}; + + /// internal container to keep of collections to be merged + std::vector m_hitCollections; + /// internal container to keep of collections to be merged + std::vector m_digiTrackHitsCollections; + /// internal container to keep of collections to be merge + std::vector m_particleCollections; + + /// Output of this tool: merged collection + DataHandle m_hitsMerged{"overlay/hits", Gaudi::DataHandle::Writer, this}; + /// Output of this tool: merged collection + DataHandle m_digiHitsMerged{"overlay/digiHits", Gaudi::DataHandle::Writer, + this}; + /// Output of this tool: merged collection + DataHandle m_particlesMerged{"overlay/particles", Gaudi::DataHandle::Writer, this}; + + /// input to this tool: signal collection + DataHandle m_hitsSignal{"overlay/signalHits", Gaudi::DataHandle::Reader, this}; + /// input to this tool: signal collection + DataHandle m_posHitsSignal{"overlay/signalDigiHits", + Gaudi::DataHandle::Reader, this}; + /// input to this tool: signal collection + DataHandle m_particlesSignal{"overlay/signalParticles", Gaudi::DataHandle::Reader, this}; + + /// offset with which the pileup event number is added to the trackID + const unsigned int m_trackIDCollectionOffset = 2.5e6; + + /// Mask to obtain system ID + const unsigned long long m_systemMask = 0xf; + + /// Private function internally used to check if a hit is a tracker hit + bool isTrackerHit(unsigned long long cellId) const; +}; + +inline bool PileupDigiTrackHitMergeTool::isTrackerHit(unsigned long long cellId) const { + return ((cellId & m_systemMask) < 5); +} + +#endif // FWCORE_PILEUPDIGITRACKERHITSMERGETOOL_H diff --git a/JugBase/src/components/PileupHitMergeTool.cpp b/JugBase/src/components/PileupHitMergeTool.cpp new file mode 100644 index 00000000..1601d2da --- /dev/null +++ b/JugBase/src/components/PileupHitMergeTool.cpp @@ -0,0 +1,119 @@ +#include "PileupHitMergeTool.h" + +#include "podio/EventStore.h" + +#include "edm4hep/CaloHitCollection.h" +#include "edm4hep/PositionedCaloHitCollection.h" +#include "edm4hep/PositionedTrackHitCollection.h" +#include "edm4hep/TrackHitCollection.h" + +typedef PileupHitMergeTool PileupCaloHitMergeTool; +typedef PileupHitMergeTool PileupTrackHitMergeTool; +DECLARE_COMPONENT(PileupTrackHitMergeTool) +DECLARE_COMPONENT_WITH_ID(PileupTrackHitMergeTool, "PileupTrackHitMergeTool") +DECLARE_COMPONENT(PileupCaloHitMergeTool) +DECLARE_COMPONENT_WITH_ID(PileupCaloHitMergeTool, "PileupCaloHitMergeTool") + +template +PileupHitMergeTool::PileupHitMergeTool(const std::string& aType, const std::string& aName, + const IInterface* aParent) + : GaudiTool(aType, aName, aParent) { + declareInterface(this); + declareProperty("signalHits", m_hitsSignal); + declareProperty("signalPositionedHits", m_posHitsSignal); + declareProperty("mergedHits", m_hitsMerged); + declareProperty("mergedPositionedHits", m_posHitsMerged); +} + +template +StatusCode PileupHitMergeTool::initialize() { + return StatusCode::SUCCESS; +} + +template +StatusCode PileupHitMergeTool::finalize() { + return StatusCode::SUCCESS; +} + +template +StatusCode PileupHitMergeTool::readPileupCollection(podio::EventStore& store) { + // local pointers, to be filled by the event store + const Hits* hitCollection; + const PositionedHits* posHitCollection; + + // get collection address and store it in container + bool hitCollectionPresent = store.get(m_pileupHitsBranchName, hitCollection); + if (hitCollectionPresent) { + m_hitCollections.push_back(hitCollection); + } else { + warning() << "No collection could be read from branch " << m_pileupHitsBranchName << endmsg; + return StatusCode::FAILURE; + } + + /// as above, for the positioned collection + bool posHitCollectionPresent = store.get(m_pileupPosHitsBranchName, posHitCollection); + if (posHitCollectionPresent) { + m_posHitCollections.push_back(posHitCollection); + } else { + warning() << "No collection could be read from branch " << m_pileupPosHitsBranchName << endmsg; + return StatusCode::FAILURE; + } + + return StatusCode::SUCCESS; +} + +template +StatusCode PileupHitMergeTool::readSignal() { + // get collection from event sture + auto collHitsSig = m_hitsSignal.get(); + auto collPosHitsSig = m_posHitsSignal.get(); + + // store them in internal container + m_hitCollections.push_back(collHitsSig); + m_posHitCollections.push_back(collPosHitsSig); + + return StatusCode::SUCCESS; +} + +template +StatusCode PileupHitMergeTool::mergeCollections() { + + // ownership given to data service at end of execute + Hits* collHitsMerged = new Hits(); + PositionedHits* collPosHitsMerged = new PositionedHits(); + + unsigned int collectionCounter = 0; + for (auto hitColl : m_hitCollections) { + // copy hits + for (const auto elem : *hitColl) { + + auto clon = elem.clone(); + // add pileup vertex counter with an offset + // i.e. for the signal event, 'bits' is just the trackID taken from geant + // for the n-th pileup event, 'bits' is the trackID + n * offset + // offset needs to be big enough to ensure uniqueness of trackID + if (elem.bits() > m_trackIDCollectionOffset) { + error() << "Event contains too many tracks to guarantee a unique trackID"; + error() << " The offset width or trackID field size needs to be adjusted!" << endmsg; + return StatusCode::FAILURE; + } + + clon.bits(clon.bits() + collectionCounter * m_trackIDCollectionOffset); + collHitsMerged->push_back(clon); + } + ++collectionCounter; + } + for (auto posHitColl : m_posHitCollections) { + // copy positioned hits + for (const auto elem : *posHitColl) { + collPosHitsMerged->push_back(elem.clone()); + } + } + + m_hitsMerged.put(collHitsMerged); + m_posHitsMerged.put(collPosHitsMerged); + + m_hitCollections.clear(); + m_posHitCollections.clear(); + return StatusCode::SUCCESS; +} diff --git a/JugBase/src/components/PileupHitMergeTool.h b/JugBase/src/components/PileupHitMergeTool.h new file mode 100644 index 00000000..4a3894ff --- /dev/null +++ b/JugBase/src/components/PileupHitMergeTool.h @@ -0,0 +1,71 @@ +#ifndef FWCORE_PILEUPTRACKERHITSMERGETOOL_H +#define FWCORE_PILEUPTRACKERHITSMERGETOOL_H + +// Gaudi +#include "GaudiAlg/GaudiTool.h" + +// FCCSW +#include "JugBase/DataHandle.h" +#include "JugBase/IEDMMergeTool.h" + +namespace fcc { +class HitCollection; +class PositionedHitCollection; +} + +/** @class PileupHitMergeTool + * + * Implemenation of the MergeTool for *Hits and *PositionedHits, templated to give versions for Tracker / Calorimeter + * While merging, this algorithm tries to keep the trackIDs unique by adding the pileup event number with an offset. + * This should be transparent, but the trackIDs will be non-consecutive. + * + * + */ + +template +class PileupHitMergeTool : public GaudiTool, virtual public IEDMMergeTool { +public: + explicit PileupHitMergeTool(const std::string& aType, const std::string& aName, const IInterface* aParent); + virtual ~PileupHitMergeTool() = default; + + virtual StatusCode initialize() override final; + + virtual StatusCode finalize() override final; + + /// fill the member containers of collection pointers in this class with a collection from the event store + virtual StatusCode readPileupCollection(podio::EventStore& store) override final; + + /// merge the collections in the member container and put them on the event store + virtual StatusCode mergeCollections() override final; + + /// fill the member container of collection pointer with a collection from the event store + virtual StatusCode readSignal() override final; + +private: + /// Name of the branch from which to read pileup collection + Gaudi::Property m_pileupHitsBranchName{this, "pileupHitsBranch", "hits", + "Name of the branch from which to read pileup collection"}; + /// Name of the branch from which to read pileup collection + Gaudi::Property m_pileupPosHitsBranchName{this, "pileupPositionedHitsBranch", "positionedHits", + "Name of the branch from which to read pileup collection"}; + + /// internal container to keep of collections to be merged + std::vector m_hitCollections; + /// internal container to keep of collections to be merged + std::vector m_posHitCollections; + + /// Output of this tool: merged collection + DataHandle m_hitsMerged{"overlay/hits", Gaudi::DataHandle::Writer, this}; + /// Output of this tool: merged collection + DataHandle m_posHitsMerged{"overlay/positionedHits", Gaudi::DataHandle::Writer, this}; + + /// input to this tool: signal collection + DataHandle m_hitsSignal{"overlay/signalHits", Gaudi::DataHandle::Reader, this}; + /// input to this tool: signal collection + DataHandle m_posHitsSignal{"overlay/signalPositionedHits", Gaudi::DataHandle::Reader, this}; + + /// offset with which the pileup event number is added to the trackID + const unsigned int m_trackIDCollectionOffset = 2 << 20; +}; + +#endif // FWCORE_PILEUPTRACKERHITSMERGETOOL_H diff --git a/JugBase/src/components/PileupOverlayAlg.cpp b/JugBase/src/components/PileupOverlayAlg.cpp new file mode 100644 index 00000000..82e86b5f --- /dev/null +++ b/JugBase/src/components/PileupOverlayAlg.cpp @@ -0,0 +1,97 @@ +#include "PileupOverlayAlg.h" + +#include "JugBase/IEDMMergeTool.h" + + + + +DECLARE_COMPONENT(PileupOverlayAlg) + +PileupOverlayAlg::PileupOverlayAlg(const std::string& aName, ISvcLocator* aSvcLoc) + : GaudiAlgorithm(aName, aSvcLoc), m_store(), m_reader() { + declareProperty("PileUpTool", m_pileUpTool); +} + +StatusCode PileupOverlayAlg::initialize() { + + for (auto& toolname : m_saveToolNames) { + m_mergeTools.push_back(tool(toolname)); + } + m_minBiasEventIndex = 0; + StatusCode sc = GaudiAlgorithm::initialize(); + IRndmGenSvc* randSvc = svc("RndmGenSvc", true); + sc = m_flatDist.initialize(randSvc, Rndm::Flat(0., 1.)); + + m_pileupFileIndex = 0; + if (m_doShuffleInputFiles) { + m_pileupFileIndex = static_cast(m_flatDist() * m_pileupFilenames.size()); + } + m_reader.openFiles({m_pileupFilenames[m_pileupFileIndex]}); + m_store.setReader(&m_reader); + return sc; +} + +StatusCode PileupOverlayAlg::finalize() { + StatusCode sc = GaudiAlgorithm::finalize(); + return sc; +} + +StatusCode PileupOverlayAlg::execute() { + unsigned nEvents = m_reader.getEntries(); + + if (!m_noSignal) { + for (auto& tool : m_mergeTools) { + auto sc = tool->readSignal(); + if (sc.isFailure()) { + return StatusCode::FAILURE; + } + } + } + + const unsigned int numPileUp = m_pileUpTool->numberOfPileUp(); + for (unsigned iev = 0; iev < numPileUp; ++iev) { + + // introduce some randomness into min bias event selection + // to avoid bias + if (m_randomizePileup == true) { + if (m_flatDist() < m_skipProbability) { + ++m_minBiasEventIndex; + } + } + m_reader.goToEvent(m_minBiasEventIndex); + + verbose() << "Reading in pileup event #" << m_minBiasEventIndex << " from pool #" << m_pileupFileIndex << " ..." + << endmsg; + for (auto& tool : m_mergeTools) { + auto sc = tool->readPileupCollection(m_store); + if (sc.isFailure()) { + return StatusCode::FAILURE; + } + } + + m_minBiasEventIndex = (m_minBiasEventIndex + 1); + if (m_minBiasEventIndex >= nEvents) { + m_minBiasEventIndex = 0; + if (m_doShuffleInputFiles) { + m_pileupFileIndex = static_cast(m_flatDist() * m_pileupFilenames.size()); + } else { + m_pileupFileIndex = (m_pileupFileIndex + 1) % m_pileupFilenames.size(); + } + m_store.clearCaches(); + m_reader.closeFiles(); + verbose() << "switching to pileup file " << m_pileupFilenames[m_pileupFileIndex] << endmsg; + m_reader.openFiles({m_pileupFilenames[m_pileupFileIndex]}); + nEvents = m_reader.getEntries(); + } + m_store.clearCaches(); + } + + for (auto& tool : m_mergeTools) { + auto sc = tool->mergeCollections(); + if (sc.isFailure()) { + return StatusCode::FAILURE; + } + } + + return StatusCode::SUCCESS; +} diff --git a/JugBase/src/components/PileupOverlayAlg.h b/JugBase/src/components/PileupOverlayAlg.h new file mode 100644 index 00000000..81d55b4b --- /dev/null +++ b/JugBase/src/components/PileupOverlayAlg.h @@ -0,0 +1,73 @@ +#ifndef FWCORE_PILEUPOVERLAYALG_H +#define FWCORE_PILEUPOVERLAYALG_H + +#include "JugBase/DataHandle.h" +#include "JugBase/IEDMMergeTool.h" +#include "JugBase/IPileUpTool.h" +#include "GaudiAlg/GaudiAlgorithm.h" +#include "GaudiKernel/RndmGenerators.h" + +#include "podio/EventStore.h" +#include "podio/ROOTReader.h" + +// forward declarations +class IEDMMergeTool; + +/*** @class PileupOverlayAlg + * Algorithm for Pileup I/O and merging + * + * The main purpose of this class is to keep a list of merge tools corresponding to different + * datatypes, open a file and let the tools read pileup data and merge it with the signal. + * It is thus flexible with regards to the data collection types, as everything is delegated to + * tools with a very general interface. One thing this algorithm does do is keep track of the + * number of pileup events (via the pileup tool) and the current position in the minimum bias pool + * ( which can be randomized by skipping events with a probability that can be set in the options file ). + * + */ + +class PileupOverlayAlg : public GaudiAlgorithm { +public: + explicit PileupOverlayAlg(const std::string&, ISvcLocator*); + virtual ~PileupOverlayAlg() = default; + virtual StatusCode initialize() override final; + virtual StatusCode execute() override final; + virtual StatusCode finalize() override final; + +private: + /// flag to determine whether to randomly skipy events when reading the pileup store + Gaudi::Property m_randomizePileup{ + this, "randomizePileup", false, + "flag to determine whether to randomly skip events when reading the pileup store"}; + /// probability of skipping next event in randomized read + Gaudi::Property m_skipProbability{this, "skipProbability", 0.5, "Probability with which to skip next event"}; + /// random number generator for randomizing pileup reads + Rndm::Numbers m_flatDist; + /// Pileup Tool to specify the number of minimum bias events + ToolHandle m_pileUpTool{"ConstPileUp/PileUpTool", this}; + + /// vector of tool interface pointers, to be replaced with ToolHandleArray + std::vector m_mergeTools; + /// names of the merge tools + Gaudi::Property> m_saveToolNames{this, "mergeTools", {}, "names of the merge tools"}; + + // list of minimum bias data readers / stores + // once the root reader supports multiple files + // this algo should take advantage of it + /// list of filenames for the minimum bias pools + Gaudi::Property> m_pileupFilenames{ + this, "pileupFilenames", {"min_bias_pool.root"}, "Name of File containing pileup events"}; + Gaudi::Property m_doShuffleInputFiles{this, "doShuffleInputFiles", false, "Shuffle list of input files for additional randomization"}; + + Gaudi::Property m_noSignal{this, "noSignal", false, "Set to true if you don't want to provide a signal collection"}; + /// store for the minimum bias file + podio::EventStore m_store; + /// reader for the minimum bias file + podio::ROOTReader m_reader; + + /// event index within current minimum bias pool + unsigned int m_minBiasEventIndex; + /// index of the current minimum bias pool + unsigned int m_pileupFileIndex; +}; + +#endif /* FWCORE_PILEUPOVERLAYALG_H */ diff --git a/JugBase/src/components/PileupParticlesMergeTool.cpp b/JugBase/src/components/PileupParticlesMergeTool.cpp new file mode 100644 index 00000000..3bb9cd10 --- /dev/null +++ b/JugBase/src/components/PileupParticlesMergeTool.cpp @@ -0,0 +1,106 @@ +#include "PileupParticlesMergeTool.h" + +#include "podio/EventStore.h" + +#include "edm4hep/GenVertexCollection.h" +#include "edm4hep/MCParticleCollection.h" + +DECLARE_COMPONENT(PileupParticlesMergeTool) + +PileupParticlesMergeTool::PileupParticlesMergeTool(const std::string& aType, const std::string& aName, + const IInterface* aParent) + : GaudiTool(aType, aName, aParent) { + declareInterface(this); + declareProperty("signalGenVertices", m_vertIn); + declareProperty("signalGenParticles", m_partIn); + declareProperty("mergedGenParticles", m_partOut); + declareProperty("mergedGenVertices", m_vertOut); +} + +StatusCode PileupParticlesMergeTool::initialize() { return StatusCode::SUCCESS; } + +StatusCode PileupParticlesMergeTool::finalize() { return StatusCode::SUCCESS; } + +StatusCode PileupParticlesMergeTool::readPileupCollection(podio::EventStore& store) { + const edm4hep::MCParticleCollection* mcParticleCollection; + const edm4hep::GenVertexCollection* genVertexCollection; + + bool mcParticleCollectionPresent = store.get(m_pileupGenParticlesBranchName, mcParticleCollection); + debug() << "size of collection about to be read " << mcParticleCollection->size() << endmsg; + if (mcParticleCollectionPresent) { + m_MCParticleCollections.push_back(mcParticleCollection); + } else { + warning() << "No collection could be read from branch " << m_pileupGenParticlesBranchName << endmsg; + } + + bool genVertexCollectionPresent = store.get(m_pileupGenVerticesBranchName, genVertexCollection); + if (genVertexCollectionPresent) { + m_GenVertexCollections.push_back(genVertexCollection); + } else { + warning() << "No collection could be read from branch " << m_pileupGenVerticesBranchName << endmsg; + } + + return StatusCode::SUCCESS; +} + +StatusCode PileupParticlesMergeTool::readSignal() { + auto collVSig = m_vertIn.get(); + auto collPSig = m_partIn.get(); + + m_MCParticleCollections.push_back(collPSig); + m_GenVertexCollections.push_back(collVSig); + + return StatusCode::SUCCESS; +} + +StatusCode PileupParticlesMergeTool::mergeCollections() { + + debug() << "merge " << m_GenVertexCollections.size() << " collections ..." << endmsg; + + // ownership given to data service at end of execute + edm4hep::MCParticleCollection* collPOut = new edm4hep::MCParticleCollection(); + edm4hep::GenVertexCollection* collVOut = new edm4hep::GenVertexCollection(); + + // need to keep track of the accumulated length + // of collections already merged + unsigned int collectionOffset = 0; + + // merge vertices + + for (auto genVertexColl : m_GenVertexCollections) { + debug() << " traverse collection of size " << genVertexColl->size() << endmsg; + for (const auto elem : *genVertexColl) { + collVOut->push_back(elem.clone()); + } + } + // merge particles, keeping the references up to date + for (unsigned int collCounter = 0; collCounter < m_MCParticleCollections.size(); ++collCounter) { + auto mcPartColl = m_MCParticleCollections[collCounter]; + for (const auto elem : *mcPartColl) { + auto newPart = edm4hep::MCParticle(elem.core()); + if (elem.startVertex().isAvailable()) { + // update reference: find new index in merged collection + // add offset since the signal particles are copied in new collection first + unsigned int newIdStart = elem.startVertex().getObjectID().index + collectionOffset; + // update startVertex + newPart.startVertex((*collVOut)[newIdStart]); + } + if (elem.endVertex().isAvailable()) { + // update reference: find new index in merged collection + // add offset since the signal particles are copied in new collection first + unsigned int newIdEnd = elem.endVertex().getObjectID().index + collectionOffset; + // update startVertex + newPart.endVertex((*collVOut)[newIdEnd]); + } + collPOut->push_back(newPart); + } + collectionOffset += m_GenVertexCollections[collCounter]->size(); + } + + m_vertOut.put(collVOut); + m_partOut.put(collPOut); + + m_MCParticleCollections.clear(); + m_GenVertexCollections.clear(); + return StatusCode::SUCCESS; +} diff --git a/JugBase/src/components/PileupParticlesMergeTool.h b/JugBase/src/components/PileupParticlesMergeTool.h new file mode 100644 index 00000000..ef8d186b --- /dev/null +++ b/JugBase/src/components/PileupParticlesMergeTool.h @@ -0,0 +1,51 @@ +#ifndef FWCORE_PILEUPPARTICLESMERGETOOL_H +#define FWCORE_PILEUPPARTICLESMERGETOOL_H + +// Gaudi +#include "GaudiAlg/GaudiTool.h" + +// FCCSW +#include "JugBase/DataHandle.h" +#include "JugBase/IEDMMergeTool.h" + +namespace fcc { +class MCParticleCollection; +class GenVertexCollection; +} + +/** @class PileupParticlesMergeTool + * + * Implementation of the EDMMergeTool interface used for pileup overlay + * of generated particles. Keeps track of the association between particles and vertices. + */ + +class PileupParticlesMergeTool : public GaudiTool, virtual public IEDMMergeTool { +public: + explicit PileupParticlesMergeTool(const std::string& aType, const std::string& aName, const IInterface* aParent); + virtual ~PileupParticlesMergeTool() = default; + + virtual StatusCode initialize() override final; + + virtual StatusCode finalize() override final; + + virtual StatusCode readPileupCollection(podio::EventStore& store) override final; + + virtual StatusCode mergeCollections() override final; + + virtual StatusCode readSignal() override final; + +private: + Gaudi::Property m_pileupGenParticlesBranchName{this, "genParticlesBranch", "genParticles", ""}; + Gaudi::Property m_pileupGenVerticesBranchName{this, "genVerticesBranch", "genVertices", ""}; + + std::vector m_MCParticleCollections; + std::vector m_GenVertexCollections; + + DataHandle m_vertOut{"overlay/allGenVertices", Gaudi::DataHandle::Writer, this}; + DataHandle m_partOut{"overlay/allGenParticles", Gaudi::DataHandle::Writer, this}; + + DataHandle m_vertIn{"overlay/signalGenVertices", Gaudi::DataHandle::Reader, this}; + DataHandle m_partIn{"overlay/signalGenParticles", Gaudi::DataHandle::Reader, this}; +}; + +#endif // FWCORE_PILEUPPARTICLESMERGETOOL_H diff --git a/JugBase/src/components/PoissonPileUp.cpp b/JugBase/src/components/PoissonPileUp.cpp new file mode 100644 index 00000000..073dfe2a --- /dev/null +++ b/JugBase/src/components/PoissonPileUp.cpp @@ -0,0 +1,46 @@ +#include "PoissonPileUp.h" +#include "GaudiKernel/IRndmGenSvc.h" + +DECLARE_COMPONENT(PoissonPileUp) + +PoissonPileUp::PoissonPileUp(const std::string& type, const std::string& name, const IInterface* parent) + : GaudiTool(type, name, parent) { + declareInterface(this); +} + +PoissonPileUp::~PoissonPileUp() { ; } + +StatusCode PoissonPileUp::initialize() { + StatusCode sc = GaudiTool::initialize(); + if (sc.isFailure()) { + return sc; + } + + IRndmGenSvc* randSvc = svc("RndmGenSvc", true); + if (randSvc == nullptr) { + return StatusCode::FAILURE; + } + + if (m_meanPileUpEvents < 0) return Error("Number of Pileup events cannot be negative!"); + sc = m_PoissonDist.initialize(randSvc, Rndm::Poisson(m_meanPileUpEvents)); + if (!sc.isSuccess()) return Error("Could not initialize Poisson random number generator"); + + if (release(randSvc).isFailure()) { + return StatusCode::FAILURE; + } + + m_currentNumPileUpEvents = m_PoissonDist(); + printPileUpCounters(); + return sc; +} + +unsigned int PoissonPileUp::numberOfPileUp() { + m_currentNumPileUpEvents = m_PoissonDist(); + return m_currentNumPileUpEvents; +} + +double PoissonPileUp::getMeanPileUp() { return m_meanPileUpEvents; } + +void PoissonPileUp::printPileUpCounters() { + info() << "Current number of pileup events: " << m_currentNumPileUpEvents << endmsg; +} diff --git a/JugBase/src/components/PoissonPileUp.h b/JugBase/src/components/PoissonPileUp.h new file mode 100644 index 00000000..b07402dc --- /dev/null +++ b/JugBase/src/components/PoissonPileUp.h @@ -0,0 +1,37 @@ +#ifndef GENERATION_POISSONPILEUP_H +#define GENERATION_POISSONPILEUP_H + +#include "JugBase/IPileUpTool.h" +#include "GaudiAlg/GaudiTool.h" +#include "GaudiKernel/RndmGenerators.h" + +/** @class PoissonPileUp + * + * Tool to generate number of pile-up events to be mixed with signal event. + * Concrete implementation of a IPileUpTool, returning a random variable + * drawn from a Poisson distribution centered around the average number + * of pileup events. + * + * @author Valentin Volkl + * @date 2016-01-18 + */ +class PoissonPileUp : public GaudiTool, virtual public Jug::Base::IPileUpTool { +public: + PoissonPileUp(const std::string& type, const std::string& name, const IInterface* parent); + virtual ~PoissonPileUp(); + virtual StatusCode initialize(); + virtual unsigned int numberOfPileUp(); + virtual double getMeanPileUp(); + virtual void printPileUpCounters(); + + +private: + /// average number of min bias events to pile on signal event. + Gaudi::Property m_meanPileUpEvents{this, "numPileUpEvents", 0, "Average number of pile-up events"}; + /// holds last realization of the random variable + unsigned int m_currentNumPileUpEvents; + /// Poisson random number generator + Rndm::Numbers m_PoissonDist; +}; + +#endif // GENERATION_POISSONPILEUP_H diff --git a/JugBase/src/components/RangePileUp.cpp b/JugBase/src/components/RangePileUp.cpp new file mode 100644 index 00000000..bcdca30d --- /dev/null +++ b/JugBase/src/components/RangePileUp.cpp @@ -0,0 +1,28 @@ +#include "RangePileUp.h" + +DECLARE_COMPONENT(RangePileUp) + +RangePileUp::RangePileUp(const std::string& type, const std::string& name, const IInterface* parent) + : GaudiTool(type, name, parent) { + declareInterface(this); +} + +RangePileUp::~RangePileUp() { ; } + +StatusCode RangePileUp::initialize() { + StatusCode sc = GaudiTool::initialize(); + if (sc.isFailure()) return sc; + return sc; +} + +unsigned int RangePileUp::numberOfPileUp() { + m_currentNumPileUpEvents = m_pileupRange[m_rangeIndex]; + m_rangeIndex = (m_rangeIndex + 1) % m_pileupRange.size(); + return m_currentNumPileUpEvents; +} + +double RangePileUp::getMeanPileUp() { return m_currentNumPileUpEvents; } + +void RangePileUp::printPileUpCounters() { + info() << "Current number of pileup events: " << m_currentNumPileUpEvents << endmsg; +} diff --git a/JugBase/src/components/RangePileUp.h b/JugBase/src/components/RangePileUp.h new file mode 100644 index 00000000..c6263913 --- /dev/null +++ b/JugBase/src/components/RangePileUp.h @@ -0,0 +1,33 @@ +#ifndef GENERATION_RANGEPILEUP_H +#define GENERATION_RANGEPILEUP_H + +#include "JugBase/IPileUpTool.h" +#include "GaudiAlg/GaudiTool.h" + +/** @class RangePileUp + * + * Tool to generate number of pile-up events to be mixed with signal event. + * Concrete implementation of a IPileUpTool, returning a random variable + * drawn from a pre-defined list - useful for detailed detector studies of + * the effects of increasing pile-up. + * + * @author Valentin Volkl + * @date 2016-01-18 + */ +class RangePileUp : public GaudiTool, virtual public Jug::Base::IPileUpTool { +public: + RangePileUp(const std::string& type, const std::string& name, const IInterface* parent); + virtual ~RangePileUp(); + virtual StatusCode initialize(); + virtual unsigned int numberOfPileUp(); + virtual double getMeanPileUp(); + virtual void printPileUpCounters(); + +private: + /// average number of min bias events to pile on signal event. + Gaudi::Property> m_pileupRange{this, "numPileUpEvents", {0}, "Average number of pile-up events"}; + unsigned int m_currentNumPileUpEvents; + unsigned int m_rangeIndex; +}; + +#endif // GENERATION_RANGEPILEUP_H -- GitLab From 3dfe1feb39bdf75bc17d7cc30a531d9d1495c600 Mon Sep 17 00:00:00 2001 From: Wouter Deconinck Date: Mon, 11 Jul 2022 19:25:13 +0200 Subject: [PATCH 2/3] Remove DigiTrackHit and Particles mergers; fix PileupHitMergeTool compilation --- .../PileupDigiTrackHitMergeTool.cpp | 162 ------------------ .../components/PileupDigiTrackHitMergeTool.h | 99 ----------- JugBase/src/components/PileupHitMergeTool.cpp | 65 ++----- JugBase/src/components/PileupHitMergeTool.h | 20 +-- .../components/PileupParticlesMergeTool.cpp | 106 ------------ .../src/components/PileupParticlesMergeTool.h | 51 ------ 6 files changed, 22 insertions(+), 481 deletions(-) delete mode 100644 JugBase/src/components/PileupDigiTrackHitMergeTool.cpp delete mode 100644 JugBase/src/components/PileupDigiTrackHitMergeTool.h delete mode 100644 JugBase/src/components/PileupParticlesMergeTool.cpp delete mode 100644 JugBase/src/components/PileupParticlesMergeTool.h diff --git a/JugBase/src/components/PileupDigiTrackHitMergeTool.cpp b/JugBase/src/components/PileupDigiTrackHitMergeTool.cpp deleted file mode 100644 index 4c88b226..00000000 --- a/JugBase/src/components/PileupDigiTrackHitMergeTool.cpp +++ /dev/null @@ -1,162 +0,0 @@ -#include "PileupDigiTrackHitMergeTool.h" - -#include "podio/EventStore.h" - -#include "edm4hep/DigiTrackHitAssociationCollection.h" -#include "edm4hep/MCParticle.h" -#include "edm4hep/MCParticleCollection.h" -#include "edm4hep/TrackHitCollection.h" - -DECLARE_COMPONENT(PileupDigiTrackHitMergeTool) -DECLARE_COMPONENT_WITH_ID(PileupDigiTrackHitMergeTool, "PileupDigiTrackHitMergeTool") - -PileupDigiTrackHitMergeTool::PileupDigiTrackHitMergeTool(const std::string& aType, const std::string& aName, - const IInterface* aParent) - : GaudiTool(aType, aName, aParent) { - declareProperty("signalPositionedHits", m_hitsSignal); - declareProperty("signalDigiHits", m_posHitsSignal); - declareProperty("signalParticles", m_particlesSignal); - declareProperty("mergedPositionedHits", m_hitsMerged); - declareProperty("mergedDigiHits", m_digiHitsMerged); - declareProperty("mergedParticles", m_particlesMerged); -} - -StatusCode PileupDigiTrackHitMergeTool::initialize() { return StatusCode::SUCCESS; } - -StatusCode PileupDigiTrackHitMergeTool::finalize() { return StatusCode::SUCCESS; } - -StatusCode PileupDigiTrackHitMergeTool::readPileupCollection(podio::EventStore& store) { - // local pointers, to be filled by the event store - const edm4hep::PositionedTrackHitCollection* hitCollection; - const edm4hep::DigiTrackHitAssociationCollection* digiHitCollection; - const edm4hep::MCParticleCollection* particleCollection; - - // get collection address and store it in container - bool hitCollectionPresent = store.get(m_pileupHitsBranchName, hitCollection); - if (hitCollectionPresent) { - m_hitCollections.push_back(hitCollection); - } else { - warning() << "No collection could be read from branch " << m_pileupHitsBranchName << endmsg; - return StatusCode::FAILURE; - } - - /// as above, for the positioned collection - bool digiHitCollectionPresent = store.get(m_pileupPosHitsBranchName, digiHitCollection); - if (digiHitCollectionPresent) { - m_digiTrackHitsCollections.push_back(digiHitCollection); - } else { - warning() << "No collection could be read from branch " << m_pileupPosHitsBranchName << endmsg; - return StatusCode::FAILURE; - } - - /// as above, for the particle collection (if wanted) - if (m_mergeParticles) { - bool particleCollectionPresent = store.get(m_pileupParticleBranchName, particleCollection); - if (particleCollectionPresent) { - m_particleCollections.push_back(particleCollection); - } else { - warning() << "No collection could be read from branch " << m_pileupParticleBranchName << endmsg; - return StatusCode::FAILURE; - } - } - - return StatusCode::SUCCESS; -} - -StatusCode PileupDigiTrackHitMergeTool::readSignal() { - // get collection from event store - auto collHitsSig = m_hitsSignal.get(); - auto collPosHitsSig = m_posHitsSignal.get(); - auto collPartSig = m_particlesSignal.get(); - // store them in internal container - m_hitCollections.push_back(collHitsSig); - m_digiTrackHitsCollections.push_back(collPosHitsSig); - // check if particles should be merged - if (m_mergeParticles) { - m_particleCollections.push_back(collPartSig); - } - - return StatusCode::SUCCESS; -} - -StatusCode PileupDigiTrackHitMergeTool::mergeCollections() { - - // ownership given to data service at end of execute - edm4hep::PositionedTrackHitCollection* collHitsMerged = new edm4hep::PositionedTrackHitCollection(); - edm4hep::DigiTrackHitAssociationCollection* collDigiHitsMerged = new edm4hep::DigiTrackHitAssociationCollection(); - edm4hep::MCParticleCollection* collParticlesMerged = new edm4hep::MCParticleCollection(); - - unsigned int collectionCounter = 0; - for (size_t j = 0; j < m_digiTrackHitsCollections.size(); j++) { - unsigned int trackIDCounter = 0; - std::unordered_map oldToNewTrackIDs; - - auto offset = collectionCounter * m_trackIDCollectionOffset; - auto digiHitColl = m_digiTrackHitsCollections.at(j); - auto hitColl = m_hitCollections.at(j); - for (int i = 0; i < digiHitColl->size(); i++) { - auto clon = hitColl->at(i).clone(); - // add pileup vertex counter with an offset - // i.e. for the signal event, 'bits' is just the trackID taken from geant - // for the n-th pileup event, 'bits' is the trackID + n * offset - // offset needs to be big enough to ensure uniqueness of trackID - if (isTrackerHit(hitColl->at(i).cellId())) { - if (trackIDCounter > m_trackIDCollectionOffset) { - error() << "Event contains too many tracks to guarantee a unique trackID"; - error() << " The offset width or trackID field size needs to be adjusted!" << endmsg; - return StatusCode::FAILURE; - } - auto search = oldToNewTrackIDs.find(clon.bits()); - if (search != oldToNewTrackIDs.end()) { - // was already present - use already created new trackID - unsigned newTrackID = (search->second + offset); - clon.bits(newTrackID); - } else { - // was not already present - create new trackID - oldToNewTrackIDs[clon.bits()] = trackIDCounter; - - unsigned newTrackID = (trackIDCounter + offset); - clon.bits(newTrackID); - - trackIDCounter++; - } - } - auto newDigiHit = digiHitColl->at(i).clone(); - newDigiHit.hit(clon); - collHitsMerged->push_back(clon); - collDigiHitsMerged->push_back(newDigiHit); - } - ///@todo find possibility in case not only tracker hits are produced - if (m_mergeParticles) { - auto partColl = m_particleCollections.at(j); - for (int i = 0; i < partColl->size(); i++) { - auto clonParticle = partColl->at(i).clone(); - auto search = oldToNewTrackIDs.find(clonParticle.bits()); - if (search != oldToNewTrackIDs.end()) { - // was not already present - create new trackID - clonParticle.bits(search->second + offset); - } else { - // was not already present - create new trackID - // this can happen if particles do not reach the tracker - oldToNewTrackIDs[clonParticle.bits()] = trackIDCounter; - - unsigned newTrackID = (trackIDCounter + offset); - clonParticle.bits(newTrackID); - - trackIDCounter++; - } - collParticlesMerged->push_back(clonParticle); - } - } - ++collectionCounter; - } - m_hitsMerged.put(collHitsMerged); - m_digiHitsMerged.put(collDigiHitsMerged); - m_particlesMerged.put(collParticlesMerged); - - m_hitCollections.clear(); - m_digiTrackHitsCollections.clear(); - m_particleCollections.clear(); - - return StatusCode::SUCCESS; -} diff --git a/JugBase/src/components/PileupDigiTrackHitMergeTool.h b/JugBase/src/components/PileupDigiTrackHitMergeTool.h deleted file mode 100644 index f59612ff..00000000 --- a/JugBase/src/components/PileupDigiTrackHitMergeTool.h +++ /dev/null @@ -1,99 +0,0 @@ -#ifndef FWCORE_PILEUPDIGITRACKERHITSMERGETOOL_H -#define FWCORE_PILEUPDIGITRACKERHITSMERGETOOL_H - -// Gaudi -#include "GaudiAlg/GaudiTool.h" - -// FCCSW -#include "JugBase/DataHandle.h" -#include "JugBase/IEDMMergeTool.h" - -namespace fcc { -class PositionedTrackHit; -class PositionedTrackHitCollection; -class DigiTrackHitAssociation; -class DigiTrackHitAssociationCollection; -class MCParticle; -class MCParticleCollection; -} - -/** @class PileupDigiTrackHitMergeTool - * - * Implemenation of the MergeTool for Digitization hits and, if requested, simulated hits. - * While merging, this algorithm tries to keep the trackIDs unique by adding the pileup event number with an offset. - * This should be transparent, but the trackIDs will be non-consecutive. - * - */ - -class PileupDigiTrackHitMergeTool : public GaudiTool, virtual public IEDMMergeTool { -public: - explicit PileupDigiTrackHitMergeTool(const std::string& aType, const std::string& aName, const IInterface* aParent); - virtual ~PileupDigiTrackHitMergeTool() = default; - - virtual StatusCode initialize() override final; - - virtual StatusCode finalize() override final; - - /// fill the member containers of collection pointers in this class with a collection from the event store - virtual StatusCode readPileupCollection(podio::EventStore& store) override final; - - /// merge the collections in the member container and put them on the event store - virtual StatusCode mergeCollections() override final; - - /// fill the member container of collection pointer with a collection from the event store - virtual StatusCode readSignal() override final; - -private: - /// Name of the branch from which to read pileup collection - Gaudi::Property m_pileupHitsBranchName{this, "pileupHitsBranch", "positionedHits", - "Name of the branch from which to read pileup collection"}; - /// Name of the branch from which to read pileup collection - Gaudi::Property m_pileupPosHitsBranchName{this, "pileupDigiHitsBranch", "digiHits", - "Name of the branch from which to read pileup collection"}; - - /// Name of the branch from which to read pileup collection - Gaudi::Property m_pileupParticleBranchName{this, "pileupParticleBranch", "simParticles", - "Name of the branch from which to read pileup collection"}; - - /// Flag indicating if particles are present and should be merged - Gaudi::Property m_mergeParticles{this, "mergeParticles", false, - "Flag indicating if particles are present and should be merged"}; - - /// internal container to keep of collections to be merged - std::vector m_hitCollections; - /// internal container to keep of collections to be merged - std::vector m_digiTrackHitsCollections; - /// internal container to keep of collections to be merge - std::vector m_particleCollections; - - /// Output of this tool: merged collection - DataHandle m_hitsMerged{"overlay/hits", Gaudi::DataHandle::Writer, this}; - /// Output of this tool: merged collection - DataHandle m_digiHitsMerged{"overlay/digiHits", Gaudi::DataHandle::Writer, - this}; - /// Output of this tool: merged collection - DataHandle m_particlesMerged{"overlay/particles", Gaudi::DataHandle::Writer, this}; - - /// input to this tool: signal collection - DataHandle m_hitsSignal{"overlay/signalHits", Gaudi::DataHandle::Reader, this}; - /// input to this tool: signal collection - DataHandle m_posHitsSignal{"overlay/signalDigiHits", - Gaudi::DataHandle::Reader, this}; - /// input to this tool: signal collection - DataHandle m_particlesSignal{"overlay/signalParticles", Gaudi::DataHandle::Reader, this}; - - /// offset with which the pileup event number is added to the trackID - const unsigned int m_trackIDCollectionOffset = 2.5e6; - - /// Mask to obtain system ID - const unsigned long long m_systemMask = 0xf; - - /// Private function internally used to check if a hit is a tracker hit - bool isTrackerHit(unsigned long long cellId) const; -}; - -inline bool PileupDigiTrackHitMergeTool::isTrackerHit(unsigned long long cellId) const { - return ((cellId & m_systemMask) < 5); -} - -#endif // FWCORE_PILEUPDIGITRACKERHITSMERGETOOL_H diff --git a/JugBase/src/components/PileupHitMergeTool.cpp b/JugBase/src/components/PileupHitMergeTool.cpp index 1601d2da..c2e5eaaf 100644 --- a/JugBase/src/components/PileupHitMergeTool.cpp +++ b/JugBase/src/components/PileupHitMergeTool.cpp @@ -2,44 +2,36 @@ #include "podio/EventStore.h" -#include "edm4hep/CaloHitCollection.h" -#include "edm4hep/PositionedCaloHitCollection.h" -#include "edm4hep/PositionedTrackHitCollection.h" -#include "edm4hep/TrackHitCollection.h" +#include "edm4hep/SimCalorimeterHitCollection.h" +#include "edm4hep/SimTrackerHitCollection.h" -typedef PileupHitMergeTool PileupCaloHitMergeTool; -typedef PileupHitMergeTool PileupTrackHitMergeTool; +typedef PileupHitMergeTool PileupCaloHitMergeTool; +typedef PileupHitMergeTool PileupTrackHitMergeTool; DECLARE_COMPONENT(PileupTrackHitMergeTool) -DECLARE_COMPONENT_WITH_ID(PileupTrackHitMergeTool, "PileupTrackHitMergeTool") DECLARE_COMPONENT(PileupCaloHitMergeTool) -DECLARE_COMPONENT_WITH_ID(PileupCaloHitMergeTool, "PileupCaloHitMergeTool") -template -PileupHitMergeTool::PileupHitMergeTool(const std::string& aType, const std::string& aName, - const IInterface* aParent) +template +PileupHitMergeTool::PileupHitMergeTool(const std::string& aType, const std::string& aName, + const IInterface* aParent) : GaudiTool(aType, aName, aParent) { - declareInterface(this); declareProperty("signalHits", m_hitsSignal); - declareProperty("signalPositionedHits", m_posHitsSignal); declareProperty("mergedHits", m_hitsMerged); - declareProperty("mergedPositionedHits", m_posHitsMerged); } -template -StatusCode PileupHitMergeTool::initialize() { +template +StatusCode PileupHitMergeTool::initialize() { return StatusCode::SUCCESS; } -template -StatusCode PileupHitMergeTool::finalize() { +template +StatusCode PileupHitMergeTool::finalize() { return StatusCode::SUCCESS; } -template -StatusCode PileupHitMergeTool::readPileupCollection(podio::EventStore& store) { +template +StatusCode PileupHitMergeTool::readPileupCollection(podio::EventStore& store) { // local pointers, to be filled by the event store const Hits* hitCollection; - const PositionedHits* posHitCollection; // get collection address and store it in container bool hitCollectionPresent = store.get(m_pileupHitsBranchName, hitCollection); @@ -50,37 +42,25 @@ StatusCode PileupHitMergeTool::readPileupCollection(podio: return StatusCode::FAILURE; } - /// as above, for the positioned collection - bool posHitCollectionPresent = store.get(m_pileupPosHitsBranchName, posHitCollection); - if (posHitCollectionPresent) { - m_posHitCollections.push_back(posHitCollection); - } else { - warning() << "No collection could be read from branch " << m_pileupPosHitsBranchName << endmsg; - return StatusCode::FAILURE; - } - return StatusCode::SUCCESS; } -template -StatusCode PileupHitMergeTool::readSignal() { +template +StatusCode PileupHitMergeTool::readSignal() { // get collection from event sture auto collHitsSig = m_hitsSignal.get(); - auto collPosHitsSig = m_posHitsSignal.get(); // store them in internal container m_hitCollections.push_back(collHitsSig); - m_posHitCollections.push_back(collPosHitsSig); return StatusCode::SUCCESS; } -template -StatusCode PileupHitMergeTool::mergeCollections() { +template +StatusCode PileupHitMergeTool::mergeCollections() { // ownership given to data service at end of execute Hits* collHitsMerged = new Hits(); - PositionedHits* collPosHitsMerged = new PositionedHits(); unsigned int collectionCounter = 0; for (auto hitColl : m_hitCollections) { @@ -92,6 +72,7 @@ StatusCode PileupHitMergeTool::mergeCollections() { // i.e. for the signal event, 'bits' is just the trackID taken from geant // for the n-th pileup event, 'bits' is the trackID + n * offset // offset needs to be big enough to ensure uniqueness of trackID +/* if (elem.bits() > m_trackIDCollectionOffset) { error() << "Event contains too many tracks to guarantee a unique trackID"; error() << " The offset width or trackID field size needs to be adjusted!" << endmsg; @@ -99,21 +80,13 @@ StatusCode PileupHitMergeTool::mergeCollections() { } clon.bits(clon.bits() + collectionCounter * m_trackIDCollectionOffset); +*/ collHitsMerged->push_back(clon); } ++collectionCounter; } - for (auto posHitColl : m_posHitCollections) { - // copy positioned hits - for (const auto elem : *posHitColl) { - collPosHitsMerged->push_back(elem.clone()); - } - } m_hitsMerged.put(collHitsMerged); - m_posHitsMerged.put(collPosHitsMerged); - m_hitCollections.clear(); - m_posHitCollections.clear(); return StatusCode::SUCCESS; } diff --git a/JugBase/src/components/PileupHitMergeTool.h b/JugBase/src/components/PileupHitMergeTool.h index 4a3894ff..22d5fcc8 100644 --- a/JugBase/src/components/PileupHitMergeTool.h +++ b/JugBase/src/components/PileupHitMergeTool.h @@ -4,25 +4,20 @@ // Gaudi #include "GaudiAlg/GaudiTool.h" -// FCCSW +// JugBase #include "JugBase/DataHandle.h" #include "JugBase/IEDMMergeTool.h" -namespace fcc { -class HitCollection; -class PositionedHitCollection; -} - /** @class PileupHitMergeTool * - * Implemenation of the MergeTool for *Hits and *PositionedHits, templated to give versions for Tracker / Calorimeter + * Implemenation of the MergeTool for *Hits, templated to give versions for Tracker / Calorimeter * While merging, this algorithm tries to keep the trackIDs unique by adding the pileup event number with an offset. * This should be transparent, but the trackIDs will be non-consecutive. * * */ -template +template class PileupHitMergeTool : public GaudiTool, virtual public IEDMMergeTool { public: explicit PileupHitMergeTool(const std::string& aType, const std::string& aName, const IInterface* aParent); @@ -45,24 +40,15 @@ private: /// Name of the branch from which to read pileup collection Gaudi::Property m_pileupHitsBranchName{this, "pileupHitsBranch", "hits", "Name of the branch from which to read pileup collection"}; - /// Name of the branch from which to read pileup collection - Gaudi::Property m_pileupPosHitsBranchName{this, "pileupPositionedHitsBranch", "positionedHits", - "Name of the branch from which to read pileup collection"}; /// internal container to keep of collections to be merged std::vector m_hitCollections; - /// internal container to keep of collections to be merged - std::vector m_posHitCollections; /// Output of this tool: merged collection DataHandle m_hitsMerged{"overlay/hits", Gaudi::DataHandle::Writer, this}; - /// Output of this tool: merged collection - DataHandle m_posHitsMerged{"overlay/positionedHits", Gaudi::DataHandle::Writer, this}; /// input to this tool: signal collection DataHandle m_hitsSignal{"overlay/signalHits", Gaudi::DataHandle::Reader, this}; - /// input to this tool: signal collection - DataHandle m_posHitsSignal{"overlay/signalPositionedHits", Gaudi::DataHandle::Reader, this}; /// offset with which the pileup event number is added to the trackID const unsigned int m_trackIDCollectionOffset = 2 << 20; diff --git a/JugBase/src/components/PileupParticlesMergeTool.cpp b/JugBase/src/components/PileupParticlesMergeTool.cpp deleted file mode 100644 index 3bb9cd10..00000000 --- a/JugBase/src/components/PileupParticlesMergeTool.cpp +++ /dev/null @@ -1,106 +0,0 @@ -#include "PileupParticlesMergeTool.h" - -#include "podio/EventStore.h" - -#include "edm4hep/GenVertexCollection.h" -#include "edm4hep/MCParticleCollection.h" - -DECLARE_COMPONENT(PileupParticlesMergeTool) - -PileupParticlesMergeTool::PileupParticlesMergeTool(const std::string& aType, const std::string& aName, - const IInterface* aParent) - : GaudiTool(aType, aName, aParent) { - declareInterface(this); - declareProperty("signalGenVertices", m_vertIn); - declareProperty("signalGenParticles", m_partIn); - declareProperty("mergedGenParticles", m_partOut); - declareProperty("mergedGenVertices", m_vertOut); -} - -StatusCode PileupParticlesMergeTool::initialize() { return StatusCode::SUCCESS; } - -StatusCode PileupParticlesMergeTool::finalize() { return StatusCode::SUCCESS; } - -StatusCode PileupParticlesMergeTool::readPileupCollection(podio::EventStore& store) { - const edm4hep::MCParticleCollection* mcParticleCollection; - const edm4hep::GenVertexCollection* genVertexCollection; - - bool mcParticleCollectionPresent = store.get(m_pileupGenParticlesBranchName, mcParticleCollection); - debug() << "size of collection about to be read " << mcParticleCollection->size() << endmsg; - if (mcParticleCollectionPresent) { - m_MCParticleCollections.push_back(mcParticleCollection); - } else { - warning() << "No collection could be read from branch " << m_pileupGenParticlesBranchName << endmsg; - } - - bool genVertexCollectionPresent = store.get(m_pileupGenVerticesBranchName, genVertexCollection); - if (genVertexCollectionPresent) { - m_GenVertexCollections.push_back(genVertexCollection); - } else { - warning() << "No collection could be read from branch " << m_pileupGenVerticesBranchName << endmsg; - } - - return StatusCode::SUCCESS; -} - -StatusCode PileupParticlesMergeTool::readSignal() { - auto collVSig = m_vertIn.get(); - auto collPSig = m_partIn.get(); - - m_MCParticleCollections.push_back(collPSig); - m_GenVertexCollections.push_back(collVSig); - - return StatusCode::SUCCESS; -} - -StatusCode PileupParticlesMergeTool::mergeCollections() { - - debug() << "merge " << m_GenVertexCollections.size() << " collections ..." << endmsg; - - // ownership given to data service at end of execute - edm4hep::MCParticleCollection* collPOut = new edm4hep::MCParticleCollection(); - edm4hep::GenVertexCollection* collVOut = new edm4hep::GenVertexCollection(); - - // need to keep track of the accumulated length - // of collections already merged - unsigned int collectionOffset = 0; - - // merge vertices - - for (auto genVertexColl : m_GenVertexCollections) { - debug() << " traverse collection of size " << genVertexColl->size() << endmsg; - for (const auto elem : *genVertexColl) { - collVOut->push_back(elem.clone()); - } - } - // merge particles, keeping the references up to date - for (unsigned int collCounter = 0; collCounter < m_MCParticleCollections.size(); ++collCounter) { - auto mcPartColl = m_MCParticleCollections[collCounter]; - for (const auto elem : *mcPartColl) { - auto newPart = edm4hep::MCParticle(elem.core()); - if (elem.startVertex().isAvailable()) { - // update reference: find new index in merged collection - // add offset since the signal particles are copied in new collection first - unsigned int newIdStart = elem.startVertex().getObjectID().index + collectionOffset; - // update startVertex - newPart.startVertex((*collVOut)[newIdStart]); - } - if (elem.endVertex().isAvailable()) { - // update reference: find new index in merged collection - // add offset since the signal particles are copied in new collection first - unsigned int newIdEnd = elem.endVertex().getObjectID().index + collectionOffset; - // update startVertex - newPart.endVertex((*collVOut)[newIdEnd]); - } - collPOut->push_back(newPart); - } - collectionOffset += m_GenVertexCollections[collCounter]->size(); - } - - m_vertOut.put(collVOut); - m_partOut.put(collPOut); - - m_MCParticleCollections.clear(); - m_GenVertexCollections.clear(); - return StatusCode::SUCCESS; -} diff --git a/JugBase/src/components/PileupParticlesMergeTool.h b/JugBase/src/components/PileupParticlesMergeTool.h deleted file mode 100644 index ef8d186b..00000000 --- a/JugBase/src/components/PileupParticlesMergeTool.h +++ /dev/null @@ -1,51 +0,0 @@ -#ifndef FWCORE_PILEUPPARTICLESMERGETOOL_H -#define FWCORE_PILEUPPARTICLESMERGETOOL_H - -// Gaudi -#include "GaudiAlg/GaudiTool.h" - -// FCCSW -#include "JugBase/DataHandle.h" -#include "JugBase/IEDMMergeTool.h" - -namespace fcc { -class MCParticleCollection; -class GenVertexCollection; -} - -/** @class PileupParticlesMergeTool - * - * Implementation of the EDMMergeTool interface used for pileup overlay - * of generated particles. Keeps track of the association between particles and vertices. - */ - -class PileupParticlesMergeTool : public GaudiTool, virtual public IEDMMergeTool { -public: - explicit PileupParticlesMergeTool(const std::string& aType, const std::string& aName, const IInterface* aParent); - virtual ~PileupParticlesMergeTool() = default; - - virtual StatusCode initialize() override final; - - virtual StatusCode finalize() override final; - - virtual StatusCode readPileupCollection(podio::EventStore& store) override final; - - virtual StatusCode mergeCollections() override final; - - virtual StatusCode readSignal() override final; - -private: - Gaudi::Property m_pileupGenParticlesBranchName{this, "genParticlesBranch", "genParticles", ""}; - Gaudi::Property m_pileupGenVerticesBranchName{this, "genVerticesBranch", "genVertices", ""}; - - std::vector m_MCParticleCollections; - std::vector m_GenVertexCollections; - - DataHandle m_vertOut{"overlay/allGenVertices", Gaudi::DataHandle::Writer, this}; - DataHandle m_partOut{"overlay/allGenParticles", Gaudi::DataHandle::Writer, this}; - - DataHandle m_vertIn{"overlay/signalGenVertices", Gaudi::DataHandle::Reader, this}; - DataHandle m_partIn{"overlay/signalGenParticles", Gaudi::DataHandle::Reader, this}; -}; - -#endif // FWCORE_PILEUPPARTICLESMERGETOOL_H -- GitLab From 87014a8c8fae9efb9f8686ce735ad4f0136c00cc Mon Sep 17 00:00:00 2001 From: Wouter Deconinck Date: Mon, 11 Jul 2022 19:25:23 +0200 Subject: [PATCH 3/3] FCCSW to JugBase --- JugBase/src/components/InputCopier.cpp | 2 +- JugTrack/src/components/TestACTSLogger.cpp | 2 +- JugTrack/src/components/TestACTSLogger.h | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/JugBase/src/components/InputCopier.cpp b/JugBase/src/components/InputCopier.cpp index 1b3e3e29..b53dc050 100644 --- a/JugBase/src/components/InputCopier.cpp +++ b/JugBase/src/components/InputCopier.cpp @@ -9,7 +9,7 @@ #include "GaudiAlg/GaudiTool.h" #include "Gaudi/Algorithm.h" -// FCCSW +// JugBase #include "JugBase/DataHandle.h" // Event Model related classes diff --git a/JugTrack/src/components/TestACTSLogger.cpp b/JugTrack/src/components/TestACTSLogger.cpp index c706b7d1..5ad60ae8 100644 --- a/JugTrack/src/components/TestACTSLogger.cpp +++ b/JugTrack/src/components/TestACTSLogger.cpp @@ -3,7 +3,7 @@ #include "TestACTSLogger.h" -// FCCSW +// JugBase //#include "DetCommon/DetUtils.h" #include "JugBase/IGeoSvc.h" diff --git a/JugTrack/src/components/TestACTSLogger.h b/JugTrack/src/components/TestACTSLogger.h index 31af7eb5..7413c393 100644 --- a/JugTrack/src/components/TestACTSLogger.h +++ b/JugTrack/src/components/TestACTSLogger.h @@ -8,7 +8,7 @@ #include "GaudiAlg/GaudiAlgorithm.h" #include "Gaudi/Property.h" -// FCCSW +// JugBase #include "JugBase/DataHandle.h" /** Logging for ACTS. -- GitLab