Combine pulses in the same CellID

src/algorithms/digi/PulseCombiner.cc

@@ -0,0 +1,146 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2025 Simon Gardner
+//
+// Combine pulses into a larger pulse if they are within a certain time of each other
+
+#include <unordered_map>
+#include <vector>
+
+#include "PulseCombiner.h"
+#include <algorithms/geo.h>
+
+
+namespace eicrecon {
+
+void PulseCombiner::init() {
+
+  // Get the detector readout and set CellID bit mask if set
+  if(!(m_cfg.readout.empty() && m_cfg.combine_field.empty())) {
+    try {
+      auto detector = algorithms::GeoSvc::instance().detector();
+      auto id_spec = detector->readout(m_cfg.readout).idSpec();
+      auto id_dec = id_spec.decoder();
+      m_detector_bitmask = 0;
+
+      for(auto & field : id_spec.fields()) {
+        // Get the field name
+        std::string field_name = field.first;
+        // Check if the field is the one we want to combine
+        m_detector_bitmask |= id_spec.field(field_name)->mask();
+        if(field_name == m_cfg.combine_field) {
+          break;
+        }
+      }
+
+    } catch (...) {
+      error("Failed set bitshift for detector {} with segmentation id {}", m_cfg.readout, m_cfg.combine_field);
+      throw std::runtime_error("Failed to load ID decoder");
+    }
+  }
+
+}
+
+void PulseCombiner::process(const PulseCombiner::Input& input,
+                            const PulseCombiner::Output& output) const {
+  const auto [inPulses] = input;
+  auto [outPulses]      = output;
+
+  // Create map containing vector of pulses from each CellID
+  std::map<uint64_t, std::vector<edm4hep::TimeSeries>> cell_pulses;
+  for (const edm4hep::TimeSeries& pulse: *inPulses) {
+    uint64_t shiftedCellID = pulse.getCellID() & m_detector_bitmask;
+    cell_pulses[shiftedCellID].push_back(pulse);
+  }
+
+  // Loop over detector elements and combine pulses
+  for (const auto& [cellID, pulses] : cell_pulses) {
+    if(pulses.size() == 1) {
+      outPulses->push_back(pulses.at(0).clone());
+      debug("CellID {} has only one pulse, no combination needed", cellID);
+    }
+    else{
+      std::vector<std::vector<edm4hep::TimeSeries>> clusters = clusterPulses(pulses);
+      for (auto cluster: clusters) {
+        // Clone the first pulse in the cluster
+        auto sum_pulse = outPulses->create();
+        sum_pulse.setCellID(cluster[0].getCellID());
+        sum_pulse.setInterval(cluster[0].getInterval());
+        sum_pulse.setTime(cluster[0].getTime());
+
+        auto newPulse = sumPulses(cluster);
+        for(auto pulse : newPulse) {
+          sum_pulse.addToAmplitude(pulse);
+        }
+      }
+      debug("CellID {} has {} pulses, combined into {} clusters", cellID, pulses.size(), clusters.size());
+    }
+
+  }
+
+} // PulseCombiner:process
+
+std::vector<std::vector<edm4hep::TimeSeries>> PulseCombiner::clusterPulses(const std::vector<edm4hep::TimeSeries> pulses) const {
+
+  // Clone the pulses array of pointers so they aren't const
+  std::vector<edm4hep::TimeSeries> ordered_pulses {pulses};
+
+  // Sort pulses by time, greaty simplifying the combination process
+  std::sort(ordered_pulses.begin(), ordered_pulses.end(), [](edm4hep::TimeSeries a, edm4hep::TimeSeries b) {
+    return a.getTime() < b.getTime();
+  });
+
+  // Create vector of pulses
+  std::vector<std::vector<edm4hep::TimeSeries>> cluster_pulses;
+  float clusterEndTime = 0;
+  bool  makeNewPulse = true;
+  // Create clusters of pulse indices which overlap with at least the minimum separation
+  for (auto pulse: ordered_pulses) {
+    float pulseStartTime = pulse.getTime();
+    float pulseEndTime = pulse.getTime() + pulse.getInterval()*pulse.getAmplitude().size();
+    if(!makeNewPulse) {
+      if (pulseStartTime < clusterEndTime + m_cfg.minimum_separation) {
+        cluster_pulses.back().push_back(pulse);
+        clusterEndTime = std::max(clusterEndTime, pulseEndTime);
+      } else {
+        makeNewPulse = true;
+      }
+    }
+    if(makeNewPulse) {
+      cluster_pulses.push_back({pulse});
+      clusterEndTime = pulseEndTime;
+      makeNewPulse = false;
+    }
+  }
+
+  return cluster_pulses;
+
+} // PulseCombiner::clusterPulses
+
+std::vector<float> PulseCombiner::sumPulses(const std::vector<edm4hep::TimeSeries> pulses) const {
+
+  // Find maximum time of pulses in cluster
+  float maxTime = 0;
+  for (auto pulse: pulses) {
+    maxTime = std::max(maxTime, pulse.getTime() + pulse.getInterval()*pulse.getAmplitude().size());
+  }
+
+  //Calculate maxTime in interval bins
+  int maxStep = std::round((maxTime - pulses[0].getTime()) / pulses[0].getInterval());
+
+  std::vector<float> newPulse(maxStep, 0.0);
+
+  for(auto pulse: pulses) {
+    //Calculate start and end of pulse in interval bins
+    int startStep = (pulse.getTime() - pulses[0].getTime())/pulse.getInterval();
+    int pulseSize = pulse.getAmplitude().size();
+    int endStep = startStep + pulseSize;
+    for(int i = startStep; i < endStep; i++) {
+      // Add pulse values to new pulse
+      newPulse[i] += pulse.getAmplitude()[i - startStep];
+    }
+  }
+
+  return newPulse;
+} // PulseCombiner::sumPulses
+
+} // namespace eicrecon

