tests/algorithms_test/digi_SiliconPulseGeneration.cc

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2025, Dmitry Kalinkin, Simon Gardner
0003
0004 #include <catch2/catch_test_macros.hpp>
0005 #include <catch2/generators/catch_generators.hpp>
0006 #include <edm4eic/EDM4eicVersion.h>
0007 #include <edm4eic/unit_system.h>
0008 #include <edm4hep/SimTrackerHitCollection.h>
0009 #if EDM4EIC_VERSION_MAJOR > 8 || (EDM4EIC_VERSION_MAJOR == 8 && EDM4EIC_VERSION_MINOR >= 1)
0010 #include <edm4eic/SimPulseCollection.h>
0011 #else
0012 #include <edm4hep/TimeSeriesCollection.h>
0013 #endif
0014 #include <podio/RelationRange.h>
0015 #include <cmath>
0016 #include <cstddef>
0017 #include <memory>
0018 #include <string>
0019 #include <tuple>
0020 #include <vector>
0021
0022 #include "algorithms/digi/SiliconPulseGeneration.h"
0023 #include "algorithms/digi/SiliconPulseGenerationConfig.h"
0024
0025 #if EDM4EIC_VERSION_MAJOR > 8 || (EDM4EIC_VERSION_MAJOR == 8 && EDM4EIC_VERSION_MINOR >= 1)
0026 using PulseType = edm4eic::SimPulse;
0027 #else
0028 using PulseType = edm4hep::TimeSeries;
0029 #endif
0030
0031 TEST_CASE("SiliconPulseGeneration generates correct number of pulses", "[SiliconPulseGeneration]") {
0032
0033   eicrecon::SiliconPulseGeneration algo("SiliconPulseGeneration");
0034   eicrecon::SiliconPulseGenerationConfig cfg;
0035   cfg.pulse_shape_function = "LandauPulse"; // Example pulse shape
0036   cfg.pulse_shape_params   = {1.0, 1.0};    // Example parameters for the pulse shape
0037   cfg.ignore_thres         = 1;
0038
0039   algo.applyConfig(cfg);
0040   algo.init();
0041
0042   SECTION("from different hits") {
0043     std::size_t nHits = GENERATE(1, 2, 3);
0044
0045     edm4hep::SimTrackerHitCollection hits_coll;
0046
0047     for (std::size_t i = 0; i < nHits; i++) {
0048       hits_coll.create(12345 + i, 10.0, 5.0); // cellID, charge, time
0049     }
0050
0051     auto pulses = std::make_unique<PulseType::collection_type>();
0052
0053     auto input  = std::make_tuple(&hits_coll);
0054     auto output = std::make_tuple(pulses.get());
0055
0056     algo.process(input, output);
0057
0058     REQUIRE(pulses->size() == nHits);
0059     REQUIRE((*pulses)[0].getCellID() == 12345);
0060     if (nHits > 1) {
0061       REQUIRE((*pulses)[1].getCellID() == 12346);
0062     }
0063     if (nHits > 2) {
0064       REQUIRE((*pulses)[2].getCellID() == 12347);
0065     }
0066   }
0067 }
0068
0069 TEST_CASE("Test the EvaluatorSvc pulse generation with a square pulse",
0070           "[SiliconPulseGeneration]") {
0071
0072   eicrecon::SiliconPulseGeneration algo("SiliconPulseGeneration");
0073   eicrecon::SiliconPulseGenerationConfig cfg;
0074
0075   // Square wave expression
0076   std::string expression = "(time >= param0 && time < param1) ? charge : 0";
0077
0078   double startTime      = 0.0 * edm4eic::unit::ns;
0079   double endTime        = 1.0 * edm4eic::unit::ns;
0080   std::size_t nTimeBins = 10;
0081   double timeStep       = (endTime - startTime) / nTimeBins;
0082
0083   cfg.pulse_shape_function = expression;
0084   cfg.pulse_shape_params   = {startTime, endTime}; // Example parameters for the square pulse
0085   cfg.ignore_thres         = 1;
0086   cfg.timestep             = timeStep;
0087   cfg.min_sampling_time    = startTime + timeStep;
0088
0089   algo.applyConfig(cfg);
0090   algo.init();
0091
0092   double charge      = 10.0 * cfg.ignore_thres;
0093   double time        = GENERATE_COPY(0.0, 0.5 * timeStep, timeStep);
0094   float rounded_time = std::floor(time / timeStep) * timeStep;
0095
0096   edm4hep::SimTrackerHitCollection hits_coll;
0097   hits_coll.create(12345, charge, time); // cellID, charge, time
0098
0099   auto pulses = std::make_unique<PulseType::collection_type>();
0100
0101   auto input  = std::make_tuple(&hits_coll);
0102   auto output = std::make_tuple(pulses.get());
0103
0104   algo.process(input, output);
0105
0106   REQUIRE(pulses->size() == 1);
0107   REQUIRE((*pulses)[0].getCellID() == 12345);
0108   REQUIRE((*pulses)[0].getTime() == rounded_time);
0109   auto amplitudes = (*pulses)[0].getAmplitude();
0110   REQUIRE(amplitudes.size() == nTimeBins); // Two time bins for the square pulse
0111   for (auto amplitude : amplitudes) {
0112     REQUIRE(amplitude == charge); // All time bins should be zero
0113   }
0114 }