tests/algorithms_test/calorimetry_CalorimeterClusterRecoCoG.cc

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2024, Sebouh Paul
0003
0004
0005 #include <catch2/catch_test_macros.hpp>            // for AssertionHandler, operator""_catch_sr, StringRef, REQUIRE, operator<, operator==, operator>, TEST_CASE
0006 #include <catch2/matchers/catch_matchers.hpp>
0007 #include <catch2/matchers/catch_matchers_floating_point.hpp>
0008 #include <edm4eic/CalorimeterHitCollection.h>      // for CalorimeterHitCollection, MutableCalorimeterHit, CalorimeterHitMutableCollectionIterator
0009 #include <edm4eic/ClusterCollection.h>
0010 #include <edm4eic/EDM4eicVersion.h>
0011 #include <edm4hep/Vector2i.h>
0012 #include <edm4hep/Vector3d.h>
0013 #if EDM4EIC_VERSION_MAJOR >= 7
0014 #include <edm4eic/MCRecoCalorimeterHitAssociationCollection.h>
0015 #endif
0016 #include <edm4eic/MCRecoClusterParticleAssociationCollection.h>
0017 #include <edm4eic/ProtoClusterCollection.h>
0018 #include <edm4eic/unit_system.h>
0019 #include <edm4hep/CaloHitContributionCollection.h>
0020 #include <edm4hep/MCParticleCollection.h>
0021 #include <edm4hep/RawCalorimeterHitCollection.h>
0022 #include <edm4hep/SimCalorimeterHitCollection.h>
0023 #include <edm4hep/Vector3f.h>
0024 #include <math.h>
0025 #include <spdlog/common.h>                         // for level_enum
0026 #include <spdlog/logger.h>                         // for logger
0027 #include <spdlog/spdlog.h>                         // for default_logger
0028 #include <memory>                                  // for allocator, unique_ptr, make_unique, shared_ptr, __shared_ptr_access
0029 #include <tuple>
0030 #include <vector>
0031
0032 #include "algorithms/calorimetry/CalorimeterClusterRecoCoG.h"        // for CalorimeterClusterRecoCoG
0033 #include "algorithms/calorimetry/CalorimeterClusterRecoCoGConfig.h"  // for CalorimeterClusterRecoCoGConfig
0034
0035 using eicrecon::CalorimeterClusterRecoCoG;
0036 using eicrecon::CalorimeterClusterRecoCoGConfig;
0037
0038 using edm4eic::CalorimeterHit;
0039
0040 TEST_CASE( "the calorimeter CoG algorithm runs", "[CalorimeterClusterRecoCoG]" ) {
0041   const float EPSILON = 1e-5;
0042
0043   CalorimeterClusterRecoCoG algo("CalorimeterClusterRecoCoG");
0044
0045   std::shared_ptr<spdlog::logger> logger = spdlog::default_logger()->clone("CalorimeterClusterRecoCoG");
0046   logger->set_level(spdlog::level::trace);
0047
0048   CalorimeterClusterRecoCoGConfig cfg;
0049   cfg.energyWeight = "log";
0050   cfg.sampFrac = 0.0203;
0051   cfg.logWeightBaseCoeffs = {5.0,0.65,0.31};
0052   cfg.logWeightBase_Eref = 50 * edm4eic::unit::GeV;
0053   cfg.longitudinalShowerInfoAvailable = true;
0054
0055
0056   algo.applyConfig(cfg);
0057   algo.init();
0058
0059   edm4hep::RawCalorimeterHitCollection rawhits_coll;
0060   edm4eic::CalorimeterHitCollection hits_coll;
0061   edm4eic::ProtoClusterCollection pclust_coll;
0062   edm4hep::SimCalorimeterHitCollection simhits_coll;
0063 #if EDM4EIC_VERSION_MAJOR >= 7
0064   edm4eic::MCRecoCalorimeterHitAssociationCollection hitassocs_coll;
0065 #endif
0066   edm4hep::CaloHitContributionCollection contribs_coll;
0067   edm4hep::MCParticleCollection mcparts_coll;
0068   auto assoc_coll = std::make_unique<edm4eic::MCRecoClusterParticleAssociationCollection>();
0069   auto clust_coll = std::make_unique<edm4eic::ClusterCollection>();
0070
0071   //create a protocluster with 3 hits
0072   auto pclust = pclust_coll.create();
0073   edm4hep::Vector3f position({0, 0, 1 * edm4eic::unit::mm});
0074
0075   auto rawhit1 = rawhits_coll.create();
0076
0077   auto hit1 = hits_coll.create();
0078   hit1.setCellID(0);
0079   hit1.setEnergy(0.1 * edm4eic::unit::GeV);
0080   hit1.setEnergyError(0);
0081   hit1.setTime(0);
0082   hit1.setTimeError(0);
0083   hit1.setPosition(position);
0084   hit1.setDimension({0,0,0});
0085   hit1.setLocal(position);
0086 #if EDM4EIC_VERSION_MAJOR >= 7
0087   hit1.setRawHit(rawhit1);
0088 #endif
0089   pclust.addToHits(hit1);
0090   pclust.addToWeights(1);
0091
0092   auto mcpart11 = mcparts_coll.create(
0093     11, // std::int32_t PDG
0094     1, // std::int32_t generatorStatus
0095     0, // std::int32_t simulatorStatus
0096     -1., // float charge
0097     0., // float time
0098     0., // double mass
0099     edm4hep::Vector3d(), // edm4hep::Vector3d vertex
0100     edm4hep::Vector3d(), // edm4hep::Vector3d endpoint
0101     edm4hep::Vector3f(), // edm4hep::Vector3f momentum
0102     edm4hep::Vector3f(), // edm4hep::Vector3f momentumAtEndpoint
0103     edm4hep::Vector3f(), // edm4hep::Vector3f spin
0104     edm4hep::Vector2i() // edm4hep::Vector2i colorFlow
0105   );
0106
0107   auto mcpart12 = mcparts_coll.create(
0108     22, // std::int32_t PDG
0109     0, // std::int32_t generatorStatus
0110     (0x1 << edm4hep::MCParticle::BITCreatedInSimulation), // std::int32_t simulatorStatus
0111     0., // float charge
0112     0., // float time
0113     0., // double mass
0114     edm4hep::Vector3d(), // edm4hep::Vector3d vertex
0115     edm4hep::Vector3d(), // edm4hep::Vector3d endpoint
0116     edm4hep::Vector3f(), // edm4hep::Vector3f momentum
0117     edm4hep::Vector3f(), // edm4hep::Vector3f momentumAtEndpoint
0118     edm4hep::Vector3f(), // edm4hep::Vector3f spin
0119     edm4hep::Vector2i() // edm4hep::Vector2i colorFlow
0120   );
0121
0122   mcpart12.addToParents(mcpart11);
0123   mcpart11.addToDaughters(mcpart12);
0124
0125   auto contrib11 = contribs_coll.create(
0126     0, // int32_t PDG
0127     0.05 * edm4eic::unit::GeV, // float energy
0128     0.0, // float time
0129     edm4hep::Vector3f() // edm4hep::Vector3f stepPosition
0130   );
0131   contrib11.setParticle(mcpart11);
0132   auto contrib12 = contribs_coll.create(
0133     0, // int32_t PDG
0134     0.05 * edm4eic::unit::GeV, // float energy
0135     0.0, // float time
0136     edm4hep::Vector3f() // edm4hep::Vector3f stepPosition
0137   );
0138   contrib12.setParticle(mcpart12);
0139
0140   auto simhit1 = simhits_coll.create();
0141   simhit1.setCellID(hit1.getCellID());
0142   simhit1.setEnergy(0.1 * edm4eic::unit::GeV);
0143   simhit1.setPosition(hit1.getPosition());
0144   simhit1.addToContributions(contrib11);
0145   simhit1.addToContributions(contrib12);
0146
0147 #if EDM4EIC_VERSION_MAJOR >= 7
0148   auto hitassoc1 = hitassocs_coll->create();
0149   hitassoc1.setRawHit(rawhit1);
0150   hitassoc1.setSimHit(simhit1);
0151 #endif
0152
0153   auto rawhit2 = rawhits_coll.create();
0154
0155   position={-1 * edm4eic::unit::mm, 0, 2 * edm4eic::unit::mm};
0156   auto hit2 = hits_coll.create();
0157   hit2.setCellID(1);
0158   hit2.setEnergy(0.1 * edm4eic::unit::GeV);
0159   hit2.setEnergyError(0);
0160   hit2.setTime(0);
0161   hit2.setTimeError(0);
0162   hit2.setPosition(position);
0163   hit2.setDimension({0,0,0});
0164   hit2.setLocal(position);
0165 #if EDM4EIC_VERSION_MAJOR >= 7
0166   hit2.setRawHit(rawhit2);
0167 #endif
0168   pclust.addToHits(hit2);
0169   pclust.addToWeights(1);
0170
0171   auto mcpart2 = mcparts_coll.create(
0172     211, // std::int32_t PDG
0173     0, // std::int32_t generatorStatus
0174     (0x1 << edm4hep::MCParticle::BITCreatedInSimulation), // std::int32_t simulatorStatus
0175     0., // float charge
0176     0., // float time
0177     0., // double mass
0178     edm4hep::Vector3d(), // edm4hep::Vector3d vertex
0179     edm4hep::Vector3d(), // edm4hep::Vector3d endpoint
0180     edm4hep::Vector3f(), // edm4hep::Vector3f momentum
0181     edm4hep::Vector3f(), // edm4hep::Vector3f momentumAtEndpoint
0182     edm4hep::Vector3f(), // edm4hep::Vector3f spin
0183     edm4hep::Vector2i() // edm4hep::Vector2i colorFlow
0184   );
0185
0186   auto contrib2 = contribs_coll.create(
0187     0, // int32_t PDG
0188     0.1 * edm4eic::unit::GeV, // float energy
0189     0.0, // float time
0190     edm4hep::Vector3f() // edm4hep::Vector3f stepPosition
0191   );
0192   contrib2.setParticle(mcpart2);
0193
0194   auto simhit2 = simhits_coll.create();
0195   simhit2.setCellID(hit2.getCellID());
0196   simhit2.setEnergy(0.1 * edm4eic::unit::GeV);
0197   simhit2.setPosition(hit2.getPosition());
0198   simhit2.addToContributions(contrib2);
0199
0200 #if EDM4EIC_VERSION_MAJOR >= 7
0201   auto hitassoc2 = hitassocs_coll->create();
0202   hitassoc2.setRawHit(rawhit2);
0203   hitassoc2.setSimHit(simhit2);
0204 #endif
0205
0206  // Constructing input and output as per the algorithm's expected signature
0207 #if EDM4EIC_VERSION_MAJOR >= 7
0208   auto input = std::make_tuple(&pclust_coll, &hitassocs_coll);
0209 #else
0210   auto input = std::make_tuple(&pclust_coll, &simhits_coll);
0211 #endif
0212   auto output = std::make_tuple(clust_coll.get(), assoc_coll.get());
0213
0214   algo.process(input, output);
0215
0216   REQUIRE(clust_coll->size() == 1);
0217   auto clust = (*clust_coll)[0];
0218   REQUIRE_THAT(clust.getIntrinsicTheta(), Catch::Matchers::WithinAbs(M_PI / 4, EPSILON));
0219   // std::abs() checks if we land on -M_PI
0220   REQUIRE_THAT(std::abs(clust.getIntrinsicPhi()), Catch::Matchers::WithinAbs(M_PI, EPSILON));
0221
0222   REQUIRE(assoc_coll->size() == 2);
0223
0224   // Half of the energy comes from mcpart11 and its daughter mcpart12
0225   REQUIRE_THAT((*assoc_coll)[0].getWeight(), Catch::Matchers::WithinAbs(0.5, EPSILON));
0226   REQUIRE((*assoc_coll)[0].getRec() == clust);
0227   REQUIRE((*assoc_coll)[0].getSim() == mcpart11);
0228
0229   // Half of the energy comes from mcpart2
0230   REQUIRE_THAT((*assoc_coll)[1].getWeight(), Catch::Matchers::WithinAbs(0.5, EPSILON));
0231   REQUIRE((*assoc_coll)[1].getRec() == clust);
0232   REQUIRE((*assoc_coll)[1].getSim() == mcpart2);
0233
0234 }