Acts/Utilities/TrackHelpers.hpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2024 CERN for the benefit of the Acts project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
0008
0009 #pragma once
0010
0011 #include "Acts/Definitions/Tolerance.hpp"
0012 #include "Acts/EventData/MultiTrajectoryHelpers.hpp"
0013 #include "Acts/EventData/TrackParameters.hpp"
0014 #include "Acts/EventData/TrackStateType.hpp"
0015 #include "Acts/Geometry/GeometryContext.hpp"
0016 #include "Acts/Propagator/StandardAborters.hpp"
0017 #include "Acts/Surfaces/Surface.hpp"
0018 #include "Acts/TrackFitting/GainMatrixSmoother.hpp"
0019 #include "Acts/Utilities/Logger.hpp"
0020 #include "Acts/Utilities/Result.hpp"
0021
0022 #include <utility>
0023
0024 namespace Acts {
0025
0026 enum class TrackExtrapolationStrategy {
0027   /// Use the first track state to reach target surface
0028   first,
0029   /// Use the last track state to reach target surface
0030   last,
0031   /// Use the first or last track state to reach target surface depending on the
0032   /// distance
0033   firstOrLast,
0034 };
0035
0036 enum class TrackExtrapolationError {
0037   CompatibleTrackStateNotFound = 1,
0038   ReferenceSurfaceUnreachable = 2,
0039 };
0040
0041 std::error_code make_error_code(TrackExtrapolationError e);
0042
0043 template <typename track_proxy_t>
0044 Result<typename track_proxy_t::ConstTrackStateProxy> findFirstMeasurementState(
0045     const track_proxy_t &track) {
0046   using TrackStateProxy = typename track_proxy_t::ConstTrackStateProxy;
0047
0048   // TODO specialize if track is forward linked
0049
0050   auto result = Result<TrackStateProxy>::failure(
0051       TrackExtrapolationError::CompatibleTrackStateNotFound);
0052
0053   for (const auto &trackState : track.trackStatesReversed()) {
0054     bool isMeasurement =
0055         trackState.typeFlags().test(TrackStateFlag::MeasurementFlag);
0056     bool isOutlier = trackState.typeFlags().test(TrackStateFlag::OutlierFlag);
0057
0058     if (isMeasurement && !isOutlier) {
0059       result = trackState;
0060     }
0061   }
0062
0063   return result;
0064 }
0065
0066 template <typename track_proxy_t>
0067 Result<typename track_proxy_t::ConstTrackStateProxy> findLastMeasurementState(
0068     const track_proxy_t &track) {
0069   using TrackStateProxy = typename track_proxy_t::ConstTrackStateProxy;
0070
0071   for (const auto &trackState : track.trackStatesReversed()) {
0072     bool isMeasurement =
0073         trackState.typeFlags().test(TrackStateFlag::MeasurementFlag);
0074     bool isOutlier = trackState.typeFlags().test(TrackStateFlag::OutlierFlag);
0075
0076     if (isMeasurement && !isOutlier) {
0077       return trackState;
0078     }
0079   }
0080
0081   return Result<TrackStateProxy>::failure(
0082       TrackExtrapolationError::CompatibleTrackStateNotFound);
0083 }
0084
0085 /// @brief Smooth a track using the gain matrix smoother
0086 ///
0087 /// @tparam track_proxy_t The track proxy type
0088 /// @tparam smoother_t The smoother type
0089 ///
0090 /// @param geoContext The geometry context
0091 /// @param track The track to smooth
0092 /// @param logger The logger
0093 /// @param smoother The smoother
0094 ///
0095 /// @return The result of the smoothing
0096 template <typename track_proxy_t, typename smoother_t = GainMatrixSmoother>
0097 Result<void> smoothTrack(
0098     const GeometryContext &geoContext, track_proxy_t &track,
0099     const Logger &logger = *getDefaultLogger("TrackSmoother", Logging::INFO),
0100     smoother_t smoother = GainMatrixSmoother()) {
0101   auto &trackContainer = track.container();
0102   auto &trackStateContainer = trackContainer.trackStateContainer();
0103
0104   auto last = findLastMeasurementState(track);
0105   if (!last.ok()) {
0106     ACTS_ERROR("no last track state found");
0107     return last.error();
0108   }
0109
0110   auto smoothingResult =
0111       smoother(geoContext, trackStateContainer, last->index(), logger);
0112
0113   if (!smoothingResult.ok()) {
0114     ACTS_ERROR("Smoothing track " << track.index() << " failed with error "
0115                                   << smoothingResult.error());
0116     return smoothingResult.error();
0117   }
0118
0119   return Result<void>::success();
0120 }
0121
0122 /// @brief Smooth tracks using the gain matrix smoother
0123 ///
0124 /// @tparam track_container_t The track container type
0125 ///
0126 /// @param geoContext The geometry context
0127 /// @param trackContainer The track container
0128 /// @param logger The logger
0129 ///
0130 /// @return The result of the smoothing
0131 template <typename track_container_t>
0132 Result<void> smoothTracks(
0133     const GeometryContext &geoContext, const track_container_t &trackContainer,
0134     const Logger &logger = *getDefaultLogger("TrackSmoother", Logging::INFO)) {
0135   Result<void> result = Result<void>::success();
0136
0137   for (const auto &track : trackContainer) {
0138     auto smoothingResult = smoothTrack(geoContext, track, logger);
0139
0140     // Only keep the first error
0141     if (!smoothingResult.ok() && result.ok()) {
0142       result = smoothingResult.error();
0143     }
0144   }
0145
0146   return result;
0147 }
0148
0149 /// @brief Find a track state for extrapolation
0150 ///
0151 /// @tparam track_proxy_t The track proxy type
0152 ///
0153 /// @param geoContext The geometry context
0154 /// @param track The track
0155 /// @param referenceSurface The reference surface
0156 /// @param strategy The extrapolation strategy
0157 /// @param logger The logger
0158 ///
0159 /// @return The result of the search containing the track state
0160 ///         and the distance to the reference surface
0161 template <typename track_proxy_t>
0162 Result<std::pair<typename track_proxy_t::ConstTrackStateProxy, double>>
0163 findTrackStateForExtrapolation(
0164     const GeometryContext &geoContext, const track_proxy_t &track,
0165     const Surface &referenceSurface, TrackExtrapolationStrategy strategy,
0166     const Logger &logger = *getDefaultLogger("TrackExtrapolation",
0167                                              Logging::INFO)) {
0168   using TrackStateProxy = typename track_proxy_t::ConstTrackStateProxy;
0169
0170   auto intersect = [&](const TrackStateProxy &state) -> SurfaceIntersection {
0171     assert(state.hasSmoothed() || state.hasFiltered());
0172
0173     FreeVector freeVector;
0174     if (state.hasSmoothed()) {
0175       freeVector = MultiTrajectoryHelpers::freeSmoothed(geoContext, state);
0176     } else {
0177       freeVector = MultiTrajectoryHelpers::freeFiltered(geoContext, state);
0178     }
0179
0180     return referenceSurface
0181         .intersect(geoContext, freeVector.template segment<3>(eFreePos0),
0182                    freeVector.template segment<3>(eFreeDir0),
0183                    BoundaryTolerance::None(), s_onSurfaceTolerance)
0184         .closest();
0185   };
0186
0187   switch (strategy) {
0188     case TrackExtrapolationStrategy::first: {
0189       ACTS_VERBOSE("looking for first track state");
0190
0191       auto first = findFirstMeasurementState(track);
0192       if (!first.ok()) {
0193         ACTS_ERROR("no first track state found");
0194         return first.error();
0195       }
0196
0197       SurfaceIntersection intersection = intersect(*first);
0198       if (!intersection.isValid()) {
0199         ACTS_ERROR("no intersection found");
0200         return Result<std::pair<TrackStateProxy, double>>::failure(
0201             TrackExtrapolationError::ReferenceSurfaceUnreachable);
0202       }
0203
0204       ACTS_VERBOSE("found intersection at " << intersection.pathLength());
0205       return std::make_pair(*first, intersection.pathLength());
0206     }
0207
0208     case TrackExtrapolationStrategy::last: {
0209       ACTS_VERBOSE("looking for last track state");
0210
0211       auto last = findLastMeasurementState(track);
0212       if (!last.ok()) {
0213         ACTS_ERROR("no last track state found");
0214         return last.error();
0215       }
0216
0217       SurfaceIntersection intersection = intersect(*last);
0218       if (!intersection.isValid()) {
0219         ACTS_ERROR("no intersection found");
0220         return Result<std::pair<TrackStateProxy, double>>::failure(
0221             TrackExtrapolationError::ReferenceSurfaceUnreachable);
0222       }
0223
0224       ACTS_VERBOSE("found intersection at " << intersection.pathLength());
0225       return std::make_pair(*last, intersection.pathLength());
0226     }
0227
0228     case TrackExtrapolationStrategy::firstOrLast: {
0229       ACTS_VERBOSE("looking for first or last track state");
0230
0231       auto first = findFirstMeasurementState(track);
0232       if (!first.ok()) {
0233         ACTS_ERROR("no first track state found");
0234         return first.error();
0235       }
0236
0237       auto last = findLastMeasurementState(track);
0238       if (!last.ok()) {
0239         ACTS_ERROR("no last track state found");
0240         return last.error();
0241       }
0242
0243       SurfaceIntersection intersectionFirst = intersect(*first);
0244       SurfaceIntersection intersectionLast = intersect(*last);
0245
0246       double absDistanceFirst = std::abs(intersectionFirst.pathLength());
0247       double absDistanceLast = std::abs(intersectionLast.pathLength());
0248
0249       if (intersectionFirst.isValid() && absDistanceFirst <= absDistanceLast) {
0250         ACTS_VERBOSE("using first track state with intersection at "
0251                      << intersectionFirst.pathLength());
0252         return std::make_pair(*first, intersectionFirst.pathLength());
0253       }
0254
0255       if (intersectionLast.isValid() && absDistanceLast <= absDistanceFirst) {
0256         ACTS_VERBOSE("using last track state with intersection at "
0257                      << intersectionLast.pathLength());
0258         return std::make_pair(*last, intersectionLast.pathLength());
0259       }
0260
0261       ACTS_ERROR("no intersection found");
0262       return Result<std::pair<TrackStateProxy, double>>::failure(
0263           TrackExtrapolationError::ReferenceSurfaceUnreachable);
0264     }
0265   }
0266
0267   // unreachable
0268   return Result<std::pair<TrackStateProxy, double>>::failure(
0269       TrackExtrapolationError::CompatibleTrackStateNotFound);
0270 }
0271
0272 /// @brief Extrapolate a track to a reference surface
0273 ///
0274 /// @tparam track_proxy_t The track proxy type
0275 /// @tparam propagator_t The propagator type
0276 /// @tparam propagator_options_t The propagator options type
0277 ///
0278 /// @param track The track which is modified in-place
0279 /// @param referenceSurface The reference surface
0280 /// @param propagator The propagator
0281 /// @param options The propagator options
0282 /// @param strategy The extrapolation strategy
0283 /// @param logger The logger
0284 ///
0285 /// @return The result of the extrapolation
0286 template <typename track_proxy_t, typename propagator_t,
0287           typename propagator_options_t>
0288 Result<void> extrapolateTrackToReferenceSurface(
0289     track_proxy_t &track, const Surface &referenceSurface,
0290     const propagator_t &propagator, propagator_options_t options,
0291     TrackExtrapolationStrategy strategy,
0292     const Logger &logger = *getDefaultLogger("TrackExtrapolation",
0293                                              Logging::INFO)) {
0294   auto findResult = findTrackStateForExtrapolation(
0295       options.geoContext, track, referenceSurface, strategy, logger);
0296
0297   if (!findResult.ok()) {
0298     ACTS_ERROR("failed to find track state for extrapolation");
0299     return findResult.error();
0300   }
0301
0302   auto &[trackState, distance] = *findResult;
0303
0304   options.direction = Direction::fromScalarZeroAsPositive(distance);
0305
0306   BoundTrackParameters parameters = track.createParametersFromState(trackState);
0307   ACTS_VERBOSE("extrapolating track to reference surface at distance "
0308                << distance << " with direction " << options.direction
0309                << " with starting parameters " << parameters);
0310
0311   auto propagateResult =
0312       propagator.template propagate<BoundTrackParameters, propagator_options_t,
0313                                     ForcedSurfaceReached>(
0314           parameters, referenceSurface, options);
0315
0316   if (!propagateResult.ok()) {
0317     ACTS_ERROR("failed to extrapolate track: " << propagateResult.error());
0318     return propagateResult.error();
0319   }
0320
0321   track.setReferenceSurface(referenceSurface.getSharedPtr());
0322   track.parameters() = propagateResult->endParameters.value().parameters();
0323   track.covariance() =
0324       propagateResult->endParameters.value().covariance().value();
0325
0326   return Result<void>::success();
0327 }
0328
0329 /// @brief Extrapolate tracks to a reference surface
0330 ///
0331 /// @tparam track_container_t The track container type
0332 /// @tparam propagator_t The propagator type
0333 /// @tparam propagator_options_t The propagator options type
0334 ///
0335 /// @param trackContainer The track container which is modified in-place
0336 /// @param referenceSurface The reference surface
0337 /// @param propagator The propagator
0338 /// @param options The propagator options
0339 /// @param strategy The extrapolation strategy
0340 /// @param logger The logger
0341 ///
0342 /// @return The result of the extrapolation
0343 template <typename track_container_t, typename propagator_t,
0344           typename propagator_options_t>
0345 Result<void> extrapolateTracksToReferenceSurface(
0346     const track_container_t &trackContainer, const Surface &referenceSurface,
0347     const propagator_t &propagator, propagator_options_t options,
0348     TrackExtrapolationStrategy strategy,
0349     const Logger &logger = *getDefaultLogger("TrackExtrapolation",
0350                                              Logging::INFO)) {
0351   Result<void> result = Result<void>::success();
0352
0353   for (const auto &track : trackContainer) {
0354     auto extrapolateResult = extrapolateTrackToReferenceSurface(
0355         track, referenceSurface, propagator, options, strategy, logger);
0356
0357     // Only keep the first error
0358     if (!extrapolateResult.ok() && result.ok()) {
0359       result = extrapolateResult.error();
0360     }
0361   }
0362
0363   return result;
0364 }
0365
0366 /// Helper function to calculate a number of track level quantities and store
0367 /// them on the track itself
0368 /// @note The input track needs to be mutable, so @c ReadOnly=false
0369 /// @tparam track_container_t the track container backend
0370 /// @tparam track_state_container_t the track state container backend
0371 /// @tparam holder_t the holder type for the track container backends
0372 /// @param track A mutable track proxy to operate on
0373 template <typename track_container_t, typename track_state_container_t,
0374           template <typename> class holder_t>
0375 void calculateTrackQuantities(
0376     Acts::TrackProxy<track_container_t, track_state_container_t, holder_t,
0377                      false>
0378         track) {
0379   track.chi2() = 0;
0380   track.nDoF() = 0;
0381
0382   track.nHoles() = 0;
0383   track.nMeasurements() = 0;
0384   track.nSharedHits() = 0;
0385   track.nOutliers() = 0;
0386
0387   for (const auto &trackState : track.trackStatesReversed()) {
0388     auto typeFlags = trackState.typeFlags();
0389
0390     if (typeFlags.test(Acts::TrackStateFlag::HoleFlag)) {
0391       track.nHoles()++;
0392     } else if (typeFlags.test(Acts::TrackStateFlag::OutlierFlag)) {
0393       track.nOutliers()++;
0394     } else if (typeFlags.test(Acts::TrackStateFlag::MeasurementFlag)) {
0395       if (typeFlags.test(Acts::TrackStateFlag::SharedHitFlag)) {
0396         track.nSharedHits()++;
0397       }
0398       track.nMeasurements()++;
0399       track.chi2() += trackState.chi2();
0400       track.nDoF() += trackState.calibratedSize();
0401     }
0402   }
0403 }
0404
0405 }  // namespace Acts
0406
0407 namespace std {
0408 // register with STL
0409 template <>
0410 struct is_error_code_enum<Acts::TrackExtrapolationError> : std::true_type {};
0411 }  // namespace std