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