src/algorithms/digi/PulseCombiner.h

@@ -0,0 +1,42 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2025 Simon Gardner
+//
+// Convert energy deposition into ADC pulses
+// ADC pulses are assumed to follow the shape of landau function
+
+#pragma once
+
+#include <algorithms/algorithm.h>
+#include <edm4hep/TimeSeriesCollection.h>
+#include <memory>
+#include <string>
+#include <string_view>
+#include <edm4eic/unit_system.h>
+
+#include "algorithms/digi/PulseCombinerConfig.h"
+#include "algorithms/interfaces/WithPodConfig.h"
+
+namespace eicrecon {
+
+using PulseCombinerAlgorithm =
+    algorithms::Algorithm<algorithms::Input<edm4hep::TimeSeriesCollection>,
+                          algorithms::Output<edm4hep::TimeSeriesCollection>>;
+
+class PulseCombiner : public PulseCombinerAlgorithm,
+                               public WithPodConfig<PulseCombinerConfig> {
+
+public:
+  PulseCombiner(std::string_view name)
+      : PulseCombinerAlgorithm{name, {"InputPulses"}, {"OutputPulses"}, {}} {}
+  virtual void init() final;
+  void process(const Input&, const Output&) const final;
+
+private:
+
+  std::vector<std::vector<edm4hep::TimeSeries>> clusterPulses(const std::vector<edm4hep::TimeSeries> pulses) const;
+  std::vector<float> sumPulses(const std::vector<edm4hep::TimeSeries> pulses) const;
+  uint64_t m_detector_bitmask = 0xFFFFFFFFFFFFFFFF;
+
+};
+
+} // namespace eicrecon

src/algorithms/digi/PulseCombinerConfig.h

@@ -0,0 +1,16 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2025 Simon Gardner
+
+#pragma once
+
+#include <edm4eic/unit_system.h>
+
+namespace eicrecon {
+
+struct PulseCombinerConfig {
+  double minimum_separation = 50 * edm4eic::unit::ns; // Minimum distance between pulses to keep separate
+  std::string readout = "";
+  std::string combine_field = "";
+};
+
+} // namespace eicrecon

