tests/algorithms_test/digi_EICROCDigitization.cc

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2024, Chun Yuen Tsang, Prithwish Tribedy
0003
0004 #include <DD4hep/Detector.h>
0005 #include <DD4hep/IDDescriptor.h>
0006 #include <DD4hep/Readout.h>
0007 #include <Evaluator/DD4hepUnits.h>
0008 #include <TMath.h>
0009 #include <algorithms/geo.h>
0010 #include <catch2/catch_test_macros.hpp> // for AssertionHandler, operator""_catch_sr, StringRef, REQUIRE, operator<, operator==, operator>, TEST_CASE
0011 #include <edm4eic/RawTrackerHitCollection.h>
0012 #include <edm4hep/RawTimeSeriesCollection.h>
0013 #include <spdlog/common.h> // for level_enum
0014 #include <spdlog/logger.h> // for logger
0015 #include <spdlog/spdlog.h> // for default_logger
0016 #include <cmath>
0017 #include <gsl/pointers>
0018 #include <memory> // for allocator, unique_ptr, make_unique, shared_ptr, __shared_ptr_access
0019 #include <string>
0020 #include <tuple>
0021 #include <utility>
0022 #include <vector>
0023
0024 #include "algorithms/digi/EICROCDigitization.h"
0025 #include "algorithms/digi/EICROCDigitizationConfig.h"
0026
0027 TEST_CASE("the Silicon charge sharing algorithm runs", "[EICROCDigitization]") {
0028   eicrecon::EICROCDigitization algo("EICROCDigitization");
0029
0030   std::shared_ptr<spdlog::logger> logger = spdlog::default_logger()->clone("EICROCDigitization");
0031   logger->set_level(spdlog::level::trace);
0032
0033   eicrecon::EICROCDigitizationConfig cfg;
0034
0035   auto detector = algorithms::GeoSvc::instance().detector();
0036   auto id_desc  = detector->readout("MockSiliconHits").idSpec();
0037
0038   cfg.tdc_bit   = 8;
0039   cfg.adc_bit   = 7;
0040   cfg.tMax      = 10 * dd4hep::ns;
0041   cfg.tdc_range = pow(2, cfg.tdc_bit);
0042   cfg.adc_range = pow(2, cfg.adc_bit);
0043   cfg.t_thres   = -cfg.adc_range * 0.1;
0044
0045   // check if max pulse height is linearly proportional to the initial Edep
0046   algo.applyConfig(cfg);
0047   algo.init();
0048
0049   SECTION("TDC vs analytic solution scan") {
0050     logger->info("Begin TDC vs analytic solution scan");
0051
0052     // NOLINTNEXTLINE(clang-analyzer-security.FloatLoopCounter)
0053     for (double time = -cfg.tMax; time <= cfg.tMax; time += cfg.tMax) {
0054       if (time < 0) {
0055         logger->info("Generation pulse at the negative time");
0056       } else if (time == 0) {
0057         logger->info("Generation pulse at the first EICROC cycle");
0058       } else {
0059         logger->info("Generation pulse at the second EICROC cycle");
0060       }
0061
0062       // test pulse with gaussian shape
0063       for (double tdc_frac = 0.4; tdc_frac < 1; tdc_frac += 0.1) {
0064         edm4hep::RawTimeSeriesCollection time_series_coll;
0065         auto rawhits_coll = std::make_unique<edm4eic::RawTrackerHitCollection>();
0066
0067         auto pulse = time_series_coll.create();
0068         auto cellID =
0069             id_desc.encode({{"system", 0}, {"module", 0}, {"sensor", 1}, {"x", 1}, {"y", 1}});
0070
0071         pulse.setCellID(cellID);
0072         pulse.setCharge(1.); // placeholder
0073         pulse.setTime(time); // placeholder
0074         pulse.setInterval(1);
0075
0076         int test_peak_TDC   = static_cast<int>(tdc_frac * cfg.tdc_range);
0077         int test_peak       = static_cast<int>(0.7 * cfg.adc_range);
0078         int test_peak_sigma = static_cast<int>(0.1 * cfg.tdc_range);
0079
0080         for (int i = 0; i < cfg.tdc_range; ++i) {
0081           int ADC =
0082               -test_peak *
0083               TMath::Exp(-0.5 * pow((i - test_peak_TDC) / static_cast<double>(test_peak_sigma), 2));
0084           pulse.addToAdcCounts(ADC);
0085         }
0086
0087         // calculate analytically when the pulse passes t_thres
0088         // ADC = amp*exp(-(TDC - mean)^2/(2sigma^2))
0089         // TDC = mean - (-2*sigma^2*ln(ADC/amp))^0.5
0090         int analytic_TDC = ceil(test_peak_TDC - sqrt(-2 * pow(test_peak_sigma, 2) *
0091                                                      TMath::Log(cfg.adc_range * 0.1 / test_peak)));
0092
0093         // Constructing input and output as per the algorithm's expected signature
0094         auto input  = std::make_tuple(&time_series_coll);
0095         auto output = std::make_tuple(rawhits_coll.get());
0096
0097         algo.process(input, output);
0098
0099         REQUIRE(rawhits_coll->size() == 1);
0100         auto hit = (*rawhits_coll)[0];
0101         REQUIRE(hit.getCellID() == cellID);
0102         REQUIRE(hit.getCharge() == test_peak);
0103         if (time < 0) {
0104           REQUIRE(hit.getTimeStamp() == analytic_TDC - cfg.tdc_range);
0105         } else if (time == 0) {
0106           REQUIRE(hit.getTimeStamp() == analytic_TDC);
0107         } else {
0108           REQUIRE(hit.getTimeStamp() == analytic_TDC + cfg.tdc_range);
0109         }
0110       }
0111     }
0112   }
0113
0114   SECTION("ADC scan") {
0115     logger->info("Begin ADC scan");
0116
0117     // test pulse with gaussian shape
0118     for (double adc_frac = 0.4; adc_frac < 1; adc_frac += 0.1) {
0119       edm4hep::RawTimeSeriesCollection time_series_coll;
0120       auto rawhits_coll = std::make_unique<edm4eic::RawTrackerHitCollection>();
0121
0122       auto pulse = time_series_coll.create();
0123       auto cellID =
0124           id_desc.encode({{"system", 0}, {"module", 0}, {"sensor", 1}, {"x", 1}, {"y", 1}});
0125
0126       pulse.setCellID(cellID);
0127       pulse.setCharge(1.); // placeholder
0128       pulse.setTime(0.);   // placeholder
0129       pulse.setInterval(1);
0130
0131       int test_peak_TDC   = static_cast<int>(0.5 * cfg.tdc_range);
0132       int test_peak       = static_cast<int>(adc_frac * cfg.adc_range);
0133       int test_peak_sigma = static_cast<int>(0.1 * cfg.tdc_range);
0134
0135       for (int i = 0; i < cfg.tdc_range; ++i) {
0136         int ADC =
0137             -test_peak *
0138             TMath::Exp(-0.5 * pow((i - test_peak_TDC) / static_cast<double>(test_peak_sigma), 2));
0139         pulse.addToAdcCounts(ADC);
0140       }
0141
0142       // Constructing input and output as per the algorithm's expected signature
0143       auto input  = std::make_tuple(&time_series_coll);
0144       auto output = std::make_tuple(rawhits_coll.get());
0145
0146       algo.process(input, output);
0147
0148       REQUIRE(rawhits_coll->size() == 1);
0149       auto hit = (*rawhits_coll)[0];
0150       REQUIRE(hit.getCellID() == cellID);
0151       REQUIRE(hit.getCharge() == test_peak);
0152     }
0153   }
0154 }