Realistic BTOF digitization

src/algorithms/digi/BTOFChargeSharing.cc

@@ -0,0 +1,180 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2024 Chun Yuen Tsang, Prithwish Tribedy
+//
+// Spread energy deposition from one strip to neighboring strips within sensor boundaries
+
+#include <DD4hep/DetElement.h>
+#include <DD4hep/Handle.h>
+#include <DD4hep/Readout.h>
+#include <DD4hep/Segmentations.h>
+#include <DD4hep/Volumes.h>
+#include <Evaluator/DD4hepUnits.h>
+#include <Math/GenVector/Cartesian3D.h>
+#include <Math/GenVector/DisplacementVector3D.h>
+#include <Math/SpecFuncMathCore.h>
+#include <TGeoManager.h>
+#include <TGeoMatrix.h>
+#include <TGeoVolume.h>
+#include <algorithms/geo.h>
+#include <edm4hep/Vector3d.h>
+#include <edm4hep/Vector3f.h>
+#include <gsl/pointers>
+#include <stdexcept>
+#include <utility>
+#include <vector>
+
+#include "BTOFChargeSharing.h"
+#include "DD4hep/Detector.h"
+#include "TMath.h"
+#include "algorithms/digi/TOFHitDigiConfig.h"
+
+
+namespace eicrecon {
+
+void BTOFChargeSharing::init() {
+  m_detector = algorithms::GeoSvc::instance().detector();;
+  m_converter = algorithms::GeoSvc::instance().cellIDPositionConverter();
+
+  auto seg       = m_detector->readout(m_cfg.readout).segmentation();
+  auto type      = seg.type();
+  if (type != "CartesianGridXY")
+    throw std::runtime_error("Unsupported segmentation type: " + type +
+                             ". BarrelTOF must use CartesianGridXY.");
+  // retrieve meaning of cellID bits
+  m_decoder = seg.decoder();
+}
+
+void BTOFChargeSharing::process(const BTOFChargeSharing::Input& input,
+                                const BTOFChargeSharing::Output& output) const {
+  const auto [simhits] = input;
+  auto [sharedHits] = output;
+  std::shared_ptr<std::vector<dd4hep::rec::CellID>> neighbors;
+
+  for(const auto & hit : *simhits) {
+    auto cellID = hit.getCellID();
+
+    if(!neighbors){
+      std::unordered_set<dd4hep::rec::CellID> dp;
+      neighbors = std::make_shared<std::vector<dd4hep::rec::CellID>>();
+      this -> _findAllNeighborsInSensor(cellID, neighbors, dp);
+    }
+
+    double edep = hit.getEDep();
+    double time = hit.getTime();
+    auto momentum = hit.getMomentum();
+    auto truePos = hit.getPosition();
+    auto localPos_hit = this -> _global2Local(dd4hep::Position(truePos.x*dd4hep::mm, truePos.y*dd4hep::mm, truePos.z*dd4hep::mm));
+
+    for(const auto neighbor : *neighbors) {
+       // integrate over neighbor area to get total energy deposition
+       auto localPos_neighbor = this -> _cell2LocalPosition(neighbor);
+       auto cellDimension = m_converter -> cellDimensions(neighbor);
+
+       double edep_cell = edep *
+                      _integralGaus(localPos_hit.x(), m_cfg.sigma_sharingx,
+                                    localPos_neighbor.x() - 0.5 * cellDimension[0],
+                                    localPos_neighbor.x() + 0.5 * cellDimension[0]) *
+                      _integralGaus(localPos_hit.y(), m_cfg.sigma_sharingy,
+                                    localPos_neighbor.y() - 0.5 * cellDimension[1],
+                                    localPos_neighbor.y() + 0.5 * cellDimension[1]);
+
+       if(edep_cell > 0) {
+         auto globalPos = m_converter -> position(neighbor);
+         auto hit =  sharedHits->create();
+         hit.setCellID(neighbor);
+         hit.setEDep(edep_cell);
+         hit.setTime(time);
+         hit.setPosition({globalPos.x(), globalPos.y(), globalPos.z()});
+         hit.setMomentum({momentum.x, momentum.y, momentum.z});
+       }
+    }
+  }
+} // BTOFChargeSharing:process
+
+
+void BTOFChargeSharing::_findAllNeighborsInSensor(
+    dd4hep::rec::CellID hitCell, std::shared_ptr<std::vector<dd4hep::rec::CellID>>& answer,
+    std::unordered_set<dd4hep::rec::CellID>& dp) const {
+  // use MST to find all neighbor within a sensor
+  // I can probably write down the formula by hand, but why do things manually when computer do
+  // everything for you?
+  const std::vector<std::pair<int, int>> searchDirs{{0, 1}, {0, -1}, {1, 0}, {-1, 0}};
+  answer->push_back(hitCell);
+  dp.insert(hitCell);
+
+  auto sensorID = this -> _getSensorID(hitCell);
+  auto xID      = m_decoder->get(hitCell, "x");
+  auto yID      = m_decoder->get(hitCell, "y");
+  for (const auto& dir : searchDirs) {
+    auto testCell = hitCell;
+    try {
+      m_decoder->set(testCell, "x", xID + dir.first);
+      m_decoder->set(testCell, "y", yID + dir.second);
+    } catch (const std::runtime_error& err) {
+      // catch overflow error
+      // ignore if invalid position ID
+      continue;
+    }
+
+    try {
+      auto pos = m_converter->position(testCell);
+      if (testCell != m_converter->cellID(pos))
+        continue;
+    } catch (const std::invalid_argument& err) {
+      // Ignore CellID that is invalid
+      continue;
+    }
+
+    // only look for cells that have not been searched
+    if (dp.find(testCell) == dp.end()) {
+      auto testSensorID = _getSensorID(testCell);
+      if (testSensorID == sensorID) {
+        // inside the same sensor
+        this->_findAllNeighborsInSensor(testCell, answer, dp);
+      }
+    }
+  }
+}
+
+const dd4hep::rec::CellID
+BTOFChargeSharing::_getSensorID(const dd4hep::rec::CellID& hitCell) const {
+  // fix x-y, what you left with are ids that corresponds to sensor info
+  // cellID may change when position changes.
+  auto sensorID = hitCell; //_converter -> cellID(_converter -> position(hitCell));
+  m_decoder->set(sensorID, "x", 0);
+  m_decoder->set(sensorID, "y", 0);
+
+  return sensorID;
+}
+
+double BTOFChargeSharing::_integralGaus(double mean, double sd, double low_lim, double up_lim) const {
+  // return integral Gauss(mean, sd) dx from x = low_lim to x = up_lim
+  // default value is set when sd = 0
+  double up = mean > up_lim? -0.5 : 0.5;
+  double low = mean > low_lim? -0.5 : 0.5;
+  if(sd > 0) {
+    up  = -0.5 * std::erf(std::sqrt(2) * (mean - up_lim) / sd);
+    low = -0.5 * std::erf(std::sqrt(2) * (mean - low_lim) / sd);
+  }
+  return up - low;
+}
+
+dd4hep::Position BTOFChargeSharing::_cell2LocalPosition(const dd4hep::rec::CellID& cell) const {
+  auto position   = m_converter -> position(cell); // global position
+  return this -> _global2Local(position);
+}
+
+dd4hep::Position BTOFChargeSharing::_global2Local(const dd4hep::Position& pos) const {
+  auto geoManager = m_detector->world().volume()->GetGeoManager();
+  auto node       = geoManager->FindNode(pos.x(), pos.y(), pos.z());
+  auto currMatrix = geoManager->GetCurrentMatrix();
+
+  double g[3], l[3];
+  pos.GetCoordinates(g);
+  currMatrix->MasterToLocal(g, l);
+  dd4hep::Position position;
+  position.SetCoordinates(l);
+  return position;
+}
+
+} // namespace eicrecon

