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