algorithms/pid/MatchToRICHPID.cc

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2022 - 2024, Sylvester Joosten, Wouter Deconinck, Dmitry Romanov, Christopher Dilks, Dmitry Kalinkin
0003
0004 #include "MatchToRICHPID.h"
0005
0006 #include <edm4eic/EDM4eicVersion.h>
0007 #include <edm4eic/TrackPoint.h>
0008 #include <edm4eic/TrackSegmentCollection.h>
0009 #include <edm4hep/MCParticle.h>
0010 #include <edm4hep/Vector3f.h>
0011 #include <edm4hep/utils/vector_utils.h>
0012 #include <podio/ObjectID.h>
0013 #include <podio/RelationRange.h>
0014 #include <podio/detail/Link.h>
0015 #include <podio/detail/LinkCollectionImpl.h>
0016 #include <algorithm>
0017 #include <cmath>
0018 #include <cstddef>
0019 #include <gsl/pointers>
0020 #include <map>
0021 #include <memory>
0022 #include <vector>
0023
0024 #include "algorithms/pid/ConvertParticleID.h"
0025 #include "algorithms/pid/MatchToRICHPIDConfig.h"
0026
0027 namespace eicrecon {
0028
0029 void MatchToRICHPID::init() {}
0030
0031 void MatchToRICHPID::process(const MatchToRICHPID::Input& input,
0032                              const MatchToRICHPID::Output& output) const {
0033   const auto [parts_in, assocs_in, drich_cherenkov_pid] = input;
0034 #if EDM4EIC_BUILD_VERSION >= EDM4EIC_VERSION(8, 7, 0)
0035   auto [parts_out, links_out, assocs_out, pids] = output;
0036 #else
0037   auto [parts_out, assocs_out, pids] = output;
0038 #endif
0039
0040   for (auto part_in : *parts_in) {
0041     auto part_out = part_in.clone();
0042
0043     // link Cherenkov PID objects
0044     auto success = linkCherenkovPID(part_out, *drich_cherenkov_pid, *pids);
0045     if (success) {
0046       trace("Previous PDG vs. CherenkovPID PDG: {:>10} vs. {:<10}", part_in.getPDG(),
0047             part_out.getParticleIDUsed().isAvailable() ? part_out.getParticleIDUsed().getPDG() : 0);
0048     }
0049
0050     for (auto assoc_in : *assocs_in) {
0051       if (assoc_in.getRec() == part_in) {
0052 #if EDM4EIC_BUILD_VERSION >= EDM4EIC_VERSION(8, 7, 0)
0053         auto link_out = links_out->create();
0054         link_out.setFrom(part_out);
0055         link_out.setTo(assoc_in.getSim());
0056         link_out.setWeight(assoc_in.getWeight());
0057 #endif
0058         auto assoc_out = assoc_in.clone();
0059         assoc_out.setRec(part_out);
0060         assocs_out->push_back(assoc_out);
0061       }
0062     }
0063
0064     parts_out->push_back(part_out);
0065   }
0066 }
0067
0068 /* link PID objects to input particle
0069      * - finds `CherenkovParticleID` object in `in_pids` associated to particle `in_part`
0070      *   by proximity matching to the associated track
0071      * - converts this `CherenkovParticleID` object's PID hypotheses to `ParticleID` objects,
0072      *   relates them to `in_part`, and adds them to the collection `out_pids` for persistency
0073      * - returns `true` iff PID objects were found and linked
0074      */
0075 bool MatchToRICHPID::linkCherenkovPID(edm4eic::MutableReconstructedParticle& in_part,
0076                                       const edm4eic::CherenkovParticleIDCollection& in_pids,
0077                                       edm4hep::ParticleIDCollection& out_pids) const {
0078
0079   // skip this particle, if neutral
0080   if (std::abs(in_part.getCharge()) < 0.001) {
0081     return false;
0082   }
0083
0084   // structure to store list of candidate matches
0085   struct ProxMatch {
0086     double match_dist;
0087     std::size_t pid_idx;
0088   };
0089   std::vector<ProxMatch> prox_match_list;
0090
0091   // get input reconstructed particle momentum angles
0092   auto in_part_p   = in_part.getMomentum();
0093   auto in_part_eta = edm4hep::utils::eta(in_part_p);
0094   auto in_part_phi = edm4hep::utils::angleAzimuthal(in_part_p);
0095   trace("Input particle: (eta,phi) = ( {:>5.4}, {:>5.4} deg )", in_part_eta,
0096         in_part_phi * 180.0 / M_PI);
0097
0098   // loop over input CherenkovParticleID objects
0099   for (std::size_t in_pid_idx = 0; in_pid_idx < in_pids.size(); in_pid_idx++) {
0100     auto in_pid = in_pids.at(in_pid_idx);
0101
0102     // get charged particle track associated to this CherenkovParticleID object
0103     auto in_track = in_pid.getChargedParticle();
0104     if (!in_track.isAvailable()) {
0105       error("found CherenkovParticleID object with no chargedParticle");
0106       return false;
0107     }
0108     if (in_track.points_size() == 0) {
0109       error("found chargedParticle for CherenkovParticleID, but it has no TrackPoints");
0110       return false;
0111     }
0112
0113     // get average momentum direction of the track's TrackPoints
0114     decltype(edm4eic::TrackPoint::momentum) in_track_p{0.0, 0.0, 0.0};
0115     for (const auto& in_track_point : in_track.getPoints()) {
0116       in_track_p = in_track_p + (in_track_point.momentum / in_track.points_size());
0117     }
0118     auto in_track_eta = edm4hep::utils::eta(in_track_p);
0119     auto in_track_phi = edm4hep::utils::angleAzimuthal(in_track_p);
0120
0121     // calculate dist(eta,phi)
0122     auto match_dist = std::hypot(in_part_eta - in_track_eta, in_part_phi - in_track_phi);
0123
0124     // check if the match is close enough: within user-specified tolerances
0125     auto match_is_close = std::abs(in_part_eta - in_track_eta) < m_cfg.etaTolerance &&
0126                           std::abs(in_part_phi - in_track_phi) < m_cfg.phiTolerance;
0127     if (match_is_close) {
0128       prox_match_list.push_back(ProxMatch{.match_dist = match_dist, .pid_idx = in_pid_idx});
0129     }
0130
0131     // logging
0132     trace("  - (eta,phi) = ( {:>5.4}, {:>5.4} deg ),  match_dist = {:<5.4}{}", in_track_eta,
0133           in_track_phi * 180.0 / M_PI, match_dist, match_is_close ? " => CLOSE!" : "");
0134
0135   } // end loop over input CherenkovParticleID objects
0136
0137   // check if at least one match was found
0138   if (prox_match_list.empty()) {
0139     trace("  => no matching CherenkovParticleID found for this particle");
0140     return false;
0141   }
0142
0143   // choose the closest matching CherenkovParticleID object corresponding to this input reconstructed particle
0144   auto closest_prox_match = *std::ranges::min_element(
0145       prox_match_list, [](ProxMatch a, ProxMatch b) { return a.match_dist < b.match_dist; });
0146   auto in_pid_matched = in_pids.at(closest_prox_match.pid_idx);
0147   trace("  => best match: match_dist = {:<5.4} at idx = {}", closest_prox_match.match_dist,
0148         closest_prox_match.pid_idx);
0149
0150   // convert `CherenkovParticleID` object's hypotheses => set of `ParticleID` objects
0151   auto out_pid_index_map = ConvertParticleID::ConvertToParticleIDs(in_pid_matched, out_pids, true);
0152   if (out_pid_index_map.empty()) {
0153     error("found CherenkovParticleID object with no hypotheses");
0154     return false;
0155   }
0156
0157   // relate matched ParticleID objects to output particle
0158   for (const auto& [out_pids_index, out_pids_id] : out_pid_index_map) {
0159     const auto& out_pid = out_pids.at(out_pids_index);
0160     if (out_pid.getObjectID().index != static_cast<int>(out_pids_id)) { // sanity check
0161       error("indexing error in `edm4eic::ParticleID` collection");
0162       return false;
0163     }
0164     in_part.addToParticleIDs(out_pid);
0165   }
0166   in_part.setParticleIDUsed(in_part.getParticleIDs().at(0)); // highest likelihood is the first
0167   in_part.setGoodnessOfPID(1); // FIXME: not used yet, aside from 0=noPID vs 1=hasPID
0168
0169   // trace logging
0170   trace("    {:.^50}", " PID results ");
0171   trace("      Hypotheses (sorted):");
0172   for (auto hyp : in_part.getParticleIDs()) {
0173     float npe =
0174         hyp.parameters_size() > 0 ? hyp.getParameters(0) : -1; // assume NPE is the first parameter
0175     trace("{:{}}{:>6}  {:>10.8}  {:>10.8}", "", 8, hyp.getPDG(), hyp.getLikelihood(), npe);
0176   }
0177   trace("    {:'^50}", "");
0178
0179   return true;
0180 }
0181 } // namespace eicrecon