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