Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-01-18 09:11:12

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