Enable calorimeter hit merging by functions

src/algorithms/calorimetry/CalorimeterHitsMerger.cc

@@ -1,5 +1,5 @@
 // SPDX-License-Identifier: LGPL-3.0-or-later
-// Copyright (C) 2022 Chao Peng, Jihee Kim, Sylvester Joosten, Whitney Armstrong, Wouter Deconinck, David Lawrence
+// Copyright (C) 2022 - 2025 Chao Peng, Jihee Kim, Sylvester Joosten, Whitney Armstrong, Wouter Deconinck, David Lawrence, Derek Anderson
 
 /*
  *  An algorithm to group readout hits from a calorimeter
@@ -12,14 +12,14 @@
 
 #include <DD4hep/Alignments.h>
 #include <DD4hep/DetElement.h>
-#include <DD4hep/IDDescriptor.h>
 #include <DD4hep/Objects.h>
 #include <DD4hep/Readout.h>
 #include <DD4hep/VolumeManager.h>
 #include <DDSegmentation/BitFieldCoder.h>
 #include <Evaluator/DD4hepUnits.h>
 #include <Math/GenVector/Cartesian3D.h>
 #include <Math/GenVector/DisplacementVector3D.h>
+#include <algorithms/service.h>
 #include <fmt/core.h>
 #include <algorithm>
 #include <cmath>
@@ -31,6 +31,7 @@
 #include <vector>
 
 #include "algorithms/calorimetry/CalorimeterHitsMergerConfig.h"
+#include "services/evaluator/EvaluatorSvc.h"
 
 namespace eicrecon {
 
@@ -41,24 +42,64 @@ void CalorimeterHitsMerger::init() {
         return;
     }
 
+    // split parameters into vectors of fields
+    // and of transformations
+    std::vector<std::string> fields;
+    std::vector<std::string> transforms;
+    for (const std::string& field_transform : m_cfg.fieldTransformations) {
+
+      const std::size_t isplit = field_transform.find_first_of(':');
+      if (isplit == std::string::npos) {
+        warning("transform '{}' ill-formatted. Format is <field>:<transformation>.", field_transform);
+      }
+
+
+      fields.emplace_back(
+        field_transform.substr(0, isplit)
+      );
+      transforms.emplace_back(
+        field_transform.substr(isplit + 1)
+      );
+    }
+
+
+    // initialize descriptor + decoders
     try {
-        auto id_desc = m_detector->readout(m_cfg.readout).idSpec();
-        id_mask = 0;
-        std::vector<std::pair<std::string, int>> ref_fields;
-        for (size_t i = 0; i < m_cfg.fields.size(); ++i) {
-            id_mask |= id_desc.field(m_cfg.fields[i])->mask();
-            // use the provided id number to find ref cell, or use 0
-            int ref = i < m_cfg.refs.size() ? m_cfg.refs[i] : 0;
-            ref_fields.emplace_back(m_cfg.fields[i], ref);
-        }
-        ref_mask = id_desc.encode(ref_fields);
+      id_desc = m_detector->readout(m_cfg.readout).idSpec();
+      id_decoder = id_desc.decoder();
+      for (const std::string& field : fields) {
+        const short index = id_decoder->index(field);
+      }
     } catch (...) {
-        auto mess = fmt::format("Failed to load ID decoder for {}", m_cfg.readout);
-        warning(mess);
+      auto mess = fmt::format("Failed to load ID decoder for {}", m_cfg.readout);
+      warning(mess);
 //        throw std::runtime_error(mess);
     }
-    id_mask = ~id_mask;
-    debug("ID mask in {:s}: {:#064b}", m_cfg.readout, id_mask);
+
+    // lambda to translate IDDescriptor fields into function parameters
+    std::function hit_transform = [this](const edm4eic::CalorimeterHit& hit) {
+      std::unordered_map<std::string, double> params;
+      for (const auto& name_field : id_desc.fields()) {
+        params.emplace(name_field.first, name_field.second->value(hit.getCellID()));
+        trace("{} = {}", name_field.first, name_field.second->value(hit.getCellID()));
+      }
+      return params;
+    };
+
+    // loop through provided readout fields
+    auto& svc = algorithms::ServiceSvc::instance();
+    for (std::size_t iField = 0; std::string& field : fields) {
+
+      // grab provided transformation and field
+      const std::string field_transform = transforms.at(iField);
+      auto name_field = id_desc.field(field);
+
+      // set transformation for each field
+      ref_maps[field] = svc.service<EvaluatorSvc>("EvaluatorSvc")->compile(field_transform, hit_transform);
+      trace("{}: using transformation '{}'", field, field_transform);
+      ++iField;
+    }  // end field loop
+
 }
 
 void CalorimeterHitsMerger::process(
@@ -69,14 +110,9 @@ void CalorimeterHitsMerger::process(
     auto [out_hits] = output;
 
     // find the hits that belong to the same group (for merging)
-    std::unordered_map<uint64_t, std::vector<std::size_t>> merge_map;
-    std::size_t ix = 0;
-    for (const auto &h : *in_hits) {
-        uint64_t id = h.getCellID() & id_mask;
-        merge_map[id].push_back(ix);
-
-        ix++;
-    }
+    MergeMap merge_map;
+    build_merge_map(in_hits, merge_map);
+    debug("Merge map built: merging {} hits into {} merged hits", in_hits->size(), merge_map.size());
 
     // sort hits by energy from large to small
     for (auto &it : merge_map) {
@@ -140,4 +176,37 @@ void CalorimeterHitsMerger::process(
     debug("Size before = {}, after = {}", in_hits->size(), out_hits->size());
 }
 
+void CalorimeterHitsMerger::build_merge_map(
+      const edm4eic::CalorimeterHitCollection* in_hits,
+      MergeMap& merge_map) const {
+
+  std::vector<RefField> ref_fields;
+  for (std::size_t iHit = 0; const auto& hit : *in_hits) {
+
+    ref_fields.clear();
+    for (std::size_t iField = 0; const auto& name_field : id_desc.fields()) {
+
+      // apply mapping to field if provided,
+      // otherwise copy value of field
+      if (ref_maps.count(name_field.first) > 0) {
+        ref_fields.push_back(
+          {name_field.first, ref_maps[name_field.first](hit)}
+        );
+      } else {
+        ref_fields.push_back(
+          {name_field.first, id_decoder->get(hit.getCellID(), name_field.first)}
+        );
+      }
+      ++iField;
+    }
+
+    // encode new cell ID and add hit to map
+    const uint64_t ref_id = id_desc.encode(ref_fields);
+    merge_map[ref_id].push_back(iHit);
+    ++iHit;
+
+  }  // end hit loop
+
+}  // end 'build_merge_map(edm4eic::CalorimeterHitsCollection*, MergeMap&)'
+
 } // namespace eicrecon

src/algorithms/calorimetry/CalorimeterHitsMerger.h

@@ -1,5 +1,5 @@
 // SPDX-License-Identifier: LGPL-3.0-or-later
-// Copyright (C) 2022 Chao Peng, Jihee Kim, Sylvester Joosten, Whitney Armstrong, Wouter Deconinck, David Lawrence
+// Copyright (C) 2022 - 2025 Chao Peng, Jihee Kim, Sylvester Joosten, Whitney Armstrong, Wouter Deconinck, David Lawrence, Derek Anderson
 
 /*
  *  An algorithm to group readout hits from a calorimeter
@@ -11,20 +11,33 @@
 #pragma once
 
 #include <DD4hep/Detector.h>
+#include <DD4hep/IDDescriptor.h>
 #include <DDRec/CellIDPositionConverter.h>
+#include <Parsers/Primitives.h>
 #include <algorithms/algorithm.h>
 #include <algorithms/geo.h>
 #include <edm4eic/CalorimeterHitCollection.h>
 #include <stdint.h>
+#include <cstddef>
+#include <functional>
 #include <gsl/pointers>
+#include <map>
 #include <string>
 #include <string_view>
+#include <unordered_map>
+#include <utility>
+#include <vector>
 
 #include "CalorimeterHitsMergerConfig.h"
 #include "algorithms/interfaces/WithPodConfig.h"
 
 namespace eicrecon {
 
+  // aliases for convenience
+  using MergeMap = std::unordered_map<uint64_t, std::vector<std::size_t>>;
+  using RefField = std::pair<std::string, int>;
+  using MapFunc = std::function<int(const edm4eic::CalorimeterHit&)>;
+
   using CalorimeterHitsMergerAlgorithm = algorithms::Algorithm<
     algorithms::Input<
       edm4eic::CalorimeterHitCollection
@@ -49,12 +62,20 @@ namespace eicrecon {
     void process(const Input&, const Output&) const final;
 
   private:
-    uint64_t id_mask{0}, ref_mask{0};
+    uint64_t ref_mask{0};
+
+  private:
+    mutable std::map<std::string, MapFunc> ref_maps;
+    dd4hep::IDDescriptor id_desc;
+    dd4hep::BitFieldCoder* id_decoder;
 
   private:
     const dd4hep::Detector* m_detector{algorithms::GeoSvc::instance().detector()};
     const dd4hep::rec::CellIDPositionConverter* m_converter{algorithms::GeoSvc::instance().cellIDPositionConverter()};
 
+  private:
+    void build_merge_map(const edm4eic::CalorimeterHitCollection* in_hits, MergeMap& merge_map) const;
+
   };
 
 } // namespace eicrecon

src/algorithms/calorimetry/CalorimeterHitsMergerConfig.h

@@ -11,8 +11,7 @@ namespace eicrecon {
   struct CalorimeterHitsMergerConfig {
 
     std::string              readout{""};
-    std::vector<std::string> fields{};
-    std::vector<int>         refs{};
+    std::vector<std::string> fieldTransformations{};
 
   };
 

src/detectors/BHCAL/BHCAL.cc

@@ -12,6 +12,7 @@
 #include "factories/calorimetry/CalorimeterClusterRecoCoG_factory.h"
 #include "factories/calorimetry/CalorimeterHitDigi_factory.h"
 #include "factories/calorimetry/CalorimeterHitReco_factory.h"
+#include "factories/calorimetry/CalorimeterHitsMerger_factory.h"
 #include "factories/calorimetry/CalorimeterIslandCluster_factory.h"
 #include "factories/calorimetry/CalorimeterTruthClustering_factory.h"
 #include "factories/calorimetry/TrackClusterMergeSplitter_factory.h"
@@ -66,6 +67,7 @@ extern "C" {
           },
           app   // TODO: Remove me once fixed
         ));
+
         app->Add(new JOmniFactoryGeneratorT<CalorimeterHitReco_factory>(
           "HcalBarrelRecHits", {"HcalBarrelRawHits"}, {"HcalBarrelRecHits"},
           {
@@ -83,10 +85,25 @@ extern "C" {
           },
           app   // TODO: Remove me once fixed
         ));
+
+        // --------------------------------------------------------------------
+        // If needed, merge adjacent phi tiles into towers. By default,
+        // NO merging will be done. This can be changed at runtime.
+        // --------------------------------------------------------------------
+        app->Add(new JOmniFactoryGeneratorT<CalorimeterHitsMerger_factory>(
+          "HcalBarrelMergedHits", {"HcalBarrelRecHits"}, {"HcalBarrelMergedHits"},
+          {
+            .readout = "HcalBarrelHits",
+            .fieldTransformations = {"phi:phi"}
+          },
+          app   // TODO: Remove me once fixed
+        ));
+
         app->Add(new JOmniFactoryGeneratorT<CalorimeterTruthClustering_factory>(
           "HcalBarrelTruthProtoClusters", {"HcalBarrelRecHits", "HcalBarrelHits"}, {"HcalBarrelTruthProtoClusters"},
           app   // TODO: Remove me once fixed
         ));
+
         app->Add(new JOmniFactoryGeneratorT<CalorimeterIslandCluster_factory>(
           "HcalBarrelIslandProtoClusters", {"HcalBarrelRecHits"}, {"HcalBarrelIslandProtoClusters"},
           {
@@ -179,5 +196,6 @@ extern "C" {
             app   // TODO: Remove me once fixed
           )
         );
+
     }
 }

src/detectors/EHCAL/EHCAL.cc

@@ -71,8 +71,7 @@ extern "C" {
           "HcalEndcapNMergedHits", {"HcalEndcapNRecHits"}, {"HcalEndcapNMergedHits"},
           {
             .readout = "HcalEndcapNHits",
-            .fields = {"layer", "slice"},
-            .refs = {4, 0}, // place merged hits at ~1 interaction length deep
+            .fieldTransformations = {"layer:4", "slice:0"}
           },
           app   // TODO: Remove me once fixed
         ));

src/detectors/FHCAL/FHCAL.cc

@@ -83,8 +83,7 @@ extern "C" {
           "HcalEndcapPInsertMergedHits", {"HcalEndcapPInsertRecHits"}, {"HcalEndcapPInsertMergedHits"},
           {
             .readout = "HcalEndcapPInsertHits",
-            .fields = {"layer", "slice"},
-            .refs = {1, 0},
+            .fieldTransformations = {"layer:1", "slice:0"}
           },
           app   // TODO: Remove me once fixed
         ));

src/factories/calorimetry/CalorimeterHitsMerger_factory.h

@@ -20,8 +20,7 @@ class CalorimeterHitsMerger_factory : public JOmniFactory<CalorimeterHitsMerger_
     PodioOutput<edm4eic::CalorimeterHit> m_hits_output {this};
 
     ParameterRef<std::string> m_readout {this, "readout", config().readout};
-    ParameterRef<std::vector<std::string>> m_fields {this, "fields", config().fields};
-    ParameterRef<std::vector<int>> m_refs {this, "refs", config().refs};
+    ParameterRef<std::vector<std::string>> m_field_transformations {this, "fieldTransformations", config().fieldTransformations};
 
     Service<AlgorithmsInit_service> m_algorithmsInit {this};
 

src/services/io/podio/JEventProcessorPODIO.cc

@@ -290,6 +290,7 @@ JEventProcessorPODIO::JEventProcessorPODIO() {
             "LFHCALSplitMergeClusterAssociations",
             "HcalBarrelRawHits",
             "HcalBarrelRecHits",
+            "HcalBarrelMergedHits",
             "HcalBarrelClusters",
             "HcalBarrelClusterAssociations",
             "HcalBarrelSplitMergeClusters",