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 <edm4eic/TrackParametersCollection.h>
0005 #include <edm4eic/TrajectoryCollection.h>
0006 #include <edm4hep/MCParticleCollection.h>
0007 #include <edm4hep/Vector3d.h>
0008 #include <edm4hep/Vector3f.h>
0009 #include <edm4hep/utils/vector_utils.h>
0010 #include <fmt/core.h>
0011 #include <podio/ObjectID.h>
0012 #include <podio/RelationRange.h>
0013 #include <cmath>
0014 #include <gsl/pointers>
0015 #include <vector>
0016
0017 #include "TracksToParticles.h"
0018
0019
0020 namespace eicrecon {
0021
0022     void TracksToParticles::init() {}
0023
0024     void TracksToParticles::process(
0025       const TracksToParticles::Input& input, const TracksToParticles::Output& output
0026     ) const {
0027         const auto [tracks, track_assocs] = input;
0028         auto [parts, part_assocs]         = output;
0029
0030         for (const auto &track: *tracks) {
0031           auto trajectory = track.getTrajectory();
0032           for (const auto &trk: trajectory.getTrackParameters()) {
0033             const auto mom = edm4hep::utils::sphericalToVector(1.0 / std::abs(trk.getQOverP()), trk.getTheta(),
0034                                                         trk.getPhi());
0035             const auto charge_rec = std::copysign(1., trk.getQOverP());
0036
0037
0038             debug("Converting track: index={:<4} momentum={:<8.3f} theta={:<8.3f} phi={:<8.2f} charge={:<4}",
0039                   trk.getObjectID().index, edm4hep::utils::magnitude(mom), edm4hep::utils::anglePolar(mom), edm4hep::utils::angleAzimuthal(mom), charge_rec);
0040
0041             auto rec_part = parts->create();
0042             rec_part.addToTracks(track);
0043             rec_part.setType(0);
0044             rec_part.setEnergy(edm4hep::utils::magnitude(mom));
0045             rec_part.setMomentum(mom);
0046             rec_part.setCharge(charge_rec);
0047             rec_part.setMass(0.);
0048             rec_part.setGoodnessOfPID(0); // assume no PID until proven otherwise
0049             // rec_part.covMatrix()  // @TODO: covariance matrix on 4-momentum
0050
0051             double max_weight = -1.;
0052             for (auto track_assoc : *track_assocs) {
0053                 if (track_assoc.getRec() == track) {
0054                     trace("Found track association: index={} -> index={}, weight={}",
0055                           track_assoc.getRec().getObjectID().index,
0056                           track_assoc.getSim().getObjectID().index,
0057                           track_assoc.getWeight());
0058                     auto part_assoc = part_assocs->create();
0059                     part_assoc.setRec(rec_part);
0060                     part_assoc.setSim(track_assoc.getSim());
0061                     part_assoc.setRecID(part_assoc.getRec().getObjectID().index);
0062                     part_assoc.setSimID(part_assoc.getSim().getObjectID().index);
0063                     part_assoc.setWeight(track_assoc.getWeight());
0064
0065                     if (max_weight < track_assoc.getWeight()) {
0066                         max_weight = track_assoc.getWeight();
0067                         edm4hep::Vector3f referencePoint = {
0068                             static_cast<float>(track_assoc.getSim().getVertex().x),
0069                             static_cast<float>(track_assoc.getSim().getVertex().y),
0070                             static_cast<float>(track_assoc.getSim().getVertex().z)}; // @TODO: not sure if vertex/reference point makes sense here
0071                         rec_part.setReferencePoint(referencePoint);
0072                     }
0073                 }
0074             }
0075           }
0076         }
0077     }
0078 }