Add LGADHitReconstruction that does clustering

src/algorithms/digi/EICROCDigitization.cc

@@ -26,9 +26,8 @@ void EICROCDigitization::process(const EICROCDigitization::Input& input,
   int adc_range = m_cfg.adc_range;
 
   for (const auto& pulse : *simhits) {
-    int tdc  = std::numeric_limits<int>::max();
-    int adc  = 0;
-    double V = 0.0;
+    int tdc = std::numeric_limits<int>::max();
+    int V   = 0;
 
     int time_bin          = 0;
     double adc_prev       = 0;
@@ -45,7 +44,7 @@ void EICROCDigitization::process(const EICROCDigitization::Input& input,
     }
 
     // limit the range of adc values
-    adc = std::min(static_cast<double>(adc_range), std::round(-V));
+    int adc = std::min(adc_range, -V);
     // only store valid hits
     if (tdc < std::numeric_limits<int>::max())
       rawhits->create(pulse.getCellID(), adc, tdc);

src/algorithms/reco/LGADHitCalibration.cc

@@ -0,0 +1,67 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2024 Chun Yuen Tsang
+
+#include "LGADHitCalibration.h"
+
+#include <Evaluator/DD4hepUnits.h>
+#include <Math/GenVector/Cartesian3D.h>
+#include <Math/GenVector/DisplacementVector3D.h>
+#include <algorithms/geo.h>
+#include <edm4eic/CovDiag3f.h>
+#include <edm4hep/Vector3f.h>
+#include <algorithm>
+#include <cmath>
+#include <gsl/pointers>
+#include <vector>
+
+#include "algorithms/reco/LGADHitCalibrationConfig.h"
+
+namespace eicrecon {
+
+void LGADHitCalibration::init() {
+  m_converter = algorithms::GeoSvc::instance().cellIDPositionConverter();
+}
+
+void LGADHitCalibration::process(const LGADHitCalibration::Input& input,
+                                 const LGADHitCalibration::Output& output) const {
+  using dd4hep::mm;
+
+  const auto [TDCADC_hits] = input;
+  auto [calibrated_hits]   = output;
+
+  for (const auto& TDCADC_hit : *TDCADC_hits) {
+
+    auto id = TDCADC_hit.getCellID();
+
+    // Get position and dimension
+    auto pos   = m_converter->position(id);
+    double ADC = TDCADC_hit.getCharge();
+    double TDC = TDCADC_hit.getTimeStamp();
+
+    // adc to charge
+    double charge = ADC * m_cfg.c_slope + m_cfg.c_intercept;
+    // TDC to time
+    float time = TDC * m_cfg.t_slope + m_cfg.t_intercept;
+
+    auto cellSize = m_converter->cellDimensions(id);
+    // sqrt(12) factor convertes ranges of uniform distribution to it's standard deviation
+    double varX = cellSize[0] / mm / std::sqrt(12.);
+    varX *= varX; // square of cell size
+    double varY = cellSize[1] / mm / std::sqrt(12.);
+    varY *= varY;
+    double varZ = cellSize.size() > 2 ? cellSize[2] / mm / std::sqrt(12.) : 0;
+    varZ *= varZ;
+
+    calibrated_hits->create(
+        id,
+        edm4hep::Vector3f{static_cast<float>(pos.x()), static_cast<float>(pos.y()),
+                          static_cast<float>(pos.z())},
+        edm4eic::CovDiag3f{varX, varY, varZ}, // should be the covariance of position
+        time,                                 // ns
+        0.0F,                                 // covariance of time
+        std::max(0., charge),                 // total ADC sum
+        0.0F);                                // Error on the energy
+  }
+}
+
+} // namespace eicrecon

src/algorithms/reco/LGADHitCalibration.h

@@ -0,0 +1,39 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2024 Chun Yuen Tsang
+
+#pragma once
+
+#include <DDRec/CellIDPositionConverter.h>
+#include <algorithms/algorithm.h>
+#include <edm4eic/RawTrackerHitCollection.h>
+#include <edm4eic/TrackerHitCollection.h>
+#include <spdlog/logger.h>
+#include <memory>
+#include <string>
+#include <string_view>
+
+#include "LGADHitCalibrationConfig.h"
+#include "algorithms/interfaces/WithPodConfig.h"
+
+namespace eicrecon {
+
+using LGADHitCalibrationAlgorithm =
+    algorithms::Algorithm<algorithms::Input<edm4eic::RawTrackerHitCollection>,
+                          algorithms::Output<edm4eic::TrackerHitCollection>>;
+
+class LGADHitCalibration : public LGADHitCalibrationAlgorithm,
+                           public WithPodConfig<LGADHitCalibrationConfig> {
+
+public:
+  LGADHitCalibration(std::string_view name)
+      : LGADHitCalibrationAlgorithm{name, {"TOFBarrelADCTDC"}, {"TOFBarrelCalHit"}, ""} {};
+
+  void init() final;
+  void process(const Input&, const Output&) const final;
+
+private:
+  std::shared_ptr<spdlog::logger> m_log;
+
+  const dd4hep::rec::CellIDPositionConverter* m_converter = nullptr;
+};
+} // namespace eicrecon

src/algorithms/reco/LGADHitCalibrationConfig.h

@@ -0,0 +1,13 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2025 Chun Yuen Tsang
+
+#pragma once
+
+namespace eicrecon {
+struct LGADHitCalibrationConfig {
+  // parameters that convert ADC to EDep
+  double c_slope = 3.86976e-7, c_intercept = 2.42716e-5;
+  // parameters that convert TDC to hit time (ns)
+  double t_slope = 0.0197305, t_intercept = 0.208047;
+};
+} // namespace eicrecon

src/algorithms/reco/LGADHitClusterAssociation.cc

@@ -0,0 +1,107 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2024 Chun Yuen Tsang
+
+#include "LGADHitClusterAssociation.h"
+
+#include <DD4hep/Alignments.h>
+#include <DD4hep/DetElement.h>
+#include <DD4hep/Detector.h>
+#include <DD4hep/Readout.h>
+#include <DD4hep/Segmentations.h>
+#include <Evaluator/DD4hepUnits.h>
+#include <TGeoMatrix.h>
+#include <algorithms/geo.h>
+#include <edm4eic/CovDiag3f.h>
+#include <edm4hep/Vector3f.h>
+#include <algorithm>
+#include <gsl/pointers>
+#include <unordered_map>
+#include <utility>
+#include <vector>
+
+#include "algorithms/reco/LGADHitClusterAssociationConfig.h"
+
+namespace eicrecon {
+
+void LGADHitClusterAssociation::init() {
+  auto detector = algorithms::GeoSvc::instance().detector();
+  auto seg      = detector->readout(m_cfg.readout).segmentation();
+  m_converter   = algorithms::GeoSvc::instance().cellIDPositionConverter();
+  m_decoder     = seg.decoder();
+}
+
+dd4hep::rec::CellID LGADHitClusterAssociation::getSensorInfos(const dd4hep::rec::CellID& id) const {
+  // CellID for BarrelTOF is composed of 6 parts
+  // system, layer, module, sensor, x, y
+  // If we fix x and y to zero, what remains will be the sensor information only
+  auto id_return = id;
+  m_decoder->set(id_return, "x", 0);
+  m_decoder->set(id_return, "y", 0);
+  return id_return;
+}
+
+edm4hep::Vector3f
+LGADHitClusterAssociation::_local2Global(const dd4hep::VolumeManagerContext* context,
+                                         const edm4hep::Vector2f& locPos) const {
+  auto nodeMatrix = context->element.nominal().worldTransformation();
+
+  double g[3], l[3];
+  l[0] = locPos.a * dd4hep::mm;
+  l[1] = locPos.b * dd4hep::mm;
+  l[2] = 0;
+  nodeMatrix.LocalToMaster(l, g);
+  return edm4hep::Vector3f{static_cast<float>(g[0] / dd4hep::mm),
+                           static_cast<float>(g[1] / dd4hep::mm),
+                           static_cast<float>(g[2] / dd4hep::mm)};
+}
+
+void LGADHitClusterAssociation::process(const LGADHitClusterAssociation::Input& input,
+                                        const LGADHitClusterAssociation::Output& output) const {
+  using dd4hep::mm;
+  const auto [meas2D_hits, raw_hits] = input;
+  auto [asso_hits]                   = output;
+
+  // group raw hit by cell ID
+  std::unordered_map<dd4hep::rec::CellID, std::vector<edm4eic::MutableTrackerHit>> cellHitMap;
+
+  for (const auto& meas2D_hit : *meas2D_hits) {
+    auto cellID = meas2D_hit.getSurface();
+    auto time   = meas2D_hit.getTime();
+
+    // sum ADC info
+    double tot_charge = 0;
+    // position info
+    auto locPos         = meas2D_hit.getLoc();
+    const auto* context = m_converter->findContext(cellID);
+    auto ave_pos        = this->_local2Global(context, locPos);
+    auto asso_hit = asso_hits->create(cellID, ave_pos, edm4eic::CovDiag3f{0., 0., 0.}, time, 0.,
+                                      tot_charge, 0.);
+    cellHitMap[getSensorInfos(asso_hit.getCellID())].push_back(asso_hit);
+  }
+
+  // sort the tracker hits by time
+  for (auto& [cellID, hits] : cellHitMap)
+    std::sort(hits.begin(), hits.end(),
+              [](const edm4eic::MutableTrackerHit& a, const edm4eic::MutableTrackerHit& b) {
+                return a.getTime() < b.getTime();
+              });
+
+  // get the associated raw_hits by the closest time
+  for (const auto& raw_hit : *raw_hits) {
+    auto time     = raw_hit.getTimeStamp() / 1e3; // ps->ns
+    auto sensorID = getSensorInfos(raw_hit.getCellID());
+    auto it       = cellHitMap.find(sensorID);
+    if (it != cellHitMap.end()) {
+      auto& hits = it->second;
+      // if t_raw - assoDeltaT < t_hit < t_raw + assoDeltaT, we make the association
+      auto itVec = std::lower_bound(
+          hits.begin(), hits.end(), time - m_cfg.assoDeltaT,
+          [](const edm4eic::MutableTrackerHit& a, double t) { return a.getTime() < t; });
+      for (; itVec != hits.end() && itVec->getTime() <= time + m_cfg.assoDeltaT; ++itVec) {
+        itVec->setRawHit(raw_hit);
+      }
+    }
+  }
+}
+
+} // namespace eicrecon

src/algorithms/reco/LGADHitClusterAssociation.h

@@ -0,0 +1,50 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2024 Chun Yuen Tsang
+
+#pragma once
+
+#include <DD4hep/VolumeManager.h>
+#include <DDRec/CellIDPositionConverter.h>
+#include <DDSegmentation/BitFieldCoder.h>
+#include <algorithms/algorithm.h>
+#include <edm4eic/Measurement2DCollection.h>
+#include <edm4eic/RawTrackerHitCollection.h>
+#include <edm4eic/TrackerHitCollection.h>
+#include <edm4hep/Vector2f.h>
+#include <edm4hep/Vector3f.h>
+#include <spdlog/logger.h>
+#include <memory>
+#include <string>
+#include <string_view>
+
+#include "LGADHitClusterAssociationConfig.h"
+#include "algorithms/interfaces/WithPodConfig.h"
+
+namespace eicrecon {
+
+using LGADHitClusterAssociationAlgorithm = algorithms::Algorithm<
+    algorithms::Input<edm4eic::Measurement2DCollection, edm4eic::RawTrackerHitCollection>,
+    algorithms::Output<edm4eic::TrackerHitCollection>>;
+
+class LGADHitClusterAssociation : public LGADHitClusterAssociationAlgorithm,
+                                  public WithPodConfig<LGADHitClusterAssociationConfig> {
+
+public:
+  LGADHitClusterAssociation(std::string_view name)
+      : LGADHitClusterAssociationAlgorithm{
+            name, {"TOFBarrelClusterHits", "TOFBarrelRawHits"}, {"TOFBarrelRecHits"}, ""} {};
+
+  void init() final;
+  void process(const Input&, const Output&) const final;
+
+private:
+  dd4hep::rec::CellID getSensorInfos(const dd4hep::rec::CellID& id) const;
+  edm4hep::Vector3f _local2Global(const dd4hep::VolumeManagerContext* context,
+                                  const edm4hep::Vector2f& locPos) const;
+
+  std::shared_ptr<spdlog::logger> m_log;
+  /// fetch sensor information from cellID
+  const dd4hep::DDSegmentation::BitFieldCoder* m_decoder  = nullptr;
+  const dd4hep::rec::CellIDPositionConverter* m_converter = nullptr;
+};
+} // namespace eicrecon

src/algorithms/reco/LGADHitClusterAssociationConfig.h

@@ -0,0 +1,13 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2025 Chun Yuen Tsang
+
+#pragma once
+
+#include <edm4eic/unit_system.h>
+
+namespace eicrecon {
+struct LGADHitClusterAssociationConfig {
+  std::string readout = "TOFBarrelHits";
+  double assoDeltaT   = 1 * edm4eic::unit::ns;
+};
+} // namespace eicrecon