Core/Propagator/MultiStepperTests.hpp

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