Enable calorimeter hit merging by functions (#1668)

### Briefly, what does this PR introduce? This PR extends the functionality of the `CalorimeterHitsMerger` algorithm. Previously, reconstructed hits could only be merged across a given field of the readout of a calorimeter. This presents a challenge for calorimeters such as the Barrel HCal where 1. our readout has no segmentation beyond just eta & phi, and 2. it could be useful to study the response of the detector as a function of how many readout channels we gang together after reconstruction. This PR addresses this point by utilizing the `EvaluatorSvc` in a manner similar to the `adjacencyMatrix`, `peakNeighbourhoodMatrix` of the `CalorimeterIslandClustering` and the `sampFrac` of `CalorimeterHitReco`. Now the user has the ability to specify an (almost) arbitrarily complex transformation for a specific field of the readout via the `fieldTransformations` parameter, which defines both the field to transform and the function to map the indices of that field onto the desired reference indices. For example: ``` app->Add(new JOmniFactoryGeneratorT<CalorimeterHitsMerger_factory>( "HcalBarrelMergedHits", {"HcalBarrelRecHits"}, {"HcalBarrelMergedHits"}, { .readout = "HcalBarrelHits", .fieldTransformations = {"phi:phi-(5*((phi/5)-floor(phi/5)))"} }, app // TODO: Remove me once fixed )); ``` Here, the `HcalBarrelMergedHits` collection will merge 5 hits (i.e. scintillator tiles for the BHCal) adjacent in phi into a one with the position and cellID of the 1st of the 5, and no hits will be merged along eta. The previous behavior of the algorithm can be recovered by simply specifying the index to be mapped onto. For example: ``` app->Add(new JOmniFactoryGeneratorT<CalorimeterHitsMerger_factory>( "HcalEndcapNMergedHits", {"HcalEndcapNRecHits"}, {"HcalEndcapNMergedHits"}, { .readout = "HcalEndcapNHits", .fieldTransformations = {"layer:4", "slice:0"} }, app // TODO: Remove me once fixed )); ``` An example script of to change a transformation (and how to update the adjacency matrix accordingly) from the command-line is provided in the snippets repo [here](https://github.com/eic/snippets/blob/main/Calorimetery/CaloDebugTools/UtilityScripts/RunEICReconWithTileMerging.rb). ### What kind of change does this PR introduce? - [ ] Bug fix (issue #__) - [x] New feature (issue #1669 ) - [ ] Documentation update - [ ] Other: __ ### Please check if this PR fulfills the following: - [ ] Tests for the changes have been added - [ ] Documentation has been added / updated - [x] Changes have been communicated to collaborators ### Does this PR introduce breaking changes? What changes might users need to make to their code? No. ### Does this PR change default behavior? No. --------- Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> Co-authored-by: Dmitry Kalinkin <[email protected]>
eic · Jan 9, 2025 · d96cbff · d96cbff
1 parent 8fa634b
commit d96cbff
Show file tree

Hide file tree

Showing 8 changed files with 139 additions and 34 deletions.
diff --git a/src/algorithms/calorimetry/CalorimeterHitsMerger.cc b/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
diff --git a/src/algorithms/calorimetry/CalorimeterHitsMerger.h b/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
diff --git a/src/algorithms/calorimetry/CalorimeterHitsMergerConfig.h b/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{};
 
   };
 

diff --git a/src/detectors/BHCAL/BHCAL.cc b/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
           )
         );
+
     }
 }
diff --git a/src/detectors/EHCAL/EHCAL.cc b/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
         ));

diff --git a/src/detectors/FHCAL/FHCAL.cc b/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
         ));

diff --git a/src/factories/calorimetry/CalorimeterHitsMerger_factory.h b/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};
 

diff --git a/src/services/io/podio/JEventProcessorPODIO.cc b/src/services/io/podio/JEventProcessorPODIO.cc
@@ -290,6 +290,7 @@ JEventProcessorPODIO::JEventProcessorPODIO() {
             "LFHCALSplitMergeClusterAssociations",
             "HcalBarrelRawHits",
             "HcalBarrelRecHits",
+            "HcalBarrelMergedHits",
             "HcalBarrelClusters",
             "HcalBarrelClusterAssociations",
             "HcalBarrelSplitMergeClusters",