src/algorithms/digi/BTOFChargeSharing.h

@@ -0,0 +1,55 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2024 Chun Yuen Tsang, Prithwish Tribedy
+//
+// Spread energy deposition from one strip to neighboring strips within sensor boundaries
+
+#pragma once
+
+#include <DD4hep/Objects.h>
+#include <DDRec/CellIDPositionConverter.h>
+#include <DDSegmentation/BitFieldCoder.h>
+#include <algorithms/algorithm.h>
+#include <edm4hep/SimTrackerHitCollection.h>
+#include <memory>
+#include <string>
+#include <string_view>
+#include <unordered_set>
+#include <vector>
+
+#include "DD4hep/Detector.h"
+#include "algorithms/digi/TOFHitDigiConfig.h"
+#include "algorithms/interfaces/WithPodConfig.h"
+
+namespace eicrecon {
+
+using BTOFChargeSharingAlgorithm =
+        algorithms::Algorithm<algorithms::Input<edm4hep::SimTrackerHitCollection>,
+                              algorithms::Output<edm4hep::SimTrackerHitCollection>>;
+
+class BTOFChargeSharing : public BTOFChargeSharingAlgorithm,
+                          public WithPodConfig<TOFHitDigiConfig> {
+
+public:
+  BTOFChargeSharing(std::string_view name) : BTOFChargeSharingAlgorithm{name,
+                                                  {"TOFBarrelHits"},
+                                                  {"TOFBarrelSharedHits"},
+                                                  ""} {};
+
+  void init() final;
+  void process(const Input&, const Output&) const final;
+protected:
+  void _findAllNeighborsInSensor(dd4hep::rec::CellID hitCell,
+                                 std::shared_ptr<std::vector<dd4hep::rec::CellID>>& answer,
+                                 std::unordered_set<dd4hep::rec::CellID>& dp) const;
+  const dd4hep::rec::CellID _getSensorID(const dd4hep::rec::CellID& hitCell) const;
+  double _integralGaus(double mean, double sd, double low_lim, double up_lim) const;
+  dd4hep::Position _cell2LocalPosition(const dd4hep::rec::CellID& cell) const;
+  dd4hep::Position _global2Local(const dd4hep::Position& pos) const;
+
+  const dd4hep::DDSegmentation::BitFieldCoder* m_decoder  = nullptr;
+  const dd4hep::Detector* m_detector                      = nullptr;
+  const dd4hep::rec::CellIDPositionConverter* m_converter = nullptr;
+
+};
+
+} // namespace eicrecon

