tests/algorithms_test/calorimetry_CalorimeterClusterShape.cc

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2024 - 2025, Sebouh Paul, Dmitry Kalinkin
0003
0004 #include <catch2/catch_test_macros.hpp>
0005 #include <catch2/matchers/catch_matchers.hpp>
0006 #include <catch2/matchers/catch_matchers_floating_point.hpp>
0007 #include <edm4eic/CalorimeterHitCollection.h>
0008 #include <edm4eic/ClusterCollection.h>
0009 #include <edm4eic/MCRecoClusterParticleAssociationCollection.h>
0010 #include <edm4eic/MCRecoClusterParticleLinkCollection.h>
0011 #include <edm4eic/unit_system.h>
0012 #include <edm4hep/MCParticleCollection.h>
0013 #include <edm4hep/Vector3f.h>
0014 #include <edm4hep/utils/vector_utils.h>
0015 #include <podio/detail/Link.h>
0016 #include <podio/detail/LinkCollectionImpl.h>
0017 #include <spdlog/common.h>
0018 #include <spdlog/logger.h>
0019 #include <spdlog/spdlog.h>
0020 #include <cmath>
0021 #include <deque>
0022 #include <memory>
0023 #include <string>
0024 #include <tuple>
0025
0026 #include "algorithms/calorimetry/CalorimeterClusterShape.h"
0027 #include "algorithms/calorimetry/CalorimeterClusterShapeConfig.h"
0028
0029 using eicrecon::CalorimeterClusterShape;
0030 using eicrecon::CalorimeterClusterShapeConfig;
0031
0032 TEST_CASE("the calorimeter CoG algorithm runs", "[CalorimeterClusterShape]") {
0033   const float EPSILON         = 1e-5;
0034   const float EXPECTED_WEIGHT = 0.123;
0035
0036   CalorimeterClusterShape algo("CalorimeterClusterShape");
0037
0038   std::shared_ptr<spdlog::logger> logger =
0039       spdlog::default_logger()->clone("CalorimeterClusterShape");
0040   logger->set_level(spdlog::level::trace);
0041
0042   CalorimeterClusterShapeConfig cfg;
0043   cfg.longitudinalShowerInfoAvailable = true;
0044
0045   algo.applyConfig(cfg);
0046   algo.init();
0047
0048   edm4eic::CalorimeterHitCollection hits_coll;
0049   edm4eic::MCRecoClusterParticleLinkCollection link_in_coll;
0050   edm4eic::ClusterCollection clust_in_coll;
0051   edm4hep::MCParticleCollection mcparts_coll;
0052   auto assoc_out_coll = std::make_unique<edm4eic::MCRecoClusterParticleAssociationCollection>();
0053   auto link_out_coll  = std::make_unique<edm4eic::MCRecoClusterParticleLinkCollection>();
0054   auto clust_out_coll = std::make_unique<edm4eic::ClusterCollection>();
0055
0056   auto hit1 = hits_coll.create();
0057   hit1.setCellID(0);
0058   hit1.setEnergy(0.1 * edm4eic::unit::GeV);
0059   hit1.setEnergyError(0);
0060   hit1.setTime(0);
0061   hit1.setTimeError(0);
0062   hit1.setPosition(edm4hep::Vector3f{0, 0, 1 * edm4eic::unit::mm});
0063   hit1.setDimension({0, 0, 0});
0064   hit1.setLocal(edm4hep::Vector3f{0, 0, 1 * edm4eic::unit::mm});
0065
0066   auto hit2 = hits_coll.create();
0067   hit2.setCellID(1);
0068   hit2.setEnergy(0.1 * edm4eic::unit::GeV);
0069   hit2.setEnergyError(0);
0070   hit2.setTime(0);
0071   hit2.setTimeError(0);
0072   hit2.setPosition(edm4hep::Vector3f{-1 * edm4eic::unit::mm, 0, 2 * edm4eic::unit::mm});
0073   hit2.setDimension({0, 0, 0});
0074   hit2.setLocal(edm4hep::Vector3f{-1 * edm4eic::unit::mm, 0, 2 * edm4eic::unit::mm});
0075
0076   // Create a cluster with 2 hits
0077   auto clust_in = clust_in_coll.create();
0078   clust_in.addToHits(hit1);
0079   clust_in.addToHitContributions(hit1.getEnergy());
0080   clust_in.addToHits(hit2);
0081   clust_in.addToHitContributions(hit2.getEnergy());
0082   clust_in.setNhits(clust_in.hits_size());
0083   clust_in.setEnergy(hit1.getEnergy() + hit2.getEnergy());
0084   clust_in.setPosition((hit1.getPosition() + hit2.getPosition()) / 2);
0085
0086   auto mcpart_in = mcparts_coll.create();
0087
0088   auto link_in = link_in_coll.create();
0089   link_in.setWeight(EXPECTED_WEIGHT);
0090   link_in.setFrom(clust_in);
0091   link_in.setTo(mcpart_in);
0092
0093   // Constructing input and output as per the algorithm's expected signature
0094   auto input  = std::make_tuple(&clust_in_coll, &link_in_coll);
0095   auto output = std::make_tuple(clust_out_coll.get(), link_out_coll.get(), assoc_out_coll.get());
0096
0097   algo.process(input, output);
0098
0099   REQUIRE(clust_out_coll->size() == 1);
0100   auto clust_out = (*clust_out_coll)[0];
0101   REQUIRE(clust_in.getNhits() == clust_out.getNhits());
0102
0103   REQUIRE_THAT(clust_out.getIntrinsicTheta(), Catch::Matchers::WithinAbs(M_PI / 4, EPSILON));
0104   // std::abs() checks if we land on -M_PI
0105   REQUIRE_THAT(std::abs(clust_out.getIntrinsicPhi()), Catch::Matchers::WithinAbs(M_PI, EPSILON));
0106
0107   REQUIRE(assoc_out_coll->size() == 1);
0108   REQUIRE((*assoc_out_coll)[0].getRec() == clust_out);
0109   REQUIRE((*assoc_out_coll)[0].getWeight() == link_in.getWeight());
0110
0111   // Validate links collection
0112   REQUIRE(link_out_coll->size() == 1);
0113
0114   // Check link from/to relationships - getFrom() should be the reconstructed cluster
0115   REQUIRE((*link_out_coll)[0].getFrom() == clust_out);
0116   REQUIRE((*link_out_coll)[0].getTo().isAvailable());
0117   REQUIRE((*link_out_coll)[0].getTo() == mcpart_in);
0118
0119   // Verify weight is propagated correctly
0120   REQUIRE((*link_out_coll)[0].getWeight() == EXPECTED_WEIGHT);
0121 }