Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2026-07-07 08:05:29

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/MeasurementHelpers.hpp"
0015 #include "Acts/EventData/MultiTrajectoryHelpers.hpp"
0016 #include "Acts/EventData/TrackContainerFrontendConcept.hpp"
0017 #include "Acts/EventData/TrackParameters.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<BoundTrackParameters, propagator_options_t,
0331                                     ForcedSurfaceReached>(
0332           parameters, referenceSurface, options);
0333 
0334   if (!propagateResult.ok()) {
0335     ACTS_DEBUG("failed to extrapolate track: " << propagateResult.error());
0336     return propagateResult.error();
0337   }
0338 
0339   track.setReferenceSurface(referenceSurface.getSharedPtr());
0340   track.parameters() = propagateResult->endParameters.value().parameters();
0341   track.covariance() =
0342       propagateResult->endParameters.value().covariance().value();
0343 
0344   return Result<void>::success();
0345 }
0346 
0347 /// @brief Extrapolate tracks to a reference surface
0348 ///
0349 /// @tparam track_container_t The track container type
0350 /// @tparam propagator_t The propagator type
0351 /// @tparam propagator_options_t The propagator options type
0352 ///
0353 /// @param trackContainer The track container which is modified in-place
0354 /// @param referenceSurface The reference surface
0355 /// @param propagator The propagator
0356 /// @param options The propagator options
0357 /// @param strategy The extrapolation strategy
0358 /// @param logger The logger
0359 ///
0360 /// @return The result of the extrapolation
0361 template <TrackContainerFrontend track_container_t, typename propagator_t,
0362           typename propagator_options_t>
0363 Result<void> extrapolateTracksToReferenceSurface(
0364     const track_container_t &trackContainer, const Surface &referenceSurface,
0365     const propagator_t &propagator, propagator_options_t options,
0366     TrackExtrapolationStrategy strategy,
0367     const Logger &logger = *getDefaultLogger("TrackExtrapolation",
0368                                              Logging::INFO)) {
0369   Result<void> result = Result<void>::success();
0370 
0371   for (const auto &track : trackContainer) {
0372     auto extrapolateResult = extrapolateTrackToReferenceSurface(
0373         track, referenceSurface, propagator, options, strategy, logger);
0374 
0375     // Only keep the first error
0376     if (!extrapolateResult.ok() && result.ok()) {
0377       result = extrapolateResult.error();
0378     }
0379   }
0380 
0381   return result;
0382 }
0383 
0384 /// Helper function to calculate a number of track level quantities and store
0385 /// them on the track itself
0386 /// @tparam track_proxy_t The track proxy type
0387 /// @param track A mutable track proxy to operate on
0388 template <TrackProxyConcept track_proxy_t>
0389 void calculateTrackQuantities(track_proxy_t track)
0390   requires(!track_proxy_t::ReadOnly)
0391 {
0392   track.chi2() = 0;
0393   track.nDoF() = 0;
0394 
0395   track.nHoles() = 0;
0396   track.nMeasurements() = 0;
0397   track.nSharedHits() = 0;
0398   track.nOutliers() = 0;
0399 
0400   for (const auto &trackState : track.trackStatesReversed()) {
0401     ConstTrackStateTypeMap typeFlags = trackState.typeFlags();
0402 
0403     if (typeFlags.isHole()) {
0404       track.nHoles()++;
0405     } else if (typeFlags.isOutlier()) {
0406       track.nOutliers()++;
0407     } else if (typeFlags.isMeasurement()) {
0408       if (typeFlags.isSharedHit()) {
0409         track.nSharedHits()++;
0410       }
0411       track.nMeasurements()++;
0412       track.chi2() += trackState.chi2();
0413       track.nDoF() += trackState.calibratedSize();
0414     }
0415   }
0416 }
0417 
0418 /// Helper function to trim track states from the front of a track
0419 /// @tparam track_proxy_t the track proxy type
0420 /// @param track the track to trim
0421 /// @param trimHoles whether to trim holes
0422 /// @param trimOutliers whether to trim outliers
0423 /// @param trimMaterial whether to trim pure material states
0424 /// @param trimOtherNoneMeasurement whether to trim other, non measurement, states
0425 template <TrackProxyConcept track_proxy_t>
0426 void trimTrackFront(track_proxy_t track, bool trimHoles, bool trimOutliers,
0427                     bool trimMaterial, bool trimOtherNoneMeasurement)
0428   requires(!track_proxy_t::ReadOnly)
0429 {
0430   using TrackStateProxy = typename track_proxy_t::TrackStateProxy;
0431 
0432   // TODO specialize if track is forward linked
0433 
0434   std::optional<TrackStateProxy> front;
0435 
0436   for (TrackStateProxy trackState : track.trackStatesReversed()) {
0437     TrackStateTypeMap typeFlags = trackState.typeFlags();
0438     bool isHole = typeFlags.isHole();
0439     bool isOutlier = typeFlags.isOutlier();
0440     bool isMaterial = typeFlags.isMaterial();
0441     bool isOtherNoneMeasurement =
0442         !typeFlags.hasMeasurement() && !isHole && !isOutlier && !isMaterial;
0443     if (trimHoles && isHole) {
0444       continue;
0445     }
0446     if (trimOutliers && isOutlier) {
0447       continue;
0448     }
0449     if (trimMaterial && isMaterial) {
0450       continue;
0451     }
0452     if (trimOtherNoneMeasurement && isOtherNoneMeasurement) {
0453       continue;
0454     }
0455 
0456     front = trackState;
0457   }
0458 
0459   if (front.has_value()) {
0460     front.value().previous() = TrackStateProxy::kInvalid;
0461   }
0462 }
0463 
0464 /// Helper function to trim track states from the back of a track
0465 /// @tparam track_proxy_t the track proxy type
0466 /// @param track the track to trim
0467 /// @param trimHoles whether to trim holes
0468 /// @param trimOutliers whether to trim outliers
0469 /// @param trimMaterial whether to trim pure material states
0470 /// @param trimOtherNoneMeasurement whether to trim other, non measurement, states
0471 template <TrackProxyConcept track_proxy_t>
0472 void trimTrackBack(track_proxy_t track, bool trimHoles, bool trimOutliers,
0473                    bool trimMaterial, bool trimOtherNoneMeasurement)
0474   requires(!track_proxy_t::ReadOnly)
0475 {
0476   using TrackStateProxy = typename track_proxy_t::TrackStateProxy;
0477 
0478   std::optional<TrackStateProxy> back;
0479 
0480   for (TrackStateProxy trackState : track.trackStatesReversed()) {
0481     back = trackState;
0482 
0483     TrackStateTypeMap typeFlags = trackState.typeFlags();
0484     bool isHole = typeFlags.isHole();
0485     bool isOutlier = typeFlags.isOutlier();
0486     bool isMaterial = typeFlags.isMaterial();
0487     bool isOtherNoneMeasurement =
0488         !typeFlags.hasMeasurement() && !isHole && !isOutlier && !isMaterial;
0489     if (trimHoles && isHole) {
0490       continue;
0491     }
0492     if (trimOutliers && isOutlier) {
0493       continue;
0494     }
0495     if (trimMaterial && isMaterial) {
0496       continue;
0497     }
0498     if (trimOtherNoneMeasurement && isOtherNoneMeasurement) {
0499       continue;
0500     }
0501 
0502     break;
0503   }
0504 
0505   if (back.has_value()) {
0506     track.tipIndex() = back.value().index();
0507   }
0508 }
0509 
0510 /// Helper function to trim track states from the front and back of a track
0511 /// @tparam track_proxy_t the track proxy type
0512 /// @param track the track to trim
0513 /// @param trimHoles whether to trim holes
0514 /// @param trimOutliers whether to trim outliers
0515 /// @param trimMaterial whether to trim pure material states
0516 /// @param trimOtherNoneMeasurement whether to trim other, non measurement, states
0517 template <TrackProxyConcept track_proxy_t>
0518 void trimTrack(track_proxy_t track, bool trimHoles, bool trimOutliers,
0519                bool trimMaterial, bool trimOtherNoneMeasurement)
0520   requires(!track_proxy_t::ReadOnly)
0521 {
0522   trimTrackFront(track, trimHoles, trimOutliers, trimMaterial,
0523                  trimOtherNoneMeasurement);
0524   trimTrackBack(track, trimHoles, trimOutliers, trimMaterial,
0525                 trimOtherNoneMeasurement);
0526 }
0527 
0528 /// Helper function to calculate the predicted residual and its covariance
0529 /// @tparam nMeasurementDim the dimension of the measurement
0530 /// @tparam track_state_proxy_t the track state proxy type
0531 /// @param trackState the track state to calculate the residual from
0532 /// @return a pair of the residual and its covariance
0533 template <std::size_t nMeasurementDim,
0534           TrackStateProxyConcept track_state_proxy_t>
0535 std::pair<Vector<nMeasurementDim>, SquareMatrix<nMeasurementDim>>
0536 calculatePredictedResidual(track_state_proxy_t trackState) {
0537   using MeasurementVector = Vector<nMeasurementDim>;
0538   using MeasurementMatrix = SquareMatrix<nMeasurementDim>;
0539 
0540   if (!trackState.hasPredicted()) {
0541     throw std::invalid_argument("track state has no predicted parameters");
0542   }
0543   if (!trackState.hasCalibrated()) {
0544     throw std::invalid_argument("track state has no calibrated parameters");
0545   }
0546 
0547   auto subspaceHelper =
0548       trackState.template projectorSubspaceHelper<nMeasurementDim>();
0549 
0550   auto measurement = trackState.template calibrated<nMeasurementDim>();
0551   auto measurementCovariance =
0552       trackState.template calibratedCovariance<nMeasurementDim>();
0553   MeasurementVector predicted =
0554       subspaceHelper.projectVector(trackState.predicted());
0555   MeasurementMatrix predictedCovariance =
0556       subspaceHelper.projectMatrix(trackState.predictedCovariance());
0557 
0558   MeasurementVector residual = measurement - predicted;
0559   MeasurementMatrix residualCovariance =
0560       measurementCovariance + predictedCovariance;
0561 
0562   return {residual, residualCovariance};
0563 }
0564 
0565 /// Helper function to calculate the filtered residual and its covariance
0566 /// @tparam nMeasurementDim the dimension of the measurement
0567 /// @tparam track_state_proxy_t the track state proxy type
0568 /// @param trackState the track state to calculate the residual from
0569 /// @return a pair of the residual and its covariance
0570 template <std::size_t nMeasurementDim,
0571           TrackStateProxyConcept track_state_proxy_t>
0572 std::pair<Vector<nMeasurementDim>, SquareMatrix<nMeasurementDim>>
0573 calculateFilteredResidual(track_state_proxy_t trackState) {
0574   using MeasurementVector = Vector<nMeasurementDim>;
0575   using MeasurementMatrix = SquareMatrix<nMeasurementDim>;
0576 
0577   if (!trackState.hasFiltered()) {
0578     throw std::invalid_argument("track state has no filtered parameters");
0579   }
0580   if (!trackState.hasCalibrated()) {
0581     throw std::invalid_argument("track state has no calibrated parameters");
0582   }
0583 
0584   auto subspaceHelper =
0585       trackState.template projectorSubspaceHelper<nMeasurementDim>();
0586 
0587   auto measurement = trackState.template calibrated<nMeasurementDim>();
0588   auto measurementCovariance =
0589       trackState.template calibratedCovariance<nMeasurementDim>();
0590   MeasurementVector filtered =
0591       subspaceHelper.projectVector(trackState.filtered());
0592   MeasurementMatrix filteredCovariance =
0593       subspaceHelper.projectMatrix(trackState.filteredCovariance());
0594 
0595   MeasurementVector residual = measurement - filtered;
0596   MeasurementMatrix residualCovariance =
0597       measurementCovariance - filteredCovariance;
0598 
0599   return {residual, residualCovariance};
0600 }
0601 
0602 /// Helper function to calculate the smoothed residual and its covariance
0603 /// @tparam nMeasurementDim the dimension of the measurement
0604 /// @tparam track_state_proxy_t the track state proxy type
0605 /// @param trackState the track state to calculate the residual from
0606 /// @return a pair of the residual and its covariance
0607 template <std::size_t nMeasurementDim,
0608           TrackStateProxyConcept track_state_proxy_t>
0609 std::pair<Vector<nMeasurementDim>, SquareMatrix<nMeasurementDim>>
0610 calculateSmoothedResidual(track_state_proxy_t trackState) {
0611   using MeasurementVector = Vector<nMeasurementDim>;
0612   using MeasurementMatrix = SquareMatrix<nMeasurementDim>;
0613 
0614   if (!trackState.hasSmoothed()) {
0615     throw std::invalid_argument("track state has no smoothed parameters");
0616   }
0617   if (!trackState.hasCalibrated()) {
0618     throw std::invalid_argument("track state has no calibrated parameters");
0619   }
0620 
0621   auto subspaceHelper =
0622       trackState.template projectorSubspaceHelper<nMeasurementDim>();
0623 
0624   auto measurement = trackState.template calibrated<nMeasurementDim>();
0625   auto measurementCovariance =
0626       trackState.template calibratedCovariance<nMeasurementDim>();
0627   MeasurementVector smoothed =
0628       subspaceHelper.projectVector(trackState.smoothed());
0629   MeasurementMatrix smoothedCovariance =
0630       subspaceHelper.projectMatrix(trackState.smoothedCovariance());
0631 
0632   MeasurementVector residual = measurement - smoothed;
0633   MeasurementMatrix residualCovariance =
0634       measurementCovariance - smoothedCovariance;
0635 
0636   return {residual, residualCovariance};
0637 }
0638 
0639 /// Helper function to calculate the predicted chi2
0640 /// @tparam track_state_proxy_t the track state proxy type
0641 /// @param trackState the track state to calculate the chi2 from
0642 /// @return the chi2
0643 template <TrackStateProxyConcept track_state_proxy_t>
0644 double calculatePredictedChi2(track_state_proxy_t trackState) {
0645   if (!trackState.hasPredicted()) {
0646     throw std::invalid_argument("track state has no predicted parameters");
0647   }
0648   if (!trackState.hasCalibrated()) {
0649     throw std::invalid_argument("track state has no calibrated parameters");
0650   }
0651 
0652   return visit_measurement(
0653       trackState.calibratedSize(),
0654       [&]<std::size_t measdim>(
0655           std::integral_constant<std::size_t, measdim>) -> double {
0656         auto [residual, residualCovariance] =
0657             calculatePredictedResidual<measdim>(trackState);
0658 
0659         return (residual.transpose() * residualCovariance.inverse() * residual)
0660             .eval()(0, 0);
0661       });
0662 }
0663 
0664 /// Helper function to calculate the filtered chi2
0665 /// @tparam track_state_proxy_t the track state proxy type
0666 /// @param trackState the track state to calculate the chi2 from
0667 /// @return the chi2
0668 template <TrackStateProxyConcept track_state_proxy_t>
0669 double calculateFilteredChi2(track_state_proxy_t trackState) {
0670   if (!trackState.hasFiltered()) {
0671     throw std::invalid_argument("track state has no filtered parameters");
0672   }
0673   if (!trackState.hasCalibrated()) {
0674     throw std::invalid_argument("track state has no calibrated parameters");
0675   }
0676 
0677   return visit_measurement(
0678       trackState.calibratedSize(),
0679       [&]<std::size_t measdim>(
0680           std::integral_constant<std::size_t, measdim>) -> double {
0681         auto [residual, residualCovariance] =
0682             calculateFilteredResidual<measdim>(trackState);
0683 
0684         return (residual.transpose() * residualCovariance.inverse() * residual)
0685             .eval()(0, 0);
0686       });
0687 }
0688 
0689 /// Helper function to calculate the smoothed chi2
0690 /// @tparam track_state_proxy_t the track state proxy type
0691 /// @param trackState the track state to calculate the chi2 from
0692 /// @return the chi2
0693 template <TrackStateProxyConcept track_state_proxy_t>
0694 double calculateSmoothedChi2(track_state_proxy_t trackState) {
0695   if (!trackState.hasSmoothed()) {
0696     throw std::invalid_argument("track state has no smoothed parameters");
0697   }
0698   if (!trackState.hasCalibrated()) {
0699     throw std::invalid_argument("track state has no calibrated parameters");
0700   }
0701 
0702   return visit_measurement(
0703       trackState.calibratedSize(),
0704       [&]<std::size_t measdim>(
0705           std::integral_constant<std::size_t, measdim>) -> double {
0706         auto [residual, residualCovariance] =
0707             calculateSmoothedResidual<measdim>(trackState);
0708 
0709         return (residual.transpose() * residualCovariance.inverse() * residual)
0710             .eval()(0, 0);
0711       });
0712 }
0713 
0714 /// Helper function to calculate the unbiased track parameters and their
0715 /// covariance (i.e. fitted track parameters with this measurement removed)
0716 /// using Eq.(12a)-Eq.(12c) of NIMA 262, 444 (1987)
0717 /// @tparam track_state_proxy_t the track state proxy type
0718 /// @param trackState the track state to calculate the unbiased parameters from
0719 /// @return a pair of the unbiased parameters and their covariance
0720 template <TrackStateProxyConcept track_state_proxy_t>
0721 std::pair<BoundVector, BoundMatrix> calculateUnbiasedParametersCovariance(
0722     track_state_proxy_t trackState) {
0723   if (!trackState.hasSmoothed()) {
0724     throw std::invalid_argument("track state has no smoothed parameters");
0725   }
0726   if (!trackState.hasCalibrated()) {
0727     throw std::invalid_argument("track state has no calibrated parameters");
0728   }
0729 
0730   return visit_measurement(
0731       trackState.calibratedSize(),
0732       [&]<std::size_t measdim>(std::integral_constant<std::size_t, measdim>) {
0733         FixedBoundSubspaceHelper<measdim> subspaceHelper =
0734             trackState.template projectorSubspaceHelper<measdim>();
0735 
0736         // TODO use subspace helper for projection instead
0737         auto H = subspaceHelper.projector();
0738         auto s = trackState.smoothed();
0739         auto C = trackState.smoothedCovariance();
0740         auto m = trackState.template calibrated<measdim>();
0741         auto V = trackState.template calibratedCovariance<measdim>();
0742         auto K =
0743             (C * H.transpose() * (H * C * H.transpose() - V).inverse()).eval();
0744         BoundVector unbiasedParamsVec = s + K * (m - H * s);
0745         BoundMatrix unbiasedParamsCov = C - K * H * C;
0746         return std::make_pair(unbiasedParamsVec, unbiasedParamsCov);
0747       });
0748 }
0749 
0750 }  // namespace Acts
0751 
0752 namespace std {
0753 // register with STL
0754 template <>
0755 struct is_error_code_enum<Acts::TrackExtrapolationError> : std::true_type {};
0756 }  // namespace std