Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2026-11-15 09:22:41

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 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 https://mozilla.org/MPL/2.0/.
0008 
0009 #pragma once
0010 
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Definitions/Direction.hpp"
0013 #include "Acts/Definitions/TrackParametrization.hpp"
0014 #include "Acts/EventData/AnyTrackStateProxy.hpp"
0015 #include "Acts/EventData/BoundTrackParameters.hpp"
0016 #include "Acts/EventData/MeasurementHelpers.hpp"
0017 #include "Acts/EventData/TrackContainerFrontendConcept.hpp"
0018 #include "Acts/EventData/TrackProxyConcept.hpp"
0019 #include "Acts/EventData/TrackStateProxyConcept.hpp"
0020 #include "Acts/EventData/TrackStateType.hpp"
0021 #include "Acts/EventData/TransformationHelpers.hpp"
0022 #include "Acts/Geometry/GeometryContext.hpp"
0023 #include "Acts/Propagator/StandardAborters.hpp"
0024 #include "Acts/Surfaces/BoundaryTolerance.hpp"
0025 #include "Acts/Surfaces/Surface.hpp"
0026 #include "Acts/TrackFitting/GainMatrixSmoother.hpp"
0027 #include "Acts/Utilities/Logger.hpp"
0028 #include "Acts/Utilities/Result.hpp"
0029 
0030 #include <optional>
0031 #include <utility>
0032 
0033 namespace Acts {
0034 
0035 /// Strategy for track extrapolation to target surface
0036 enum class TrackExtrapolationStrategy {
0037   /// Use the first track state to reach target surface
0038   first,
0039   /// Use the last track state to reach target surface
0040   last,
0041   /// Use the first or last track state to reach target surface depending on the
0042   /// distance
0043   firstOrLast,
0044 };
0045 
0046 /// Error codes for track extrapolation operations
0047 /// @ingroup errors
0048 enum class TrackExtrapolationError {
0049   /// Did not find a compatible track state
0050   CompatibleTrackStateNotFound = 1,
0051   /// Provided reference surface is unreachable
0052   ReferenceSurfaceUnreachable = 2,
0053 };
0054 
0055 /// Create error code from TrackExtrapolationError
0056 /// @param e The error code enum value
0057 /// @return Standard error code
0058 std::error_code make_error_code(TrackExtrapolationError e);
0059 
0060 /// Find the first measurement state in a track
0061 /// @param track The track to search
0062 /// @return Result containing the first measurement state proxy or error
0063 template <TrackProxyConcept track_proxy_t>
0064 Result<typename track_proxy_t::ConstTrackStateProxy> findFirstMeasurementState(
0065     const track_proxy_t &track) {
0066   using TrackStateProxy = typename track_proxy_t::ConstTrackStateProxy;
0067 
0068   // TODO specialize if track is forward linked
0069 
0070   std::optional<TrackStateProxy> firstMeasurementOpt;
0071 
0072   for (const auto &trackState : track.trackStatesReversed()) {
0073     if (trackState.typeFlags().isMeasurement()) {
0074       firstMeasurementOpt = trackState;
0075     }
0076   }
0077 
0078   if (firstMeasurementOpt) {
0079     return Result<TrackStateProxy>::success(*firstMeasurementOpt);
0080   }
0081 
0082   return Result<TrackStateProxy>::failure(
0083       TrackExtrapolationError::CompatibleTrackStateNotFound);
0084 }
0085 
0086 /// Find the last measurement state in a track
0087 /// @param track The track to search
0088 /// @return Result containing the last measurement state proxy or error
0089 template <TrackProxyConcept track_proxy_t>
0090 Result<typename track_proxy_t::ConstTrackStateProxy> findLastMeasurementState(
0091     const track_proxy_t &track) {
0092   using TrackStateProxy = typename track_proxy_t::ConstTrackStateProxy;
0093 
0094   for (const auto &trackState : track.trackStatesReversed()) {
0095     if (trackState.typeFlags().isMeasurement()) {
0096       return TrackStateProxy{trackState};
0097     }
0098   }
0099 
0100   return Result<TrackStateProxy>::failure(
0101       TrackExtrapolationError::CompatibleTrackStateNotFound);
0102 }
0103 
0104 /// @brief Smooth a track using the gain matrix smoother
0105 ///
0106 /// @tparam track_proxy_t The track proxy type
0107 /// @tparam smoother_t The smoother type
0108 ///
0109 /// @param geoContext The geometry context
0110 /// @param track The track to smooth
0111 /// @param logger The logger
0112 /// @param smoother The smoother
0113 ///
0114 /// @return The result of the smoothing
0115 template <TrackProxyConcept track_proxy_t,
0116           typename smoother_t = GainMatrixSmoother>
0117 Result<void> smoothTrack(
0118     const GeometryContext &geoContext, track_proxy_t &track,
0119     const Logger &logger = *getDefaultLogger("TrackSmoother", Logging::INFO),
0120     smoother_t smoother = GainMatrixSmoother()) {
0121   auto &trackContainer = track.container();
0122   auto &trackStateContainer = trackContainer.trackStateContainer();
0123 
0124   auto last = findLastMeasurementState(track);
0125   if (!last.ok()) {
0126     ACTS_DEBUG("no last track state found");
0127     return last.error();
0128   }
0129 
0130   auto smoothingResult =
0131       smoother(geoContext, trackStateContainer, last->index(), logger);
0132 
0133   if (!smoothingResult.ok()) {
0134     ACTS_DEBUG("Smoothing track " << track.index() << " failed with error "
0135                                   << smoothingResult.error());
0136     return smoothingResult.error();
0137   }
0138 
0139   return Result<void>::success();
0140 }
0141 
0142 /// @brief Smooth tracks using the gain matrix smoother
0143 ///
0144 /// @tparam track_container_t The track container type
0145 ///
0146 /// @param geoContext The geometry context
0147 /// @param trackContainer The track container
0148 /// @param logger The logger
0149 ///
0150 /// @return The result of the smoothing
0151 template <TrackContainerFrontend track_container_t>
0152 Result<void> smoothTracks(
0153     const GeometryContext &geoContext, const track_container_t &trackContainer,
0154     const Logger &logger = *getDefaultLogger("TrackSmoother", Logging::INFO)) {
0155   Result<void> result = Result<void>::success();
0156 
0157   for (const auto &track : trackContainer) {
0158     auto smoothingResult = smoothTrack(geoContext, track, logger);
0159 
0160     // Only keep the first error
0161     if (!smoothingResult.ok() && result.ok()) {
0162       result = smoothingResult.error();
0163     }
0164   }
0165 
0166   return result;
0167 }
0168 
0169 /// @brief Find a track state for extrapolation
0170 ///
0171 /// @tparam track_proxy_t The track proxy type
0172 ///
0173 /// @param geoContext The geometry context
0174 /// @param track The track
0175 /// @param referenceSurface The reference surface
0176 /// @param strategy The extrapolation strategy
0177 /// @param logger The logger
0178 ///
0179 /// @return The result of the search containing the track state
0180 ///         and the distance to the reference surface
0181 template <TrackProxyConcept track_proxy_t>
0182 Result<std::pair<typename track_proxy_t::ConstTrackStateProxy, double>>
0183 findTrackStateForExtrapolation(
0184     const GeometryContext &geoContext, const track_proxy_t &track,
0185     const Surface &referenceSurface, TrackExtrapolationStrategy strategy,
0186     const Logger &logger = *getDefaultLogger("TrackExtrapolation",
0187                                              Logging::INFO)) {
0188   using TrackStateProxy = typename track_proxy_t::ConstTrackStateProxy;
0189 
0190   // Intersect the reference surface with the trajectory at a track state.
0191   // Returns `std::nullopt` if the state carries no parameters at all and can
0192   // therefore not be started from; an invalid intersection means the state has
0193   // parameters but does not reach the reference surface.
0194   auto intersect =
0195       [&](const TrackStateProxy &state) -> std::optional<Intersection3D> {
0196     if (!state.hasSmoothed() && !state.hasFiltered() && !state.hasPredicted()) {
0197       return std::nullopt;
0198     }
0199 
0200     // `parameters` picks smoothed over filtered over predicted, which is what
0201     // `TrackProxy::createParametersFromState` starts the propagation from. The
0202     // distance has to be measured on the same parameters.
0203     const FreeVector freeVector = transformBoundToFreeParameters(
0204         state.referenceSurface(), geoContext, state.parameters());
0205 
0206     return referenceSurface
0207         .intersect(geoContext, freeVector.template segment<3>(eFreePos0),
0208                    freeVector.template segment<3>(eFreeDir0),
0209                    BoundaryTolerance::None(), s_onSurfaceTolerance)
0210         .closest();
0211   };
0212 
0213   switch (strategy) {
0214     case TrackExtrapolationStrategy::first: {
0215       ACTS_VERBOSE("looking for first track state");
0216 
0217       auto first = findFirstMeasurementState(track);
0218       if (!first.ok()) {
0219         ACTS_DEBUG("no first track state found");
0220         return first.error();
0221       }
0222 
0223       std::optional<Intersection3D> intersection = intersect(*first);
0224       if (!intersection.has_value()) {
0225         ACTS_DEBUG("first track state carries no parameters");
0226         return Result<std::pair<TrackStateProxy, double>>::failure(
0227             TrackExtrapolationError::CompatibleTrackStateNotFound);
0228       }
0229       if (!intersection->isValid()) {
0230         ACTS_DEBUG("no intersection found");
0231         return Result<std::pair<TrackStateProxy, double>>::failure(
0232             TrackExtrapolationError::ReferenceSurfaceUnreachable);
0233       }
0234 
0235       ACTS_VERBOSE("found intersection at " << intersection->pathLength());
0236       return std::pair(*first, intersection->pathLength());
0237     }
0238 
0239     case TrackExtrapolationStrategy::last: {
0240       ACTS_VERBOSE("looking for last track state");
0241 
0242       auto last = findLastMeasurementState(track);
0243       if (!last.ok()) {
0244         ACTS_DEBUG("no last track state found");
0245         return last.error();
0246       }
0247 
0248       std::optional<Intersection3D> intersection = intersect(*last);
0249       if (!intersection.has_value()) {
0250         ACTS_DEBUG("last track state carries no parameters");
0251         return Result<std::pair<TrackStateProxy, double>>::failure(
0252             TrackExtrapolationError::CompatibleTrackStateNotFound);
0253       }
0254       if (!intersection->isValid()) {
0255         ACTS_DEBUG("no intersection found");
0256         return Result<std::pair<TrackStateProxy, double>>::failure(
0257             TrackExtrapolationError::ReferenceSurfaceUnreachable);
0258       }
0259 
0260       ACTS_VERBOSE("found intersection at " << intersection->pathLength());
0261       return std::pair(*last, intersection->pathLength());
0262     }
0263 
0264     case TrackExtrapolationStrategy::firstOrLast: {
0265       ACTS_VERBOSE("looking for first or last track state");
0266 
0267       auto first = findFirstMeasurementState(track);
0268       if (!first.ok()) {
0269         ACTS_DEBUG("no first track state found");
0270         return first.error();
0271       }
0272 
0273       auto last = findLastMeasurementState(track);
0274       if (!last.ok()) {
0275         ACTS_DEBUG("no last track state found");
0276         return last.error();
0277       }
0278 
0279       std::optional<Intersection3D> intersectionFirstOpt = intersect(*first);
0280       std::optional<Intersection3D> intersectionLastOpt = intersect(*last);
0281 
0282       if (!intersectionFirstOpt.has_value() &&
0283           !intersectionLastOpt.has_value()) {
0284         ACTS_DEBUG("neither first nor last track state carries parameters");
0285         return Result<std::pair<TrackStateProxy, double>>::failure(
0286             TrackExtrapolationError::CompatibleTrackStateNotFound);
0287       }
0288 
0289       // an end without parameters cannot be started from, so it loses the
0290       // comparison below through the infinite path length of an invalid
0291       // intersection
0292       Intersection3D intersectionFirst =
0293           intersectionFirstOpt.value_or(Intersection3D::Invalid());
0294       Intersection3D intersectionLast =
0295           intersectionLastOpt.value_or(Intersection3D::Invalid());
0296 
0297       double absDistanceFirst = std::abs(intersectionFirst.pathLength());
0298       double absDistanceLast = std::abs(intersectionLast.pathLength());
0299 
0300       if (intersectionFirst.isValid() && absDistanceFirst <= absDistanceLast) {
0301         ACTS_VERBOSE("using first track state with intersection at "
0302                      << intersectionFirst.pathLength());
0303         return std::pair(*first, intersectionFirst.pathLength());
0304       }
0305 
0306       if (intersectionLast.isValid() && absDistanceLast <= absDistanceFirst) {
0307         ACTS_VERBOSE("using last track state with intersection at "
0308                      << intersectionLast.pathLength());
0309         return std::pair(*last, intersectionLast.pathLength());
0310       }
0311 
0312       ACTS_DEBUG("no intersection found");
0313       return Result<std::pair<TrackStateProxy, double>>::failure(
0314           TrackExtrapolationError::ReferenceSurfaceUnreachable);
0315     }
0316   }
0317 
0318   // unreachable
0319   return Result<std::pair<TrackStateProxy, double>>::failure(
0320       TrackExtrapolationError::CompatibleTrackStateNotFound);
0321 }
0322 
0323 /// @brief Extrapolate a track to a reference surface
0324 ///
0325 /// @tparam track_proxy_t The track proxy type
0326 /// @tparam propagator_t The propagator type
0327 /// @tparam propagator_options_t The propagator options type
0328 ///
0329 /// @param track The track which is modified in-place
0330 /// @param referenceSurface The reference surface
0331 /// @param propagator The propagator
0332 /// @param options The propagator options
0333 /// @param strategy The extrapolation strategy
0334 /// @param logger The logger
0335 ///
0336 /// @return The result of the extrapolation
0337 template <TrackProxyConcept track_proxy_t, typename propagator_t,
0338           typename propagator_options_t>
0339 Result<void> extrapolateTrackToReferenceSurface(
0340     track_proxy_t &track, const Surface &referenceSurface,
0341     const propagator_t &propagator, propagator_options_t options,
0342     TrackExtrapolationStrategy strategy,
0343     const Logger &logger = *getDefaultLogger("TrackExtrapolation",
0344                                              Logging::INFO)) {
0345   auto findResult = findTrackStateForExtrapolation(
0346       options.geoContext, track, referenceSurface, strategy, logger);
0347 
0348   if (!findResult.ok()) {
0349     ACTS_DEBUG("failed to find track state for extrapolation");
0350     return findResult.error();
0351   }
0352 
0353   auto &[trackState, distance] = *findResult;
0354 
0355   options.direction = Direction::fromScalarZeroAsPositive(distance);
0356 
0357   BoundTrackParameters parameters = track.createParametersFromState(trackState);
0358   ACTS_VERBOSE("extrapolating track to reference surface at distance "
0359                << distance << " with direction " << options.direction
0360                << " with starting parameters " << parameters);
0361 
0362   auto propagateResult =
0363       propagator.template propagate<propagator_options_t, ForcedSurfaceReached>(
0364           parameters, referenceSurface, options);
0365 
0366   if (!propagateResult.ok()) {
0367     ACTS_DEBUG("failed to extrapolate track: " << propagateResult.error());
0368     return propagateResult.error();
0369   }
0370 
0371   track.setReferenceSurface(referenceSurface.getSharedPtr());
0372   track.parameters() = propagateResult->endParameters.value().parameters();
0373   track.covariance() =
0374       propagateResult->endParameters.value().covariance().value();
0375 
0376   return Result<void>::success();
0377 }
0378 
0379 /// @brief Extrapolate tracks to a reference surface
0380 ///
0381 /// @tparam track_container_t The track container type
0382 /// @tparam propagator_t The propagator type
0383 /// @tparam propagator_options_t The propagator options type
0384 ///
0385 /// @param trackContainer The track container which is modified in-place
0386 /// @param referenceSurface The reference surface
0387 /// @param propagator The propagator
0388 /// @param options The propagator options
0389 /// @param strategy The extrapolation strategy
0390 /// @param logger The logger
0391 ///
0392 /// @return The result of the extrapolation
0393 template <TrackContainerFrontend track_container_t, typename propagator_t,
0394           typename propagator_options_t>
0395 Result<void> extrapolateTracksToReferenceSurface(
0396     const track_container_t &trackContainer, const Surface &referenceSurface,
0397     const propagator_t &propagator, propagator_options_t options,
0398     TrackExtrapolationStrategy strategy,
0399     const Logger &logger = *getDefaultLogger("TrackExtrapolation",
0400                                              Logging::INFO)) {
0401   Result<void> result = Result<void>::success();
0402 
0403   for (const auto &track : trackContainer) {
0404     auto extrapolateResult = extrapolateTrackToReferenceSurface(
0405         track, referenceSurface, propagator, options, strategy, logger);
0406 
0407     // Only keep the first error
0408     if (!extrapolateResult.ok() && result.ok()) {
0409       result = extrapolateResult.error();
0410     }
0411   }
0412 
0413   return result;
0414 }
0415 
0416 /// Helper function to calculate a number of track level quantities and store
0417 /// them on the track itself
0418 /// @tparam track_proxy_t The track proxy type
0419 /// @param track A mutable track proxy to operate on
0420 template <TrackProxyConcept track_proxy_t>
0421 void calculateTrackQuantities(track_proxy_t track)
0422   requires(!track_proxy_t::ReadOnly)
0423 {
0424   track.chi2() = 0;
0425   track.nDoF() = 0;
0426 
0427   track.nHoles() = 0;
0428   track.nMeasurements() = 0;
0429   track.nSharedHits() = 0;
0430   track.nSplitHits() = 0;
0431   track.nOutliers() = 0;
0432 
0433   for (const auto &trackState : track.trackStatesReversed()) {
0434     ConstTrackStateTypeMap typeFlags = trackState.typeFlags();
0435 
0436     if (typeFlags.isHole()) {
0437       track.nHoles()++;
0438     } else if (typeFlags.isOutlier()) {
0439       track.nOutliers()++;
0440     } else if (typeFlags.isMeasurement()) {
0441       if (typeFlags.isSharedHit()) {
0442         track.nSharedHits()++;
0443       }
0444       if (typeFlags.isSplitHit()) {
0445         track.nSplitHits()++;
0446       }
0447       track.nMeasurements()++;
0448       track.chi2() += trackState.chi2();
0449       track.nDoF() += trackState.calibratedSize();
0450     }
0451   }
0452 }
0453 
0454 /// Helper function to trim track states from the front of a track
0455 /// @tparam track_proxy_t the track proxy type
0456 /// @param track the track to trim
0457 /// @param trimHoles whether to trim holes
0458 /// @param trimOutliers whether to trim outliers
0459 /// @param trimMaterial whether to trim pure material states
0460 /// @param trimOtherNoneMeasurement whether to trim other, non measurement, states
0461 template <TrackProxyConcept track_proxy_t>
0462 void trimTrackFront(track_proxy_t track, bool trimHoles, bool trimOutliers,
0463                     bool trimMaterial, bool trimOtherNoneMeasurement)
0464   requires(!track_proxy_t::ReadOnly)
0465 {
0466   using TrackStateProxy = typename track_proxy_t::TrackStateProxy;
0467 
0468   // TODO specialize if track is forward linked
0469 
0470   std::optional<TrackStateProxy> front;
0471 
0472   for (TrackStateProxy trackState : track.trackStatesReversed()) {
0473     TrackStateTypeMap typeFlags = trackState.typeFlags();
0474     bool isHole = typeFlags.isHole();
0475     bool isOutlier = typeFlags.isOutlier();
0476     bool isMaterial = typeFlags.isMaterial();
0477     bool isOtherNoneMeasurement =
0478         !typeFlags.hasMeasurement() && !isHole && !isOutlier && !isMaterial;
0479     if (trimHoles && isHole) {
0480       continue;
0481     }
0482     if (trimOutliers && isOutlier) {
0483       continue;
0484     }
0485     if (trimMaterial && isMaterial) {
0486       continue;
0487     }
0488     if (trimOtherNoneMeasurement && isOtherNoneMeasurement) {
0489       continue;
0490     }
0491 
0492     front = trackState;
0493   }
0494 
0495   if (front.has_value()) {
0496     front.value().previous() = TrackStateProxy::kInvalid;
0497   }
0498 }
0499 
0500 /// Helper function to trim track states from the back of a track
0501 /// @tparam track_proxy_t the track proxy type
0502 /// @param track the track to trim
0503 /// @param trimHoles whether to trim holes
0504 /// @param trimOutliers whether to trim outliers
0505 /// @param trimMaterial whether to trim pure material states
0506 /// @param trimOtherNoneMeasurement whether to trim other, non measurement, states
0507 template <TrackProxyConcept track_proxy_t>
0508 void trimTrackBack(track_proxy_t track, bool trimHoles, bool trimOutliers,
0509                    bool trimMaterial, bool trimOtherNoneMeasurement)
0510   requires(!track_proxy_t::ReadOnly)
0511 {
0512   using TrackStateProxy = typename track_proxy_t::TrackStateProxy;
0513 
0514   std::optional<TrackStateProxy> back;
0515 
0516   for (TrackStateProxy trackState : track.trackStatesReversed()) {
0517     back = trackState;
0518 
0519     TrackStateTypeMap typeFlags = trackState.typeFlags();
0520     bool isHole = typeFlags.isHole();
0521     bool isOutlier = typeFlags.isOutlier();
0522     bool isMaterial = typeFlags.isMaterial();
0523     bool isOtherNoneMeasurement =
0524         !typeFlags.hasMeasurement() && !isHole && !isOutlier && !isMaterial;
0525     if (trimHoles && isHole) {
0526       continue;
0527     }
0528     if (trimOutliers && isOutlier) {
0529       continue;
0530     }
0531     if (trimMaterial && isMaterial) {
0532       continue;
0533     }
0534     if (trimOtherNoneMeasurement && isOtherNoneMeasurement) {
0535       continue;
0536     }
0537 
0538     break;
0539   }
0540 
0541   if (back.has_value()) {
0542     track.tipIndex() = back.value().index();
0543   }
0544 }
0545 
0546 /// Helper function to trim track states from the front and back of a track
0547 /// @tparam track_proxy_t the track proxy type
0548 /// @param track the track to trim
0549 /// @param trimHoles whether to trim holes
0550 /// @param trimOutliers whether to trim outliers
0551 /// @param trimMaterial whether to trim pure material states
0552 /// @param trimOtherNoneMeasurement whether to trim other, non measurement, states
0553 template <TrackProxyConcept track_proxy_t>
0554 void trimTrack(track_proxy_t track, bool trimHoles, bool trimOutliers,
0555                bool trimMaterial, bool trimOtherNoneMeasurement)
0556   requires(!track_proxy_t::ReadOnly)
0557 {
0558   trimTrackFront(track, trimHoles, trimOutliers, trimMaterial,
0559                  trimOtherNoneMeasurement);
0560   trimTrackBack(track, trimHoles, trimOutliers, trimMaterial,
0561                 trimOtherNoneMeasurement);
0562 }
0563 
0564 /// Helper function to calculate the predicted residual and its covariance
0565 /// @tparam nMeasurementDim the dimension of the measurement
0566 /// @tparam track_state_proxy_t the track state proxy type
0567 /// @param trackState the track state to calculate the residual from
0568 /// @return a pair of the residual and its covariance
0569 template <std::size_t nMeasurementDim,
0570           TrackStateProxyConcept track_state_proxy_t>
0571 std::pair<Vector<nMeasurementDim>, SquareMatrix<nMeasurementDim>>
0572 calculatePredictedResidual(track_state_proxy_t trackState) {
0573   using MeasurementVector = Vector<nMeasurementDim>;
0574   using MeasurementMatrix = SquareMatrix<nMeasurementDim>;
0575 
0576   if (!trackState.hasPredicted()) {
0577     throw std::invalid_argument("track state has no predicted parameters");
0578   }
0579   if (!trackState.hasCalibrated()) {
0580     throw std::invalid_argument("track state has no calibrated parameters");
0581   }
0582 
0583   auto subspaceHelper =
0584       trackState.template projectorSubspaceHelper<nMeasurementDim>();
0585 
0586   auto measurement = trackState.template calibrated<nMeasurementDim>();
0587   auto measurementCovariance =
0588       trackState.template calibratedCovariance<nMeasurementDim>();
0589   MeasurementVector predicted =
0590       subspaceHelper.projectVector(trackState.predicted());
0591   MeasurementMatrix predictedCovariance =
0592       subspaceHelper.projectMatrix(trackState.predictedCovariance());
0593 
0594   MeasurementVector residual = measurement - predicted;
0595   MeasurementMatrix residualCovariance =
0596       measurementCovariance + predictedCovariance;
0597 
0598   return {residual, residualCovariance};
0599 }
0600 
0601 /// Helper function to calculate the filtered residual and its covariance
0602 /// @tparam nMeasurementDim the dimension of the measurement
0603 /// @tparam track_state_proxy_t the track state proxy type
0604 /// @param trackState the track state to calculate the residual from
0605 /// @return a pair of the residual and its covariance
0606 template <std::size_t nMeasurementDim,
0607           TrackStateProxyConcept track_state_proxy_t>
0608 std::pair<Vector<nMeasurementDim>, SquareMatrix<nMeasurementDim>>
0609 calculateFilteredResidual(track_state_proxy_t trackState) {
0610   using MeasurementVector = Vector<nMeasurementDim>;
0611   using MeasurementMatrix = SquareMatrix<nMeasurementDim>;
0612 
0613   if (!trackState.hasFiltered()) {
0614     throw std::invalid_argument("track state has no filtered parameters");
0615   }
0616   if (!trackState.hasCalibrated()) {
0617     throw std::invalid_argument("track state has no calibrated parameters");
0618   }
0619 
0620   auto subspaceHelper =
0621       trackState.template projectorSubspaceHelper<nMeasurementDim>();
0622 
0623   auto measurement = trackState.template calibrated<nMeasurementDim>();
0624   auto measurementCovariance =
0625       trackState.template calibratedCovariance<nMeasurementDim>();
0626   MeasurementVector filtered =
0627       subspaceHelper.projectVector(trackState.filtered());
0628   MeasurementMatrix filteredCovariance =
0629       subspaceHelper.projectMatrix(trackState.filteredCovariance());
0630 
0631   MeasurementVector residual = measurement - filtered;
0632   MeasurementMatrix residualCovariance =
0633       measurementCovariance - filteredCovariance;
0634 
0635   return {residual, residualCovariance};
0636 }
0637 
0638 /// Helper function to calculate the smoothed residual and its covariance
0639 /// @tparam nMeasurementDim the dimension of the measurement
0640 /// @tparam track_state_proxy_t the track state proxy type
0641 /// @param trackState the track state to calculate the residual from
0642 /// @return a pair of the residual and its covariance
0643 template <std::size_t nMeasurementDim,
0644           TrackStateProxyConcept track_state_proxy_t>
0645 std::pair<Vector<nMeasurementDim>, SquareMatrix<nMeasurementDim>>
0646 calculateSmoothedResidual(track_state_proxy_t trackState) {
0647   using MeasurementVector = Vector<nMeasurementDim>;
0648   using MeasurementMatrix = SquareMatrix<nMeasurementDim>;
0649 
0650   if (!trackState.hasSmoothed()) {
0651     throw std::invalid_argument("track state has no smoothed parameters");
0652   }
0653   if (!trackState.hasCalibrated()) {
0654     throw std::invalid_argument("track state has no calibrated parameters");
0655   }
0656 
0657   auto subspaceHelper =
0658       trackState.template projectorSubspaceHelper<nMeasurementDim>();
0659 
0660   auto measurement = trackState.template calibrated<nMeasurementDim>();
0661   auto measurementCovariance =
0662       trackState.template calibratedCovariance<nMeasurementDim>();
0663   MeasurementVector smoothed =
0664       subspaceHelper.projectVector(trackState.smoothed());
0665   MeasurementMatrix smoothedCovariance =
0666       subspaceHelper.projectMatrix(trackState.smoothedCovariance());
0667 
0668   MeasurementVector residual = measurement - smoothed;
0669   MeasurementMatrix residualCovariance =
0670       measurementCovariance - smoothedCovariance;
0671 
0672   return {residual, residualCovariance};
0673 }
0674 
0675 /// Helper function to calculate the predicted chi2
0676 /// @tparam track_state_proxy_t the track state proxy type
0677 /// @param trackState the track state to calculate the chi2 from
0678 /// @return the chi2
0679 template <TrackStateProxyConcept track_state_proxy_t>
0680 double calculatePredictedChi2(track_state_proxy_t trackState) {
0681   if (!trackState.hasPredicted()) {
0682     throw std::invalid_argument("track state has no predicted parameters");
0683   }
0684   if (!trackState.hasCalibrated()) {
0685     throw std::invalid_argument("track state has no calibrated parameters");
0686   }
0687 
0688   return visit_measurement(
0689       trackState.calibratedSize(),
0690       [&]<std::size_t measdim>(
0691           std::integral_constant<std::size_t, measdim>) -> double {
0692         auto [residual, residualCovariance] =
0693             calculatePredictedResidual<measdim>(trackState);
0694 
0695         return (residual.transpose() * residualCovariance.inverse() * residual)
0696             .eval()(0, 0);
0697       });
0698 }
0699 
0700 /// Helper function to calculate the filtered chi2
0701 /// @tparam track_state_proxy_t the track state proxy type
0702 /// @param trackState the track state to calculate the chi2 from
0703 /// @return the chi2
0704 template <TrackStateProxyConcept track_state_proxy_t>
0705 double calculateFilteredChi2(track_state_proxy_t trackState) {
0706   if (!trackState.hasFiltered()) {
0707     throw std::invalid_argument("track state has no filtered parameters");
0708   }
0709   if (!trackState.hasCalibrated()) {
0710     throw std::invalid_argument("track state has no calibrated parameters");
0711   }
0712 
0713   return visit_measurement(
0714       trackState.calibratedSize(),
0715       [&]<std::size_t measdim>(
0716           std::integral_constant<std::size_t, measdim>) -> double {
0717         auto [residual, residualCovariance] =
0718             calculateFilteredResidual<measdim>(trackState);
0719 
0720         return (residual.transpose() * residualCovariance.inverse() * residual)
0721             .eval()(0, 0);
0722       });
0723 }
0724 
0725 /// Helper function to calculate the smoothed chi2
0726 /// @tparam track_state_proxy_t the track state proxy type
0727 /// @param trackState the track state to calculate the chi2 from
0728 /// @return the chi2
0729 template <TrackStateProxyConcept track_state_proxy_t>
0730 double calculateSmoothedChi2(track_state_proxy_t trackState) {
0731   if (!trackState.hasSmoothed()) {
0732     throw std::invalid_argument("track state has no smoothed parameters");
0733   }
0734   if (!trackState.hasCalibrated()) {
0735     throw std::invalid_argument("track state has no calibrated parameters");
0736   }
0737 
0738   return visit_measurement(
0739       trackState.calibratedSize(),
0740       [&]<std::size_t measdim>(
0741           std::integral_constant<std::size_t, measdim>) -> double {
0742         auto [residual, residualCovariance] =
0743             calculateSmoothedResidual<measdim>(trackState);
0744 
0745         return (residual.transpose() * residualCovariance.inverse() * residual)
0746             .eval()(0, 0);
0747       });
0748 }
0749 
0750 /// Helper function to calculate the unbiased track parameters and their
0751 /// covariance (i.e. fitted track parameters with this measurement removed)
0752 /// using Eq.(12a)-Eq.(12c) of NIMA 262, 444 (1987)
0753 /// @tparam track_state_proxy_t the track state proxy type
0754 /// @param trackState the track state to calculate the unbiased parameters from
0755 /// @return a pair of the unbiased parameters and their covariance
0756 /// @deprecated Instantiating this template is very expensive in compiler memory
0757 ///   (it expands the Eigen expression templates over all measurement
0758 ///   dimensions) and it does so in every calling translation unit. Prefer the
0759 ///   non-template overload taking a type-erased @c AnyConstTrackStateProxy,
0760 ///   which is compiled once in the Acts core library:
0761 ///   @code
0762 ///   calculateUnbiasedParametersCovariance(Acts::AnyConstTrackStateProxy{state});
0763 ///   @endcode
0764 template <TrackStateProxyConcept track_state_proxy_t>
0765 [[deprecated(
0766     "Use calculateUnbiasedParametersCovariance(const AnyConstTrackStateProxy&) "
0767     "instead; the templated form instantiates expensive Eigen code in every "
0768     "translation unit.")]]
0769 std::pair<BoundVector, BoundMatrix> calculateUnbiasedParametersCovariance(
0770     track_state_proxy_t trackState) {
0771   // Explicitly select the non-template overload taking a type-erased proxy.
0772   // A plain call here would re-resolve to this very template (the wrapped
0773   // AnyConstTrackStateProxy satisfies TrackStateProxyConcept and is an exact
0774   // by-value match), causing infinite recursion and a self-deprecation error.
0775   std::pair<BoundVector, BoundMatrix> (&impl)(const AnyConstTrackStateProxy &) =
0776       calculateUnbiasedParametersCovariance;
0777   return impl(AnyConstTrackStateProxy{trackState});
0778 }
0779 
0780 /// Calculate the unbiased track parameters and their covariance for a
0781 /// type-erased track state proxy. See the templated overload above for the
0782 /// underlying formula.
0783 ///
0784 /// This is the preferred entry point. It is not a template, so the (very
0785 /// expensive) Eigen expression templates are instantiated exactly once, in the
0786 /// Acts core library (TrackHelpersUnbiased.cpp), instead of in every calling
0787 /// translation unit. Callers holding a concrete track state proxy wrap it:
0788 /// @code
0789 /// calculateUnbiasedParametersCovariance(Acts::AnyConstTrackStateProxy{state});
0790 /// @endcode
0791 /// @param trackState the (type-erased) track state to calculate from
0792 /// @return a pair of the unbiased parameters and their covariance
0793 std::pair<BoundVector, BoundMatrix> calculateUnbiasedParametersCovariance(
0794     const AnyConstTrackStateProxy &trackState);
0795 
0796 }  // namespace Acts
0797 
0798 namespace std {
0799 // register with STL
0800 template <>
0801 struct is_error_code_enum<Acts::TrackExtrapolationError> : std::true_type {};
0802 }  // namespace std