Core/TrackFitting/FitterTestsCommon.hpp

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 #pragma once
0010
0011 #include <boost/test/unit_test.hpp>
0012
0013 #include "Acts/Definitions/TrackParametrization.hpp"
0014 #include "Acts/Definitions/Units.hpp"
0015 #include "Acts/EventData/MultiTrajectory.hpp"
0016 #include "Acts/EventData/ProxyAccessor.hpp"
0017 #include "Acts/EventData/SourceLink.hpp"
0018 #include "Acts/EventData/VectorMultiTrajectory.hpp"
0019 #include "Acts/EventData/VectorTrackContainer.hpp"
0020 #include "Acts/EventData/detail/TestSourceLink.hpp"
0021 #include "Acts/Geometry/TrackingGeometry.hpp"
0022 #include "Acts/MagneticField/ConstantBField.hpp"
0023 #include "Acts/MagneticField/MagneticFieldContext.hpp"
0024 #include "Acts/Propagator/Navigator.hpp"
0025 #include "Acts/Propagator/Propagator.hpp"
0026 #include "Acts/Propagator/StraightLineStepper.hpp"
0027 #include "Acts/Surfaces/CurvilinearSurface.hpp"
0028 #include "Acts/Tests/CommonHelpers/CubicTrackingGeometry.hpp"
0029 #include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
0030 #include "Acts/Tests/CommonHelpers/MeasurementsCreator.hpp"
0031 #include "Acts/TrackFitting/detail/KalmanGlobalCovariance.hpp"
0032 #include "Acts/Utilities/CalibrationContext.hpp"
0033 #include "Acts/Utilities/Logger.hpp"
0034
0035 #include <iterator>
0036
0037 using namespace Acts::UnitLiterals;
0038 using namespace Acts::Test;
0039
0040 /// Find outliers using plain distance for testing purposes.
0041 ///
0042 /// In a real setup, the outlier classification can be much more involved, e.g.
0043 /// by computing the weighted distance/ local chi2 value. Here, the purpose is
0044 /// to test that the basic principle works using simplified, synthetic data.
0045 /// Thus, the simplest possible implementation should do.
0046 struct TestOutlierFinder {
0047   double distanceMax = std::numeric_limits<double>::max();
0048
0049   /// Classify a measurement as a valid one or an outlier.
0050   ///
0051   /// @tparam track_state_t Type of the track state
0052   /// @param state The track state to classify
0053   /// @retval False if the measurement is not an outlier
0054   /// @retval True if the measurement is an outlier
0055   template <typename traj_t>
0056   bool operator()(typename traj_t::ConstTrackStateProxy state) const {
0057     // can't determine an outlier w/o a measurement or predicted parameters
0058     if (!state.hasCalibrated() || !state.hasPredicted()) {
0059       return false;
0060     }
0061     auto subspaceHelper = state.projectorSubspaceHelper();
0062     auto projector =
0063         subspaceHelper.fullProjector()
0064             .topLeftCorner(state.calibratedSize(), Acts::eBoundSize)
0065             .eval();
0066     auto residuals =
0067         (state.effectiveCalibrated() - projector * state.predicted()).eval();
0068     auto distance = residuals.norm();
0069     return (distanceMax <= distance);
0070   }
0071 };
0072
0073 /// Determine if the smoothing of a track should be done with or without reverse
0074 /// filtering
0075 struct TestReverseFilteringLogic {
0076   double momentumMax = std::numeric_limits<double>::max();
0077
0078   /// Classify a measurement as a valid one or an outlier.
0079   ///
0080   /// @param trackState The trackState of the last measurement
0081   /// @retval False if we don't use the reverse filtering for the smoothing of the track
0082   /// @retval True if we use the reverse filtering for the smoothing of the track
0083   template <typename traj_t>
0084   bool operator()(typename traj_t::ConstTrackStateProxy state) const {
0085     // can't determine an outlier w/o a measurement or predicted parameters
0086     auto momentum = std::abs(1 / state.filtered()[Acts::eBoundQOverP]);
0087     std::cout << "momentum : " << momentum << std::endl;
0088     return (momentum <= momentumMax);
0089   }
0090 };
0091
0092 // Construct a straight-line propagator.
0093 auto makeStraightPropagator(std::shared_ptr<const Acts::TrackingGeometry> geo) {
0094   Acts::Navigator::Config cfg{std::move(geo)};
0095   cfg.resolvePassive = false;
0096   cfg.resolveMaterial = true;
0097   cfg.resolveSensitive = true;
0098   Acts::Navigator navigator(
0099       cfg, Acts::getDefaultLogger("Navigator", Acts::Logging::INFO));
0100   Acts::StraightLineStepper stepper;
0101   return Acts::Propagator<Acts::StraightLineStepper, Acts::Navigator>(
0102       stepper, std::move(navigator));
0103 }
0104
0105 // Construct a propagator using a constant magnetic field along z.
0106 template <typename stepper_t>
0107 auto makeConstantFieldPropagator(
0108     std::shared_ptr<const Acts::TrackingGeometry> geo, double bz) {
0109   Acts::Navigator::Config cfg{std::move(geo)};
0110   cfg.resolvePassive = false;
0111   cfg.resolveMaterial = true;
0112   cfg.resolveSensitive = true;
0113   Acts::Navigator navigator(
0114       cfg, Acts::getDefaultLogger("Navigator", Acts::Logging::INFO));
0115   auto field =
0116       std::make_shared<Acts::ConstantBField>(Acts::Vector3(0.0, 0.0, bz));
0117   stepper_t stepper(std::move(field));
0118   return Acts::Propagator<decltype(stepper), Acts::Navigator>(
0119       std::move(stepper), std::move(navigator));
0120 }
0121
0122 // Put all this in a struct to avoid that all these objects are exposed as
0123 // global objects in the header
0124 struct FitterTester {
0125   using Rng = std::default_random_engine;
0126
0127   // Context objects
0128   Acts::GeometryContext geoCtx;
0129   Acts::MagneticFieldContext magCtx;
0130   Acts::CalibrationContext calCtx;
0131
0132   // detector geometry
0133   CubicTrackingGeometry geometryStore{geoCtx};
0134   std::shared_ptr<const Acts::TrackingGeometry> geometry = geometryStore();
0135
0136   Acts::detail::Test::TestSourceLink::SurfaceAccessor surfaceAccessor{
0137       *geometry};
0138
0139   // expected number of measurements for the given detector
0140   constexpr static std::size_t nMeasurements = 6u;
0141
0142   // detector resolutions
0143   MeasurementResolution resPixel = {MeasurementType::eLoc01, {25_um, 50_um}};
0144   MeasurementResolution resStrip0 = {MeasurementType::eLoc0, {100_um}};
0145   MeasurementResolution resStrip1 = {MeasurementType::eLoc1, {150_um}};
0146   MeasurementResolutionMap resolutions = {
0147       {Acts::GeometryIdentifier().withVolume(2), resPixel},
0148       {Acts::GeometryIdentifier().withVolume(3).withLayer(2), resStrip0},
0149       {Acts::GeometryIdentifier().withVolume(3).withLayer(4), resStrip1},
0150       {Acts::GeometryIdentifier().withVolume(3).withLayer(6), resStrip0},
0151       {Acts::GeometryIdentifier().withVolume(3).withLayer(8), resStrip1},
0152   };
0153
0154   // simulation propagator
0155   Acts::Propagator<Acts::StraightLineStepper, Acts::Navigator> simPropagator =
0156       makeStraightPropagator(geometry);
0157
0158   static std::vector<Acts::SourceLink> prepareSourceLinks(
0159       const std::vector<Acts::detail::Test::TestSourceLink>& sourceLinks) {
0160     std::vector<Acts::SourceLink> result;
0161     std::transform(sourceLinks.begin(), sourceLinks.end(),
0162                    std::back_inserter(result),
0163                    [](const auto& sl) { return Acts::SourceLink{sl}; });
0164     return result;
0165   }
0166
0167   //////////////////////////
0168   // The testing functions
0169   //////////////////////////
0170
0171   template <typename fitter_t, typename fitter_options_t, typename parameters_t>
0172   void test_ZeroFieldNoSurfaceForward(const fitter_t& fitter,
0173                                       fitter_options_t options,
0174                                       const parameters_t& start, Rng& rng,
0175                                       const bool expected_reversed,
0176                                       const bool expected_smoothed,
0177                                       const bool doDiag) const {
0178     auto measurements = createMeasurements(simPropagator, geoCtx, magCtx, start,
0179                                            resolutions, rng);
0180
0181     auto sourceLinks = prepareSourceLinks(measurements.sourceLinks);
0182     BOOST_REQUIRE_EQUAL(sourceLinks.size(), nMeasurements);
0183
0184     // this is the default option. set anyway for consistency
0185     options.referenceSurface = nullptr;
0186
0187     Acts::ConstProxyAccessor<bool> reversed{"reversed"};
0188     Acts::ConstProxyAccessor<bool> smoothed{"smoothed"};
0189
0190     auto doTest = [&](bool diag) {
0191       Acts::TrackContainer tracks{Acts::VectorTrackContainer{},
0192                                   Acts::VectorMultiTrajectory{}};
0193       if (diag) {
0194         tracks.addColumn<bool>("reversed");
0195         tracks.addColumn<bool>("smoothed");
0196
0197         BOOST_CHECK(tracks.hasColumn("reversed"));
0198         BOOST_CHECK(tracks.hasColumn("smoothed"));
0199       }
0200
0201       auto res = fitter.fit(sourceLinks.begin(), sourceLinks.end(), start,
0202                             options, tracks);
0203       BOOST_REQUIRE(res.ok());
0204
0205       const auto track = res.value();
0206       BOOST_CHECK_NE(track.tipIndex(), Acts::MultiTrajectoryTraits::kInvalid);
0207       BOOST_CHECK(!track.hasReferenceSurface());
0208       BOOST_CHECK_EQUAL(track.nMeasurements(), sourceLinks.size());
0209       BOOST_CHECK_EQUAL(track.nHoles(), 0u);
0210
0211       if (diag) {
0212         // check the output status flags
0213         BOOST_CHECK_EQUAL(reversed(track), expected_reversed);
0214         BOOST_CHECK_EQUAL(smoothed(track), expected_smoothed);
0215       }
0216     };
0217
0218     if (doDiag) {
0219       doTest(true);
0220     }  // with reversed & smoothed columns
0221     doTest(false);  // without the extra columns
0222   }
0223
0224   template <typename fitter_t, typename fitter_options_t, typename parameters_t>
0225   void test_ZeroFieldWithSurfaceForward(const fitter_t& fitter,
0226                                         fitter_options_t options,
0227                                         const parameters_t& start, Rng& rng,
0228                                         const bool expected_reversed,
0229                                         const bool expected_smoothed,
0230                                         const bool doDiag) const {
0231     auto measurements = createMeasurements(simPropagator, geoCtx, magCtx, start,
0232                                            resolutions, rng);
0233     auto sourceLinks = prepareSourceLinks(measurements.sourceLinks);
0234     BOOST_REQUIRE_EQUAL(sourceLinks.size(), nMeasurements);
0235
0236     // initial fitter options configured for backward filtering mode
0237     // backward filtering requires a reference surface
0238     options.referenceSurface = &start.referenceSurface();
0239     // this is the default option. set anyway for consistency
0240     options.propagatorPlainOptions.direction = Acts::Direction::Forward();
0241
0242     Acts::TrackContainer tracks{Acts::VectorTrackContainer{},
0243                                 Acts::VectorMultiTrajectory{}};
0244     tracks.addColumn<bool>("reversed");
0245     tracks.addColumn<bool>("smoothed");
0246
0247     auto res = fitter.fit(sourceLinks.begin(), sourceLinks.end(), start,
0248                           options, tracks);
0249     BOOST_REQUIRE(res.ok());
0250
0251     const auto& track = res.value();
0252     BOOST_CHECK_NE(track.tipIndex(), Acts::MultiTrajectoryTraits::kInvalid);
0253     BOOST_CHECK(track.hasReferenceSurface());
0254     BOOST_CHECK_EQUAL(track.nMeasurements(), sourceLinks.size());
0255     BOOST_CHECK_EQUAL(track.nHoles(), 0u);
0256
0257     BOOST_CHECK(tracks.hasColumn("reversed"));
0258     BOOST_CHECK(tracks.hasColumn("smoothed"));
0259
0260     Acts::ConstProxyAccessor<bool> reversed{"reversed"};
0261     Acts::ConstProxyAccessor<bool> smoothed{"smoothed"};
0262
0263     // check the output status flags
0264     if (doDiag) {
0265       BOOST_CHECK_EQUAL(smoothed(track), expected_smoothed);
0266       BOOST_CHECK_EQUAL(reversed(track), expected_reversed);
0267     }
0268
0269     // count the number of `smoothed` states
0270     if (expected_reversed && expected_smoothed) {
0271       std::size_t nSmoothed = 0;
0272       for (const auto ts : track.trackStatesReversed()) {
0273         nSmoothed += ts.hasSmoothed();
0274       }
0275       BOOST_CHECK_EQUAL(nSmoothed, sourceLinks.size());
0276     }
0277   }
0278
0279   template <typename fitter_t, typename fitter_options_t, typename parameters_t>
0280   void test_ZeroFieldWithSurfaceBackward(const fitter_t& fitter,
0281                                          fitter_options_t options,
0282                                          const parameters_t& start, Rng& rng,
0283                                          const bool expected_reversed,
0284                                          const bool expected_smoothed,
0285                                          const bool doDiag) const {
0286     auto measurements = createMeasurements(simPropagator, geoCtx, magCtx, start,
0287                                            resolutions, rng);
0288     auto sourceLinks = prepareSourceLinks(measurements.sourceLinks);
0289     BOOST_REQUIRE_EQUAL(sourceLinks.size(), nMeasurements);
0290
0291     // create a track near the tracker exit for outward->inward filtering
0292     Acts::Vector4 posOuter = start.fourPosition(geoCtx);
0293     posOuter[Acts::ePos0] = 3_m;
0294     Acts::BoundTrackParameters startOuter =
0295         Acts::BoundTrackParameters::createCurvilinear(
0296             posOuter, start.direction(), start.qOverP(), start.covariance(),
0297             Acts::ParticleHypothesis::pion());
0298
0299     options.referenceSurface = &startOuter.referenceSurface();
0300     options.propagatorPlainOptions.direction = Acts::Direction::Backward();
0301
0302     Acts::TrackContainer tracks{Acts::VectorTrackContainer{},
0303                                 Acts::VectorMultiTrajectory{}};
0304     tracks.addColumn<bool>("reversed");
0305     tracks.addColumn<bool>("smoothed");
0306
0307     auto res = fitter.fit(sourceLinks.begin(), sourceLinks.end(), startOuter,
0308                           options, tracks);
0309     BOOST_CHECK(res.ok());
0310
0311     const auto& track = res.value();
0312     BOOST_CHECK_NE(track.tipIndex(), Acts::MultiTrajectoryTraits::kInvalid);
0313     BOOST_CHECK(track.hasReferenceSurface());
0314     BOOST_CHECK_EQUAL(track.nMeasurements(), sourceLinks.size());
0315     BOOST_CHECK_EQUAL(track.nHoles(), 0u);
0316
0317     Acts::ConstProxyAccessor<bool> reversed{"reversed"};
0318     Acts::ConstProxyAccessor<bool> smoothed{"smoothed"};
0319     // check the output status flags
0320     if (doDiag) {
0321       BOOST_CHECK_EQUAL(smoothed(track), expected_smoothed);
0322       BOOST_CHECK_EQUAL(reversed(track), expected_reversed);
0323     }
0324
0325     // count the number of `smoothed` states
0326     if (expected_reversed && expected_smoothed) {
0327       std::size_t nSmoothed = 0;
0328       for (const auto ts : track.trackStatesReversed()) {
0329         nSmoothed += ts.hasSmoothed();
0330       }
0331       BOOST_CHECK_EQUAL(nSmoothed, sourceLinks.size());
0332     }
0333   }
0334
0335   template <typename fitter_t, typename fitter_options_t, typename parameters_t>
0336   void test_ZeroFieldWithSurfaceAtExit(const fitter_t& fitter,
0337                                        fitter_options_t options,
0338                                        const parameters_t& start, Rng& rng,
0339                                        const bool expected_reversed,
0340                                        const bool expected_smoothed,
0341                                        const bool doDiag) const {
0342     auto measurements = createMeasurements(simPropagator, geoCtx, magCtx, start,
0343                                            resolutions, rng);
0344     auto sourceLinks = prepareSourceLinks(measurements.sourceLinks);
0345     BOOST_REQUIRE_EQUAL(sourceLinks.size(), nMeasurements);
0346
0347     // create a boundless target surface near the tracker exit
0348     Acts::Vector3 center(3._m, 0., 0.);
0349     Acts::Vector3 normal(1., 0., 0.);
0350     std::shared_ptr<Acts::PlaneSurface> targetSurface =
0351         Acts::CurvilinearSurface(center, normal).planeSurface();
0352
0353     options.referenceSurface = targetSurface.get();
0354
0355     Acts::TrackContainer tracks{Acts::VectorTrackContainer{},
0356                                 Acts::VectorMultiTrajectory{}};
0357     tracks.addColumn<bool>("reversed");
0358     tracks.addColumn<bool>("smoothed");
0359
0360     auto res = fitter.fit(sourceLinks.begin(), sourceLinks.end(), start,
0361                           options, tracks);
0362     BOOST_REQUIRE(res.ok());
0363
0364     const auto& track = res.value();
0365     BOOST_CHECK_NE(track.tipIndex(), Acts::MultiTrajectoryTraits::kInvalid);
0366     BOOST_CHECK(track.hasReferenceSurface());
0367     BOOST_CHECK_EQUAL(track.nMeasurements(), sourceLinks.size());
0368     BOOST_CHECK_EQUAL(track.nHoles(), 0u);
0369
0370     Acts::ConstProxyAccessor<bool> reversed{"reversed"};
0371     Acts::ConstProxyAccessor<bool> smoothed{"smoothed"};
0372
0373     // check the output status flags
0374     if (doDiag) {
0375       BOOST_CHECK_EQUAL(smoothed(track), expected_smoothed);
0376       BOOST_CHECK_EQUAL(reversed(track), expected_reversed);
0377     }
0378   }
0379
0380   template <typename fitter_t, typename fitter_options_t, typename parameters_t>
0381   void test_ZeroFieldShuffled(const fitter_t& fitter, fitter_options_t options,
0382                               const parameters_t& start, Rng& rng,
0383                               const bool expected_reversed,
0384                               const bool expected_smoothed,
0385                               const bool doDiag) const {
0386     auto measurements = createMeasurements(simPropagator, geoCtx, magCtx, start,
0387                                            resolutions, rng);
0388     auto sourceLinks = prepareSourceLinks(measurements.sourceLinks);
0389     BOOST_REQUIRE_EQUAL(sourceLinks.size(), nMeasurements);
0390
0391     options.referenceSurface = &start.referenceSurface();
0392
0393     Acts::BoundVector parameters = Acts::BoundVector::Zero();
0394
0395     Acts::TrackContainer tracks{Acts::VectorTrackContainer{},
0396                                 Acts::VectorMultiTrajectory{}};
0397     tracks.addColumn<bool>("reversed");
0398     tracks.addColumn<bool>("smoothed");
0399
0400     Acts::ConstProxyAccessor<bool> reversed{"reversed"};
0401     Acts::ConstProxyAccessor<bool> smoothed{"smoothed"};
0402
0403     // fit w/ all hits in order
0404     {
0405       auto res = fitter.fit(sourceLinks.begin(), sourceLinks.end(), start,
0406                             options, tracks);
0407       BOOST_REQUIRE(res.ok());
0408
0409       const auto& track = res.value();
0410       BOOST_CHECK_NE(track.tipIndex(), Acts::MultiTrajectoryTraits::kInvalid);
0411       BOOST_CHECK_EQUAL(track.nMeasurements(), sourceLinks.size());
0412       BOOST_REQUIRE(track.hasReferenceSurface());
0413       parameters = track.parameters();
0414       BOOST_CHECK_EQUAL(track.nHoles(), 0u);
0415
0416       // check the output status flags
0417       if (doDiag) {
0418         BOOST_CHECK_EQUAL(smoothed(track), expected_smoothed);
0419         BOOST_CHECK_EQUAL(reversed(track), expected_reversed);
0420       }
0421     }
0422     // fit w/ all hits in random order
0423     {
0424       decltype(sourceLinks) shuffledSourceLinks = sourceLinks;
0425       std::shuffle(shuffledSourceLinks.begin(), shuffledSourceLinks.end(), rng);
0426       auto res = fitter.fit(shuffledSourceLinks.begin(),
0427                             shuffledSourceLinks.end(), start, options, tracks);
0428       BOOST_REQUIRE(res.ok());
0429
0430       const auto& track = res.value();
0431       BOOST_CHECK_NE(track.tipIndex(), Acts::MultiTrajectoryTraits::kInvalid);
0432       BOOST_REQUIRE(track.hasReferenceSurface());
0433       // check consistency w/ un-shuffled measurements
0434       CHECK_CLOSE_ABS(track.parameters(), parameters, 1e-5);
0435       BOOST_CHECK_EQUAL(track.nMeasurements(), sourceLinks.size());
0436       // check the output status flags
0437       if (doDiag) {
0438         BOOST_CHECK_EQUAL(smoothed(track), expected_smoothed);
0439         BOOST_CHECK_EQUAL(reversed(track), expected_reversed);
0440       }
0441     }
0442   }
0443
0444   template <typename fitter_t, typename fitter_options_t, typename parameters_t>
0445   void test_ZeroFieldWithHole(const fitter_t& fitter, fitter_options_t options,
0446                               const parameters_t& start, Rng& rng,
0447                               const bool expected_reversed,
0448                               const bool expected_smoothed,
0449                               const bool doDiag) const {
0450     auto measurements = createMeasurements(simPropagator, geoCtx, magCtx, start,
0451                                            resolutions, rng);
0452     auto sourceLinks = prepareSourceLinks(measurements.sourceLinks);
0453     BOOST_REQUIRE_EQUAL(sourceLinks.size(), nMeasurements);
0454
0455     Acts::TrackContainer tracks{Acts::VectorTrackContainer{},
0456                                 Acts::VectorMultiTrajectory{}};
0457     tracks.addColumn<bool>("reversed");
0458     tracks.addColumn<bool>("smoothed");
0459
0460     Acts::ConstProxyAccessor<bool> reversed{"reversed"};
0461     Acts::ConstProxyAccessor<bool> smoothed{"smoothed"};
0462
0463     // always keep the first and last measurement. leaving those in seems to not
0464     // count the respective surfaces as holes.
0465     for (std::size_t i = 1u; (i + 1u) < sourceLinks.size(); ++i) {
0466       // remove the i-th measurement
0467       auto withHole = sourceLinks;
0468       withHole.erase(std::next(withHole.begin(), i));
0469       BOOST_REQUIRE_EQUAL(withHole.size() + 1u, sourceLinks.size());
0470       BOOST_TEST_INFO("Removed measurement " << i);
0471
0472       auto res =
0473           fitter.fit(withHole.begin(), withHole.end(), start, options, tracks);
0474       BOOST_REQUIRE(res.ok());
0475
0476       const auto& track = res.value();
0477       BOOST_CHECK_NE(track.tipIndex(), Acts::MultiTrajectoryTraits::kInvalid);
0478       BOOST_REQUIRE(!track.hasReferenceSurface());
0479       BOOST_CHECK_EQUAL(track.nMeasurements(), withHole.size());
0480       // check the output status flags
0481       if (doDiag) {
0482         BOOST_CHECK_EQUAL(smoothed(track), expected_smoothed);
0483         BOOST_CHECK_EQUAL(reversed(track), expected_reversed);
0484       }
0485       BOOST_CHECK_EQUAL(track.nHoles(), 1u);
0486     }
0487     BOOST_CHECK_EQUAL(tracks.size(), sourceLinks.size() - 2);
0488   }
0489
0490   template <typename fitter_t, typename fitter_options_t, typename parameters_t>
0491   void test_ZeroFieldWithOutliers(const fitter_t& fitter,
0492                                   fitter_options_t options,
0493                                   const parameters_t& start, Rng& rng,
0494                                   const bool expected_reversed,
0495                                   const bool expected_smoothed,
0496                                   const bool doDiag) const {
0497     auto measurements = createMeasurements(simPropagator, geoCtx, magCtx, start,
0498                                            resolutions, rng);
0499     auto sourceLinks = prepareSourceLinks(measurements.sourceLinks);
0500     auto outlierSourceLinks =
0501         prepareSourceLinks(measurements.outlierSourceLinks);
0502     BOOST_REQUIRE_EQUAL(sourceLinks.size(), nMeasurements);
0503     BOOST_REQUIRE_EQUAL(outlierSourceLinks.size(), nMeasurements);
0504
0505     Acts::TrackContainer tracks{Acts::VectorTrackContainer{},
0506                                 Acts::VectorMultiTrajectory{}};
0507     tracks.addColumn<bool>("reversed");
0508     tracks.addColumn<bool>("smoothed");
0509
0510     Acts::ConstProxyAccessor<bool> reversed{"reversed"};
0511     Acts::ConstProxyAccessor<bool> smoothed{"smoothed"};
0512
0513     for (std::size_t i = 0; i < sourceLinks.size(); ++i) {
0514       // replace the i-th measurement with an outlier
0515       auto withOutlier = sourceLinks;
0516       withOutlier[i] = outlierSourceLinks[i];
0517       BOOST_REQUIRE_EQUAL(withOutlier.size(), sourceLinks.size());
0518       BOOST_TEST_INFO("Replaced measurement " << i << " with outlier");
0519
0520       auto res = fitter.fit(withOutlier.begin(), withOutlier.end(), start,
0521                             options, tracks);
0522       BOOST_REQUIRE(res.ok());
0523
0524       const auto& track = res.value();
0525       BOOST_CHECK_NE(track.tipIndex(), Acts::MultiTrajectoryTraits::kInvalid);
0526       // count the number of outliers
0527       std::size_t nOutliers = 0;
0528       for (const auto state : track.trackStatesReversed()) {
0529         nOutliers += state.typeFlags().test(Acts::TrackStateFlag::OutlierFlag);
0530       }
0531       BOOST_CHECK_EQUAL(nOutliers, 1u);
0532       BOOST_REQUIRE(!track.hasReferenceSurface());
0533       BOOST_CHECK_EQUAL(track.nMeasurements(), withOutlier.size() - 1u);
0534       // check the output status flags
0535       if (doDiag) {
0536         BOOST_CHECK_EQUAL(smoothed(track), expected_smoothed);
0537         BOOST_CHECK_EQUAL(reversed(track), expected_reversed);
0538       }
0539       BOOST_CHECK_EQUAL(track.nHoles(), 0u);
0540     }
0541     BOOST_CHECK_EQUAL(tracks.size(), sourceLinks.size());
0542   }
0543
0544   template <typename fitter_t, typename fitter_options_t, typename parameters_t>
0545   void test_ZeroFieldWithReverseFiltering(const fitter_t& fitter,
0546                                           fitter_options_t options,
0547                                           const parameters_t& start, Rng& rng,
0548                                           const bool expected_reversed,
0549                                           const bool expected_smoothed,
0550                                           const bool doDiag) const {
0551     auto measurements = createMeasurements(simPropagator, geoCtx, magCtx, start,
0552                                            resolutions, rng);
0553
0554     Acts::TrackContainer tracks{Acts::VectorTrackContainer{},
0555                                 Acts::VectorMultiTrajectory{}};
0556     tracks.addColumn<bool>("reversed");
0557     tracks.addColumn<bool>("smoothed");
0558
0559     Acts::ConstProxyAccessor<bool> reversed{"reversed"};
0560     Acts::ConstProxyAccessor<bool> smoothed{"smoothed"};
0561
0562     auto sourceLinks = prepareSourceLinks(measurements.sourceLinks);
0563
0564     const auto& outlierSourceLinks = measurements.outlierSourceLinks;
0565     BOOST_REQUIRE_EQUAL(sourceLinks.size(), nMeasurements);
0566     BOOST_REQUIRE_EQUAL(outlierSourceLinks.size(), nMeasurements);
0567
0568     // create a boundless target surface near the tracker entry
0569     Acts::Vector3 center(-3._m, 0., 0.);
0570     Acts::Vector3 normal(1., 0., 0.);
0571     std::shared_ptr<Acts::PlaneSurface> targetSurface =
0572         Acts::CurvilinearSurface(center, normal).planeSurface();
0573
0574     options.referenceSurface = targetSurface.get();
0575
0576     auto res = fitter.fit(sourceLinks.begin(), sourceLinks.end(), start,
0577                           options, tracks);
0578     BOOST_REQUIRE(res.ok());
0579     const auto& track = res.value();
0580
0581     // Track of 1 GeV with a threshold set at 0.1 GeV, reversed filtering should
0582     // not be used
0583     if (doDiag) {
0584       BOOST_CHECK_EQUAL(smoothed(track), expected_smoothed);
0585       BOOST_CHECK_EQUAL(reversed(track), expected_reversed);
0586     }
0587   }
0588
0589   // TODO this is not really Kalman fitter specific. is probably better tested
0590   // with a synthetic trajectory.
0591   template <typename fitter_t, typename fitter_options_t, typename parameters_t>
0592   void test_GlobalCovariance(const fitter_t& fitter, fitter_options_t options,
0593                              const parameters_t& start, Rng& rng) const {
0594     auto measurements = createMeasurements(simPropagator, geoCtx, magCtx, start,
0595                                            resolutions, rng);
0596     auto sourceLinks = prepareSourceLinks(measurements.sourceLinks);
0597     BOOST_REQUIRE_EQUAL(sourceLinks.size(), nMeasurements);
0598
0599     Acts::TrackContainer tracks{Acts::VectorTrackContainer{},
0600                                 Acts::VectorMultiTrajectory{}};
0601
0602     auto res = fitter.fit(sourceLinks.begin(), sourceLinks.end(), start,
0603                           options, tracks);
0604     BOOST_REQUIRE(res.ok());
0605
0606     // Calculate global track parameters covariance matrix
0607     const auto& track = res.value();
0608     auto [trackParamsCov, stateRowIndices] =
0609         Acts::detail::globalTrackParametersCovariance(
0610             tracks.trackStateContainer(), track.tipIndex());
0611     BOOST_CHECK_EQUAL(trackParamsCov.rows(),
0612                       sourceLinks.size() * Acts::eBoundSize);
0613     BOOST_CHECK_EQUAL(stateRowIndices.size(), sourceLinks.size());
0614     // Each smoothed track state will have eBoundSize rows/cols in the global
0615     // covariance. stateRowIndices is a map of the starting row/index with the
0616     // state tip as the key. Thus, the last track state (i.e. the state
0617     // corresponding track.tipIndex()) has a starting row/index =
0618     // eBoundSize * (nMeasurements - 1), i.e. 6*(6-1) = 30.
0619     BOOST_CHECK_EQUAL(stateRowIndices.at(track.tipIndex()),
0620                       Acts::eBoundSize * (nMeasurements - 1));
0621   }
0622 };