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