Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-05-12 08:02:40

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