Back to home page

EIC code displayed by LXR

 
 

    


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