Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2026-08-25 08:16:39

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 res = stepper.boundState(state.stepping, *targetReached.surface);
0408           if (!res.ok()) {
0409             ACTS_DEBUG("Error while acquiring bound state for target surface: "
0410                        << res.error() << " " << res.error().message());
0411             return res.error();
0412           } else {
0413             const auto& [boundParams, jacobian, pathLength] = *res;
0414             result.fittedParameters = boundParams;
0415           }
0416         }
0417 
0418         result.finished = true;
0419       }
0420 
0421       return Result<void>::success();
0422     }
0423 
0424     template <typename propagator_state_t, typename stepper_t,
0425               typename navigator_t>
0426     bool checkAbort(propagator_state_t& /*state*/, const stepper_t& /*stepper*/,
0427                     const navigator_t& /*navigator*/, const result_type& result,
0428                     const Logger& /*logger*/) const {
0429       return result.finished;
0430     }
0431 
0432     /// @brief Kalman actor operation: update
0433     ///
0434     /// @tparam propagator_state_t is the type of Propagator state
0435     /// @tparam stepper_t Type of the stepper
0436     /// @tparam navigator_t Type of the navigator
0437     ///
0438     /// @param surface The surface where the update happens
0439     /// @param state The mutable propagator state object
0440     /// @param stepper The stepper in use
0441     /// @param navigator The navigator in use
0442     /// @param result The mutable result state object
0443     template <typename propagator_state_t, typename stepper_t,
0444               typename navigator_t>
0445     Result<void> filter(const Surface& surface, propagator_state_t& state,
0446                         const stepper_t& stepper, const navigator_t& navigator,
0447                         result_type& result) const {
0448       const bool precedingMeasurementExists = result.measurementStates > 0;
0449       const bool surfaceIsSensitive = surface.isSensitive();
0450       const bool surfaceHasMaterial = surface.hasMaterial();
0451 
0452       // Try to find the surface in the measurement surfaces
0453       const auto sourceLinkIt = inputMeasurements.find(&surface);
0454       if (sourceLinkIt != inputMeasurements.end()) {
0455         // Screen output message
0456         ACTS_VERBOSE("Measurement surface " << surface.geometryId()
0457                                             << " detected.");
0458         // Transport the covariance to the surface
0459         stepper.transportCovarianceToBound(state.stepping, surface,
0460                                            freeToBoundCorrection);
0461 
0462         // Update state and stepper with pre material effects
0463         const Result<detail::PointwiseMaterialEffects>
0464             materialInteractionPreRes = detail::performMaterialInteraction(
0465                 state, stepper, surface,
0466                 detail::determineMaterialUpdateMode(
0467                     state, navigator, MaterialUpdateMode::PreUpdate),
0468                 NoiseUpdateMode::addNoise, multipleScattering, energyLoss,
0469                 logger());
0470         if (!materialInteractionPreRes.ok()) {
0471           ACTS_DEBUG("Material interaction failed during filter: "
0472                      << materialInteractionPreRes.error().message());
0473           return materialInteractionPreRes.error();
0474         }
0475 
0476         // Create a track state with the desired components. Note that the
0477         // filtered parameters are deliberately not allocated here: they are
0478         // only added once the Kalman update is about to write them, so that
0479         // the calibrator and the outlier finder cannot observe allocated but
0480         // uninitialized filtered parameters via `parameters()`.
0481         TrackStatePropMask mask = TrackStatePropMask::Predicted |
0482                                   TrackStatePropMask::Jacobian |
0483                                   TrackStatePropMask::Calibrated;
0484         typename traj_t::TrackStateProxy trackStateProxy =
0485             result.fittedStates->makeTrackState(mask, result.lastTrackIndex);
0486 
0487         typename traj_t::ConstTrackStateProxy trackStateProxyConst{
0488             trackStateProxy};
0489 
0490         // Set the trackStateProxy components with the state from the ongoing
0491         // propagation
0492         trackStateProxy.setReferenceSurface(surface.getSharedPtr());
0493         // Bind the transported state to the current surface
0494         auto res = stepper.boundState(state.stepping, surface, false,
0495                                       freeToBoundCorrection);
0496         if (!res.ok()) {
0497           ACTS_DEBUG("Propagate to surface " << surface.geometryId()
0498                                              << " failed: " << res.error());
0499           return res.error();
0500         }
0501         const auto& [boundParams, jacobian, pathLength] = *res;
0502 
0503         // Fill the track state
0504         trackStateProxy.predicted() = boundParams.parameters();
0505         trackStateProxy.predictedCovariance() = state.stepping.cov;
0506 
0507         trackStateProxy.jacobian() = jacobian;
0508         trackStateProxy.pathLength() = pathLength;
0509 
0510         // We have predicted parameters, so calibrate the uncalibrated input
0511         // measurement
0512         extensions.calibrator(state.geoContext, *calibrationContext,
0513                               sourceLinkIt->second, trackStateProxy);
0514 
0515         // Get and set the type flags
0516         auto typeFlags = trackStateProxy.typeFlags();
0517         typeFlags.setHasParameters();
0518         if (surface.hasMaterial()) {
0519           typeFlags.setHasMaterial();
0520         }
0521 
0522         // Check if the state is an outlier.
0523         // If not:
0524         // - run Kalman update
0525         // - tag it as a measurement
0526         // - update the stepping state.
0527         // Else, just tag it as an outlier
0528         if (!extensions.outlierFinder(trackStateProxyConst)) {
0529           // Allocate the filtered parameters right before they are written
0530           trackStateProxy.addComponents(TrackStatePropMask::Filtered);
0531           // Run Kalman update
0532           auto updateRes =
0533               extensions.updater(state.geoContext, trackStateProxy, logger());
0534           if (!updateRes.ok()) {
0535             ACTS_DEBUG("Update step failed: " << updateRes.error());
0536             return updateRes.error();
0537           }
0538           // Set the measurement type flag
0539           typeFlags.setIsMeasurement();
0540         } else {
0541           ACTS_VERBOSE(
0542               "Filtering step successful. But measurement is determined "
0543               "to be an outlier. Stepping state is not updated.");
0544           // Set the outlier type flag
0545           typeFlags.setIsOutlier();
0546           trackStateProxy.shareFrom(trackStateProxy,
0547                                     TrackStatePropMask::Predicted,
0548                                     TrackStatePropMask::Filtered);
0549         }
0550 
0551         result.lastTrackIndex = trackStateProxy.index();
0552 
0553         // Update the stepper if it is not an outlier
0554         if (trackStateProxy.typeFlags().isMeasurement()) {
0555           // Update the stepping state with filtered parameters
0556           ACTS_VERBOSE("Filtering step successful, updated parameters are:\n"
0557                        << trackStateProxy.filtered().transpose());
0558           // update stepping state using filtered parameters after kalman
0559           stepper.update(state.stepping,
0560                          MultiTrajectoryHelpers::freeFiltered(
0561                              state.options.geoContext, trackStateProxy),
0562                          trackStateProxy.filtered(),
0563                          trackStateProxy.filteredCovariance(), surface);
0564           // We count the state with measurement
0565           ++result.measurementStates;
0566         }
0567 
0568         // Update state and stepper with post material effects
0569         const Result<detail::PointwiseMaterialEffects>
0570             materialInteractionPostRes = detail::performMaterialInteraction(
0571                 state, stepper, surface,
0572                 detail::determineMaterialUpdateMode(
0573                     state, navigator, MaterialUpdateMode::PostUpdate),
0574                 NoiseUpdateMode::addNoise, multipleScattering, energyLoss,
0575                 logger());
0576         if (!materialInteractionPostRes.ok()) {
0577           ACTS_DEBUG("Material interaction failed during filter: "
0578                      << materialInteractionPostRes.error().message());
0579           return materialInteractionPostRes.error();
0580         }
0581         // We count the processed state
0582         ++result.processedStates;
0583         // Update the number of holes count only when encountering a
0584         // measurement
0585         result.measurementHoles = result.missedActiveSurfaces.size();
0586         // Since we encountered a measurement update the lastMeasurementIndex to
0587         // the lastTrackIndex.
0588         result.lastMeasurementIndex = result.lastTrackIndex;
0589 
0590       } else if ((precedingMeasurementExists && surfaceIsSensitive) ||
0591                  surfaceHasMaterial) {
0592         // We only create track states here if there is already measurement
0593         // detected or if the surface has material (no holes before the first
0594         // measurement)
0595 
0596         // Create a track state with the desired components
0597         TrackStatePropMask mask =
0598             TrackStatePropMask::Predicted | TrackStatePropMask::Jacobian;
0599         typename traj_t::TrackStateProxy trackStateProxy =
0600             result.fittedStates->makeTrackState(mask, result.lastTrackIndex);
0601 
0602         // Set the trackStateProxy components with the state from the ongoing
0603         // propagation
0604         trackStateProxy.setReferenceSurface(surface.getSharedPtr());
0605         // Bind the transported state to the current surface
0606         auto res = stepper.boundState(state.stepping, surface, true,
0607                                       freeToBoundCorrection);
0608         if (!res.ok()) {
0609           return res.error();
0610         }
0611         const auto& [boundParams, jacobian, pathLength] = *res;
0612 
0613         // Fill the track state
0614         trackStateProxy.predicted() = boundParams.parameters();
0615         trackStateProxy.predictedCovariance() = state.stepping.cov;
0616 
0617         trackStateProxy.jacobian() = jacobian;
0618         trackStateProxy.pathLength() = pathLength;
0619 
0620         // Set the filtered parameter index to be the same with predicted
0621         // parameter
0622         trackStateProxy.shareFrom(trackStateProxy,
0623                                   TrackStatePropMask::Predicted,
0624                                   TrackStatePropMask::Filtered);
0625 
0626         // Set the track state flags
0627         auto typeFlags = trackStateProxy.typeFlags();
0628         typeFlags.setHasParameters();
0629 
0630         if (surfaceHasMaterial) {
0631           typeFlags.setHasMaterial();
0632         }
0633 
0634         if (surfaceIsSensitive && precedingMeasurementExists) {
0635           ACTS_VERBOSE("Detected hole on " << surface.geometryId());
0636           // If the surface is sensitive, set the hole type flag
0637           typeFlags.setIsHole();
0638         } else if (surfaceIsSensitive) {
0639           ACTS_VERBOSE("Skip hole (no preceding measurements) on surface "
0640                        << surface.geometryId());
0641         } else if (surfaceHasMaterial) {
0642           ACTS_VERBOSE("Detected in-sensitive surface "
0643                        << surface.geometryId());
0644         }
0645 
0646         result.lastTrackIndex = trackStateProxy.index();
0647 
0648         if (trackStateProxy.typeFlags().isHole()) {
0649           // Count the missed surface
0650           result.missedActiveSurfaces.push_back(&surface);
0651         }
0652 
0653         ++result.processedStates;
0654 
0655         // Update state and stepper with (possible) material effects
0656         const Result<detail::PointwiseMaterialEffects> materialInteractionRes =
0657             detail::performMaterialInteraction(
0658                 state, stepper, surface,
0659                 detail::determineMaterialUpdateMode(
0660                     state, navigator, MaterialUpdateMode::FullUpdate),
0661                 NoiseUpdateMode::addNoise, multipleScattering, energyLoss,
0662                 logger());
0663         if (!materialInteractionRes.ok()) {
0664           ACTS_DEBUG("Material interaction failed during filter: "
0665                      << materialInteractionRes.error().message());
0666           return materialInteractionRes.error();
0667         }
0668       }
0669 
0670       return Result<void>::success();
0671     }
0672   };
0673 
0674  public:
0675   /// Fit implementation of the forward filter, calls the
0676   /// the filter and smoother/reversed filter
0677   ///
0678   /// @tparam source_link_iterator_t Iterator type used to pass source links
0679   /// @tparam track_container_t Type of the track container
0680   ///
0681   /// @param it Begin iterator for the fittable uncalibrated measurements
0682   /// @param end End iterator for the fittable uncalibrated measurements
0683   /// @param sParameters The initial track parameters
0684   /// @param kfOptions KalmanOptions steering the fit
0685   /// @param trackContainer Input track container storage to append into
0686   /// @note The input measurements are given in the form of @c SourceLink s.
0687   /// It's the calibrators job to turn them into calibrated measurements used in
0688   /// the fit.
0689   ///
0690   /// @return the output as an output track
0691   template <typename source_link_iterator_t,
0692             TrackContainerFrontend track_container_t>
0693   Result<typename track_container_t::TrackProxy> fit(
0694       source_link_iterator_t it, source_link_iterator_t end,
0695       const BoundTrackParameters& sParameters,
0696       const KalmanFitterOptions<traj_t>& kfOptions,
0697       track_container_t& trackContainer) const {
0698     return fit_impl(it, end, sParameters, kfOptions, nullptr, trackContainer);
0699   }
0700 
0701   /// Fit implementation of the forward filter, calls the
0702   /// the filter and smoother/reversed filter
0703   ///
0704   /// @tparam source_link_iterator_t Iterator type used to pass source links
0705   /// @tparam track_container_t Type of the track container
0706   ///
0707   /// @param it Begin iterator for the fittable uncalibrated measurements
0708   /// @param end End iterator for the fittable uncalibrated measurements
0709   /// @param sParameters The initial track parameters
0710   /// @param kfOptions KalmanOptions steering the fit
0711   /// @param sSequence surface sequence used to initialize a DirectNavigator
0712   /// @param trackContainer Input track container storage to append into
0713   /// @note The input measurements are given in the form of @c SourceLinks.
0714   /// It's
0715   /// @c calibrator_t's job to turn them into calibrated measurements used in
0716   /// the fit.
0717   ///
0718   /// @return the output as an output track
0719   template <typename source_link_iterator_t,
0720             TrackContainerFrontend track_container_t>
0721   Result<typename track_container_t::TrackProxy> fit(
0722       source_link_iterator_t it, source_link_iterator_t end,
0723       const BoundTrackParameters& sParameters,
0724       const KalmanFitterOptions<traj_t>& kfOptions,
0725       const std::vector<const Surface*>& sSequence,
0726       track_container_t& trackContainer) const
0727     requires(isDirectNavigator)
0728   {
0729     return fit_impl(it, end, sParameters, kfOptions, &sSequence,
0730                     trackContainer);
0731   }
0732 
0733  private:
0734   template <typename source_link_iterator_t>
0735   auto make_propagator_options(source_link_iterator_t it,
0736                                source_link_iterator_t end,
0737                                const KalmanFitterOptions<traj_t>& kfOptions,
0738                                const std::vector<const Surface*>* sSequence,
0739                                const Surface* targetSurface,
0740                                bool reverseDirection) const {
0741     using KalmanActor = Actor;
0742     using Actors = ActorList<KalmanActor>;
0743     using PropagatorOptions = typename propagator_t::template Options<Actors>;
0744 
0745     const std::size_t nMeasurements = std::distance(it, end);
0746 
0747     // To be able to find measurements later, we put them into a map
0748     // We need to copy input SourceLinks anyway, so the map can own them.
0749     ACTS_VERBOSE("Preparing " << nMeasurements << " input measurements");
0750     std::unordered_map<const Surface*, SourceLink> inputMeasurements;
0751     for (; it != end; ++it) {
0752       SourceLink sl = *it;
0753       const Surface* surface = kfOptions.extensions.surfaceAccessor(sl);
0754       inputMeasurements.try_emplace(surface, std::move(sl));
0755     }
0756 
0757     // Create relevant options for the propagation options
0758     PropagatorOptions propagatorOptions(kfOptions.geoContext,
0759                                         kfOptions.magFieldContext);
0760 
0761     // Set the trivial propagator options
0762     propagatorOptions.setPlainOptions(kfOptions.propagatorPlainOptions);
0763 
0764     if (reverseDirection) {
0765       propagatorOptions.direction = propagatorOptions.direction.invert();
0766     }
0767 
0768     if constexpr (!isDirectNavigator) {
0769       // Add the measurement surface as external surface to navigator.
0770       // We will try to hit those surface by ignoring boundary checks.
0771       for (const auto& [surface, _] : inputMeasurements) {
0772         propagatorOptions.navigation.appendExternalSurface(*surface);
0773       }
0774     } else {
0775       assert(sSequence != nullptr &&
0776              "DirectNavigator requires a surface sequence for KalmanFitter");
0777       // Set the surface sequence
0778       propagatorOptions.navigation.externalSurfaces = *sSequence;
0779     }
0780 
0781     // Catch the actor and set the measurements
0782     auto& kalmanActor = propagatorOptions.actorList.template get<KalmanActor>();
0783     kalmanActor.inputMeasurements = std::move(inputMeasurements);
0784     kalmanActor.targetReached.surface = targetSurface;
0785     kalmanActor.multipleScattering = kfOptions.multipleScattering;
0786     kalmanActor.energyLoss = kfOptions.energyLoss;
0787     kalmanActor.freeToBoundCorrection = kfOptions.freeToBoundCorrection;
0788     kalmanActor.calibrationContext = &kfOptions.calibrationContext.get();
0789     kalmanActor.extensions = kfOptions.extensions;
0790     kalmanActor.actorLogger = m_actorLogger.get();
0791 
0792     return propagatorOptions;
0793   }
0794 
0795   template <typename propagator_options_t,
0796             TrackContainerFrontend track_container_t>
0797   auto filter_impl(const BoundTrackParameters& sParameters,
0798                    const propagator_options_t& propagatorOptions,
0799                    track_container_t& trackContainer) const
0800       -> Result<typename track_container_t::TrackProxy> {
0801     auto propagatorState = m_propagator.makeState(propagatorOptions);
0802 
0803     auto propagatorInitResult =
0804         m_propagator.initialize(propagatorState, sParameters);
0805     if (!propagatorInitResult.ok()) {
0806       ACTS_DEBUG("Propagation initialization failed: "
0807                  << propagatorInitResult.error());
0808       return propagatorInitResult.error();
0809     }
0810 
0811     auto& kalmanResult =
0812         propagatorState.template get<KalmanFitterResult<traj_t>>();
0813     kalmanResult.fittedStates = &trackContainer.trackStateContainer();
0814 
0815     // Run the fitter
0816     auto result = m_propagator.propagate(propagatorState);
0817 
0818     if (!result.ok()) {
0819       ACTS_DEBUG("Propagation failed: " << result.error());
0820       return result.error();
0821     }
0822 
0823     /// It could happen that the fit ends in zero measurement states.
0824     /// The result gets meaningless so such case is regarded as fit failure.
0825     if (!kalmanResult.measurementStates) {
0826       ACTS_DEBUG("KalmanFilter failed: No measurement states found");
0827       return KalmanFitterError::NoMeasurementFound;
0828     }
0829 
0830     auto track = trackContainer.makeTrack();
0831     track.tipIndex() = kalmanResult.lastMeasurementIndex;
0832     if (kalmanResult.fittedParameters) {
0833       const auto& params = kalmanResult.fittedParameters.value();
0834       track.parameters() = params.parameters();
0835       track.covariance() = params.covariance().value();
0836       track.setReferenceSurface(params.referenceSurface().getSharedPtr());
0837     }
0838 
0839     calculateTrackQuantities(track);
0840 
0841     track.linkForward();
0842 
0843     return track;
0844   }
0845 
0846   /// Common fit implementation
0847   ///
0848   /// @tparam source_link_iterator_t Iterator type used to pass source links
0849   /// @tparam track_container_t Type of the track container
0850   ///
0851   /// @param it Begin iterator for the fittable uncalibrated measurements
0852   /// @param end End iterator for the fittable uncalibrated measurements
0853   /// @param sParameters The initial track parameters
0854   /// @param kfOptions KalmanOptions steering the fit
0855   /// @param sSequence surface sequence used to initialize a DirectNavigator
0856   /// @param trackContainer Input track container storage to append into
0857   ///
0858   /// @return the output as an output track
0859   template <typename source_link_iterator_t,
0860             TrackContainerFrontend track_container_t>
0861   auto fit_impl(source_link_iterator_t it, source_link_iterator_t end,
0862                 const BoundTrackParameters& sParameters,
0863                 const KalmanFitterOptions<traj_t>& kfOptions,
0864                 const std::vector<const Surface*>* sSequence,
0865                 track_container_t& trackContainer) const
0866       -> Result<typename track_container_t::TrackProxy> {
0867     using TrackProxy = typename track_container_t::TrackProxy;
0868     using TrackStateProxy = typename track_container_t::TrackStateProxy;
0869 
0870     auto forwardPropagatorOptions =
0871         make_propagator_options(it, end, kfOptions, sSequence, nullptr, false);
0872 
0873     auto forwardFilterResult =
0874         filter_impl(sParameters, forwardPropagatorOptions, trackContainer);
0875 
0876     if (!forwardFilterResult.ok()) {
0877       ACTS_DEBUG("KalmanFilter failed: "
0878                  << forwardFilterResult.error() << ", "
0879                  << forwardFilterResult.error().message());
0880       return forwardFilterResult.error();
0881     }
0882 
0883     TrackProxy forwardTrack = forwardFilterResult.value();
0884 
0885     TrackStateProxy firstMeasurementState =
0886         trackContainer.trackStateContainer().getTrackState(
0887             findFirstMeasurementState(forwardTrack).value().index());
0888     TrackStateProxy lastMeasurementState =
0889         trackContainer.trackStateContainer().getTrackState(
0890             findLastMeasurementState(forwardTrack).value().index());
0891     lastMeasurementState.shareFrom(lastMeasurementState,
0892                                    TrackStatePropMask::Filtered,
0893                                    TrackStatePropMask::Smoothed);
0894 
0895     TrackProxy track = forwardTrack;
0896 
0897     const bool doReverseFilter =
0898         kfOptions.reverseFiltering ||
0899         kfOptions.extensions.reverseFilteringLogic(
0900             typename traj_t::ConstTrackStateProxy(lastMeasurementState));
0901     if (doReverseFilter) {
0902       ACTS_VERBOSE("Smooth track by reversed filtering");
0903 
0904       auto reverseStartParameters = forwardTrack.createParametersFromState(
0905           typename traj_t::ConstTrackStateProxy(lastMeasurementState));
0906       reverseStartParameters.covariance().value() *=
0907           kfOptions.reverseFilteringCovarianceScaling;
0908       auto reversePropagatorOptions = make_propagator_options(
0909           it, end, kfOptions, sSequence, kfOptions.referenceSurface, true);
0910       auto reverseFilterResult = filter_impl(
0911           reverseStartParameters, reversePropagatorOptions, trackContainer);
0912 
0913       if (!reverseFilterResult.ok()) {
0914         ACTS_DEBUG("Reversed KalmanFilter failed: "
0915                    << reverseFilterResult.error() << ", "
0916                    << reverseFilterResult.error().message());
0917         return reverseFilterResult.error();
0918       }
0919 
0920       TrackProxy reverseTrack = reverseFilterResult.value();
0921 
0922       TrackStateProxy reverseLastMeasurementState =
0923           trackContainer.trackStateContainer().getTrackState(
0924               findLastMeasurementState(reverseTrack).value().index());
0925 
0926       if (&firstMeasurementState.referenceSurface() !=
0927           &reverseLastMeasurementState.referenceSurface()) {
0928         ACTS_DEBUG(
0929             "Inconsistent reference surfaces between forward and "
0930             "reversed filtered tracks");
0931         return Result<TrackProxy>::failure(
0932             KalmanFitterError::InconsistentTrackStates);
0933       }
0934       firstMeasurementState.shareFrom(reverseLastMeasurementState,
0935                                       TrackStatePropMask::Filtered,
0936                                       TrackStatePropMask::Smoothed);
0937 
0938       if (reverseTrack.hasReferenceSurface()) {
0939         track.parameters() = reverseTrack.parameters();
0940         track.covariance() = reverseTrack.covariance();
0941         track.setReferenceSurface(
0942             reverseTrack.referenceSurface().getSharedPtr());
0943       }
0944 
0945       trackContainer.removeTrack(reverseTrack.index());
0946     } else {
0947       ACTS_VERBOSE("Smooth track directly without reversed filtering");
0948 
0949       auto smoothRes = kfOptions.extensions.smoother(
0950           kfOptions.geoContext, trackContainer.trackStateContainer(),
0951           forwardTrack.tipIndex(), logger());
0952       if (!smoothRes.ok()) {
0953         ACTS_DEBUG("Smoothing step failed: " << smoothRes.error() << ", "
0954                                              << smoothRes.error().message());
0955         return smoothRes.error();
0956       }
0957     }
0958 
0959     if (!track.hasReferenceSurface() && kfOptions.referenceSurface != nullptr) {
0960       typename propagator_t::template Options<> extrapolationOptions(
0961           kfOptions.geoContext, kfOptions.magFieldContext);
0962       auto extrapolationResult = extrapolateTrackToReferenceSurface(
0963           track, *kfOptions.referenceSurface, m_propagator,
0964           extrapolationOptions, kfOptions.referenceSurfaceStrategy, logger());
0965 
0966       if (!extrapolationResult.ok()) {
0967         ACTS_DEBUG("Extrapolation to reference surface failed: "
0968                    << extrapolationResult.error() << ", "
0969                    << extrapolationResult.error().message());
0970         return extrapolationResult.error();
0971       }
0972     }
0973 
0974     if (trackContainer.hasColumn(hashString("smoothed"))) {
0975       track.template component<bool, hashString("smoothed")>() =
0976           !doReverseFilter;
0977     }
0978     if (trackContainer.hasColumn(hashString("reversed"))) {
0979       track.template component<bool, hashString("reversed")>() =
0980           doReverseFilter;
0981     }
0982 
0983     return track;
0984   }
0985 };
0986 
0987 /// @}
0988 
0989 }  // namespace Acts