tests/algorithms_test/digi_LGADPulseGeneration.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 <algorithms/geo.h>
0009 #include <catch2/catch_test_macros.hpp> // for AssertionHandler, operator""_catch_sr, StringRef, REQUIRE, operator<, operator==, operator>, TEST_CASE
0010 #include <edm4hep/RawTimeSeriesCollection.h>
0011 #include <edm4hep/SimTrackerHitCollection.h>
0012 #include <fmt/core.h>
0013 #include <podio/RelationRange.h>
0014 #include <spdlog/common.h> // for level_enum
0015 #include <spdlog/logger.h> // for logger
0016 #include <spdlog/spdlog.h> // for default_logger
0017 #include <algorithm>
0018 #include <cmath>
0019 #include <gsl/pointers>
0020 #include <limits>
0021 #include <memory> // for allocator, unique_ptr, make_unique, shared_ptr, __shared_ptr_access
0022 #include <tuple>
0023 #include <utility>
0024 #include <vector>
0025
0026 #include "TF1.h"
0027 #include "TGraphErrors.h"
0028 #include "algorithms/digi/LGADPulseGeneration.h"
0029 #include "algorithms/digi/LGADPulseGenerationConfig.h"
0030
0031 TEST_CASE("the LGAD charge sharing algorithm runs", "[LGADPulseGeneration]") {
0032   const float EPSILON = 1e-5;
0033   eicrecon::LGADPulseGenerationConfig cfg;
0034   cfg.gain         = 113;
0035   cfg.Vm           = 1e-4 * dd4hep::GeV;
0036   cfg.ignore_thres = 1e-4 / 5;
0037   cfg.sigma_analog = 0.293951 * dd4hep::ns;
0038   cfg.adc_bit      = 8;
0039   cfg.adc_range    = pow(2, cfg.adc_bit);
0040
0041   eicrecon::LGADPulseGeneration algo("LGADPulseGeneration");
0042
0043   std::shared_ptr<spdlog::logger> logger = spdlog::default_logger()->clone("LGADPulseGeneration");
0044   logger->set_level(spdlog::level::trace);
0045
0046
0047   auto detector = algorithms::GeoSvc::instance().detector();
0048   auto id_desc  = detector->readout("MockLGADHits").idSpec();
0049
0050   SECTION("Pulse height linearlity test") {
0051     // check if max pulse height is linearly proportional to the initial Edep
0052     algo.applyConfig(cfg);
0053     algo.init();
0054
0055     TGraphErrors graph;
0056     for (double edep = 0; edep <= 1e-4; edep += 1e-4 / 9) {
0057       edm4hep::SimTrackerHitCollection rawhits_coll;
0058       auto time_series_coll = std::make_unique<edm4hep::RawTimeSeriesCollection>();
0059
0060       auto hit = rawhits_coll.create();
0061       auto cellID =
0062           id_desc.encode({{"system", 0}, {"module", 0}, {"sensor", 1}, {"x", 1}, {"y", 1}});
0063
0064       hit.setCellID(cellID);
0065       hit.setEDep(
0066           edep); // in GeV. Since Vm = 1e-4*gain, EDep = 1e-4 GeV corresponds to ADC = max_adc
0067       hit.setTime(1.5 * dd4hep::ns); // in ns
0068
0069       // Constructing input and output as per the algorithm's expected signature
0070       auto input  = std::make_tuple(&rawhits_coll);
0071       auto output = std::make_tuple(time_series_coll.get());
0072
0073       algo.process(input, output);
0074
0075       if (edep < 1e-4 / 5)
0076         REQUIRE(time_series_coll->size() == 0);
0077       else {
0078         REQUIRE(time_series_coll->size() == 1);
0079         auto min_adc = std::numeric_limits<int>::max();
0080         for(const auto& pulse : (*time_series_coll)) {
0081           REQUIRE(pulse.getCellID() == cellID);
0082           auto adcs    = pulse.getAdcCounts();
0083           for (const auto adc : adcs)
0084             min_adc = std::min(min_adc, adc);
0085         }
0086         int npt = graph.GetN();
0087         graph.SetPoint(npt, edep, min_adc);
0088         graph.SetPointError(npt, 0, 0.5);
0089         // make sure when energy deposition = Vm, ADC reaches max value
0090         if (edep == 1e-4)
0091           REQUIRE(min_adc == -cfg.adc_range + 1);
0092       }
0093     }
0094
0095     // test linearlity
0096     TF1 tf1("tf1", "pol1", 0, 1e-4);
0097     graph.Fit(&tf1, "R0");
0098     // slope can't be consistent with zero
0099     REQUIRE(!(tf1.GetParameter(1) - tf1.GetParError(1) < 0 && 0 < tf1.GetParameter(1) + tf1.GetParError(1) ));
0100     double chi2_dof = tf1.GetChisquare() / tf1.GetNDF();
0101     logger->info("Chi-square/dof value for Edep vs min-adc = {}", chi2_dof);
0102     REQUIRE(chi2_dof < 2);
0103   }
0104
0105   SECTION("Pulse timing linearlity test") {
0106     // check if max pulse height is linearly proportional to the initial Edep
0107     algo.applyConfig(cfg);
0108     algo.init();
0109
0110     TGraphErrors graph;
0111     std::vector<double> times;
0112
0113     // test within the same EICROC cycle
0114     for (double time = 0; time < 12; time += 1.) times.push_back(time);
0115     // test multiple EICROC cycle
0116     for (double time = 10; time < 101; time += 25.) times.push_back(time);
0117     // test negative time
0118     times.push_back(-10);
0119
0120     for (double time : times) {
0121       edm4hep::SimTrackerHitCollection rawhits_coll;
0122       auto time_series_coll = std::make_unique<edm4hep::RawTimeSeriesCollection>();
0123
0124       auto hit = rawhits_coll.create();
0125       auto cellID =
0126           id_desc.encode({{"system", 0}, {"module", 0}, {"sensor", 1}, {"x", 1}, {"y", 1}});
0127
0128       hit.setCellID(cellID);
0129       hit.setEDep(
0130           0.5e-4 * dd4hep::GeV); // in GeV. Since Vm = 1e-4*gain, EDep = 1e-4 GeV corresponds to ADC = max_adc
0131       hit.setTime(time); // in ns
0132
0133       // Constructing input and output as per the algorithm's expected signature
0134       auto input  = std::make_tuple(&rawhits_coll);
0135       auto output = std::make_tuple(time_series_coll.get());
0136
0137       algo.process(input, output);
0138
0139       REQUIRE(time_series_coll->size() == 1);
0140       auto min_adc          = std::numeric_limits<int>::max();
0141       int time_bin = 0;
0142       for(const auto& pulse: *time_series_coll) {
0143         //auto pulse = (*time_series_coll)[0];
0144         REQUIRE(pulse.getCellID() == cellID);
0145
0146         auto adcs             = pulse.getAdcCounts();
0147         for (unsigned int i = 0; i < adcs.size(); ++i) {
0148           auto adc = adcs[i];
0149           if (adc < min_adc)
0150             time_bin = i + pulse.getTime()/cfg.tMax*cfg.tdc_range;
0151           min_adc = std::min(min_adc, adc);
0152         }
0153       }
0154       int npt = graph.GetN();
0155       graph.SetPoint(npt, time, time_bin);
0156       graph.SetPointError(npt, 0, 0.5);
0157     }
0158
0159     // test linearlity
0160     TF1 tf1("tf1", "pol1", *std::min_element(times.begin(), times.end()), *std::max_element(times.begin(), times.end()));
0161     graph.Fit(&tf1, "R0");
0162     // slope can't be consistent with zero
0163     REQUIRE(!(tf1.GetParameter(1) - tf1.GetParError(1) < 0 && 0 < tf1.GetParameter(1) + tf1.GetParError(1) ));
0164     double chi2_dof = tf1.GetChisquare() / tf1.GetNDF();
0165     logger->info("Chi-square/dof value for time vs TDC-bin = {}", chi2_dof);
0166     REQUIRE(chi2_dof < 2);
0167   }
0168 }