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