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