tests/algorithms_test/tracking_MPGDHitReconstruction.cc

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2026 ePIC Collaboration
0003
0004 #include <DD4hep/Detector.h>
0005 #include <DD4hep/IDDescriptor.h>
0006 #include <DD4hep/Readout.h>
0007 #include <DDSegmentation/BitFieldCoder.h>
0008 #include <algorithms/geo.h>
0009 #include <catch2/catch_approx.hpp>
0010 #include <catch2/catch_test_macros.hpp>
0011 #include <edm4eic/RawTrackerHitCollection.h>
0012 #include <edm4eic/TrackerHitCollection.h>
0013 #include <gsl/pointers>
0014 #include <string>
0015
0016 #include "algorithms/tracking/MPGDHitReconstruction.h"
0017 #include "algorithms/tracking/MPGDHitReconstructionConfig.h"
0018
0019 using eicrecon::MPGDHitReconstruction;
0020 using eicrecon::MPGDHitReconstructionConfig;
0021
0022 // Helper: build a CellID from field values using the IDDescriptor encoder.
0023 // system=3 matches the mock geometry's addPhysVolID("system", 3).
0024 // strip=1 → p-strip (coordinate in x field at offset 32)
0025 // strip=2 → n-strip (coordinate in y field at offset 48)
0026 static dd4hep::DDSegmentation::CellID makeCellID(const dd4hep::IDDescriptor& desc, int system,
0027                                                  int layer, int module, int strip, int x, int y) {
0028   auto encoder                       = desc.decoder();
0029   dd4hep::DDSegmentation::CellID cid = 0;
0030   encoder->set(cid, "system", system);
0031   encoder->set(cid, "layer", layer);
0032   encoder->set(cid, "module", module);
0033   encoder->set(cid, "strip", strip);
0034   encoder->set(cid, "x", x);
0035   encoder->set(cid, "y", y);
0036   return cid;
0037 }
0038
0039 static dd4hep::IDDescriptor getMPGDIdDesc() {
0040   return algorithms::GeoSvc::instance().detector()->readout("MockMPGDHits").idSpec();
0041 }
0042
0043 TEST_CASE("MPGDHitReconstruction: empty input produces empty output", "[MPGDHitReconstruction]") {
0044   MPGDHitReconstruction algo("test_empty");
0045   MPGDHitReconstructionConfig cfg;
0046   cfg.readout        = "MockMPGDHits";
0047   cfg.timeResolution = 10.0;
0048   algo.applyConfig(cfg);
0049   algo.init();
0050
0051   edm4eic::RawTrackerHitCollection raw_hits;
0052   edm4eic::TrackerHitCollection rec_hits;
0053
0054   algo.process({&raw_hits}, {&rec_hits});
0055   REQUIRE(rec_hits.size() == 0);
0056 }
0057
0058 TEST_CASE("MPGDHitReconstruction: single p-strip hit produces one cluster",
0059           "[MPGDHitReconstruction]") {
0060   MPGDHitReconstruction algo("test_single_p");
0061   MPGDHitReconstructionConfig cfg;
0062   cfg.readout        = "MockMPGDHits";
0063   cfg.timeResolution = 10.0;
0064   algo.applyConfig(cfg);
0065   algo.init();
0066
0067   auto id_desc     = getMPGDIdDesc();
0068   const int pStrip = 1;
0069
0070   edm4eic::RawTrackerHitCollection raw_hits;
0071   edm4eic::TrackerHitCollection rec_hits;
0072
0073   auto cellID = makeCellID(id_desc, 3, 0, 0, pStrip, 10, 0);
0074   raw_hits.create(cellID, /*charge=*/500, /*timeStamp=*/100000);
0075
0076   algo.process({&raw_hits}, {&rec_hits});
0077
0078   REQUIRE(rec_hits.size() == 1);
0079   // Single-hit cluster: charge is the input charge (converted to GeV)
0080   CHECK(rec_hits[0].getEdep() > 0);
0081   // CellID of the reconstructed hit should match the max-charge hit
0082   CHECK(rec_hits[0].getCellID() == cellID);
0083   // Raw hit association
0084   CHECK(rec_hits[0].getRawHit().getCellID() == cellID);
0085 }
0086
0087 TEST_CASE("MPGDHitReconstruction: single n-strip hit produces one cluster",
0088           "[MPGDHitReconstruction]") {
0089   MPGDHitReconstruction algo("test_single_n");
0090   MPGDHitReconstructionConfig cfg;
0091   cfg.readout        = "MockMPGDHits";
0092   cfg.timeResolution = 10.0;
0093   algo.applyConfig(cfg);
0094   algo.init();
0095
0096   auto id_desc     = getMPGDIdDesc();
0097   const int nStrip = 2;
0098
0099   edm4eic::RawTrackerHitCollection raw_hits;
0100   edm4eic::TrackerHitCollection rec_hits;
0101
0102   auto cellID = makeCellID(id_desc, 3, 0, 0, nStrip, 0, 15);
0103   raw_hits.create(cellID, /*charge=*/700, /*timeStamp=*/200000);
0104
0105   algo.process({&raw_hits}, {&rec_hits});
0106
0107   REQUIRE(rec_hits.size() == 1);
0108   CHECK(rec_hits[0].getCellID() == cellID);
0109   CHECK(rec_hits[0].getRawHit().getCellID() == cellID);
0110 }
0111
0112 TEST_CASE("MPGDHitReconstruction: two adjacent p-strip hits cluster together",
0113           "[MPGDHitReconstruction]") {
0114   MPGDHitReconstruction algo("test_adj_p");
0115   MPGDHitReconstructionConfig cfg;
0116   cfg.readout        = "MockMPGDHits";
0117   cfg.timeResolution = 10.0;
0118   algo.applyConfig(cfg);
0119   algo.init();
0120
0121   auto id_desc     = getMPGDIdDesc();
0122   const int pStrip = 1;
0123
0124   edm4eic::RawTrackerHitCollection raw_hits;
0125   edm4eic::TrackerHitCollection rec_hits;
0126
0127   // Two adjacent p-strip hits at x=10 and x=11 (same system/layer/module/strip)
0128   auto cid0 = makeCellID(id_desc, 3, 0, 0, pStrip, 10, 0);
0129   auto cid1 = makeCellID(id_desc, 3, 0, 0, pStrip, 11, 0);
0130   raw_hits.create(cid0, /*charge=*/600, /*timeStamp=*/100000);
0131   raw_hits.create(cid1, /*charge=*/400, /*timeStamp=*/110000);
0132
0133   algo.process({&raw_hits}, {&rec_hits});
0134
0135   // Adjacent hits in the same subvolume and strip type → 1 cluster
0136   REQUIRE(rec_hits.size() == 1);
0137   // Total charge = 600 + 400 = 1000 (converted to GeV: 1000 / 1e6)
0138   CHECK(rec_hits[0].getEdep() > 0);
0139   // CellID from the max-charge hit
0140   CHECK(rec_hits[0].getCellID() == cid0);
0141   // Timing from the max-charge hit (600 > 400)
0142   CHECK(rec_hits[0].getTime() > 0);
0143 }
0144
0145 TEST_CASE("MPGDHitReconstruction: two adjacent n-strip hits cluster together",
0146           "[MPGDHitReconstruction]") {
0147   MPGDHitReconstruction algo("test_adj_n");
0148   MPGDHitReconstructionConfig cfg;
0149   cfg.readout        = "MockMPGDHits";
0150   cfg.timeResolution = 10.0;
0151   algo.applyConfig(cfg);
0152   algo.init();
0153
0154   auto id_desc     = getMPGDIdDesc();
0155   const int nStrip = 2;
0156
0157   edm4eic::RawTrackerHitCollection raw_hits;
0158   edm4eic::TrackerHitCollection rec_hits;
0159
0160   auto cid0 = makeCellID(id_desc, 3, 0, 0, nStrip, 0, 20);
0161   auto cid1 = makeCellID(id_desc, 3, 0, 0, nStrip, 0, 21);
0162   raw_hits.create(cid0, /*charge=*/300, /*timeStamp=*/100000);
0163   raw_hits.create(cid1, /*charge=*/500, /*timeStamp=*/150000);
0164
0165   algo.process({&raw_hits}, {&rec_hits});
0166
0167   REQUIRE(rec_hits.size() == 1);
0168   // CellID from the max-charge hit (500 > 300)
0169   CHECK(rec_hits[0].getCellID() == cid1);
0170 }
0171
0172 TEST_CASE("MPGDHitReconstruction: separated hits produce separate clusters",
0173           "[MPGDHitReconstruction]") {
0174   MPGDHitReconstruction algo("test_separated");
0175   MPGDHitReconstructionConfig cfg;
0176   cfg.readout        = "MockMPGDHits";
0177   cfg.timeResolution = 10.0;
0178   algo.applyConfig(cfg);
0179   algo.init();
0180
0181   auto id_desc     = getMPGDIdDesc();
0182   const int pStrip = 1;
0183
0184   edm4eic::RawTrackerHitCollection raw_hits;
0185   edm4eic::TrackerHitCollection rec_hits;
0186
0187   // Two p-strip hits at x=5 and x=50 — far apart, not adjacent
0188   auto cid0 = makeCellID(id_desc, 3, 0, 0, pStrip, 5, 0);
0189   auto cid1 = makeCellID(id_desc, 3, 0, 0, pStrip, 50, 0);
0190   raw_hits.create(cid0, /*charge=*/500, /*timeStamp=*/100000);
0191   raw_hits.create(cid1, /*charge=*/500, /*timeStamp=*/100000);
0192
0193   algo.process({&raw_hits}, {&rec_hits});
0194
0195   // Non-adjacent → 2 separate clusters
0196   REQUIRE(rec_hits.size() == 2);
0197 }
0198
0199 TEST_CASE("MPGDHitReconstruction: p-strip and n-strip hits produce separate clusters",
0200           "[MPGDHitReconstruction]") {
0201   MPGDHitReconstruction algo("test_pn_sep");
0202   MPGDHitReconstructionConfig cfg;
0203   cfg.readout        = "MockMPGDHits";
0204   cfg.timeResolution = 10.0;
0205   algo.applyConfig(cfg);
0206   algo.init();
0207
0208   auto id_desc     = getMPGDIdDesc();
0209   const int pStrip = 1;
0210   const int nStrip = 2;
0211
0212   edm4eic::RawTrackerHitCollection raw_hits;
0213   edm4eic::TrackerHitCollection rec_hits;
0214
0215   // One p-strip and one n-strip hit — different strip types → separate clusters
0216   auto cidP = makeCellID(id_desc, 3, 0, 0, pStrip, 10, 0);
0217   auto cidN = makeCellID(id_desc, 3, 0, 0, nStrip, 0, 10);
0218   raw_hits.create(cidP, /*charge=*/500, /*timeStamp=*/100000);
0219   raw_hits.create(cidN, /*charge=*/500, /*timeStamp=*/100000);
0220
0221   algo.process({&raw_hits}, {&rec_hits});
0222
0223   REQUIRE(rec_hits.size() == 2);
0224 }
0225
0226 TEST_CASE("MPGDHitReconstruction: three-hit cluster sums charge correctly",
0227           "[MPGDHitReconstruction]") {
0228   MPGDHitReconstruction algo("test_3hit");
0229   MPGDHitReconstructionConfig cfg;
0230   cfg.readout        = "MockMPGDHits";
0231   cfg.timeResolution = 10.0;
0232   algo.applyConfig(cfg);
0233   algo.init();
0234
0235   auto id_desc     = getMPGDIdDesc();
0236   const int pStrip = 1;
0237
0238   edm4eic::RawTrackerHitCollection raw_hits;
0239   edm4eic::TrackerHitCollection rec_hits;
0240
0241   // Three adjacent p-strip hits at x=10,11,12
0242   auto cid0 = makeCellID(id_desc, 3, 0, 0, pStrip, 10, 0);
0243   auto cid1 = makeCellID(id_desc, 3, 0, 0, pStrip, 11, 0);
0244   auto cid2 = makeCellID(id_desc, 3, 0, 0, pStrip, 12, 0);
0245   raw_hits.create(cid0, /*charge=*/200, /*timeStamp=*/100000);
0246   raw_hits.create(cid1, /*charge=*/800, /*timeStamp=*/110000);
0247   raw_hits.create(cid2, /*charge=*/300, /*timeStamp=*/120000);
0248
0249   algo.process({&raw_hits}, {&rec_hits});
0250
0251   REQUIRE(rec_hits.size() == 1);
0252   // Charge is (200+800+300)/1e6 GeV = 0.0013
0253   float expected_edep = (200.0f + 800.0f + 300.0f) / 1.0e6f;
0254   CHECK(rec_hits[0].getEdep() == Catch::Approx(expected_edep).epsilon(0.01));
0255   // CellID from max-charge hit (800 at cid1)
0256   CHECK(rec_hits[0].getCellID() == cid1);
0257 }
0258
0259 TEST_CASE("MPGDHitReconstruction: cluster timing from max-charge hit", "[MPGDHitReconstruction]") {
0260   MPGDHitReconstruction algo("test_timing");
0261   MPGDHitReconstructionConfig cfg;
0262   cfg.readout        = "MockMPGDHits";
0263   cfg.timeResolution = 10.0;
0264   algo.applyConfig(cfg);
0265   algo.init();
0266
0267   auto id_desc     = getMPGDIdDesc();
0268   const int pStrip = 1;
0269
0270   edm4eic::RawTrackerHitCollection raw_hits;
0271   edm4eic::TrackerHitCollection rec_hits;
0272
0273   auto cid0 = makeCellID(id_desc, 3, 0, 0, pStrip, 10, 0);
0274   auto cid1 = makeCellID(id_desc, 3, 0, 0, pStrip, 11, 0);
0275   // Hit 1 has higher charge → its timestamp determines cluster time
0276   raw_hits.create(cid0, /*charge=*/300, /*timeStamp=*/100000);
0277   raw_hits.create(cid1, /*charge=*/700, /*timeStamp=*/200000);
0278
0279   algo.process({&raw_hits}, {&rec_hits});
0280
0281   REQUIRE(rec_hits.size() == 1);
0282   // Time = max-charge timestamp / 1000.0 (conversion to ns)
0283   float expected_time = 200000.0f / 1000.0f;
0284   CHECK(rec_hits[0].getTime() == Catch::Approx(expected_time).epsilon(0.01));
0285 }
0286
0287 TEST_CASE("MPGDHitReconstruction: out-of-order input is sorted and clustered",
0288           "[MPGDHitReconstruction]") {
0289   MPGDHitReconstruction algo("test_sort");
0290   MPGDHitReconstructionConfig cfg;
0291   cfg.readout        = "MockMPGDHits";
0292   cfg.timeResolution = 10.0;
0293   algo.applyConfig(cfg);
0294   algo.init();
0295
0296   auto id_desc     = getMPGDIdDesc();
0297   const int pStrip = 1;
0298
0299   edm4eic::RawTrackerHitCollection raw_hits;
0300   edm4eic::TrackerHitCollection rec_hits;
0301
0302   // Insert in reverse coordinate order: x=11 first, x=10 second
0303   auto cid1 = makeCellID(id_desc, 3, 0, 0, pStrip, 11, 0);
0304   auto cid0 = makeCellID(id_desc, 3, 0, 0, pStrip, 10, 0);
0305   raw_hits.create(cid1, /*charge=*/400, /*timeStamp=*/110000);
0306   raw_hits.create(cid0, /*charge=*/600, /*timeStamp=*/100000);
0307
0308   algo.process({&raw_hits}, {&rec_hits});
0309
0310   // Still clusters together (algorithm sorts internally)
0311   REQUIRE(rec_hits.size() == 1);
0312   // CellID from max-charge hit (600 at cid0)
0313   CHECK(rec_hits[0].getCellID() == cid0);
0314 }
0315
0316 TEST_CASE("MPGDHitReconstruction: mixed p/n strips with adjacency", "[MPGDHitReconstruction]") {
0317   MPGDHitReconstruction algo("test_mixed");
0318   MPGDHitReconstructionConfig cfg;
0319   cfg.readout        = "MockMPGDHits";
0320   cfg.timeResolution = 10.0;
0321   algo.applyConfig(cfg);
0322   algo.init();
0323
0324   auto id_desc     = getMPGDIdDesc();
0325   const int pStrip = 1;
0326   const int nStrip = 2;
0327
0328   edm4eic::RawTrackerHitCollection raw_hits;
0329   edm4eic::TrackerHitCollection rec_hits;
0330
0331   // Two adjacent p-strips + one n-strip
0332   auto cidP0 = makeCellID(id_desc, 3, 0, 0, pStrip, 10, 0);
0333   auto cidP1 = makeCellID(id_desc, 3, 0, 0, pStrip, 11, 0);
0334   auto cidN0 = makeCellID(id_desc, 3, 0, 0, nStrip, 0, 10);
0335   raw_hits.create(cidP0, /*charge=*/500, /*timeStamp=*/100000);
0336   raw_hits.create(cidP1, /*charge=*/500, /*timeStamp=*/100000);
0337   raw_hits.create(cidN0, /*charge=*/500, /*timeStamp=*/100000);
0338
0339   algo.process({&raw_hits}, {&rec_hits});
0340
0341   // p-strip cluster (2 hits) + n-strip cluster (1 hit) = 2 clusters
0342   REQUIRE(rec_hits.size() == 2);
0343 }
0344
0345 TEST_CASE("MPGDHitReconstruction: different modules produce separate clusters",
0346           "[MPGDHitReconstruction]") {
0347   MPGDHitReconstruction algo("test_modules");
0348   MPGDHitReconstructionConfig cfg;
0349   cfg.readout        = "MockMPGDHits";
0350   cfg.timeResolution = 10.0;
0351   algo.applyConfig(cfg);
0352   algo.init();
0353
0354   auto id_desc     = getMPGDIdDesc();
0355   const int pStrip = 1;
0356
0357   edm4eic::RawTrackerHitCollection raw_hits;
0358   edm4eic::TrackerHitCollection rec_hits;
0359
0360   // Same strip type, same coordinates, but different modules
0361   auto cid0 = makeCellID(id_desc, 3, 0, /*module=*/0, pStrip, 10, 0);
0362   auto cid1 = makeCellID(id_desc, 3, 0, /*module=*/1, pStrip, 10, 0);
0363   raw_hits.create(cid0, /*charge=*/500, /*timeStamp=*/100000);
0364   raw_hits.create(cid1, /*charge=*/500, /*timeStamp=*/100000);
0365
0366   algo.process({&raw_hits}, {&rec_hits});
0367
0368   // Different modules → different subvolumes → 2 separate clusters
0369   REQUIRE(rec_hits.size() == 2);
0370 }