File indexing completed on 2025-07-13 07:50:23
0001
0002
0003
0004
0005
0006
0007
0008
0009 #pragma once
0010
0011
0012 #include "Acts/Utilities/detail/ReferenceWrapperAnyCompat.hpp"
0013
0014 #include "Acts/EventData/MultiTrajectory.hpp"
0015 #include "Acts/EventData/MultiTrajectoryHelpers.hpp"
0016 #include "Acts/EventData/SourceLink.hpp"
0017 #include "Acts/EventData/TrackParameters.hpp"
0018 #include "Acts/EventData/TrackStateType.hpp"
0019 #include "Acts/EventData/VectorMultiTrajectory.hpp"
0020 #include "Acts/Geometry/GeometryContext.hpp"
0021 #include "Acts/MagneticField/MagneticFieldContext.hpp"
0022 #include "Acts/Propagator/ActorList.hpp"
0023 #include "Acts/Propagator/DirectNavigator.hpp"
0024 #include "Acts/Propagator/PropagatorOptions.hpp"
0025 #include "Acts/Propagator/StandardAborters.hpp"
0026 #include "Acts/Propagator/detail/PointwiseMaterialInteraction.hpp"
0027 #include "Acts/TrackFitting/KalmanFitterError.hpp"
0028 #include "Acts/TrackFitting/detail/KalmanUpdateHelpers.hpp"
0029 #include "Acts/TrackFitting/detail/VoidFitterComponents.hpp"
0030 #include "Acts/Utilities/CalibrationContext.hpp"
0031 #include "Acts/Utilities/Delegate.hpp"
0032 #include "Acts/Utilities/Logger.hpp"
0033 #include "Acts/Utilities/Result.hpp"
0034 #include "Acts/Utilities/TrackHelpers.hpp"
0035
0036 #include <functional>
0037 #include <limits>
0038 #include <map>
0039 #include <memory>
0040 #include <type_traits>
0041
0042 namespace Acts {
0043
0044 enum class KalmanFitterTargetSurfaceStrategy {
0045
0046 first,
0047
0048 last,
0049
0050
0051 firstOrLast,
0052 };
0053
0054
0055 template <typename traj_t>
0056 struct KalmanFitterExtensions {
0057 using TrackStateProxy = typename traj_t::TrackStateProxy;
0058 using ConstTrackStateProxy = typename traj_t::ConstTrackStateProxy;
0059 using Parameters = typename TrackStateProxy::Parameters;
0060
0061 using Calibrator =
0062 Delegate<void(const GeometryContext&, const CalibrationContext&,
0063 const SourceLink&, TrackStateProxy)>;
0064
0065 using Smoother = Delegate<Result<void>(const GeometryContext&, traj_t&,
0066 std::size_t, const Logger&)>;
0067
0068 using Updater = Delegate<Result<void>(const GeometryContext&, TrackStateProxy,
0069 const Logger&)>;
0070
0071 using OutlierFinder = Delegate<bool(ConstTrackStateProxy)>;
0072
0073 using ReverseFilteringLogic = Delegate<bool(ConstTrackStateProxy)>;
0074
0075
0076
0077
0078 Calibrator calibrator;
0079
0080
0081 Updater updater;
0082
0083
0084 Smoother smoother;
0085
0086
0087
0088 OutlierFinder outlierFinder;
0089
0090
0091
0092 ReverseFilteringLogic reverseFilteringLogic;
0093
0094
0095 SourceLinkSurfaceAccessor surfaceAccessor;
0096
0097
0098 KalmanFitterExtensions() {
0099 calibrator.template connect<&detail::voidFitterCalibrator<traj_t>>();
0100 updater.template connect<&detail::voidFitterUpdater<traj_t>>();
0101 smoother.template connect<&detail::voidFitterSmoother<traj_t>>();
0102 outlierFinder.template connect<&detail::voidOutlierFinder<traj_t>>();
0103 reverseFilteringLogic
0104 .template connect<&detail::voidReverseFilteringLogic<traj_t>>();
0105 surfaceAccessor.connect<&detail::voidSurfaceAccessor>();
0106 }
0107 };
0108
0109
0110
0111
0112 template <typename traj_t>
0113 struct KalmanFitterOptions {
0114
0115
0116
0117
0118
0119
0120
0121
0122
0123
0124
0125
0126
0127 KalmanFitterOptions(const GeometryContext& gctx,
0128 const MagneticFieldContext& mctx,
0129 std::reference_wrapper<const CalibrationContext> cctx,
0130 KalmanFitterExtensions<traj_t> extensions_,
0131 const PropagatorPlainOptions& pOptions,
0132 const Surface* rSurface = nullptr,
0133 bool mScattering = true, bool eLoss = true,
0134 bool rFiltering = false, double rfScaling = 1.0,
0135 const FreeToBoundCorrection& freeToBoundCorrection_ =
0136 FreeToBoundCorrection(false))
0137 : geoContext(gctx),
0138 magFieldContext(mctx),
0139 calibrationContext(cctx),
0140 extensions(extensions_),
0141 propagatorPlainOptions(pOptions),
0142 referenceSurface(rSurface),
0143 multipleScattering(mScattering),
0144 energyLoss(eLoss),
0145 reversedFiltering(rFiltering),
0146 reversedFilteringCovarianceScaling(rfScaling),
0147 freeToBoundCorrection(freeToBoundCorrection_) {}
0148
0149 KalmanFitterOptions() = delete;
0150
0151
0152 std::reference_wrapper<const GeometryContext> geoContext;
0153
0154 std::reference_wrapper<const MagneticFieldContext> magFieldContext;
0155
0156 std::reference_wrapper<const CalibrationContext> calibrationContext;
0157
0158 KalmanFitterExtensions<traj_t> extensions;
0159
0160
0161 PropagatorPlainOptions propagatorPlainOptions;
0162
0163
0164 const Surface* referenceSurface = nullptr;
0165
0166
0167 KalmanFitterTargetSurfaceStrategy referenceSurfaceStrategy =
0168 KalmanFitterTargetSurfaceStrategy::firstOrLast;
0169
0170
0171 bool multipleScattering = true;
0172
0173
0174 bool energyLoss = true;
0175
0176
0177
0178 bool reversedFiltering = false;
0179
0180
0181
0182
0183
0184 double reversedFilteringCovarianceScaling = 1.0;
0185
0186
0187
0188 FreeToBoundCorrection freeToBoundCorrection;
0189 };
0190
0191 template <typename traj_t>
0192 struct KalmanFitterResult {
0193
0194 traj_t* fittedStates{nullptr};
0195
0196
0197
0198
0199
0200 std::size_t lastMeasurementIndex = Acts::MultiTrajectoryTraits::kInvalid;
0201
0202
0203
0204
0205
0206 std::size_t lastTrackIndex = Acts::MultiTrajectoryTraits::kInvalid;
0207
0208
0209 std::optional<BoundTrackParameters> fittedParameters;
0210
0211
0212 std::size_t measurementStates = 0;
0213
0214
0215
0216
0217
0218 std::size_t measurementHoles = 0;
0219
0220
0221 std::size_t processedStates = 0;
0222
0223
0224 bool smoothed = false;
0225
0226
0227 bool reversed = false;
0228
0229
0230 bool finished = false;
0231
0232
0233 std::vector<const Surface*> missedActiveSurfaces;
0234
0235
0236 std::vector<const Surface*> passedAgainSurfaces;
0237
0238 Result<void> result{Result<void>::success()};
0239 };
0240
0241
0242
0243
0244
0245
0246
0247
0248
0249
0250
0251
0252
0253
0254
0255
0256
0257
0258
0259
0260
0261 template <typename propagator_t, typename traj_t>
0262 class KalmanFitter {
0263
0264 using KalmanNavigator = typename propagator_t::Navigator;
0265
0266
0267 static constexpr bool isDirectNavigator =
0268 std::is_same_v<KalmanNavigator, DirectNavigator>;
0269
0270 public:
0271 explicit KalmanFitter(propagator_t pPropagator,
0272 std::unique_ptr<const Logger> _logger =
0273 getDefaultLogger("KalmanFitter", Logging::INFO))
0274 : m_propagator(std::move(pPropagator)),
0275 m_logger{std::move(_logger)},
0276 m_actorLogger{m_logger->cloneWithSuffix("Actor")} {}
0277
0278 private:
0279
0280 propagator_t m_propagator;
0281
0282
0283 std::unique_ptr<const Logger> m_logger;
0284 std::unique_ptr<const Logger> m_actorLogger;
0285
0286 const Logger& logger() const { return *m_logger; }
0287
0288
0289
0290
0291
0292
0293
0294
0295
0296 template <typename parameters_t>
0297 class Actor {
0298 public:
0299
0300 using result_type = KalmanFitterResult<traj_t>;
0301
0302
0303 SurfaceReached targetReached{std::numeric_limits<double>::lowest()};
0304
0305
0306 KalmanFitterTargetSurfaceStrategy targetSurfaceStrategy =
0307 KalmanFitterTargetSurfaceStrategy::firstOrLast;
0308
0309
0310 const std::map<GeometryIdentifier, SourceLink>* inputMeasurements = nullptr;
0311
0312
0313 bool multipleScattering = true;
0314
0315
0316 bool energyLoss = true;
0317
0318
0319 bool reversedFiltering = false;
0320
0321
0322 double reversedFilteringCovarianceScaling = 1.0;
0323
0324
0325
0326 FreeToBoundCorrection freeToBoundCorrection;
0327
0328
0329 std::shared_ptr<traj_t> outputStates;
0330
0331
0332 const Logger* actorLogger{nullptr};
0333
0334
0335 const Logger& logger() const { return *actorLogger; }
0336
0337 KalmanFitterExtensions<traj_t> extensions;
0338
0339
0340 const CalibrationContext* calibrationContext{nullptr};
0341
0342
0343
0344
0345
0346
0347
0348
0349
0350
0351
0352 template <typename propagator_state_t, typename stepper_t,
0353 typename navigator_t>
0354 void act(propagator_state_t& state, const stepper_t& stepper,
0355 const navigator_t& navigator, result_type& result,
0356 const Logger& ) const {
0357 assert(result.fittedStates && "No MultiTrajectory set");
0358
0359 if (result.finished) {
0360 return;
0361 }
0362
0363 ACTS_VERBOSE("KalmanFitter step at pos: "
0364 << stepper.position(state.stepping).transpose()
0365 << " dir: " << stepper.direction(state.stepping).transpose()
0366 << " momentum: "
0367 << stepper.absoluteMomentum(state.stepping));
0368
0369
0370
0371 auto surface = navigator.currentSurface(state.navigation);
0372 std::string direction = state.options.direction.toString();
0373 if (surface != nullptr) {
0374
0375
0376
0377
0378
0379
0380
0381 if (!result.smoothed && !result.reversed) {
0382 ACTS_VERBOSE("Perform " << direction << " filter step");
0383 auto res = filter(surface, state, stepper, navigator, result);
0384 if (!res.ok()) {
0385 ACTS_ERROR("Error in " << direction << " filter: " << res.error());
0386 result.result = res.error();
0387 }
0388 }
0389 if (result.reversed) {
0390 ACTS_VERBOSE("Perform " << direction << " filter step");
0391 auto res = reversedFilter(surface, state, stepper, navigator, result);
0392 if (!res.ok()) {
0393 ACTS_ERROR("Error in " << direction << " filter: " << res.error());
0394 result.result = res.error();
0395 }
0396 }
0397 }
0398
0399
0400
0401
0402
0403 if (!result.smoothed && !result.reversed) {
0404 if (result.measurementStates == inputMeasurements->size() ||
0405 (result.measurementStates > 0 &&
0406 navigator.navigationBreak(state.navigation))) {
0407
0408 result.missedActiveSurfaces.resize(result.measurementHoles);
0409
0410 typename traj_t::ConstTrackStateProxy trackStateProxy{
0411 result.fittedStates->getTrackState(result.lastMeasurementIndex)};
0412 if (reversedFiltering ||
0413 extensions.reverseFilteringLogic(trackStateProxy)) {
0414
0415
0416
0417 ACTS_VERBOSE("Reverse navigation direction.");
0418 auto res = reverse(state, stepper, navigator, result);
0419 if (!res.ok()) {
0420 ACTS_ERROR("Error in reversing navigation: " << res.error());
0421 result.result = res.error();
0422 }
0423 } else {
0424
0425
0426
0427
0428 ACTS_VERBOSE("Finalize/run smoothing");
0429 auto res = finalize(state, stepper, navigator, result);
0430 if (!res.ok()) {
0431 ACTS_ERROR("Error in finalize: " << res.error());
0432 result.result = res.error();
0433 }
0434 }
0435 }
0436 }
0437
0438
0439
0440
0441 if (result.smoothed || result.reversed) {
0442 if (result.smoothed) {
0443
0444
0445 materialInteractor(navigator.currentSurface(state.navigation), state,
0446 stepper, navigator,
0447 MaterialUpdateStage::FullUpdate);
0448 }
0449
0450 if (targetReached.surface == nullptr) {
0451
0452
0453
0454 if (result.reversed) {
0455 ACTS_ERROR(
0456 "The target surface needed for aborting reversed propagation "
0457 "is not provided");
0458 result.result =
0459 Result<void>(KalmanFitterError::ReversePropagationFailed);
0460 } else {
0461 ACTS_VERBOSE(
0462 "No target surface set. Completing without fitted track "
0463 "parameter");
0464
0465 result.finished = true;
0466 }
0467 } else if (targetReached.checkAbort(state, stepper, navigator,
0468 logger())) {
0469 ACTS_VERBOSE("Completing with fitted track parameter");
0470
0471 auto res = stepper.boundState(state.stepping, *targetReached.surface,
0472 true, freeToBoundCorrection);
0473 if (!res.ok()) {
0474 ACTS_ERROR("Error in " << direction << " filter: " << res.error());
0475 result.result = res.error();
0476 return;
0477 }
0478 auto& fittedState = *res;
0479
0480 result.fittedParameters = std::get<BoundTrackParameters>(fittedState);
0481
0482
0483 if (result.reversed) {
0484 result.fittedStates->applyBackwards(
0485 result.lastMeasurementIndex, [&](auto trackState) {
0486 auto fSurface = &trackState.referenceSurface();
0487 if (!rangeContainsValue(result.passedAgainSurfaces,
0488 fSurface)) {
0489
0490
0491 trackState.unset(TrackStatePropMask::Smoothed);
0492 trackState.typeFlags().set(TrackStateFlag::OutlierFlag);
0493 }
0494 });
0495 }
0496
0497 result.finished = true;
0498 }
0499 }
0500 }
0501
0502 template <typename propagator_state_t, typename stepper_t,
0503 typename navigator_t>
0504 bool checkAbort(propagator_state_t& , const stepper_t& ,
0505 const navigator_t& , const result_type& result,
0506 const Logger& ) const {
0507 return (!result.result.ok() || result.finished);
0508 }
0509
0510
0511
0512
0513
0514
0515
0516
0517
0518
0519
0520 template <typename propagator_state_t, typename stepper_t,
0521 typename navigator_t>
0522 Result<void> reverse(propagator_state_t& state, const stepper_t& stepper,
0523 const navigator_t& navigator,
0524 result_type& result) const {
0525
0526 if (result.lastMeasurementIndex ==
0527 Acts::MultiTrajectoryTraits::kInvalid) {
0528 ACTS_ERROR("No point to reverse for a track without measurements.");
0529 return KalmanFitterError::ReversePropagationFailed;
0530 }
0531
0532
0533 result.reversed = true;
0534
0535
0536 state.options.direction = state.options.direction.invert();
0537
0538
0539
0540 state.options.pathLimit =
0541 state.options.direction * std::abs(state.options.pathLimit);
0542
0543
0544
0545 auto st = result.fittedStates->getTrackState(result.lastMeasurementIndex);
0546
0547
0548 stepper.initialize(
0549 state.stepping, st.filtered(),
0550 reversedFilteringCovarianceScaling * st.filteredCovariance(),
0551 stepper.particleHypothesis(state.stepping), st.referenceSurface());
0552
0553
0554 st.smoothed() = st.filtered();
0555 st.smoothedCovariance() = st.filteredCovariance();
0556 result.passedAgainSurfaces.push_back(&st.referenceSurface());
0557
0558
0559
0560
0561 state.navigation.options.startSurface = &st.referenceSurface();
0562 state.navigation.options.targetSurface = nullptr;
0563 auto navInitRes = navigator.initialize(
0564 state.navigation, stepper.position(state.stepping),
0565 stepper.direction(state.stepping), state.options.direction);
0566 if (!navInitRes.ok()) {
0567 return navInitRes.error();
0568 }
0569
0570
0571
0572 materialInteractor(navigator.currentSurface(state.navigation), state,
0573 stepper, navigator, MaterialUpdateStage::FullUpdate);
0574
0575 return Result<void>::success();
0576 }
0577
0578
0579
0580
0581
0582
0583
0584
0585
0586
0587
0588
0589 template <typename propagator_state_t, typename stepper_t,
0590 typename navigator_t>
0591 Result<void> filter(const Surface* surface, propagator_state_t& state,
0592 const stepper_t& stepper, const navigator_t& navigator,
0593 result_type& result) const {
0594 const bool precedingMeasurementExists = result.measurementStates > 0;
0595 const bool surfaceIsSensitive =
0596 surface->associatedDetectorElement() != nullptr;
0597 const bool surfaceHasMaterial = surface->surfaceMaterial() != nullptr;
0598
0599
0600 auto sourceLinkIt = inputMeasurements->find(surface->geometryId());
0601 if (sourceLinkIt != inputMeasurements->end()) {
0602
0603 ACTS_VERBOSE("Measurement surface " << surface->geometryId()
0604 << " detected.");
0605
0606 stepper.transportCovarianceToBound(state.stepping, *surface,
0607 freeToBoundCorrection);
0608
0609
0610 materialInteractor(surface, state, stepper, navigator,
0611 MaterialUpdateStage::PreUpdate);
0612
0613
0614
0615 auto trackStateProxyRes = detail::kalmanHandleMeasurement(
0616 *calibrationContext, state, stepper, extensions, *surface,
0617 sourceLinkIt->second, *result.fittedStates, result.lastTrackIndex,
0618 false, logger());
0619
0620 if (!trackStateProxyRes.ok()) {
0621 return trackStateProxyRes.error();
0622 }
0623
0624 const auto& trackStateProxy = *trackStateProxyRes;
0625 result.lastTrackIndex = trackStateProxy.index();
0626
0627
0628 if (trackStateProxy.typeFlags().test(
0629 Acts::TrackStateFlag::MeasurementFlag)) {
0630
0631 ACTS_VERBOSE("Filtering step successful, updated parameters are:\n"
0632 << trackStateProxy.filtered().transpose());
0633
0634 stepper.update(state.stepping,
0635 MultiTrajectoryHelpers::freeFiltered(
0636 state.options.geoContext, trackStateProxy),
0637 trackStateProxy.filtered(),
0638 trackStateProxy.filteredCovariance(), *surface);
0639
0640 ++result.measurementStates;
0641 }
0642
0643
0644 materialInteractor(surface, state, stepper, navigator,
0645 MaterialUpdateStage::PostUpdate);
0646
0647 ++result.processedStates;
0648
0649
0650 result.measurementHoles = result.missedActiveSurfaces.size();
0651
0652
0653 result.lastMeasurementIndex = result.lastTrackIndex;
0654
0655 } else if ((precedingMeasurementExists && surfaceIsSensitive) ||
0656 surfaceHasMaterial) {
0657
0658
0659
0660 auto trackStateProxyRes = detail::kalmanHandleNoMeasurement(
0661 state.stepping, stepper, *surface, *result.fittedStates,
0662 result.lastTrackIndex, true, logger(), precedingMeasurementExists,
0663 freeToBoundCorrection);
0664
0665 if (!trackStateProxyRes.ok()) {
0666 return trackStateProxyRes.error();
0667 }
0668
0669 const auto& trackStateProxy = *trackStateProxyRes;
0670 result.lastTrackIndex = trackStateProxy.index();
0671
0672 if (trackStateProxy.typeFlags().test(TrackStateFlag::HoleFlag)) {
0673
0674 result.missedActiveSurfaces.push_back(surface);
0675 }
0676
0677 ++result.processedStates;
0678
0679
0680 materialInteractor(surface, state, stepper, navigator,
0681 MaterialUpdateStage::FullUpdate);
0682 }
0683 return Result<void>::success();
0684 }
0685
0686
0687
0688
0689
0690
0691
0692
0693
0694
0695
0696
0697 template <typename propagator_state_t, typename stepper_t,
0698 typename navigator_t>
0699 Result<void> reversedFilter(const Surface* surface,
0700 propagator_state_t& state,
0701 const stepper_t& stepper,
0702 const navigator_t& navigator,
0703 result_type& result) const {
0704
0705 auto sourceLinkIt = inputMeasurements->find(surface->geometryId());
0706 if (sourceLinkIt != inputMeasurements->end()) {
0707
0708 ACTS_VERBOSE("Measurement surface "
0709 << surface->geometryId()
0710 << " detected in reversed propagation.");
0711
0712
0713
0714 if (result.reversed &&
0715 surface == navigator.startSurface(state.navigation)) {
0716 materialInteractor(surface, state, stepper, navigator,
0717 MaterialUpdateStage::FullUpdate);
0718 return Result<void>::success();
0719 }
0720
0721
0722 stepper.transportCovarianceToBound(state.stepping, *surface,
0723 freeToBoundCorrection);
0724
0725
0726 materialInteractor(surface, state, stepper, navigator,
0727 MaterialUpdateStage::PreUpdate);
0728
0729 auto fittedStates = *result.fittedStates;
0730
0731
0732
0733 TrackStatePropMask mask =
0734 TrackStatePropMask::Predicted | TrackStatePropMask::Filtered |
0735 TrackStatePropMask::Smoothed | TrackStatePropMask::Jacobian |
0736 TrackStatePropMask::Calibrated;
0737 const std::size_t currentTrackIndex = fittedStates.addTrackState(
0738 mask, Acts::MultiTrajectoryTraits::kInvalid);
0739
0740
0741 typename traj_t::TrackStateProxy trackStateProxy =
0742 fittedStates.getTrackState(currentTrackIndex);
0743
0744
0745
0746 {
0747 trackStateProxy.setReferenceSurface(surface->getSharedPtr());
0748
0749 auto res = stepper.boundState(state.stepping, *surface, false,
0750 freeToBoundCorrection);
0751 if (!res.ok()) {
0752 return res.error();
0753 }
0754 const auto& [boundParams, jacobian, pathLength] = *res;
0755
0756
0757 trackStateProxy.predicted() = boundParams.parameters();
0758 trackStateProxy.predictedCovariance() = state.stepping.cov;
0759
0760 trackStateProxy.jacobian() = jacobian;
0761 trackStateProxy.pathLength() = pathLength;
0762 }
0763
0764
0765
0766 extensions.calibrator(state.geoContext, *calibrationContext,
0767 sourceLinkIt->second, trackStateProxy);
0768
0769
0770 auto updateRes =
0771 extensions.updater(state.geoContext, trackStateProxy, logger());
0772 if (!updateRes.ok()) {
0773 ACTS_ERROR("Backward update step failed: " << updateRes.error());
0774 return updateRes.error();
0775 } else {
0776
0777 ACTS_VERBOSE(
0778 "Backward filtering step successful, updated parameters are:\n"
0779 << trackStateProxy.filtered().transpose());
0780
0781
0782 result.fittedStates->applyBackwards(
0783 result.lastMeasurementIndex, [&](auto trackState) {
0784 auto fSurface = &trackState.referenceSurface();
0785 if (fSurface == surface) {
0786 result.passedAgainSurfaces.push_back(surface);
0787 trackState.smoothed() = trackStateProxy.filtered();
0788 trackState.smoothedCovariance() =
0789 trackStateProxy.filteredCovariance();
0790 return false;
0791 }
0792 return true;
0793 });
0794
0795
0796
0797
0798 stepper.update(state.stepping,
0799 MultiTrajectoryHelpers::freeFiltered(
0800 state.options.geoContext, trackStateProxy),
0801 trackStateProxy.filtered(),
0802 trackStateProxy.filteredCovariance(), *surface);
0803
0804
0805 materialInteractor(surface, state, stepper, navigator,
0806 MaterialUpdateStage::PostUpdate);
0807 }
0808 } else if (surface->associatedDetectorElement() != nullptr ||
0809 surface->surfaceMaterial() != nullptr) {
0810
0811 if (surface->associatedDetectorElement() != nullptr) {
0812 ACTS_VERBOSE("Detected hole on " << surface->geometryId()
0813 << " in reversed filtering");
0814 if (state.stepping.covTransport) {
0815 stepper.transportCovarianceToBound(state.stepping, *surface);
0816 }
0817 } else if (surface->surfaceMaterial() != nullptr) {
0818 ACTS_VERBOSE("Detected in-sensitive surface "
0819 << surface->geometryId() << " in reversed filtering");
0820 if (state.stepping.covTransport) {
0821 stepper.transportCovarianceToCurvilinear(state.stepping);
0822 }
0823 }
0824
0825
0826 stepper.setIdentityJacobian(state.stepping);
0827 if (surface->surfaceMaterial() != nullptr) {
0828
0829 materialInteractor(surface, state, stepper, navigator,
0830 MaterialUpdateStage::FullUpdate);
0831 }
0832 }
0833
0834 return Result<void>::success();
0835 }
0836
0837
0838
0839
0840
0841
0842
0843
0844
0845
0846
0847
0848
0849 template <typename propagator_state_t, typename stepper_t,
0850 typename navigator_t>
0851 void materialInteractor(const Surface* surface, propagator_state_t& state,
0852 const stepper_t& stepper,
0853 const navigator_t& navigator,
0854 const MaterialUpdateStage& updateStage) const {
0855
0856 if (!surface) {
0857 ACTS_VERBOSE(
0858 "Surface is nullptr. Cannot be used for material interaction");
0859 return;
0860 }
0861
0862
0863 bool hasMaterial = false;
0864
0865 if (surface && surface->surfaceMaterial()) {
0866
0867 detail::PointwiseMaterialInteraction interaction(surface, state,
0868 stepper);
0869
0870 if (interaction.evaluateMaterialSlab(state, navigator, updateStage)) {
0871
0872 hasMaterial = true;
0873
0874
0875 interaction.evaluatePointwiseMaterialInteraction(multipleScattering,
0876 energyLoss);
0877
0878
0879 ACTS_VERBOSE("Material effects on surface: "
0880 << surface->geometryId()
0881 << " at update stage: " << updateStage << " are :");
0882 ACTS_VERBOSE("eLoss = "
0883 << interaction.Eloss << ", "
0884 << "variancePhi = " << interaction.variancePhi << ", "
0885 << "varianceTheta = " << interaction.varianceTheta
0886 << ", "
0887 << "varianceQoverP = " << interaction.varianceQoverP);
0888
0889
0890 interaction.updateState(state, stepper, addNoise);
0891 }
0892 }
0893
0894 if (!hasMaterial) {
0895
0896 ACTS_VERBOSE("No material effects on surface: " << surface->geometryId()
0897 << " at update stage: "
0898 << updateStage);
0899 }
0900 }
0901
0902
0903
0904
0905
0906
0907
0908
0909
0910
0911
0912 template <typename propagator_state_t, typename stepper_t,
0913 typename navigator_t>
0914 Result<void> finalize(propagator_state_t& state, const stepper_t& stepper,
0915 const navigator_t& navigator,
0916 result_type& result) const {
0917
0918 result.smoothed = true;
0919
0920
0921
0922 std::size_t firstStateIndex = result.lastMeasurementIndex;
0923
0924 std::size_t nStates = 0;
0925 result.fittedStates->applyBackwards(
0926 result.lastMeasurementIndex, [&](auto st) {
0927 bool isMeasurement =
0928 st.typeFlags().test(TrackStateFlag::MeasurementFlag);
0929 bool isOutlier = st.typeFlags().test(TrackStateFlag::OutlierFlag);
0930
0931
0932
0933 if (isMeasurement && !isOutlier) {
0934 firstStateIndex = st.index();
0935 }
0936 nStates++;
0937 });
0938
0939
0940 if (nStates == 0) {
0941 ACTS_ERROR("Smoothing for a track without measurements.");
0942 return KalmanFitterError::SmoothFailed;
0943 }
0944
0945 ACTS_VERBOSE("Apply smoothing on " << nStates
0946 << " filtered track states.");
0947
0948
0949 auto smoothRes =
0950 extensions.smoother(state.geoContext, *result.fittedStates,
0951 result.lastMeasurementIndex, logger());
0952 if (!smoothRes.ok()) {
0953 ACTS_ERROR("Smoothing step failed: " << smoothRes.error());
0954 return smoothRes.error();
0955 }
0956
0957
0958 if (targetReached.surface == nullptr) {
0959 return Result<void>::success();
0960 }
0961
0962
0963 auto firstCreatedState =
0964 result.fittedStates->getTrackState(firstStateIndex);
0965 auto lastCreatedMeasurement =
0966 result.fittedStates->getTrackState(result.lastMeasurementIndex);
0967
0968
0969 auto target = [&](const FreeVector& freeVector) -> SurfaceIntersection {
0970 return targetReached.surface
0971 ->intersect(
0972 state.geoContext, freeVector.segment<3>(eFreePos0),
0973 state.options.direction * freeVector.segment<3>(eFreeDir0),
0974 BoundaryTolerance::None(), state.options.surfaceTolerance)
0975 .closest();
0976 };
0977
0978
0979
0980 auto firstParams = MultiTrajectoryHelpers::freeSmoothed(
0981 state.options.geoContext, firstCreatedState);
0982 auto lastParams = MultiTrajectoryHelpers::freeSmoothed(
0983 state.options.geoContext, lastCreatedMeasurement);
0984
0985
0986 const auto firstIntersection = target(firstParams);
0987 const auto lastIntersection = target(lastParams);
0988
0989
0990
0991
0992
0993
0994
0995
0996
0997 bool useFirstTrackState = true;
0998 switch (targetSurfaceStrategy) {
0999 case KalmanFitterTargetSurfaceStrategy::first:
1000 useFirstTrackState = true;
1001 break;
1002 case KalmanFitterTargetSurfaceStrategy::last:
1003 useFirstTrackState = false;
1004 break;
1005 case KalmanFitterTargetSurfaceStrategy::firstOrLast:
1006 useFirstTrackState = std::abs(firstIntersection.pathLength()) <=
1007 std::abs(lastIntersection.pathLength());
1008 break;
1009 default:
1010 ACTS_ERROR("Unknown target surface strategy");
1011 return KalmanFitterError::SmoothFailed;
1012 }
1013 bool reverseDirection = false;
1014 if (useFirstTrackState) {
1015 stepper.initialize(state.stepping, firstCreatedState.smoothed(),
1016 firstCreatedState.smoothedCovariance(),
1017 stepper.particleHypothesis(state.stepping),
1018 firstCreatedState.referenceSurface());
1019 reverseDirection = firstIntersection.pathLength() < 0;
1020 } else {
1021 stepper.initialize(state.stepping, lastCreatedMeasurement.smoothed(),
1022 lastCreatedMeasurement.smoothedCovariance(),
1023 stepper.particleHypothesis(state.stepping),
1024 lastCreatedMeasurement.referenceSurface());
1025 reverseDirection = lastIntersection.pathLength() < 0;
1026 }
1027
1028 if (reverseDirection) {
1029 ACTS_VERBOSE(
1030 "Reverse navigation direction after smoothing for reaching the "
1031 "target surface");
1032 state.options.direction = state.options.direction.invert();
1033 }
1034 const auto& surface = useFirstTrackState
1035 ? firstCreatedState.referenceSurface()
1036 : lastCreatedMeasurement.referenceSurface();
1037
1038 ACTS_VERBOSE(
1039 "Smoothing successful, updating stepping state to smoothed "
1040 "parameters at surface "
1041 << surface.geometryId() << ". Prepared to reach the target surface.");
1042
1043
1044
1045
1046 auto navigationOptions = state.navigation.options;
1047 navigationOptions.startSurface = &surface;
1048 navigationOptions.targetSurface = nullptr;
1049 state.navigation = navigator.makeState(navigationOptions);
1050 auto navInitRes = navigator.initialize(
1051 state.navigation, stepper.position(state.stepping),
1052 stepper.direction(state.stepping), state.options.direction);
1053 if (!navInitRes.ok()) {
1054 return navInitRes.error();
1055 }
1056
1057 return Result<void>::success();
1058 }
1059 };
1060
1061 public:
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080 template <typename source_link_iterator_t, typename start_parameters_t,
1081 typename parameters_t = BoundTrackParameters,
1082 TrackContainerFrontend track_container_t>
1083 Result<typename track_container_t::TrackProxy> fit(
1084 source_link_iterator_t it, source_link_iterator_t end,
1085 const start_parameters_t& sParameters,
1086 const KalmanFitterOptions<traj_t>& kfOptions,
1087 track_container_t& trackContainer) const
1088 requires(!isDirectNavigator)
1089 {
1090
1091
1092 ACTS_VERBOSE("Preparing " << std::distance(it, end)
1093 << " input measurements");
1094 std::map<GeometryIdentifier, SourceLink> inputMeasurements;
1095
1096 for (; it != end; ++it) {
1097 SourceLink sl = *it;
1098 const Surface* surface = kfOptions.extensions.surfaceAccessor(sl);
1099
1100 auto geoId = surface->geometryId();
1101 inputMeasurements.emplace(geoId, std::move(sl));
1102 }
1103
1104
1105 using KalmanActor = Actor<parameters_t>;
1106
1107 using KalmanResult = typename KalmanActor::result_type;
1108 using Actors = ActorList<KalmanActor>;
1109 using PropagatorOptions = typename propagator_t::template Options<Actors>;
1110
1111
1112 PropagatorOptions propagatorOptions(kfOptions.geoContext,
1113 kfOptions.magFieldContext);
1114
1115
1116 propagatorOptions.setPlainOptions(kfOptions.propagatorPlainOptions);
1117
1118
1119
1120 for (const auto& [surfaceId, _] : inputMeasurements) {
1121 propagatorOptions.navigation.insertExternalSurface(surfaceId);
1122 }
1123
1124
1125 auto& kalmanActor = propagatorOptions.actorList.template get<KalmanActor>();
1126 kalmanActor.inputMeasurements = &inputMeasurements;
1127 kalmanActor.targetReached.surface = kfOptions.referenceSurface;
1128 kalmanActor.targetSurfaceStrategy = kfOptions.referenceSurfaceStrategy;
1129 kalmanActor.multipleScattering = kfOptions.multipleScattering;
1130 kalmanActor.energyLoss = kfOptions.energyLoss;
1131 kalmanActor.reversedFiltering = kfOptions.reversedFiltering;
1132 kalmanActor.reversedFilteringCovarianceScaling =
1133 kfOptions.reversedFilteringCovarianceScaling;
1134 kalmanActor.freeToBoundCorrection = kfOptions.freeToBoundCorrection;
1135 kalmanActor.calibrationContext = &kfOptions.calibrationContext.get();
1136 kalmanActor.extensions = kfOptions.extensions;
1137 kalmanActor.actorLogger = m_actorLogger.get();
1138
1139 return fit_impl<start_parameters_t, PropagatorOptions, KalmanResult,
1140 track_container_t>(sParameters, propagatorOptions,
1141 trackContainer);
1142 }
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164 template <typename source_link_iterator_t, typename start_parameters_t,
1165 typename parameters_t = BoundTrackParameters,
1166 TrackContainerFrontend track_container_t>
1167 Result<typename track_container_t::TrackProxy> fit(
1168 source_link_iterator_t it, source_link_iterator_t end,
1169 const start_parameters_t& sParameters,
1170 const KalmanFitterOptions<traj_t>& kfOptions,
1171 const std::vector<const Surface*>& sSequence,
1172 track_container_t& trackContainer) const
1173 requires(isDirectNavigator)
1174 {
1175
1176
1177 ACTS_VERBOSE("Preparing " << std::distance(it, end)
1178 << " input measurements");
1179 std::map<GeometryIdentifier, SourceLink> inputMeasurements;
1180 for (; it != end; ++it) {
1181 SourceLink sl = *it;
1182 const Surface* surface = kfOptions.extensions.surfaceAccessor(sl);
1183
1184 auto geoId = surface->geometryId();
1185 inputMeasurements.emplace(geoId, std::move(sl));
1186 }
1187
1188
1189 using KalmanActor = Actor<parameters_t>;
1190
1191 using KalmanResult = typename KalmanActor::result_type;
1192 using Actors = ActorList<KalmanActor>;
1193 using PropagatorOptions = typename propagator_t::template Options<Actors>;
1194
1195
1196 PropagatorOptions propagatorOptions(kfOptions.geoContext,
1197 kfOptions.magFieldContext);
1198
1199
1200 propagatorOptions.setPlainOptions(kfOptions.propagatorPlainOptions);
1201
1202
1203 auto& kalmanActor = propagatorOptions.actorList.template get<KalmanActor>();
1204 kalmanActor.inputMeasurements = &inputMeasurements;
1205 kalmanActor.targetReached.surface = kfOptions.referenceSurface;
1206 kalmanActor.targetSurfaceStrategy = kfOptions.referenceSurfaceStrategy;
1207 kalmanActor.multipleScattering = kfOptions.multipleScattering;
1208 kalmanActor.energyLoss = kfOptions.energyLoss;
1209 kalmanActor.reversedFiltering = kfOptions.reversedFiltering;
1210 kalmanActor.reversedFilteringCovarianceScaling =
1211 kfOptions.reversedFilteringCovarianceScaling;
1212 kalmanActor.freeToBoundCorrection = kfOptions.freeToBoundCorrection;
1213 kalmanActor.calibrationContext = &kfOptions.calibrationContext.get();
1214 kalmanActor.extensions = kfOptions.extensions;
1215 kalmanActor.actorLogger = m_actorLogger.get();
1216
1217
1218 propagatorOptions.navigation.surfaces = sSequence;
1219
1220 return fit_impl<start_parameters_t, PropagatorOptions, KalmanResult,
1221 track_container_t>(sParameters, propagatorOptions,
1222 trackContainer);
1223 }
1224
1225 private:
1226
1227
1228
1229
1230
1231
1232
1233
1234
1235
1236
1237
1238
1239 template <typename start_parameters_t, typename propagator_options_t,
1240 typename kalman_result_t, TrackContainerFrontend track_container_t>
1241 auto fit_impl(const start_parameters_t& sParameters,
1242 const propagator_options_t& propagatorOptions,
1243 track_container_t& trackContainer) const
1244 -> Result<typename track_container_t::TrackProxy> {
1245 auto propagatorState = m_propagator.makeState(propagatorOptions);
1246
1247 auto propagatorInitResult =
1248 m_propagator.initialize(propagatorState, sParameters);
1249 if (!propagatorInitResult.ok()) {
1250 ACTS_ERROR("Propagation initialization failed: "
1251 << propagatorInitResult.error());
1252 return propagatorInitResult.error();
1253 }
1254
1255 auto& kalmanResult =
1256 propagatorState.template get<KalmanFitterResult<traj_t>>();
1257 kalmanResult.fittedStates = &trackContainer.trackStateContainer();
1258
1259
1260 auto result = m_propagator.propagate(propagatorState);
1261
1262 if (!result.ok()) {
1263 ACTS_ERROR("Propagation failed: " << result.error());
1264 return result.error();
1265 }
1266
1267
1268
1269 if (kalmanResult.result.ok() && !kalmanResult.measurementStates) {
1270 kalmanResult.result = Result<void>(KalmanFitterError::NoMeasurementFound);
1271 }
1272
1273 if (!kalmanResult.result.ok()) {
1274 ACTS_ERROR("KalmanFilter failed: "
1275 << kalmanResult.result.error() << ", "
1276 << kalmanResult.result.error().message());
1277 return kalmanResult.result.error();
1278 }
1279
1280 auto track = trackContainer.makeTrack();
1281 track.tipIndex() = kalmanResult.lastMeasurementIndex;
1282 if (kalmanResult.fittedParameters) {
1283 const auto& params = kalmanResult.fittedParameters.value();
1284 track.parameters() = params.parameters();
1285 track.covariance() = params.covariance().value();
1286 track.setReferenceSurface(params.referenceSurface().getSharedPtr());
1287 }
1288
1289 calculateTrackQuantities(track);
1290
1291 if (trackContainer.hasColumn(hashString("smoothed"))) {
1292 track.template component<bool, hashString("smoothed")>() =
1293 kalmanResult.smoothed;
1294 }
1295
1296 if (trackContainer.hasColumn(hashString("reversed"))) {
1297 track.template component<bool, hashString("reversed")>() =
1298 kalmanResult.reversed;
1299 }
1300
1301
1302 return track;
1303 }
1304 };
1305
1306 }