Add CalorimeterCALOROCCalibration

src/algorithms/calorimetry/CalorimeterCALOROCReco.h

@@ -0,0 +1,84 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2026 Chun Yuen Tsang, Minho Kim
+
+#pragma once
+
+#include <DD4hep/DetElement.h>
+#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 <edm4eic/MCRecoCalorimeterHitAssociationCollection.h>
+#include <edm4eic/MCRecoCalorimeterHitLinkCollection.h>
+#include <edm4eic/RawCALOROCHitCollection.h>
+#include <edm4hep/RawCalorimeterHitCollection.h>
+#include <edm4hep/SimCalorimeterHitCollection.h>
+#include <cstddef>
+#include <functional>
+#include <gsl/pointers>
+#include <map>
+#include <string>
+#include <string_view>
+#include <vector>
+
+#include "CalorimeterCALOROCRecoConfig.h"
+#include "algorithms/interfaces/WithPodConfig.h"
+
+namespace eicrecon {
+
+using CalorimeterCALOROCRecoAlgorithm = algorithms::Algorithm<
+    algorithms::Input<edm4hep::SimCalorimeterHitCollection, edm4eic::RawCALOROCHitCollection,
+                      edm4hep::SimCalorimeterHitCollection, edm4eic::RawCALOROCHitCollection>,
+    algorithms::Output<edm4eic::CalorimeterHitCollection, edm4hep::RawCalorimeterHitCollection,
+                       edm4eic::MCRecoCalorimeterHitLinkCollection,
+                       edm4eic::MCRecoCalorimeterHitAssociationCollection>>;
+
+class CalorimeterCALOROCReco : public CalorimeterCALOROCRecoAlgorithm,
+                               public WithPodConfig<CalorimeterCALOROCRecoConfig> {
+
+public:
+  CalorimeterCALOROCReco(std::string_view name)
+      : CalorimeterCALOROCRecoAlgorithm{
+            name,
+            {"inputNpeHitPCollection", "inputADCPCollection", "inputNpeHitNCollection",
+             "inputADCNCollection"},
+            {"outputRecHitCollection", "outputRawHitCollection", "outputRawLink", "outputRawAssoc"},
+            "Reconstruct hit from CALOROC ADC data with MC truth from NpeHits."} {}
+
+  void init() final;
+  void process(const Input&, const Output&) const final;
+
+private:
+  std::function<double(const edm4hep::RawCalorimeterHit& h)> sampFrac;
+
+  dd4hep::IDDescriptor id_spec;
+  dd4hep::BitFieldCoder* id_dec = nullptr;
+
+  double m_reference_z_p, m_reference_z_n;
+  std::size_t sector_idx{0}, layer_idx{0};
+
+  mutable bool warned_unsupported_segmentation = false;
+
+  dd4hep::DetElement m_local;
+  std::size_t local_mask = ~static_cast<std::size_t>(0), gpos_mask = static_cast<std::size_t>(0);
+
+  std::map<std::vector<int>, double> m_edep_to_npe_lut{};
+  std::vector<std::size_t> m_field_idxs{};
+
+  double _energyCor(double referencePos, double energy, double z) const;
+  double _sumADC(const edm4eic::RawCALOROCHit& ADC) const;
+  double _toa(const edm4eic::RawCALOROCHit& ADC) const;
+  double _timeWalkCorrection(double toa, double ADC) const;
+
+  const algorithms::GeoSvc& m_geo = algorithms::GeoSvc::instance();
+
+private:
+  const dd4hep::Detector* m_detector{algorithms::GeoSvc::instance().detector()};
+  const dd4hep::rec::CellIDPositionConverter* m_converter{
+      algorithms::GeoSvc::instance().cellIDPositionConverter()};
+};
+
+} // namespace eicrecon

src/algorithms/calorimetry/CalorimeterCALOROCRecoConfig.h