src/algorithms/digi/TOFHitDigiConfig.h

@@ -0,0 +1,40 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2024 Souvik Paul
+
+#pragma once
+
+#include <DD4hep/DD4hepUnits.h>
+
+namespace eicrecon {
+
+struct TOFHitDigiConfig {
+  // single hit energy deposition threshold
+  double threshold{1.0 * dd4hep::keV};
+  double tRes = 0.1; /// TODO 8 of what units??? Same TODO in juggler. Probably [ns]
+                     // digitization settings
+  double resolutionTDC{1};
+  double resolutionADC{1};
+
+  // Parameters of AC-LGAD signal generation
+  double gain           = 80;
+  double risetime       = 0.45; // 0.02; //in ns
+  double sigma_analog   = 0.293951;
+  double sigma_sharingx = 0.1;
+  double sigma_sharingy = 0.5;
+  double Vm             = -1e-4 * dd4hep::GeV; // Vm = voltage maximum. When EDep = 1e-4 GeV, voltage corresponds to ADC = adc_max
+  double t_thres        = 0.1 * Vm;
+  double ignore_thres   = 0.01 * Vm; // If EDep below this value, digitization for the cell will be ignored. Speed up calculation
+                                     //
+  double tMin     = 0.1;
+  double tMax     = 25;// 25 ns is the period of 40MHz EIC clock
+  int total_time  = ceil(tMax - tMin);
+  int adc_bit     = 8;
+  int tdc_bit     = 10;
+
+  int adc_range = pow(2, adc_bit);
+  int tdc_range = pow(2, tdc_bit);
+
+  std::string readout = "TOFBarrelHits";
+};
+
+} // namespace eicrecon

