Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2026-05-19 07:34:01

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