diff --git a/src/algorithms/digi/BTOFChargeSharing.cc b/src/algorithms/digi/BTOFChargeSharing.cc new file mode 100644 index 0000000000..5704098347 --- /dev/null +++ b/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 +#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 "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> 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..7757bd17cb --- /dev/null +++ b/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 +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "DD4hep/Detector.h" +#include "algorithms/digi/TOFHitDigiConfig.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/TOFHitDigiConfig.h b/src/algorithms/digi/TOFHitDigiConfig.h new file mode 100644 index 0000000000..51697b9707 --- /dev/null +++ b/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 + +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 diff --git a/src/algorithms/digi/TOFPulseDigitization.cc b/src/algorithms/digi/TOFPulseDigitization.cc new file mode 100644 index 0000000000..550f463f46 --- /dev/null +++ b/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 +#include +#include +#include +#include +#include +#include + +#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::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(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(adc_range), round(-V)); + // only store valid hits + if (tdc < std::numeric_limits::max()) + rawhits->create(pulse.getCellID(), adc, tdc); + //----------------------------------------------------------- + + } +} // TOFPulseDigitization:process +} // namespace eicrecon diff --git a/src/algorithms/digi/TOFPulseDigitization.h b/src/algorithms/digi/TOFPulseDigitization.h new file mode 100644 index 0000000000..cf99a8eca4 --- /dev/null +++ b/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 +#include +#include +#include +#include + +#include "algorithms/digi/TOFHitDigiConfig.h" +#include "algorithms/interfaces/WithPodConfig.h" + +namespace eicrecon { + +using TOFPulseDigitizationAlgorithm = + algorithms::Algorithm, + algorithms::Output>; + +class TOFPulseDigitization : public TOFPulseDigitizationAlgorithm, + public WithPodConfig { + +public: + TOFPulseDigitization(std::string_view name) : + TOFPulseDigitizationAlgorithm{name, + {"TOFBarrelPulse"}, + {"TOFBarrelADCTDC"}, + {}} {} + void init() {}; + void process(const Input&, const Output&) const final; +}; + +} // namespace eicrecon diff --git a/src/algorithms/digi/TOFPulseGeneration.cc b/src/algorithms/digi/TOFPulseGeneration.cc new file mode 100644 index 0000000000..aca923fca5 --- /dev/null +++ b/src/algorithms/digi/TOFPulseGeneration.cc @@ -0,0 +1,90 @@ +// 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 "TMath.h" +#include "TOFPulseGeneration.h" +#include "algorithms/digi/TOFHitDigiConfig.h" + + +namespace eicrecon { + +double TOFPulseGeneration::_Landau(double x, double mean, double std) const { + double C = -113.755; + return C * TMath::Landau(x, mean, std, kTRUE); +} + +void TOFPulseGeneration::process(const TOFPulseGeneration::Input& input, + const TOFPulseGeneration::Output& output) const { + const auto [simhits] = input; + auto [rawADCs] = output; + + double Vm = m_cfg.Vm; + double tMin = m_cfg.tMin; + double tMax = m_cfg.tMax; + int adc_range = m_cfg.adc_range; + int tdc_range = m_cfg.tdc_range; + int nBins = 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 + std::unordered_map> adc_sum; + double interval = (tMax - tMin) / (nBins - 1); + + for (const auto& hit : *simhits) { + auto cellID = hit.getCellID(); + double sum_charge = 0.0; + double mpv_analog = 0.0; + + double time = hit.getTime(); + double charge = hit.getEDep(); + // reduce computation power by not simulating low-charge hits + if(charge < m_cfg.ignore_thres) continue; + + auto& ADCs = adc_sum[cellID]; + if(ADCs.size() == 0) ADCs.resize(nBins, 0); + + mpv_analog = time + m_cfg.risetime; + + double landau_min = this -> _Landau(0, mpv_analog, m_cfg.sigma_analog); + // find minimum of the landau function + // minimum = peak because prefactor is negative + for (int j = 0; j < nBins; ++j) { + double x = tMin + j * interval; + landau_min = std::min(landau_min, this -> _Landau(x, mpv_analog, m_cfg.sigma_analog)); + } + + double scalingFactor = -1. / Vm / m_cfg.gain / landau_min * adc_range; + + for (int j = 0; j < nBins; ++j) { + double x = tMin + j * interval; + double y = -1 * charge * m_cfg.gain * this -> _Landau(x, mpv_analog, m_cfg.sigma_analog) * scalingFactor; + ADCs[j] += y;; + } + + } + + // convert vector of ADC values to RawTimeSeries + for(const auto &[cellID, ADCs] : adc_sum) { + auto time_series = rawADCs -> create(); + time_series.setCellID(cellID); + time_series.setTime(1.); // placeholder. Don't know what to assign when there are two or more hits + time_series.setCharge(1.); // placeholder. Don't know what to assign when there are two or more hits + time_series.setInterval(interval); + + for(const auto ADC : ADCs) + time_series.addToAdcCounts(ADC); + } +} // TOFPulseGeneration:process +} // namespace eicrecon diff --git a/src/algorithms/digi/TOFPulseGeneration.h b/src/algorithms/digi/TOFPulseGeneration.h new file mode 100644 index 0000000000..20f9c910cb --- /dev/null +++ b/src/algorithms/digi/TOFPulseGeneration.h @@ -0,0 +1,41 @@ +// 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 "algorithms/digi/TOFHitDigiConfig.h" +#include "algorithms/interfaces/WithPodConfig.h" + +namespace eicrecon { + +using TOFPulseGenerationAlgorithm = + algorithms::Algorithm, + algorithms::Output>; + +class TOFPulseGeneration : public TOFPulseGenerationAlgorithm, + public WithPodConfig { + +public: + TOFPulseGeneration(std::string_view name) + : TOFPulseGenerationAlgorithm{name, + {"TOFBarrelSharedHits"}, + {"TOFBarrelPulse"}, + {}} {} + void init() {}; + void process(const Input&, const Output&) const final; + +protected: + double _Landau(double x, double mean, double std) const; +}; + +} // namespace eicrecon diff --git a/src/detectors/BTOF/BTOF.cc b/src/detectors/BTOF/BTOF.cc index 7c030ba7db..71e86826fa 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/TOFPulseGeneration_factory.h" +#include "factories/digi/TOFPulseDigitization_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..1927a0a93d --- /dev/null +++ b/src/factories/digi/BTOFChargeSharing_factory.h @@ -0,0 +1,42 @@ +// 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/TOFPulseDigitization_factory.h b/src/factories/digi/TOFPulseDigitization_factory.h new file mode 100644 index 0000000000..510ac3f183 --- /dev/null +++ b/src/factories/digi/TOFPulseDigitization_factory.h @@ -0,0 +1,41 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2024 Chun Yuen Tsang + +#pragma once + +#include "extensions/jana/JOmniFactory.h" + +#include "algorithms/digi/TOFPulseDigitization.h" +#include + +namespace eicrecon { + +class TOFPulseDigitization_factory : public JOmniFactory { +public: + using AlgoT = eicrecon::TOFPulseDigitization; +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/TOFPulseGeneration_factory.h b/src/factories/digi/TOFPulseGeneration_factory.h new file mode 100644 index 0000000000..cd3956f46b --- /dev/null +++ b/src/factories/digi/TOFPulseGeneration_factory.h @@ -0,0 +1,47 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2024 Chun Yuen Tsang + +#pragma once + +#include "extensions/jana/JOmniFactory.h" + +#include "algorithms/digi/TOFPulseGeneration.h" +#include + +namespace eicrecon { + +class TOFPulseGeneration_factory : public JOmniFactory { +public: + using AlgoT = eicrecon::TOFPulseGeneration; +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_tMin{this, "tMin", config().tMin}; + 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() { + 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/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..b6b2531ca5 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_TOFPulseGeneration.cc + digi_TOFPulseDigitization.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..1d5a839e28 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 readoutTOF(std::string("MockTOFHits")); + dd4hep::IDDescriptor id_desc_TOF("MockTOFHits", "system:8,layer:4,module:12,sensor:10,x:40:-8,y:-16"); + //Create segmentation with 1x1 mm pixels + dd4hep::Segmentation segmentation_TOF("CartesianGridXY","TOFHitsSeg", id_desc_tracker.decoder()); + readoutTOF.setIDDescriptor(id_desc_TOF); + readoutTOF.setSegmentation(segmentation_TOF); + detector->add(id_desc_TOF); + detector->add(readoutTOF); + m_detector = std::move(detector); auto& serviceSvc = algorithms::ServiceSvc::instance(); diff --git a/src/tests/algorithms_test/digi_TOFPulseDigitization.cc b/src/tests/algorithms_test/digi_TOFPulseDigitization.cc new file mode 100644 index 0000000000..f431086530 --- /dev/null +++ b/src/tests/algorithms_test/digi_TOFPulseDigitization.cc @@ -0,0 +1,150 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2024, Chun Yuen Tsang, Prithwish Tribedy + +#include // for AssertionHandler, operator""_catch_sr, StringRef, REQUIRE, operator<, operator==, operator>, TEST_CASE +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include // for allocator, unique_ptr, make_unique, shared_ptr, __shared_ptr_access +#include // for level_enum +#include // for logger +#include // for default_logger +#include +#include + +#include "algorithms/digi/TOFHitDigiConfig.h" +#include "algorithms/digi/TOFPulseDigitization.h" +#include +#include +#include +#include +#include + +#include "TF1.h" +#include "TGraphErrors.h" + +TEST_CASE("the BTOF charge sharing algorithm runs", "[TOFPulseDigitization]") { + const float EPSILON = 1e-5; + + eicrecon::TOFPulseDigitization algo("TOFPulseDigitization"); + + std::shared_ptr logger = spdlog::default_logger()->clone("TOFPulseDigitization"); + logger->set_level(spdlog::level::trace); + + eicrecon::TOFHitDigiConfig cfg; + cfg.readout = "MockTOFHits"; + + auto detector = algorithms::GeoSvc::instance().detector(); + auto id_desc = detector->readout(cfg.readout).idSpec(); + + cfg.gain = 10; + cfg.Vm = -1e-4; + cfg.ignore_thres = 1e-4 / 5; + cfg.t_thres = cfg.Vm * 0.1; + cfg.tdc_bit = 8; + cfg.adc_bit = 7; + cfg.tdc_range = pow(2, cfg.tdc_bit); + cfg.adc_range = pow(2, cfg.adc_bit); + + // check if max pulse height is linearly proportional to the initial Edep + algo.applyConfig(cfg); + algo.init(); + + SECTION("TDC vs analytic solution scan") { + + // 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(1.); // 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); + REQUIRE(hit.getTimeStamp() == analytic_TDC); + } + } + + SECTION("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(1.); // 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_TOFPulseGeneration.cc b/src/tests/algorithms_test/digi_TOFPulseGeneration.cc new file mode 100644 index 0000000000..b09bae1bb2 --- /dev/null +++ b/src/tests/algorithms_test/digi_TOFPulseGeneration.cc @@ -0,0 +1,158 @@ +// SPDX-License-Identifier: LGPL-3.0-or-later +// Copyright (C) 2024, Chun Yuen Tsang, Prithwish Tribedy + +#include // for AssertionHandler, operator""_catch_sr, StringRef, REQUIRE, operator<, operator==, operator>, TEST_CASE +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include // for allocator, unique_ptr, make_unique, shared_ptr, __shared_ptr_access +#include // for level_enum +#include // for logger +#include // for default_logger +#include +#include + +#include "algorithms/digi/TOFHitDigiConfig.h" +#include "algorithms/digi/TOFPulseGeneration.h" +#include +#include +#include +#include + +#include "TF1.h" +#include "TGraphErrors.h" + +TEST_CASE("the BTOF charge sharing algorithm runs", "[TOFPulseGeneration]") { + const float EPSILON = 1e-5; + + eicrecon::TOFPulseGeneration algo("TOFPulseGeneration"); + + std::shared_ptr logger = spdlog::default_logger()->clone("TOFPulseGeneration"); + logger->set_level(spdlog::level::trace); + + eicrecon::TOFHitDigiConfig cfg; + cfg.readout = "MockTOFHits"; + + auto detector = algorithms::GeoSvc::instance().detector(); + auto id_desc = detector->readout(cfg.readout).idSpec(); + + cfg.gain = 10; + cfg.Vm = -1e-4; + cfg.ignore_thres = 1e-4 / 5; + + 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); // 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 pulse = (*time_series_coll)[0]; + REQUIRE(pulse.getCellID() == cellID); + + auto adcs = pulse.getAdcCounts(); + auto min_adc = std::numeric_limits::max(); + 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"); + 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; + for (double time = 0; time < 12; time += 1.) { + 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); // 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 pulse = (*time_series_coll)[0]; + REQUIRE(pulse.getCellID() == cellID); + + auto adcs = pulse.getAdcCounts(); + auto min_adc = std::numeric_limits::max(); + unsigned int time_bin = 0; + for (unsigned int i = 0; i < adcs.size(); ++i) { + auto adc = adcs[i]; + if (adc < min_adc) + time_bin = i; + 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", 0, 12); + graph.Fit(&tf1, "R0"); + double chi2_dof = tf1.GetChisquare() / tf1.GetNDF(); + logger->info("Chi-square/dof value for time vs TDC-bin = {}", chi2_dof); + REQUIRE(chi2_dof < 2); + } +}