Core/Propagator/EigenStepperTests.cpp

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