@@ -0,0 +1,86 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2026 Chun Yuen Tsang, Minho Kim
+
+#pragma once
+
+#include <cstdint>
+#include <iostream>
+#include <string>
+#include <vector>
+
+namespace eicrecon {
+
+struct CalorimeterCALOROCRecoConfig {
+  // readout fields
+  std::string readout{""};
+  std::string layerField{""};
+  std::string sectorField{""};
+
+  // name of detelment or fields to find the local detector (for global->local transform)
+  // if nothing is provided, the lowest level DetElement (from cellID) will be used
+  std::string localDetElement{""};
+  std::vector<std::string> localDetFields{};
+  std::string maskPos{""};
+  std::vector<std::string> maskPosFields{};
+
+  std::string edep_to_npe_filename{""};
+  std::vector<std::string> edep_to_npe_fields{};
+  enum class ProxyType { sum = 0, simpson = 1, templateFit = 2 } proxy_type = ProxyType::sum;
+
+  bool timeWalkCor{false};
+  bool useNpeHitPos{false};
+
+  uint16_t highGainDR{0}; // high gain dynamic range
+  double gainRatio{1};    // gain ratio of low gain to high gain, it should be < 1
+
+  std::string attenuationReferencePositionNamePos{""};
+  std::string attenuationReferencePositionNameNeg{""};
+
+  // parameters for attenuation function
+  // [0] * exp(-|z_ref - z| / [1]) + (1 - [0]) * exp(-|z_ref - z| / [2])
+  // specified in edm4eic::units where dimensionfull
+  std::vector<double> attenuationParameters{0};
+
+  std::vector<double> timeWalkCorrectionParameters{0};
+  std::vector<double> lightSpeedParameters{0};
+
+  // calibration parameters
+  double slope     = 1 / 9.41e-2;
+  double intercept = 0;
+};
+
+std::istream& operator>>(std::istream& in, CalorimeterCALOROCRecoConfig::ProxyType& proxyType) {
+  std::string s;
+  in >> s;
+  // stringifying the enums causes them to be converted to integers before conversion to strings
+  if (s == "sum" or s == "0") {
+    proxyType = CalorimeterCALOROCRecoConfig::ProxyType::sum;
+  } else if (s == "simpson" or s == "1") {
+    proxyType = CalorimeterCALOROCRecoConfig::ProxyType::simpson;
+  } else if (s == "templateFit" or s == "2") {
+    proxyType = CalorimeterCALOROCRecoConfig::ProxyType::templateFit;
+  } else {
+    in.setstate(std::ios::failbit); // Set the fail bit if the input is not valid
+  }
+
+  return in;
+}
+std::ostream& operator<<(std::ostream& out,
+                         const CalorimeterCALOROCRecoConfig::ProxyType& proxyType) {
+  switch (proxyType) {
+  case CalorimeterCALOROCRecoConfig::ProxyType::sum:
+    out << "sum";
+    break;
+  case CalorimeterCALOROCRecoConfig::ProxyType::simpson:
+    out << "simpson";
+    break;
+  case CalorimeterCALOROCRecoConfig::ProxyType::templateFit:
+    out << "templateFit";
+    break;
+  default:
+    out.setstate(std::ios::failbit);
+  }
+  return out;
+}
+
+} // namespace eicrecon

src/algorithms/calorimetry/EdepToNpeConversion.cc

