tests/algorithms_test/particle_flow_CaloRemnantCombiner.cc

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2026, ePIC Collaboration
0003
0004 #include <algorithms/logger.h>
0005 #include <catch2/catch_test_macros.hpp>
0006 #include <edm4eic/ClusterCollection.h>
0007 #include <edm4eic/ReconstructedParticleCollection.h>
0008 #include <edm4hep/Vector3f.h>
0009 #include <cmath>
0010 #include <cstddef>
0011 #include <memory>
0012
0013 #include "algorithms/particle_flow/CaloRemnantCombiner.h"
0014 #include "algorithms/particle_flow/CaloRemnantCombinerConfig.h"
0015
0016 namespace {
0017
0018 //! Create a cluster at a given (eta, phi) with a given energy at radius r
0019 edm4eic::MutableCluster make_cluster(edm4eic::ClusterCollection& coll, float energy, float eta,
0020                                      float phi, float r = 2000.F) {
0021   auto cluster = coll.create();
0022   cluster.setEnergy(energy);
0023   cluster.setPosition(
0024       {r * std::cos(phi) / std::cosh(eta), r * std::sin(phi) / std::cosh(eta), r * std::tanh(eta)});
0025   return cluster;
0026 }
0027
0028 } // namespace
0029
0030 TEST_CASE("the CaloRemnantCombiner algorithm runs", "[CaloRemnantCombiner]") {
0031   eicrecon::CaloRemnantCombiner algo("test");
0032
0033   eicrecon::CaloRemnantCombinerConfig cfg;
0034   cfg.ecalDeltaR = 0.03;
0035   cfg.hcalDeltaR = 0.15;
0036
0037   algo.applyConfig(cfg);
0038   algo.level(algorithms::LogLevel::kDebug);
0039   algo.init();
0040
0041   auto ecal       = std::make_unique<edm4eic::ClusterCollection>();
0042   auto hcal       = std::make_unique<edm4eic::ClusterCollection>();
0043   auto candidates = std::make_unique<edm4eic::ReconstructedParticleCollection>();
0044
0045   SECTION("empty input produces empty output") {
0046     algo.process({ecal.get(), hcal.get()}, {candidates.get()});
0047
0048     REQUIRE(candidates->size() == 0);
0049   }
0050
0051   SECTION("single ecal cluster produces one candidate") {
0052     auto cluster = make_cluster(*ecal, 1.0F, 0.5F, 0.5F);
0053
0054     algo.process({ecal.get(), hcal.get()}, {candidates.get()});
0055
0056     REQUIRE(candidates->size() == 1);
0057     REQUIRE(candidates->at(0).clusters_size() == 1);
0058     REQUIRE(candidates->at(0).getClusters(0) == cluster);
0059   }
0060
0061   SECTION("single hcal cluster produces one candidate") {
0062     auto cluster = make_cluster(*hcal, 1.0F, 0.5F, 0.5F);
0063
0064     algo.process({ecal.get(), hcal.get()}, {candidates.get()});
0065
0066     REQUIRE(candidates->size() == 1);
0067     REQUIRE(candidates->at(0).clusters_size() == 1);
0068     REQUIRE(candidates->at(0).getClusters(0) == cluster);
0069   }
0070
0071   SECTION("hcal cluster within hcalDeltaR of ecal seed is merged") {
0072     auto ecal_cluster = make_cluster(*ecal, 1.0F, 0.5F, 0.5F);
0073     auto hcal_cluster = make_cluster(*hcal, 2.0F, 0.5F + 0.1F, 0.5F); // dR = 0.1 < 0.15
0074
0075     algo.process({ecal.get(), hcal.get()}, {candidates.get()});
0076
0077     REQUIRE(candidates->size() == 1);
0078     REQUIRE(candidates->at(0).clusters_size() == 2);
0079     REQUIRE(candidates->at(0).getClusters(0) == ecal_cluster);
0080     REQUIRE(candidates->at(0).getClusters(1) == hcal_cluster);
0081   }
0082
0083   SECTION("hcal cluster outside hcalDeltaR of ecal seed becomes its own candidate") {
0084     auto ecal_cluster = make_cluster(*ecal, 1.0F, 0.5F, 0.5F);
0085     auto hcal_cluster = make_cluster(*hcal, 2.0F, 0.5F + 0.5F, 0.5F); // dR = 0.5 > 0.15
0086
0087     algo.process({ecal.get(), hcal.get()}, {candidates.get()});
0088
0089     REQUIRE(candidates->size() == 2);
0090     // ecal-seeded candidate comes first
0091     REQUIRE(candidates->at(0).clusters_size() == 1);
0092     REQUIRE(candidates->at(0).getClusters(0) == ecal_cluster);
0093     // hcal-seeded candidate follows
0094     REQUIRE(candidates->at(1).clusters_size() == 1);
0095     REQUIRE(candidates->at(1).getClusters(0) == hcal_cluster);
0096   }
0097
0098   SECTION("nearby ecal clusters are merged into one candidate seeded by highest energy") {
0099     auto low  = make_cluster(*ecal, 1.0F, 0.5F + 0.02F, 0.5F); // dR = 0.02 < 0.03
0100     auto high = make_cluster(*ecal, 5.0F, 0.5F, 0.5F);
0101
0102     algo.process({ecal.get(), hcal.get()}, {candidates.get()});
0103
0104     REQUIRE(candidates->size() == 1);
0105     REQUIRE(candidates->at(0).clusters_size() == 2);
0106     // clusters are attached in descending energy order (seed first)
0107     REQUIRE(candidates->at(0).getClusters(0) == high);
0108     REQUIRE(candidates->at(0).getClusters(1) == low);
0109   }
0110
0111   SECTION("distant ecal clusters produce separate candidates") {
0112     auto low  = make_cluster(*ecal, 1.0F, -0.5F, 0.5F);
0113     auto high = make_cluster(*ecal, 5.0F, 0.5F, 0.5F);
0114
0115     algo.process({ecal.get(), hcal.get()}, {candidates.get()});
0116
0117     REQUIRE(candidates->size() == 2);
0118     // highest energy cluster is seeded first
0119     REQUIRE(candidates->at(0).clusters_size() == 1);
0120     REQUIRE(candidates->at(0).getClusters(0) == high);
0121     REQUIRE(candidates->at(1).clusters_size() == 1);
0122     REQUIRE(candidates->at(1).getClusters(0) == low);
0123   }
0124
0125   SECTION("nearby hcal clusters are merged into one candidate in phase 2") {
0126     auto high = make_cluster(*hcal, 5.0F, 0.5F, 0.5F);
0127     auto low  = make_cluster(*hcal, 1.0F, 0.5F, 0.5F + 0.1F); // dR = 0.1 < 0.15
0128
0129     algo.process({ecal.get(), hcal.get()}, {candidates.get()});
0130
0131     REQUIRE(candidates->size() == 1);
0132     REQUIRE(candidates->at(0).clusters_size() == 2);
0133     REQUIRE(candidates->at(0).getClusters(0) == high);
0134     REQUIRE(candidates->at(0).getClusters(1) == low);
0135   }
0136
0137   SECTION("clusters are merged across the phi = +/- pi boundary") {
0138     constexpr float pi = M_PI;
0139
0140     auto seed  = make_cluster(*ecal, 5.0F, 0.5F, pi - 0.01F);
0141     auto other = make_cluster(*ecal, 1.0F, 0.5F, -pi + 0.01F); // wrapped dphi = 0.02 < 0.03
0142
0143     algo.process({ecal.get(), hcal.get()}, {candidates.get()});
0144
0145     REQUIRE(candidates->size() == 1);
0146     REQUIRE(candidates->at(0).clusters_size() == 2);
0147     REQUIRE(candidates->at(0).getClusters(0) == seed);
0148     REQUIRE(candidates->at(0).getClusters(1) == other);
0149   }
0150
0151   SECTION("every input cluster ends up in exactly one candidate") {
0152     make_cluster(*ecal, 1.0F, -1.0F, 0.0F);
0153     make_cluster(*ecal, 2.0F, 0.0F, 1.0F);
0154     make_cluster(*ecal, 3.0F, 1.0F, 2.0F);
0155     make_cluster(*hcal, 1.5F, -1.0F, 0.0F);
0156     make_cluster(*hcal, 2.5F, 2.0F, -2.0F);
0157
0158     algo.process({ecal.get(), hcal.get()}, {candidates.get()});
0159
0160     std::size_t n_clusters = 0;
0161     for (const auto& candidate : *candidates) {
0162       n_clusters += candidate.clusters_size();
0163     }
0164     REQUIRE(n_clusters == ecal->size() + hcal->size());
0165   }
0166 }