algorithms/digi/PulseCombiner.cc

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2025-2026 Simon Gardner, Minho Kim
0003 //
0004 // Combine pulses into a larger pulse if they are within a certain time of each other
0005
0006 #include <DD4hep/Detector.h>
0007 #include <DD4hep/IDDescriptor.h>
0008 #include <DD4hep/Readout.h>
0009 #include <DDSegmentation/BitFieldCoder.h>
0010 #include <algorithms/geo.h>
0011 #include <edm4hep/MCParticle.h>
0012 #include <edm4hep/SimCalorimeterHit.h>
0013 #include <edm4hep/SimTrackerHit.h>
0014 #include <edm4hep/Vector3f.h>
0015 #include <podio/RelationRange.h>
0016 #include <algorithm>
0017 #include <cmath>
0018 #include <cstddef>
0019 #include <gsl/pointers>
0020 #include <map>
0021 #include <numeric>
0022 #include <stdexcept>
0023 #include <tuple>
0024 #include <utility>
0025 #include <vector>
0026
0027 #include "PulseCombiner.h"
0028
0029 namespace eicrecon {
0030
0031 void PulseCombiner::init() {
0032
0033   // Get the detector readout and set CellID bit mask if set
0034   if (!(m_cfg.readout.empty() && m_cfg.combine_field.empty())) {
0035     try {
0036       auto detector      = algorithms::GeoSvc::instance().detector();
0037       auto id_spec       = detector->readout(m_cfg.readout).idSpec();
0038       m_detector_bitmask = 0;
0039
0040       for (const auto& field : id_spec.fields()) {
0041         // Get the field name
0042         std::string field_name = field.first;
0043         // Check if the field is the one we want to combine
0044         m_detector_bitmask |= id_spec.field(field_name)->mask();
0045         if (field_name == m_cfg.combine_field) {
0046           break;
0047         }
0048       }
0049
0050     } catch (...) {
0051       error("Failed set bitshift for detector {} with segmentation id {}", m_cfg.readout,
0052             m_cfg.combine_field);
0053       throw std::runtime_error("Failed to load ID decoder");
0054     }
0055   }
0056 }
0057
0058 void PulseCombiner::process(const PulseCombiner::Input& input,
0059                             const PulseCombiner::Output& output) const {
0060   const auto [inPulses] = input;
0061   auto [outPulses]      = output;
0062
0063   // Create map containing vector of pulses from each CellID
0064   std::map<uint64_t, std::vector<PulseType>> cell_pulses;
0065   for (const PulseType& pulse : *inPulses) {
0066     uint64_t shiftedCellID = pulse.getCellID() & m_detector_bitmask;
0067     cell_pulses[shiftedCellID].push_back(pulse);
0068   }
0069
0070   // Loop over detector elements and combine pulses
0071   for (const auto& [cellID, pulses] : cell_pulses) {
0072     if (pulses.size() == 1) {
0073       outPulses->push_back(pulses.at(0).clone());
0074       debug("CellID {} has only one pulse, no combination needed", cellID);
0075     } else {
0076       // Order the pulses by time and group those that are close in time into clusters
0077       std::vector<std::vector<PulseType>> clusters = clusterPulses(pulses);
0078       for (const auto& cluster : clusters) {
0079         // Clone the first pulse in the cluster
0080         auto sum_pulse = outPulses->create();
0081         sum_pulse.setCellID(cluster[0].getCellID());
0082         sum_pulse.setInterval(cluster[0].getInterval());
0083         sum_pulse.setTime(cluster[0].getTime());
0084
0085         // Sum the amplitudes of the pulses in the cluster.
0086         // The pulses must be time-ordered, which clusterPulses() has already done
0087         auto newPulse = sumTimeOrderedPulses(cluster);
0088         for (auto pulse : newPulse) {
0089           sum_pulse.addToAmplitude(pulse);
0090         }
0091
0092         // Sum the pulse array
0093         float integral = std::accumulate(newPulse.begin(), newPulse.end(), 0.0F);
0094         sum_pulse.setIntegral(integral);
0095         sum_pulse.setPosition(edm4hep::Vector3f(
0096             cluster[0].getPosition().x, cluster[0].getPosition().y, cluster[0].getPosition().z));
0097         for (const auto& pulse : cluster) {
0098           sum_pulse.addToPulses(pulse);
0099           for (auto particle : pulse.getParticles()) {
0100             sum_pulse.addToParticles(particle);
0101           }
0102           for (auto hit : pulse.getTrackerHits()) {
0103             sum_pulse.addToTrackerHits(hit);
0104           }
0105           for (auto hit : pulse.getCalorimeterHits()) {
0106             sum_pulse.addToCalorimeterHits(hit);
0107           }
0108         }
0109       }
0110       debug("CellID {} has {} pulses, combined into {} clusters", cellID, pulses.size(),
0111             clusters.size());
0112     }
0113   }
0114
0115 } // PulseCombiner:process
0116
0117 std::vector<std::vector<PulseType>>
0118 PulseCombiner::clusterPulses(const std::vector<PulseType>& pulses) const {
0119
0120   // Copied so they can be sorted
0121   std::vector<PulseType> ordered_pulses{pulses};
0122
0123   // Sort pulses by time, greaty simplifying the combination process
0124   std::ranges::sort(ordered_pulses, [](const PulseType& a, const PulseType& b) {
0125     return a.getTime() < b.getTime();
0126   });
0127
0128   // Create vector of pulses
0129   std::vector<std::vector<PulseType>> cluster_pulses;
0130   float clusterEndTime = 0;
0131   bool makeNewPulse    = true;
0132   // Create clusters of pulse indices which overlap with at least the minimum separation
0133   for (const auto& pulse : ordered_pulses) {
0134     float pulseStartTime = pulse.getTime();
0135     float pulseEndTime   = pulse.getTime() + pulse.getInterval() * pulse.getAmplitude().size();
0136     if (!makeNewPulse) {
0137       if (pulseStartTime < clusterEndTime + m_cfg.minimum_separation) {
0138         cluster_pulses.back().push_back(pulse);
0139         clusterEndTime = std::max(clusterEndTime, pulseEndTime);
0140       } else {
0141         makeNewPulse = true;
0142       }
0143     }
0144     if (makeNewPulse) {
0145       cluster_pulses.push_back({pulse});
0146       clusterEndTime = pulseEndTime;
0147       makeNewPulse   = false;
0148     }
0149   }
0150
0151   return cluster_pulses;
0152
0153 } // PulseCombiner::clusterPulses
0154
0155 std::vector<float> PulseCombiner::sumTimeOrderedPulses(const std::vector<PulseType>& pulses) {
0156
0157   const float startTime = pulses[0].getTime();
0158   const float interval  = pulses[0].getInterval();
0159
0160   std::vector<float> newPulse;
0161
0162   for (const auto& pulse : pulses) {
0163     const auto startStep =
0164         static_cast<std::size_t>(std::round((pulse.getTime() - startTime) / interval));
0165     const auto& amplitude = pulse.getAmplitude();
0166
0167     // Extend the combined pulse so that it can hold this pulse's contribution
0168     newPulse.resize(std::max(newPulse.size(), startStep + amplitude.size()), 0.0);
0169
0170     for (std::size_t i = 0; i < amplitude.size(); ++i) {
0171       newPulse[startStep + i] += amplitude[i];
0172     }
0173   }
0174
0175   return newPulse;
0176 } // PulseCombiner::sumTimeOrderedPulses
0177
0178 } // namespace eicrecon