src/algorithms/digi/SiliconPulseGenerationConfig.h

@@ -1,5 +1,5 @@
 // SPDX-License-Identifier: LGPL-3.0-or-later
-// Copyright (C) 2024 Souvik Paul
+// Copyright (C) 2025 Simon Gardner
 
 #pragma once
 

src/detectors/LOWQ2/LOWQ2.cc

@@ -1,10 +1,12 @@
-// SPDX-License-Identifier: LGPL-3.0-or-later
-// Copyright (C) 2023 - 2025, Simon Gardner
+// Copyright 2023-2025, Simon Gardner
+// Subject to the terms in the LICENSE file found in the top-level directory.
+//
+//
 
-#include <JANA/JApplication.h>
 #include <edm4eic/EDM4eicVersion.h>
+#include <JANA/JApplication.h>
+#include <edm4eic/RawTrackerHit.h>
 #include <edm4eic/TrackSegment.h>
-#include <edm4eic/TrackerHit.h>
 #include <edm4eic/unit_system.h>
 #include <fmt/core.h>
 #include <math.h>
@@ -16,15 +18,15 @@
 #include "algorithms/interfaces/WithPodConfig.h"
 #include "algorithms/meta/SubDivideFunctors.h"
 #include "extensions/jana/JOmniFactoryGeneratorT.h"
-#include "factories/digi/PulseNoise_factory.h"
-#include "factories/digi/SiliconPulseGeneration_factory.h"
 #include "factories/digi/SiliconTrackerDigi_factory.h"
+#include "factories/digi/SiliconPulseGeneration_factory.h"
+#include "factories/digi/PulseCombiner_factory.h"
+#include "factories/digi/PulseNoise_factory.h"
 #include "factories/fardetectors/FarDetectorLinearProjection_factory.h"
 #include "factories/fardetectors/FarDetectorLinearTracking_factory.h"
-#include "factories/tracking/TrackerHitReconstruction_factory.h"
 #if EDM4EIC_VERSION_MAJOR >= 8
-#include "factories/fardetectors/FarDetectorTransportationPostML_factory.h"
 #include "factories/fardetectors/FarDetectorTransportationPreML_factory.h"
+#include "factories/fardetectors/FarDetectorTransportationPostML_factory.h"
 #endif
 #include "factories/fardetectors/FarDetectorMLReconstruction_factory.h"
 #include "factories/fardetectors/FarDetectorTrackerCluster_factory.h"
@@ -54,16 +56,27 @@ extern "C" {
       app
     ));
 
+    // Combine pulses into larger pulses
+    app->Add(new JOmniFactoryGeneratorT<PulseCombiner_factory>(
+      "TaggerTrackerPulseCombiner",
+      {"TaggerTrackerHitPulses"},
+      {"TaggerTrackerCombinedPulses"},
+      {
+          .minimum_separation = 25 * edm4eic::unit::ns,
+      },
+      app
+    ));
+
     // Add noise to pulses
     app->Add(new JOmniFactoryGeneratorT<PulseNoise_factory>(
       "TaggerTrackerPulseNoise",
-      {"TaggerTrackerHitPulses"},
-      {"TaggerTrackerHitPulsesWithNoise"},
+      {"TaggerTrackerCombinedPulses"},
+      {"TaggerTrackerCombinedPulsesWithNoise"},
       {
           .poles = 5,
           .variance = 1.0,
           .alpha = 0.5,
-          .scale = 500.0,
+          .scale = 0.000002,
       },
       app
     ));
@@ -85,17 +98,6 @@ extern "C" {
          app
     ));
 
