algorithms/tracking/ActsToTracks.cc

0001 // SPDX-License-Identifier: LGPL-3.0-or-later
0002 // Copyright (C) 2024 - 2025 Whitney Armstrong, Wouter Deconinck, Dmitry Romanov, Shujie Li, Dmitry Kalinkin
0003
0004 #include <Acts/Definitions/Algebra.hpp>
0005 #include <Acts/Definitions/TrackParametrization.hpp>
0006 #include <Acts/EventData/MultiTrajectoryHelpers.hpp>
0007 #include <Acts/EventData/ParticleHypothesis.hpp>
0008 #include <Acts/EventData/ProxyAccessor.hpp>
0009 #include <Acts/EventData/SourceLink.hpp>
0010 #include <Acts/EventData/TrackProxy.hpp>
0011 #include <Acts/EventData/VectorMultiTrajectory.hpp>
0012 #include <Acts/Geometry/GeometryContext.hpp>
0013 #include <Acts/Geometry/GeometryIdentifier.hpp>
0014 #include <Acts/Surfaces/Surface.hpp>
0015 #include <Acts/Utilities/UnitVectors.hpp>
0016 #include <ActsExamples/EventData/IndexSourceLink.hpp>
0017 #include <ActsExamples/EventData/Track.hpp>
0018 #include <edm4eic/Cov6f.h>
0019 #include <edm4eic/RawTrackerHit.h>
0020 #include <edm4eic/TrackerHit.h>
0021 #include <edm4hep/MCParticleCollection.h>
0022 #include <edm4hep/SimTrackerHit.h>
0023 #include <edm4hep/Vector2f.h>
0024 #include <edm4hep/Vector3f.h>
0025 #include <edm4hep/utils/vector_utils.h>
0026 #include <podio/ObjectID.h>
0027 #include <podio/RelationRange.h>
0028 #include <podio/detail/Link.h>
0029 #include <podio/detail/LinkCollectionImpl.h>
0030 #include <any>
0031 #include <array>
0032 #include <cmath>
0033 #include <cstddef>
0034 #include <map>
0035 #include <memory>
0036 #include <numeric>
0037 #include <tuple>
0038 #include <utility>
0039 #include <vector>
0040
0041 #include "ActsToTracks.h"
0042 #include "extensions/edm4eic/EDM4eicToActs.h"
0043
0044 namespace eicrecon {
0045
0046 // Custom comparator for MCParticle that uses deterministic ObjectID-based comparison
0047 // instead of podio's default memory-address-based comparison
0048 namespace {
0049   struct MCParticleCompare {
0050     bool operator()(const edm4hep::MCParticle& p_a, const edm4hep::MCParticle& p_b) const {
0051       // Compare particles by ObjectID for deterministic ordering
0052       auto id_a = p_a.getObjectID();
0053       auto id_b = p_b.getObjectID();
0054       if (id_a.collectionID != id_b.collectionID) {
0055         return id_a.collectionID < id_b.collectionID;
0056       }
0057       return id_a.index < id_b.index;
0058     }
0059   };
0060 } // namespace
0061
0062 void ActsToTracks::init() {}
0063
0064 void ActsToTracks::process(const Input& input, const Output& output) const {
0065   const auto [meas2Ds, track_seeds, acts_track_states, acts_tracks, raw_hit_assocs] = input;
0066   auto [trajectories, track_parameters, tracks, tracks_links, tracks_assoc]         = output;
0067
0068   // Create accessor for seed number dynamic column
0069   Acts::ConstProxyAccessor<unsigned int> seedNumber("seed");
0070
0071   // Construct ActsExamples::ConstTrackContainer from underlying containers
0072   auto trackStateContainer = std::make_shared<Acts::ConstVectorMultiTrajectory>(*acts_track_states);
0073   auto trackContainer      = std::make_shared<Acts::ConstVectorTrackContainer>(*acts_tracks);
0074   ActsExamples::ConstTrackContainer acts_track_container(trackContainer, trackStateContainer);
0075
0076   // Loop over tracks
0077   for (const auto& track : acts_track_container) {
0078     // Collect the trajectory summary info
0079     auto trajectoryState = Acts::MultiTrajectoryHelpers::trajectoryState(
0080         acts_track_container.trackStateContainer(), track.tipIndex());
0081
0082     // Create trajectory
0083     auto trajectory = trajectories->create();
0084     trajectory.setNMeasurements(trajectoryState.nMeasurements);
0085     trajectory.setNStates(trajectoryState.nStates);
0086     trajectory.setNOutliers(trajectoryState.nOutliers);
0087     trajectory.setNHoles(trajectoryState.nHoles);
0088     trajectory.setNSharedHits(trajectoryState.nSharedHits);
0089
0090     // Set the seed that was used to obtain this track
0091     unsigned int iseed = seedNumber(track);
0092     if (iseed < track_seeds->size()) {
0093       trajectory.setSeed((*track_seeds)[iseed]);
0094     } else {
0095       warning("ActsToTracks: seed index {} is out of bounds (track_seeds size = {}), seed will not "
0096               "be set for this trajectory",
0097               iseed, track_seeds->size());
0098     }
0099
0100     debug("trajectory state, measurement, outlier, hole: {} {} {} {}", trajectoryState.nStates,
0101           trajectoryState.nMeasurements, trajectoryState.nOutliers, trajectoryState.nHoles);
0102
0103     for (const auto& measurementChi2 : trajectoryState.measurementChi2) {
0104       trajectory.addToMeasurementChi2(measurementChi2);
0105     }
0106
0107     for (const auto& outlierChi2 : trajectoryState.outlierChi2) {
0108       trajectory.addToOutlierChi2(outlierChi2);
0109     }
0110
0111     // Get the fitted track parameter
0112     const auto& parameter  = track.parameters();
0113     const auto& covariance = track.covariance();
0114
0115     auto pars = track_parameters->create();
0116     pars.setType(0); // type: track head --> 0
0117     pars.setLoc({static_cast<float>(parameter[Acts::eBoundLoc0]),
0118                  static_cast<float>(parameter[Acts::eBoundLoc1])});
0119     pars.setTheta(static_cast<float>(parameter[Acts::eBoundTheta]));
0120     pars.setPhi(static_cast<float>(parameter[Acts::eBoundPhi]));
0121     pars.setQOverP(static_cast<float>(parameter[Acts::eBoundQOverP]));
0122     pars.setTime(static_cast<float>(parameter[Acts::eBoundTime]));
0123     edm4eic::Cov6f cov;
0124     for (std::size_t i = 0; const auto& [a, x] : edm4eic_indexed_units) {
0125       for (std::size_t j = 0; const auto& [b, y] : edm4eic_indexed_units) {
0126         // FIXME why not pars.getCovariance()(i,j) = covariance(a,b) / x / y;
0127         cov(i, j) = covariance(a, b) / x / y;
0128         ++j;
0129       }
0130       ++i;
0131     }
0132     pars.setCovariance(cov);
0133
0134     trajectory.addToTrackParameters(pars);
0135
0136     // Fill tracks
0137     auto track_out = tracks->create();
0138     track_out.setType( // Flag that defines the type of track
0139         pars.getType());
0140
0141     // Compute 3D position from perigee local coordinates via localToGlobal.
0142     // A default GeometryContext is sufficient here: the perigee surface is
0143     // defined purely by its center point and carries no alignment data.
0144     const Acts::Vector2 localPos{parameter[Acts::eBoundLoc0], parameter[Acts::eBoundLoc1]};
0145     const Acts::Vector3 direction =
0146         Acts::makeDirectionFromPhiTheta(parameter[Acts::eBoundPhi], parameter[Acts::eBoundTheta]);
0147     const Acts::Vector3 globalPos = track.referenceSurface().localToGlobal(
0148 #if Acts_VERSION_MAJOR >= 45
0149         Acts::GeometryContext::dangerouslyDefaultConstruct(),
0150 #else
0151         Acts::GeometryContext{},
0152 #endif
0153         localPos, direction);
0154     track_out.setPosition( // Track 3-position at the perigee [mm]
0155         edm4hep::Vector3f{static_cast<float>(globalPos.x()), static_cast<float>(globalPos.y()),
0156                           static_cast<float>(globalPos.z())});
0157
0158     // Compute Cartesian momentum from spherical parameters.
0159     const double qOverP = parameter[Acts::eBoundQOverP];
0160     double p_abs        = 0.0;
0161     if (std::isfinite(qOverP) && qOverP != 0.0) {
0162       p_abs = std::abs(1.0 / qOverP);
0163     } else {
0164       warning("ActsToTracks: track has qOverP={}, which yields non-finite momentum; setting "
0165               "momentum to zero",
0166               qOverP);
0167     }
0168     const double p = p_abs;
0169     track_out.setMomentum( // Track 3-momentum at the perigee [GeV]
0170         edm4hep::utils::sphericalToVector(p, parameter[Acts::eBoundTheta],
0171                                           parameter[Acts::eBoundPhi]));
0172
0173     track_out.setPositionMomentumCovariance( // Covariance matrix in basis [x,y,z,px,py,pz]
0174         edm4eic::Cov6f());
0175     track_out.setTime( // Track time at the perigee [ns]
0176         static_cast<float>(parameter[Acts::eBoundTime]));
0177     track_out.setTimeError( // Error on the track perigee time
0178         sqrt(static_cast<float>(covariance(Acts::eBoundTime, Acts::eBoundTime))));
0179     const double charge = // Particle charge (0 if qOverP is invalid or zero)
0180         (std::isfinite(qOverP) && qOverP != 0.0) ? std::copysign(1.0, qOverP) : 0.0;
0181     track_out.setCharge(charge);
0182     track_out.setChi2(trajectoryState.chi2Sum); // Total chi2
0183     track_out.setNdf(trajectoryState.NDF);      // Number of degrees of freedom
0184     track_out.setPdg(                           // PDG particle ID hypothesis
0185         track.particleHypothesis().absolutePdg());
0186     track_out.setTrajectory(trajectory); // Trajectory of this track
0187
0188     // Determine track association with MCParticle, weighted by number of used measurements
0189     std::map<edm4hep::MCParticle, double, MCParticleCompare> mcparticle_weight_by_hit_count;
0190
0191     // save measurement2d to good measurements or outliers according to srclink index
0192     // fix me: ideally, this should be integrated into multitrajectoryhelper
0193     // fix me: should say "OutlierMeasurements" instead of "OutlierHits" etc
0194     for (const auto& state : track.trackStatesReversed()) {
0195       auto geoID     = state.referenceSurface().geometryId().value();
0196       auto typeFlags = state.typeFlags();
0197
0198       // find the associated hit (2D measurement) with state sourcelink index
0199       // fix me: calibrated or not?
0200       if (state.hasUncalibratedSourceLink()) {
0201
0202         std::size_t srclink_index =
0203             state.getUncalibratedSourceLink().template get<ActsExamples::IndexSourceLink>().index();
0204
0205         // no hit on this state/surface, skip
0206 #if Acts_VERSION_MAJOR >= 45
0207         if (typeFlags.isHole()) {
0208 #else
0209         if (typeFlags.test(Acts::TrackStateFlag::HoleFlag)) {
0210 #endif
0211           debug("No hit found on geo id={}", geoID);
0212
0213         } else {
0214           auto meas2D = (*meas2Ds)[srclink_index];
0215 #if Acts_VERSION_MAJOR >= 45
0216           if (typeFlags.isOutlier()) {
0217 #else
0218           if (typeFlags.test(Acts::TrackStateFlag::OutlierFlag)) {
0219 #endif
0220             trajectory.addToOutliers_deprecated(meas2D);
0221             debug("Outlier on geo id={}, index={}, loc={},{}", geoID, srclink_index,
0222                   meas2D.getLoc().a, meas2D.getLoc().b);
0223 #if Acts_VERSION_MAJOR >= 45
0224           } else if (typeFlags.isMeasurement()) {
0225 #else
0226           } else if (typeFlags.test(Acts::TrackStateFlag::MeasurementFlag)) {
0227 #endif
0228             track_out.addToMeasurements(meas2D);
0229             trajectory.addToMeasurements_deprecated(meas2D);
0230             debug("Measurement on geo id={}, index={}, loc={},{}", geoID, srclink_index,
0231                   meas2D.getLoc().a, meas2D.getLoc().b);
0232
0233             // Determine track associations if hit associations provided
0234             // FIXME: not able to check whether optional inputs were provided
0235             //if (raw_hit_assocs->has_value()) {
0236             for (const auto& hit : meas2D.getHits()) {
0237               auto raw_hit = hit.getRawHit();
0238               for (const auto raw_hit_assoc : *raw_hit_assocs) {
0239                 if (raw_hit_assoc.getRawHit() == raw_hit) {
0240                   auto sim_hit     = raw_hit_assoc.getSimHit();
0241                   auto mc_particle = sim_hit.getParticle();
0242                   mcparticle_weight_by_hit_count[mc_particle]++;
0243                 }
0244               }
0245             }
0246             //}
0247           }
0248         }
0249       }
0250     }
0251
0252     // Store track associations if hit associations provided
0253     // FIXME: not able to check whether optional inputs were provided
0254     //if (raw_hit_assocs->has_value()) {
0255     double total_weight = std::accumulate(
0256         mcparticle_weight_by_hit_count.begin(), mcparticle_weight_by_hit_count.end(), 0,
0257         [](const double sum, const auto& i) { return sum + i.second; });
0258     for (const auto& [mcparticle, weight] : mcparticle_weight_by_hit_count) {
0259       double normalized_weight = weight / total_weight;
0260       auto track_link          = tracks_links->create();
0261       track_link.setFrom(track_out);
0262       track_link.setTo(mcparticle);
0263       track_link.setWeight(normalized_weight);
0264       auto track_assoc = tracks_assoc->create();
0265       track_assoc.setRec(track_out);
0266       track_assoc.setSim(mcparticle);
0267       track_assoc.setWeight(normalized_weight);
0268       debug("track {}: mcparticle {} weight {}", track_out.id().index, mcparticle.id().index,
0269             normalized_weight);
0270     }
0271     //}
0272   }
0273 }
0274
0275 } // namespace eicrecon