Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2026-05-22 08:00: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/TrackParameters.hpp"
0015 #include "Acts/EventData/TrackStatePropMask.hpp"
0016 #include "Acts/EventData/Types.hpp"
0017 #include "Acts/Geometry/GeometryContext.hpp"
0018 #include "Acts/MagneticField/MagneticFieldContext.hpp"
0019 #include "Acts/Propagator/ActorList.hpp"
0020 #include "Acts/Propagator/ConstrainedStep.hpp"
0021 #include "Acts/Propagator/PropagatorState.hpp"
0022 #include "Acts/Propagator/StandardAborters.hpp"
0023 #include "Acts/Propagator/detail/LoopProtection.hpp"
0024 #include "Acts/Propagator/detail/PointwiseMaterialInteraction.hpp"
0025 #include "Acts/TrackFinding/CombinatorialKalmanFilterError.hpp"
0026 #include "Acts/TrackFinding/CombinatorialKalmanFilterExtensions.hpp"
0027 #include "Acts/Utilities/CalibrationContext.hpp"
0028 #include "Acts/Utilities/Logger.hpp"
0029 #include "Acts/Utilities/Result.hpp"
0030 
0031 #include <functional>
0032 #include <limits>
0033 #include <memory>
0034 #include <type_traits>
0035 
0036 namespace Acts {
0037 
0038 /// @addtogroup track_finding
0039 /// @{
0040 
0041 /// Combined options for the combinatorial Kalman filter.
0042 ///
0043 /// @tparam source_link_iterator_t Type of the source link iterator
0044 /// @tparam track_container_t Type of the track container
0045 template <typename track_container_t>
0046 struct CombinatorialKalmanFilterOptions {
0047   /// Type alias for track state container backend
0048   using TrackStateContainerBackend =
0049       typename track_container_t::TrackStateContainerBackend;
0050   /// Type alias for track state proxy from the container
0051   using TrackStateProxy = typename track_container_t::TrackStateProxy;
0052 
0053   /// PropagatorOptions with context
0054   ///
0055   /// @param gctx The geometry context for this track finding/fitting
0056   /// @param mctx The magnetic context for this track finding/fitting
0057   /// @param cctx The calibration context for this track finding/fitting
0058   /// @param extensions_ The extension struct
0059   /// @param pOptions The plain propagator options
0060   /// @param mScattering Whether to include multiple scattering
0061   /// @param eLoss Whether to include energy loss
0062   CombinatorialKalmanFilterOptions(
0063       const GeometryContext& gctx, const MagneticFieldContext& mctx,
0064       std::reference_wrapper<const CalibrationContext> cctx,
0065       CombinatorialKalmanFilterExtensions<track_container_t> extensions_,
0066       const PropagatorPlainOptions& pOptions, bool mScattering = true,
0067       bool eLoss = true)
0068       : geoContext(gctx),
0069         magFieldContext(mctx),
0070         calibrationContext(cctx),
0071         extensions(extensions_),
0072         propagatorPlainOptions(pOptions),
0073         multipleScattering(mScattering),
0074         energyLoss(eLoss) {}
0075 
0076   /// Contexts are required and the options must not be default-constructible.
0077   CombinatorialKalmanFilterOptions() = delete;
0078 
0079   /// Context object for the geometry
0080   std::reference_wrapper<const GeometryContext> geoContext;
0081   /// Context object for the magnetic field
0082   std::reference_wrapper<const MagneticFieldContext> magFieldContext;
0083   /// context object for the calibration
0084   std::reference_wrapper<const CalibrationContext> calibrationContext;
0085 
0086   /// The filter extensions
0087   CombinatorialKalmanFilterExtensions<track_container_t> extensions;
0088 
0089   /// The trivial propagator options
0090   PropagatorPlainOptions propagatorPlainOptions;
0091 
0092   /// The target surface
0093   /// @note This is useful if the filtering should be terminated at a
0094   ///       certain surface
0095   const Surface* targetSurface = nullptr;
0096 
0097   /// Whether to consider multiple scattering.
0098   bool multipleScattering = true;
0099 
0100   /// Whether to consider energy loss.
0101   bool energyLoss = true;
0102 
0103   /// Skip the pre propagation call. This effectively skips the first surface
0104   /// @note This is useful if the first surface should not be considered in a second reverse pass
0105   bool skipPrePropagationUpdate = false;
0106 };
0107 
0108 /// Result container for the combinatorial Kalman filter actor.
0109 ///
0110 /// @tparam track_container_t Type of the track container storing results
0111 template <typename track_container_t>
0112 struct CombinatorialKalmanFilterResult {
0113   /// Track state container backend type
0114   using TrackStateContainerBackend =
0115       typename track_container_t::TrackStateContainerBackend;
0116   /// Track proxy type
0117   using TrackProxy = typename track_container_t::TrackProxy;
0118   /// Track state proxy type
0119   using TrackStateProxy = typename track_container_t::TrackStateProxy;
0120 
0121   /// The track container to store the found tracks
0122   track_container_t* tracks{nullptr};
0123 
0124   /// Fitted states that the actor has handled.
0125   TrackStateContainerBackend* trackStates{nullptr};
0126 
0127   /// Indices into `tracks` which mark active branches
0128   std::vector<TrackProxy> activeBranches;
0129 
0130   /// Indices into `tracks` which mark active branches
0131   std::vector<TrackProxy> collectedTracks;
0132 
0133   /// Track state candidates buffer which can be used by the track state creator
0134   std::vector<TrackStateProxy> trackStateCandidates;
0135 
0136   /// Indicator if track finding has been done
0137   bool finished = false;
0138 
0139   /// Path limit aborter
0140   PathLimitReached pathLimitReached;
0141 };
0142 
0143 /// Combinatorial Kalman filter to find tracks.
0144 ///
0145 /// @tparam propagator_t Type of the propagator
0146 ///
0147 /// The CombinatorialKalmanFilter contains an Actor and a Sequencer sub-class.
0148 /// The Sequencer has to be part of the Navigator of the Propagator in order to
0149 /// initialize and provide the measurement surfaces.
0150 ///
0151 /// The Actor is part of the Propagation call and does the Kalman update.
0152 /// Updater and Calibrator are given to the Actor for further use:
0153 /// - The Updater is the implemented kalman updater formalism, it
0154 ///   runs via a visitor pattern through the measurements.
0155 ///
0156 /// Measurements are not required to be ordered for the
0157 /// CombinatorialKalmanFilter, measurement ordering needs to be figured out by
0158 /// the navigation of the propagator.
0159 ///
0160 /// The void components are provided mainly for unit testing.
0161 ///
0162 template <typename propagator_t, typename track_container_t>
0163 class CombinatorialKalmanFilter {
0164  public:
0165   /// Default constructor is deleted
0166   CombinatorialKalmanFilter() = delete;
0167 
0168   /// Constructor with propagator and logging level
0169   /// @param pPropagator The propagator used for the track finding
0170   /// @param _logger The logger for messages
0171   explicit CombinatorialKalmanFilter(propagator_t pPropagator,
0172                                      std::unique_ptr<const Logger> _logger =
0173                                          getDefaultLogger("CKF", Logging::INFO))
0174       : m_propagator(std::move(pPropagator)),
0175         m_logger(std::move(_logger)),
0176         m_actorLogger{m_logger->cloneWithSuffix("Actor")},
0177         m_updaterLogger{m_logger->cloneWithSuffix("Updater")} {}
0178 
0179  private:
0180   using BoundState = std::tuple<BoundTrackParameters, BoundMatrix, double>;
0181   using TrackStateContainerBackend =
0182       typename track_container_t::TrackStateContainerBackend;
0183   using TrackProxy = typename track_container_t::TrackProxy;
0184   using TrackStateProxy = typename track_container_t::TrackStateProxy;
0185 
0186   /// The propagator for the transport and material update
0187   propagator_t m_propagator;
0188 
0189   std::unique_ptr<const Logger> m_logger;
0190   std::shared_ptr<const Logger> m_actorLogger;
0191   std::shared_ptr<const Logger> m_updaterLogger;
0192 
0193   const Logger& logger() const { return *m_logger; }
0194 
0195   /// @brief Propagator Actor plugin for the CombinatorialKalmanFilter
0196   ///
0197   /// The CombinatorialKalmanFilter Actor does not rely on the measurements to
0198   /// be sorted along the track.
0199   class Actor {
0200    public:
0201     using BoundState = std::tuple<BoundTrackParameters, BoundMatrix, double>;
0202     /// Broadcast the result_type
0203     using result_type = CombinatorialKalmanFilterResult<track_container_t>;
0204 
0205     using BranchStopperResult = CombinatorialKalmanFilterBranchStopperResult;
0206 
0207     /// The target surface aborter
0208     SurfaceReached targetReached{std::numeric_limits<double>::lowest()};
0209 
0210     /// Whether to consider multiple scattering.
0211     bool multipleScattering = true;
0212 
0213     /// Whether to consider energy loss.
0214     bool energyLoss = true;
0215 
0216     /// Skip the pre propagation call. This effectively skips the first surface
0217     bool skipPrePropagationUpdate = false;
0218 
0219     /// Calibration context for the finding run
0220     const CalibrationContext* calibrationContextPtr{nullptr};
0221 
0222     CombinatorialKalmanFilterExtensions<track_container_t> extensions;
0223 
0224     /// End of world aborter
0225     EndOfWorldReached endOfWorldReached;
0226 
0227     /// Volume constraint aborter
0228     VolumeConstraintAborter volumeConstraintAborter;
0229 
0230     /// Actor logger instance
0231     const Logger* actorLogger{nullptr};
0232     /// Updater logger instance
0233     const Logger* updaterLogger{nullptr};
0234 
0235     const Logger& logger() const { return *actorLogger; }
0236 
0237     /// @brief CombinatorialKalmanFilter actor operation
0238     ///
0239     /// @tparam propagator_state_t Type of the Propagator state
0240     /// @tparam stepper_t Type of the stepper
0241     ///
0242     /// @param state is the mutable propagator state object
0243     /// @param stepper is the stepper in use
0244     /// @param navigator is the navigator in use
0245     /// @param result is the mutable result state object
0246     template <typename propagator_state_t, typename stepper_t,
0247               typename navigator_t>
0248     Result<void> act(propagator_state_t& state, const stepper_t& stepper,
0249                      const navigator_t& navigator, result_type& result,
0250                      const Logger& /*logger*/) const {
0251       ACTS_VERBOSE("CKF Actor called");
0252 
0253       assert(result.trackStates && "No MultiTrajectory set");
0254 
0255       if (state.stage == PropagatorStage::prePropagation &&
0256           skipPrePropagationUpdate) {
0257         ACTS_VERBOSE("Skip pre-propagation update (first surface)");
0258         return Result<void>::success();
0259       }
0260       if (state.stage == PropagatorStage::postPropagation) {
0261         ACTS_VERBOSE("Skip post-propagation action");
0262         return Result<void>::success();
0263       }
0264 
0265       ACTS_VERBOSE("CombinatorialKalmanFilter step");
0266 
0267       assert(!result.activeBranches.empty() && "No active branches");
0268       assert(!result.finished && "Should never reach this when finished");
0269 
0270       // Initialize path limit reached aborter
0271       if (result.pathLimitReached.internalLimit ==
0272           std::numeric_limits<double>::max()) {
0273         detail::setupLoopProtection(state, stepper, result.pathLimitReached,
0274                                     true, logger());
0275       }
0276 
0277       // Update:
0278       // - Waiting for a current surface
0279       if (const Surface* surface = navigator.currentSurface(state.navigation);
0280           surface != nullptr) {
0281         // There are three scenarios:
0282         // 1) The surface is in the measurement map
0283         // -> Select source links
0284         // -> Perform the kalman update for selected non-outlier source links
0285         // -> Add track states in multitrajectory. Multiple states mean branch
0286         // splitting.
0287         // -> Call branch stopper to justify each branch
0288         // -> If there is non-outlier state, update stepper information
0289         // 2) The surface is not in the measurement map but with material or is
0290         // an active surface
0291         // -> Add a hole or passive material state in multitrajectory
0292         // -> Call branch stopper to justify the branch
0293         // 3) The surface is neither in the measurement map nor with material
0294         // -> Do nothing
0295         ACTS_VERBOSE("Perform filter step");
0296         auto res = filter(*surface, state, stepper, navigator, result);
0297         if (!res.ok()) {
0298           ACTS_DEBUG("Error in filter: " << res.error().message());
0299           return res.error();
0300         }
0301 
0302         if (result.finished) {
0303           ACTS_VERBOSE("CKF Actor returns after filter step");
0304           return Result<void>::success();
0305         }
0306       }
0307 
0308       assert(!result.activeBranches.empty() && "No active branches");
0309 
0310       const bool isEndOfWorldReached =
0311           endOfWorldReached.checkAbort(state, stepper, navigator, logger());
0312       const bool isVolumeConstraintReached = volumeConstraintAborter.checkAbort(
0313           state, stepper, navigator, logger());
0314       const bool isPathLimitReached = result.pathLimitReached.checkAbort(
0315           state, stepper, navigator, logger());
0316       const bool isTargetReached =
0317           targetReached.checkAbort(state, stepper, navigator, logger());
0318       if (isEndOfWorldReached || isVolumeConstraintReached ||
0319           isPathLimitReached || isTargetReached) {
0320         if (isEndOfWorldReached) {
0321           ACTS_VERBOSE("End of world reached");
0322         } else if (isVolumeConstraintReached) {
0323           ACTS_VERBOSE("Volume constraint reached");
0324         } else if (isPathLimitReached) {
0325           ACTS_VERBOSE("Path limit reached");
0326         } else if (isTargetReached) {
0327           ACTS_VERBOSE("Target surface reached");
0328 
0329           // Bind the parameter to the target surface
0330           auto res = stepper.boundState(state.stepping, *targetReached.surface);
0331           if (!res.ok()) {
0332             ACTS_DEBUG("Error while acquiring bound state for target surface: "
0333                        << res.error() << " " << res.error().message());
0334             return res.error();
0335           }
0336 
0337           const auto& [boundParams, jacobian, pathLength] = *res;
0338           auto currentBranch = result.activeBranches.back();
0339           // Assign the fitted parameters
0340           currentBranch.parameters() = boundParams.parameters();
0341           currentBranch.covariance() = *boundParams.covariance();
0342           currentBranch.setReferenceSurface(
0343               boundParams.referenceSurface().getSharedPtr());
0344 
0345           stepper.releaseStepSize(state.stepping,
0346                                   ConstrainedStep::Type::Navigator);
0347         }
0348 
0349         // Record the active branch and remove it from the list
0350         storeLastActiveBranch(result);
0351         result.activeBranches.pop_back();
0352 
0353         // Reset propagation state to track state at next active branch
0354         auto resetRes = reset(state, stepper, navigator, result);
0355         if (!resetRes.ok()) {
0356           return resetRes.error();
0357         }
0358       }
0359 
0360       return Result<void>::success();
0361     }
0362 
0363     template <typename propagator_state_t, typename stepper_t,
0364               typename navigator_t>
0365     bool checkAbort(propagator_state_t& /*state*/, const stepper_t& /*stepper*/,
0366                     const navigator_t& /*navigator*/, const result_type& result,
0367                     const Logger& /*logger*/) const {
0368       return result.finished;
0369     }
0370 
0371     /// @brief CombinatorialKalmanFilter actor operation: reset propagation
0372     ///
0373     /// @tparam propagator_state_t Type of Propagator state
0374     /// @tparam stepper_t Type of the stepper
0375     /// @tparam navigator_t Type of the navigator
0376     ///
0377     /// @param state is the mutable propagator state object
0378     /// @param stepper is the stepper in use
0379     /// @param navigator is the navigator in use
0380     /// @param result is the mutable result state object
0381     template <typename propagator_state_t, typename stepper_t,
0382               typename navigator_t>
0383     Result<void> reset(propagator_state_t& state, const stepper_t& stepper,
0384                        const navigator_t& navigator,
0385                        result_type& result) const {
0386       if (result.activeBranches.empty()) {
0387         ACTS_VERBOSE("Stop CKF with " << result.collectedTracks.size()
0388                                       << " found tracks");
0389         result.finished = true;
0390 
0391         return Result<void>::success();
0392       }
0393 
0394       auto currentBranch = result.activeBranches.back();
0395       auto currentState = currentBranch.outermostTrackState();
0396 
0397       ACTS_VERBOSE("Propagation jumps to branch with tip = "
0398                    << currentBranch.tipIndex());
0399 
0400       // Reset the stepping state
0401       stepper.initialize(state.stepping, currentState.filtered(),
0402                          currentState.filteredCovariance(),
0403                          stepper.particleHypothesis(state.stepping),
0404                          currentState.referenceSurface());
0405 
0406       // Reset the navigation state
0407       // Set targetSurface to nullptr for forward filtering
0408       state.navigation.options.startSurface = &currentState.referenceSurface();
0409       state.navigation.options.targetSurface = nullptr;
0410       auto navInitRes = navigator.initialize(
0411           state.navigation, stepper.position(state.stepping),
0412           stepper.direction(state.stepping), state.options.direction);
0413       if (!navInitRes.ok()) {
0414         ACTS_DEBUG("Navigation initialization failed: " << navInitRes.error());
0415         return navInitRes.error();
0416       }
0417 
0418       // No Kalman filtering for the starting surface, but still need
0419       // to consider the material effects here
0420       detail::performMaterialInteraction(
0421           state, stepper, navigator,
0422           detail::determineMaterialUpdateMode(state, navigator,
0423                                               MaterialUpdateMode::PostUpdate),
0424           NoiseUpdateMode::addNoise, multipleScattering, energyLoss, logger());
0425 
0426       // Set path limit based on loop protection
0427       detail::setupLoopProtection(state, stepper, result.pathLimitReached, true,
0428                                   logger());
0429 
0430       // Set path limit based on target surface
0431       targetReached.checkAbort(state, stepper, navigator, logger());
0432 
0433       return Result<void>::success();
0434     }
0435 
0436     /// @brief CombinatorialKalmanFilter actor operation:
0437     /// - filtering for all measurement(s) on surface
0438     /// - store selected track states in multiTrajectory
0439     /// - update propagator state to the (last) selected track state
0440     ///
0441     /// @tparam propagator_state_t Type of the 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       using PM = TrackStatePropMask;
0456 
0457       bool isSensitive = surface.isSensitive();
0458       bool hasMaterial = surface.surfaceMaterial() != nullptr;
0459       bool isMaterialOnly = hasMaterial && !isSensitive;
0460       bool expectMeasurements = isSensitive;
0461 
0462       if (isSensitive) {
0463         ACTS_VERBOSE("Measurement surface " << surface.geometryId()
0464                                             << " detected.");
0465       } else if (isMaterialOnly) {
0466         ACTS_VERBOSE("Material surface " << surface.geometryId()
0467                                          << " detected.");
0468       } else {
0469         ACTS_VERBOSE("Passive surface " << surface.geometryId()
0470                                         << " detected.");
0471         return Result<void>::success();
0472       }
0473 
0474       // Transport the covariance to the surface
0475       if (isMaterialOnly) {
0476         stepper.transportCovarianceToCurvilinear(state.stepping);
0477       } else {
0478         stepper.transportCovarianceToBound(state.stepping, surface);
0479       }
0480 
0481       // Update state and stepper with pre material effects
0482       detail::performMaterialInteraction(
0483           state, stepper, navigator,
0484           detail::determineMaterialUpdateMode(state, navigator,
0485                                               MaterialUpdateMode::PreUpdate),
0486           NoiseUpdateMode::addNoise, multipleScattering, energyLoss, logger());
0487 
0488       // Bind the transported state to the current surface
0489       auto boundStateRes = stepper.boundState(state.stepping, surface, false);
0490       if (!boundStateRes.ok()) {
0491         return boundStateRes.error();
0492       }
0493       auto& boundState = *boundStateRes;
0494       auto& [boundParams, jacobian, pathLength] = boundState;
0495       boundParams.covariance() = state.stepping.cov;
0496 
0497       auto currentBranch = result.activeBranches.back();
0498       TrackIndexType prevTip = currentBranch.tipIndex();
0499 
0500       using TrackStatesResult = Result<CkfTypes::BranchVector<TrackIndexType>>;
0501       TrackStatesResult tsRes = TrackStatesResult::success({});
0502       if (isSensitive) {
0503         // extend trajectory with measurements associated to the current surface
0504         // which may create extra trajectory branches if more than one
0505         // measurement is selected.
0506         tsRes = extensions.createTrackStates(
0507             state.geoContext, *calibrationContextPtr, surface, boundState,
0508             prevTip, result.trackStateCandidates, *result.trackStates,
0509             logger());
0510       }
0511 
0512       if (tsRes.ok() && !(*tsRes).empty()) {
0513         const CkfTypes::BranchVector<TrackIndexType>& newTrackStateList =
0514             *tsRes;
0515         Result<unsigned int> procRes =
0516             processNewTrackStates(state.geoContext, newTrackStateList, result);
0517         if (!procRes.ok()) {
0518           ACTS_DEBUG("Processing of selected track states failed: "
0519                      << procRes.error().message());
0520           return procRes.error();
0521         }
0522         unsigned int nBranchesOnSurface = *procRes;
0523 
0524         if (nBranchesOnSurface == 0) {
0525           ACTS_VERBOSE("All branches on surface " << surface.geometryId()
0526                                                   << " have been stopped");
0527 
0528           reset(state, stepper, navigator, result);
0529 
0530           return Result<void>::success();
0531         }
0532 
0533         // `currentBranch` is invalidated after `processNewTrackStates`
0534         currentBranch = result.activeBranches.back();
0535         prevTip = currentBranch.tipIndex();
0536       } else {
0537         if (!tsRes.ok()) {
0538           if (static_cast<CombinatorialKalmanFilterError>(
0539                   tsRes.error().value()) ==
0540               CombinatorialKalmanFilterError::NoMeasurementExpected) {
0541             // recoverable error returned by track state creator
0542             expectMeasurements = false;
0543           } else {
0544             ACTS_DEBUG("Track state creation failed on surface "
0545                        << surface.geometryId() << ": " << tsRes.error());
0546             return tsRes.error();
0547           }
0548         }
0549 
0550         if (expectMeasurements) {
0551           ACTS_VERBOSE("Detected hole after measurement selection on surface "
0552                        << surface.geometryId());
0553         }
0554 
0555         auto stateMask = PM::Predicted | PM::Jacobian;
0556 
0557         // Add a hole or material track state to the multitrajectory
0558         TrackIndexType currentTip =
0559             addNonSourcelinkState(stateMask, boundState, result, isSensitive,
0560                                   expectMeasurements, prevTip);
0561         currentBranch.tipIndex() = currentTip;
0562         auto currentState = currentBranch.outermostTrackState();
0563         if (expectMeasurements) {
0564           currentBranch.nHoles()++;
0565         }
0566 
0567         BranchStopperResult branchStopperResult =
0568             extensions.branchStopper(currentBranch, currentState);
0569 
0570         // Check the branch
0571         if (branchStopperResult == BranchStopperResult::Continue) {
0572           // Remembered the active branch and its state
0573         } else {
0574           // No branch on this surface
0575           if (branchStopperResult == BranchStopperResult::StopAndKeep) {
0576             storeLastActiveBranch(result);
0577           }
0578           // Remove the branch from list
0579           result.activeBranches.pop_back();
0580 
0581           // Branch on the surface has been stopped - reset
0582           ACTS_VERBOSE("Branch on surface " << surface.geometryId()
0583                                             << " has been stopped");
0584 
0585           reset(state, stepper, navigator, result);
0586 
0587           return Result<void>::success();
0588         }
0589       }
0590 
0591       auto currentState = currentBranch.outermostTrackState();
0592 
0593       if (currentState.typeFlags().isOutlier()) {
0594         // We don't need to update the stepper given an outlier state
0595         ACTS_VERBOSE("Outlier state detected on surface "
0596                      << surface.geometryId());
0597       } else if (currentState.typeFlags().isMeasurement()) {
0598         // If there are measurement track states on this surface
0599         // Update stepping state using filtered parameters of last track
0600         // state on this surface
0601         stepper.update(state.stepping,
0602                        MultiTrajectoryHelpers::freeFiltered(
0603                            state.options.geoContext, currentState),
0604                        currentState.filtered(),
0605                        currentState.filteredCovariance(), surface);
0606         ACTS_VERBOSE("Stepping state is updated with filtered parameter:");
0607         ACTS_VERBOSE("-> " << currentState.filtered().transpose()
0608                            << " of track state with tip = "
0609                            << currentState.index());
0610       }
0611 
0612       // Update state and stepper with post material effects
0613       detail::performMaterialInteraction(
0614           state, stepper, navigator,
0615           detail::determineMaterialUpdateMode(state, navigator,
0616                                               MaterialUpdateMode::PostUpdate),
0617           NoiseUpdateMode::addNoise, multipleScattering, energyLoss, logger());
0618 
0619       return Result<void>::success();
0620     }
0621 
0622     /// Process new, incompomplete track states and set the filtered state
0623     ///
0624     /// @note will process the given list of new states, run the updater
0625     ///     or share the predicted state for states flagged as outliers
0626     ///     and add them to the list of active branches
0627     ///
0628     /// @param gctx The geometry context for this track finding/fitting
0629     /// @param newTrackStateList index list of new track states
0630     /// @param result which contains among others the new states, and the list of active branches
0631     /// @return the number of newly added branches or an error
0632     Result<unsigned int> processNewTrackStates(
0633         const GeometryContext& gctx,
0634         const CkfTypes::BranchVector<TrackIndexType>& newTrackStateList,
0635         result_type& result) const {
0636       using PM = TrackStatePropMask;
0637 
0638       unsigned int nBranchesOnSurface = 0;
0639 
0640       auto rootBranch = result.activeBranches.back();
0641 
0642       // Build the new branches by forking the root branch. Reverse the order
0643       // to process the best candidate first
0644       CkfTypes::BranchVector<TrackProxy> newBranches;
0645       for (auto it = newTrackStateList.rbegin(); it != newTrackStateList.rend();
0646            ++it) {
0647         // Keep the root branch as the first branch, make a copy for the
0648         // others
0649         auto shallowCopy = [&] {
0650           auto sc = rootBranch.container().makeTrack();
0651           sc.copyFromShallow(rootBranch);
0652           return sc;
0653         };
0654         auto newBranch =
0655             (it == newTrackStateList.rbegin()) ? rootBranch : shallowCopy();
0656         newBranch.tipIndex() = *it;
0657         newBranches.push_back(newBranch);
0658       }
0659 
0660       // Remove the root branch
0661       result.activeBranches.pop_back();
0662 
0663       // Update and select from the new branches
0664       for (TrackProxy newBranch : newBranches) {
0665         auto trackState = newBranch.outermostTrackState();
0666         TrackStateTypeMap typeFlags = trackState.typeFlags();
0667 
0668         if (typeFlags.isOutlier()) {
0669           // No Kalman update for outlier
0670           // Set the filtered parameter index to be the same with predicted
0671           // parameter
0672           trackState.shareFrom(PM::Predicted, PM::Filtered);
0673           // Increment number of outliers
0674           newBranch.nOutliers()++;
0675         } else if (typeFlags.isMeasurement()) {
0676           // Kalman update
0677           auto updateRes = extensions.updater(gctx, trackState, *updaterLogger);
0678           if (!updateRes.ok()) {
0679             ACTS_DEBUG("Update step failed: " << updateRes.error().message());
0680             return updateRes.error();
0681           }
0682           ACTS_VERBOSE("Appended measurement track state with tip = "
0683                        << newBranch.tipIndex());
0684           // Increment number of measurements
0685           newBranch.nMeasurements()++;
0686           newBranch.nDoF() += trackState.calibratedSize();
0687           newBranch.chi2() += trackState.chi2();
0688         } else {
0689           ACTS_WARNING("Cannot handle this track state flags");
0690           continue;
0691         }
0692 
0693         result.activeBranches.push_back(newBranch);
0694 
0695         BranchStopperResult branchStopperResult =
0696             extensions.branchStopper(newBranch, trackState);
0697 
0698         // Check if need to stop this branch
0699         if (branchStopperResult == BranchStopperResult::Continue) {
0700           // Record the number of branches on surface
0701           nBranchesOnSurface++;
0702         } else {
0703           // Record the number of stopped branches
0704           if (branchStopperResult == BranchStopperResult::StopAndKeep) {
0705             storeLastActiveBranch(result);
0706           }
0707           // Remove the branch from list
0708           result.activeBranches.pop_back();
0709         }
0710       }
0711 
0712       return nBranchesOnSurface;
0713     }
0714 
0715     /// @brief CombinatorialKalmanFilter actor operation: add a hole or material track state
0716     ///
0717     /// @param stateMask The bitmask that instructs which components to allocate
0718     /// @param boundState The bound state on current surface
0719     /// @param result is the mutable result state object and which to leave invalid
0720     /// @param isSensitive The surface is sensitive or passive
0721     /// @param expectMeasurements True if measurements where expected for this surface
0722     /// @param prevTip The index of the previous state
0723     ///
0724     /// @return The tip of added state
0725     TrackIndexType addNonSourcelinkState(TrackStatePropMask stateMask,
0726                                          const BoundState& boundState,
0727                                          result_type& result, bool isSensitive,
0728                                          bool expectMeasurements,
0729                                          TrackIndexType prevTip) const {
0730       using PM = TrackStatePropMask;
0731 
0732       // Add a track state
0733       auto trackStateProxy =
0734           result.trackStates->makeTrackState(stateMask, prevTip);
0735       ACTS_VERBOSE("Create "
0736                    << (isSensitive
0737                            ? (expectMeasurements ? "Hole"
0738                                                  : "noMeasurementExpected")
0739                            : "Material")
0740                    << " output track state #" << trackStateProxy.index()
0741                    << " with mask: " << stateMask);
0742 
0743       const auto& [boundParams, jacobian, pathLength] = boundState;
0744       // Fill the track state
0745       trackStateProxy.predicted() = boundParams.parameters();
0746       trackStateProxy.predictedCovariance() = boundParams.covariance().value();
0747       trackStateProxy.jacobian() = jacobian;
0748       trackStateProxy.pathLength() = pathLength;
0749       // Set the surface
0750       trackStateProxy.setReferenceSurface(
0751           boundParams.referenceSurface().getSharedPtr());
0752 
0753       // Set the track state flags
0754       auto typeFlags = trackStateProxy.typeFlags();
0755       if (trackStateProxy.referenceSurface().surfaceMaterial() != nullptr) {
0756         typeFlags.setHasMaterial();
0757       }
0758       typeFlags.setHasParameters();
0759       if (isSensitive) {
0760         if (expectMeasurements) {
0761           typeFlags.setIsHole();
0762         } else {
0763           typeFlags.setHasNoExpectedHit();
0764         }
0765       }
0766 
0767       // Set the filtered parameter index to be the same with predicted
0768       // parameter
0769       trackStateProxy.shareFrom(PM::Predicted, PM::Filtered);
0770 
0771       return trackStateProxy.index();
0772     }
0773 
0774     void storeLastActiveBranch(result_type& result) const {
0775       auto currentBranch = result.activeBranches.back();
0776       TrackIndexType currentTip = currentBranch.tipIndex();
0777 
0778       ACTS_VERBOSE("Storing track "
0779                    << currentBranch.index() << " with tip index " << currentTip
0780                    << ". nMeasurements = " << currentBranch.nMeasurements()
0781                    << ", nOutliers = " << currentBranch.nOutliers()
0782                    << ", nHoles = " << currentBranch.nHoles());
0783 
0784       result.collectedTracks.push_back(currentBranch);
0785     }
0786   };
0787 
0788   /// Void path limit reached aborter to replace the default since the path
0789   /// limit is handled in the CKF actor internally.
0790   struct StubPathLimitReached {
0791     double internalLimit{};
0792 
0793     template <typename propagator_state_t, typename stepper_t,
0794               typename navigator_t>
0795     bool checkAbort(propagator_state_t& /*state*/, const stepper_t& /*stepper*/,
0796                     const navigator_t& /*navigator*/,
0797                     const Logger& /*logger*/) const {
0798       return false;
0799     }
0800   };
0801 
0802  public:
0803   /// Combinatorial Kalman Filter implementation, calls the Kalman filter
0804   ///
0805   /// @param initialParameters The initial track parameters
0806   /// @param tfOptions CombinatorialKalmanFilterOptions steering the track
0807   ///                  finding
0808   /// @param trackContainer Track container in which to store the results
0809   /// @param rootBranch The track to be used as the root branch
0810   ///
0811   /// @note The input measurements are given in the form of @c SourceLinks.
0812   ///       It's @c calibrator_t's job to turn them into calibrated measurements
0813   ///       used in the track finding.
0814   ///
0815   /// @return a container of track finding result for all the initial track
0816   /// parameters
0817   auto findTracks(
0818       const BoundTrackParameters& initialParameters,
0819       const CombinatorialKalmanFilterOptions<track_container_t>& tfOptions,
0820       track_container_t& trackContainer,
0821       typename track_container_t::TrackProxy rootBranch) const
0822       -> Result<std::vector<
0823           typename std::decay_t<decltype(trackContainer)>::TrackProxy>> {
0824     // Create the ActorList
0825     using CombinatorialKalmanFilterActor = Actor;
0826     using Actors = ActorList<CombinatorialKalmanFilterActor>;
0827 
0828     // Create relevant options for the propagation options
0829     using PropagatorOptions = typename propagator_t::template Options<Actors>;
0830     PropagatorOptions propOptions(tfOptions.geoContext,
0831                                   tfOptions.magFieldContext);
0832 
0833     // Set the trivial propagator options
0834     propOptions.setPlainOptions(tfOptions.propagatorPlainOptions);
0835 
0836     // Catch the actor
0837     auto& combKalmanActor =
0838         propOptions.actorList.template get<CombinatorialKalmanFilterActor>();
0839     combKalmanActor.targetReached.surface = tfOptions.targetSurface;
0840     combKalmanActor.multipleScattering = tfOptions.multipleScattering;
0841     combKalmanActor.energyLoss = tfOptions.energyLoss;
0842     combKalmanActor.skipPrePropagationUpdate =
0843         tfOptions.skipPrePropagationUpdate;
0844     combKalmanActor.actorLogger = m_actorLogger.get();
0845     combKalmanActor.updaterLogger = m_updaterLogger.get();
0846     combKalmanActor.calibrationContextPtr = &tfOptions.calibrationContext.get();
0847 
0848     // copy delegates to calibrator, updater, branch stopper
0849     combKalmanActor.extensions = tfOptions.extensions;
0850 
0851     auto propState =
0852         m_propagator
0853             .template makeState<PropagatorOptions, StubPathLimitReached>(
0854                 propOptions);
0855 
0856     auto initResult = m_propagator.template initialize<
0857         decltype(propState), BoundTrackParameters, StubPathLimitReached>(
0858         propState, initialParameters);
0859     if (!initResult.ok()) {
0860       ACTS_DEBUG("Propagation initialization failed: " << initResult.error());
0861       return initResult.error();
0862     }
0863 
0864     auto& r =
0865         propState
0866             .template get<CombinatorialKalmanFilterResult<track_container_t>>();
0867     r.tracks = &trackContainer;
0868     r.trackStates = &trackContainer.trackStateContainer();
0869 
0870     // make sure the right particle hypothesis is set on the root branch
0871     rootBranch.setParticleHypothesis(initialParameters.particleHypothesis());
0872 
0873     r.activeBranches.push_back(rootBranch);
0874 
0875     auto propagationResult = m_propagator.propagate(propState);
0876 
0877     auto result = m_propagator.makeResult(
0878         std::move(propState), propagationResult, propOptions, false);
0879 
0880     if (!result.ok()) {
0881       ACTS_DEBUG("Propagation failed: " << result.error() << " "
0882                                         << result.error().message()
0883                                         << " with the initial parameters: \n"
0884                                         << initialParameters.parameters());
0885       return result.error();
0886     }
0887 
0888     auto& propRes = *result;
0889 
0890     // Get the result of the CombinatorialKalmanFilter
0891     auto combKalmanResult =
0892         std::move(propRes.template get<
0893                   CombinatorialKalmanFilterResult<track_container_t>>());
0894 
0895     // Check if track finding finished properly
0896     if (!combKalmanResult.finished) {
0897       ACTS_DEBUG("CombinatorialKalmanFilter failed: "
0898                  << "Propagation reached max steps "
0899                  << "with the initial parameters: "
0900                  << initialParameters.parameters().transpose());
0901       return CombinatorialKalmanFilterError::PropagationReachesMaxSteps;
0902     }
0903 
0904     return std::move(combKalmanResult.collectedTracks);
0905   }
0906 
0907   /// Combinatorial Kalman Filter implementation, calls the Kalman filter
0908   ///
0909   /// @param initialParameters The initial track parameters
0910   /// @param tfOptions CombinatorialKalmanFilterOptions steering the track
0911   ///                  finding
0912   /// @param trackContainer Track container in which to store the results
0913   /// @note The input measurements are given in the form of @c SourceLinks.
0914   ///       It's @c calibrator_t's job to turn them into calibrated measurements
0915   ///       used in the track finding.
0916   ///
0917   /// @return a container of track finding result for all the initial track
0918   /// parameters
0919   auto findTracks(
0920       const BoundTrackParameters& initialParameters,
0921       const CombinatorialKalmanFilterOptions<track_container_t>& tfOptions,
0922       track_container_t& trackContainer) const
0923       -> Result<std::vector<
0924           typename std::decay_t<decltype(trackContainer)>::TrackProxy>> {
0925     auto rootBranch = trackContainer.makeTrack();
0926     return findTracks(initialParameters, tfOptions, trackContainer, rootBranch);
0927   }
0928 };  // namespace Acts
0929 
0930 /// @}
0931 
0932 }  // namespace Acts