tests/algorithms_test/reco_ClustersToParticles.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 <catch2/matchers/catch_matchers.hpp>
0007 #include <catch2/matchers/catch_matchers_floating_point.hpp>
0008 #include <edm4eic/ClusterCollection.h>
0009 #include <edm4eic/EDM4eicVersion.h>
0010 #include <edm4eic/MCRecoClusterParticleAssociationCollection.h>
0011 #include <edm4eic/MCRecoParticleAssociationCollection.h>
0012 #include <edm4eic/ReconstructedParticleCollection.h>
0013 #include <edm4hep/EDM4hepVersion.h>
0014 #include <edm4hep/MCParticleCollection.h>
0015 #include <edm4hep/Vector3f.h>
0016 #include <podio/detail/Link.h>
0017 #include <podio/detail/LinkCollectionImpl.h>
0018 #if EDM4EIC_BUILD_VERSION >= EDM4EIC_VERSION(8, 7, 0)
0019 #include <edm4eic/MCRecoParticleLinkCollection.h>
0020 #endif
0021 #if EDM4HEP_BUILD_VERSION < EDM4HEP_VERSION(0, 99, 2)
0022 #include <edm4hep/Vector2i.h>
0023 #endif
0024 #include <edm4hep/Vector3d.h>
0025 #include <cmath>
0026 #include <deque>
0027 #include <memory>
0028
0029 #include "algorithms/reco/ClustersToParticles.h"
0030 #include "algorithms/reco/ClustersToParticlesConfig.h"
0031
0032 using eicrecon::ClustersToParticles;
0033 using eicrecon::ClustersToParticlesConfig;
0034
0035 TEST_CASE("the ClustersToParticles algorithm runs", "[ClustersToParticles]") {
0036   const float EPSILON = 1e-5;
0037
0038   ClustersToParticles algo("ClustersToParticles");
0039   algo.level(algorithms::LogLevel::kDebug);
0040
0041   ClustersToParticlesConfig cfg;
0042   algo.applyConfig(cfg);
0043   algo.init();
0044
0045   // Create input: two clusters at different positions with different energies,
0046   // plus one MC association on the first cluster
0047   edm4eic::ClusterCollection clusters;
0048
0049   auto cluster1 = clusters.create();
0050   cluster1.setEnergy(5.0);
0051   cluster1.setPosition({0.0f, 0.0f, 100.0f});
0052
0053   auto cluster2 = clusters.create();
0054   cluster2.setEnergy(10.0);
0055   // Cluster at 45 degrees in x-z plane
0056   cluster2.setPosition({1.0f, 0.0f, 1.0f});
0057
0058   // Create an MC particle and associate it with cluster1
0059   edm4hep::MCParticleCollection mcparts;
0060   auto mcpart = mcparts.create(22,                  // std::int32_t PDG
0061                                1,                   // std::int32_t generatorStatus
0062                                0,                   // std::int32_t simulatorStatus
0063                                0.,                  // float charge
0064                                0.,                  // float time
0065                                0.,                  // double mass
0066                                edm4hep::Vector3d(), // edm4hep::Vector3d vertex
0067                                edm4hep::Vector3d(), // edm4hep::Vector3d endpoint
0068 #if EDM4HEP_BUILD_VERSION < EDM4HEP_VERSION(0, 99, 1)
0069                                edm4hep::Vector3f(), // edm4hep::Vector3f momentum
0070                                edm4hep::Vector3f(), // edm4hep::Vector3f momentumAtEndpoint
0071 #else
0072                                edm4hep::Vector3d(), // edm4hep::Vector3d momentum
0073                                edm4hep::Vector3d(), // edm4hep::Vector3d momentumAtEndpoint
0074 #endif
0075 #if EDM4HEP_BUILD_VERSION < EDM4HEP_VERSION(0, 99, 3)
0076                                edm4hep::Vector3f() // edm4hep::Vector3f spin
0077 #else
0078                                9 // int32_t helicity (9 if unset)
0079 #endif
0080 #if EDM4HEP_BUILD_VERSION < EDM4HEP_VERSION(0, 99, 2)
0081                                ,
0082                                edm4hep::Vector2i() // edm4hep::Vector2i colorFlow
0083 #endif
0084   );
0085
0086   edm4eic::MCRecoClusterParticleAssociationCollection cluster_assocs;
0087   auto cluster_assoc = cluster_assocs.create();
0088   cluster_assoc.setRec(cluster1);
0089   cluster_assoc.setSim(mcpart);
0090   cluster_assoc.setWeight(0.9f);
0091
0092   // Run algorithm
0093   auto parts       = std::make_unique<edm4eic::ReconstructedParticleCollection>();
0094   auto part_assocs = std::make_unique<edm4eic::MCRecoParticleAssociationCollection>();
0095 #if EDM4EIC_BUILD_VERSION >= EDM4EIC_VERSION(8, 7, 0)
0096   auto part_links = std::make_unique<edm4eic::MCRecoParticleLinkCollection>();
0097   algo.process({&clusters, &cluster_assocs}, {parts.get(), part_links.get(), part_assocs.get()});
0098 #else
0099   algo.process({&clusters, &cluster_assocs}, {parts.get(), part_assocs.get()});
0100 #endif
0101
0102   // Two clusters in, two particles out
0103   REQUIRE(parts->size() == 2);
0104
0105   // --- First particle: from cluster1 at (0,0,100) with E=5 ---
0106   auto part1 = (*parts)[0];
0107   REQUIRE(part1.getPDG() == 22);
0108   REQUIRE_THAT(part1.getMass(), Catch::Matchers::WithinAbs(0.0, EPSILON));
0109   REQUIRE_THAT(part1.getCharge(), Catch::Matchers::WithinAbs(0.0, EPSILON));
0110   REQUIRE_THAT(part1.getEnergy(), Catch::Matchers::WithinAbs(5.0, EPSILON));
0111   // For massless particle |p| = E, directed along z
0112   REQUIRE_THAT(part1.getMomentum().x, Catch::Matchers::WithinAbs(0.0, EPSILON));
0113   REQUIRE_THAT(part1.getMomentum().y, Catch::Matchers::WithinAbs(0.0, EPSILON));
0114   REQUIRE_THAT(part1.getMomentum().z, Catch::Matchers::WithinAbs(5.0, EPSILON));
0115
0116   // --- Second particle: from cluster2 at (1,0,1) with E=10 ---
0117   auto part2 = (*parts)[1];
0118   REQUIRE_THAT(part2.getEnergy(), Catch::Matchers::WithinAbs(10.0, EPSILON));
0119   // Position (1,0,1) -> unit vector (1/sqrt2, 0, 1/sqrt2), |p| = E = 10
0120   double expected_comp = 10.0 / std::sqrt(2.0);
0121   REQUIRE_THAT(part2.getMomentum().x, Catch::Matchers::WithinAbs(expected_comp, 1e-4));
0122   REQUIRE_THAT(part2.getMomentum().y, Catch::Matchers::WithinAbs(0.0, EPSILON));
0123   REQUIRE_THAT(part2.getMomentum().z, Catch::Matchers::WithinAbs(expected_comp, 1e-4));
0124
0125   // --- Association: cluster1's MC association is propagated ---
0126   REQUIRE(part_assocs->size() == 1);
0127   auto assoc = (*part_assocs)[0];
0128   REQUIRE(assoc.getRec() == part1);
0129   REQUIRE(assoc.getSim() == mcpart);
0130   REQUIRE_THAT(assoc.getWeight(), Catch::Matchers::WithinAbs(0.9, EPSILON));
0131
0132 #if EDM4EIC_BUILD_VERSION >= EDM4EIC_VERSION(8, 7, 0)
0133   REQUIRE(part_links->size() == 1);
0134   auto link = (*part_links)[0];
0135   REQUIRE(link.getFrom() == part1);
0136   REQUIRE(link.getTo() == mcpart);
0137   REQUIRE_THAT(link.getWeight(), Catch::Matchers::WithinAbs(0.9, EPSILON));
0138 #endif
0139 }