Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2026-10-24 08:25:21 

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2026 Chun Yuen Tsang, Minho Kim
0003 
0004 #include <DD4hep/Detector.h>
0005 #include <DD4hep/IDDescriptor.h>
0006 #include <DD4hep/Readout.h>
0007 #include <algorithms/geo.h>
0008 #include <catch2/catch_test_macros.hpp>
0009 #include <catch2/matchers/catch_matchers.hpp>
0010 #include <catch2/matchers/catch_matchers_floating_point.hpp>
0011 #include <edm4eic/CALOROC1BSample.h>
0012 #include <edm4eic/CalorimeterHitCollection.h>
0013 #include <edm4hep/MCParticleCollection.h>
0014 #include <edm4eic/MCRecoCalorimeterHitAssociationCollection.h>
0015 #include <edm4eic/MCRecoCalorimeterHitLinkCollection.h>
0016 #include <edm4eic/RawCALOROCHitCollection.h>
0017 #include <edm4hep/CaloHitContributionCollection.h>
0018 #include <edm4hep/RawCalorimeterHitCollection.h>
0019 #include <edm4hep/SimCalorimeterHitCollection.h>
0020 #include <edm4hep/Vector3f.h>
0021 #include <podio/detail/LinkCollectionIterator.h>
0022 #include <stdint.h>
0023 #include <algorithm>
0024 #include <deque>
0025 #include <filesystem>
0026 #include <fstream>
0027 #include <gsl/pointers>
0028 #include <initializer_list>
0029 #include <iostream>
0030 #include <memory>
0031 #include <string>
0032 #include <utility>
0033 #include <vector>
0034 
0035 #include "algorithms/calorimetry/CalorimeterCALOROCReco.h"
0036 #include "algorithms/calorimetry/CalorimeterCALOROCRecoConfig.h"
0037 
0038 // Helper: build a RawCALOROCHit with one B-sample
0039 // timeStamp and samplePhase control the TOA formula; highGainADC/lowGainADC control amplitude
0040 static edm4eic::MutableRawCALOROCHit make_raw_hit(edm4eic::RawCALOROCHitCollection& coll,
0041                                                   uint64_t cellID, int32_t timeStamp,
0042                                                   int32_t samplePhase, uint16_t highGainADC,
0043                                                   uint16_t lowGainADC, uint16_t timeOfArrival) {
0044   auto hit = coll.create();
0045   hit.setCellID(cellID);
0046   hit.setTimeStamp(timeStamp);
0047   hit.setSamplePhase(samplePhase);
0048   edm4eic::CALOROC1BSample sample;
0049   sample.highGainADC   = highGainADC;
0050   sample.lowGainADC    = lowGainADC;
0051   sample.timeOfArrival = timeOfArrival;
0052   hit.addToBSamples(sample);
0053   return hit;
0054 }
0055 
0056 // Helper: build an NpeHit (SimCalorimeterHit) with given contributions
0057 static edm4hep::MutableSimCalorimeterHit
0058 make_npe_hit(edm4hep::SimCalorimeterHitCollection& hitColl,
0059              edm4hep::CaloHitContributionCollection& contribColl, uint64_t cellID,
0060              edm4hep::Vector3f position, std::initializer_list<float> contrib_energies,
0061              std::initializer_list<edm4hep::MCParticle> particles = {}) {
0062   REQUIRE((particles.size() == 0 || particles.size() == contrib_energies.size()));
0063 
0064   float total = 0.f;
0065   for (auto e : contrib_energies)
0066     total += e;
0067   auto hit = hitColl.create();
0068   hit.setCellID(cellID);
0069   hit.setEnergy(total);
0070   hit.setPosition(position);
0071   auto particle = particles.begin();
0072   for (auto e : contrib_energies) {
0073     auto c = contribColl.create();
0074     c.setEnergy(e);
0075     c.setTime(0.f);
0076     c.setStepPosition({0.f, 0.f, 0.f});
0077     if (particle != particles.end())
0078       c.setParticle(*particle++);
0079     hit.addToContributions(c);
0080   }
0081   return hit;
0082 }
0083 
0084 // Write a minimal LUT file with one entry: layer=0 sector=0 -> factor=1.0
0085 // Returns the file path as a string
0086 static std::string write_lut_file(const std::string& path) {
0087   std::ofstream f(path);
0088   f << "0 0 1.0\n";
0089   return path;
0090 }
0091 
0092 // Convenience: build the standard config used in most tests
0093 static eicrecon::CalorimeterCALOROCRecoConfig make_cfg(const std::string& lut_path) {
0094   eicrecon::CalorimeterCALOROCRecoConfig cfg;
0095   cfg.readout              = "MockCALOROCHits";
0096   cfg.layerField           = "layer";
0097   cfg.sectorField          = "sector";
0098   cfg.localDetFields       = {"system"};
0099   cfg.edep_to_npe_filename = lut_path;
0100   cfg.edep_to_npe_fields   = {"layer", "sector"};
0101 
0102   // Identity-like attenuation: factor ≈ 1 for all reasonable z
0103   cfg.attenuationParameters = {1.0, 1.0e9, 1.0e9};
0104 
0105   cfg.attenuationReferencePositionNamePos = "MockCALOROC_RefPosP";
0106   cfg.attenuationReferencePositionNameNeg = "MockCALOROC_RefPosN";
0107 
0108   cfg.slope     = 1.0;
0109   cfg.intercept = 0.0;
0110 
0111   cfg.highGainDR = 1000; // saturation threshold
0112   cfg.gainRatio  = 1.0;  // high/low gain ratio = 1 for simplicity
0113 
0114   // Light speed parameters: zpos = dt * 1.0 + 0.0
0115   cfg.lightSpeedParameters = {1.0, 0.0};
0116 
0117   cfg.timeWalkCor                  = false; // disabled by default; enabled per test-case
0118   cfg.timeWalkCorrectionParameters = {0.0, 1.0, 0.0, 1.0};
0119   // correction = toa - (1.0 * (ADCsum - 0)^1 + 0) = toa - ADCsum
0120 
0121   cfg.useNpeHitPos = false; // use timing-based z by default
0122   cfg.proxy_type   = eicrecon::CalorimeterCALOROCRecoConfig::ProxyType::sum;
0123 
0124   return cfg;
0125 }
0126 
0127 // ── TOA formula (for expected-value computation in tests) ──────────────────────
0128 // toa = samplePhase*(25/1042) + (timeStamp + idx_toa)*25 - timeOfArrival*(25/1024)
0129 static double expected_toa(int samplePhase, int timeStamp, int idx_toa, uint16_t toa_val) {
0130   return samplePhase * (25.0 / 1042.0) + (timeStamp + idx_toa) * 25.0 - toa_val * (25.0 / 1024.0);
0131 }
0132 
0133 // ────────────────────────────────────────────────────────────────────────────────
0134 // Test 1: Amplitude from ADC sum (no saturation)
0135 // ────────────────────────────────────────────────────────────────────────────────
0136 TEST_CASE("CalorimeterCALOROCReco: ADC amplitude in low-signal regime",
0137           "[CalorimeterCALOROCReco][Amplitude]") {
0138 
0139   // Write LUT
0140   const std::string lut_path =
0141       std::filesystem::temp_directory_path() / "caloroc_test_lut_amplitude.txt";
0142   write_lut_file(lut_path);
0143 
0144   auto cfg = make_cfg(lut_path);
0145 
0146   eicrecon::CalorimeterCALOROCReco algo("CALOROCRecoAmplitude");
0147   algo.applyConfig(cfg);
0148   algo.init();
0149 
0150   auto detector = algorithms::GeoSvc::instance().detector();
0151   auto id_desc  = detector->readout("MockCALOROCHits").idSpec();
0152 
0153   // cellID with layer=0, sector=0 (matches LUT entry)
0154   uint64_t cellID =
0155       id_desc.encode({{"system", 4}, {"layer", 0}, {"sector", 0}, {"x", 0}, {"y", 0}});
0156 
0157   // P-side: highGainADC=100, below highGainDR=1000 -> npeP = 100/1.0 = 100
0158   // N-side: highGainADC=200 -> npeN = 200
0159   edm4eic::RawCALOROCHitCollection adcP_coll, adcN_coll;
0160   // Symmetric timing on both sides: tP = tN = -12.5, dt = 0, zpos = 0.
0161   // Timing is irrelevant for this test; both sides use identical timeStamp=0, toa=512.
0162   make_raw_hit(adcP_coll, cellID, /*timeStamp=*/0, /*samplePhase=*/0,
0163                /*highGainADC=*/100, /*lowGainADC=*/50, /*timeOfArrival=*/512);
0164   make_raw_hit(adcN_coll, cellID, /*timeStamp=*/0, /*samplePhase=*/0,
0165                /*highGainADC=*/200, /*lowGainADC=*/100, /*timeOfArrival=*/512);
0166 
0167   edm4hep::SimCalorimeterHitCollection npeP_coll, npeN_coll;
0168   edm4hep::CaloHitContributionCollection contribs_coll;
0169   make_npe_hit(npeP_coll, contribs_coll, cellID, {0.f, 0.f, 0.f}, {60.f});
0170   make_npe_hit(npeN_coll, contribs_coll, cellID, {0.f, 0.f, 0.f}, {40.f});
0171 
0172   auto recohits = std::make_unique<edm4eic::CalorimeterHitCollection>();
0173   auto rawhits  = std::make_unique<edm4hep::RawCalorimeterHitCollection>();
0174   edm4eic::MCRecoCalorimeterHitLinkCollection links;
0175   edm4eic::MCRecoCalorimeterHitAssociationCollection assocs;
0176 
0177   algo.process({&npeP_coll, &adcP_coll, &npeN_coll, &adcN_coll},
0178                {recohits.get(), rawhits.get(), &links, &assocs});
0179 
0180   const int64_t actual_amplitude = (*rawhits)[0].getAmplitude();
0181   std::cout << "--- ADC amplitude (no saturation) ---" << "\n";
0182   std::cout << "P-side highGainADC = 100, gainRatio = 1.0" << "\n";
0183   std::cout << "Expected amplitude (npeP = highGainADC / gainRatio) = 100" << "\n";
0184   std::cout << "Number of rawhits produced = " << rawhits->size() << "\n";
0185   std::cout << "Actual rawhit amplitude    = " << actual_amplitude << "\n";
0186   REQUIRE(rawhits->size() == 1);
0187   REQUIRE(actual_amplitude == 100);
0188 }
0189 
0190 // ────────────────────────────────────────────────────────────────────────────────
0191 // Test 2: Amplitude uses lowGainADC when highGain saturates
0192 // ────────────────────────────────────────────────────────────────────────────────
0193 TEST_CASE("CalorimeterCALOROCReco: ADC amplitude switches to low-gain on saturation",
0194           "[CalorimeterCALOROCReco][Amplitude]") {
0195 
0196   const std::string lut_path =
0197       std::filesystem::temp_directory_path() / "caloroc_test_lut_saturation.txt";
0198   write_lut_file(lut_path);
0199 
0200   auto cfg = make_cfg(lut_path);
0201 
0202   eicrecon::CalorimeterCALOROCReco algo("CALOROCRecoSaturation");
0203   algo.applyConfig(cfg);
0204   algo.init();
0205 
0206   auto detector = algorithms::GeoSvc::instance().detector();
0207   auto id_desc  = detector->readout("MockCALOROCHits").idSpec();
0208   uint64_t cellID =
0209       id_desc.encode({{"system", 4}, {"layer", 0}, {"sector", 0}, {"x", 0}, {"y", 0}});
0210 
0211   edm4eic::RawCALOROCHitCollection adcP_coll, adcN_coll;
0212   // Symmetric timing on both sides: timing is not under test here.
0213   // P-side: highGainADC=1000 >= highGainDR=1000 -> use lowGainADC=75 directly
0214   make_raw_hit(adcP_coll, cellID, 0, 0, /*highGainADC=*/1000, /*lowGainADC=*/75, /*toa=*/512);
0215   // N-side: also saturated -> lowGainADC=150
0216   make_raw_hit(adcN_coll, cellID, 0, 0, /*highGainADC=*/1000, /*lowGainADC=*/150, /*toa=*/512);
0217 
0218   edm4hep::SimCalorimeterHitCollection npeP_coll, npeN_coll;
0219   edm4hep::CaloHitContributionCollection contribs_coll;
0220   make_npe_hit(npeP_coll, contribs_coll, cellID, {0.f, 0.f, 0.f}, {60.f});
0221   make_npe_hit(npeN_coll, contribs_coll, cellID, {0.f, 0.f, 0.f}, {40.f});
0222 
0223   auto recohits = std::make_unique<edm4eic::CalorimeterHitCollection>();
0224   auto rawhits  = std::make_unique<edm4hep::RawCalorimeterHitCollection>();
0225   edm4eic::MCRecoCalorimeterHitLinkCollection links;
0226   edm4eic::MCRecoCalorimeterHitAssociationCollection assocs;
0227 
0228   algo.process({&npeP_coll, &adcP_coll, &npeN_coll, &adcN_coll},
0229                {recohits.get(), rawhits.get(), &links, &assocs});
0230 
0231   const int64_t actual_amplitude = (*rawhits)[0].getAmplitude();
0232   std::cout << "--- ADC amplitude (high-gain saturation) ---" << "\n";
0233   std::cout << "P-side highGainADC = 1000 >= highGainDR = 1000 -> switch to lowGainADC = 75"
0234             << "\n";
0235   std::cout << "Expected amplitude = 75" << "\n";
0236   std::cout << "Number of rawhits produced = " << rawhits->size() << "\n";
0237   std::cout << "Actual rawhit amplitude    = " << actual_amplitude << "\n";
0238   REQUIRE(rawhits->size() == 1);
0239   REQUIRE(actual_amplitude == 75);
0240 }
0241 
0242 // ────────────────────────────────────────────────────────────────────────────────
0243 // Test 3: z-position from delta-t
0244 // ────────────────────────────────────────────────────────────────────────────────
0245 TEST_CASE("CalorimeterCALOROCReco: z-position from timing difference",
0246           "[CalorimeterCALOROCReco][ZPosition]") {
0247 
0248   const std::string lut_path = std::filesystem::temp_directory_path() / "caloroc_test_lut_zpos.txt";
0249   write_lut_file(lut_path);
0250 
0251   auto cfg         = make_cfg(lut_path);
0252   cfg.timeWalkCor  = false;
0253   cfg.useNpeHitPos = false;
0254 
0255   eicrecon::CalorimeterCALOROCReco algo("CALOROCRecoZPos");
0256   algo.applyConfig(cfg);
0257   algo.init();
0258 
0259   auto detector = algorithms::GeoSvc::instance().detector();
0260   auto id_desc  = detector->readout("MockCALOROCHits").idSpec();
0261   uint64_t cellID =
0262       id_desc.encode({{"system", 4}, {"layer", 0}, {"sector", 0}, {"x", 0}, {"y", 0}});
0263 
0264   // P-side: timeStamp=0, samplePhase=0, toa=512 -> tP = 0 + 0*25 - 512*(25/1024) = -12.5
0265   // N-side: timeStamp=1, samplePhase=0, toa=512 -> tN = 0 + 1*25 - 12.5 = 12.5
0266   // dt = tN - tP = 12.5 - (-12.5) = 25.0
0267   // zpos = 25.0 * 1.0 + 0.0 = 25.0
0268   edm4eic::RawCALOROCHitCollection adcP_coll, adcN_coll;
0269   make_raw_hit(adcP_coll, cellID, 0, 0, 100, 50, 512);
0270   make_raw_hit(adcN_coll, cellID, 1, 0, 200, 100, 512);
0271 
0272   edm4hep::SimCalorimeterHitCollection npeP_coll, npeN_coll;
0273   edm4hep::CaloHitContributionCollection contribs_coll;
0274   make_npe_hit(npeP_coll, contribs_coll, cellID, {0.f, 0.f, 0.f}, {60.f});
0275   make_npe_hit(npeN_coll, contribs_coll, cellID, {0.f, 0.f, 0.f}, {40.f});
0276 
0277   auto recohits = std::make_unique<edm4eic::CalorimeterHitCollection>();
0278   auto rawhits  = std::make_unique<edm4hep::RawCalorimeterHitCollection>();
0279   edm4eic::MCRecoCalorimeterHitLinkCollection links;
0280   edm4eic::MCRecoCalorimeterHitAssociationCollection assocs;
0281 
0282   algo.process({&npeP_coll, &adcP_coll, &npeN_coll, &adcN_coll},
0283                {recohits.get(), rawhits.get(), &links, &assocs});
0284 
0285   const double tP_expected   = expected_toa(0, 0, 0, 512); // -12.5
0286   const double tN_expected   = expected_toa(0, 1, 0, 512); //  12.5
0287   const double dt_expected   = tN_expected - tP_expected;  //  25.0
0288   const double zpos_expected = dt_expected * 1.0 + 0.0;    //  25.0
0289   const double actual_z      = (*recohits)[0].getPosition().z;
0290   std::cout << "--- z-position from timing difference ---" << "\n";
0291   std::cout << "P-side: timeStamp=0, samplePhase=0, timeOfArrival=512" << "\n";
0292   std::cout << "N-side: timeStamp=1, samplePhase=0, timeOfArrival=512" << "\n";
0293   std::cout << "tP (expected)   = " << tP_expected << " ns" << "\n";
0294   std::cout << "tN (expected)   = " << tN_expected << " ns" << "\n";
0295   std::cout << "dt = tN - tP    = " << dt_expected << " ns" << "\n";
0296   std::cout << "zpos = dt * lightSpeed[0] + lightSpeed[1] = " << zpos_expected << " mm" << "\n";
0297   std::cout << "Number of recohits = " << recohits->size() << "\n";
0298   std::cout << "Actual recohit z   = " << actual_z << " mm" << "\n";
0299   REQUIRE(recohits->size() == 1);
0300   REQUIRE_THAT(actual_z, Catch::Matchers::WithinAbs(zpos_expected, 1e-4));
0301 }
0302 
0303 // ────────────────────────────────────────────────────────────────────────────────
0304 // Test 4: z-position from NpeHit (useNpeHitPos=true)
0305 // ────────────────────────────────────────────────────────────────────────────────
0306 TEST_CASE("CalorimeterCALOROCReco: z-position read from NpeHit position",
0307           "[CalorimeterCALOROCReco][ZPosition]") {
0308 
0309   const std::string lut_path =
0310       std::filesystem::temp_directory_path() / "caloroc_test_lut_npepos.txt";
0311   write_lut_file(lut_path);
0312 
0313   auto cfg         = make_cfg(lut_path);
0314   cfg.useNpeHitPos = true;
0315 
0316   eicrecon::CalorimeterCALOROCReco algo("CALOROCRecoNpePos");
0317   algo.applyConfig(cfg);
0318   algo.init();
0319 
0320   auto detector = algorithms::GeoSvc::instance().detector();
0321   auto id_desc  = detector->readout("MockCALOROCHits").idSpec();
0322   uint64_t cellID =
0323       id_desc.encode({{"system", 4}, {"layer", 0}, {"sector", 0}, {"x", 0}, {"y", 0}});
0324 
0325   edm4eic::RawCALOROCHitCollection adcP_coll, adcN_coll;
0326   make_raw_hit(adcP_coll, cellID, 0, 0, 100, 50, 512);
0327   make_raw_hit(adcN_coll, cellID, 1, 0, 200, 100, 512);
0328 
0329   edm4hep::SimCalorimeterHitCollection npeP_coll, npeN_coll;
0330   edm4hep::CaloHitContributionCollection contribs_coll;
0331   // NpeHitP has z = 42.0 mm — this should be used as zpos
0332   make_npe_hit(npeP_coll, contribs_coll, cellID, {0.f, 0.f, 42.f}, {60.f});
0333   make_npe_hit(npeN_coll, contribs_coll, cellID, {0.f, 0.f, 0.f}, {40.f});
0334 
0335   auto recohits = std::make_unique<edm4eic::CalorimeterHitCollection>();
0336   auto rawhits  = std::make_unique<edm4hep::RawCalorimeterHitCollection>();
0337   edm4eic::MCRecoCalorimeterHitLinkCollection links;
0338   edm4eic::MCRecoCalorimeterHitAssociationCollection assocs;
0339 
0340   algo.process({&npeP_coll, &adcP_coll, &npeN_coll, &adcN_coll},
0341                {recohits.get(), rawhits.get(), &links, &assocs});
0342 
0343   const double npe_z    = 42.0;
0344   const double actual_z = (*recohits)[0].getPosition().z;
0345   std::cout << "--- z-position from NpeHit position (useNpeHitPos=true) ---" << "\n";
0346   std::cout << "NpeHitP position.z = " << npe_z << " mm  (set explicitly)" << "\n";
0347   std::cout << "Expected recohit z = " << npe_z << " mm" << "\n";
0348   std::cout << "Number of recohits = " << recohits->size() << "\n";
0349   std::cout << "Actual recohit z   = " << actual_z << " mm" << "\n";
0350   REQUIRE(recohits->size() == 1);
0351   REQUIRE_THAT(actual_z, Catch::Matchers::WithinAbs(npe_z, 1e-4));
0352 }
0353 
0354 // ────────────────────────────────────────────────────────────────────────────────
0355 // Test 5: Time walk correction shifts the reconstructed z
0356 // ────────────────────────────────────────────────────────────────────────────────
0357 TEST_CASE("CalorimeterCALOROCReco: time walk correction changes z from delta-t",
0358           "[CalorimeterCALOROCReco][TimeWalk]") {
0359 
0360   const std::string lut_path =
0361       std::filesystem::temp_directory_path() / "caloroc_test_lut_timewalk.txt";
0362   write_lut_file(lut_path);
0363 
0364   auto detector = algorithms::GeoSvc::instance().detector();
0365   auto id_desc  = detector->readout("MockCALOROCHits").idSpec();
0366   uint64_t cellID =
0367       id_desc.encode({{"system", 4}, {"layer", 0}, {"sector", 0}, {"x", 0}, {"y", 0}});
0368 
0369   edm4eic::RawCALOROCHitCollection adcP_coll, adcN_coll;
0370   // P-side ADC sum = 100 (highGainADC=100 < 1000, gainRatio=1)
0371   // N-side ADC sum = 200
0372   make_raw_hit(adcP_coll, cellID, 0, 0, 100, 50, 512);
0373   make_raw_hit(adcN_coll, cellID, 1, 0, 200, 100, 512);
0374 
0375   edm4hep::SimCalorimeterHitCollection npeP_coll, npeN_coll;
0376   edm4hep::CaloHitContributionCollection contribs_coll;
0377   make_npe_hit(npeP_coll, contribs_coll, cellID, {0.f, 0.f, 0.f}, {60.f});
0378   make_npe_hit(npeN_coll, contribs_coll, cellID, {0.f, 0.f, 0.f}, {40.f});
0379 
0380   // Without time walk correction:
0381   // tP = -12.5, tN = 12.5, dt = 25.0, zpos = 25.0
0382   auto cfg_nocor        = make_cfg(lut_path);
0383   cfg_nocor.timeWalkCor = false;
0384   {
0385     eicrecon::CalorimeterCALOROCReco algo("CALOROCRecoNoTWC");
0386     algo.applyConfig(cfg_nocor);
0387     algo.init();
0388 
0389     auto recohits = std::make_unique<edm4eic::CalorimeterHitCollection>();
0390     auto rawhits  = std::make_unique<edm4hep::RawCalorimeterHitCollection>();
0391     edm4eic::MCRecoCalorimeterHitLinkCollection links;
0392     edm4eic::MCRecoCalorimeterHitAssociationCollection assocs;
0393     algo.process({&npeP_coll, &adcP_coll, &npeN_coll, &adcN_coll},
0394                  {recohits.get(), rawhits.get(), &links, &assocs});
0395 
0396     const double tP_raw   = expected_toa(0, 0, 0, 512); // -12.5
0397     const double tN_raw   = expected_toa(0, 1, 0, 512); //  12.5
0398     const double dt_raw   = tN_raw - tP_raw;
0399     const double zpos_raw = dt_raw * 1.0 + 0.0;
0400     const double actual_z = (*recohits)[0].getPosition().z;
0401     std::cout << "--- Time walk correction OFF ---" << "\n";
0402     std::cout << "tP (raw)         = " << tP_raw << " ns" << "\n";
0403     std::cout << "tN (raw)         = " << tN_raw << " ns" << "\n";
0404     std::cout << "dt (raw)         = " << dt_raw << " ns" << "\n";
0405     std::cout << "Expected zpos    = " << zpos_raw << " mm" << "\n";
0406     std::cout << "Actual recohit z = " << actual_z << " mm" << "\n";
0407     REQUIRE(recohits->size() == 1);
0408     REQUIRE_THAT(actual_z, Catch::Matchers::WithinAbs(25.0, 1e-4));
0409   }
0410 
0411   // With time walk correction (params = {0, 1, 0, 1}):
0412   // corrected tP = -12.5 - (1.0*(100-0)^1 + 0) = -112.5
0413   // corrected tN = 12.5  - (1.0*(200-0)^1 + 0) = -187.5
0414   // dt_corrected = -187.5 - (-112.5) = -75.0
0415   // zpos_corrected = -75.0 (within [-1000, 1000], no clamping)
0416   auto cfg_cor                         = make_cfg(lut_path);
0417   cfg_cor.timeWalkCor                  = true;
0418   cfg_cor.timeWalkCorrectionParameters = {0.0, 1.0, 0.0, 1.0};
0419   {
0420     eicrecon::CalorimeterCALOROCReco algo("CALOROCRecoWithTWC");
0421     algo.applyConfig(cfg_cor);
0422     algo.init();
0423 
0424     auto recohits = std::make_unique<edm4eic::CalorimeterHitCollection>();
0425     auto rawhits  = std::make_unique<edm4hep::RawCalorimeterHitCollection>();
0426     edm4eic::MCRecoCalorimeterHitLinkCollection links;
0427     edm4eic::MCRecoCalorimeterHitAssociationCollection assocs;
0428     algo.process({&npeP_coll, &adcP_coll, &npeN_coll, &adcN_coll},
0429                  {recohits.get(), rawhits.get(), &links, &assocs});
0430 
0431     const double tP_raw  = expected_toa(0, 0, 0, 512); // -12.5
0432     const double tN_raw  = expected_toa(0, 1, 0, 512); //  12.5
0433     const double adcSumP = 100.0;                      // highGainADC=100, gainRatio=1
0434     const double adcSumN = 200.0;                      // highGainADC=200, gainRatio=1
0435     // correction = toa - (p[1] * (ADCsum - p[2])^p[3] + p[0])
0436     //            = toa - (1.0  * (ADCsum -    0)^1.0  + 0.0)
0437     //            = toa - ADCsum
0438     const double tP_corrected   = tP_raw - adcSumP;            // -12.5 - 100 = -112.5
0439     const double tN_corrected   = tN_raw - adcSumN;            //  12.5 - 200 = -187.5
0440     const double dt_corrected   = tN_corrected - tP_corrected; // -75.0
0441     const double zpos_corrected = dt_corrected * 1.0 + 0.0;    // -75.0
0442     const double actual_z       = (*recohits)[0].getPosition().z;
0443     std::cout << "--- Time walk correction ON ---" << "\n";
0444     std::cout << "timeWalkCorrectionParameters = {0, 1, 0, 1}" << "\n";
0445     std::cout << "tP (raw)         = " << tP_raw << " ns" << "\n";
0446     std::cout << "tN (raw)         = " << tN_raw << " ns" << "\n";
0447     std::cout << "ADC sum P-side   = " << adcSumP << "\n";
0448     std::cout << "ADC sum N-side   = " << adcSumN << "\n";
0449     std::cout << "tP (corrected)   = " << tP_corrected << " ns" << "\n";
0450     std::cout << "tN (corrected)   = " << tN_corrected << " ns" << "\n";
0451     std::cout << "dt (corrected)   = " << dt_corrected << " ns" << "\n";
0452     std::cout << "Expected zpos    = " << zpos_corrected << " mm" << "\n";
0453     std::cout << "Actual recohit z = " << actual_z << " mm" << "\n";
0454     REQUIRE(recohits->size() == 1);
0455     REQUIRE_THAT(actual_z, Catch::Matchers::WithinAbs(-75.0, 1e-4));
0456   }
0457 }
0458 
0459 // ────────────────────────────────────────────────────────────────────────────────
0460 // Test 6: MC truth links – correct counts and energy-normalized weights
0461 // ────────────────────────────────────────────────────────────────────────────────
0462 TEST_CASE("CalorimeterCALOROCReco: MC truth link weights are energy-normalized",
0463           "[CalorimeterCALOROCReco][MCLinks]") {
0464 
0465   const std::string lut_path =
0466       std::filesystem::temp_directory_path() / "caloroc_test_lut_links.txt";
0467   write_lut_file(lut_path);
0468 
0469   auto cfg = make_cfg(lut_path);
0470 
0471   eicrecon::CalorimeterCALOROCReco algo("CALOROCRecoLinks");
0472   algo.applyConfig(cfg);
0473   algo.init();
0474 
0475   auto detector = algorithms::GeoSvc::instance().detector();
0476   auto id_desc  = detector->readout("MockCALOROCHits").idSpec();
0477   uint64_t cellID =
0478       id_desc.encode({{"system", 4}, {"layer", 0}, {"sector", 0}, {"x", 0}, {"y", 0}});
0479 
0480   edm4eic::RawCALOROCHitCollection adcP_coll, adcN_coll;
0481   // Symmetric timing: timing irrelevant here, both sides identical to keep zpos=0.
0482   make_raw_hit(adcP_coll, cellID, 0, 0, 100, 50, 512);
0483   make_raw_hit(adcN_coll, cellID, 0, 0, 200, 100, 512);
0484 
0485   // NpeHitP has 2 contributions: 30 and 70 -> total 100
0486   // NpeHitN has 1 contribution: 50
0487   // Total edep = 30 + 70 + 50 = 150
0488   edm4hep::SimCalorimeterHitCollection npeP_coll, npeN_coll;
0489   edm4hep::CaloHitContributionCollection contribs_coll;
0490   edm4hep::MCParticleCollection particles;
0491   auto particleP1 = particles.create();
0492   auto particleP2 = particles.create();
0493   auto particleN1 = particles.create();
0494   make_npe_hit(npeP_coll, contribs_coll, cellID, {0.f, 0.f, 0.f}, {30.f, 70.f},
0495                {particleP1, particleP2});
0496   make_npe_hit(npeN_coll, contribs_coll, cellID, {0.f, 0.f, 0.f}, {50.f}, {particleN1});
0497 
0498   auto recohits = std::make_unique<edm4eic::CalorimeterHitCollection>();
0499   auto rawhits  = std::make_unique<edm4hep::RawCalorimeterHitCollection>();
0500   edm4eic::MCRecoCalorimeterHitLinkCollection links;
0501   edm4eic::MCRecoCalorimeterHitAssociationCollection assocs;
0502 
0503   algo.process({&npeP_coll, &adcP_coll, &npeN_coll, &adcN_coll},
0504                {recohits.get(), rawhits.get(), &links, &assocs});
0505 
0506   const double edepP1 = 30.0, edepP2 = 70.0, edepN1 = 50.0;
0507   const double total_edep = edepP1 + edepP2 + edepN1;
0508   std::vector<double> weights;
0509   for (const auto& lnk : links)
0510     weights.push_back(lnk.getWeight());
0511   std::sort(weights.begin(), weights.end());
0512   double total_weight = 0.;
0513   for (auto w : weights)
0514     total_weight += w;
0515 
0516   std::cout << "--- MC truth link weights ---" << "\n";
0517   std::cout << "NpeHitP contributions: " << edepP1 << ", " << edepP2 << "\n";
0518   std::cout << "NpeHitN contributions: " << edepN1 << "\n";
0519   std::cout << "Total edep = " << total_edep << "\n";
0520   std::cout << "Expected weights (sorted): " << edepP1 / total_edep << ", " << edepN1 / total_edep
0521             << ", " << edepP2 / total_edep << "\n";
0522   std::cout << "Number of links  = " << links.size() << "\n";
0523   std::cout << "Number of assocs = " << assocs.size() << "\n";
0524   std::cout << "Sum of weights   = " << total_weight << "\n";
0525   if (weights.size() == 3) {
0526     std::cout << "Actual weights (sorted): " << weights[0] << ", " << weights[1] << ", "
0527               << weights[2] << "\n";
0528   }
0529 
0530   REQUIRE(links.size() == 3);
0531   REQUIRE(assocs.size() == 3);
0532   REQUIRE_THAT(total_weight, Catch::Matchers::WithinAbs(1.0, 1e-5));
0533   REQUIRE_THAT(weights[0], Catch::Matchers::WithinAbs(30.0 / 150.0, 1e-5));
0534   REQUIRE_THAT(weights[1], Catch::Matchers::WithinAbs(50.0 / 150.0, 1e-5));
0535   REQUIRE_THAT(weights[2], Catch::Matchers::WithinAbs(70.0 / 150.0, 1e-5));
0536 }
0537 
0538 // ────────────────────────────────────────────────────────────────────────────────
0539 // Test 7: Average time stored in rawhit timestamp
0540 // ────────────────────────────────────────────────────────────────────────────────
0541 TEST_CASE("CalorimeterCALOROCReco: rawhit timestamp is average of P and N TOA",
0542           "[CalorimeterCALOROCReco][Timing]") {
0543 
0544   const std::string lut_path = std::filesystem::temp_directory_path() / "caloroc_test_lut_time.txt";
0545   write_lut_file(lut_path);
0546 
0547   auto cfg        = make_cfg(lut_path);
0548   cfg.timeWalkCor = false;
0549 
0550   eicrecon::CalorimeterCALOROCReco algo("CALOROCRecoTime");
0551   algo.applyConfig(cfg);
0552   algo.init();
0553 
0554   auto detector = algorithms::GeoSvc::instance().detector();
0555   auto id_desc  = detector->readout("MockCALOROCHits").idSpec();
0556   uint64_t cellID =
0557       id_desc.encode({{"system", 4}, {"layer", 0}, {"sector", 0}, {"x", 0}, {"y", 0}});
0558 
0559   // tP = expected_toa(0, 0, 0, 512) = -12.5
0560   // tN = expected_toa(0, 1, 0, 512) =  12.5
0561   // time = 0.5 * (-12.5 + 12.5) = 0.0
0562   edm4eic::RawCALOROCHitCollection adcP_coll, adcN_coll;
0563   make_raw_hit(adcP_coll, cellID, 0, 0, 100, 50, 512);
0564   make_raw_hit(adcN_coll, cellID, 1, 0, 200, 100, 512);
0565 
0566   edm4hep::SimCalorimeterHitCollection npeP_coll, npeN_coll;
0567   edm4hep::CaloHitContributionCollection contribs_coll;
0568   make_npe_hit(npeP_coll, contribs_coll, cellID, {0.f, 0.f, 0.f}, {60.f});
0569   make_npe_hit(npeN_coll, contribs_coll, cellID, {0.f, 0.f, 0.f}, {40.f});
0570 
0571   auto recohits = std::make_unique<edm4eic::CalorimeterHitCollection>();
0572   auto rawhits  = std::make_unique<edm4hep::RawCalorimeterHitCollection>();
0573   edm4eic::MCRecoCalorimeterHitLinkCollection links;
0574   edm4eic::MCRecoCalorimeterHitAssociationCollection assocs;
0575 
0576   algo.process({&npeP_coll, &adcP_coll, &npeN_coll, &adcN_coll},
0577                {recohits.get(), rawhits.get(), &links, &assocs});
0578 
0579   const double tP          = expected_toa(0, 0, 0, 512);
0580   const double tN          = expected_toa(0, 1, 0, 512);
0581   const double expected    = 0.5 * (tP + tN);
0582   const double actual_time = static_cast<double>((*rawhits)[0].getTimeStamp());
0583   std::cout << "--- Rawhit timestamp = average of P and N TOA ---" << "\n";
0584   std::cout << "P-side: timeStamp=0, samplePhase=0, timeOfArrival=512" << "\n";
0585   std::cout << "N-side: timeStamp=1, samplePhase=0, timeOfArrival=512" << "\n";
0586   std::cout << "tP (formula)       = " << tP << " ns" << "\n";
0587   std::cout << "tN (formula)       = " << tN << " ns" << "\n";
0588   std::cout << "Expected timestamp = 0.5*(tP+tN) = " << expected << " ns" << "\n";
0589   std::cout << "Number of rawhits  = " << rawhits->size() << "\n";
0590   std::cout << "Actual timestamp   = " << actual_time << " ns" << "\n";
0591   REQUIRE(rawhits->size() == 1);
0592   REQUIRE_THAT(actual_time, Catch::Matchers::WithinAbs(expected, 1e-4));
0593 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!