File indexing completed on 2024-09-27 07:03:08
0001
0002
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/CherenkovParticleIDCollection.h>
0008 #include <edm4eic/CherenkovParticleIDHypothesis.h>
0009 #include <edm4eic/TrackSegmentCollection.h>
0010 #include <edm4hep/Vector2f.h>
0011 #include <podio/RelationRange.h>
0012 #include <spdlog/common.h>
0013 #include <spdlog/logger.h>
0014 #include <spdlog/spdlog.h>
0015 #include <cmath>
0016 #include <memory>
0017 #include <stdexcept>
0018 #include <vector>
0019
0020 #include "algorithms/pid/MergeParticleID.h"
0021 #include "algorithms/pid/MergeParticleIDConfig.h"
0022
0023 TEST_CASE("the PID MergeParticleID algorithm runs", "[MergeParticleID]") {
0024
0025 const float EPSILON = 1e-5;
0026
0027 std::shared_ptr<spdlog::logger> logger = spdlog::default_logger()->clone("MergeParticleID");
0028 logger->set_level(spdlog::level::debug);
0029
0030
0031
0032
0033
0034
0035 auto tracks = std::make_unique<edm4eic::TrackSegmentCollection>();
0036 for(int i=0; i<2; i++)
0037 tracks->create();
0038
0039
0040
0041
0042
0043
0044
0045 edm4eic::MutableCherenkovParticleID pid;
0046 edm4eic::CherenkovParticleIDHypothesis pid_hyp;
0047
0048
0049 auto coll_cherenkov_1 = std::make_unique<edm4eic::CherenkovParticleIDCollection>();
0050
0051 pid = coll_cherenkov_1->create();
0052 pid.setChargedParticle(tracks->at(0));
0053 pid.setNpe(10);
0054 pid.setRefractiveIndex(1.05);
0055 pid.setPhotonEnergy(3e-9);
0056 pid_hyp.PDG = 211;
0057 pid_hyp.npe = 10;
0058 pid_hyp.weight = 100;
0059 pid.addToHypotheses(pid_hyp);
0060 pid_hyp.PDG = 321;
0061 pid_hyp.npe = 8;
0062 pid_hyp.weight = 80;
0063 pid.addToHypotheses(pid_hyp);
0064 for(int i=0; i<=pid.getNpe(); i++)
0065 pid.addToThetaPhiPhotons(edm4hep::Vector2f{M_PI/4, 0});
0066
0067 pid = coll_cherenkov_1->create();
0068 pid.setChargedParticle(tracks->at(1));
0069 pid.setNpe(11);
0070 pid.setRefractiveIndex(1.06);
0071 pid.setPhotonEnergy(4e-9);
0072 pid_hyp.PDG = 321;
0073 pid_hyp.npe = 10;
0074 pid_hyp.weight = 100;
0075 pid.addToHypotheses(pid_hyp);
0076 pid_hyp.PDG = 211;
0077 pid_hyp.npe = 10;
0078 pid_hyp.weight = 80;
0079 pid.addToHypotheses(pid_hyp);
0080 for(int i=0; i<=pid.getNpe(); i++)
0081 pid.addToThetaPhiPhotons(edm4hep::Vector2f{-M_PI/4, 0});
0082
0083
0084 auto coll_cherenkov_2 = std::make_unique<edm4eic::CherenkovParticleIDCollection>();
0085
0086 pid = coll_cherenkov_2->create();
0087 pid.setChargedParticle(tracks->at(1));
0088 pid.setNpe(4);
0089 pid.setRefractiveIndex(1.5);
0090 pid.setPhotonEnergy(3e-9);
0091 pid_hyp.PDG = 211;
0092 pid_hyp.npe = 3;
0093 pid_hyp.weight = 200;
0094 pid.addToHypotheses(pid_hyp);
0095 pid_hyp.PDG = 321;
0096 pid_hyp.npe = 2;
0097 pid_hyp.weight = 90;
0098 pid.addToHypotheses(pid_hyp);
0099 for(int i=0; i<=pid.getNpe(); i++)
0100 pid.addToThetaPhiPhotons(edm4hep::Vector2f{M_PI/4, 0});
0101
0102 pid = coll_cherenkov_2->create();
0103 pid.setChargedParticle(tracks->at(0));
0104 pid.setNpe(6);
0105 pid.setRefractiveIndex(1.3);
0106 pid.setPhotonEnergy(4e-9);
0107 pid_hyp.PDG = 321;
0108 pid_hyp.npe = 4;
0109 pid_hyp.weight = 70;
0110 pid.addToHypotheses(pid_hyp);
0111 pid_hyp.PDG = 211;
0112 pid_hyp.npe = 1;
0113 pid_hyp.weight = 80;
0114 pid.addToHypotheses(pid_hyp);
0115 pid_hyp.PDG = 11;
0116 pid_hyp.npe = 1;
0117 pid_hyp.weight = 60;
0118 pid.addToHypotheses(pid_hyp);
0119 for(int i=0; i<=pid.getNpe(); i++)
0120 pid.addToThetaPhiPhotons(edm4hep::Vector2f{-M_PI/4, 0});
0121
0122
0123
0124
0125 std::vector<gsl::not_null<const edm4eic::CherenkovParticleIDCollection*>> coll_cherenkov_list = {
0126 coll_cherenkov_1.get(),
0127 coll_cherenkov_2.get()
0128 };
0129
0130 auto find_cherenkov_pid_for_track = [] (const auto& coll, auto trk) {
0131 for(auto obj : *coll) {
0132 if(obj.getChargedParticle().id() == trk.id())
0133 return obj;
0134 }
0135 FAIL("ERROR: cannot find CherenkovParticleID given track");
0136 if(coll->size() == 0)
0137 throw std::runtime_error("empty collection used in pid_MergeParticleID::find_cherenkov_pid_for_track");
0138 return coll->at(0);
0139 };
0140
0141
0142 SECTION("merge CherenkovParticleID: add hypothesis weights") {
0143
0144 eicrecon::MergeParticleID algo("test");
0145 eicrecon::MergeParticleIDConfig cfg;
0146 cfg.mergeMode = eicrecon::MergeParticleIDConfig::kAddWeights;
0147 algo.applyConfig(cfg);
0148 algo.init(logger);
0149
0150 auto result = std::make_unique<edm4eic::CherenkovParticleIDCollection>();
0151 algo.process({coll_cherenkov_list}, {result.get()});
0152 auto pid_0 = find_cherenkov_pid_for_track(result, tracks->at(0));
0153 auto pid_1 = find_cherenkov_pid_for_track(result, tracks->at(1));
0154
0155 REQUIRE_THAT( pid_0.getNpe(), Catch::Matchers::WithinAbs(10+6, EPSILON) );
0156 REQUIRE_THAT( pid_1.getNpe(), Catch::Matchers::WithinAbs(11+4, EPSILON) );
0157
0158 REQUIRE_THAT( pid_0.getRefractiveIndex(), Catch::Matchers::WithinAbs((10*1.05+6*1.3)/(10+6), EPSILON) );
0159 REQUIRE_THAT( pid_1.getRefractiveIndex(), Catch::Matchers::WithinAbs((11*1.06+4*1.5)/(11+4), EPSILON) );
0160
0161 REQUIRE_THAT( pid_0.getPhotonEnergy(), Catch::Matchers::WithinAbs((10*3e-9+6*4e-9)/(10+6), EPSILON) );
0162 REQUIRE_THAT( pid_1.getPhotonEnergy(), Catch::Matchers::WithinAbs((11*4e-9+4*3e-9)/(11+4), EPSILON) );
0163
0164 REQUIRE(pid_0.hypotheses_size() == 3);
0165 for(auto hyp : pid_0.getHypotheses()) {
0166 switch(hyp.PDG) {
0167 case 211:
0168 REQUIRE_THAT( hyp.npe, Catch::Matchers::WithinAbs(10+1, EPSILON) );
0169 REQUIRE_THAT( hyp.weight, Catch::Matchers::WithinAbs(100+80, EPSILON) );
0170 break;
0171 case 321:
0172 REQUIRE_THAT( hyp.npe, Catch::Matchers::WithinAbs(8+4, EPSILON) );
0173 REQUIRE_THAT( hyp.weight, Catch::Matchers::WithinAbs(80+70, EPSILON) );
0174 break;
0175 case 11:
0176 REQUIRE_THAT( hyp.npe, Catch::Matchers::WithinAbs(1, EPSILON) );
0177 REQUIRE_THAT( hyp.weight, Catch::Matchers::WithinAbs(60, EPSILON) );
0178 break;
0179 default:
0180 FAIL("untested PDG hypothesis");
0181 }
0182 }
0183
0184 REQUIRE(pid_1.hypotheses_size() == 2);
0185 for(auto hyp : pid_1.getHypotheses()) {
0186 switch(hyp.PDG) {
0187 case 211:
0188 REQUIRE_THAT( hyp.npe, Catch::Matchers::WithinAbs(10+3, EPSILON) );
0189 REQUIRE_THAT( hyp.weight, Catch::Matchers::WithinAbs(80+200, EPSILON) );
0190 break;
0191 case 321:
0192 REQUIRE_THAT( hyp.npe, Catch::Matchers::WithinAbs(10+2, EPSILON) );
0193 REQUIRE_THAT( hyp.weight, Catch::Matchers::WithinAbs(100+90, EPSILON) );
0194 break;
0195 default:
0196 FAIL("untested PDG hypothesis");
0197 }
0198 }
0199
0200 }
0201
0202
0203
0204
0205
0206
0207 SECTION("merge CherenkovParticleID: multiply hypothesis weights") {
0208
0209 eicrecon::MergeParticleID algo("test");
0210 eicrecon::MergeParticleIDConfig cfg;
0211 cfg.mergeMode = eicrecon::MergeParticleIDConfig::kMultiplyWeights;
0212 algo.applyConfig(cfg);
0213 algo.init(logger);
0214
0215 auto result = std::make_unique<edm4eic::CherenkovParticleIDCollection>();
0216 algo.process({coll_cherenkov_list}, {result.get()});
0217 auto pid_0 = find_cherenkov_pid_for_track(result, tracks->at(0));
0218 auto pid_1 = find_cherenkov_pid_for_track(result, tracks->at(1));
0219
0220 REQUIRE_THAT( pid_0.getNpe(), Catch::Matchers::WithinAbs(10+6, EPSILON) );
0221 REQUIRE_THAT( pid_1.getNpe(), Catch::Matchers::WithinAbs(11+4, EPSILON) );
0222
0223 REQUIRE_THAT( pid_0.getRefractiveIndex(), Catch::Matchers::WithinAbs((10*1.05+6*1.3)/(10+6), EPSILON) );
0224 REQUIRE_THAT( pid_1.getRefractiveIndex(), Catch::Matchers::WithinAbs((11*1.06+4*1.5)/(11+4), EPSILON) );
0225
0226 REQUIRE_THAT( pid_0.getPhotonEnergy(), Catch::Matchers::WithinAbs((10*3e-9+6*4e-9)/(10+6), EPSILON) );
0227 REQUIRE_THAT( pid_1.getPhotonEnergy(), Catch::Matchers::WithinAbs((11*4e-9+4*3e-9)/(11+4), EPSILON) );
0228
0229 for(auto hyp : pid_0.getHypotheses()) {
0230 switch(hyp.PDG) {
0231 case 211:
0232 REQUIRE_THAT( hyp.npe, Catch::Matchers::WithinAbs(10+1, EPSILON) );
0233 REQUIRE_THAT( hyp.weight, Catch::Matchers::WithinAbs(100*80, EPSILON) );
0234 break;
0235 case 321:
0236 REQUIRE_THAT( hyp.npe, Catch::Matchers::WithinAbs(8+4, EPSILON) );
0237 REQUIRE_THAT( hyp.weight, Catch::Matchers::WithinAbs(80*70, EPSILON) );
0238 break;
0239 case 11:
0240 REQUIRE_THAT( hyp.npe, Catch::Matchers::WithinAbs(1, EPSILON) );
0241 REQUIRE_THAT( hyp.weight, Catch::Matchers::WithinAbs(60, EPSILON) );
0242 break;
0243 default:
0244 FAIL("untested PDG hypothesis");
0245 }
0246 }
0247
0248 for(auto hyp : pid_1.getHypotheses()) {
0249 switch(hyp.PDG) {
0250 case 211:
0251 REQUIRE_THAT( hyp.npe, Catch::Matchers::WithinAbs(10+3, EPSILON) );
0252 REQUIRE_THAT( hyp.weight, Catch::Matchers::WithinAbs(80*200, EPSILON) );
0253 break;
0254 case 321:
0255 REQUIRE_THAT( hyp.npe, Catch::Matchers::WithinAbs(10+2, EPSILON) );
0256 REQUIRE_THAT( hyp.weight, Catch::Matchers::WithinAbs(100*90, EPSILON) );
0257 break;
0258 default:
0259 FAIL("untested PDG hypothesis");
0260 }
0261 }
0262
0263 }
0264
0265 }