src/algorithms/digi/TOFPulseDigitization.cc

@@ -0,0 +1,65 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2024 Souvik Paul, Chun Yuen Tsang, Prithwish Tribedy
+// Special Acknowledgement: Kolja Kauder
+//
+// Convert ADC pulses from TOFPulseGeneration into ADC and TDC values
+
+#include <math.h>
+#include <podio/RelationRange.h>
+#include <stdlib.h>
+#include <algorithm>
+#include <gsl/pointers>
+#include <limits>
+#include <vector>
+
+#include "TOFPulseDigitization.h"
+#include "algorithms/digi/TOFHitDigiConfig.h"
+
+
+namespace eicrecon {
+
+void TOFPulseDigitization::process(const TOFPulseDigitization::Input& input,
+                                   const TOFPulseDigitization::Output& output) const {
+  const auto [simhits] = input;
+  auto [rawhits] = output;
+
+  double thres = m_cfg.t_thres;
+  // Vm in unit of GeV. When Edep = Vm, ADC = cfg.adc_range-1
+  double Vm = m_cfg.Vm;
+  int adc_range = m_cfg.adc_range;
+
+  // normalized time threshold
+  // convert threshold EDep to voltage
+  double norm_threshold = -thres * adc_range / Vm;
+
+  for(const auto& pulse : *simhits) {
+    double intersectionX = 0.0;
+    int tdc              = std::numeric_limits<int>::max();
+    int adc              = 0;
+    double V             = 0.0;
+
+    int time_bin = 0;
+    double adc_prev = 0;
+    double time_interval = pulse.getInterval();
+    auto adcs = pulse.getAdcCounts();
+    for (const auto adc : adcs) {
+      if (adc_prev >= norm_threshold && adc <= norm_threshold) {
+        intersectionX = time_bin*time_interval + time_interval * (norm_threshold - adc_prev) / (adc - adc_prev);
+        tdc = static_cast<int>(intersectionX / time_interval);
+      }
+      if (abs(adc) > abs(V)) // To get peak of the Analog signal
+        V = adc;
+      adc_prev = adc;
+      ++time_bin;
+    }
+
+    // limit the range of adc values
+    adc = std::min(static_cast<double>(adc_range), round(-V));
+    // only store valid hits
+    if (tdc < std::numeric_limits<int>::max())
+      rawhits->create(pulse.getCellID(), adc, tdc);
+    //-----------------------------------------------------------
+
+  }
+} // TOFPulseDigitization:process
+} // namespace eicrecon

src/algorithms/digi/TOFPulseDigitization.h

@@ -0,0 +1,37 @@
+// SPDX-License-Identifier: LGPL-3.0-or-later
+// Copyright (C) 2024 Souvik Paul, Chun Yuen Tsang, Prithwish Tribedy
+// Special Acknowledgement: Kolja Kauder
+//
+// Convert ADC pulses from TOFPulseGeneration into ADC and TDC values
+
+#pragma once
+
+#include <algorithms/algorithm.h>
+#include <edm4eic/RawTrackerHitCollection.h>
+#include <edm4hep/RawTimeSeriesCollection.h>
+#include <string>
+#include <string_view>
+
+#include "algorithms/digi/TOFHitDigiConfig.h"
+#include "algorithms/interfaces/WithPodConfig.h"
+
+namespace eicrecon {
+
+using TOFPulseDigitizationAlgorithm =
+        algorithms::Algorithm<algorithms::Input<edm4hep::RawTimeSeriesCollection>,
+                              algorithms::Output<edm4eic::RawTrackerHitCollection>>;
+
+class TOFPulseDigitization : public TOFPulseDigitizationAlgorithm,
+                             public WithPodConfig<TOFHitDigiConfig> {
+
+public:
+  TOFPulseDigitization(std::string_view name) :
+          TOFPulseDigitizationAlgorithm{name,
+                                        {"TOFBarrelPulse"},
+                                        {"TOFBarrelADCTDC"},
+                                        {}} {}
+  void init() {};
+  void process(const Input&, const Output&) const final;
+};
+
+} // namespace eicrecon