diff --git a/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.cc b/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.cc index 5afc5288f7..ba7d7a0ac2 100644 --- a/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.cc +++ b/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.cc @@ -31,6 +31,7 @@ #include "CalorimeterClusterRecoCoG.h" #include "algorithms/calorimetry/CalorimeterClusterRecoCoGConfig.h" +#include "MCTools.h" namespace eicrecon { @@ -219,7 +220,7 @@ void CalorimeterClusterRecoCoG::associate( // -------------------------------------------------------------------- // grab primary responsible for contribution & increment relevant sum // -------------------------------------------------------------------- - edm4hep::MCParticle primary = get_primary(contrib); + edm4hep::MCParticle primary = MCTools::lookup_primary(contrib); mapMCParToContrib[primary] += contrib.getEnergy(); trace("Identified primary: id = {}, pid = {}, total energy = {}, contributed = {}", @@ -251,22 +252,4 @@ void CalorimeterClusterRecoCoG::associate( } } -edm4hep::MCParticle -CalorimeterClusterRecoCoG::get_primary(const edm4hep::CaloHitContribution& contrib) { - // get contributing particle - const auto contributor = contrib.getParticle(); - - // walk back through parents to find primary - // - TODO finalize primary selection. This - // can be improved!! - edm4hep::MCParticle primary = contributor; - while (primary.parents_size() > 0) { - if (primary.getGeneratorStatus() != 0) { - break; - } - primary = primary.getParents(0); - } - return primary; -} - } // namespace eicrecon diff --git a/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.h b/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.h index ebadcdd685..0c8562909c 100644 --- a/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.h +++ b/src/algorithms/calorimetry/CalorimeterClusterRecoCoG.h @@ -81,7 +81,6 @@ class CalorimeterClusterRecoCoG : public CalorimeterClusterRecoCoGAlgorithm, void associate(const edm4eic::Cluster& cl, const edm4eic::MCRecoCalorimeterHitAssociationCollection* mchitassociations, edm4eic::MCRecoClusterParticleAssociationCollection* assocs) const; - static edm4hep::MCParticle get_primary(const edm4hep::CaloHitContribution& contrib); }; } // namespace eicrecon diff --git a/src/algorithms/calorimetry/ImagingClusterReco.cc b/src/algorithms/calorimetry/ImagingClusterReco.cc index 3cd3379c3d..0f296e8f9f 100644 --- a/src/algorithms/calorimetry/ImagingClusterReco.cc +++ b/src/algorithms/calorimetry/ImagingClusterReco.cc @@ -30,6 +30,7 @@ #include "algorithms/calorimetry/ClusterTypes.h" #include "algorithms/calorimetry/ImagingClusterRecoConfig.h" +#include "MCTools.h" namespace eicrecon { @@ -277,7 +278,7 @@ void ImagingClusterReco::associate_mc_particles( // -------------------------------------------------------------------- // grab primary responsible for contribution & increment relevant sum // -------------------------------------------------------------------- - edm4hep::MCParticle primary = get_primary(contrib); + edm4hep::MCParticle primary = MCTools::lookup_primary(contrib); mapMCParToContrib[primary] += contrib.getEnergy(); trace("Identified primary: id = {}, pid = {}, total energy = {}, contributed = {}", @@ -309,22 +310,4 @@ void ImagingClusterReco::associate_mc_particles( } } -edm4hep::MCParticle -ImagingClusterReco::get_primary(const edm4hep::CaloHitContribution& contrib) const { - // get contributing particle - const auto contributor = contrib.getParticle(); - - // walk back through parents to find primary - // - TODO finalize primary selection. This - // can be improved!! - edm4hep::MCParticle primary = contributor; - while (primary.parents_size() > 0) { - if (primary.getGeneratorStatus() != 0) { - break; - } - primary = primary.getParents(0); - } - return primary; -} - } // namespace eicrecon diff --git a/src/algorithms/calorimetry/ImagingClusterReco.h b/src/algorithms/calorimetry/ImagingClusterReco.h index d5d0d349a4..f7d6e3cf42 100644 --- a/src/algorithms/calorimetry/ImagingClusterReco.h +++ b/src/algorithms/calorimetry/ImagingClusterReco.h @@ -76,8 +76,6 @@ class ImagingClusterReco : public ImagingClusterRecoAlgorithm, const edm4eic::Cluster& cl, const edm4eic::MCRecoCalorimeterHitAssociationCollection* mchitassociations, edm4eic::MCRecoClusterParticleAssociationCollection* assocs) const; - - edm4hep::MCParticle get_primary(const edm4hep::CaloHitContribution& contrib) const; }; } // namespace eicrecon diff --git a/src/algorithms/calorimetry/MCTools.h b/src/algorithms/calorimetry/MCTools.h new file mode 100644 index 0000000000..70c92b49a8 --- /dev/null +++ b/src/algorithms/calorimetry/MCTools.h @@ -0,0 +1,40 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2025 Minho Kim, Sylvester Joosten, Derek Anderson, Wouter Deconinck +#pragma once + +// +// @TODO should be migrated to a shared utility function in edm4xxx +// + +#include +#include + +namespace eicrecon::MCTools { + +// Lookup primary MCParticle +// we stop looking if we find the parent has status 1 +// so we don't merge e.g. radiative photons with the primary electron as this prevents us +// from properly linking the clusters back to the event geometry. Note that we also +// enforce for this case that no steps back towards the primary were taken to avoid +// storing the first pair of calorimetric showers that start inside the tracking volume. +// Hence, this algorithm will return: +// - Contribution came from primary: primary +// - Contribution came from immediate daughter of primary which has no children -> daughter +// - All other cases (i.e. early showers, multi-radiation): primary +inline edm4hep::MCParticle lookup_primary(const edm4hep::CaloHitContribution& contrib) { + const auto contributor = contrib.getParticle(); + + edm4hep::MCParticle primary = contributor; + size_t steps_taken = 0; // The number of steps taken looking for the primary + while (primary.parents_size() > 0) { + auto parent = primary.getParents(0); + if (primary.getGeneratorStatus() != 0 || + (parent.getGeneratorStatus() != 0 && steps_taken == 0)) { + break; + } + primary = parent; + steps_taken += 1; + } + return primary; +} +} // namespace eicrecon::MCTools diff --git a/src/algorithms/calorimetry/SimCalorimeterHitProcessor.cc b/src/algorithms/calorimetry/SimCalorimeterHitProcessor.cc index ab27984b24..2c1804f126 100644 --- a/src/algorithms/calorimetry/SimCalorimeterHitProcessor.cc +++ b/src/algorithms/calorimetry/SimCalorimeterHitProcessor.cc @@ -27,6 +27,7 @@ #include #include "algorithms/calorimetry/SimCalorimeterHitProcessorConfig.h" +#include "MCTools.h" using namespace dd4hep; @@ -70,20 +71,6 @@ template <> struct hash> { // unnamed namespace for internal utility namespace { -// Lookup primary MCParticle @TODO this should be a shared utiliy function in the edm4xxx -// libraries -edm4hep::MCParticle lookup_primary(const edm4hep::CaloHitContribution& contrib) { - const auto contributor = contrib.getParticle(); - - edm4hep::MCParticle primary = contributor; - while (primary.parents_size() > 0) { - if (primary.getGeneratorStatus() != 0) { - break; - } - primary = primary.getParents(0); - } - return primary; -} class HitContributionAccumulator { private: float m_energy{0}; @@ -195,7 +182,7 @@ void SimCalorimeterHitProcessor::process(const SimCalorimeterHitProcessor::Input m_attenuationReferencePosition ? get_attenuation(ih.getPosition().z) : 1.; // Use primary particle (traced back through parents) to group contributions for (const auto& contrib : ih.getContributions()) { - edm4hep::MCParticle primary = lookup_primary(contrib); + edm4hep::MCParticle primary = MCTools::lookup_primary(contrib); const double propagationTime = m_attenuationReferencePosition ? std::abs(m_attenuationReferencePosition.value() - ih.getPosition().z) *