Core/Propagator/MultiStepperTests.cpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2018-2020 CERN for the benefit of the Acts project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
0008
0009 #include <boost/test/unit_test.hpp>
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Definitions/Direction.hpp"
0013 #include "Acts/Definitions/TrackParametrization.hpp"
0014 #include "Acts/EventData/Charge.hpp"
0015 #include "Acts/EventData/GenericBoundTrackParameters.hpp"
0016 #include "Acts/EventData/MultiComponentTrackParameters.hpp"
0017 #include "Acts/EventData/TrackParameters.hpp"
0018 #include "Acts/EventData/detail/CorrectedTransformationFreeToBound.hpp"
0019 #include "Acts/Geometry/GeometryContext.hpp"
0020 #include "Acts/MagneticField/ConstantBField.hpp"
0021 #include "Acts/MagneticField/MagneticFieldContext.hpp"
0022 #include "Acts/MagneticField/NullBField.hpp"
0023 #include "Acts/Propagator/DefaultExtension.hpp"
0024 #include "Acts/Propagator/EigenStepper.hpp"
0025 #include "Acts/Propagator/MultiEigenStepperLoop.hpp"
0026 #include "Acts/Propagator/Navigator.hpp"
0027 #include "Acts/Propagator/Propagator.hpp"
0028 #include "Acts/Propagator/StepperExtensionList.hpp"
0029 #include "Acts/Surfaces/PlaneSurface.hpp"
0030 #include "Acts/Surfaces/Surface.hpp"
0031 #include "Acts/Utilities/Helpers.hpp"
0032 #include "Acts/Utilities/Intersection.hpp"
0033 #include "Acts/Utilities/Logger.hpp"
0034 #include "Acts/Utilities/Result.hpp"
0035
0036 #include <algorithm>
0037 #include <array>
0038 #include <cmath>
0039 #include <cstddef>
0040 #include <memory>
0041 #include <optional>
0042 #include <stdexcept>
0043 #include <tuple>
0044 #include <type_traits>
0045 #include <utility>
0046 #include <vector>
0047
0048 namespace Acts {
0049 struct MultiStepperSurfaceReached;
0050 }  // namespace Acts
0051
0052 using namespace Acts;
0053 using namespace Acts::VectorHelpers;
0054
0055 /////////////////////////////////////////////////////
0056 // Some useful global objects, typedefs and structs
0057 /////////////////////////////////////////////////////
0058 const MagneticFieldContext magCtx;
0059 const GeometryContext geoCtx;
0060
0061 using MultiStepperLoop =
0062     MultiEigenStepperLoop<StepperExtensionList<DefaultExtension>>;
0063 using SingleStepper = EigenStepper<StepperExtensionList<DefaultExtension>>;
0064
0065 const double defaultStepSize = 123.;
0066 const auto defaultNDir = Direction::Backward;
0067
0068 const auto defaultBField =
0069     std::make_shared<ConstantBField>(Vector3(1., 2.5, 33.33));
0070 const auto defaultNullBField = std::make_shared<NullBField>();
0071
0072 const auto particleHypothesis = ParticleHypothesis::pion();
0073
0074 struct Options {
0075   double stepTolerance = 1e-4;
0076   double stepSizeCutOff = 0.0;
0077   std::size_t maxRungeKuttaStepTrials = 10;
0078   Direction direction = defaultNDir;
0079   const Acts::Logger &logger = Acts::getDummyLogger();
0080 };
0081
0082 struct MockNavigator {};
0083
0084 static constexpr MockNavigator mockNavigator;
0085
0086 struct Navigation {};
0087
0088 template <typename stepper_state_t>
0089 struct DummyPropState {
0090   stepper_state_t &stepping;
0091   Options options;
0092   Navigation navigation;
0093   GeometryContext geoContext;
0094
0095   DummyPropState(Direction direction, stepper_state_t &ss)
0096       : stepping(ss),
0097         options(Options{}),
0098         navigation(Navigation{}),
0099         geoContext(geoCtx) {
0100     options.direction = direction;
0101   }
0102 };
0103
0104 template <typename T>
0105 using components_t = typename T::components;
0106
0107 // Makes random bound parameters and covariance and a plane surface at {0,0,0}
0108 // with normal {1,0,0}. Optionally some external fixed bound parameters can be
0109 // supplied
0110 auto makeDefaultBoundPars(bool cov = true, std::size_t n = 4,
0111                           std::optional<BoundVector> ext_pars = std::nullopt) {
0112   std::vector<std::tuple<double, BoundVector, std::optional<BoundSquareMatrix>>>
0113       cmps;
0114   using Opt = std::optional<BoundSquareMatrix>;
0115
0116   auto make_random_sym_matrix = []() {
0117     auto c = BoundSquareMatrix::Random().eval();
0118     c *= c.transpose();
0119     return c;
0120   };
0121
0122   for (auto i = 0ul; i < n; ++i) {
0123     cmps.push_back({1. / n, ext_pars ? *ext_pars : BoundVector::Random(),
0124                     cov ? Opt{make_random_sym_matrix()} : Opt{}});
0125   }
0126
0127   auto surface = Acts::Surface::makeShared<Acts::PlaneSurface>(
0128       Vector3::Zero(), Vector3{1., 0., 0.});
0129
0130   return MultiComponentBoundTrackParameters(surface, cmps, particleHypothesis);
0131 }
0132
0133 //////////////////////////////////////////////////////
0134 /// Test the construction of the MultiStepper::State
0135 //////////////////////////////////////////////////////
0136 template <typename multi_stepper_t, bool Cov>
0137 void test_multi_stepper_state() {
0138   using MultiState = typename multi_stepper_t::State;
0139   using MultiStepper = multi_stepper_t;
0140
0141   constexpr std::size_t N = 4;
0142   const auto multi_pars = makeDefaultBoundPars(Cov, N, BoundVector::Ones());
0143
0144   MultiState state(geoCtx, magCtx, defaultBField, multi_pars, defaultStepSize);
0145
0146   MultiStepper ms(defaultBField);
0147
0148   BOOST_CHECK_EQUAL(N, ms.numberComponents(state));
0149
0150   // Test the result & compare with the input/test for reasonable members
0151   auto const_iterable = ms.constComponentIterable(state);
0152   for (const auto cmp : const_iterable) {
0153     BOOST_CHECK_EQUAL(cmp.jacTransport(), FreeMatrix::Identity());
0154     BOOST_CHECK_EQUAL(cmp.derivative(), FreeVector::Zero());
0155     if constexpr (!Cov) {
0156       BOOST_CHECK_EQUAL(cmp.jacToGlobal(), BoundToFreeMatrix::Zero());
0157       BOOST_CHECK_EQUAL(cmp.cov(), BoundSquareMatrix::Zero());
0158     }
0159   }
0160
0161   BOOST_CHECK_EQUAL(state.pathAccumulated, 0.);
0162   for (const auto cmp : const_iterable) {
0163     BOOST_CHECK_EQUAL(cmp.pathAccumulated(), 0.);
0164   }
0165
0166   // covTransport in the MultiEigenStepperLoop is redundant and
0167   // thus not part of the interface. However, we want to check them for
0168   // consistency.
0169   if constexpr (Acts::Concepts::exists<components_t, MultiState>) {
0170     BOOST_CHECK(!state.covTransport);
0171     for (const auto &cmp : state.components) {
0172       BOOST_CHECK_EQUAL(cmp.state.covTransport, Cov);
0173     }
0174   }
0175 }
0176
0177 BOOST_AUTO_TEST_CASE(multi_stepper_state_no_cov) {
0178   test_multi_stepper_state<MultiStepperLoop, false>();
0179 }
0180
0181 template <typename multi_stepper_t>
0182 void test_multi_stepper_state_invalid() {
0183   using MultiState = typename multi_stepper_t::State;
0184
0185   // Empty component vector
0186   const auto multi_pars = makeDefaultBoundPars(false, 0);
0187
0188   BOOST_CHECK_THROW(
0189       MultiState(geoCtx, magCtx, defaultBField, multi_pars, defaultStepSize),
0190       std::invalid_argument);
0191 }
0192
0193 BOOST_AUTO_TEST_CASE(multi_eigen_stepper_state_invalid) {
0194   test_multi_stepper_state_invalid<MultiStepperLoop>();
0195 }
0196
0197 ////////////////////////////////////////////////////////////////////////
0198 // Compare the Multi-Stepper against the Eigen-Stepper for consistency
0199 ////////////////////////////////////////////////////////////////////////
0200 template <typename multi_stepper_t>
0201 void test_multi_stepper_vs_eigen_stepper() {
0202   using MultiState = typename multi_stepper_t::State;
0203   using MultiStepper = multi_stepper_t;
0204
0205   const BoundVector pars = BoundVector::Ones();
0206   const BoundSquareMatrix cov = BoundSquareMatrix::Identity();
0207
0208   std::vector<std::tuple<double, BoundVector, std::optional<BoundSquareMatrix>>>
0209       cmps(4, {0.25, pars, cov});
0210
0211   auto surface = Acts::Surface::makeShared<Acts::PlaneSurface>(
0212       Vector3::Zero(), Vector3::Ones().normalized());
0213
0214   MultiComponentBoundTrackParameters multi_pars(surface, cmps,
0215                                                 particleHypothesis);
0216   BoundTrackParameters single_pars(surface, pars, cov, particleHypothesis);
0217
0218   MultiState multi_state(geoCtx, magCtx, defaultBField, multi_pars,
0219                          defaultStepSize);
0220   SingleStepper::State single_state(geoCtx, defaultBField->makeCache(magCtx),
0221                                     single_pars, defaultStepSize);
0222
0223   MultiStepper multi_stepper(defaultBField);
0224   SingleStepper single_stepper(defaultBField);
0225
0226   for (auto cmp : multi_stepper.componentIterable(multi_state)) {
0227     cmp.status() = Acts::Intersection3D::Status::reachable;
0228   }
0229
0230   // Do some steps and check that the results match
0231   for (int i = 0; i < 10; ++i) {
0232     // Single stepper
0233     auto single_prop_state = DummyPropState(defaultNDir, single_state);
0234     auto single_result = single_stepper.step(single_prop_state, mockNavigator);
0235     single_stepper.transportCovarianceToCurvilinear(single_state);
0236
0237     // Multi stepper;
0238     auto multi_prop_state = DummyPropState(defaultNDir, multi_state);
0239     auto multi_result = multi_stepper.step(multi_prop_state, mockNavigator);
0240     multi_stepper.transportCovarianceToCurvilinear(multi_state);
0241
0242     // Check equality
0243     BOOST_REQUIRE(multi_result.ok());
0244     BOOST_REQUIRE_EQUAL(multi_result.ok(), single_result.ok());
0245
0246     BOOST_CHECK_EQUAL(*single_result, *multi_result);
0247
0248     for (const auto cmp : multi_stepper.constComponentIterable(multi_state)) {
0249       BOOST_CHECK_EQUAL(cmp.pars(), single_state.pars);
0250       BOOST_CHECK_EQUAL(cmp.cov(), single_state.cov);
0251       BOOST_CHECK_EQUAL(cmp.jacTransport(), single_state.jacTransport);
0252       BOOST_CHECK_EQUAL(cmp.jacToGlobal(), single_state.jacToGlobal);
0253       BOOST_CHECK_EQUAL(cmp.derivative(), single_state.derivative);
0254       BOOST_CHECK_EQUAL(cmp.pathAccumulated(), single_state.pathAccumulated);
0255     }
0256   }
0257 }
0258
0259 BOOST_AUTO_TEST_CASE(multi_eigen_vs_single_eigen) {
0260   test_multi_stepper_vs_eigen_stepper<MultiStepperLoop>();
0261 }
0262
0263 /////////////////////////////
0264 // Test stepsize accessors
0265 /////////////////////////////
0266
0267 // TODO do this later, when we introduce the MultiEigenStepperSIMD, which there
0268 // needs new interfaces...
0269
0270 ////////////////////////////////////////////////////
0271 // Test the modifying accessors to the components
0272 ////////////////////////////////////////////////////
0273 template <typename multi_stepper_t>
0274 void test_components_modifying_accessors() {
0275   using MultiState = typename multi_stepper_t::State;
0276   using MultiStepper = multi_stepper_t;
0277
0278   const auto multi_pars = makeDefaultBoundPars();
0279
0280   MultiState mutable_multi_state(geoCtx, magCtx, defaultBField, multi_pars,
0281                                  defaultStepSize);
0282   const MultiState const_multi_state(geoCtx, magCtx, defaultBField, multi_pars,
0283                                      defaultStepSize);
0284
0285   MultiStepper multi_stepper(defaultBField);
0286
0287   auto modify = [&](const auto &projector) {
0288     // Here test the mutable overloads of the mutable iterable
0289     for (auto cmp : multi_stepper.componentIterable(mutable_multi_state)) {
0290       using type = std::decay_t<decltype(projector(cmp))>;
0291       if constexpr (std::is_enum_v<type>) {
0292         projector(cmp) =
0293             static_cast<type>(static_cast<int>(projector(cmp)) + 1);
0294       } else {
0295         projector(cmp) *= 2.0;
0296       }
0297     }
0298   };
0299
0300   auto check = [&](const auto &projector) {
0301     // Here test the const-member functions of the mutable iterable
0302     auto mutable_state_iterable =
0303         multi_stepper.componentIterable(mutable_multi_state);
0304     // Here test the const iterable
0305     auto const_state_iterable =
0306         multi_stepper.constComponentIterable(const_multi_state);
0307
0308     auto mstate_it = mutable_state_iterable.begin();
0309     auto cstate_it = const_state_iterable.begin();
0310     for (; cstate_it != const_state_iterable.end(); ++mstate_it, ++cstate_it) {
0311       const auto mstate_cmp = *mstate_it;
0312       auto cstate_cmp = *cstate_it;
0313
0314       using type = std::decay_t<decltype(projector(mstate_cmp))>;
0315
0316       if constexpr (std::is_arithmetic_v<type>) {
0317         BOOST_CHECK_CLOSE(projector(mstate_cmp), 2.0 * projector(cstate_cmp),
0318                           1.e-8);
0319       } else if constexpr (std::is_enum_v<type>) {
0320         BOOST_CHECK_EQUAL(static_cast<int>(projector(mstate_cmp)),
0321                           1 + static_cast<int>(projector(cstate_cmp)));
0322       } else {
0323         BOOST_CHECK(
0324             projector(mstate_cmp).isApprox(2.0 * projector(cstate_cmp), 1.e-8));
0325       }
0326     }
0327   };
0328
0329   const auto projectors = std::make_tuple(
0330       [](auto &cmp) -> decltype(auto) { return cmp.status(); },
0331       [](auto &cmp) -> decltype(auto) { return cmp.pathAccumulated(); },
0332       [](auto &cmp) -> decltype(auto) { return cmp.weight(); },
0333       [](auto &cmp) -> decltype(auto) { return cmp.pars(); },
0334       [](auto &cmp) -> decltype(auto) { return cmp.cov(); },
0335       [](auto &cmp) -> decltype(auto) { return cmp.jacTransport(); },
0336       [](auto &cmp) -> decltype(auto) { return cmp.derivative(); },
0337       [](auto &cmp) -> decltype(auto) { return cmp.jacobian(); },
0338       [](auto &cmp) -> decltype(auto) { return cmp.jacToGlobal(); });
0339
0340   std::apply(
0341       [&](const auto &...projs) {
0342         // clang-format off
0343         ( [&]() { modify(projs); check(projs); }(), ...);
0344         // clang-format on
0345       },
0346       projectors);
0347 }
0348
0349 BOOST_AUTO_TEST_CASE(multi_eigen_component_iterable_with_modification) {
0350   test_components_modifying_accessors<MultiStepperLoop>();
0351 }
0352
0353 /////////////////////////////////////////////
0354 // Test if the surface status update works
0355 /////////////////////////////////////////////
0356 template <typename multi_stepper_t>
0357 void test_multi_stepper_surface_status_update() {
0358   using MultiState = typename multi_stepper_t::State;
0359   using MultiStepper = multi_stepper_t;
0360
0361   auto start_surface = Acts::Surface::makeShared<Acts::PlaneSurface>(
0362       Vector3::Zero(), Vector3{1.0, 0.0, 0.0});
0363
0364   auto right_surface = Acts::Surface::makeShared<Acts::PlaneSurface>(
0365       Vector3{1.0, 0.0, 0.0}, Vector3{1.0, 0.0, 0.0});
0366
0367   std::vector<std::tuple<double, BoundVector, std::optional<BoundSquareMatrix>>>
0368       cmps(2, {0.5, BoundVector::Zero(), std::nullopt});
0369   std::get<BoundVector>(cmps[0])[eBoundTheta] = M_PI_2;
0370   std::get<BoundVector>(cmps[1])[eBoundTheta] = -M_PI_2;
0371   std::get<BoundVector>(cmps[0])[eBoundQOverP] = 1.0;
0372   std::get<BoundVector>(cmps[1])[eBoundQOverP] = 1.0;
0373
0374   MultiComponentBoundTrackParameters multi_pars(start_surface, cmps,
0375                                                 particleHypothesis);
0376
0377   BOOST_REQUIRE(std::get<1>(multi_pars[0])
0378                     .direction()
0379                     .isApprox(Vector3{1.0, 0.0, 0.0}, 1.e-10));
0380   BOOST_REQUIRE(std::get<1>(multi_pars[1])
0381                     .direction()
0382                     .isApprox(Vector3{-1.0, 0.0, 0.0}, 1.e-10));
0383
0384   MultiState multi_state(geoCtx, magCtx, defaultNullBField, multi_pars,
0385                          defaultStepSize);
0386   SingleStepper::State single_state(
0387       geoCtx, defaultNullBField->makeCache(magCtx), std::get<1>(multi_pars[0]),
0388       defaultStepSize);
0389
0390   MultiStepper multi_stepper(defaultNullBField);
0391   SingleStepper single_stepper(defaultNullBField);
0392
0393   // Update surface status and check
0394   {
0395     auto status = multi_stepper.updateSurfaceStatus(multi_state, *right_surface,
0396                                                     0, Direction::Forward,
0397                                                     BoundaryCheck(false));
0398
0399     BOOST_CHECK_EQUAL(status, Intersection3D::Status::reachable);
0400
0401     auto cmp_iterable = multi_stepper.constComponentIterable(multi_state);
0402
0403     BOOST_CHECK_EQUAL((*cmp_iterable.begin()).status(),
0404                       Intersection3D::Status::reachable);
0405     BOOST_CHECK_EQUAL((*(++cmp_iterable.begin())).status(),
0406                       Intersection3D::Status::missed);
0407   }
0408
0409   // Step forward now
0410   {
0411     auto multi_prop_state = DummyPropState(Direction::Forward, multi_state);
0412     multi_stepper.step(multi_prop_state, mockNavigator);
0413
0414     // Single stepper
0415     auto single_prop_state = DummyPropState(Direction::Forward, single_state);
0416     single_stepper.step(single_prop_state, mockNavigator);
0417   }
0418
0419   // Update surface status and check again
0420   {
0421     auto status = multi_stepper.updateSurfaceStatus(multi_state, *right_surface,
0422                                                     0, Direction::Forward,
0423                                                     BoundaryCheck(false));
0424
0425     BOOST_CHECK_EQUAL(status, Intersection3D::Status::onSurface);
0426
0427     auto cmp_iterable = multi_stepper.constComponentIterable(multi_state);
0428
0429     BOOST_CHECK_EQUAL((*cmp_iterable.begin()).status(),
0430                       Intersection3D::Status::onSurface);
0431     BOOST_CHECK_EQUAL((*(++cmp_iterable.begin())).status(),
0432                       Intersection3D::Status::missed);
0433   }
0434
0435   // Start surface should be unreachable
0436   {
0437     auto status = multi_stepper.updateSurfaceStatus(multi_state, *start_surface,
0438                                                     0, Direction::Forward,
0439                                                     BoundaryCheck(false));
0440
0441     BOOST_CHECK_EQUAL(status, Intersection3D::Status::unreachable);
0442
0443     auto cmp_iterable = multi_stepper.constComponentIterable(multi_state);
0444
0445     BOOST_CHECK_EQUAL((*cmp_iterable.begin()).status(),
0446                       Intersection3D::Status::unreachable);
0447     BOOST_CHECK_EQUAL((*(++cmp_iterable.begin())).status(),
0448                       Intersection3D::Status::unreachable);
0449   }
0450 }
0451
0452 BOOST_AUTO_TEST_CASE(test_surface_status_and_cmpwise_bound_state) {
0453   test_multi_stepper_surface_status_update<MultiStepperLoop>();
0454 }
0455
0456 //////////////////////////////////
0457 // Test Bound state computations
0458 //////////////////////////////////
0459 template <typename multi_stepper_t>
0460 void test_component_bound_state() {
0461   using MultiState = typename multi_stepper_t::State;
0462   using MultiStepper = multi_stepper_t;
0463
0464   auto start_surface = Acts::Surface::makeShared<Acts::PlaneSurface>(
0465       Vector3::Zero(), Vector3{1.0, 0.0, 0.0});
0466
0467   auto right_surface = Acts::Surface::makeShared<Acts::PlaneSurface>(
0468       Vector3{1.0, 0.0, 0.0}, Vector3{1.0, 0.0, 0.0});
0469
0470   std::vector<std::tuple<double, BoundVector, std::optional<BoundSquareMatrix>>>
0471       cmps(2, {0.5, BoundVector::Zero(), std::nullopt});
0472   std::get<BoundVector>(cmps[0])[eBoundTheta] = M_PI_2;
0473   std::get<BoundVector>(cmps[1])[eBoundTheta] = -M_PI_2;
0474   std::get<BoundVector>(cmps[0])[eBoundQOverP] = 1.0;
0475   std::get<BoundVector>(cmps[1])[eBoundQOverP] = 1.0;
0476
0477   MultiComponentBoundTrackParameters multi_pars(start_surface, cmps,
0478                                                 particleHypothesis);
0479
0480   BOOST_REQUIRE(std::get<1>(multi_pars[0])
0481                     .direction()
0482                     .isApprox(Vector3{1.0, 0.0, 0.0}, 1.e-10));
0483   BOOST_REQUIRE(std::get<1>(multi_pars[1])
0484                     .direction()
0485                     .isApprox(Vector3{-1.0, 0.0, 0.0}, 1.e-10));
0486
0487   MultiState multi_state(geoCtx, magCtx, defaultNullBField, multi_pars,
0488                          defaultStepSize);
0489   SingleStepper::State single_state(
0490       geoCtx, defaultNullBField->makeCache(magCtx), std::get<1>(multi_pars[0]),
0491       defaultStepSize);
0492
0493   MultiStepper multi_stepper(defaultNullBField);
0494   SingleStepper single_stepper(defaultNullBField);
0495
0496   // Step forward now
0497   {
0498     multi_stepper.updateSurfaceStatus(multi_state, *right_surface, 0,
0499                                       Direction::Forward, BoundaryCheck(false));
0500     auto multi_prop_state = DummyPropState(Direction::Forward, multi_state);
0501     multi_stepper.step(multi_prop_state, mockNavigator);
0502
0503     // Single stepper
0504     single_stepper.updateSurfaceStatus(single_state, *right_surface, 0,
0505                                        Direction::Forward,
0506                                        BoundaryCheck(false));
0507     auto single_prop_state = DummyPropState(Direction::Forward, single_state);
0508     single_stepper.step(single_prop_state, mockNavigator);
0509   }
0510
0511   // Check component-wise bound-state
0512   {
0513     auto single_bound_state = single_stepper.boundState(
0514         single_state, *right_surface, true, FreeToBoundCorrection(false));
0515     BOOST_REQUIRE(single_bound_state.ok());
0516
0517     auto cmp_iterable = multi_stepper.componentIterable(multi_state);
0518
0519     auto ok_bound_state =
0520         (*cmp_iterable.begin())
0521             .boundState(*right_surface, true, FreeToBoundCorrection(false));
0522     BOOST_REQUIRE(ok_bound_state.ok());
0523     BOOST_CHECK(*single_bound_state == *ok_bound_state);
0524
0525     auto failed_bound_state =
0526         (*(++cmp_iterable.begin()))
0527             .boundState(*right_surface, true, FreeToBoundCorrection(false));
0528     BOOST_CHECK(!failed_bound_state.ok());
0529   }
0530 }
0531
0532 BOOST_AUTO_TEST_CASE(test_component_wise_bound_state) {
0533   test_component_bound_state<MultiStepperLoop>();
0534 }
0535
0536 template <typename multi_stepper_t>
0537 void test_combined_bound_state_function() {
0538   using MultiState = typename multi_stepper_t::State;
0539   using MultiStepper = multi_stepper_t;
0540
0541   auto surface = Acts::Surface::makeShared<Acts::PlaneSurface>(
0542       Vector3::Zero(), Vector3{1.0, 0.0, 0.0});
0543
0544   // Use Ones() here, so that the angles are in correct range
0545   const auto pars = BoundVector::Ones().eval();
0546   const auto cov = []() {
0547     auto c = BoundSquareMatrix::Random().eval();
0548     c *= c.transpose();
0549     return c;
0550   }();
0551
0552   std::vector<std::tuple<double, BoundVector, std::optional<BoundSquareMatrix>>>
0553       cmps(4, {0.25, pars, cov});
0554
0555   MultiComponentBoundTrackParameters multi_pars(surface, cmps,
0556                                                 particleHypothesis);
0557   MultiState multi_state(geoCtx, magCtx, defaultBField, multi_pars,
0558                          defaultStepSize);
0559   MultiStepper multi_stepper(defaultBField);
0560
0561   auto res = multi_stepper.boundState(multi_state, *surface, true,
0562                                       FreeToBoundCorrection(false));
0563
0564   BOOST_REQUIRE(res.ok());
0565
0566   const auto [bound_pars, jacobian, pathLength] = *res;
0567
0568   BOOST_CHECK_EQUAL(jacobian, decltype(jacobian)::Zero());
0569   BOOST_CHECK_EQUAL(pathLength, 0.0);
0570   BOOST_CHECK(bound_pars.parameters().isApprox(pars, 1.e-8));
0571   BOOST_CHECK(bound_pars.covariance()->isApprox(cov, 1.e-8));
0572 }
0573
0574 BOOST_AUTO_TEST_CASE(test_combined_bound_state) {
0575   test_combined_bound_state_function<MultiStepperLoop>();
0576 }
0577
0578 //////////////////////////////////////////////////
0579 // Test the combined curvilinear state function
0580 //////////////////////////////////////////////////
0581 template <typename multi_stepper_t>
0582 void test_combined_curvilinear_state_function() {
0583   using MultiState = typename multi_stepper_t::State;
0584   using MultiStepper = multi_stepper_t;
0585
0586   auto surface = Acts::Surface::makeShared<Acts::PlaneSurface>(
0587       Vector3::Zero(), Vector3{1.0, 0.0, 0.0});
0588
0589   // Use Ones() here, so that the angles are in correct range
0590   const auto pars = BoundVector::Ones().eval();
0591   const auto cov = []() {
0592     auto c = BoundSquareMatrix::Random().eval();
0593     c *= c.transpose();
0594     return c;
0595   }();
0596
0597   std::vector<std::tuple<double, BoundVector, std::optional<BoundSquareMatrix>>>
0598       cmps(4, {0.25, pars, cov});
0599   BoundTrackParameters check_pars(surface, pars, cov, particleHypothesis);
0600
0601   MultiComponentBoundTrackParameters multi_pars(surface, cmps,
0602                                                 particleHypothesis);
0603   MultiState multi_state(geoCtx, magCtx, defaultBField, multi_pars,
0604                          defaultStepSize);
0605   MultiStepper multi_stepper(defaultBField);
0606
0607   const auto [curv_pars, jac, pathLength] =
0608       multi_stepper.curvilinearState(multi_state);
0609
0610   BOOST_CHECK(
0611       curv_pars.fourPosition(multi_state.geoContext)
0612           .isApprox(check_pars.fourPosition(multi_state.geoContext), 1.e-8));
0613   BOOST_CHECK(curv_pars.direction().isApprox(check_pars.direction(), 1.e-8));
0614   BOOST_CHECK_CLOSE(curv_pars.absoluteMomentum(), check_pars.absoluteMomentum(),
0615                     1.e-8);
0616   BOOST_CHECK_CLOSE(curv_pars.charge(), check_pars.charge(), 1.e-8);
0617 }
0618
0619 BOOST_AUTO_TEST_CASE(test_curvilinear_state) {
0620   test_combined_curvilinear_state_function<MultiStepperLoop>();
0621 }
0622
0623 ////////////////////////////////////
0624 // Test single component interface
0625 ////////////////////////////////////
0626
0627 template <typename multi_stepper_t>
0628 void test_single_component_interface_function() {
0629   using MultiState = typename multi_stepper_t::State;
0630   using MultiStepper = multi_stepper_t;
0631
0632   std::vector<std::tuple<double, BoundVector, std::optional<BoundSquareMatrix>>>
0633       cmps;
0634   for (int i = 0; i < 4; ++i) {
0635     cmps.push_back({0.25, BoundVector::Random(), BoundSquareMatrix::Random()});
0636   }
0637
0638   auto surface = Acts::Surface::makeShared<Acts::PlaneSurface>(
0639       Vector3::Zero(), Vector3::Ones().normalized());
0640
0641   MultiComponentBoundTrackParameters multi_pars(surface, cmps,
0642                                                 particleHypothesis);
0643
0644   MultiState multi_state(geoCtx, magCtx, defaultBField, multi_pars,
0645                          defaultStepSize);
0646
0647   MultiStepper multi_stepper(defaultBField);
0648
0649   DummyPropState multi_prop_state(defaultNDir, multi_state);
0650
0651   // Check at least some properties at the moment
0652   auto check = [&](auto cmp) {
0653     auto sstepper = cmp.singleStepper(multi_stepper);
0654     auto &sstepping = cmp.singleState(multi_prop_state).stepping;
0655
0656     BOOST_CHECK_EQUAL(sstepper.position(sstepping),
0657                       cmp.pars().template segment<3>(eFreePos0));
0658     BOOST_CHECK_EQUAL(sstepper.direction(sstepping),
0659                       cmp.pars().template segment<3>(eFreeDir0));
0660     BOOST_CHECK_EQUAL(sstepper.time(sstepping), cmp.pars()[eFreeTime]);
0661     BOOST_CHECK_CLOSE(sstepper.qOverP(sstepping), cmp.pars()[eFreeQOverP],
0662                       1.e-8);
0663   };
0664
0665   for (const auto cmp : multi_stepper.constComponentIterable(multi_state)) {
0666     check(cmp);
0667   }
0668
0669   for (auto cmp : multi_stepper.componentIterable(multi_state)) {
0670     check(cmp);
0671   }
0672 }
0673
0674 BOOST_AUTO_TEST_CASE(test_single_component_interface) {
0675   test_single_component_interface_function<MultiStepperLoop>();
0676 }
0677
0678 //////////////////////////////
0679 // Remove and add components
0680 //////////////////////////////
0681
0682 template <typename multi_stepper_t>
0683 void remove_add_components_function() {
0684   using MultiState = typename multi_stepper_t::State;
0685   using MultiStepper = multi_stepper_t;
0686
0687   const auto multi_pars = makeDefaultBoundPars(4);
0688
0689   MultiState multi_state(geoCtx, magCtx, defaultBField, multi_pars,
0690                          defaultStepSize);
0691
0692   MultiStepper multi_stepper(defaultBField);
0693
0694   {
0695     BoundTrackParameters pars(multi_pars.referenceSurface().getSharedPtr(),
0696                               BoundVector::Ones(), std::nullopt,
0697                               particleHypothesis);
0698     multi_stepper.addComponent(multi_state, pars, 0.0);
0699   }
0700
0701   BOOST_CHECK_EQUAL(multi_stepper.numberComponents(multi_state),
0702                     multi_pars.components().size() + 1);
0703
0704   multi_stepper.clearComponents(multi_state);
0705
0706   BOOST_CHECK_EQUAL(multi_stepper.numberComponents(multi_state), 0);
0707 }
0708
0709 BOOST_AUTO_TEST_CASE(remove_add_components_test) {
0710   remove_add_components_function<MultiStepperLoop>();
0711 }
0712
0713 //////////////////////////////////////////////////
0714 // Instantiate a Propagator with the MultiStepper
0715 //////////////////////////////////////////////////
0716
0717 template <typename multi_stepper_t>
0718 void propagator_instatiation_test_function() {
0719   auto bField = std::make_shared<NullBField>();
0720   multi_stepper_t multi_stepper(bField);
0721   Propagator<multi_stepper_t, Navigator> propagator(
0722       std::move(multi_stepper), Navigator{Navigator::Config{}});
0723
0724   auto surface = Acts::Surface::makeShared<Acts::PlaneSurface>(
0725       Vector3::Zero(), Vector3{1.0, 0.0, 0.0});
0726   PropagatorOptions options(geoCtx, magCtx);
0727
0728   std::vector<std::tuple<double, BoundVector, std::optional<BoundSquareMatrix>>>
0729       cmps(4, {0.25, BoundVector::Ones().eval(),
0730                BoundSquareMatrix::Identity().eval()});
0731   MultiComponentBoundTrackParameters pars(surface, cmps, particleHypothesis);
0732
0733   // This only checks that this compiles, not that it runs without errors
0734   // @TODO: Add test that checks the target aborter works correctly
0735
0736   // Instantiate with target
0737   using type_a =
0738       decltype(propagator.template propagate<decltype(pars), decltype(options),
0739                                              MultiStepperSurfaceReached>(
0740           pars, *surface, options));
0741
0742   // Instantiate without target
0743   using tybe_b = decltype(propagator.propagate(pars, options));
0744 }
0745
0746 BOOST_AUTO_TEST_CASE(propagator_instatiation_test) {
0747   propagator_instatiation_test_function<MultiStepperLoop>();
0748 }