Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2026-09-16 08:19:30

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