diff --git a/src/algorithms/digi/BTOFChargeSharing.cc b/src/algorithms/digi/BTOFChargeSharing.cc new file mode 100644 index 0000000000..9210cf6334 --- /dev/null +++ b/src/algorithms/digi/BTOFChargeSharing.cc @@ -0,0 +1,179 @@ +// 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 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "BTOFChargeSharing.h" +#include "DD4hep/Detector.h" +#include "algorithms/digi/BTOFChargeSharingConfig.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> neighbors; + + for (const auto& hit : *simhits) { + auto cellID = hit.getCellID(); + + if (!neighbors) { + std::unordered_set dp; + neighbors = std::make_shared>(); + 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>& answer, + std::unordered_set& 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> 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 diff --git a/src/algorithms/digi/BTOFChargeSharing.h b/src/algorithms/digi/BTOFChargeSharing.h new file mode 100644 index 0000000000..ff7c954249 --- /dev/null +++ b/src/algorithms/digi/BTOFChargeSharing.h @@ -0,0 +1,53 @@ +// 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 +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "DD4hep/Detector.h" +#include "algorithms/digi/BTOFChargeSharingConfig.h" +#include "algorithms/interfaces/WithPodConfig.h" + +namespace eicrecon { + +using BTOFChargeSharingAlgorithm = + algorithms::Algorithm, + algorithms::Output>; + +class BTOFChargeSharing : public BTOFChargeSharingAlgorithm, + public WithPodConfig { + +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>& answer, + std::unordered_set& 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 diff --git a/src/algorithms/digi/BTOFChargeSharingConfig.h b/src/algorithms/digi/BTOFChargeSharingConfig.h new file mode 100644 index 0000000000..0d42fa70a9 --- /dev/null +++ b/src/algorithms/digi/BTOFChargeSharingConfig.h @@ -0,0 +1,18 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2024 Souvik Paul + +#pragma once + +#include + +namespace eicrecon { + +struct BTOFChargeSharingConfig { + // Parameters of AC-LGAD signal generation + double sigma_sharingx = 0.1 * dd4hep::cm; + double sigma_sharingy = 0.5 * dd4hep::cm; + + std::string readout = "TOFBarrelHits"; +}; + +} // namespace eicrecon diff --git a/src/algorithms/digi/LGADChargeSharingConfig.h b/src/algorithms/digi/LGADChargeSharingConfig.h new file mode 100644 index 0000000000..b7cb9f4641 --- /dev/null +++ b/src/algorithms/digi/LGADChargeSharingConfig.h @@ -0,0 +1,39 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2024 Souvik Paul + +#pragma once + +#include + +namespace eicrecon { + +struct TOFHitDigiConfig { + // Parameters of AC-LGAD signal generation + double gain = 113.755; + double risetime = 0.45 * dd4hep::ns; + double sigma_analog = 0.293951 * dd4hep::ns; + double sigma_sharingx = 0.1 * dd4hep::cm; + double sigma_sharingy = 0.5 * dd4hep::cm; + 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; // TDC value = time when pulse exceed t_thres + double ignore_thres = 0.001 * Vm; // If EDep below this value, digitization for the cell will be + // ignored. Speed up calculation + int adc_bit = 8; + int tdc_bit = 10; + + // total number of TDC/ADC values + // Since digitizer starts at zero, max ADC value = adc_range - 1 + // Similar for TDC + int adc_range = pow(2, adc_bit); + int tdc_range = pow(2, tdc_bit); + + // period of the sensor clock. Time internal to sensor will all be digitized to integer multiple + // of tInterval + double tMax = 25 * dd4hep::ns; // 25 ns is the period of 40MHz EIC clock + double tInterval = tMax / (tdc_range - 1); + + std::string readout = "TOFBarrelHits"; +}; + +} // namespace eicrecon diff --git a/src/algorithms/digi/LGADPulseDigitization.cc b/src/algorithms/digi/LGADPulseDigitization.cc new file mode 100644 index 0000000000..ba426165df --- /dev/null +++ b/src/algorithms/digi/LGADPulseDigitization.cc @@ -0,0 +1,59 @@ +// 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 LGADPulseGeneration into ADC and TDC values + +#include +#include +#include +#include +#include +#include +#include + +#include "LGADPulseDigitization.h" +#include "algorithms/digi/LGADPulseDigitizationConfig.h" + +namespace eicrecon { + +void LGADPulseDigitization::process(const LGADPulseDigitization::Input& input, + const LGADPulseDigitization::Output& output) const { + const auto [simhits] = input; + auto [rawhits] = output; + + double thres = -m_cfg.t_thres; + int adc_range = m_cfg.adc_range; + + for (const auto& pulse : *simhits) { + double intersectionX = 0.0; + int tdc = std::numeric_limits::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(); + double n_EICROC_cycle = static_cast(std::floor(pulse.getTime()/m_cfg.tMax + 1e-3)); + for (const auto adc : adcs) { + if (adc_prev >= thres && adc <= thres) { + intersectionX = time_bin * time_interval + + time_interval * (thres - adc_prev) / (adc - adc_prev); + tdc = static_cast(intersectionX / time_interval) + n_EICROC_cycle * m_cfg.tdc_range; + } + 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(adc_range), round(-V)); + // only store valid hits + if (tdc < std::numeric_limits::max()) + rawhits->create(pulse.getCellID(), adc, tdc); + //----------------------------------------------------------- + } +} // LGADPulseDigitization:process +} // namespace eicrecon diff --git a/src/algorithms/digi/LGADPulseDigitization.h b/src/algorithms/digi/LGADPulseDigitization.h new file mode 100644 index 0000000000..09ee182b9e --- /dev/null +++ b/src/algorithms/digi/LGADPulseDigitization.h @@ -0,0 +1,34 @@ +// 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 LGADPulseGeneration into ADC and TDC values + +#pragma once + +#include +#include +#include +#include +#include + +#include "algorithms/digi/LGADPulseDigitizationConfig.h" +#include "algorithms/interfaces/WithPodConfig.h" + +namespace eicrecon { + +using LGADPulseDigitizationAlgorithm = + algorithms::Algorithm, + algorithms::Output>; + +class LGADPulseDigitization : public LGADPulseDigitizationAlgorithm, + public WithPodConfig { + +public: + LGADPulseDigitization(std::string_view name) + : LGADPulseDigitizationAlgorithm{name, {"LGADPulse"}, {"ADCTDCOutput"}, {}} {} + void init(){}; + void process(const Input&, const Output&) const final; +}; + +} // namespace eicrecon diff --git a/src/algorithms/digi/LGADPulseDigitizationConfig.h b/src/algorithms/digi/LGADPulseDigitizationConfig.h new file mode 100644 index 0000000000..498e907139 --- /dev/null +++ b/src/algorithms/digi/LGADPulseDigitizationConfig.h @@ -0,0 +1,26 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2024 Souvik Paul + +#pragma once + +#include + +namespace eicrecon { + +struct LGADPulseDigitizationConfig { + int adc_bit = 8; + int tdc_bit = 10; + // total number of TDC/ADC values + // Since digitizer starts at zero, max ADC value = adc_range - 1 + // Similar for TDC + int adc_range = pow(2, adc_bit); + int tdc_range = pow(2, tdc_bit); + + + double t_thres = 0.1 * adc_range; // TDC value = time when pulse exceed t_thres + // period of the sensor clock. Time internal to sensor will all be digitized to integer multiple + // of tInterval + double tMax = 25 * dd4hep::ns; // 25 ns is the period of 40MHz EIC clock +}; + +} // namespace eicrecon diff --git a/src/algorithms/digi/LGADPulseGeneration.cc b/src/algorithms/digi/LGADPulseGeneration.cc new file mode 100644 index 0000000000..63004822d7 --- /dev/null +++ b/src/algorithms/digi/LGADPulseGeneration.cc @@ -0,0 +1,107 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2024 Souvik Paul, Chun Yuen Tsang, Prithwish Tribedy +// Special Acknowledgement: Kolja Kauder +// +// Convert energy deposition into ADC pulses +// ADC pulses are assumed to follow the shape of landau function + +#include +#include +#include +#include +#include +#include +#include + +#include "LGADPulseGeneration.h" +#include "TMath.h" +#include "algorithms/digi/LGADPulseGenerationConfig.h" + +namespace eicrecon { + +LandauPulse::LandauPulse(double gain, double Vm, double sigma_analog, double adc_range) { + m_gain = gain; + m_sigma_analog = sigma_analog; + + // calculation of the extreme values for Landau distribution can be found on lin 514-520 of + // https://root.cern.ch/root/html524/src/TMath.cxx.html#fsokrB Landau reaches minimum for mpv = + // 0 and sigma = 1 at x = -0.22278 + const double x_when_landau_min = -0.22278; + const double landau_min = -gain * TMath::Landau(x_when_landau_min, 0, 1, kTRUE) / m_sigma_analog; + m_scalingFactor = 1. / Vm / landau_min * adc_range; +} + +double LandauPulse::Eval(double time, double hit_time, double charge) { + return charge * m_gain * TMath::Landau(time, hit_time, m_sigma_analog, kTRUE) * m_scalingFactor; +} + +void LGADPulseGeneration::_FillADCArray(AdcArray& adc_sum, double charge, double mpv_analog, int n_EICROC_cycle, dd4hep::rec::CellID cellID) const { + double Vm = m_cfg.Vm; + double t = 0; + double tMax = m_cfg.tMax; + double interval = m_cfg.tInterval; + int adc_range = m_cfg.adc_range; + int nBins = m_cfg.tdc_range; + + // amplitude has to be negative + // because voltage is negative + // fetch the corresponding array + auto& ADCs = adc_sum[cellID][n_EICROC_cycle]; + if (ADCs.size() == 0) + ADCs.resize(nBins, 0); + + // keep filling the array until added amplitude < ignore_thres + for (unsigned int j = 0; j <= ADCs.size(); ++j, t += interval) { + double amplitude = m_pulse -> Eval(t, mpv_analog, charge); + if(std::fabs(amplitude) > std::fabs(m_cfg.ignore_thres * adc_range/m_cfg.Vm)){ + if(j >= ADCs.size()) { + // pulse has to be saved in the next clock cycle + this -> _FillADCArray(adc_sum, charge, mpv_analog - tMax, n_EICROC_cycle+1, cellID); + } else ADCs[j] += amplitude; + } + } + +} + +void LGADPulseGeneration::process(const LGADPulseGeneration::Input& input, + const LGADPulseGeneration::Output& output) const { + const auto [simhits] = input; + auto [rawADCs] = output; + + int tdc_range = m_cfg.tdc_range; + + // signal sum + // NOTE: we take the cellID of the most energetic hit in this group so it is a real cellID from an + // MC hit + AdcArray adc_sum; + + for (const auto& hit : *simhits) { + auto cellID = hit.getCellID(); + + double time = hit.getTime() * dd4hep::ns; + double charge = hit.getEDep(); + // reduce computation power by not simulating low-charge hits + if (charge < m_cfg.ignore_thres) + continue; + + int n_EICROC_cycle = static_cast(std::floor(time/m_cfg.tMax + 1e-3)); + double time_in_cycle = time - n_EICROC_cycle * m_cfg.tMax; + double mpv_analog = time_in_cycle + m_cfg.risetime; + this -> _FillADCArray(adc_sum, charge, mpv_analog, n_EICROC_cycle, cellID); + } + + // convert vector of ADC values to RawTimeSeries + for (const auto& [cellID, nCycleADCs] : adc_sum) { + for (const auto& [nCycle, ADCs] : nCycleADCs) { + auto time_series = rawADCs->create(); + time_series.setCellID(cellID); + time_series.setTime(nCycle * m_cfg.tMax); + time_series.setCharge(1.); // placeholder. Don't know what to assign when there are two or more hits + time_series.setInterval(m_cfg.tInterval); + + for (const auto ADC : ADCs) + time_series.addToAdcCounts(ADC); + } + } +} // LGADPulseGeneration:process +} // namespace eicrecon diff --git a/src/algorithms/digi/LGADPulseGeneration.h b/src/algorithms/digi/LGADPulseGeneration.h new file mode 100644 index 0000000000..e78e05ecb3 --- /dev/null +++ b/src/algorithms/digi/LGADPulseGeneration.h @@ -0,0 +1,69 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2024 Souvik Paul, Chun Yuen Tsang, Prithwish Tribedy +// Special Acknowledgement: Kolja Kauder +// +// Convert energy deposition into ADC pulses +// ADC pulses are assumed to follow the shape of landau function + +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "algorithms/digi/LGADPulseGenerationConfig.h" +#include "algorithms/interfaces/WithPodConfig.h" + +namespace eicrecon { + +using LGADPulseGenerationAlgorithm = + algorithms::Algorithm, + algorithms::Output>; + +class PulseShape { +public: + // return voltage in units of ADC value + // ranges from 0 to adc_range + // charge is in unit GeV, it's the EDep in the detector + virtual double Eval(double time, double hit_time, double charge) = 0; + virtual ~PulseShape() {}; +}; + +// default pulse shape is Landau +class LandauPulse : public PulseShape { +public: + LandauPulse(double gain, double Vm, double sigma_analog, double adc_range); + double Eval(double time, double hit_time, double charge); +private: + double m_gain; + double m_sigma_analog; + double m_scalingFactor; +}; + + +class LGADPulseGeneration : public LGADPulseGenerationAlgorithm, + public WithPodConfig { +// The key pair is +typedef std::unordered_map>> AdcArray; + +public: + LGADPulseGeneration(std::string_view name, std::unique_ptr&& pulse) + : LGADPulseGenerationAlgorithm{name, {"RawHits"}, {"OutputPulses"}, {}}, m_pulse(std::move(pulse)) {} + virtual void init(){}; + void process(const Input&, const Output&) const final; + +protected: + double _Landau(double amp, double x, double mean, double std) const; + void _FillADCArray(AdcArray& adc_sum, double charge, double mpv_analog, + int n_EICROC_cycle, dd4hep::rec::CellID cellID) const; + std::unique_ptr m_pulse; +}; + +} // namespace eicrecon diff --git a/src/algorithms/digi/LGADPulseGenerationConfig.h b/src/algorithms/digi/LGADPulseGenerationConfig.h new file mode 100644 index 0000000000..e46d6b85a1 --- /dev/null +++ b/src/algorithms/digi/LGADPulseGenerationConfig.h @@ -0,0 +1,35 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2024 Souvik Paul + +#pragma once + +#include + +namespace eicrecon { + +struct LGADPulseGenerationConfig { + // Parameters of AC-LGAD signal generation + double gain = 113.755; + double risetime = 0.45 * dd4hep::ns; + double sigma_analog = 0.293951 * dd4hep::ns; + double Vm = -1e-4 * dd4hep::GeV; // Vm = voltage maximum. When EDep = 1e-4 GeV, voltage + // corresponds to ADC = adc_max + double ignore_thres = 0.001 * Vm; // If EDep below this value, digitization for the cell will be + // ignored. Speed up calculation + int adc_bit = 8; + int tdc_bit = 10; + + // total number of TDC/ADC values + // Since digitizer starts at zero, max ADC value = adc_range - 1 + // Similar for TDC + int adc_range = pow(2, adc_bit); + int tdc_range = pow(2, tdc_bit); + + // period of the sensor clock. Time internal to sensor will all be digitized to integer multiple + // of tInterval + double tMax = 25 * dd4hep::ns; // 25 ns is the period of 40MHz EIC clock + double tInterval = tMax / (tdc_range - 1); + +}; + +} // namespace eicrecon diff --git a/src/detectors/BTOF/BTOF.cc b/src/detectors/BTOF/BTOF.cc index 7c030ba7db..2a2a41ec24 100644 --- a/src/detectors/BTOF/BTOF.cc +++ b/src/detectors/BTOF/BTOF.cc @@ -18,90 +18,115 @@ #include "extensions/jana/JOmniFactoryGeneratorT.h" #include "factories/digi/SiliconTrackerDigi_factory.h" #include "factories/tracking/TrackerHitReconstruction_factory.h" +#include "factories/digi/BTOFChargeSharing_factory.h" +#include "factories/digi/LGADPulseGeneration_factory.h" +#include "factories/digi/LGADPulseDigitization_factory.h" #include "global/pid_lut/PIDLookup_factory.h" #include "services/geometry/dd4hep/DD4hep_service.h" extern "C" { -void InitPlugin(JApplication *app) { - InitJANAPlugin(app); +void InitPlugin(JApplication* app) { + InitJANAPlugin(app); - using namespace eicrecon; + using namespace eicrecon; - // Digitization - app->Add(new JOmniFactoryGeneratorT( + // Digitization + app->Add(new JOmniFactoryGeneratorT( + "TOFBarrelRawHits", + { + "TOFBarrelHits" + }, + { "TOFBarrelRawHits", - { - "TOFBarrelHits" - }, - { - "TOFBarrelRawHits", - "TOFBarrelRawHitAssociations" - }, - { - .threshold = 6.0 * dd4hep::keV, - .timeResolution = 0.025, // [ns] - }, - app - )); + "TOFBarrelRawHitAssociations" + }, + { + .threshold = 6.0 * dd4hep::keV, + .timeResolution = 0.025, // [ns] + }, + app + )); - // Convert raw digitized hits into hits with geometry info (ready for tracking) - app->Add(new JOmniFactoryGeneratorT( - "TOFBarrelRecHits", - {"TOFBarrelRawHits"}, // Input data collection tags - {"TOFBarrelRecHits"}, // Output data tag - { - .timeResolution = 10, - }, - app - )); // Hit reco default config for factories + // Convert raw digitized hits into hits with geometry info (ready for tracking) + app->Add(new JOmniFactoryGeneratorT( + "TOFBarrelRecHits", + {"TOFBarrelRawHits"}, // Input data collection tags + {"TOFBarrelRecHits"}, // Output data tag + { + .timeResolution = 10, + }, + app + )); // Hit reco default config for factories - int BarrelTOF_ID = 0; - try { - auto detector = app->GetService()->detector(); - BarrelTOF_ID = detector->constant("BarrelTOF_ID"); - } catch(const std::runtime_error&) { - // Nothing - } - PIDLookupConfig pid_cfg { - .filename="calibrations/tof.lut", - .system=BarrelTOF_ID, - .pdg_values={11, 211, 321, 2212}, - .charge_values={1}, - .momentum_edges={0.0, 0.3, 0.6, 0.9, 1.2, 1.5, 1.8, 2.1, 2.4, 2.7, 3.0, 3.3, 3.6, 3.9, 4.2, 4.5, 4.8, 5.1, 5.4, 5.7, 6.0}, - .polar_edges={2.50, 10.95, 19.40, 27.85, 36.30, 44.75, 53.20, 61.65, 70.10, 78.55, 87.00, 95.45, 103.90, 112.35, 120.80, 129.25, 137.70, 146.15, 154.60}, - .azimuthal_binning={0., 360., 360.}, // lower, upper, step - .momentum_bin_centers_in_lut=true, - .polar_bin_centers_in_lut=true, - }; + app->Add(new JOmniFactoryGeneratorT( + "BTOFChargeSharing", + {"TOFBarrelHits"}, + {"TOFBarrelSharedHits"}, + {}, + app + )); - app->Add(new JOmniFactoryGeneratorT( - "CombinedTOFTruthSeededLUTPID", - { + app->Add(new JOmniFactoryGeneratorT( + "BTOFPulseGeneration", + {"TOFBarrelSharedHits"}, + {"TOFBarrelPulse"}, + {}, + app + )); + + app->Add(new JOmniFactoryGeneratorT( + "BTOFPulseDigitization", + {"TOFBarrelPulse"}, + {"TOFBarrelADCTDC"}, + {}, + app + )); + + int BarrelTOF_ID = 0; + try { + auto detector = app->GetService()->detector(); + BarrelTOF_ID = detector->constant("BarrelTOF_ID"); + } catch (const std::runtime_error&) { + // Nothing + } + PIDLookupConfig pid_cfg{ + .filename = "calibrations/tof.lut", + .system = BarrelTOF_ID, + .pdg_values = {11, 211, 321, 2212}, + .charge_values = {1}, + .momentum_edges = {0.0, 0.3, 0.6, 0.9, 1.2, 1.5, 1.8, 2.1, 2.4, 2.7, 3.0, + 3.3, 3.6, 3.9, 4.2, 4.5, 4.8, 5.1, 5.4, 5.7, 6.0}, + .polar_edges = {2.50, 10.95, 19.40, 27.85, 36.30, 44.75, 53.20, 61.65, 70.10, 78.55, 87.00, + 95.45, 103.90, 112.35, 120.80, 129.25, 137.70, 146.15, 154.60}, + .azimuthal_binning = {0., 360., 360.}, // lower, upper, step + .momentum_bin_centers_in_lut = true, + .polar_bin_centers_in_lut = true, + }; + + app->Add(new JOmniFactoryGeneratorT( + "CombinedTOFTruthSeededLUTPID", + { "ReconstructedTruthSeededChargedWithPFRICHPIDParticles", "ReconstructedTruthSeededChargedWithPFRICHPIDParticleAssociations", - }, - { + }, + { "ReconstructedTruthSeededChargedWithPFRICHTOFPIDParticles", "ReconstructedTruthSeededChargedWithPFRICHTOFPIDParticleAssociations", "CombinedTOFTruthSeededParticleIDs", - }, - pid_cfg, - app - )); + }, + pid_cfg, app)); - app->Add(new JOmniFactoryGeneratorT( - "CombinedTOFLUTPID", - { + app->Add(new JOmniFactoryGeneratorT( + "CombinedTOFLUTPID", + { "ReconstructedChargedWithPFRICHPIDParticles", "ReconstructedChargedWithPFRICHPIDParticleAssociations", - }, - { + }, + { "ReconstructedChargedWithPFRICHTOFPIDParticles", "ReconstructedChargedWithPFRICHTOFPIDParticleAssociations", "CombinedTOFParticleIDs", - }, - pid_cfg, - app - )); + }, + pid_cfg, app)); } } // extern "C" diff --git a/src/factories/digi/BTOFChargeSharing_factory.h b/src/factories/digi/BTOFChargeSharing_factory.h new file mode 100644 index 0000000000..801d7454ef --- /dev/null +++ b/src/factories/digi/BTOFChargeSharing_factory.h @@ -0,0 +1,43 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2024 Chun Yuen Tsang + +#pragma once + +#include "extensions/jana/JOmniFactory.h" + +#include "algorithms/digi/BTOFChargeSharing.h" +#include + +namespace eicrecon { + +class BTOFChargeSharing_factory : public JOmniFactory { +public: + using AlgoT = eicrecon::BTOFChargeSharing; + +private: + std::unique_ptr m_algo; + + PodioInput m_in_sim_track{this}; + + PodioOutput m_out_reco_particles{this}; + + ParameterRef m_sigma_sharingx{this, "sigmaSharingX", config().sigma_sharingx}; + ParameterRef m_sigma_sharingy{this, "sigmaSharingY", config().sigma_sharingy}; + + Service m_algorithmsInit{this}; + +public: + void Configure() { + m_algo = std::make_unique(GetPrefix()); + m_algo->level(static_cast(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_sim_track()}, {m_out_reco_particles().get()}); + } +}; +} // namespace eicrecon diff --git a/src/factories/digi/LGADPulseDigitization_factory.h b/src/factories/digi/LGADPulseDigitization_factory.h new file mode 100644 index 0000000000..e3ef225cdd --- /dev/null +++ b/src/factories/digi/LGADPulseDigitization_factory.h @@ -0,0 +1,43 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2024 Chun Yuen Tsang + +#pragma once + +#include "extensions/jana/JOmniFactory.h" + +#include "algorithms/digi/LGADPulseDigitization.h" +#include + +namespace eicrecon { + +class LGADPulseDigitization_factory + : public JOmniFactory { +public: + using AlgoT = eicrecon::LGADPulseDigitization; + +private: + std::unique_ptr m_algo; + + PodioInput m_in_sim_track{this}; + + PodioOutput m_out_reco_particles{this}; + + ParameterRef m_t_thres{this, "tThreshold", config().t_thres}; + + Service m_algorithmsInit{this}; + +public: + void Configure() { + m_algo = std::make_unique(GetPrefix()); + m_algo->level(static_cast(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_sim_track()}, {m_out_reco_particles().get()}); + } +}; +} // namespace eicrecon diff --git a/src/factories/digi/LGADPulseGeneration_factory.h b/src/factories/digi/LGADPulseGeneration_factory.h new file mode 100644 index 0000000000..4659282bb2 --- /dev/null +++ b/src/factories/digi/LGADPulseGeneration_factory.h @@ -0,0 +1,50 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2024 Chun Yuen Tsang + +#pragma once + +#include "extensions/jana/JOmniFactory.h" + +#include "algorithms/digi/LGADPulseGeneration.h" +#include + +namespace eicrecon { + +class LGADPulseGeneration_factory + : public JOmniFactory { +public: + using AlgoT = eicrecon::LGADPulseGeneration; + +private: + std::unique_ptr m_algo; + + PodioInput m_in_sim_track{this}; + + PodioOutput m_out_reco_particles{this}; + + ParameterRef m_Vm{this, "Vm", config().Vm}; + ParameterRef m_tMax{this, "tMax", config().tMax}; + ParameterRef m_adc_range{this, "adcRange", config().adc_range}; + ParameterRef m_tdc_range{this, "tdcRange", config().tdc_range}; + ParameterRef m_ignore_thres{this, "ignoreThreshold", config().ignore_thres}; + + Service m_algorithmsInit{this}; + +public: + void Configure() { + std::unique_ptr landau = std::make_unique(config().gain, config().Vm, + config().sigma_analog, config().adc_range); + m_algo = std::make_unique(GetPrefix(), std::move(landau)); + m_algo->level(static_cast(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_sim_track()}, {m_out_reco_particles().get()}); + } +}; + +} // namespace eicrecon diff --git a/src/services/io/podio/JEventProcessorPODIO.cc b/src/services/io/podio/JEventProcessorPODIO.cc index 2575c326a9..6df43a4c5a 100644 --- a/src/services/io/podio/JEventProcessorPODIO.cc +++ b/src/services/io/podio/JEventProcessorPODIO.cc @@ -92,6 +92,7 @@ JEventProcessorPODIO::JEventProcessorPODIO() { "TOFEndcapRawHits", "TOFBarrelHits", + "TOFBarrelADCTDC", "TOFEndcapHits", "TOFBarrelRawHitAssociations", diff --git a/src/tests/algorithms_test/CMakeLists.txt b/src/tests/algorithms_test/CMakeLists.txt index 57ecd29c88..add56f0616 100644 --- a/src/tests/algorithms_test/CMakeLists.txt +++ b/src/tests/algorithms_test/CMakeLists.txt @@ -11,6 +11,8 @@ add_executable( calorimetry_CalorimeterHitDigi.cc calorimetry_CalorimeterClusterRecoCoG.cc calorimetry_HEXPLIT.cc + digi_LGADPulseGeneration.cc + digi_LGADPulseDigitization.cc pid_MergeTracks.cc pid_MergeParticleID.cc pid_lut_PIDLookup.cc diff --git a/src/tests/algorithms_test/algorithmsInit.cc b/src/tests/algorithms_test/algorithmsInit.cc index 357e928001..bb673c717e 100644 --- a/src/tests/algorithms_test/algorithmsInit.cc +++ b/src/tests/algorithms_test/algorithmsInit.cc @@ -44,6 +44,15 @@ class algorithmsInitListener : public Catch::EventListenerBase { detector->add(id_desc_tracker); detector->add(readoutTracker); + dd4hep::Readout readoutLGAD(std::string("MockLGADHits")); + dd4hep::IDDescriptor id_desc_LGAD("MockLGADHits", "system:8,layer:4,module:12,sensor:10,x:40:-8,y:-16"); + //Create segmentation with 1x1 mm pixels + dd4hep::Segmentation segmentation_LGAD("CartesianGridXY","LGADHitsSeg", id_desc_tracker.decoder()); + readoutLGAD.setIDDescriptor(id_desc_LGAD); + readoutLGAD.setSegmentation(segmentation_LGAD); + detector->add(id_desc_LGAD); + detector->add(readoutLGAD); + m_detector = std::move(detector); auto& serviceSvc = algorithms::ServiceSvc::instance(); diff --git a/src/tests/algorithms_test/digi_LGADPulseDigitization.cc b/src/tests/algorithms_test/digi_LGADPulseDigitization.cc new file mode 100644 index 0000000000..5f11c0858c --- /dev/null +++ b/src/tests/algorithms_test/digi_LGADPulseDigitization.cc @@ -0,0 +1,145 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2024, Chun Yuen Tsang, Prithwish Tribedy + +#include +#include +#include +#include +#include +#include +#include // for AssertionHandler, operator""_catch_sr, StringRef, REQUIRE, operator<, operator==, operator>, TEST_CASE +#include +#include +#include // for level_enum +#include // for logger +#include // for default_logger +#include +#include +#include // for allocator, unique_ptr, make_unique, shared_ptr, __shared_ptr_access +#include +#include + +#include "algorithms/digi/LGADPulseDigitization.h" +#include "algorithms/digi/LGADPulseDigitizationConfig.h" + +TEST_CASE("the LGAD charge sharing algorithm runs", "[LGADPulseDigitization]") { + const float EPSILON = 1e-5; + + eicrecon::LGADPulseDigitization algo("LGADPulseDigitization"); + + std::shared_ptr logger = spdlog::default_logger()->clone("LGADPulseDigitization"); + logger->set_level(spdlog::level::trace); + + eicrecon::LGADPulseDigitizationConfig cfg; + + auto detector = algorithms::GeoSvc::instance().detector(); + auto id_desc = detector->readout("MockLGADHits").idSpec(); + + cfg.tdc_bit = 8; + cfg.adc_bit = 7; + cfg.tMax = 10 * dd4hep::ns; + cfg.tdc_range = pow(2, cfg.tdc_bit); + cfg.adc_range = pow(2, cfg.adc_bit); + cfg.t_thres = cfg.adc_range * 0.1; + + // check if max pulse height is linearly proportional to the initial Edep + algo.applyConfig(cfg); + algo.init(); + + SECTION("TDC vs analytic solution scan") { + logger->info("Begin TDC vs analytic solution scan"); + + for(double time = -cfg.tMax; time <= cfg.tMax; time += cfg.tMax) { + if(time < 0) logger->info("Generation pulse at the negative time"); + else if(time == 0) logger->info("Generation pulse at the first EICROC cycle"); + else logger->info("Generation pulse at the second EICROC cycle"); + + // test pulse with gaussian shape + for (double tdc_frac = 0.4; tdc_frac < 1; tdc_frac += 0.1) { + edm4hep::RawTimeSeriesCollection time_series_coll; + auto rawhits_coll = std::make_unique(); + + auto pulse = time_series_coll.create(); + auto cellID = + id_desc.encode({{"system", 0}, {"module", 0}, {"sensor", 1}, {"x", 1}, {"y", 1}}); + + pulse.setCellID(cellID); + pulse.setCharge(1.); // placeholder + pulse.setTime(time); // placeholder + pulse.setInterval(1); + + int test_peak_TDC = static_cast(tdc_frac * cfg.tdc_range); + int test_peak = static_cast(0.7 * cfg.adc_range); + int test_peak_sigma = static_cast(0.1 * cfg.tdc_range); + + for (int i = 0; i < cfg.tdc_range; ++i) { + int ADC = + -test_peak * + TMath::Exp(-0.5 * pow((i - test_peak_TDC) / static_cast(test_peak_sigma), 2)); + pulse.addToAdcCounts(ADC); + } + + // calculate analytically when the pulse passes t_thres + // ADC = amp*exp(-(TDC - mean)^2/(2sigma^2)) + // TDC = mean - (-2*sigma^2*ln(ADC/amp))^0.5 + int analytic_TDC = ceil(test_peak_TDC - sqrt(-2 * pow(test_peak_sigma, 2) * + TMath::Log(cfg.adc_range * 0.1 / test_peak))); + + // Constructing input and output as per the algorithm's expected signature + auto input = std::make_tuple(&time_series_coll); + auto output = std::make_tuple(rawhits_coll.get()); + + algo.process(input, output); + + REQUIRE(rawhits_coll->size() == 1); + auto hit = (*rawhits_coll)[0]; + REQUIRE(hit.getCellID() == cellID); + REQUIRE(hit.getCharge() == test_peak); + if(time < 0) REQUIRE(hit.getTimeStamp() == analytic_TDC - cfg.tdc_range); + else if(time == 0) REQUIRE(hit.getTimeStamp() == analytic_TDC); + else REQUIRE(hit.getTimeStamp() == analytic_TDC + cfg.tdc_range); + } + } + } + + SECTION("ADC scan") { + logger->info("Begin ADC scan"); + + // test pulse with gaussian shape + for (double adc_frac = 0.4; adc_frac < 1; adc_frac += 0.1) { + edm4hep::RawTimeSeriesCollection time_series_coll; + auto rawhits_coll = std::make_unique(); + + auto pulse = time_series_coll.create(); + auto cellID = + id_desc.encode({{"system", 0}, {"module", 0}, {"sensor", 1}, {"x", 1}, {"y", 1}}); + + pulse.setCellID(cellID); + pulse.setCharge(1.); // placeholder + pulse.setTime(0.); // placeholder + pulse.setInterval(1); + + int test_peak_TDC = static_cast(0.5 * cfg.tdc_range); + int test_peak = static_cast(adc_frac * cfg.adc_range); + int test_peak_sigma = static_cast(0.1 * cfg.tdc_range); + + for (int i = 0; i < cfg.tdc_range; ++i) { + int ADC = + -test_peak * + TMath::Exp(-0.5 * pow((i - test_peak_TDC) / static_cast(test_peak_sigma), 2)); + pulse.addToAdcCounts(ADC); + } + + // Constructing input and output as per the algorithm's expected signature + auto input = std::make_tuple(&time_series_coll); + auto output = std::make_tuple(rawhits_coll.get()); + + algo.process(input, output); + + REQUIRE(rawhits_coll->size() == 1); + auto hit = (*rawhits_coll)[0]; + REQUIRE(hit.getCellID() == cellID); + REQUIRE(hit.getCharge() == test_peak); + } + } +} diff --git a/src/tests/algorithms_test/digi_LGADPulseGeneration.cc b/src/tests/algorithms_test/digi_LGADPulseGeneration.cc new file mode 100644 index 0000000000..e4ff7d8f04 --- /dev/null +++ b/src/tests/algorithms_test/digi_LGADPulseGeneration.cc @@ -0,0 +1,171 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2024, Chun Yuen Tsang, Prithwish Tribedy + +#include +#include +#include +#include +#include +#include // for AssertionHandler, operator""_catch_sr, StringRef, REQUIRE, operator<, operator==, operator>, TEST_CASE +#include +#include +#include +#include +#include // for level_enum +#include // for logger +#include // for default_logger +#include +#include +#include +#include +#include // for allocator, unique_ptr, make_unique, shared_ptr, __shared_ptr_access +#include +#include +#include + +#include "TF1.h" +#include "TGraphErrors.h" +#include "algorithms/digi/LGADPulseGeneration.h" +#include "algorithms/digi/LGADPulseGenerationConfig.h" + +TEST_CASE("the LGAD charge sharing algorithm runs", "[LGADPulseGeneration]") { + const float EPSILON = 1e-5; + eicrecon::LGADPulseGenerationConfig cfg; + cfg.gain = 113; + cfg.Vm = 1e-4 * dd4hep::GeV; + cfg.ignore_thres = 1e-4 / 5; + cfg.sigma_analog = 0.293951 * dd4hep::ns; + cfg.adc_bit = 8; + cfg.adc_range = pow(2, cfg.adc_bit); + + std::unique_ptr landau = std::make_unique(cfg.gain, cfg.Vm, + cfg.sigma_analog, cfg.adc_range); + + eicrecon::LGADPulseGeneration algo("LGADPulseGeneration", std::move(landau)); + + std::shared_ptr logger = spdlog::default_logger()->clone("LGADPulseGeneration"); + logger->set_level(spdlog::level::trace); + + + auto detector = algorithms::GeoSvc::instance().detector(); + auto id_desc = detector->readout("MockLGADHits").idSpec(); + + SECTION("Pulse height linearlity test") { + // check if max pulse height is linearly proportional to the initial Edep + algo.applyConfig(cfg); + algo.init(); + + TGraphErrors graph; + for (double edep = 0; edep <= 1e-4; edep += 1e-4 / 9) { + edm4hep::SimTrackerHitCollection rawhits_coll; + auto time_series_coll = std::make_unique(); + + auto hit = rawhits_coll.create(); + auto cellID = + id_desc.encode({{"system", 0}, {"module", 0}, {"sensor", 1}, {"x", 1}, {"y", 1}}); + + hit.setCellID(cellID); + hit.setEDep( + edep); // in GeV. Since Vm = 1e-4*gain, EDep = 1e-4 GeV corresponds to ADC = max_adc + hit.setTime(1.5 * dd4hep::ns); // in ns + + // Constructing input and output as per the algorithm's expected signature + auto input = std::make_tuple(&rawhits_coll); + auto output = std::make_tuple(time_series_coll.get()); + + algo.process(input, output); + + if (edep < 1e-4 / 5) + REQUIRE(time_series_coll->size() == 0); + else { + REQUIRE(time_series_coll->size() == 1); + auto min_adc = std::numeric_limits::max(); + for(const auto& pulse : (*time_series_coll)) { + REQUIRE(pulse.getCellID() == cellID); + auto adcs = pulse.getAdcCounts(); + for (const auto adc : adcs) + min_adc = std::min(min_adc, adc); + } + int npt = graph.GetN(); + graph.SetPoint(npt, edep, min_adc); + graph.SetPointError(npt, 0, 0.5); + // make sure when energy deposition = Vm, ADC reaches max value + if (edep == 1e-4) + REQUIRE(min_adc == -cfg.adc_range + 1); + } + } + + // test linearlity + TF1 tf1("tf1", "pol1", 0, 1e-4); + graph.Fit(&tf1, "R0"); + // slope can't be consistent with zero + REQUIRE(!(tf1.GetParameter(1) - tf1.GetParError(1) < 0 && 0 < tf1.GetParameter(1) + tf1.GetParError(1) )); + double chi2_dof = tf1.GetChisquare() / tf1.GetNDF(); + logger->info("Chi-square/dof value for Edep vs min-adc = {}", chi2_dof); + REQUIRE(chi2_dof < 2); + } + + SECTION("Pulse timing linearlity test") { + // check if max pulse height is linearly proportional to the initial Edep + algo.applyConfig(cfg); + algo.init(); + + TGraphErrors graph; + std::vector times; + + // test within the same EICROC cycle + for (double time = 0; time < 12; time += 1.) times.push_back(time); + // test multiple EICROC cycle + for (double time = 10; time < 101; time += 25.) times.push_back(time); + // test negative time + times.push_back(-10); + + for (double time : times) { + edm4hep::SimTrackerHitCollection rawhits_coll; + auto time_series_coll = std::make_unique(); + + auto hit = rawhits_coll.create(); + auto cellID = + id_desc.encode({{"system", 0}, {"module", 0}, {"sensor", 1}, {"x", 1}, {"y", 1}}); + + hit.setCellID(cellID); + hit.setEDep( + 0.5e-4 * dd4hep::GeV); // in GeV. Since Vm = 1e-4*gain, EDep = 1e-4 GeV corresponds to ADC = max_adc + hit.setTime(time); // in ns + + // Constructing input and output as per the algorithm's expected signature + auto input = std::make_tuple(&rawhits_coll); + auto output = std::make_tuple(time_series_coll.get()); + + algo.process(input, output); + + REQUIRE(time_series_coll->size() == 1); + auto min_adc = std::numeric_limits::max(); + int time_bin = 0; + for(const auto& pulse: *time_series_coll) { + //auto pulse = (*time_series_coll)[0]; + REQUIRE(pulse.getCellID() == cellID); + + auto adcs = pulse.getAdcCounts(); + for (unsigned int i = 0; i < adcs.size(); ++i) { + auto adc = adcs[i]; + if (adc < min_adc) + time_bin = i + pulse.getTime()/cfg.tMax*cfg.tdc_range; + min_adc = std::min(min_adc, adc); + } + } + int npt = graph.GetN(); + graph.SetPoint(npt, time, time_bin); + graph.SetPointError(npt, 0, 0.5); + } + + // test linearlity + TF1 tf1("tf1", "pol1", *std::min_element(times.begin(), times.end()), *std::max_element(times.begin(), times.end())); + graph.Fit(&tf1, "R0"); + // slope can't be consistent with zero + REQUIRE(!(tf1.GetParameter(1) - tf1.GetParError(1) < 0 && 0 < tf1.GetParameter(1) + tf1.GetParError(1) )); + double chi2_dof = tf1.GetChisquare() / tf1.GetNDF(); + logger->info("Chi-square/dof value for time vs TDC-bin = {}", chi2_dof); + REQUIRE(chi2_dof < 2); + } +}