HEXPLIT: adjust timing cut to respect t0

src/algorithms/calorimetry/HEXPLIT.cc

@@ -1,4 +1,4 @@
-// Copyright (C) 2023 Sebouh Paul
+// Copyright (C) 2023 - 2026 Sebouh Paul, ePIC Collaboration
 // SPDX-License-Identifier: LGPL-3.0-or-later
 
 // References:
@@ -11,7 +11,11 @@
 #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 <podio/ObjectID.h>
+#include <unordered_map>
 #include <algorithm>
 #include <cmath>
 #include <cstdlib>
@@ -99,19 +103,27 @@ 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;
+  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;
+
+  // Per-event t0 cache: populated lazily by get_t0(), keyed on hit ObjectID.
+  T0Cache t0_cache;
+  std::optional<podio::LinkNavigator<edm4eic::MCRecoCalorimeterHitLinkCollection>> nav;
+  if (mchitlinks != nullptr && !mchitlinks->empty()) {
+    nav.emplace(*mchitlinks);
+  }
 
   auto volman = m_detector->volumeManager();
 
   for (const auto& hit : *hits) {
+    const auto t0 = nav ? get_t0(hit, *nav, t0_cache) : std::nullopt;
     //skip hits that do not pass E and t cuts
-    if (hit.getEnergy() < Emin || hit.getTime() > tmax) {
+    if (hit.getEnergy() < Emin || (t0 && (hit.getTime() - *t0) > max_dt)) {
       continue;
     }
 
@@ -130,7 +142,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 ? get_t0(other_hit, *nav, t0_cache) : std::nullopt;
+      if (other_hit.getEnergy() < Emin ||
+          (other_t0 && (other_hit.getTime() - *other_t0) > max_dt)) {
         continue;
       }
       //difference in transverse position (in units of side lengths)
@@ -221,4 +235,39 @@ 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,
+                T0Cache& cache) {
+  const auto id = hit.getObjectID();
+  if (auto it = cache.find(id); it != cache.end()) {
+    return it->second;
+  }
+  std::optional<double> t0;
+  for (const auto& [simhit, weight] : nav.getLinked(hit.getRawHit())) {
+    for (const auto& contrib : simhit.getContributions()) {
+      t0 = get_primary(contrib).getTime();
+      break;
+    }
+    break;
+  }
+  cache.emplace(id, t0);
+  return t0;
+}
+
 } // namespace eicrecon

src/algorithms/calorimetry/HEXPLIT.h

@@ -13,7 +13,14 @@
 #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 <podio/ObjectID.h>
+#include <unordered_map>
 #include <string>      // for basic_string
 #include <string_view> // for string_view
 #include <vector>
@@ -23,9 +30,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 +110,16 @@ 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);
+
+  using T0Cache = std::unordered_map<podio::ObjectID, std::optional<double>>;
+
+  static std::optional<double>
+  get_t0(const edm4eic::CalorimeterHit& hit,
+         const podio::LinkNavigator<edm4eic::MCRecoCalorimeterHitLinkCollection>& nav,
+         T0Cache& cache);
+
 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()});
   }
 };