Core/Propagator/EigenStepperTests.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/Definitions/Units.hpp"
0016 #include "Acts/EventData/GenericBoundTrackParameters.hpp"
0017 #include "Acts/EventData/ParticleHypothesis.hpp"
0018 #include "Acts/EventData/TrackParameters.hpp"
0019 #include "Acts/EventData/TransformationHelpers.hpp"
0020 #include "Acts/Geometry/BoundarySurfaceT.hpp"
0021 #include "Acts/Geometry/CuboidVolumeBuilder.hpp"
0022 #include "Acts/Geometry/GeometryContext.hpp"
0023 #include "Acts/Geometry/TrackingGeometry.hpp"
0024 #include "Acts/Geometry/TrackingGeometryBuilder.hpp"
0025 #include "Acts/Geometry/TrackingVolume.hpp"
0026 #include "Acts/MagneticField/ConstantBField.hpp"
0027 #include "Acts/MagneticField/MagneticFieldContext.hpp"
0028 #include "Acts/MagneticField/MagneticFieldProvider.hpp"
0029 #include "Acts/MagneticField/NullBField.hpp"
0030 #include "Acts/Material/HomogeneousSurfaceMaterial.hpp"
0031 #include "Acts/Material/HomogeneousVolumeMaterial.hpp"
0032 #include "Acts/Material/MaterialSlab.hpp"
0033 #include "Acts/Propagator/ActorList.hpp"
0034 #include "Acts/Propagator/ConstrainedStep.hpp"
0035 #include "Acts/Propagator/EigenStepper.hpp"
0036 #include "Acts/Propagator/EigenStepperDenseExtension.hpp"
0037 #include "Acts/Propagator/EigenStepperError.hpp"
0038 #include "Acts/Propagator/MaterialInteractor.hpp"
0039 #include "Acts/Propagator/Navigator.hpp"
0040 #include "Acts/Propagator/Propagator.hpp"
0041 #include "Acts/Surfaces/BoundaryTolerance.hpp"
0042 #include "Acts/Surfaces/CurvilinearSurface.hpp"
0043 #include "Acts/Surfaces/PlaneSurface.hpp"
0044 #include "Acts/Surfaces/RectangleBounds.hpp"
0045 #include "Acts/Surfaces/Surface.hpp"
0046 #include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
0047 #include "Acts/Tests/CommonHelpers/PredefinedMaterials.hpp"
0048 #include "Acts/Utilities/Logger.hpp"
0049 #include "Acts/Utilities/Result.hpp"
0050 #include "Acts/Utilities/UnitVectors.hpp"
0051
0052 #include <cmath>
0053 #include <limits>
0054 #include <memory>
0055 #include <numbers>
0056 #include <optional>
0057 #include <string>
0058 #include <type_traits>
0059 #include <utility>
0060 #include <vector>
0061
0062 using namespace Acts::UnitLiterals;
0063 using Acts::VectorHelpers::makeVector4;
0064
0065 namespace Acts::Test {
0066
0067 using Covariance = BoundSquareMatrix;
0068
0069 static constexpr auto eps = 3 * std::numeric_limits<double>::epsilon();
0070
0071 // Create a test context
0072 GeometryContext tgContext = GeometryContext();
0073 MagneticFieldContext mfContext = MagneticFieldContext();
0074
0075 /// @brief Aborter for the case that a particle leaves the detector or reaches
0076 /// a custom made threshold.
0077 ///
0078 struct EndOfWorld {
0079   /// Maximum value in x-direction of the detector
0080   double maxX = 1_m;
0081
0082   /// @brief Main call operator for the abort operation
0083   ///
0084   /// @tparam propagator_state_t State of the propagator
0085   /// @tparam stepper_t Type of the stepper
0086   /// @tparam navigator_t Type of the navigator
0087   ///
0088   /// @param [in] state State of the propagation
0089   /// @param [in] stepper Stepper of the propagation
0090   /// @param [in] navigator Navigator of the propagation
0091   ///
0092   /// @return Boolean statement if the particle is still in the detector
0093   template <typename propagator_state_t, typename stepper_t,
0094             typename navigator_t>
0095   bool checkAbort(propagator_state_t& state, const stepper_t& stepper,
0096                   const navigator_t& /*navigator*/,
0097                   const Logger& /*logger*/) const {
0098     const double tolerance = state.options.surfaceTolerance;
0099     if (maxX - std::abs(stepper.position(state.stepping).x()) <= tolerance ||
0100         std::abs(stepper.position(state.stepping).y()) >= 0.5_m ||
0101         std::abs(stepper.position(state.stepping).z()) >= 0.5_m) {
0102       return true;
0103     }
0104     return false;
0105   }
0106 };
0107
0108 ///
0109 /// @brief Data collector while propagation
0110 ///
0111 struct StepCollector {
0112   ///
0113   /// @brief Data container for result analysis
0114   ///
0115   struct this_result {
0116     // Position of the propagator after each step
0117     std::vector<Vector3> position;
0118     // Momentum of the propagator after each step
0119     std::vector<Vector3> momentum;
0120   };
0121
0122   using result_type = this_result;
0123
0124   /// @brief Main call operator for the action list. It stores the data for
0125   /// analysis afterwards
0126   ///
0127   /// @tparam propagator_state_t Type of the propagator state
0128   /// @tparam stepper_t Type of the stepper
0129   /// @tparam navigator_t Type of the navigator
0130   ///
0131   /// @param [in] state State of the propagator
0132   /// @param [in] stepper Stepper of the propagation
0133   /// @param [in] navigator Navigator of the propagation
0134   /// @param [out] result Struct which is filled with the data
0135   template <typename propagator_state_t, typename stepper_t,
0136             typename navigator_t>
0137   void act(propagator_state_t& state, const stepper_t& stepper,
0138            const navigator_t& /*navigator*/, result_type& result,
0139            const Logger& /*logger*/) const {
0140     result.position.push_back(stepper.position(state.stepping));
0141     result.momentum.push_back(stepper.momentum(state.stepping));
0142   }
0143
0144   template <typename propagator_state_t, typename stepper_t,
0145             typename navigator_t>
0146   bool checkAbort(propagator_state_t& /*state*/, const stepper_t& /*stepper*/,
0147                   const navigator_t& /*navigator*/, result_type& /*result*/,
0148                   const Logger& /*logger*/) const {
0149     return false;
0150   }
0151 };
0152
0153 /// These tests are aiming to test whether the state setup is working properly
0154 BOOST_AUTO_TEST_CASE(eigen_stepper_state_test) {
0155   // Set up some variables
0156   auto bField = std::make_shared<ConstantBField>(Vector3(1., 2.5, 33.33));
0157
0158   Vector3 pos(1., 2., 3.);
0159   Vector3 dir(4., 5., 6.);
0160   double time = 7.;
0161   double absMom = 8.;
0162   double charge = -1.;
0163
0164   EigenStepper<>::Options esOptions(tgContext, mfContext);
0165
0166   // Test charged parameters without covariance matrix
0167   BoundTrackParameters cp = BoundTrackParameters::createCurvilinear(
0168       makeVector4(pos, time), dir, charge / absMom, std::nullopt,
0169       ParticleHypothesis::pion());
0170   EigenStepper<> es(bField);
0171   EigenStepper<>::State esState = es.makeState(esOptions);
0172   es.initialize(esState, cp);
0173
0174   // Test the result & compare with the input/test for reasonable members
0175   BOOST_CHECK_EQUAL(esState.jacToGlobal, BoundToFreeMatrix::Zero());
0176   BOOST_CHECK_EQUAL(esState.jacTransport, FreeMatrix::Identity());
0177   BOOST_CHECK_EQUAL(esState.derivative, FreeVector::Zero());
0178   BOOST_CHECK(!esState.covTransport);
0179   BOOST_CHECK_EQUAL(esState.cov, Covariance::Zero());
0180   BOOST_CHECK_EQUAL(esState.pathAccumulated, 0.);
0181   BOOST_CHECK_EQUAL(esState.previousStepSize, 0.);
0182
0183   // Test without charge and covariance matrix
0184   BoundTrackParameters ncp = BoundTrackParameters::createCurvilinear(
0185       makeVector4(pos, time), dir, 1 / absMom, std::nullopt,
0186       ParticleHypothesis::pion0());
0187   esOptions = EigenStepper<>::Options(tgContext, mfContext);
0188   es.initialize(esState, ncp);
0189   BOOST_CHECK_EQUAL(es.charge(esState), 0.);
0190
0191   // Test with covariance matrix
0192   Covariance cov = 8. * Covariance::Identity();
0193   ncp = BoundTrackParameters::createCurvilinear(makeVector4(pos, time), dir,
0194                                                 1 / absMom, cov,
0195                                                 ParticleHypothesis::pion0());
0196   es.initialize(esState, ncp);
0197   BOOST_CHECK_NE(esState.jacToGlobal, BoundToFreeMatrix::Zero());
0198   BOOST_CHECK(esState.covTransport);
0199   BOOST_CHECK_EQUAL(esState.cov, cov);
0200 }
0201
0202 /// These tests are aiming to test the functions of the EigenStepper
0203 /// The numerical correctness of the stepper is tested in the integration tests
0204 BOOST_AUTO_TEST_CASE(eigen_stepper_test) {
0205   // Set up some variables for the state
0206   Direction navDir = Direction::Backward();
0207   double stepSize = 123.;
0208   auto bField = std::make_shared<ConstantBField>(Vector3(1., 2.5, 33.33));
0209   auto bCache = bField->makeCache(mfContext);
0210
0211   // Construct the parameters
0212   Vector3 pos(1., 2., 3.);
0213   Vector3 dir = Vector3(4., 5., 6.).normalized();
0214   double time = 7.;
0215   double absMom = 8.;
0216   double charge = -1.;
0217   Covariance cov = 8. * Covariance::Identity();
0218   BoundTrackParameters cp = BoundTrackParameters::createCurvilinear(
0219       makeVector4(pos, time), dir, charge / absMom, cov,
0220       ParticleHypothesis::pion());
0221
0222   EigenStepper<>::Options esOptions(tgContext, mfContext);
0223   esOptions.maxStepSize = stepSize;
0224   esOptions.initialStepSize = 10_m;
0225
0226   // Build the stepper and the state
0227   EigenStepper<> es(bField);
0228   EigenStepper<>::State esState = es.makeState(esOptions);
0229   es.initialize(esState, cp);
0230
0231   // Test the getters
0232   CHECK_CLOSE_ABS(es.position(esState), pos, eps);
0233   CHECK_CLOSE_ABS(es.direction(esState), dir, eps);
0234   CHECK_CLOSE_ABS(es.absoluteMomentum(esState), absMom, eps);
0235   CHECK_CLOSE_ABS(es.charge(esState), charge, eps);
0236   CHECK_CLOSE_ABS(es.time(esState), time, eps);
0237   BOOST_CHECK_EQUAL(es.getField(esState, pos).value(),
0238                     bField->getField(pos, bCache).value());
0239
0240   // Step size modifies
0241   const std::string originalStepSize = esState.stepSize.toString();
0242
0243   es.updateStepSize(esState, -1337., ConstrainedStep::Type::Navigator);
0244   BOOST_CHECK_EQUAL(esState.previousStepSize, stepSize);
0245   BOOST_CHECK_EQUAL(esState.stepSize.value(), -1337.);
0246
0247   es.releaseStepSize(esState, ConstrainedStep::Type::Navigator);
0248   BOOST_CHECK_EQUAL(esState.stepSize.value(), stepSize);
0249   BOOST_CHECK_EQUAL(es.outputStepSize(esState), originalStepSize);
0250
0251   // Test the curvilinear state construction
0252   auto curvState = es.curvilinearState(esState);
0253   auto curvPars = std::get<0>(curvState);
0254   CHECK_CLOSE_ABS(curvPars.position(tgContext), cp.position(tgContext), eps);
0255   CHECK_CLOSE_ABS(curvPars.momentum(), cp.momentum(), 10e-6);
0256   CHECK_CLOSE_ABS(curvPars.charge(), cp.charge(), eps);
0257   CHECK_CLOSE_ABS(curvPars.time(), cp.time(), eps);
0258   BOOST_CHECK(curvPars.covariance().has_value());
0259   BOOST_CHECK_NE(*curvPars.covariance(), cov);
0260   CHECK_CLOSE_COVARIANCE(std::get<1>(curvState),
0261                          BoundMatrix(BoundMatrix::Identity()), eps);
0262   CHECK_CLOSE_ABS(std::get<2>(curvState), 0., eps);
0263
0264   // Test the update method
0265   Vector3 newPos(2., 4., 8.);
0266   Vector3 newMom(3., 9., 27.);
0267   double newTime(321.);
0268   es.update(esState, newPos, newMom.normalized(), charge / newMom.norm(),
0269             newTime);
0270   BOOST_CHECK_EQUAL(es.position(esState), newPos);
0271   BOOST_CHECK_EQUAL(es.direction(esState), newMom.normalized());
0272   BOOST_CHECK_EQUAL(es.absoluteMomentum(esState), newMom.norm());
0273   BOOST_CHECK_EQUAL(es.charge(esState), charge);
0274   BOOST_CHECK_EQUAL(es.time(esState), newTime);
0275
0276   // The covariance transport
0277   esState.cov = cov;
0278   es.transportCovarianceToCurvilinear(esState);
0279   BOOST_CHECK_NE(esState.cov, cov);
0280   BOOST_CHECK_NE(esState.jacToGlobal, BoundToFreeMatrix::Zero());
0281   BOOST_CHECK_EQUAL(esState.jacTransport, FreeMatrix::Identity());
0282   BOOST_CHECK_EQUAL(esState.derivative, FreeVector::Zero());
0283
0284   // Perform a step without and with covariance transport
0285   esState.cov = cov;
0286
0287   esState.covTransport = false;
0288   es.step(esState, navDir, nullptr).value();
0289   CHECK_CLOSE_COVARIANCE(esState.cov, cov, eps);
0290   BOOST_CHECK_NE(es.position(esState).norm(), newPos.norm());
0291   BOOST_CHECK_NE(es.direction(esState), newMom.normalized());
0292   BOOST_CHECK_EQUAL(es.charge(esState), charge);
0293   BOOST_CHECK_LT(es.time(esState), newTime);
0294   BOOST_CHECK_EQUAL(esState.derivative, FreeVector::Zero());
0295   BOOST_CHECK_EQUAL(esState.jacTransport, FreeMatrix::Identity());
0296
0297   esState.covTransport = true;
0298   es.step(esState, navDir, nullptr).value();
0299   CHECK_CLOSE_COVARIANCE(esState.cov, cov, eps);
0300   BOOST_CHECK_NE(es.position(esState).norm(), newPos.norm());
0301   BOOST_CHECK_NE(es.direction(esState), newMom.normalized());
0302   BOOST_CHECK_EQUAL(es.charge(esState), charge);
0303   BOOST_CHECK_LT(es.time(esState), newTime);
0304   BOOST_CHECK_NE(esState.derivative, FreeVector::Zero());
0305   BOOST_CHECK_NE(esState.jacTransport, FreeMatrix::Identity());
0306
0307   /// Test the state reset
0308   // Construct the parameters
0309   Vector3 pos2(1.5, -2.5, 3.5);
0310   Vector3 dir2 = Vector3(4.5, -5.5, 6.5).normalized();
0311   double time2 = 7.5;
0312   double absMom2 = 8.5;
0313   double charge2 = 1.;
0314   BoundSquareMatrix cov2 = 8.5 * Covariance::Identity();
0315   BoundTrackParameters cp2 = BoundTrackParameters::createCurvilinear(
0316       makeVector4(pos2, time2), dir2, charge2 / absMom2, cov2,
0317       ParticleHypothesis::pion());
0318   FreeVector freeParams = transformBoundToFreeParameters(
0319       cp2.referenceSurface(), tgContext, cp2.parameters());
0320   navDir = Direction::Forward();
0321
0322   auto copyState = [&](auto& field, const auto& state) {
0323     using field_t = std::decay_t<decltype(field)>;
0324     std::decay_t<decltype(state)> copy = es.makeState(esOptions);
0325     es.initialize(copy, cp);
0326     copy.pars = state.pars;
0327     copy.covTransport = state.covTransport;
0328     copy.cov = state.cov;
0329     copy.jacobian = state.jacobian;
0330     copy.jacToGlobal = state.jacToGlobal;
0331     copy.jacTransport = state.jacTransport;
0332     copy.derivative = state.derivative;
0333     copy.pathAccumulated = state.pathAccumulated;
0334     copy.stepSize = state.stepSize;
0335     copy.previousStepSize = state.previousStepSize;
0336
0337     copy.fieldCache = MagneticFieldProvider::Cache(
0338         std::in_place_type<typename field_t::Cache>,
0339         state.fieldCache.template as<typename field_t::Cache>());
0340
0341     copy.extension = state.extension;
0342     copy.stepData = state.stepData;
0343
0344     return copy;
0345   };
0346
0347   // Reset all possible parameters
0348   EigenStepper<>::State esStateCopy = copyState(*bField, esState);
0349   BOOST_CHECK(cp2.covariance().has_value());
0350   es.initialize(esStateCopy, cp2.parameters(), *cp2.covariance(),
0351                 cp2.particleHypothesis(), cp2.referenceSurface());
0352   // Test all components
0353   BOOST_CHECK_NE(esStateCopy.jacToGlobal, BoundToFreeMatrix::Zero());
0354   BOOST_CHECK_NE(esStateCopy.jacToGlobal, esState.jacToGlobal);
0355   BOOST_CHECK_EQUAL(esStateCopy.jacTransport, FreeMatrix::Identity());
0356   BOOST_CHECK_EQUAL(esStateCopy.derivative, FreeVector::Zero());
0357   BOOST_CHECK(esStateCopy.covTransport);
0358   BOOST_CHECK_EQUAL(esStateCopy.cov, cov2);
0359   BOOST_CHECK_EQUAL(es.position(esStateCopy),
0360                     freeParams.template segment<3>(eFreePos0));
0361   BOOST_CHECK_EQUAL(es.direction(esStateCopy),
0362                     freeParams.template segment<3>(eFreeDir0).normalized());
0363   BOOST_CHECK_EQUAL(es.absoluteMomentum(esStateCopy),
0364                     std::abs(1. / freeParams[eFreeQOverP]));
0365   BOOST_CHECK_EQUAL(es.charge(esStateCopy), -es.charge(esState));
0366   BOOST_CHECK_EQUAL(es.time(esStateCopy), freeParams[eFreeTime]);
0367   BOOST_CHECK_EQUAL(esStateCopy.pathAccumulated, 0.);
0368   BOOST_CHECK_EQUAL(esStateCopy.stepSize.value(), stepSize);
0369   BOOST_CHECK_EQUAL(esStateCopy.previousStepSize, 0.);
0370
0371   /// Repeat with surface related methods
0372   std::shared_ptr<PlaneSurface> plane =
0373       CurvilinearSurface(pos, dir.normalized()).planeSurface();
0374   auto bp = BoundTrackParameters::create(
0375                 tgContext, plane, makeVector4(pos, time), dir, charge / absMom,
0376                 cov, ParticleHypothesis::pion())
0377                 .value();
0378   es.initialize(esState, bp);
0379
0380   // Test the intersection in the context of a surface
0381   std::shared_ptr<PlaneSurface> targetSurface =
0382       CurvilinearSurface(pos + navDir * 2. * dir, dir).planeSurface();
0383   es.updateSurfaceStatus(esState, *targetSurface, 0, navDir,
0384                          BoundaryTolerance::Infinite(), s_onSurfaceTolerance,
0385                          ConstrainedStep::Type::Navigator);
0386   CHECK_CLOSE_ABS(esState.stepSize.value(ConstrainedStep::Type::Navigator),
0387                   navDir * 2., eps);
0388
0389   // Test the step size modification in the context of a surface
0390   es.updateStepSize(esState,
0391                     targetSurface
0392                         ->intersect(tgContext, es.position(esState),
0393                                     navDir * es.direction(esState),
0394                                     BoundaryTolerance::Infinite())
0395                         .closest(),
0396                     navDir, ConstrainedStep::Type::Navigator);
0397   CHECK_CLOSE_ABS(esState.stepSize.value(), 2., eps);
0398   esState.stepSize.setUser(navDir * stepSize);
0399   es.releaseStepSize(esState, ConstrainedStep::Type::Navigator);
0400   es.updateStepSize(esState,
0401                     targetSurface
0402                         ->intersect(tgContext, es.position(esState),
0403                                     navDir * es.direction(esState),
0404                                     BoundaryTolerance::Infinite())
0405                         .closest(),
0406                     navDir, ConstrainedStep::Type::Navigator);
0407   CHECK_CLOSE_ABS(esState.stepSize.value(), 2., eps);
0408
0409   // Test the bound state construction
0410   auto boundState = es.boundState(esState, *plane).value();
0411   auto boundPars = std::get<0>(boundState);
0412   CHECK_CLOSE_ABS(boundPars.position(tgContext), bp.position(tgContext), eps);
0413   CHECK_CLOSE_ABS(boundPars.momentum(), bp.momentum(), 1e-7);
0414   CHECK_CLOSE_ABS(boundPars.charge(), bp.charge(), eps);
0415   CHECK_CLOSE_ABS(boundPars.time(), bp.time(), eps);
0416   BOOST_CHECK(boundPars.covariance().has_value());
0417   BOOST_CHECK_NE(*boundPars.covariance(), cov);
0418   CHECK_CLOSE_COVARIANCE(std::get<1>(boundState),
0419                          BoundMatrix(BoundMatrix::Identity()), eps);
0420   CHECK_CLOSE_ABS(std::get<2>(boundState), 0., eps);
0421
0422   // Transport the covariance in the context of a surface
0423   es.transportCovarianceToBound(esState, *plane);
0424   BOOST_CHECK_NE(esState.cov, cov);
0425   BOOST_CHECK_NE(esState.jacToGlobal, BoundToFreeMatrix::Zero());
0426   BOOST_CHECK_EQUAL(esState.jacTransport, FreeMatrix::Identity());
0427   BOOST_CHECK_EQUAL(esState.derivative, FreeVector::Zero());
0428
0429   // Update in context of a surface
0430   freeParams = transformBoundToFreeParameters(bp.referenceSurface(), tgContext,
0431                                               bp.parameters());
0432
0433   es.update(esState, freeParams, bp.parameters(), 2 * (*bp.covariance()),
0434             *plane);
0435   CHECK_CLOSE_OR_SMALL(es.position(esState), pos, eps, eps);
0436   CHECK_CLOSE_OR_SMALL(es.direction(esState), dir, eps, eps);
0437   CHECK_CLOSE_REL(es.absoluteMomentum(esState), absMom, eps);
0438   BOOST_CHECK_EQUAL(es.charge(esState), charge);
0439   CHECK_CLOSE_OR_SMALL(es.time(esState), time, eps, eps);
0440   CHECK_CLOSE_COVARIANCE(esState.cov, Covariance(2. * cov), eps);
0441
0442   // Test a case where no step size adjustment is required
0443   esState.options.stepTolerance = 2. * 4.4258e+09;
0444   double h0 = esState.stepSize.value();
0445   es.step(esState, navDir, nullptr);
0446   CHECK_CLOSE_ABS(h0, esState.stepSize.value(), eps);
0447
0448   // Produce some errors
0449   auto nBfield = std::make_shared<NullBField>();
0450   EigenStepper<> nes(nBfield);
0451   EigenStepper<>::Options nesOptions(tgContext, mfContext);
0452   EigenStepper<>::State nesState = nes.makeState(nesOptions);
0453   nes.initialize(nesState, cp);
0454   auto nEsState = copyState(*nBfield, nesState);
0455   // Test that we can reach the minimum step size
0456   nEsState.options.stepTolerance = 1e-21;
0457   nEsState.options.stepSizeCutOff = 1e20;
0458   auto res = nes.step(nEsState, navDir, nullptr);
0459   BOOST_CHECK(!res.ok());
0460   BOOST_CHECK_EQUAL(res.error(), EigenStepperError::StepSizeStalled);
0461
0462   // Test that the number of trials exceeds
0463   nEsState.options.stepSizeCutOff = 0.;
0464   nEsState.options.maxRungeKuttaStepTrials = 0.;
0465   res = nes.step(nEsState, navDir, nullptr);
0466   BOOST_CHECK(!res.ok());
0467   BOOST_CHECK_EQUAL(res.error(), EigenStepperError::StepSizeAdjustmentFailed);
0468 }
0469
0470 /// @brief This function tests the EigenStepper with the EigenStepperDefaultExtension and
0471 /// the EigenStepperDenseExtension. The focus of this tests lies in
0472 /// the choosing of the right extension for the individual use case. This is
0473 /// performed with three different detectors:
0474 /// a) Pure vacuum -> DefaultExtension needs to act
0475 /// b) Pure Be -> DenseEnvironmentExtension needs to act
0476 /// c) Vacuum - Be - Vacuum -> Both should act and switch during the propagation
0477
0478 // Test case a). The DenseEnvironmentExtension should state that it is not
0479 // valid in this case.
0480 BOOST_AUTO_TEST_CASE(step_extension_vacuum_test) {
0481   CuboidVolumeBuilder cvb;
0482   CuboidVolumeBuilder::VolumeConfig vConf;
0483   vConf.position = {0.5_m, 0., 0.};
0484   vConf.length = {1_m, 1_m, 1_m};
0485   CuboidVolumeBuilder::Config conf;
0486   conf.volumeCfg.push_back(vConf);
0487   conf.position = {0.5_m, 0., 0.};
0488   conf.length = {1_m, 1_m, 1_m};
0489
0490   // Build detector
0491   cvb.setConfig(conf);
0492   TrackingGeometryBuilder::Config tgbCfg;
0493   tgbCfg.trackingVolumeBuilders.push_back(
0494       [=](const auto& context, const auto& inner, const auto& vb) {
0495         return cvb.trackingVolume(context, inner, vb);
0496       });
0497   TrackingGeometryBuilder tgb(tgbCfg);
0498   std::shared_ptr<const TrackingGeometry> vacuum =
0499       tgb.trackingGeometry(tgContext);
0500
0501   // Build navigator
0502   Navigator naviVac({vacuum, true, true, true});
0503
0504   // Set initial parameters for the particle track
0505   Covariance cov = Covariance::Identity();
0506   const Vector3 startDir = makeDirectionFromPhiTheta(0_degree, 90_degree);
0507   const Vector3 startMom = 1_GeV * startDir;
0508   const BoundTrackParameters sbtp = BoundTrackParameters::createCurvilinear(
0509       Vector4::Zero(), startDir, 1_e / 1_GeV, cov, ParticleHypothesis::pion());
0510
0511   using Stepper = EigenStepper<EigenStepperDenseExtension>;
0512   using Propagator = Propagator<Stepper, Navigator>;
0513   using PropagatorOptions =
0514       Propagator::Options<ActorList<StepCollector, EndOfWorld>>;
0515
0516   // Set options for propagator
0517   PropagatorOptions propOpts(tgContext, mfContext);
0518   propOpts.maxSteps = 100;
0519   propOpts.stepping.maxStepSize = 1.5_m;
0520
0521   // Build stepper and propagator
0522   auto bField = std::make_shared<ConstantBField>(Vector3(0., 0., 0.));
0523   Stepper es(bField);
0524   Propagator prop(es, naviVac);
0525
0526   // Launch and collect results
0527   const auto& result = prop.propagate(sbtp, propOpts).value();
0528   const StepCollector::this_result& stepResult =
0529       result.get<typename StepCollector::result_type>();
0530
0531   // Check that the propagation happened without interactions
0532   for (const auto& pos : stepResult.position) {
0533     CHECK_SMALL(pos.y(), 1_um);
0534     CHECK_SMALL(pos.z(), 1_um);
0535     if (pos == stepResult.position.back()) {
0536       CHECK_CLOSE_ABS(pos.x(), 1_m, 1_um);
0537     }
0538   }
0539   for (const auto& mom : stepResult.momentum) {
0540     CHECK_CLOSE_ABS(mom, startMom, 1_keV);
0541   }
0542
0543   using DefStepper = EigenStepper<EigenStepperDenseExtension>;
0544   using DefPropagator = Acts::Propagator<DefStepper, Navigator>;
0545   using DefPropagatorOptions =
0546       DefPropagator::Options<ActorList<StepCollector, EndOfWorld>>;
0547
0548   // Set options for propagator
0549   DefPropagatorOptions propOptsDef(tgContext, mfContext);
0550   propOptsDef.maxSteps = 100;
0551   propOptsDef.stepping.maxStepSize = 1.5_m;
0552
0553   DefStepper esDef(bField);
0554   DefPropagator propDef(esDef, naviVac);
0555
0556   // Launch and collect results
0557   const auto& resultDef = propDef.propagate(sbtp, propOptsDef).value();
0558   const StepCollector::this_result& stepResultDef =
0559       resultDef.get<typename StepCollector::result_type>();
0560
0561   // Check that the right extension was chosen
0562   // If chosen correctly, the number of elements should be identical
0563   BOOST_CHECK_EQUAL(stepResult.position.size(), stepResultDef.position.size());
0564   for (unsigned int i = 0; i < stepResult.position.size(); i++) {
0565     CHECK_CLOSE_ABS(stepResult.position[i], stepResultDef.position[i], 1_um);
0566   }
0567   BOOST_CHECK_EQUAL(stepResult.momentum.size(), stepResultDef.momentum.size());
0568   for (unsigned int i = 0; i < stepResult.momentum.size(); i++) {
0569     CHECK_CLOSE_ABS(stepResult.momentum[i], stepResultDef.momentum[i], 1_keV);
0570   }
0571 }
0572 // Test case b). The DefaultExtension should state that it is invalid here.
0573 BOOST_AUTO_TEST_CASE(step_extension_material_test) {
0574   CuboidVolumeBuilder cvb;
0575   CuboidVolumeBuilder::VolumeConfig vConf;
0576   vConf.position = {0.5_m, 0., 0.};
0577   vConf.length = {1_m, 1_m, 1_m};
0578   vConf.volumeMaterial =
0579       std::make_shared<HomogeneousVolumeMaterial>(makeBeryllium());
0580   CuboidVolumeBuilder::Config conf;
0581   conf.volumeCfg.push_back(vConf);
0582   conf.position = {0.5_m, 0., 0.};
0583   conf.length = {1_m, 1_m, 1_m};
0584
0585   // Build detector
0586   cvb.setConfig(conf);
0587   TrackingGeometryBuilder::Config tgbCfg;
0588   tgbCfg.trackingVolumeBuilders.push_back(
0589       [=](const auto& context, const auto& inner, const auto& vb) {
0590         return cvb.trackingVolume(context, inner, vb);
0591       });
0592   TrackingGeometryBuilder tgb(tgbCfg);
0593   std::shared_ptr<const TrackingGeometry> material =
0594       tgb.trackingGeometry(tgContext);
0595
0596   // Build navigator
0597   Navigator naviMat({material, true, true, true});
0598
0599   // Set initial parameters for the particle track
0600   Covariance cov = Covariance::Identity();
0601   const Vector3 startDir = makeDirectionFromPhiTheta(0_degree, 90_degree);
0602   const Vector3 startMom = 5_GeV * startDir;
0603   const BoundTrackParameters sbtp = BoundTrackParameters::createCurvilinear(
0604       Vector4::Zero(), startDir, 1_e / 5_GeV, cov, ParticleHypothesis::pion());
0605
0606   using Stepper = EigenStepper<EigenStepperDenseExtension>;
0607   using Propagator = Propagator<Stepper, Navigator>;
0608   using PropagatorOptions =
0609       Propagator::Options<ActorList<StepCollector, EndOfWorld>>;
0610
0611   // Set options for propagator
0612   PropagatorOptions propOpts(tgContext, mfContext);
0613   propOpts.maxSteps = 10000;
0614   propOpts.stepping.maxStepSize = 1.5_m;
0615
0616   // Build stepper and propagator
0617   auto bField = std::make_shared<ConstantBField>(Vector3(0., 0., 0.));
0618   Stepper es(bField);
0619   Propagator prop(es, naviMat,
0620                   Acts::getDefaultLogger("Propagator", Acts::Logging::VERBOSE));
0621
0622   // Launch and collect results
0623   const auto& result = prop.propagate(sbtp, propOpts).value();
0624   const StepCollector::this_result& stepResult =
0625       result.get<typename StepCollector::result_type>();
0626
0627   // Check that there occurred interaction
0628   for (const auto& pos : stepResult.position) {
0629     CHECK_SMALL(pos.y(), 1_um);
0630     CHECK_SMALL(pos.z(), 1_um);
0631     if (pos == stepResult.position.front()) {
0632       CHECK_SMALL(pos.x(), 1_um);
0633     } else {
0634       BOOST_CHECK_GT(std::abs(pos.x()), 1_um);
0635     }
0636   }
0637   for (const auto& mom : stepResult.momentum) {
0638     CHECK_SMALL(mom.y(), 1_keV);
0639     CHECK_SMALL(mom.z(), 1_keV);
0640     if (mom == stepResult.momentum.front()) {
0641       CHECK_CLOSE_ABS(mom.x(), 5_GeV, 1_keV);
0642     } else {
0643       BOOST_CHECK_LT(mom.x(), 5_GeV);
0644     }
0645   }
0646
0647   using DenseStepper = EigenStepper<EigenStepperDenseExtension>;
0648   using DensePropagator = Acts::Propagator<DenseStepper, Navigator>;
0649   using DensePropagatorOptions =
0650       DensePropagator::Options<ActorList<StepCollector, EndOfWorld>>;
0651
0652   // Rebuild and check the choice of extension
0653   // Set options for propagator
0654   DensePropagatorOptions propOptsDense(tgContext, mfContext);
0655   propOptsDense.maxSteps = 1000;
0656   propOptsDense.stepping.maxStepSize = 1.5_m;
0657
0658   // Build stepper and propagator
0659   DenseStepper esDense(bField);
0660   DensePropagator propDense(esDense, naviMat);
0661
0662   // Launch and collect results
0663   const auto& resultDense = propDense.propagate(sbtp, propOptsDense).value();
0664   const StepCollector::this_result& stepResultDense =
0665       resultDense.get<typename StepCollector::result_type>();
0666
0667   // Check that the right extension was chosen
0668   // If chosen correctly, the number of elements should be identical
0669   BOOST_CHECK_EQUAL(stepResult.position.size(),
0670                     stepResultDense.position.size());
0671   for (unsigned int i = 0; i < stepResult.position.size(); i++) {
0672     CHECK_CLOSE_ABS(stepResult.position[i], stepResultDense.position[i], 1_um);
0673   }
0674   BOOST_CHECK_EQUAL(stepResult.momentum.size(),
0675                     stepResultDense.momentum.size());
0676   for (unsigned int i = 0; i < stepResult.momentum.size(); i++) {
0677     CHECK_CLOSE_ABS(stepResult.momentum[i], stepResultDense.momentum[i], 1_keV);
0678   }
0679
0680   ////////////////////////////////////////////////////////////////////
0681
0682   // Re-launch the configuration with magnetic field
0683   bField->setField(Vector3{0., 1_T, 0.});
0684   Stepper esB(bField);
0685   Propagator propB(esB, naviMat);
0686
0687   const auto& resultB = propB.propagate(sbtp, propOptsDense).value();
0688   const StepCollector::this_result& stepResultB =
0689       resultB.get<typename StepCollector::result_type>();
0690
0691   // Check that there occurred interaction
0692   for (const auto& pos : stepResultB.position) {
0693     if (pos == stepResultB.position.front()) {
0694       CHECK_SMALL(pos, 1_um);
0695     } else {
0696       BOOST_CHECK_GT(std::abs(pos.x()), 1_um);
0697       CHECK_SMALL(pos.y(), 1_um);
0698       BOOST_CHECK_GT(std::abs(pos.z()), 0.125_um);
0699     }
0700   }
0701   for (const auto& mom : stepResultB.momentum) {
0702     if (mom == stepResultB.momentum.front()) {
0703       CHECK_CLOSE_ABS(mom, startMom, 1_keV);
0704     } else {
0705       BOOST_CHECK_NE(mom.x(), 5_GeV);
0706       CHECK_SMALL(mom.y(), 1_keV);
0707       BOOST_CHECK_NE(mom.z(), 0.);
0708     }
0709   }
0710 }
0711 // Test case c). Both should be involved in their part of the detector
0712 BOOST_AUTO_TEST_CASE(step_extension_vacmatvac_test) {
0713   CuboidVolumeBuilder cvb;
0714   CuboidVolumeBuilder::VolumeConfig vConfVac1;
0715   vConfVac1.position = {0.5_m, 0., 0.};
0716   vConfVac1.length = {1_m, 1_m, 1_m};
0717   vConfVac1.name = "First vacuum volume";
0718   CuboidVolumeBuilder::VolumeConfig vConfMat;
0719   vConfMat.position = {1.5_m, 0., 0.};
0720   vConfMat.length = {1_m, 1_m, 1_m};
0721   vConfMat.volumeMaterial =
0722       std::make_shared<const HomogeneousVolumeMaterial>(makeBeryllium());
0723   vConfMat.name = "Material volume";
0724   CuboidVolumeBuilder::VolumeConfig vConfVac2;
0725   vConfVac2.position = {2.5_m, 0., 0.};
0726   vConfVac2.length = {1_m, 1_m, 1_m};
0727   vConfVac2.name = "Second vacuum volume";
0728   CuboidVolumeBuilder::Config conf;
0729   conf.volumeCfg = {vConfVac1, vConfMat, vConfVac2};
0730   conf.position = {1.5_m, 0., 0.};
0731   conf.length = {3_m, 1_m, 1_m};
0732
0733   // Build detector
0734   cvb.setConfig(conf);
0735   TrackingGeometryBuilder::Config tgbCfg;
0736   tgbCfg.trackingVolumeBuilders.push_back(
0737       [=](const auto& context, const auto& inner, const auto& vb) {
0738         return cvb.trackingVolume(context, inner, vb);
0739       });
0740   TrackingGeometryBuilder tgb(tgbCfg);
0741   std::shared_ptr<const TrackingGeometry> det = tgb.trackingGeometry(tgContext);
0742
0743   // Build navigator
0744   Navigator naviDet({det, true, true, true});
0745
0746   // Set initial parameters for the particle track
0747   BoundTrackParameters sbtp = BoundTrackParameters::createCurvilinear(
0748       Vector4::Zero(), 0_degree, 90_degree, 1_e / 5_GeV, Covariance::Identity(),
0749       ParticleHypothesis::pion());
0750
0751   using Stepper = EigenStepper<EigenStepperDenseExtension>;
0752   using Propagator = Acts::Propagator<Stepper, Navigator>;
0753   using PropagatorOptions =
0754       Propagator::Options<ActorList<StepCollector, EndOfWorld>>;
0755
0756   // Set options for propagator
0757   PropagatorOptions propOpts(tgContext, mfContext);
0758   propOpts.actorList.get<EndOfWorld>().maxX = 3_m;
0759   propOpts.maxSteps = 1000;
0760   propOpts.stepping.maxStepSize = 1.5_m;
0761
0762   // Build stepper and propagator
0763   auto bField = std::make_shared<ConstantBField>(Vector3(0., 1_T, 0.));
0764   Stepper es(bField);
0765   Propagator prop(es, naviDet);
0766
0767   // Launch and collect results
0768   const auto& result = prop.propagate(sbtp, propOpts).value();
0769   const StepCollector::this_result& stepResult =
0770       result.get<typename StepCollector::result_type>();
0771
0772   // Manually set the extensions for each step and propagate through each
0773   // volume by propagation to the boundaries
0774   // Collect boundaries
0775   std::vector<Surface const*> surs;
0776   std::vector<std::shared_ptr<const BoundarySurfaceT<TrackingVolume>>>
0777       boundaries = det->lowestTrackingVolume(tgContext, {0.5_m, 0., 0.})
0778                        ->boundarySurfaces();
0779   for (auto& b : boundaries) {
0780     if (b->surfaceRepresentation().center(tgContext).x() == 1_m) {
0781       surs.push_back(&(b->surfaceRepresentation()));
0782       break;
0783     }
0784   }
0785   boundaries =
0786       det->lowestTrackingVolume(tgContext, {1.5_m, 0., 0.})->boundarySurfaces();
0787   for (auto& b : boundaries) {
0788     if (b->surfaceRepresentation().center(tgContext).x() == 2_m) {
0789       surs.push_back(&(b->surfaceRepresentation()));
0790       break;
0791     }
0792   }
0793   boundaries =
0794       det->lowestTrackingVolume(tgContext, {2.5_m, 0., 0.})->boundarySurfaces();
0795   for (auto& b : boundaries) {
0796     if (b->surfaceRepresentation().center(tgContext).x() == 3_m) {
0797       surs.push_back(&(b->surfaceRepresentation()));
0798       break;
0799     }
0800   }
0801
0802   // Build launcher through vacuum
0803   // Set options for propagator
0804
0805   using DefStepper = EigenStepper<EigenStepperDenseExtension>;
0806   using DefPropagator = Acts::Propagator<DefStepper, Navigator>;
0807   using DefPropagatorOptions =
0808       DefPropagator::Options<ActorList<StepCollector, EndOfWorld>>;
0809
0810   DefPropagatorOptions propOptsDef(tgContext, mfContext);
0811   propOptsDef.actorList.get<EndOfWorld>().maxX = 3_m;
0812   propOptsDef.maxSteps = 1000;
0813   propOptsDef.stepping.maxStepSize = 1.5_m;
0814
0815   // Build stepper and propagator
0816   DefStepper esDef(bField);
0817   DefPropagator propDef(esDef, naviDet);
0818
0819   // Launch and collect results
0820   const auto& resultDef = propDef.propagate(sbtp, propOptsDef).value();
0821   const StepCollector::this_result& stepResultDef =
0822       resultDef.get<typename StepCollector::result_type>();
0823
0824   // Check the exit situation of the first volume
0825   std::pair<Vector3, Vector3> endParams, endParamsControl;
0826   for (unsigned int i = 0; i < stepResultDef.position.size(); i++) {
0827     if (1_m - stepResultDef.position[i].x() < 1e-4) {
0828       endParams =
0829           std::make_pair(stepResultDef.position[i], stepResultDef.momentum[i]);
0830       break;
0831     }
0832   }
0833   for (unsigned int i = 0; i < stepResult.position.size(); i++) {
0834     if (1_m - stepResult.position[i].x() < 1e-4) {
0835       endParamsControl =
0836           std::make_pair(stepResult.position[i], stepResult.momentum[i]);
0837       break;
0838     }
0839   }
0840
0841   CHECK_CLOSE_ABS(endParams.first, endParamsControl.first, 1_um);
0842   CHECK_CLOSE_ABS(endParams.second, endParamsControl.second, 1_um);
0843
0844   CHECK_CLOSE_ABS(endParams.first.x(), endParamsControl.first.x(), 1e-5);
0845   CHECK_CLOSE_ABS(endParams.first.y(), endParamsControl.first.y(), 1e-5);
0846   CHECK_CLOSE_ABS(endParams.first.z(), endParamsControl.first.z(), 1e-5);
0847   CHECK_CLOSE_ABS(endParams.second.x(), endParamsControl.second.x(), 1e-5);
0848   CHECK_CLOSE_ABS(endParams.second.y(), endParamsControl.second.y(), 1e-5);
0849   CHECK_CLOSE_ABS(endParams.second.z(), endParamsControl.second.z(), 1e-5);
0850
0851   // Build launcher through material
0852   // Set initial parameters for the particle track by using the result of the
0853   // first volume
0854
0855   using DenseStepper = EigenStepper<EigenStepperDenseExtension>;
0856   using DensePropagator = Acts::Propagator<DenseStepper, Navigator>;
0857   using DensePropagatorOptions =
0858       DensePropagator::Options<ActorList<StepCollector, EndOfWorld>>;
0859
0860   // Set options for propagator
0861   DensePropagatorOptions propOptsDense(tgContext, mfContext);
0862   propOptsDense.actorList.get<EndOfWorld>().maxX = 3_m;
0863   propOptsDense.maxSteps = 1000;
0864   propOptsDense.stepping.maxStepSize = 1.5_m;
0865
0866   // Build stepper and propagator
0867   DenseStepper esDense(bField);
0868   DensePropagator propDense(esDense, naviDet);
0869
0870   // Launch and collect results
0871   const auto& resultDense = propDense.propagate(sbtp, propOptsDense).value();
0872   const StepCollector::this_result& stepResultDense =
0873       resultDense.get<typename StepCollector::result_type>();
0874
0875   // Check the exit situation of the second volume
0876   for (unsigned int i = 0; i < stepResultDense.position.size(); i++) {
0877     if (1_m - stepResultDense.position[i].x() < 1e-4) {
0878       endParams = std::make_pair(stepResultDense.position[i],
0879                                  stepResultDense.momentum[i]);
0880       break;
0881     }
0882   }
0883   for (unsigned int i = 0; i < stepResult.position.size(); i++) {
0884     if (1_m - stepResult.position[i].x() < 1e-4) {
0885       endParamsControl =
0886           std::make_pair(stepResult.position[i], stepResult.momentum[i]);
0887       break;
0888     }
0889   }
0890
0891   CHECK_CLOSE_ABS(endParams.first, endParamsControl.first, 1_um);
0892   CHECK_CLOSE_ABS(endParams.second, endParamsControl.second, 1_um);
0893 }
0894
0895 // Test case a). The DenseEnvironmentExtension should state that it is not
0896 // valid in this case.
0897 BOOST_AUTO_TEST_CASE(step_extension_trackercalomdt_test) {
0898   double rotationAngle = std::numbers::pi / 2.;
0899   Vector3 xPos(cos(rotationAngle), 0., sin(rotationAngle));
0900   Vector3 yPos(0., 1., 0.);
0901   Vector3 zPos(-sin(rotationAngle), 0., cos(rotationAngle));
0902   MaterialSlab matProp(makeBeryllium(), 0.5_mm);
0903
0904   CuboidVolumeBuilder cvb;
0905   CuboidVolumeBuilder::SurfaceConfig sConf1;
0906   sConf1.position = Vector3(0.3_m, 0., 0.);
0907   sConf1.rotation.col(0) = xPos;
0908   sConf1.rotation.col(1) = yPos;
0909   sConf1.rotation.col(2) = zPos;
0910   sConf1.rBounds =
0911       std::make_shared<const RectangleBounds>(RectangleBounds(0.5_m, 0.5_m));
0912   sConf1.surMat = std::shared_ptr<const ISurfaceMaterial>(
0913       new HomogeneousSurfaceMaterial(matProp));
0914   sConf1.thickness = 1._mm;
0915   CuboidVolumeBuilder::LayerConfig lConf1;
0916   lConf1.surfaceCfg = {sConf1};
0917
0918   CuboidVolumeBuilder::SurfaceConfig sConf2;
0919   sConf2.position = Vector3(0.6_m, 0., 0.);
0920   sConf2.rotation.col(0) = xPos;
0921   sConf2.rotation.col(1) = yPos;
0922   sConf2.rotation.col(2) = zPos;
0923   sConf2.rBounds =
0924       std::make_shared<const RectangleBounds>(RectangleBounds(0.5_m, 0.5_m));
0925   sConf2.surMat = std::shared_ptr<const ISurfaceMaterial>(
0926       new HomogeneousSurfaceMaterial(matProp));
0927   sConf2.thickness = 1._mm;
0928   CuboidVolumeBuilder::LayerConfig lConf2;
0929   lConf2.surfaceCfg = {sConf2};
0930
0931   CuboidVolumeBuilder::VolumeConfig muConf1;
0932   muConf1.position = {2.3_m, 0., 0.};
0933   muConf1.length = {20._cm, 20._cm, 20._cm};
0934   muConf1.volumeMaterial =
0935       std::make_shared<HomogeneousVolumeMaterial>(makeBeryllium());
0936   muConf1.name = "MDT1";
0937   CuboidVolumeBuilder::VolumeConfig muConf2;
0938   muConf2.position = {2.7_m, 0., 0.};
0939   muConf2.length = {20._cm, 20._cm, 20._cm};
0940   muConf2.volumeMaterial =
0941       std::make_shared<HomogeneousVolumeMaterial>(makeBeryllium());
0942   muConf2.name = "MDT2";
0943
0944   CuboidVolumeBuilder::VolumeConfig vConf1;
0945   vConf1.position = {0.5_m, 0., 0.};
0946   vConf1.length = {1._m, 1._m, 1._m};
0947   vConf1.layerCfg = {lConf1, lConf2};
0948   vConf1.name = "Tracker";
0949   CuboidVolumeBuilder::VolumeConfig vConf2;
0950   vConf2.position = {1.5_m, 0., 0.};
0951   vConf2.length = {1._m, 1._m, 1._m};
0952   vConf2.volumeMaterial =
0953       std::make_shared<HomogeneousVolumeMaterial>(makeBeryllium());
0954   vConf2.name = "Calorimeter";
0955   CuboidVolumeBuilder::VolumeConfig vConf3;
0956   vConf3.position = {2.5_m, 0., 0.};
0957   vConf3.length = {1._m, 1._m, 1._m};
0958   vConf3.volumeCfg = {muConf1, muConf2};
0959   vConf3.name = "Muon system";
0960   CuboidVolumeBuilder::Config conf;
0961   conf.volumeCfg = {vConf1, vConf2, vConf3};
0962   conf.position = {1.5_m, 0., 0.};
0963   conf.length = {3._m, 1._m, 1._m};
0964
0965   // Build detector
0966   cvb.setConfig(conf);
0967   TrackingGeometryBuilder::Config tgbCfg;
0968   tgbCfg.trackingVolumeBuilders.push_back(
0969       [=](const auto& context, const auto& inner, const auto& vb) {
0970         return cvb.trackingVolume(context, inner, vb);
0971       });
0972   TrackingGeometryBuilder tgb(tgbCfg);
0973   std::shared_ptr<const TrackingGeometry> detector =
0974       tgb.trackingGeometry(tgContext);
0975
0976   // Build navigator
0977   Navigator naviVac({detector, true, true, true});
0978
0979   // Set initial parameters for the particle track
0980   BoundTrackParameters sbtp = BoundTrackParameters::createCurvilinear(
0981       Vector4::Zero(), 0_degree, 90_degree, 1_e / 1_GeV, Covariance::Identity(),
0982       ParticleHypothesis::pion());
0983
0984   using Stepper = EigenStepper<EigenStepperDenseExtension>;
0985   using Propagator = Acts::Propagator<Stepper, Navigator>;
0986   using PropagatorOptions = Propagator::Options<
0987       ActorList<StepCollector, MaterialInteractor, EndOfWorld>>;
0988
0989   // Set options for propagator
0990   PropagatorOptions propOpts(tgContext, mfContext);
0991   propOpts.actorList.get<EndOfWorld>().maxX = 3._m;
0992   propOpts.maxSteps = 10000;
0993
0994   // Build stepper and propagator
0995   auto bField = std::make_shared<ConstantBField>(Vector3(0., 0., 0.));
0996   Stepper es(bField);
0997   Propagator prop(es, naviVac);
0998
0999   // Launch and collect results
1000   const auto& result = prop.propagate(sbtp, propOpts).value();
1001   const StepCollector::this_result& stepResult =
1002       result.get<typename StepCollector::result_type>();
1003
1004   // Test that momentum changes only occurred at the right detector parts
1005   double lastMomentum = stepResult.momentum[0].x();
1006   for (unsigned int i = 0; i < stepResult.position.size(); i++) {
1007     // Test for changes
1008     if ((stepResult.position[i].x() > 0.3_m &&
1009          stepResult.position[i].x() < 0.6_m) ||
1010         (stepResult.position[i].x() > 0.6_m &&
1011          stepResult.position[i].x() <= 1._m) ||
1012         (stepResult.position[i].x() > 1._m &&
1013          stepResult.position[i].x() <= 2._m) ||
1014         (stepResult.position[i].x() > 2.2_m &&
1015          stepResult.position[i].x() <= 2.4_m) ||
1016         (stepResult.position[i].x() > 2.6_m &&
1017          stepResult.position[i].x() <= 2.8_m)) {
1018       BOOST_CHECK_LE(stepResult.momentum[i].x(), lastMomentum);
1019       lastMomentum = stepResult.momentum[i].x();
1020     } else
1021     // Test the absence of momentum loss
1022     {
1023       if (stepResult.position[i].x() < 0.3_m ||
1024           (stepResult.position[i].x() > 2._m &&
1025            stepResult.position[i].x() <= 2.2_m) ||
1026           (stepResult.position[i].x() > 2.4_m &&
1027            stepResult.position[i].x() <= 2.6_m) ||
1028           (stepResult.position[i].x() > 2.8_m &&
1029            stepResult.position[i].x() <= 3._m)) {
1030         BOOST_CHECK_EQUAL(stepResult.momentum[i].x(), lastMomentum);
1031       }
1032     }
1033   }
1034 }
1035
1036 }  // namespace Acts::Test