TrackFitting/src/KalmanFitterFunction.cpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 #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/Navigator.hpp"
0020 #include "Acts/Propagator/Propagator.hpp"
0021 #include "Acts/Propagator/SympyStepper.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::SympyStepper;
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 = std::abs(1 / trackState.filtered()[Acts::eBoundQOverP]);
0066     return (momentum <= momentumThreshold);
0067   }
0068 };
0069
0070 struct SimpleOutlierFinder {
0071   double chi2Cut = std::numeric_limits<double>::infinity();
0072
0073   bool isOutlier(
0074       Acts::VectorMultiTrajectory::ConstTrackStateProxy trackState) const {
0075     double chi2 = Acts::calculatePredictedChi2(trackState);
0076     return chi2 > chi2Cut;
0077   }
0078 };
0079
0080 using namespace ActsExamples;
0081
0082 struct KalmanFitterFunctionImpl final : public TrackFitterFunction {
0083   Fitter fitter;
0084   DirectFitter directFitter;
0085
0086   Acts::GainMatrixUpdater kfUpdater;
0087   Acts::GainMatrixSmoother kfSmoother;
0088   SimpleReverseFilteringLogic reverseFilteringLogic;
0089   SimpleOutlierFinder outlierFinder;
0090
0091   bool multipleScattering = false;
0092   bool energyLoss = false;
0093   Acts::FreeToBoundCorrection freeToBoundCorrection;
0094
0095   IndexSourceLink::SurfaceAccessor slSurfaceAccessor;
0096
0097   KalmanFitterFunctionImpl(Fitter&& f, DirectFitter&& df,
0098                            const Acts::TrackingGeometry& trkGeo)
0099       : fitter(std::move(f)),
0100         directFitter(std::move(df)),
0101         slSurfaceAccessor{trkGeo} {}
0102
0103   template <typename calibrator_t>
0104   auto makeKfOptions(const GeneralFitterOptions& options,
0105                      const calibrator_t& calibrator) const {
0106     Acts::KalmanFitterExtensions<Acts::VectorMultiTrajectory> extensions;
0107     extensions.updater.connect<
0108         &Acts::GainMatrixUpdater::operator()<Acts::VectorMultiTrajectory>>(
0109         &kfUpdater);
0110     extensions.smoother.connect<
0111         &Acts::GainMatrixSmoother::operator()<Acts::VectorMultiTrajectory>>(
0112         &kfSmoother);
0113     extensions.reverseFilteringLogic
0114         .connect<&SimpleReverseFilteringLogic::doBackwardFiltering>(
0115             &reverseFilteringLogic);
0116     extensions.outlierFinder.connect<&SimpleOutlierFinder::isOutlier>(
0117         &outlierFinder);
0118
0119     Acts::KalmanFitterOptions<Acts::VectorMultiTrajectory> kfOptions(
0120         options.geoContext, options.magFieldContext, options.calibrationContext,
0121         extensions, options.propOptions, &(*options.referenceSurface));
0122
0123     kfOptions.referenceSurfaceStrategy =
0124         Acts::KalmanFitterTargetSurfaceStrategy::first;
0125     kfOptions.multipleScattering = multipleScattering;
0126     kfOptions.energyLoss = energyLoss;
0127     kfOptions.freeToBoundCorrection = freeToBoundCorrection;
0128     kfOptions.extensions.calibrator.connect<&calibrator_t::calibrate>(
0129         &calibrator);
0130
0131     if (options.doRefit) {
0132       kfOptions.extensions.surfaceAccessor
0133           .connect<&RefittingCalibrator::accessSurface>();
0134     } else {
0135       kfOptions.extensions.surfaceAccessor
0136           .connect<&IndexSourceLink::SurfaceAccessor::operator()>(
0137               &slSurfaceAccessor);
0138     }
0139
0140     return kfOptions;
0141   }
0142
0143   TrackFitterResult operator()(const std::vector<Acts::SourceLink>& sourceLinks,
0144                                const TrackParameters& initialParameters,
0145                                const GeneralFitterOptions& options,
0146                                const MeasurementCalibratorAdapter& calibrator,
0147                                TrackContainer& tracks) const override {
0148     const auto kfOptions = makeKfOptions(options, calibrator);
0149     return fitter.fit(sourceLinks.begin(), sourceLinks.end(), initialParameters,
0150                       kfOptions, tracks);
0151   }
0152
0153   TrackFitterResult operator()(
0154       const std::vector<Acts::SourceLink>& sourceLinks,
0155       const TrackParameters& initialParameters,
0156       const GeneralFitterOptions& options,
0157       const RefittingCalibrator& calibrator,
0158       const std::vector<const Acts::Surface*>& surfaceSequence,
0159       TrackContainer& tracks) const override {
0160     const auto kfOptions = makeKfOptions(options, calibrator);
0161     return directFitter.fit(sourceLinks.begin(), sourceLinks.end(),
0162                             initialParameters, kfOptions, surfaceSequence,
0163                             tracks);
0164   }
0165 };
0166
0167 }  // namespace
0168
0169 std::shared_ptr<ActsExamples::TrackFitterFunction>
0170 ActsExamples::makeKalmanFitterFunction(
0171     std::shared_ptr<const Acts::TrackingGeometry> trackingGeometry,
0172     std::shared_ptr<const Acts::MagneticFieldProvider> magneticField,
0173     bool multipleScattering, bool energyLoss,
0174     double reverseFilteringMomThreshold,
0175     Acts::FreeToBoundCorrection freeToBoundCorrection, double chi2Cut,
0176     const Acts::Logger& logger) {
0177   // Stepper should be copied into the fitters
0178   const Stepper stepper(std::move(magneticField));
0179
0180   // Standard fitter
0181   const auto& geo = *trackingGeometry;
0182   Acts::Navigator::Config cfg{std::move(trackingGeometry)};
0183   cfg.resolvePassive = false;
0184   cfg.resolveMaterial = true;
0185   cfg.resolveSensitive = true;
0186   Acts::Navigator navigator(cfg, logger.cloneWithSuffix("Navigator"));
0187   Propagator propagator(stepper, std::move(navigator),
0188                         logger.cloneWithSuffix("Propagator"));
0189   Fitter trackFitter(std::move(propagator), logger.cloneWithSuffix("Fitter"));
0190
0191   // Direct fitter
0192   Acts::DirectNavigator directNavigator{
0193       logger.cloneWithSuffix("DirectNavigator")};
0194   DirectPropagator directPropagator(stepper, std::move(directNavigator),
0195                                     logger.cloneWithSuffix("DirectPropagator"));
0196   DirectFitter directTrackFitter(std::move(directPropagator),
0197                                  logger.cloneWithSuffix("DirectFitter"));
0198
0199   // build the fitter function. owns the fitter object.
0200   auto fitterFunction = std::make_shared<KalmanFitterFunctionImpl>(
0201       std::move(trackFitter), std::move(directTrackFitter), geo);
0202   fitterFunction->multipleScattering = multipleScattering;
0203   fitterFunction->energyLoss = energyLoss;
0204   fitterFunction->reverseFilteringLogic.momentumThreshold =
0205       reverseFilteringMomThreshold;
0206   fitterFunction->freeToBoundCorrection = freeToBoundCorrection;
0207   fitterFunction->outlierFinder.chi2Cut = chi2Cut;
0208
0209   return fitterFunction;
0210 }