algorithms/digi/PulseCombiner.cc

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