-    // Convert raw digitized hits into hits with geometry info (ready for tracking)
-    app->Add(new JOmniFactoryGeneratorT<TrackerHitReconstruction_factory>(
-      "TaggerTrackerRecHits",
-      {"TaggerTrackerRawHits"},
-      {"TaggerTrackerRecHits"},
-      {
-          .timeResolution = 2,
-      },
-      app
-    ));
-
     // Divide collection based on geometry segmentation labels
     // This should really be done before digitization as summing hits in the same cell couldn't even be mixed between layers. At the moment just prep for clustering.
     std::string readout = "TaggerTrackerHits";
@@ -113,15 +115,15 @@ extern "C" {
       moduleClusterTags.push_back({});
       for(int lay_id : layerIDs){
         geometryDivisions.push_back({mod_id,lay_id});
-        geometryDivisionCollectionNames.push_back(fmt::format("TaggerTrackerM{}L{}RecHits",mod_id,lay_id));
+        geometryDivisionCollectionNames.push_back(fmt::format("TaggerTrackerM{}L{}RawHits",mod_id,lay_id));
         outputClusterCollectionNames.push_back(fmt::format("TaggerTrackerM{}L{}ClusterPositions",mod_id,lay_id));
         moduleClusterTags.back().push_back(outputClusterCollectionNames.back());
       }
     }
 
-    app->Add(new JOmniFactoryGeneratorT<SubDivideCollection_factory<edm4eic::TrackerHit>>(
+    app->Add(new JOmniFactoryGeneratorT<SubDivideCollection_factory<edm4eic::RawTrackerHit>>(
          "TaggerTrackerSplitHits",
-         {"TaggerTrackerRecHits"},
+         {"TaggerTrackerRawHits"},
          geometryDivisionCollectionNames,
          {
           .function = GeometrySplit{geometryDivisions,readout,geometryLabels},

src/factories/digi/PulseCombiner_factory.h

@@ -0,0 +1,46 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2025 Simon Gardner
+
+#pragma once
+
+#include "extensions/jana/JOmniFactory.h"
+
+#include "algorithms/digi/PulseCombiner.h"
+#include <iostream>
+
+namespace eicrecon {
+
+class PulseCombiner_factory
+    : public JOmniFactory<PulseCombiner_factory, PulseCombinerConfig> {
+public:
+  using AlgoT = eicrecon::PulseCombiner;
+
+private:
+  std::unique_ptr<AlgoT> m_algo;
+
+  PodioInput<edm4hep::TimeSeries> m_in_pulses{this};
+
+  PodioOutput<edm4hep::TimeSeries> m_out_pulses{this};
+
+  ParameterRef<double>      m_minimum_separation{this, "minimumSeperation", config().minimum_separation};
+  ParameterRef<std::string> m_readout{this, "readout", config().readout};
+  ParameterRef<std::string> m_combine_field{this, "combineField", config().combine_field};
+
+  Service<AlgorithmsInit_service> m_algorithmsInit{this};
+
+public:
+  void Configure() {
+    m_algo = std::make_unique<eicrecon::PulseCombiner>(GetPrefix());
+    m_algo->level(static_cast<algorithms::LogLevel>(logger()->level()));
+    m_algo->applyConfig(config());
+    m_algo->init();
+  }
+
+  void ChangeRun(int64_t run_number) {}
+
+  void Process(int64_t run_number, uint64_t event_number) {
+    m_algo->process({m_in_pulses()}, {m_out_pulses().get()});
+  }
+};
+
+} // namespace eicrecon

src/services/io/podio/JEventProcessorPODIO.cc

@@ -137,7 +137,8 @@ JEventProcessorPODIO::JEventProcessorPODIO() {
             // LOWQ2 hits
             "TaggerTrackerHits",
             "TaggerTrackerHitPulses",
-            "TaggerTrackerHitPulsesWithNoise",
+            "TaggerTrackerCombinedPulses",
+            "TaggerTrackerCombinedPulsesWithNoise",
             "TaggerTrackerRawHits",
             "TaggerTrackerRawHitAssociations",
             "TaggerTrackerM1L0ClusterPositions",