Core/Propagator/EigenStepperTests.cpp

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