File indexing completed on 2026-08-25 08:16:39
0001
0002
0003
0004
0005
0006
0007
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
0040
0041
0042
0043 template <typename traj_t>
0044 struct KalmanFitterExtensions {
0045
0046 using TrackStateProxy = typename traj_t::TrackStateProxy;
0047
0048 using ConstTrackStateProxy = typename traj_t::ConstTrackStateProxy;
0049
0050 using Parameters = typename TrackStateProxy::Parameters;
0051
0052
0053 using Calibrator =
0054 Delegate<void(const GeometryContext&, const CalibrationContext&,
0055 const SourceLink&, TrackStateProxy)>;
0056
0057
0058 using Updater = Delegate<Result<void>(const GeometryContext&, TrackStateProxy,
0059 const Logger&)>;
0060
0061
0062 using OutlierFinder = Delegate<bool(ConstTrackStateProxy)>;
0063
0064
0065 using ReverseFilteringLogic = Delegate<bool(ConstTrackStateProxy)>;
0066
0067
0068 using Smoother = Delegate<Result<void>(const GeometryContext&, traj_t&,
0069 std::size_t, const Logger&)>;
0070
0071
0072 SourceLinkSurfaceAccessor surfaceAccessor;
0073
0074
0075
0076
0077 Calibrator calibrator;
0078
0079
0080 Updater updater;
0081
0082
0083
0084 OutlierFinder outlierFinder;
0085
0086
0087
0088 ReverseFilteringLogic reverseFilteringLogic;
0089
0090
0091 Smoother smoother;
0092
0093
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
0106
0107
0108 template <typename traj_t>
0109 struct KalmanFitterOptions {
0110
0111
0112
0113
0114
0115
0116
0117
0118
0119
0120
0121
0122
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
0146 std::reference_wrapper<const GeometryContext> geoContext;
0147
0148 std::reference_wrapper<const MagneticFieldContext> magFieldContext;
0149
0150 std::reference_wrapper<const CalibrationContext> calibrationContext;
0151
0152
0153 KalmanFitterExtensions<traj_t> extensions;
0154
0155
0156 PropagatorPlainOptions propagatorPlainOptions;
0157
0158
0159 const Surface* referenceSurface = nullptr;
0160
0161
0162 TrackExtrapolationStrategy referenceSurfaceStrategy =
0163 TrackExtrapolationStrategy::firstOrLast;
0164
0165
0166 bool multipleScattering = true;
0167
0168
0169 bool energyLoss = true;
0170
0171
0172
0173 bool reverseFiltering = false;
0174
0175
0176
0177
0178
0179
0180
0181 double reverseFilteringCovarianceScaling = 100.0;
0182
0183
0184
0185 FreeToBoundCorrection freeToBoundCorrection;
0186 };
0187
0188
0189 template <typename traj_t>
0190 struct KalmanFitterResult {
0191
0192 traj_t* fittedStates{nullptr};
0193
0194
0195
0196
0197
0198 std::size_t lastMeasurementIndex = kTrackIndexInvalid;
0199
0200
0201
0202
0203
0204 std::size_t lastTrackIndex = kTrackIndexInvalid;
0205
0206
0207 std::optional<BoundTrackParameters> fittedParameters;
0208
0209
0210 std::size_t measurementStates = 0;
0211
0212
0213
0214
0215
0216 std::size_t measurementHoles = 0;
0217
0218
0219 std::size_t processedStates = 0;
0220
0221
0222 bool finished = false;
0223
0224
0225 std::vector<const Surface*> missedActiveSurfaces;
0226
0227
0228 PathLimitReached pathLimitReached;
0229 };
0230
0231
0232
0233
0234
0235
0236
0237
0238
0239
0240
0241
0242
0243
0244
0245
0246
0247
0248
0249
0250
0251
0252 template <typename propagator_t, typename traj_t>
0253 class KalmanFitter {
0254
0255 using KalmanNavigator = typename propagator_t::Navigator;
0256
0257
0258 static constexpr bool isDirectNavigator =
0259 std::is_same_v<KalmanNavigator, DirectNavigator>;
0260
0261 public:
0262
0263
0264
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
0274 propagator_t m_propagator;
0275
0276
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
0283
0284
0285
0286 class Actor {
0287 public:
0288
0289 using result_type = KalmanFitterResult<traj_t>;
0290
0291
0292 SurfaceReached targetReached{std::numeric_limits<double>::lowest()};
0293
0294
0295 std::unordered_map<const Surface*, SourceLink> inputMeasurements;
0296
0297
0298 bool multipleScattering = true;
0299
0300
0301 bool energyLoss = true;
0302
0303
0304
0305 FreeToBoundCorrection freeToBoundCorrection;
0306
0307
0308 std::shared_ptr<traj_t> outputStates;
0309
0310 KalmanFitterExtensions<traj_t> extensions;
0311
0312
0313 const CalibrationContext* calibrationContext{nullptr};
0314
0315
0316 EndOfWorldReached endOfWorldReached;
0317
0318
0319 VolumeConstraintAborter volumeConstraintAborter;
0320
0321
0322 const Logger* actorLogger{nullptr};
0323
0324
0325 const Logger& logger() const { return *actorLogger; }
0326
0327
0328
0329
0330
0331
0332
0333
0334
0335
0336
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& ) 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
0355 if (result.pathLimitReached.internalLimit ==
0356 std::numeric_limits<double>::max()) {
0357 detail::setupLoopProtection(state, stepper, result.pathLimitReached,
0358 true, logger());
0359 }
0360
0361
0362
0363 const Surface* surface = navigator.currentSurface(state.navigation);
0364 if (surface != nullptr) {
0365
0366
0367
0368
0369
0370
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
0382
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
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& , const stepper_t& ,
0427 const navigator_t& , const result_type& result,
0428 const Logger& ) const {
0429 return result.finished;
0430 }
0431
0432
0433
0434
0435
0436
0437
0438
0439
0440
0441
0442
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
0453 const auto sourceLinkIt = inputMeasurements.find(&surface);
0454 if (sourceLinkIt != inputMeasurements.end()) {
0455
0456 ACTS_VERBOSE("Measurement surface " << surface.geometryId()
0457 << " detected.");
0458
0459 stepper.transportCovarianceToBound(state.stepping, surface,
0460 freeToBoundCorrection);
0461
0462
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
0477
0478
0479
0480
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
0491
0492 trackStateProxy.setReferenceSurface(surface.getSharedPtr());
0493
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
0504 trackStateProxy.predicted() = boundParams.parameters();
0505 trackStateProxy.predictedCovariance() = state.stepping.cov;
0506
0507 trackStateProxy.jacobian() = jacobian;
0508 trackStateProxy.pathLength() = pathLength;
0509
0510
0511
0512 extensions.calibrator(state.geoContext, *calibrationContext,
0513 sourceLinkIt->second, trackStateProxy);
0514
0515
0516 auto typeFlags = trackStateProxy.typeFlags();
0517 typeFlags.setHasParameters();
0518 if (surface.hasMaterial()) {
0519 typeFlags.setHasMaterial();
0520 }
0521
0522
0523
0524
0525
0526
0527
0528 if (!extensions.outlierFinder(trackStateProxyConst)) {
0529
0530 trackStateProxy.addComponents(TrackStatePropMask::Filtered);
0531
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
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
0545 typeFlags.setIsOutlier();
0546 trackStateProxy.shareFrom(trackStateProxy,
0547 TrackStatePropMask::Predicted,
0548 TrackStatePropMask::Filtered);
0549 }
0550
0551 result.lastTrackIndex = trackStateProxy.index();
0552
0553
0554 if (trackStateProxy.typeFlags().isMeasurement()) {
0555
0556 ACTS_VERBOSE("Filtering step successful, updated parameters are:\n"
0557 << trackStateProxy.filtered().transpose());
0558
0559 stepper.update(state.stepping,
0560 MultiTrajectoryHelpers::freeFiltered(
0561 state.options.geoContext, trackStateProxy),
0562 trackStateProxy.filtered(),
0563 trackStateProxy.filteredCovariance(), surface);
0564
0565 ++result.measurementStates;
0566 }
0567
0568
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
0582 ++result.processedStates;
0583
0584
0585 result.measurementHoles = result.missedActiveSurfaces.size();
0586
0587
0588 result.lastMeasurementIndex = result.lastTrackIndex;
0589
0590 } else if ((precedingMeasurementExists && surfaceIsSensitive) ||
0591 surfaceHasMaterial) {
0592
0593
0594
0595
0596
0597 TrackStatePropMask mask =
0598 TrackStatePropMask::Predicted | TrackStatePropMask::Jacobian;
0599 typename traj_t::TrackStateProxy trackStateProxy =
0600 result.fittedStates->makeTrackState(mask, result.lastTrackIndex);
0601
0602
0603
0604 trackStateProxy.setReferenceSurface(surface.getSharedPtr());
0605
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
0614 trackStateProxy.predicted() = boundParams.parameters();
0615 trackStateProxy.predictedCovariance() = state.stepping.cov;
0616
0617 trackStateProxy.jacobian() = jacobian;
0618 trackStateProxy.pathLength() = pathLength;
0619
0620
0621
0622 trackStateProxy.shareFrom(trackStateProxy,
0623 TrackStatePropMask::Predicted,
0624 TrackStatePropMask::Filtered);
0625
0626
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
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
0650 result.missedActiveSurfaces.push_back(&surface);
0651 }
0652
0653 ++result.processedStates;
0654
0655
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
0676
0677
0678
0679
0680
0681
0682
0683
0684
0685
0686
0687
0688
0689
0690
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
0702
0703
0704
0705
0706
0707
0708
0709
0710
0711
0712
0713
0714
0715
0716
0717
0718
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
0748
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
0758 PropagatorOptions propagatorOptions(kfOptions.geoContext,
0759 kfOptions.magFieldContext);
0760
0761
0762 propagatorOptions.setPlainOptions(kfOptions.propagatorPlainOptions);
0763
0764 if (reverseDirection) {
0765 propagatorOptions.direction = propagatorOptions.direction.invert();
0766 }
0767
0768 if constexpr (!isDirectNavigator) {
0769
0770
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
0778 propagatorOptions.navigation.externalSurfaces = *sSequence;
0779 }
0780
0781
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
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
0824
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
0847
0848
0849
0850
0851
0852
0853
0854
0855
0856
0857
0858
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 }