@@ -0,0 +1,142 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2026 Minho Kim
+
+#include <DD4hep/Detector.h>
+#include <DD4hep/Readout.h>
+#include <DDSegmentation/BitFieldCoder.h>
+#include <fmt/format.h>
+#include <fmt/ranges.h>
+#include <cstddef>
+#include <fstream>
+#include <random>
+#include <sstream>
+#include <stdexcept>
+#include <string>
+#include <tuple>
+#include <utility>
+#include <vector>
+
+#include "EdepToNpeConversion.h"
+
+namespace eicrecon {
+
+void EdepToNpeConversion::init() {
+
+  // A nonzero constant factor takes precedence if it is configured together with LUT
+  if (m_cfg.edep_to_npe > 0) {
+    if (!m_cfg.edep_to_npe_filename.empty()) {
+      info("Constant edep-to-npe factor {} is set; ignoring LUT file {}", m_cfg.edep_to_npe,
+           m_cfg.edep_to_npe_filename);
+    }
+    return;
+  }
+
+  // Otherwise a complete LUT configuration is required
+  if (m_cfg.readout.empty() || m_cfg.edep_to_npe_fields.empty() ||
+      m_cfg.edep_to_npe_filename.empty()) {
+    throw std::runtime_error(
+        "No edep-to-npe conversion configured: set either a constant edep_to_npe factor or all "
+        "of (readout, edep_to_npe_fields, edep_to_npe_filename)");
+  }
+
+  // Get the cellID decoder and the indices of the LUT key fields
+  try {
+    m_id_spec = m_detector->readout(m_cfg.readout).idSpec();
+  } catch (...) {
+    throw std::runtime_error(fmt::format("Failed to get idSpec for readout {}", m_cfg.readout));
+  }
+  m_id_dec = m_id_spec.decoder();
+  if (m_id_dec == nullptr) {
+    throw std::runtime_error(fmt::format("Failed to get ID decoder for readout {}", m_cfg.readout));
+  }
+  for (const auto& field : m_cfg.edep_to_npe_fields) {
+    try {
+      m_field_idxs.push_back(m_id_dec->index(field));
+    } catch (...) {
+      throw std::runtime_error(
+          fmt::format("Field {} not found in idSpec of readout {}", field, m_cfg.readout));
+    }
+  }
+
+  // Load the LUT and parses each line into a lookup key and a conversion factor
+  std::string filename = m_cfg.edep_to_npe_filename;
+  std::ifstream infile(filename);
+  if (!infile) {
+    throw std::runtime_error(fmt::format("Unable to open LUT file: {}", filename));
+  }
+  std::string line;
+  std::size_t lineno = 0;
+  while (std::getline(infile, line)) {
+    lineno++;
+    if (line.empty()) {
+      throw std::runtime_error(
+          fmt::format("Empty line in LUT file {} at line {}", filename, lineno));
+    }
+    std::istringstream iss(line);
+    std::vector<int> key(m_cfg.edep_to_npe_fields.size());
+    double factor;
+    for (auto& value : key) {
+      if (!(iss >> value)) {
+        throw std::runtime_error(fmt::format("Malformed LUT file {} at line {}", filename, lineno));
+      }
+    }
+    if (!(iss >> factor)) {
+      throw std::runtime_error(fmt::format("Malformed LUT file {} at line {}", filename, lineno));
+    }
+    if (!m_edep_to_npe_lut.emplace(std::move(key), factor).second) {
+      throw std::runtime_error(
+          fmt::format("Duplicate key in LUT file {} at line {}", filename, lineno));
+    }
+  }
+  if (m_edep_to_npe_lut.empty()) {
+    throw std::runtime_error(fmt::format("LUT file {} contains no entries", filename));
+  }
+} // EdepToNpeConversion:init
+
+void EdepToNpeConversion::process(const EdepToNpeConversion::Input& input,
+                                  const EdepToNpeConversion::Output& output) const {
+  const auto [headers, inhits] = input;
+  auto [outhits]               = output;
+
+  auto seed = m_uid.getUniqueID(*headers, name());
+  std::mt19937 generator(seed);
+
+  for (const auto& hit : *inhits) {
+    // Edep-to-Npe conversion & Apply Poisson smearing
+    const double mean_npe = hit.getEnergy() * get_edep_to_npe_factor(hit);
+    long npe              = 0;
+    if (mean_npe > 0) {
+      std::poisson_distribution<long> poisson(mean_npe);
+      npe = poisson(generator);
+    }
+
+    if (npe == 0) {
+      continue;
+    }
+
+    auto out_hit = hit.clone();
+    out_hit.setEnergy(static_cast<float>(npe));
+    outhits->push_back(out_hit);
+  }
+
+} // EdepToNpeConversion:process
+
+double EdepToNpeConversion::get_edep_to_npe_factor(const edm4hep::SimCalorimeterHit& hit) const {
+  if (m_cfg.edep_to_npe > 0) {
+    return m_cfg.edep_to_npe;
+  }
+  std::vector<int> key;
+  key.reserve(m_field_idxs.size());
+  for (const auto idx : m_field_idxs) {
+    key.push_back(static_cast<int>(m_id_dec->get(hit.getCellID(), idx)));
+  }
+  const auto it = m_edep_to_npe_lut.find(key);
+  if (it == m_edep_to_npe_lut.end()) {
+    throw std::runtime_error(
+        fmt::format("No edep-to-npe factor for cellID {:#x} (fields [{}] = [{}])", hit.getCellID(),
+                    fmt::join(m_cfg.edep_to_npe_fields, ", "), fmt::join(key, ", ")));
+  }
+  return it->second;
+}
+
+} // namespace eicrecon

