Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2026-10-24 08:17:47

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/Common.hpp"
0012 #include "Acts/EventData/MultiTrajectory.hpp"
0013 #include "Acts/EventData/MultiTrajectoryHelpers.hpp"
0014 #include "Acts/EventData/SourceLink.hpp"
0015 #include "Acts/EventData/VectorMultiTrajectory.hpp"
0016 #include "Acts/EventData/detail/CorrectedTransformationFreeToBound.hpp"
0017 #include "Acts/Geometry/GeometryContext.hpp"
0018 #include "Acts/MagneticField/MagneticFieldContext.hpp"
0019 #include "Acts/Propagator/ActorList.hpp"
0020 #include "Acts/Propagator/DirectNavigator.hpp"
0021 #include "Acts/Propagator/PropagatorOptions.hpp"
0022 #include "Acts/Propagator/StandardAborters.hpp"
0023 #include "Acts/Propagator/detail/LoopProtection.hpp"
0024 #include "Acts/Propagator/detail/PointwiseMaterialInteraction.hpp"
0025 #include "Acts/TrackFitting/KalmanFitterError.hpp"
0026 #include "Acts/TrackFitting/detail/VoidFitterComponents.hpp"
0027 #include "Acts/Utilities/CalibrationContext.hpp"
0028 #include "Acts/Utilities/Delegate.hpp"
0029 #include "Acts/Utilities/Logger.hpp"
0030 #include "Acts/Utilities/Result.hpp"
0031 #include "Acts/Utilities/TrackHelpers.hpp"
0032 
0033 #include <memory>
0034 #include <unordered_map>
0035 #include <vector>
0036 
0037 namespace Acts {
0038 
0039 /// @addtogroup track_fitting
0040 /// @{
0041 
0042 /// Extension struct which holds delegates to customise the KF behavior
0043 template <typename traj_t>
0044 struct KalmanFitterExtensions {
0045   /// Type alias for track state proxy from trajectory
0046   using TrackStateProxy = typename traj_t::TrackStateProxy;
0047   /// Type alias for const track state proxy from trajectory
0048   using ConstTrackStateProxy = typename traj_t::ConstTrackStateProxy;
0049   /// Type alias for track parameters from track state proxy
0050   using Parameters = typename TrackStateProxy::Parameters;
0051 
0052   /// Type alias for measurement calibrator delegate
0053   using Calibrator =
0054       Delegate<void(const GeometryContext&, const CalibrationContext&,
0055                     const SourceLink&, TrackStateProxy)>;
0056 
0057   /// Type alias for Kalman filter update delegate
0058   using Updater = Delegate<Result<void>(const GeometryContext&, TrackStateProxy,
0059                                         const Logger&)>;
0060 
0061   /// Type alias for outlier detection delegate
0062   using OutlierFinder = Delegate<bool(ConstTrackStateProxy)>;
0063 
0064   /// Type alias for reverse filtering decision delegate
0065   using ReverseFilteringLogic = Delegate<bool(ConstTrackStateProxy)>;
0066 
0067   /// Type alias for track smoothing delegate
0068   using Smoother = Delegate<Result<void>(const GeometryContext&, traj_t&,
0069                                          std::size_t, const Logger&)>;
0070 
0071   /// Retrieves the associated surface from a source link
0072   SourceLinkSurfaceAccessor surfaceAccessor;
0073 
0074   /// The Calibrator is a dedicated calibration algorithm that allows
0075   /// to calibrate measurements using track information, this could be
0076   /// e.g. sagging for wires, module deformations, etc.
0077   Calibrator calibrator;
0078 
0079   /// The updater incorporates measurement information into the track parameters
0080   Updater updater;
0081 
0082   /// Determines whether a measurement is supposed to be considered as an
0083   /// outlier
0084   OutlierFinder outlierFinder;
0085 
0086   /// Decides whether the smoothing stage uses linearized transport or full
0087   /// reverse propagation
0088   ReverseFilteringLogic reverseFilteringLogic;
0089 
0090   /// The smoother back-propagates measurement information along the track
0091   Smoother smoother;
0092 
0093   /// Default constructor which connects the default void components
0094   KalmanFitterExtensions() {
0095     surfaceAccessor.connect<&detail::voidSurfaceAccessor>();
0096     calibrator.template connect<&detail::voidFitterCalibrator<traj_t>>();
0097     updater.template connect<&detail::voidFitterUpdater<traj_t>>();
0098     outlierFinder.template connect<&detail::voidOutlierFinder<traj_t>>();
0099     reverseFilteringLogic
0100         .template connect<&detail::voidReverseFilteringLogic<traj_t>>();
0101     smoother.template connect<&detail::voidFitterSmoother<traj_t>>();
0102   }
0103 };
0104 
0105 /// Combined options for the Kalman fitter.
0106 ///
0107 /// @tparam traj_t The trajectory type
0108 template <typename traj_t>
0109 struct KalmanFitterOptions {
0110   /// PropagatorOptions with context.
0111   ///
0112   /// @param gctx The geometry context for this fit
0113   /// @param mctx The magnetic context for this fit
0114   /// @param cctx The calibration context for this fit
0115   /// @param extensions_ The KF extensions
0116   /// @param pOptions The plain propagator options
0117   /// @param tSurface The target surface for the fit
0118   /// @param mScattering Whether to include multiple scattering
0119   /// @param eLoss Whether to include energy loss
0120   /// @param rFiltering Whether to run reversed filtering
0121   /// @param rfScaling Scaling factor for covariance at input of reversed filtering
0122   /// @param freeToBoundCorrection_ Correction for non-linearity effect during transform from free to bound
0123   KalmanFitterOptions(const GeometryContext& gctx,
0124                       const MagneticFieldContext& mctx,
0125                       std::reference_wrapper<const CalibrationContext> cctx,
0126                       KalmanFitterExtensions<traj_t> extensions_,
0127                       const PropagatorPlainOptions& pOptions,
0128                       const Surface* tSurface = nullptr,
0129                       bool mScattering = true, bool eLoss = true,
0130                       bool rFiltering = false, double rfScaling = 1.0,
0131                       const FreeToBoundCorrection& freeToBoundCorrection_ =
0132                           FreeToBoundCorrection(false))
0133       : geoContext(gctx),
0134         magFieldContext(mctx),
0135         calibrationContext(cctx),
0136         extensions(extensions_),
0137         propagatorPlainOptions(pOptions),
0138         referenceSurface(tSurface),
0139         multipleScattering(mScattering),
0140         energyLoss(eLoss),
0141         reverseFiltering(rFiltering),
0142         reverseFilteringCovarianceScaling(rfScaling),
0143         freeToBoundCorrection(freeToBoundCorrection_) {}
0144 
0145   /// Context object for the geometry
0146   std::reference_wrapper<const GeometryContext> geoContext;
0147   /// Context object for the magnetic field
0148   std::reference_wrapper<const MagneticFieldContext> magFieldContext;
0149   /// context object for the calibration
0150   std::reference_wrapper<const CalibrationContext> calibrationContext;
0151 
0152   /// Extensions for calibration and outlier finding
0153   KalmanFitterExtensions<traj_t> extensions;
0154 
0155   /// The trivial propagator options
0156   PropagatorPlainOptions propagatorPlainOptions;
0157 
0158   /// The reference surface
0159   const Surface* referenceSurface = nullptr;
0160 
0161   /// Strategy to propagate to reference surface
0162   TrackExtrapolationStrategy referenceSurfaceStrategy =
0163       TrackExtrapolationStrategy::firstOrLast;
0164 
0165   /// Whether to consider multiple scattering
0166   bool multipleScattering = true;
0167 
0168   /// Whether to consider energy loss
0169   bool energyLoss = true;
0170 
0171   /// Whether to run filtering in reversed direction overwrite the
0172   /// ReverseFilteringLogic
0173   bool reverseFiltering = false;
0174 
0175   /// Factor by which the covariance of the input of the reversed filtering is
0176   /// scaled. This is only used in the backward filtering (if reverseFiltering
0177   /// is true or if the ReverseFilteringLogic return true for the track of
0178   /// interest).
0179   /// Note that the default value is not tuned and might need adjustment for
0180   /// different use cases.
0181   double reverseFilteringCovarianceScaling = 100.0;
0182 
0183   /// Whether to include non-linear correction during global to local
0184   /// transformation
0185   FreeToBoundCorrection freeToBoundCorrection;
0186 };
0187 
0188 /// Result payload returned by the Kalman fitter.
0189 template <typename traj_t>
0190 struct KalmanFitterResult {
0191   /// Fitted states that the actor has handled.
0192   traj_t* fittedStates{nullptr};
0193 
0194   /// This is the index of the 'tip' of the track stored in multitrajectory.
0195   /// This corresponds to the last measurement state in the multitrajectory.
0196   /// Since this KF only stores one trajectory, it is unambiguous.
0197   /// TrackTraits::kInvalid is the start of a trajectory.
0198   std::size_t lastMeasurementIndex = kTrackIndexInvalid;
0199 
0200   /// This is the index of the 'tip' of the states stored in multitrajectory.
0201   /// This corresponds to the last state in the multitrajectory.
0202   /// Since this KF only stores one trajectory, it is unambiguous.
0203   /// TrackTraits::kInvalid is the start of a trajectory.
0204   std::size_t lastTrackIndex = kTrackIndexInvalid;
0205 
0206   /// The optional Parameters at the provided surface
0207   std::optional<BoundTrackParameters> fittedParameters;
0208 
0209   /// Counter for states with non-outlier measurements
0210   std::size_t measurementStates = 0;
0211 
0212   /// Counter for measurements holes
0213   /// A hole correspond to a surface with an associated detector element with no
0214   /// associated measurement. Holes are only taken into account if they are
0215   /// between the first and last measurements.
0216   std::size_t measurementHoles = 0;
0217 
0218   /// Counter for handled states
0219   std::size_t processedStates = 0;
0220 
0221   /// Indicator if track fitting has been done
0222   bool finished = false;
0223 
0224   /// Measurement surfaces without hits
0225   std::vector<const Surface*> missedActiveSurfaces;
0226 
0227   /// Path limit aborter
0228   PathLimitReached pathLimitReached;
0229 };
0230 
0231 /// Kalman fitter implementation.
0232 ///
0233 /// @tparam propagator_t Type of the propagation class, usually an instance of
0234 ///         @ref Acts::Propagator
0235 ///
0236 /// The Kalman filter contains an Actor and a Sequencer sub-class.
0237 /// The Sequencer has to be part of the Navigator of the Propagator
0238 /// in order to initialize and provide the measurement surfaces.
0239 ///
0240 /// The Actor is part of the Propagation call and does the Kalman update
0241 /// and eventually the smoothing.  Updater, Smoother and Calibrator are
0242 /// given to the Actor for further use:
0243 /// - The Updater is the implemented kalman updater formalism, it
0244 ///   runs via a visitor pattern through the measurements.
0245 /// - The Smoother is called at the end of the filtering by the Actor.
0246 ///
0247 /// Measurements are not required to be ordered for the KalmanFilter,
0248 /// measurement ordering needs to be figured out by the navigation of
0249 /// the propagator.
0250 ///
0251 /// The void components are provided mainly for unit testing.
0252 template <typename propagator_t, typename traj_t>
0253 class KalmanFitter {
0254   /// The navigator type
0255   using KalmanNavigator = typename propagator_t::Navigator;
0256 
0257   /// The navigator has DirectNavigator type or not
0258   static constexpr bool isDirectNavigator =
0259       std::is_same_v<KalmanNavigator, DirectNavigator>;
0260 
0261  public:
0262   /// Constructor with propagator and logger
0263   /// @param pPropagator Propagator instance for track propagation
0264   /// @param _logger Logger for diagnostic output
0265   explicit KalmanFitter(propagator_t pPropagator,
0266                         std::unique_ptr<const Logger> _logger =
0267                             getDefaultLogger("KalmanFitter", Logging::INFO))
0268       : m_propagator(std::move(pPropagator)),
0269         m_logger{std::move(_logger)},
0270         m_actorLogger{m_logger->cloneWithSuffix("Actor")} {}
0271 
0272  private:
0273   /// The propagator for the transport and material update
0274   propagator_t m_propagator;
0275 
0276   /// The logger instance
0277   std::unique_ptr<const Logger> m_logger;
0278   std::unique_ptr<const Logger> m_actorLogger;
0279 
0280   const Logger& logger() const { return *m_logger; }
0281 
0282   /// @brief Propagator Actor plugin for the KalmanFilter
0283   ///
0284   /// The KalmanActor does not rely on the measurements to be
0285   /// sorted along the track.
0286   class Actor {
0287    public:
0288     /// Broadcast the result_type
0289     using result_type = KalmanFitterResult<traj_t>;
0290 
0291     /// The target surface aborter
0292     SurfaceReached targetReached{std::numeric_limits<double>::lowest()};
0293 
0294     /// Allows retrieving measurements for a surface
0295     std::unordered_map<const Surface*, SourceLink> inputMeasurements;
0296 
0297     /// Whether to consider multiple scattering.
0298     bool multipleScattering = true;
0299 
0300     /// Whether to consider energy loss.
0301     bool energyLoss = true;
0302 
0303     /// Whether to include non-linear correction during global to local
0304     /// transformation
0305     FreeToBoundCorrection freeToBoundCorrection;
0306 
0307     /// Input MultiTrajectory
0308     std::shared_ptr<traj_t> outputStates;
0309 
0310     KalmanFitterExtensions<traj_t> extensions;
0311 
0312     /// Calibration context for the fit
0313     const CalibrationContext* calibrationContext{nullptr};
0314 
0315     /// End of world aborter
0316     EndOfWorldReached endOfWorldReached;
0317 
0318     /// Volume constraint aborter
0319     VolumeConstraintAborter volumeConstraintAborter;
0320 
0321     /// The logger instance
0322     const Logger* actorLogger{nullptr};
0323 
0324     /// Logger helper
0325     const Logger& logger() const { return *actorLogger; }
0326 
0327     /// @brief Kalman actor operation
0328     ///
0329     /// @tparam propagator_state_t is the type of Propagator state
0330     /// @tparam stepper_t Type of the stepper
0331     /// @tparam navigator_t Type of the navigator
0332     ///
0333     /// @param state is the mutable propagator state object
0334     /// @param stepper The stepper in use
0335     /// @param navigator The navigator in use
0336     /// @param result is the mutable result state object
0337     template <typename propagator_state_t, typename stepper_t,
0338               typename navigator_t>
0339     Result<void> act(propagator_state_t& state, const stepper_t& stepper,
0340                      const navigator_t& navigator, result_type& result,
0341                      const Logger& /*logger*/) const {
0342       assert(result.fittedStates && "No MultiTrajectory set");
0343 
0344       if (result.finished) {
0345         return Result<void>::success();
0346       }
0347 
0348       ACTS_VERBOSE("KalmanFitter step at pos: "
0349                    << stepper.position(state.stepping).transpose()
0350                    << " dir: " << stepper.direction(state.stepping).transpose()
0351                    << " momentum: "
0352                    << stepper.absoluteMomentum(state.stepping));
0353 
0354       // Initialize path limit reached aborter
0355       if (result.pathLimitReached.internalLimit ==
0356           std::numeric_limits<double>::max()) {
0357         detail::setupLoopProtection(state, stepper, result.pathLimitReached,
0358                                     true, logger());
0359       }
0360 
0361       // Update:
0362       // - Waiting for a current surface
0363       const Surface* surface = navigator.currentSurface(state.navigation);
0364       if (surface != nullptr) {
0365         // Check if the surface is in the measurement map
0366         // -> Get the measurement / calibrate
0367         // -> Create the predicted state
0368         // -> Check outlier behavior, if non-outlier:
0369         // -> Perform the kalman update
0370         // -> Fill track state information & update stepper information
0371 
0372         ACTS_VERBOSE("Perform " << state.options.direction << " filter step");
0373         auto res = filter(*surface, state, stepper, navigator, result);
0374         if (!res.ok()) {
0375           ACTS_DEBUG("Error in " << state.options.direction
0376                                  << " filter: " << res.error());
0377           return res.error();
0378         }
0379       }
0380 
0381       // Finalization:
0382       // when all track states have been handled or an aborter is triggered
0383       const bool isTrackComplete =
0384           result.measurementStates == inputMeasurements.size();
0385       const bool isEndOfWorldReached =
0386           endOfWorldReached.checkAbort(state, stepper, navigator, logger());
0387       const bool isVolumeConstraintReached = volumeConstraintAborter.checkAbort(
0388           state, stepper, navigator, logger());
0389       const bool isPathLimitReached = result.pathLimitReached.checkAbort(
0390           state, stepper, navigator, logger());
0391       const bool isTargetReached =
0392           targetReached.checkAbort(state, stepper, navigator, logger());
0393       if (isTrackComplete || isEndOfWorldReached || isVolumeConstraintReached ||
0394           isPathLimitReached || isTargetReached) {
0395         ACTS_VERBOSE(
0396             "Finalizing Kalman fit: "
0397             << (isTrackComplete ? "track complete; " : "")
0398             << (isEndOfWorldReached ? "end of world reached; " : "")
0399             << (isVolumeConstraintReached ? "volume constraint reached; " : "")
0400             << (isPathLimitReached ? "path limit reached; " : "")
0401             << (isTargetReached ? "target surface reached; " : ""));
0402 
0403         if (isTargetReached) {
0404           ACTS_VERBOSE("Setting fitted parameters at target surface");
0405 
0406           // Bind the parameter to the target surface
0407           auto transportRes =
0408               stepper.transportToBound(state.stepping, *targetReached.surface);
0409           if (!transportRes.ok()) {
0410             ACTS_DEBUG("Error while transporting to the target surface: "
0411                        << transportRes.error() << " "
0412                        << transportRes.error().message());
0413             return transportRes.error();
0414           }
0415           auto res =
0416               stepper.boundParameters(state.stepping, *targetReached.surface);
0417           if (!res.ok()) {
0418             ACTS_DEBUG("Error while binding to the target surface: "
0419                        << res.error() << " " << res.error().message());
0420             return res.error();
0421           }
0422           result.fittedParameters = std::move(*res);
0423         }
0424 
0425         result.finished = true;
0426       }
0427 
0428       return Result<void>::success();
0429     }
0430 
0431     template <typename propagator_state_t, typename stepper_t,
0432               typename navigator_t>
0433     bool checkAbort(propagator_state_t& /*state*/, const stepper_t& /*stepper*/,
0434                     const navigator_t& /*navigator*/, const result_type& result,
0435                     const Logger& /*logger*/) const {
0436       return result.finished;
0437     }
0438 
0439     /// @brief Kalman actor operation: update
0440     ///
0441     /// @tparam propagator_state_t is the type of Propagator state
0442     /// @tparam stepper_t Type of the stepper
0443     /// @tparam navigator_t Type of the navigator
0444     ///
0445     /// @param surface The surface where the update happens
0446     /// @param state The mutable propagator state object
0447     /// @param stepper The stepper in use
0448     /// @param navigator The navigator in use
0449     /// @param result The mutable result state object
0450     template <typename propagator_state_t, typename stepper_t,
0451               typename navigator_t>
0452     Result<void> filter(const Surface& surface, propagator_state_t& state,
0453                         const stepper_t& stepper, const navigator_t& navigator,
0454                         result_type& result) const {
0455       const bool precedingMeasurementExists = result.measurementStates > 0;
0456       const bool surfaceIsSensitive = surface.isSensitive();
0457       const bool surfaceHasMaterial = surface.hasMaterial();
0458 
0459       // Try to find the surface in the measurement surfaces
0460       const auto sourceLinkIt = inputMeasurements.find(&surface);
0461       if (sourceLinkIt != inputMeasurements.end()) {
0462         // Screen output message
0463         ACTS_VERBOSE("Measurement surface " << surface.geometryId()
0464                                             << " detected.");
0465         // Transport the covariance to the surface
0466         auto transportRes = stepper.transportToBound(state.stepping, surface,
0467                                                      freeToBoundCorrection);
0468         if (!transportRes.ok()) {
0469           return transportRes.error();
0470         }
0471 
0472         // Update state and stepper with pre material effects
0473         const Result<detail::PointwiseMaterialEffects>
0474             materialInteractionPreRes = detail::performMaterialInteraction(
0475                 state, stepper, surface,
0476                 detail::determineMaterialUpdateMode(
0477                     state, navigator, MaterialUpdateMode::PreUpdate),
0478                 NoiseUpdateMode::addNoise, multipleScattering, energyLoss,
0479                 logger());
0480         if (!materialInteractionPreRes.ok()) {
0481           ACTS_DEBUG("Material interaction failed during filter: "
0482                      << materialInteractionPreRes.error().message());
0483           return materialInteractionPreRes.error();
0484         }
0485 
0486         // Create a track state with the desired components. Note that the
0487         // filtered parameters are deliberately not allocated here: they are
0488         // only added once the Kalman update is about to write them, so that
0489         // the calibrator and the outlier finder cannot observe allocated but
0490         // uninitialized filtered parameters via `parameters()`.
0491         TrackStatePropMask mask = TrackStatePropMask::Predicted |
0492                                   TrackStatePropMask::Jacobian |
0493                                   TrackStatePropMask::Calibrated;
0494         typename traj_t::TrackStateProxy trackStateProxy =
0495             result.fittedStates->makeTrackState(mask, result.lastTrackIndex);
0496 
0497         typename traj_t::ConstTrackStateProxy trackStateProxyConst{
0498             trackStateProxy};
0499 
0500         // Set the trackStateProxy components with the state from the ongoing
0501         // propagation
0502         trackStateProxy.setReferenceSurface(surface.getSharedPtr());
0503         // Bind the transported state to the current surface
0504         auto res = stepper.boundParameters(state.stepping, surface);
0505         if (!res.ok()) {
0506           ACTS_DEBUG("Propagate to surface " << surface.geometryId()
0507                                              << " failed: " << res.error());
0508           return res.error();
0509         }
0510 
0511         // Fill the track state
0512         trackStateProxy.predicted() = res->parameters();
0513         trackStateProxy.predictedCovariance() =
0514             stepper.covariance(state.stepping).value();
0515 
0516         trackStateProxy.jacobian() = *transportRes;
0517         trackStateProxy.pathLength() = stepper.pathLength(state.stepping);
0518 
0519         // We have predicted parameters, so calibrate the uncalibrated input
0520         // measurement
0521         extensions.calibrator(state.geoContext, *calibrationContext,
0522                               sourceLinkIt->second, trackStateProxy);
0523 
0524         // Get and set the type flags
0525         auto typeFlags = trackStateProxy.typeFlags();
0526         typeFlags.setHasParameters();
0527         if (surface.hasMaterial()) {
0528           typeFlags.setHasMaterial();
0529         }
0530 
0531         // Check if the state is an outlier.
0532         // If not:
0533         // - run Kalman update
0534         // - tag it as a measurement
0535         // - update the stepping state.
0536         // Else, just tag it as an outlier
0537         if (!extensions.outlierFinder(trackStateProxyConst)) {
0538           // Allocate the filtered parameters right before they are written
0539           trackStateProxy.addComponents(TrackStatePropMask::Filtered);
0540           // Run Kalman update
0541           auto updateRes =
0542               extensions.updater(state.geoContext, trackStateProxy, logger());
0543           if (!updateRes.ok()) {
0544             ACTS_DEBUG("Update step failed: " << updateRes.error());
0545             return updateRes.error();
0546           }
0547           // Set the measurement type flag
0548           typeFlags.setIsMeasurement();
0549         } else {
0550           ACTS_VERBOSE(
0551               "Filtering step successful. But measurement is determined "
0552               "to be an outlier. Stepping state is not updated.");
0553           // Set the outlier type flag
0554           typeFlags.setIsOutlier();
0555           trackStateProxy.shareFrom(trackStateProxy,
0556                                     TrackStatePropMask::Predicted,
0557                                     TrackStatePropMask::Filtered);
0558         }
0559 
0560         result.lastTrackIndex = trackStateProxy.index();
0561 
0562         // Update the stepper if it is not an outlier
0563         if (trackStateProxy.typeFlags().isMeasurement()) {
0564           // Update the stepping state with filtered parameters
0565           ACTS_VERBOSE("Filtering step successful, updated parameters are:\n"
0566                        << trackStateProxy.filtered().transpose());
0567           // update stepping state using filtered parameters after kalman
0568           stepper.update(state.stepping,
0569                          MultiTrajectoryHelpers::freeFiltered(
0570                              state.options.geoContext, trackStateProxy),
0571                          trackStateProxy.filtered(),
0572                          trackStateProxy.filteredCovariance(), surface);
0573           // We count the state with measurement
0574           ++result.measurementStates;
0575         }
0576 
0577         // Update state and stepper with post material effects
0578         const Result<detail::PointwiseMaterialEffects>
0579             materialInteractionPostRes = detail::performMaterialInteraction(
0580                 state, stepper, surface,
0581                 detail::determineMaterialUpdateMode(
0582                     state, navigator, MaterialUpdateMode::PostUpdate),
0583                 NoiseUpdateMode::addNoise, multipleScattering, energyLoss,
0584                 logger());
0585         if (!materialInteractionPostRes.ok()) {
0586           ACTS_DEBUG("Material interaction failed during filter: "
0587                      << materialInteractionPostRes.error().message());
0588           return materialInteractionPostRes.error();
0589         }
0590         // We count the processed state
0591         ++result.processedStates;
0592         // Update the number of holes count only when encountering a
0593         // measurement
0594         result.measurementHoles = result.missedActiveSurfaces.size();
0595         // Since we encountered a measurement update the lastMeasurementIndex to
0596         // the lastTrackIndex.
0597         result.lastMeasurementIndex = result.lastTrackIndex;
0598 
0599       } else if ((precedingMeasurementExists && surfaceIsSensitive) ||
0600                  surfaceHasMaterial) {
0601         // We only create track states here if there is already measurement
0602         // detected or if the surface has material (no holes before the first
0603         // measurement)
0604 
0605         // Create a track state with the desired components
0606         TrackStatePropMask mask =
0607             TrackStatePropMask::Predicted | TrackStatePropMask::Jacobian;
0608         typename traj_t::TrackStateProxy trackStateProxy =
0609             result.fittedStates->makeTrackState(mask, result.lastTrackIndex);
0610 
0611         // Set the trackStateProxy components with the state from the ongoing
0612         // propagation
0613         trackStateProxy.setReferenceSurface(surface.getSharedPtr());
0614         // Transport the covariance to the surface and bind the state to it
0615         auto transportRes = stepper.transportToBound(state.stepping, surface,
0616                                                      freeToBoundCorrection);
0617         if (!transportRes.ok()) {
0618           return transportRes.error();
0619         }
0620         auto res = stepper.boundParameters(state.stepping, surface);
0621         if (!res.ok()) {
0622           return res.error();
0623         }
0624 
0625         // Fill the track state
0626         trackStateProxy.predicted() = res->parameters();
0627         trackStateProxy.predictedCovariance() =
0628             stepper.covariance(state.stepping).value();
0629 
0630         trackStateProxy.jacobian() = *transportRes;
0631         trackStateProxy.pathLength() = stepper.pathLength(state.stepping);
0632 
0633         // Set the filtered parameter index to be the same with predicted
0634         // parameter
0635         trackStateProxy.shareFrom(trackStateProxy,
0636                                   TrackStatePropMask::Predicted,
0637                                   TrackStatePropMask::Filtered);
0638 
0639         // Set the track state flags
0640         auto typeFlags = trackStateProxy.typeFlags();
0641         typeFlags.setHasParameters();
0642 
0643         if (surfaceHasMaterial) {
0644           typeFlags.setHasMaterial();
0645         }
0646 
0647         if (surfaceIsSensitive && precedingMeasurementExists) {
0648           ACTS_VERBOSE("Detected hole on " << surface.geometryId());
0649           // If the surface is sensitive, set the hole type flag
0650           typeFlags.setIsHole();
0651         } else if (surfaceIsSensitive) {
0652           ACTS_VERBOSE("Skip hole (no preceding measurements) on surface "
0653                        << surface.geometryId());
0654         } else if (surfaceHasMaterial) {
0655           ACTS_VERBOSE("Detected in-sensitive surface "
0656                        << surface.geometryId());
0657         }
0658 
0659         result.lastTrackIndex = trackStateProxy.index();
0660 
0661         if (trackStateProxy.typeFlags().isHole()) {
0662           // Count the missed surface
0663           result.missedActiveSurfaces.push_back(&surface);
0664         }
0665 
0666         ++result.processedStates;
0667 
0668         // Update state and stepper with (possible) material effects
0669         const Result<detail::PointwiseMaterialEffects> materialInteractionRes =
0670             detail::performMaterialInteraction(
0671                 state, stepper, surface,
0672                 detail::determineMaterialUpdateMode(
0673                     state, navigator, MaterialUpdateMode::FullUpdate),
0674                 NoiseUpdateMode::addNoise, multipleScattering, energyLoss,
0675                 logger());
0676         if (!materialInteractionRes.ok()) {
0677           ACTS_DEBUG("Material interaction failed during filter: "
0678                      << materialInteractionRes.error().message());
0679           return materialInteractionRes.error();
0680         }
0681       }
0682 
0683       return Result<void>::success();
0684     }
0685   };
0686 
0687  public:
0688   /// Fit implementation of the forward filter, calls the
0689   /// the filter and smoother/reversed filter
0690   ///
0691   /// @tparam source_link_iterator_t Iterator type used to pass source links
0692   /// @tparam track_container_t Type of the track container
0693   ///
0694   /// @param it Begin iterator for the fittable uncalibrated measurements
0695   /// @param end End iterator for the fittable uncalibrated measurements
0696   /// @param sParameters The initial track parameters
0697   /// @param kfOptions KalmanOptions steering the fit
0698   /// @param trackContainer Input track container storage to append into
0699   /// @note The input measurements are given in the form of @c SourceLink s.
0700   /// It's the calibrators job to turn them into calibrated measurements used in
0701   /// the fit.
0702   ///
0703   /// @return the output as an output track
0704   template <typename source_link_iterator_t,
0705             TrackContainerFrontend track_container_t>
0706   Result<typename track_container_t::TrackProxy> fit(
0707       source_link_iterator_t it, source_link_iterator_t end,
0708       const BoundTrackParameters& sParameters,
0709       const KalmanFitterOptions<traj_t>& kfOptions,
0710       track_container_t& trackContainer) const {
0711     return fit_impl(it, end, sParameters, kfOptions, nullptr, trackContainer);
0712   }
0713 
0714   /// Fit implementation of the forward filter, calls the
0715   /// the filter and smoother/reversed filter
0716   ///
0717   /// @tparam source_link_iterator_t Iterator type used to pass source links
0718   /// @tparam track_container_t Type of the track container
0719   ///
0720   /// @param it Begin iterator for the fittable uncalibrated measurements
0721   /// @param end End iterator for the fittable uncalibrated measurements
0722   /// @param sParameters The initial track parameters
0723   /// @param kfOptions KalmanOptions steering the fit
0724   /// @param sSequence surface sequence used to initialize a DirectNavigator
0725   /// @param trackContainer Input track container storage to append into
0726   /// @note The input measurements are given in the form of @c SourceLinks.
0727   /// It's
0728   /// @c calibrator_t's job to turn them into calibrated measurements used in
0729   /// the fit.
0730   ///
0731   /// @return the output as an output track
0732   template <typename source_link_iterator_t,
0733             TrackContainerFrontend track_container_t>
0734   Result<typename track_container_t::TrackProxy> fit(
0735       source_link_iterator_t it, source_link_iterator_t end,
0736       const BoundTrackParameters& sParameters,
0737       const KalmanFitterOptions<traj_t>& kfOptions,
0738       const std::vector<const Surface*>& sSequence,
0739       track_container_t& trackContainer) const
0740     requires(isDirectNavigator)
0741   {
0742     return fit_impl(it, end, sParameters, kfOptions, &sSequence,
0743                     trackContainer);
0744   }
0745 
0746  private:
0747   template <typename source_link_iterator_t>
0748   auto make_propagator_options(source_link_iterator_t it,
0749                                source_link_iterator_t end,
0750                                const KalmanFitterOptions<traj_t>& kfOptions,
0751                                const std::vector<const Surface*>* sSequence,
0752                                const Surface* targetSurface,
0753                                bool reverseDirection) const {
0754     using KalmanActor = Actor;
0755     using Actors = ActorList<KalmanActor>;
0756     using PropagatorOptions = typename propagator_t::template Options<Actors>;
0757 
0758     const std::size_t nMeasurements = std::distance(it, end);
0759 
0760     // To be able to find measurements later, we put them into a map
0761     // We need to copy input SourceLinks anyway, so the map can own them.
0762     ACTS_VERBOSE("Preparing " << nMeasurements << " input measurements");
0763     std::unordered_map<const Surface*, SourceLink> inputMeasurements;
0764     for (; it != end; ++it) {
0765       SourceLink sl = *it;
0766       const Surface* surface = kfOptions.extensions.surfaceAccessor(sl);
0767       inputMeasurements.try_emplace(surface, std::move(sl));
0768     }
0769 
0770     // Create relevant options for the propagation options
0771     PropagatorOptions propagatorOptions(kfOptions.geoContext,
0772                                         kfOptions.magFieldContext);
0773 
0774     // Set the trivial propagator options
0775     propagatorOptions.setPlainOptions(kfOptions.propagatorPlainOptions);
0776 
0777     if (reverseDirection) {
0778       propagatorOptions.direction = propagatorOptions.direction.invert();
0779     }
0780 
0781     if constexpr (!isDirectNavigator) {
0782       // Relax the bounds check, so the navigator targets a measurement
0783       // surface even where the track misses it
0784       for (const auto& [surface, _] : inputMeasurements) {
0785         propagatorOptions.navigation.registerMeasurementSurface(*surface);
0786       }
0787     } else {
0788       assert(sSequence != nullptr &&
0789              "DirectNavigator requires a surface sequence for KalmanFitter");
0790       // Set the surface sequence
0791       propagatorOptions.navigation.surfaceSequence = *sSequence;
0792     }
0793 
0794     // Catch the actor and set the measurements
0795     auto& kalmanActor = propagatorOptions.actorList.template get<KalmanActor>();
0796     kalmanActor.inputMeasurements = std::move(inputMeasurements);
0797     kalmanActor.targetReached.surface = targetSurface;
0798     kalmanActor.multipleScattering = kfOptions.multipleScattering;
0799     kalmanActor.energyLoss = kfOptions.energyLoss;
0800     kalmanActor.freeToBoundCorrection = kfOptions.freeToBoundCorrection;
0801     kalmanActor.calibrationContext = &kfOptions.calibrationContext.get();
0802     kalmanActor.extensions = kfOptions.extensions;
0803     kalmanActor.actorLogger = m_actorLogger.get();
0804 
0805     return propagatorOptions;
0806   }
0807 
0808   template <typename propagator_options_t,
0809             TrackContainerFrontend track_container_t>
0810   auto filter_impl(const BoundTrackParameters& sParameters,
0811                    const propagator_options_t& propagatorOptions,
0812                    track_container_t& trackContainer) const
0813       -> Result<typename track_container_t::TrackProxy> {
0814     auto propagatorState = m_propagator.makeState(propagatorOptions);
0815 
0816     auto propagatorInitResult =
0817         m_propagator.initialize(propagatorState, sParameters, nullptr);
0818     if (!propagatorInitResult.ok()) {
0819       ACTS_DEBUG("Propagation initialization failed: "
0820                  << propagatorInitResult.error());
0821       return propagatorInitResult.error();
0822     }
0823 
0824     auto& kalmanResult =
0825         propagatorState.template get<KalmanFitterResult<traj_t>>();
0826     kalmanResult.fittedStates = &trackContainer.trackStateContainer();
0827 
0828     // Run the fitter
0829     auto result = m_propagator.propagate(propagatorState);
0830 
0831     if (!result.ok()) {
0832       ACTS_DEBUG("Propagation failed: " << result.error());
0833       return result.error();
0834     }
0835 
0836     /// It could happen that the fit ends in zero measurement states.
0837     /// The result gets meaningless so such case is regarded as fit failure.
0838     if (!kalmanResult.measurementStates) {
0839       ACTS_DEBUG("KalmanFilter failed: No measurement states found");
0840       return KalmanFitterError::NoMeasurementFound;
0841     }
0842 
0843     auto track = trackContainer.makeTrack();
0844     track.tipIndex() = kalmanResult.lastMeasurementIndex;
0845     track.setParticleHypothesis(sParameters.particleHypothesis());
0846     if (kalmanResult.fittedParameters) {
0847       const auto& params = kalmanResult.fittedParameters.value();
0848       track.parameters() = params.parameters();
0849       track.covariance() = params.covariance().value();
0850       track.setReferenceSurface(params.referenceSurface().getSharedPtr());
0851     }
0852 
0853     calculateTrackQuantities(track);
0854 
0855     track.linkForward();
0856 
0857     return track;
0858   }
0859 
0860   /// Common fit implementation
0861   ///
0862   /// @tparam source_link_iterator_t Iterator type used to pass source links
0863   /// @tparam track_container_t Type of the track container
0864   ///
0865   /// @param it Begin iterator for the fittable uncalibrated measurements
0866   /// @param end End iterator for the fittable uncalibrated measurements
0867   /// @param sParameters The initial track parameters
0868   /// @param kfOptions KalmanOptions steering the fit
0869   /// @param sSequence surface sequence used to initialize a DirectNavigator
0870   /// @param trackContainer Input track container storage to append into
0871   ///
0872   /// @return the output as an output track
0873   template <typename source_link_iterator_t,
0874             TrackContainerFrontend track_container_t>
0875   auto fit_impl(source_link_iterator_t it, source_link_iterator_t end,
0876                 const BoundTrackParameters& sParameters,
0877                 const KalmanFitterOptions<traj_t>& kfOptions,
0878                 const std::vector<const Surface*>* sSequence,
0879                 track_container_t& trackContainer) const
0880       -> Result<typename track_container_t::TrackProxy> {
0881     using TrackProxy = typename track_container_t::TrackProxy;
0882     using TrackStateProxy = typename track_container_t::TrackStateProxy;
0883 
0884     auto forwardPropagatorOptions =
0885         make_propagator_options(it, end, kfOptions, sSequence, nullptr, false);
0886 
0887     auto forwardFilterResult =
0888         filter_impl(sParameters, forwardPropagatorOptions, trackContainer);
0889 
0890     if (!forwardFilterResult.ok()) {
0891       ACTS_DEBUG("KalmanFilter failed: "
0892                  << forwardFilterResult.error() << ", "
0893                  << forwardFilterResult.error().message());
0894       return forwardFilterResult.error();
0895     }
0896 
0897     TrackProxy forwardTrack = forwardFilterResult.value();
0898 
0899     TrackStateProxy firstMeasurementState =
0900         trackContainer.trackStateContainer().getTrackState(
0901             findFirstMeasurementState(forwardTrack).value().index());
0902     TrackStateProxy lastMeasurementState =
0903         trackContainer.trackStateContainer().getTrackState(
0904             findLastMeasurementState(forwardTrack).value().index());
0905     lastMeasurementState.shareFrom(lastMeasurementState,
0906                                    TrackStatePropMask::Filtered,
0907                                    TrackStatePropMask::Smoothed);
0908 
0909     TrackProxy track = forwardTrack;
0910 
0911     const bool doReverseFilter =
0912         kfOptions.reverseFiltering ||
0913         kfOptions.extensions.reverseFilteringLogic(
0914             typename traj_t::ConstTrackStateProxy(lastMeasurementState));
0915     if (doReverseFilter) {
0916       ACTS_VERBOSE("Smooth track by reversed filtering");
0917 
0918       auto reverseStartParameters = forwardTrack.createParametersFromState(
0919           typename traj_t::ConstTrackStateProxy(lastMeasurementState));
0920       reverseStartParameters.covariance().value() *=
0921           kfOptions.reverseFilteringCovarianceScaling;
0922       auto reversePropagatorOptions = make_propagator_options(
0923           it, end, kfOptions, sSequence, kfOptions.referenceSurface, true);
0924       auto reverseFilterResult = filter_impl(
0925           reverseStartParameters, reversePropagatorOptions, trackContainer);
0926 
0927       if (!reverseFilterResult.ok()) {
0928         ACTS_DEBUG("Reversed KalmanFilter failed: "
0929                    << reverseFilterResult.error() << ", "
0930                    << reverseFilterResult.error().message());
0931         return reverseFilterResult.error();
0932       }
0933 
0934       TrackProxy reverseTrack = reverseFilterResult.value();
0935 
0936       TrackStateProxy reverseLastMeasurementState =
0937           trackContainer.trackStateContainer().getTrackState(
0938               findLastMeasurementState(reverseTrack).value().index());
0939 
0940       if (&firstMeasurementState.referenceSurface() !=
0941           &reverseLastMeasurementState.referenceSurface()) {
0942         ACTS_DEBUG(
0943             "Inconsistent reference surfaces between forward and "
0944             "reversed filtered tracks");
0945         return Result<TrackProxy>::failure(
0946             KalmanFitterError::InconsistentTrackStates);
0947       }
0948       firstMeasurementState.shareFrom(reverseLastMeasurementState,
0949                                       TrackStatePropMask::Filtered,
0950                                       TrackStatePropMask::Smoothed);
0951 
0952       if (reverseTrack.hasReferenceSurface()) {
0953         track.parameters() = reverseTrack.parameters();
0954         track.covariance() = reverseTrack.covariance();
0955         track.setReferenceSurface(
0956             reverseTrack.referenceSurface().getSharedPtr());
0957       }
0958 
0959       trackContainer.removeTrack(reverseTrack.index());
0960     } else {
0961       ACTS_VERBOSE("Smooth track directly without reversed filtering");
0962 
0963       auto smoothRes = kfOptions.extensions.smoother(
0964           kfOptions.geoContext, trackContainer.trackStateContainer(),
0965           forwardTrack.tipIndex(), logger());
0966       if (!smoothRes.ok()) {
0967         ACTS_DEBUG("Smoothing step failed: " << smoothRes.error() << ", "
0968                                              << smoothRes.error().message());
0969         return smoothRes.error();
0970       }
0971     }
0972 
0973     if (!track.hasReferenceSurface() && kfOptions.referenceSurface != nullptr) {
0974       typename propagator_t::template Options<> extrapolationOptions(
0975           kfOptions.geoContext, kfOptions.magFieldContext);
0976       auto extrapolationResult = extrapolateTrackToReferenceSurface(
0977           track, *kfOptions.referenceSurface, m_propagator,
0978           extrapolationOptions, kfOptions.referenceSurfaceStrategy, logger());
0979 
0980       if (!extrapolationResult.ok()) {
0981         ACTS_DEBUG("Extrapolation to reference surface failed: "
0982                    << extrapolationResult.error() << ", "
0983                    << extrapolationResult.error().message());
0984         return extrapolationResult.error();
0985       }
0986     }
0987 
0988     if (trackContainer.hasColumn(hashString("smoothed"))) {
0989       track.template component<bool, hashString("smoothed")>() =
0990           !doReverseFilter;
0991     }
0992     if (trackContainer.hasColumn(hashString("reversed"))) {
0993       track.template component<bool, hashString("reversed")>() =
0994           doReverseFilter;
0995     }
0996 
0997     return track;
0998   }
0999 };
1000 
1001 /// @}
1002 
1003 }  // namespace Acts