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