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