HEXPLIT: adjust timing cut to respect t0

src/algorithms/calorimetry/HEXPLIT.cc

@@ -11,7 +11,9 @@
 #include <Evaluator/DD4hepUnits.h>
 #include <Math/GenVector/Cartesian3D.h>
 #include <Math/GenVector/DisplacementVector3D.h>
+#include <edm4eic/MCRecoCalorimeterHitLinkCollection.h>
 #include <edm4hep/Vector3f.h>
+#include <podio/LinkNavigator.h>
 #include <algorithm>
 #include <cmath>
 #include <cstdlib>
@@ -99,19 +101,25 @@ void HEXPLIT::init() {
 
 void HEXPLIT::process(const HEXPLIT::Input& input, const HEXPLIT::Output& output) const {
 
-  const auto [hits]  = input;
-  auto [subcellHits] = output;
+  const auto [hits, mchitlinks] = input;
+  auto [subcellHits]            = output;
 
-  double MIP   = m_cfg.MIP / dd4hep::GeV;
-  double delta = m_cfg.delta_in_MIPs * MIP;
-  double Emin  = m_cfg.Emin_in_MIPs * MIP;
-  double tmax  = m_cfg.tmax / dd4hep::ns;
+  std::optional<podio::LinkNavigator<edm4eic::MCRecoCalorimeterHitLinkCollection>> nav;
+  if (mchitlinks != nullptr && !mchitlinks->empty()) {
+    nav.emplace(*mchitlinks);
+  }
+
+  double MIP    = m_cfg.MIP / dd4hep::GeV;
+  double delta  = m_cfg.delta_in_MIPs * MIP;
+  double Emin   = m_cfg.Emin_in_MIPs * MIP;
+  double max_dt = m_cfg.max_time_to_truth_t0 / dd4hep::ns;
 
   auto volman = m_detector->volumeManager();
 
   for (const auto& hit : *hits) {
     //skip hits that do not pass E and t cuts
-    if (hit.getEnergy() < Emin || hit.getTime() > tmax) {
+    const auto t0 = nav.has_value() ? get_t0(hit, *nav) : std::nullopt;
+    if (hit.getEnergy() < Emin || (t0.has_value() && (hit.getTime() - *t0) > max_dt)) {
       continue;
     }
 
@@ -130,7 +138,9 @@ void HEXPLIT::process(const HEXPLIT::Input& input, const HEXPLIT::Output& output
       if (dz > 2 || dz == 0) {
         continue;
       }
-      if (other_hit.getEnergy() < Emin || other_hit.getTime() > tmax) {
+      const auto other_t0 = nav.has_value() ? get_t0(other_hit, *nav) : std::nullopt;
+      if (other_hit.getEnergy() < Emin ||
+          (other_t0.has_value() && (other_hit.getTime() - *other_t0) > max_dt)) {
         continue;
       }
       //difference in transverse position (in units of side lengths)
@@ -221,4 +231,30 @@ void HEXPLIT::process(const HEXPLIT::Input& input, const HEXPLIT::Output& output
   }
 }
 
+edm4hep::MCParticle HEXPLIT::get_primary(const edm4hep::CaloHitContribution& contrib) {
+  return get_primary(contrib.getParticle());
+}
+
+edm4hep::MCParticle HEXPLIT::get_primary(const edm4hep::MCParticle& particle) {
+  edm4hep::MCParticle primary = particle;
+  while (primary.parents_size() > 0) {
+    if (primary.getGeneratorStatus() != 0) {
+      break;
+    }
+    primary = primary.getParents(0);
+  }
+  return primary;
+}
+
+std::optional<double>
+HEXPLIT::get_t0(const edm4eic::CalorimeterHit& hit,
+                const podio::LinkNavigator<edm4eic::MCRecoCalorimeterHitLinkCollection>& nav) {
+  for (const auto& [simhit, weight] : nav.getLinked(hit.getRawHit())) {
+    for (const auto& contrib : simhit.getContributions()) {
+      return get_primary(contrib).getTime();
+    }
+  }
+  return std::nullopt;
+}
+
 } // namespace eicrecon

src/algorithms/calorimetry/HEXPLIT.h

@@ -13,7 +13,12 @@
 #include <algorithms/algorithm.h>
 #include <algorithms/geo.h>
 #include <edm4eic/CalorimeterHitCollection.h>
+#include <edm4eic/MCRecoCalorimeterHitLinkCollection.h>
+#include <edm4hep/CaloHitContribution.h>
+#include <edm4hep/MCParticle.h>
 #include <gsl/pointers>
+#include <optional>
+#include <podio/LinkNavigator.h>
 #include <string>      // for basic_string
 #include <string_view> // for string_view
 #include <vector>
@@ -23,9 +28,10 @@
 
 namespace eicrecon {
 
-using HEXPLITAlgorithm =
-    algorithms::Algorithm<algorithms::Input<const edm4eic::CalorimeterHitCollection>,
-                          algorithms::Output<edm4eic::CalorimeterHitCollection>>;
+using HEXPLITAlgorithm = algorithms::Algorithm<
+    algorithms::Input<const edm4eic::CalorimeterHitCollection,
+                      std::optional<edm4eic::MCRecoCalorimeterHitLinkCollection>>,
+    algorithms::Output<edm4eic::CalorimeterHitCollection>>;
 
 class HEXPLIT : public HEXPLITAlgorithm, public WithPodConfig<HEXPLITConfig> {
 
@@ -102,6 +108,13 @@ class HEXPLIT : public HEXPLITAlgorithm, public WithPodConfig<HEXPLITConfig> {
 
   stagger_pattern stag = stag_H4;
 
+  static edm4hep::MCParticle get_primary(const edm4hep::CaloHitContribution& contrib);
+  static edm4hep::MCParticle get_primary(const edm4hep::MCParticle& particle);
+
+  static std::optional<double>
+  get_t0(const edm4eic::CalorimeterHit& hit,
+         const podio::LinkNavigator<edm4eic::MCRecoCalorimeterHitLinkCollection>& nav);
+
 private:
   const dd4hep::Detector* m_detector{algorithms::GeoSvc::instance().detector()};
 };

src/algorithms/calorimetry/HEXPLITConfig.h

@@ -9,7 +9,7 @@ struct HEXPLITConfig {
   double MIP{472. * dd4hep::keV};
   double Emin_in_MIPs{0.1};
   double delta_in_MIPs{0.01};
-  double tmax{325 * dd4hep::ns};
+  double max_time_to_truth_t0{325 * dd4hep::ns};
   enum StaggerType { H4 = 0, H3 = 1, S2 = 2 } stag_type = H4;
 };
 

src/detectors/FHCAL/FHCAL.cc

@@ -85,12 +85,12 @@ void InitPlugin(JApplication* app) {
       ));
 
   app->Add(new JOmniFactoryGeneratorT<HEXPLIT_factory>(
-      "HcalEndcapPInsertSubcellHits", {"HcalEndcapPInsertRecHits"},
+      "HcalEndcapPInsertSubcellHits", {"HcalEndcapPInsertRecHits", "HcalEndcapPInsertRawHitLinks"},
       {"HcalEndcapPInsertSubcellHits"},
       {
-          .MIP          = 480. * dd4hep::keV,
-          .Emin_in_MIPs = 0.5,
-          .tmax         = 162 * dd4hep::ns, //150 ns + (z at front face)/(speed of light)
+          .MIP                  = 480. * dd4hep::keV,
+          .Emin_in_MIPs         = 0.5,
+          .max_time_to_truth_t0 = 162 * dd4hep::ns, //150 ns + (z at front face)/(speed of light)
       },
       app // TODO: Remove me once fixed
       ));

src/detectors/ZDC/ZDC.cc

@@ -175,14 +175,14 @@ void InitPlugin(JApplication* app) {
       ));
 
   app->Add(new JOmniFactoryGeneratorT<HEXPLIT_factory>(
-      "HcalFarForwardZDCSubcellHits", {"HcalFarForwardZDCRecHits"},
+      "HcalFarForwardZDCSubcellHits", {"HcalFarForwardZDCRecHits", "HcalFarForwardZDCRawHitLinks"},
       {"HcalFarForwardZDCSubcellHits"},
       {
-          .MIP           = 630. * dd4hep::keV,
-          .Emin_in_MIPs  = 0.5,
-          .delta_in_MIPs = 0.01,
-          .tmax          = 269 * dd4hep::ns,
-          .stag_type     = HEXPLITConfig::StaggerType::S2,
+          .MIP                  = 630. * dd4hep::keV,
+          .Emin_in_MIPs         = 0.5,
+          .delta_in_MIPs        = 0.01,
+          .max_time_to_truth_t0 = 269 * dd4hep::ns, //150 ns + (z at front face)/(speed of light)
+          .stag_type            = HEXPLITConfig::StaggerType::S2,
       },
       app // TODO: Remove me once fixed
       ));

src/factories/calorimetry/HEXPLIT_factory.h

@@ -3,6 +3,8 @@
 
 #pragma once
 
+#include <edm4eic/MCRecoCalorimeterHitLinkCollection.h>
+
 #include "algorithms/calorimetry/HEXPLIT.h"
 #include "services/algorithms_init/AlgorithmsInit_service.h"
 #include "extensions/jana/JOmniFactory.h"
@@ -16,12 +18,14 @@ class HEXPLIT_factory : public JOmniFactory<HEXPLIT_factory, HEXPLITConfig> {
 private:
   std::unique_ptr<AlgoT> m_algo;
   PodioInput<edm4eic::CalorimeterHit> m_rec_hits_input{this};
+  PodioInput<edm4eic::MCRecoCalorimeterHitLink, true> m_mchitlinks_input{this};
   PodioOutput<edm4eic::CalorimeterHit> m_subcell_hits_output{this};
 
   ParameterRef<double> m_MIP{this, "MIP", config().MIP};
   ParameterRef<double> m_Emin_in_MIPs{this, "Emin_in_MIPs", config().Emin_in_MIPs};
   ParameterRef<double> m_delta_in_MIPs{this, "delta_in_MIPs", config().delta_in_MIPs};
-  ParameterRef<double> m_tmax{this, "tmax", config().tmax};
+  ParameterRef<double> m_max_time_to_truth_t0{this, "max_time_to_truth_t0",
+                                              config().max_time_to_truth_t0};
 
   Service<AlgorithmsInit_service> m_algorithmsInit{this};
 
@@ -34,7 +38,7 @@ class HEXPLIT_factory : public JOmniFactory<HEXPLIT_factory, HEXPLITConfig> {
   }
 
   void Process(int32_t /* run_number */, uint64_t /* event_number */) {
-    m_algo->process({m_rec_hits_input()}, {m_subcell_hits_output().get()});
+    m_algo->process({m_rec_hits_input(), m_mchitlinks_input()}, {m_subcell_hits_output().get()});
   }
 };
 

src/tests/algorithms_test/calorimetry_HEXPLIT.cc

@@ -35,8 +35,8 @@ TEST_CASE("the subcell-splitting algorithm runs", "[HEXPLIT]") {
   logger->set_level(spdlog::level::trace);
 
   HEXPLITConfig cfg;
-  cfg.MIP  = 472. * dd4hep::keV;
-  cfg.tmax = 1000. * dd4hep::ns;
+  cfg.MIP                  = 472. * dd4hep::keV;
+  cfg.max_time_to_truth_t0 = 1000. * dd4hep::ns;
 
   auto detector = algorithms::GeoSvc::instance().detector();
   auto id_desc  = detector->readout("MockCalorimeterHits").idSpec();
@@ -79,7 +79,7 @@ TEST_CASE("the subcell-splitting algorithm runs", "[HEXPLIT]") {
   }
 
   auto subcellhits_coll = std::make_unique<edm4eic::CalorimeterHitCollection>();
-  algo.process({&hits_coll}, {subcellhits_coll.get()});
+  algo.process({&hits_coll, nullptr}, {subcellhits_coll.get()});
 
   //the number of subcell hits should be equal to the
   //number of subcells per cell (12) times the number of cells (5)
@@ -101,3 +101,55 @@ TEST_CASE("the subcell-splitting algorithm runs", "[HEXPLIT]") {
   // is in the subcell where the other hits overlap
   REQUIRE((*subcellhits_coll)[35].getEnergy() / E[2] > 0.95);
 }
+
+struct HEXPLITFixture {
+  HEXPLIT algo{"HEXPLIT"};
+  double side_length   = 31.3 * dd4hep::mm;
+  double layer_spacing = 25.1 * dd4hep::mm;
+  double thickness     = 3 * dd4hep::mm;
+  edm4hep::Vector3f dimension;
+  uint64_t cellID;
+
+  HEXPLITFixture() {
+    std::shared_ptr<spdlog::logger> logger = spdlog::default_logger()->clone("HEXPLIT");
+    logger->set_level(spdlog::level::warn);
+
+    HEXPLITConfig cfg;
+    cfg.MIP                  = 472. * dd4hep::keV;
+    cfg.max_time_to_truth_t0 = 500. * dd4hep::ns;
+    algo.applyConfig(cfg);
+    algo.init();
+
+    dimension = edm4hep::Vector3f(2 * side_length, std::numbers::sqrt3 * side_length, thickness);
+    auto id_desc =
+        algorithms::GeoSvc::instance().detector()->readout("MockCalorimeterHits").idSpec();
+    cellID = id_desc.encode({{"system", 255}, {"x", 0}, {"y", 0}});
+  }
+
+  // 5 hits in consecutive layers at the given time, same positions as the main test.
+  edm4eic::CalorimeterHitCollection make_hits(float time_ns) {
+    std::array<double, 5> layer = {0, 1, 2, 3, 4};
+    std::array<double, 5> x     = {0, 0.75 * side_length, 0, 0.75 * side_length, 0};
+    std::array<double, 5> y     = {
+        std::numbers::sqrt3 / 2 * side_length, -0.25 * std::numbers::sqrt3 * side_length, 0,
+        0.25 * std::numbers::sqrt3 * side_length, std::numbers::sqrt3 / 2 * side_length};
+
+    edm4eic::CalorimeterHitCollection hits;
+    for (std::size_t i = 0; i < 5; i++) {
+      hits.create(cellID, 50 * dd4hep::MeV, 0.f, time_ns, 0.f,
+                  edm4hep::Vector3f(x[i], y[i], layer[i] * layer_spacing), dimension, 0, layer[i],
+                  edm4hep::Vector3f(x[i], y[i], layer[i] * layer_spacing));
+    }
+    return hits;
+  }
+};
+
+TEST_CASE("HEXPLIT timing cut: no MCParticles skips the cut", "[HEXPLIT]") {
+  HEXPLITFixture f;
+  // Hits are 1000 ns late — well outside the 500 ns window.
+  // Without a link collection the cut is skipped entirely.
+  auto hits = f.make_hits(1000.f);
+  auto out  = std::make_unique<edm4eic::CalorimeterHitCollection>();
+  f.algo.process({&hits, nullptr}, {out.get()});
+  REQUIRE(out->size() == 60);
+}