tests/algorithms_test/calorimetry_CalorimeterClusterRecoCoG.cc

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