tests/algorithms_test/pid_MergeParticleID.cc

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2023, Christopher Dilks
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   // charged particle input data
0031   //----------------------------------------------------------
0032   /* just create two unique tracks; they don't have to be filled with data,
0033    * since `MergeParticleID` matches tracks by ID
0034    */
0035   auto tracks = std::make_unique<edm4eic::TrackSegmentCollection>();
0036   for(int i=0; i<2; i++)
0037     tracks->create();
0038
0039   // Cherenkov PID inputs
0040   //----------------------------------------------------------
0041   /* 2 collections, each with 2 elements corresponding to 2 charged particles;
0042    * 2 or 3 possible PID hypotheses for each
0043    */
0044
0045   edm4eic::MutableCherenkovParticleID    pid;
0046   edm4eic::CherenkovParticleIDHypothesis pid_hyp;
0047
0048   // collection 1 ------
0049   auto coll_cherenkov_1 = std::make_unique<edm4eic::CherenkovParticleIDCollection>();
0050
0051   pid = coll_cherenkov_1->create();
0052   pid.setChargedParticle(tracks->at(0)); // track 0 first
0053   pid.setNpe(10);
0054   pid.setRefractiveIndex(1.05);
0055   pid.setPhotonEnergy(3e-9);
0056   pid_hyp.PDG    = 211; // pion hypothesis first
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; // kaon hypothesis first
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   // collection 2 ------
0084   auto coll_cherenkov_2 = std::make_unique<edm4eic::CherenkovParticleIDCollection>();
0085
0086   pid = coll_cherenkov_2->create();
0087   pid.setChargedParticle(tracks->at(1)); // track 1 first
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; // an additional hypothesis
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   // Cherenkov PID tests
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   // additive weights
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   // multiplicative weights
0203   /* NOTE: the only difference from additive weights is the resulting weight values, but it
0204    * is still wise to re-test everything, so that any changes in the algorithm require corresponding
0205    * changes in these tests
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 }