Core/TrackFitting/FitterTestsCommon.hpp

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