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 ///
0089 /// @param geoContext The geometry context
0090 /// @param track The track to smooth
0091 /// @param logger The logger
0092 ///
0093 /// @return The result of the smoothing
0094 template <typename track_proxy_t>
0095 Result<void> smoothTrack(
0096     const GeometryContext &geoContext, track_proxy_t &track,
0097     const Logger &logger = *getDefaultLogger("TrackSmoother", Logging::INFO)) {
0098   Acts::GainMatrixSmoother smoother;
0099
0100   auto &trackContainer = track.container();
0101   auto &trackStateContainer = trackContainer.trackStateContainer();
0102
0103   auto last = findLastMeasurementState(track);
0104   if (!last.ok()) {
0105     ACTS_ERROR("no last track state found");
0106     return last.error();
0107   }
0108
0109   auto smoothingResult =
0110       smoother(geoContext, trackStateContainer, last->index(), logger);
0111
0112   if (!smoothingResult.ok()) {
0113     ACTS_ERROR("Smoothing track " << track.index() << " failed with error "
0114                                   << smoothingResult.error());
0115     return smoothingResult.error();
0116   }
0117
0118   return Result<void>::success();
0119 }
0120
0121 /// @brief Smooth tracks using the gain matrix smoother
0122 ///
0123 /// @tparam track_container_t The track container type
0124 ///
0125 /// @param geoContext The geometry context
0126 /// @param trackContainer The track container
0127 /// @param logger The logger
0128 ///
0129 /// @return The result of the smoothing
0130 template <typename track_container_t>
0131 Result<void> smoothTracks(
0132     const GeometryContext &geoContext, const track_container_t &trackContainer,
0133     const Logger &logger = *getDefaultLogger("TrackSmoother", Logging::INFO)) {
0134   Result<void> result = Result<void>::success();
0135
0136   for (const auto &track : trackContainer) {
0137     auto smoothingResult = smoothTrack(geoContext, track, logger);
0138
0139     // Only keep the first error
0140     if (!smoothingResult.ok() && result.ok()) {
0141       result = smoothingResult.error();
0142     }
0143   }
0144
0145   return result;
0146 }
0147
0148 /// @brief Find a track state for extrapolation
0149 ///
0150 /// @tparam track_proxy_t The track proxy type
0151 ///
0152 /// @param geoContext The geometry context
0153 /// @param track The track
0154 /// @param referenceSurface The reference surface
0155 /// @param strategy The extrapolation strategy
0156 /// @param logger The logger
0157 ///
0158 /// @return The result of the search containing the track state
0159 ///         and the distance to the reference surface
0160 template <typename track_proxy_t>
0161 Result<std::pair<typename track_proxy_t::ConstTrackStateProxy, double>>
0162 findTrackStateForExtrapolation(
0163     const GeometryContext &geoContext, const track_proxy_t &track,
0164     const Surface &referenceSurface, TrackExtrapolationStrategy strategy,
0165     const Logger &logger = *getDefaultLogger("TrackExtrapolation",
0166                                              Logging::INFO)) {
0167   using TrackStateProxy = typename track_proxy_t::ConstTrackStateProxy;
0168
0169   auto intersect = [&](const TrackStateProxy &state) -> SurfaceIntersection {
0170     assert(state.hasSmoothed() || state.hasFiltered());
0171
0172     FreeVector freeVector;
0173     if (state.hasSmoothed()) {
0174       freeVector = MultiTrajectoryHelpers::freeSmoothed(geoContext, state);
0175     } else {
0176       freeVector = MultiTrajectoryHelpers::freeFiltered(geoContext, state);
0177     }
0178
0179     return referenceSurface
0180         .intersect(geoContext, freeVector.template segment<3>(eFreePos0),
0181                    freeVector.template segment<3>(eFreeDir0),
0182                    BoundaryCheck(true), s_onSurfaceTolerance)
0183         .closest();
0184   };
0185
0186   switch (strategy) {
0187     case TrackExtrapolationStrategy::first: {
0188       ACTS_VERBOSE("looking for first track state");
0189
0190       auto first = findFirstMeasurementState(track);
0191       if (!first.ok()) {
0192         ACTS_ERROR("no first track state found");
0193         return first.error();
0194       }
0195
0196       SurfaceIntersection intersection = intersect(*first);
0197       if (!intersection) {
0198         ACTS_ERROR("no intersection found");
0199         return Result<std::pair<TrackStateProxy, double>>::failure(
0200             TrackExtrapolationError::ReferenceSurfaceUnreachable);
0201       }
0202
0203       ACTS_VERBOSE("found intersection at " << intersection.pathLength());
0204       return std::make_pair(*first, intersection.pathLength());
0205     }
0206
0207     case TrackExtrapolationStrategy::last: {
0208       ACTS_VERBOSE("looking for last track state");
0209
0210       auto last = findLastMeasurementState(track);
0211       if (!last.ok()) {
0212         ACTS_ERROR("no last track state found");
0213         return last.error();
0214       }
0215
0216       SurfaceIntersection intersection = intersect(*last);
0217       if (!intersection) {
0218         ACTS_ERROR("no intersection found");
0219         return Result<std::pair<TrackStateProxy, double>>::failure(
0220             TrackExtrapolationError::ReferenceSurfaceUnreachable);
0221       }
0222
0223       ACTS_VERBOSE("found intersection at " << intersection.pathLength());
0224       return std::make_pair(*last, intersection.pathLength());
0225     }
0226
0227     case TrackExtrapolationStrategy::firstOrLast: {
0228       ACTS_VERBOSE("looking for first or last track state");
0229
0230       auto first = findFirstMeasurementState(track);
0231       if (!first.ok()) {
0232         ACTS_ERROR("no first track state found");
0233         return first.error();
0234       }
0235
0236       auto last = findLastMeasurementState(track);
0237       if (!last.ok()) {
0238         ACTS_ERROR("no last track state found");
0239         return last.error();
0240       }
0241
0242       SurfaceIntersection intersectionFirst = intersect(*first);
0243       SurfaceIntersection intersectionLast = intersect(*last);
0244
0245       double absDistanceFirst = std::abs(intersectionFirst.pathLength());
0246       double absDistanceLast = std::abs(intersectionLast.pathLength());
0247
0248       if (intersectionFirst && absDistanceFirst <= absDistanceLast) {
0249         ACTS_VERBOSE("using first track state with intersection at "
0250                      << intersectionFirst.pathLength());
0251         return std::make_pair(*first, intersectionFirst.pathLength());
0252       }
0253
0254       if (intersectionLast && absDistanceLast <= absDistanceFirst) {
0255         ACTS_VERBOSE("using last track state with intersection at "
0256                      << intersectionLast.pathLength());
0257         return std::make_pair(*last, intersectionLast.pathLength());
0258       }
0259
0260       ACTS_ERROR("no intersection found");
0261       return Result<std::pair<TrackStateProxy, double>>::failure(
0262           TrackExtrapolationError::ReferenceSurfaceUnreachable);
0263     }
0264   }
0265
0266   // unreachable
0267   return Result<std::pair<TrackStateProxy, double>>::failure(
0268       TrackExtrapolationError::CompatibleTrackStateNotFound);
0269 }
0270
0271 /// @brief Extrapolate a track to a reference surface
0272 ///
0273 /// @tparam track_proxy_t The track proxy type
0274 /// @tparam propagator_t The propagator type
0275 /// @tparam propagator_options_t The propagator options type
0276 ///
0277 /// @param track The track which is modified in-place
0278 /// @param referenceSurface The reference surface
0279 /// @param propagator The propagator
0280 /// @param options The propagator options
0281 /// @param strategy The extrapolation strategy
0282 /// @param logger The logger
0283 ///
0284 /// @return The result of the extrapolation
0285 template <typename track_proxy_t, typename propagator_t,
0286           typename propagator_options_t>
0287 Result<void> extrapolateTrackToReferenceSurface(
0288     track_proxy_t &track, const Surface &referenceSurface,
0289     const propagator_t &propagator, propagator_options_t options,
0290     TrackExtrapolationStrategy strategy,
0291     const Logger &logger = *getDefaultLogger("TrackExtrapolation",
0292                                              Logging::INFO)) {
0293   auto findResult = findTrackStateForExtrapolation(
0294       options.geoContext, track, referenceSurface, strategy, logger);
0295
0296   if (!findResult.ok()) {
0297     ACTS_ERROR("failed to find track state for extrapolation");
0298     return findResult.error();
0299   }
0300
0301   auto &[trackState, distance] = *findResult;
0302
0303   options.direction = Direction::fromScalarZeroAsPositive(distance);
0304
0305   BoundTrackParameters parameters = track.createParametersFromState(trackState);
0306   ACTS_VERBOSE("extrapolating track to reference surface at distance "
0307                << distance << " with direction " << options.direction
0308                << " with starting parameters " << parameters);
0309
0310   auto propagateResult =
0311       propagator.template propagate<BoundTrackParameters, propagator_options_t,
0312                                     ForcedSurfaceReached>(
0313           parameters, referenceSurface, options);
0314
0315   if (!propagateResult.ok()) {
0316     ACTS_ERROR("failed to extrapolate track: " << propagateResult.error());
0317     return propagateResult.error();
0318   }
0319
0320   track.setReferenceSurface(referenceSurface.getSharedPtr());
0321   track.parameters() = propagateResult->endParameters.value().parameters();
0322   track.covariance() =
0323       propagateResult->endParameters.value().covariance().value();
0324
0325   return Result<void>::success();
0326 }
0327
0328 /// @brief Extrapolate tracks to a reference surface
0329 ///
0330 /// @tparam track_container_t The track container type
0331 /// @tparam propagator_t The propagator type
0332 /// @tparam propagator_options_t The propagator options type
0333 ///
0334 /// @param trackContainer The track container which is modified in-place
0335 /// @param referenceSurface The reference surface
0336 /// @param propagator The propagator
0337 /// @param options The propagator options
0338 /// @param strategy The extrapolation strategy
0339 /// @param logger The logger
0340 ///
0341 /// @return The result of the extrapolation
0342 template <typename track_container_t, typename propagator_t,
0343           typename propagator_options_t>
0344 Result<void> extrapolateTracksToReferenceSurface(
0345     const track_container_t &trackContainer, const Surface &referenceSurface,
0346     const propagator_t &propagator, propagator_options_t options,
0347     TrackExtrapolationStrategy strategy,
0348     const Logger &logger = *getDefaultLogger("TrackExtrapolation",
0349                                              Logging::INFO)) {
0350   Result<void> result = Result<void>::success();
0351
0352   for (const auto &track : trackContainer) {
0353     auto extrapolateResult = extrapolateTrackToReferenceSurface(
0354         track, referenceSurface, propagator, options, strategy, logger);
0355
0356     // Only keep the first error
0357     if (!extrapolateResult.ok() && result.ok()) {
0358       result = extrapolateResult.error();
0359     }
0360   }
0361
0362   return result;
0363 }
0364
0365 }  // namespace Acts
0366
0367 namespace std {
0368 // register with STL
0369 template <>
0370 struct is_error_code_enum<Acts::TrackExtrapolationError> : std::true_type {};
0371 }  // namespace std