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