algorithms/calorimetry/EdepToNpeConversion.cc

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2026 Minho Kim
0003
0004 #include <DD4hep/Detector.h>
0005 #include <DD4hep/Readout.h>
0006 #include <DDSegmentation/BitFieldCoder.h>
0007 #include <fmt/format.h>
0008 #include <fmt/ranges.h>
0009 #include <cstddef>
0010 #include <fstream>
0011 #include <random>
0012 #include <sstream>
0013 #include <stdexcept>
0014 #include <string>
0015 #include <tuple>
0016 #include <utility>
0017 #include <vector>
0018
0019 #include "EdepToNpeConversion.h"
0020
0021 namespace eicrecon {
0022
0023 void EdepToNpeConversion::init() {
0024
0025   // A nonzero constant factor takes precedence if it is configured together with LUT
0026   if (m_cfg.edep_to_npe > 0) {
0027     if (!m_cfg.edep_to_npe_filename.empty()) {
0028       info("Constant edep-to-npe factor {} is set; ignoring LUT file {}", m_cfg.edep_to_npe,
0029            m_cfg.edep_to_npe_filename);
0030     }
0031     return;
0032   }
0033
0034   // Otherwise a complete LUT configuration is required
0035   if (m_cfg.readout.empty() || m_cfg.edep_to_npe_fields.empty() ||
0036       m_cfg.edep_to_npe_filename.empty()) {
0037     throw std::runtime_error(
0038         "No edep-to-npe conversion configured: set either a constant edep_to_npe factor or all "
0039         "of (readout, edep_to_npe_fields, edep_to_npe_filename)");
0040   }
0041
0042   // Get the cellID decoder and the indices of the LUT key fields
0043   try {
0044     m_id_spec = m_detector->readout(m_cfg.readout).idSpec();
0045   } catch (...) {
0046     throw std::runtime_error(fmt::format("Failed to get idSpec for readout {}", m_cfg.readout));
0047   }
0048   m_id_dec = m_id_spec.decoder();
0049   if (m_id_dec == nullptr) {
0050     throw std::runtime_error(fmt::format("Failed to get ID decoder for readout {}", m_cfg.readout));
0051   }
0052   for (const auto& field : m_cfg.edep_to_npe_fields) {
0053     try {
0054       m_field_idxs.push_back(m_id_dec->index(field));
0055     } catch (...) {
0056       throw std::runtime_error(
0057           fmt::format("Field {} not found in idSpec of readout {}", field, m_cfg.readout));
0058     }
0059   }
0060
0061   // Load the LUT and parses each line into a lookup key and a conversion factor
0062   std::string filename = m_cfg.edep_to_npe_filename;
0063   std::ifstream infile(filename);
0064   if (!infile) {
0065     throw std::runtime_error(fmt::format("Unable to open LUT file: {}", filename));
0066   }
0067   std::string line;
0068   std::size_t lineno = 0;
0069   while (std::getline(infile, line)) {
0070     lineno++;
0071     if (line.empty()) {
0072       throw std::runtime_error(
0073           fmt::format("Empty line in LUT file {} at line {}", filename, lineno));
0074     }
0075     std::istringstream iss(line);
0076     std::vector<int> key(m_cfg.edep_to_npe_fields.size());
0077     double factor;
0078     for (auto& value : key) {
0079       if (!(iss >> value)) {
0080         throw std::runtime_error(fmt::format("Malformed LUT file {} at line {}", filename, lineno));
0081       }
0082     }
0083     if (!(iss >> factor)) {
0084       throw std::runtime_error(fmt::format("Malformed LUT file {} at line {}", filename, lineno));
0085     }
0086     if (!m_edep_to_npe_lut.emplace(std::move(key), factor).second) {
0087       throw std::runtime_error(
0088           fmt::format("Duplicate key in LUT file {} at line {}", filename, lineno));
0089     }
0090   }
0091   if (m_edep_to_npe_lut.empty()) {
0092     throw std::runtime_error(fmt::format("LUT file {} contains no entries", filename));
0093   }
0094 } // EdepToNpeConversion:init
0095
0096 void EdepToNpeConversion::process(const EdepToNpeConversion::Input& input,
0097                                   const EdepToNpeConversion::Output& output) const {
0098   const auto [headers, inhits] = input;
0099   auto [outhits]               = output;
0100
0101   auto seed = m_uid.getUniqueID(*headers, name());
0102   std::mt19937 generator(seed);
0103
0104   for (const auto& hit : *inhits) {
0105     // Edep-to-Npe conversion & Apply Poisson smearing
0106     const double mean_npe = hit.getEnergy() * get_edep_to_npe_factor(hit);
0107     unsigned long npe     = 0;
0108     if (mean_npe > 0) {
0109       std::poisson_distribution<unsigned long> poisson(mean_npe);
0110       npe = poisson(generator);
0111     }
0112
0113     if (npe == 0) {
0114       continue;
0115     }
0116
0117     auto out_hit = hit.clone();
0118     out_hit.setEnergy(static_cast<float>(npe));
0119     outhits->push_back(out_hit);
0120   }
0121
0122 } // EdepToNpeConversion:process
0123
0124 double EdepToNpeConversion::get_edep_to_npe_factor(const edm4hep::SimCalorimeterHit& hit) const {
0125   if (m_cfg.edep_to_npe > 0) {
0126     return m_cfg.edep_to_npe;
0127   }
0128   std::vector<int> key;
0129   key.reserve(m_field_idxs.size());
0130   for (const auto idx : m_field_idxs) {
0131     key.push_back(static_cast<int>(m_id_dec->get(hit.getCellID(), idx)));
0132   }
0133   const auto it = m_edep_to_npe_lut.find(key);
0134   if (it == m_edep_to_npe_lut.end()) {
0135     throw std::runtime_error(
0136         fmt::format("No edep-to-npe factor for cellID {:#x} (fields [{}] = [{}])", hit.getCellID(),
0137                     fmt::join(m_cfg.edep_to_npe_fields, ", "), fmt::join(key, ", ")));
0138   }
0139   return it->second;
0140 }
0141
0142 } // namespace eicrecon