src/algorithms/calorimetry/EdepToNpeConversion.h

@@ -0,0 +1,51 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2026 Minho Kim
+
+#pragma once
+
+#include <DD4hep/Detector.h>
+#include <DD4hep/IDDescriptor.h>
+#include <Parsers/Primitives.h>
+#include <algorithms/algorithm.h>
+#include <algorithms/geo.h>
+#include <edm4hep/EventHeaderCollection.h>
+#include <edm4hep/SimCalorimeterHitCollection.h>
+#include <cstddef>
+#include <gsl/pointers>
+#include <map>
+#include <string>
+#include <string_view>
+#include <vector>
+
+#include "algorithms/calorimetry/EdepToNpeConversionConfig.h"
+#include "algorithms/interfaces/UniqueIDGenSvc.h"
+#include "algorithms/interfaces/WithPodConfig.h"
+
+namespace eicrecon {
+
+using EdepToNpeConversionAlgorithm = algorithms::Algorithm<
+    algorithms::Input<edm4hep::EventHeaderCollection, edm4hep::SimCalorimeterHitCollection>,
+    algorithms::Output<edm4hep::SimCalorimeterHitCollection>>;
+
+class EdepToNpeConversion : public EdepToNpeConversionAlgorithm,
+                            public WithPodConfig<EdepToNpeConversionConfig> {
+
+public:
+  EdepToNpeConversion(std::string_view name)
+      : EdepToNpeConversionAlgorithm{name, {"EventHeader", "inputHits"}, {"outputHits"}, {}} {}
+  void init() final;
+  void process(const Input&, const Output&) const final;
+
+private:
+  double get_edep_to_npe_factor(const edm4hep::SimCalorimeterHit& hit) const;
+
+  const dd4hep::Detector* m_detector{algorithms::GeoSvc::instance().detector()};
+  const algorithms::UniqueIDGenSvc& m_uid = algorithms::UniqueIDGenSvc::instance();
+
+  dd4hep::IDDescriptor m_id_spec;
+  const dd4hep::BitFieldCoder* m_id_dec = nullptr;
+  std::vector<std::size_t> m_field_idxs{};
+  std::map<std::vector<int>, double> m_edep_to_npe_lut{};
+};
+
+} // namespace eicrecon

src/algorithms/calorimetry/EdepToNpeConversionConfig.h

@@ -0,0 +1,20 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2026 Minho Kim
+
+#pragma once
+
+#include <string>
+#include <vector>
+
+namespace eicrecon {
+
+struct EdepToNpeConversionConfig {
+  std::string readout{};
+  std::vector<std::string>
+      edep_to_npe_fields{}; // Fields the edep-to-npe conversion factor depends on
+  std::string
+      edep_to_npe_filename{}; // Lookup table name to get the field-dependent edep-to-npe conversion factors
+  double edep_to_npe{}; // Constant edep-to-npe conversion factor
+};
+
+} // namespace eicrecon