TrackFitting/src/KalmanFitterFunction.cpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2019-2021 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 http://mozilla.org/MPL/2.0/.
0008
0009 #include "Acts/Definitions/Direction.hpp"
0010 #include "Acts/Definitions/TrackParametrization.hpp"
0011 #include "Acts/EventData/MultiTrajectory.hpp"
0012 #include "Acts/EventData/TrackContainer.hpp"
0013 #include "Acts/EventData/TrackStatePropMask.hpp"
0014 #include "Acts/EventData/VectorMultiTrajectory.hpp"
0015 #include "Acts/EventData/VectorTrackContainer.hpp"
0016 #include "Acts/EventData/detail/CorrectedTransformationFreeToBound.hpp"
0017 #include "Acts/Geometry/GeometryIdentifier.hpp"
0018 #include "Acts/Propagator/DirectNavigator.hpp"
0019 #include "Acts/Propagator/EigenStepper.hpp"
0020 #include "Acts/Propagator/Navigator.hpp"
0021 #include "Acts/Propagator/Propagator.hpp"
0022 #include "Acts/TrackFitting/GainMatrixSmoother.hpp"
0023 #include "Acts/TrackFitting/GainMatrixUpdater.hpp"
0024 #include "Acts/TrackFitting/KalmanFitter.hpp"
0025 #include "Acts/Utilities/Delegate.hpp"
0026 #include "Acts/Utilities/Logger.hpp"
0027 #include "ActsExamples/EventData/IndexSourceLink.hpp"
0028 #include "ActsExamples/EventData/MeasurementCalibration.hpp"
0029 #include "ActsExamples/EventData/Track.hpp"
0030 #include "ActsExamples/TrackFitting/RefittingCalibrator.hpp"
0031 #include "ActsExamples/TrackFitting/TrackFitterFunction.hpp"
0032
0033 #include <algorithm>
0034 #include <cmath>
0035 #include <functional>
0036 #include <memory>
0037 #include <utility>
0038 #include <vector>
0039
0040 namespace Acts {
0041 class MagneticFieldProvider;
0042 class SourceLink;
0043 class Surface;
0044 class TrackingGeometry;
0045 }  // namespace Acts
0046
0047 namespace {
0048
0049 using Stepper = Acts::EigenStepper<>;
0050 using Propagator = Acts::Propagator<Stepper, Acts::Navigator>;
0051 using Fitter = Acts::KalmanFitter<Propagator, Acts::VectorMultiTrajectory>;
0052 using DirectPropagator = Acts::Propagator<Stepper, Acts::DirectNavigator>;
0053 using DirectFitter =
0054     Acts::KalmanFitter<DirectPropagator, Acts::VectorMultiTrajectory>;
0055
0056 using TrackContainer =
0057     Acts::TrackContainer<Acts::VectorTrackContainer,
0058                          Acts::VectorMultiTrajectory, std::shared_ptr>;
0059
0060 struct SimpleReverseFilteringLogic {
0061   double momentumThreshold = 0;
0062
0063   bool doBackwardFiltering(
0064       Acts::VectorMultiTrajectory::ConstTrackStateProxy trackState) const {
0065     auto momentum = fabs(1 / trackState.filtered()[Acts::eBoundQOverP]);
0066     return (momentum <= momentumThreshold);
0067   }
0068 };
0069
0070 using namespace ActsExamples;
0071
0072 struct KalmanFitterFunctionImpl final : public TrackFitterFunction {
0073   Fitter fitter;
0074   DirectFitter directFitter;
0075
0076   Acts::GainMatrixUpdater kfUpdater;
0077   Acts::GainMatrixSmoother kfSmoother;
0078   SimpleReverseFilteringLogic reverseFilteringLogic;
0079
0080   bool multipleScattering = false;
0081   bool energyLoss = false;
0082   Acts::FreeToBoundCorrection freeToBoundCorrection;
0083
0084   IndexSourceLink::SurfaceAccessor slSurfaceAccessor;
0085
0086   KalmanFitterFunctionImpl(Fitter&& f, DirectFitter&& df,
0087                            const Acts::TrackingGeometry& trkGeo)
0088       : fitter(std::move(f)),
0089         directFitter(std::move(df)),
0090         slSurfaceAccessor{trkGeo} {}
0091
0092   template <typename calibrator_t>
0093   auto makeKfOptions(const GeneralFitterOptions& options,
0094                      const calibrator_t& calibrator) const {
0095     Acts::KalmanFitterExtensions<Acts::VectorMultiTrajectory> extensions;
0096     extensions.updater.connect<
0097         &Acts::GainMatrixUpdater::operator()<Acts::VectorMultiTrajectory>>(
0098         &kfUpdater);
0099     extensions.smoother.connect<
0100         &Acts::GainMatrixSmoother::operator()<Acts::VectorMultiTrajectory>>(
0101         &kfSmoother);
0102     extensions.reverseFilteringLogic
0103         .connect<&SimpleReverseFilteringLogic::doBackwardFiltering>(
0104             &reverseFilteringLogic);
0105
0106     Acts::KalmanFitterOptions<Acts::VectorMultiTrajectory> kfOptions(
0107         options.geoContext, options.magFieldContext, options.calibrationContext,
0108         extensions, options.propOptions, &(*options.referenceSurface));
0109
0110     kfOptions.referenceSurfaceStrategy =
0111         Acts::KalmanFitterTargetSurfaceStrategy::first;
0112     kfOptions.multipleScattering = multipleScattering;
0113     kfOptions.energyLoss = energyLoss;
0114     kfOptions.freeToBoundCorrection = freeToBoundCorrection;
0115     kfOptions.extensions.calibrator.connect<&calibrator_t::calibrate>(
0116         &calibrator);
0117     kfOptions.extensions.surfaceAccessor
0118         .connect<&IndexSourceLink::SurfaceAccessor::operator()>(
0119             &slSurfaceAccessor);
0120
0121     return kfOptions;
0122   }
0123
0124   TrackFitterResult operator()(const std::vector<Acts::SourceLink>& sourceLinks,
0125                                const TrackParameters& initialParameters,
0126                                const GeneralFitterOptions& options,
0127                                const MeasurementCalibratorAdapter& calibrator,
0128                                TrackContainer& tracks) const override {
0129     const auto kfOptions = makeKfOptions(options, calibrator);
0130     return fitter.fit(sourceLinks.begin(), sourceLinks.end(), initialParameters,
0131                       kfOptions, tracks);
0132   }
0133
0134   TrackFitterResult operator()(
0135       const std::vector<Acts::SourceLink>& sourceLinks,
0136       const TrackParameters& initialParameters,
0137       const GeneralFitterOptions& options,
0138       const RefittingCalibrator& calibrator,
0139       const std::vector<const Acts::Surface*>& surfaceSequence,
0140       TrackContainer& tracks) const override {
0141     const auto kfOptions = makeKfOptions(options, calibrator);
0142     return directFitter.fit(sourceLinks.begin(), sourceLinks.end(),
0143                             initialParameters, kfOptions, surfaceSequence,
0144                             tracks);
0145   }
0146 };
0147
0148 }  // namespace
0149
0150 std::shared_ptr<ActsExamples::TrackFitterFunction>
0151 ActsExamples::makeKalmanFitterFunction(
0152     std::shared_ptr<const Acts::TrackingGeometry> trackingGeometry,
0153     std::shared_ptr<const Acts::MagneticFieldProvider> magneticField,
0154     bool multipleScattering, bool energyLoss,
0155     double reverseFilteringMomThreshold,
0156     Acts::FreeToBoundCorrection freeToBoundCorrection,
0157     const Acts::Logger& logger) {
0158   // Stepper should be copied into the fitters
0159   const Stepper stepper(std::move(magneticField));
0160
0161   // Standard fitter
0162   const auto& geo = *trackingGeometry;
0163   Acts::Navigator::Config cfg{std::move(trackingGeometry)};
0164   cfg.resolvePassive = false;
0165   cfg.resolveMaterial = true;
0166   cfg.resolveSensitive = true;
0167   Acts::Navigator navigator(cfg, logger.cloneWithSuffix("Navigator"));
0168   Propagator propagator(stepper, std::move(navigator),
0169                         logger.cloneWithSuffix("Propagator"));
0170   Fitter trackFitter(std::move(propagator), logger.cloneWithSuffix("Fitter"));
0171
0172   // Direct fitter
0173   Acts::DirectNavigator directNavigator{
0174       logger.cloneWithSuffix("DirectNavigator")};
0175   DirectPropagator directPropagator(stepper, std::move(directNavigator),
0176                                     logger.cloneWithSuffix("DirectPropagator"));
0177   DirectFitter directTrackFitter(std::move(directPropagator),
0178                                  logger.cloneWithSuffix("DirectFitter"));
0179
0180   // build the fitter function. owns the fitter object.
0181   auto fitterFunction = std::make_shared<KalmanFitterFunctionImpl>(
0182       std::move(trackFitter), std::move(directTrackFitter), geo);
0183   fitterFunction->multipleScattering = multipleScattering;
0184   fitterFunction->energyLoss = energyLoss;
0185   fitterFunction->reverseFilteringLogic.momentumThreshold =
0186       reverseFilteringMomThreshold;
0187   fitterFunction->freeToBoundCorrection = freeToBoundCorrection;
0188
0189   return fitterFunction;
0190 }