Core/Propagator/MultiStepperTests.cpp

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