algorithms/tracking/TracksToParticles.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 <algorithms/logger.h>
0005 #include <edm4eic/TrackParametersCollection.h>
0006 #include <edm4eic/TrajectoryCollection.h>
0007 #include <edm4hep/ParticleIDCollection.h>
0008 #include <edm4hep/Vector3d.h>
0009 #include <edm4hep/Vector3f.h>
0010 #include <edm4hep/utils/vector_utils.h>
0011 #include <fmt/core.h>
0012 #include <podio/ObjectID.h>
0013 #include <podio/RelationRange.h>
0014 #include <stdint.h>
0015 #include <cmath>
0016 #include <cstddef>
0017 #include <gsl/pointers>
0018 #include <limits>
0019 #include <mutex>
0020 #include <vector>
0021
0022 #include "TracksToParticles.h"
0023 #include "TracksToParticlesConfig.h"
0024
0025
0026 namespace eicrecon {
0027
0028     void TracksToParticles::init() {}
0029
0030     void TracksToParticles::process(
0031       const TracksToParticles::Input& input, const TracksToParticles::Output& output
0032     ) const {
0033         const auto [mc_particles, tracks] = input;
0034         auto [parts, assocs]              = output;
0035
0036         const double sinPhiOver2Tolerance = sin(0.5 * m_cfg.phiTolerance);
0037         tracePhiToleranceOnce(sinPhiOver2Tolerance, m_cfg.phiTolerance);
0038
0039         std::vector<bool> mc_prt_is_consumed(mc_particles->size(), false);         // MCParticle is already consumed flag
0040
0041         for (const auto &track: *tracks) {
0042           auto trajectory = track.getTrajectory();
0043           for (const auto &trk: trajectory.getTrackParameters()) {
0044             const auto mom = edm4hep::utils::sphericalToVector(1.0 / std::abs(trk.getQOverP()), trk.getTheta(),
0045                                                         trk.getPhi());
0046             const auto charge_rec = std::copysign(1., trk.getQOverP());
0047
0048
0049             debug("Match:  [id]   [mom]   [theta]  [phi]    [charge]  [PID]");
0050             debug(" Track : {:<4} {:<8.3f} {:<8.3f} {:<8.2f} {:<4}",
0051                          trk.getObjectID().index, edm4hep::utils::magnitude(mom), edm4hep::utils::anglePolar(mom), edm4hep::utils::angleAzimuthal(mom), charge_rec);
0052
0053             // utility variables for matching
0054             int best_match = -1;
0055             double best_delta = std::numeric_limits<double>::max();
0056             for (size_t ip = 0; ip < mc_particles->size(); ++ip) {
0057                 const auto &mc_part = (*mc_particles)[ip];
0058                 const auto &p = mc_part.getMomentum();
0059
0060                 trace("  MCParticle with id={:<4} mom={:<8.3f} charge={}", mc_part.getObjectID().index,
0061                              edm4hep::utils::magnitude(p), mc_part.getCharge());
0062
0063                 // Check if used
0064                 if (mc_prt_is_consumed[ip]) {
0065                     trace("    Ignoring. Particle is already used");
0066                     continue;
0067                 }
0068
0069                 // Check if non-primary
0070                 if (mc_part.getGeneratorStatus() > 1) {
0071                     trace("    Ignoring. GeneratorStatus > 1 => Non-primary particle");
0072                     continue;
0073                 }
0074
0075                 // Check if neutral
0076                 if (mc_part.getCharge() == 0) {
0077                     trace("    Ignoring. Neutral particle");
0078                     continue;
0079                 }
0080
0081                 // Check opposite charge
0082                 if (mc_part.getCharge() * charge_rec < 0) {
0083                     trace("    Ignoring. Opposite charge particle");
0084                     continue;
0085                 }
0086
0087                 const auto p_mag = edm4hep::utils::magnitude(p);
0088                 const auto p_phi = edm4hep::utils::angleAzimuthal(p);
0089                 const auto p_eta = edm4hep::utils::eta(p);
0090                 const double dp_rel = std::abs((edm4hep::utils::magnitude(mom) - p_mag) / p_mag);
0091                 // check the tolerance for sin(dphi/2) to avoid the hemisphere problem and allow
0092                 // for phi rollovers
0093                 const double dsphi = std::abs(sin(0.5 * (edm4hep::utils::angleAzimuthal(mom) - p_phi)));
0094                 const double deta = std::abs((edm4hep::utils::eta(mom) - p_eta));
0095
0096                 bool is_matching = dp_rel < m_cfg.momentumRelativeTolerance &&
0097                                    deta < m_cfg.etaTolerance &&
0098                                    dsphi < sinPhiOver2Tolerance;
0099
0100                 // Matching kinematics with the static variables doesn't work at low angles and within beam divergence
0101                 // TODO - Maybe reconsider variables used or divide into regions
0102                 // Backward going
0103                 if ((p_eta < -5) && (edm4hep::utils::eta(mom) < -5)) {
0104                   is_matching = true;
0105                 }
0106                 // Forward going
0107                 if ((p_eta >  5) && (edm4hep::utils::eta(mom) >  5)) {
0108                   is_matching = true;
0109                 }
0110
0111                 trace("    Decision: {}  dp: {:.4f} < {}  &&  d_eta: {:.6f} < {}  && d_sin_phi: {:.4e} < {:.4e} ",
0112                              is_matching? "Matching":"Ignoring",
0113                              dp_rel, m_cfg.momentumRelativeTolerance,
0114                              deta, m_cfg.etaTolerance,
0115                              dsphi, sinPhiOver2Tolerance);
0116
0117                 if (is_matching) {
0118                     const double delta =
0119                             std::hypot(dp_rel / m_cfg.momentumRelativeTolerance, deta / m_cfg.etaTolerance,
0120                                        dsphi / sinPhiOver2Tolerance);
0121                     if (delta < best_delta) {
0122                         best_match = ip;
0123                         best_delta = delta;
0124                         trace("    Is the best match now");
0125                     }
0126                 }
0127             }
0128             auto rec_part = parts->create();
0129             rec_part.addToTracks(track);
0130             auto referencePoint = rec_part.getReferencePoint();
0131             // float time          = 0;
0132             float mass = 0;
0133             if (best_match >= 0) {
0134                 trace("Best match is found and is: {}", best_match);
0135                 mc_prt_is_consumed[best_match] = true;
0136                 const auto &best_mc_part = (*mc_particles)[best_match];
0137                 referencePoint = {
0138                         static_cast<float>(best_mc_part.getVertex().x),
0139                         static_cast<float>(best_mc_part.getVertex().y),
0140                         static_cast<float>(best_mc_part.getVertex().z)}; // @TODO: not sure if vertex/reference point makes sense here
0141                 // time                 = mcpart.getTime();
0142                 mass = best_mc_part.getMass();
0143             }
0144
0145             rec_part.setType(static_cast<int16_t>(best_match >= 0 ? 0 : -1)); // @TODO: determine type codes
0146             rec_part.setEnergy((float) std::hypot(edm4hep::utils::magnitude(mom), mass));
0147             rec_part.setMomentum(mom);
0148             rec_part.setReferencePoint(referencePoint);
0149             rec_part.setCharge(charge_rec);
0150             rec_part.setMass(mass);
0151             rec_part.setGoodnessOfPID(0); // assume no PID until proven otherwise
0152             // rec_part.covMatrix()  // @TODO: covariance matrix on 4-momentum
0153
0154             // Also write MC <--> truth particle association if match was found
0155             if (best_match >= 0) {
0156                 auto rec_assoc = assocs->create();
0157                 rec_assoc.setRecID(rec_part.getObjectID().index);
0158                 rec_assoc.setSimID((*mc_particles)[best_match].getObjectID().index);
0159                 rec_assoc.setWeight(1);
0160                 rec_assoc.setRec(rec_part);
0161                 auto sim = (*mc_particles)[best_match];
0162                 rec_assoc.setSim(sim);
0163
0164                 if (level() <= algorithms::LogLevel::kDebug) {
0165
0166                     const auto &mcpart = (*mc_particles)[best_match];
0167                     const auto &p = mcpart.getMomentum();
0168                     const auto p_mag = edm4hep::utils::magnitude(p);
0169                     const auto p_phi = edm4hep::utils::angleAzimuthal(p);
0170                     const auto p_theta = edm4hep::utils::anglePolar(p);
0171                     debug(" MCPart: {:<4} {:<8.3f} {:<8.3f} {:<8.2f} {:<6}",
0172                                  mcpart.getObjectID().index, p_mag, p_theta, p_phi, mcpart.getCharge(),
0173                                  mcpart.getPDG());
0174
0175                     debug(" Assoc: id={} SimId={} RecId={}",
0176                                  rec_assoc.getObjectID().index, rec_assoc.getSim().getObjectID().index, rec_assoc.getSim().getObjectID().index);
0177
0178                     trace(" Assoc PDGs: sim.PDG | rec.PDG | rec.particleIDUsed.PDG = {:^6} | {:^6} | {:^6}",
0179                                  rec_assoc.getSim().getPDG(),
0180                                  rec_assoc.getRec().getPDG(),
0181                                  rec_assoc.getRec().getParticleIDUsed().isAvailable() ? rec_assoc.getRec().getParticleIDUsed().getPDG() : 0);
0182
0183
0184                 }
0185             }
0186             else {
0187                 debug(" MCPart: Did not find a good match");
0188             }
0189           }
0190         }
0191     }
0192
0193     void TracksToParticles::tracePhiToleranceOnce(const double sinPhiOver2Tolerance, double phiTolerance) const {
0194         // This function is called once to print tolerances useful for tracing
0195         static std::once_flag do_it_once;
0196         std::call_once(do_it_once, [this, sinPhiOver2Tolerance, phiTolerance]() {
0197             trace("m_cfg.phiTolerance: {:<8.4f} => sinPhiOver2Tolerance: {:<8.4f}", sinPhiOver2Tolerance, phiTolerance);
0198         });
0199     }
0200 }