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 Chao Peng, Jihee Kim, Sylvester Joosten, Whitney Armstrong, Wouter Deconinck, David Lawrence, Derek Anderson
 
 /*
  *  An algorithm to group readout hits from a calorimeter
@@ -12,16 +12,17 @@
 
 #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 <edm4eic/CalorimeterHit.h>
 #include <Evaluator/DD4hepUnits.h>
 #include <Math/GenVector/Cartesian3D.h>
 #include <Math/GenVector/DisplacementVector3D.h>
 #include <fmt/core.h>
 #include <algorithm>
+#include <algorithms/service.h>
 #include <cmath>
 #include <cstddef>
 #include <gsl/pointers>
@@ -31,6 +32,7 @@
 #include <vector>
 
 #include "algorithms/calorimetry/CalorimeterHitsMergerConfig.h"
+#include "services/evaluator/EvaluatorSvc.h"
 
 namespace eicrecon {
 
@@ -41,24 +43,57 @@ void CalorimeterHitsMerger::init() {
         return;
     }
 
+    // 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 auto& field : m_cfg.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);
+
+    // if no field-by-field mappings provided, initialize relevant bitmasks
+    // otherwise intialize relevant functionals
+    if (m_cfg.mappings.empty()) {
+      id_mask = 0;
+      std::vector<RefField> 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_mask = ~id_mask;
+      debug("ID mask in {:s}: {:#064b}", m_cfg.readout, id_mask);
+    } else {
+
+      // lambda to translate IDDescriptor fields into function parameters
+      std::function hit_to_map = [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;
+      };
+
+      // intialize functions
+      //   - NOTE this assumes provided fields are 1-to-1!
+      auto& svc = algorithms::ServiceSvc::instance();
+      for (std::size_t iMap = 0; const auto& mapping : m_cfg.mappings) {
+        if (iMap < m_cfg.fields.size()) {
+          ref_maps[m_cfg.fields.at(iMap)] = svc.service<EvaluatorSvc>("EvaluatorSvc")->compile(mapping, hit_to_map);
+          trace("Mapping for field {} = {}", m_cfg.fields.at(iMap), mapping);
+        }
+        ++iMap;
+      }  // end loop over mappings
+    }
+
 }
 
 void CalorimeterHitsMerger::process(
@@ -69,14 +104,17 @@ 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) {
+    MergeMap merge_map;
+    if (m_cfg.mappings.empty()) {
+      for (std::size_t ix = 0; const auto &h : *in_hits) {
         uint64_t id = h.getCellID() & id_mask;
         merge_map[id].push_back(ix);
-
         ix++;
+      }
+    } else {
+      build_map_via_funcs(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 +178,44 @@ void CalorimeterHitsMerger::process(
     debug("Size before = {}, after = {}", in_hits->size(), out_hits->size());
 }
 
+void CalorimeterHitsMerger::build_map_via_funcs(
+  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) {
+
+    // make sure vector is clear
+    ref_fields.clear();
+    for (std::size_t iField = 0; const auto& name_field : id_desc.fields()) {
+
+      // check if field has associated mapping
+      const bool foundMapping = (
+        std::find(m_cfg.fields.begin(), m_cfg.fields.end(), name_field.first) != m_cfg.fields.end()
+      );
+
+      // if mapping provided for field, apply it
+      // otherwise just copy index
+      if (foundMapping) {
+        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;
+    }
+    const uint64_t ref_id = id_desc.encode(ref_fields);
+
+    // add hit to appropriate group
+    merge_map[ref_id].push_back(iHit);
+    ++iHit;
+
+  }  // end hit loop
+
+}  // end 'build_map_via_funcs(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 Chao Peng, Jihee Kim, Sylvester Joosten, Whitney Armstrong, Wouter Deconinck, David Lawrence, Derek Anderson
 
 /*
  *  An algorithm to group readout hits from a calorimeter
@@ -10,21 +10,32 @@
 
 #pragma once
 
+#include <DD4hep/BitFieldCoder.h>
 #include <DD4hep/Detector.h>
+#include <DD4hep/IDDescriptor.h>
 #include <DDRec/CellIDPositionConverter.h>
 #include <algorithms/algorithm.h>
 #include <algorithms/geo.h>
 #include <edm4eic/CalorimeterHitCollection.h>
 #include <stdint.h>
+#include <functional>
 #include <gsl/pointers>
+#include <map>
 #include <string>
 #include <string_view>
+#include <unordered_map>
+#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
@@ -51,10 +62,18 @@ namespace eicrecon {
   private:
     uint64_t id_mask{0}, 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_map_via_funcs(const edm4eic::CalorimeterHitCollection* in_hits, MergeMap& merge_map) const;
+
   };
 
 } // namespace eicrecon

src/algorithms/calorimetry/CalorimeterHitsMergerConfig.h

@@ -13,6 +13,7 @@ namespace eicrecon {
     std::string              readout{""};
     std::vector<std::string> fields{};
     std::vector<int>         refs{};
+    std::vector<std::string> mappings{};
 
   };
 

src/detectors/BHCAL/BHCAL.cc

@@ -12,10 +12,66 @@
 #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"
 
+namespace eicrecon {
+
+  // --------------------------------------------------------------------------
+  // Helper method to generate appropriate phi map
+  // --------------------------------------------------------------------------
+  std::string MakePhiMap(const std::size_t nMerge = 1) {
+
+    // convert num to merge to string
+    const std::string sNum = std::to_string(nMerge);
+
+    // construct mapping
+    std::string map  = "";
+    if (nMerge > 1) {
+      map = "phi-(" + sNum + "*((phi/" + sNum + ")-floor(phi/" + sNum + ")))";
+    } else {
+      map = "phi";
+    }
+    return map;
+
+  }  // end 'MakePhiMap(std::size_t)'
+
+  // --------------------------------------------------------------------------
+  // Helper method to generate appropriate adjacency matrix
+  // --------------------------------------------------------------------------
+  std::string MakeAdjacencyMatrix(const std::size_t nMerge = 1) {
+
+    // set up checks at wraparound.
+    // magic constants:
+    //   320 - number of tiles per row
+    const std::string sMerge    = std::to_string(nMerge);
+    const std::string wrapDef   = "(abs(phi_1 - phi_2) == (320 - 1))";
+    const std::string wrapMerge = "(abs(320 - abs(phi_1 - phi_2)) <= " + sMerge + ")";
+
+    // combine strings into horizontal adjacency conditions
+    const std::string phiAdjacent  = "(abs(phi_1 - phi_2) == " + sMerge + ")";
+    const std::string wrapAdjacent = nMerge > 1 ? wrapMerge : wrapDef;
+
+    // put everything together
+    std::string matrix("");
+    matrix += "(";
+    // check for vertically adjacent tiles
+    matrix += "  ( (abs(eta_1 - eta_2) == 1) && (abs(phi_1 - phi_2) == 0) ) ||";
+    // check for horizontally adjacent tiles
+    matrix += "  ( (abs(eta_1 - eta_2) == 0) && " + phiAdjacent + " ) ||";
+    // check for horizontally adjacent tiles at wraparound
+    matrix += "  ( (abs(eta_1 - eta_2) == 0) && " + wrapAdjacent + " )";
+    matrix += ") == 1";
+    return matrix;
+
+  }  // end 'MakeAdjacencyMatrix(std::size_t)'
+
+}
+
+
+
 extern "C" {
 
     void InitPlugin(JApplication *app) {
@@ -31,17 +87,15 @@ extern "C" {
         decltype(CalorimeterHitDigiConfig::pedSigmaADC)   HcalBarrel_pedSigmaADC = 2;
         decltype(CalorimeterHitDigiConfig::resolutionTDC) HcalBarrel_resolutionTDC = 1 * dd4hep::picosecond;
 
-        // Set default adjacency matrix. Magic constants:
-        //  320 - number of tiles per row
-        decltype(CalorimeterIslandClusterConfig::adjacencyMatrix) HcalBarrel_adjacencyMatrix =
-          "("
-          // check for vertically adjacent tiles
-          "  ( (abs(eta_1 - eta_2) == 1) && (abs(phi_1 - phi_2) == 0) ) ||"
-          // check for horizontally adjacent tiles
-          "  ( (abs(eta_1 - eta_2) == 0) && (abs(phi_1 - phi_2) == 1) ) ||"
-          // check for horizontally adjacent tiles at wraparound
-          "  ( (abs(eta_1 - eta_2) == 0) && (abs(phi_1 - phi_2) == (320 - 1)) )"
-          ") == 1";
+        // ---------------------------------------------------------------------
+        // if needed, merge adjacent tiles into a tower and adjust maps/matrices
+        // ---------------------------------------------------------------------
+        const std::size_t nPhiToMerge = 1;
+        const std::string phiMap      = MakePhiMap(nPhiToMerge);
+        const std::string adjMatrix   = MakeAdjacencyMatrix(nPhiToMerge);
+
+        // set matrix
+        decltype(CalorimeterIslandClusterConfig::adjacencyMatrix) HcalBarrel_adjacencyMatrix = adjMatrix;
 
         app->Add(new JOmniFactoryGeneratorT<CalorimeterHitDigi_factory>(
           "HcalBarrelRawHits",
@@ -66,6 +120,7 @@ extern "C" {
           },
           app   // TODO: Remove me once fixed
         ));
+
         app->Add(new JOmniFactoryGeneratorT<CalorimeterHitReco_factory>(
           "HcalBarrelRecHits", {"HcalBarrelRawHits"}, {"HcalBarrelRecHits"},
           {
@@ -83,12 +138,24 @@ extern "C" {
           },
           app   // TODO: Remove me once fixed
         ));
+
+        app->Add(new JOmniFactoryGeneratorT<CalorimeterHitsMerger_factory>(
+          "HcalBarrelMergedHits", {"HcalBarrelRecHits"}, {"HcalBarrelMergedHits"},
+          {
+            .readout = "HcalBarrelHits",
+            .fields = {"phi"},
+            .mappings = {phiMap}
+          },
+          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"},
+          "HcalBarrelIslandProtoClusters", {"HcalBarrelMergedHits"}, {"HcalBarrelIslandProtoClusters"},
           {
             .adjacencyMatrix = HcalBarrel_adjacencyMatrix,
             .readout = "HcalBarrelHits",
@@ -179,5 +246,6 @@ extern "C" {
             app   // TODO: Remove me once fixed
           )
         );
+
     }
 }

src/factories/calorimetry/CalorimeterHitsMerger_factory.h

@@ -22,6 +22,7 @@ class CalorimeterHitsMerger_factory : public JOmniFactory<CalorimeterHitsMerger_
     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_mappings {this, "mappings", config().mappings};
 
     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",