Stop looking for primaries for simple radiative photon-like geometries

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

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

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

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

src/algorithms/calorimetry/MCTools.h

@@ -0,0 +1,41 @@
+// 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 utiliy function in edm4xxx
+//
+
+#include <edm4hep/MCParticleCollection.h>
+#include <edm4hep/CaloHitContributionCollection.h>
+
+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 and has no childern -> daughter
+//  - All other cases (i.e. early showers, multi-radiation): primary
+// libraries
+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

src/algorithms/calorimetry/SimCalorimeterHitProcessor.cc

@@ -27,6 +27,7 @@
 #include <vector>
 
 #include "algorithms/calorimetry/SimCalorimeterHitProcessorConfig.h"
+#include "MCTools.h"
 
 using namespace dd4hep;
 
@@ -70,20 +71,6 @@ template <> struct hash<std::tuple<edm4hep::MCParticle, uint64_t, int>> {
 
 // 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) *