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