Core/Propagator/SympyStepperTests.cpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 #include <boost/test/unit_test.hpp>
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Definitions/Direction.hpp"
0013 #include "Acts/Definitions/Tolerance.hpp"
0014 #include "Acts/Definitions/TrackParametrization.hpp"
0015 #include "Acts/Definitions/Units.hpp"
0016 #include "Acts/EventData/GenericBoundTrackParameters.hpp"
0017 #include "Acts/EventData/ParticleHypothesis.hpp"
0018 #include "Acts/EventData/TrackParameters.hpp"
0019 #include "Acts/EventData/TransformationHelpers.hpp"
0020 #include "Acts/Geometry/GeometryContext.hpp"
0021 #include "Acts/MagneticField/ConstantBField.hpp"
0022 #include "Acts/MagneticField/MagneticFieldContext.hpp"
0023 #include "Acts/MagneticField/MagneticFieldProvider.hpp"
0024 #include "Acts/Propagator/ConstrainedStep.hpp"
0025 #include "Acts/Propagator/SympyStepper.hpp"
0026 #include "Acts/Surfaces/BoundaryTolerance.hpp"
0027 #include "Acts/Surfaces/CurvilinearSurface.hpp"
0028 #include "Acts/Surfaces/PlaneSurface.hpp"
0029 #include "Acts/Utilities/Logger.hpp"
0030 #include "Acts/Utilities/Result.hpp"
0031 #include "ActsTests/CommonHelpers/FloatComparisons.hpp"
0032
0033 #include <cmath>
0034 #include <limits>
0035 #include <memory>
0036 #include <optional>
0037 #include <string>
0038 #include <type_traits>
0039 #include <utility>
0040 #include <vector>
0041
0042 using namespace Acts;
0043 using namespace Acts::UnitLiterals;
0044 using Acts::VectorHelpers::makeVector4;
0045
0046 namespace ActsTests {
0047
0048 using Covariance = BoundSquareMatrix;
0049
0050 static constexpr auto eps = 3 * std::numeric_limits<double>::epsilon();
0051
0052 // Create a test context
0053 GeometryContext tgContext = GeometryContext();
0054 MagneticFieldContext mfContext = MagneticFieldContext();
0055
0056 /// @brief Aborter for the case that a particle leaves the detector or reaches
0057 /// a custom made threshold.
0058 ///
0059 struct EndOfWorld {
0060   /// Maximum value in x-direction of the detector
0061   double maxX = 1_m;
0062
0063   /// @brief Constructor
0064   EndOfWorld() = default;
0065
0066   /// @brief Main call operator for the abort operation
0067   ///
0068   /// @tparam propagator_state_t State of the propagator
0069   /// @tparam stepper_t Type of the stepper
0070   /// @tparam navigator_t Type of the navigator
0071   ///
0072   /// @param [in] state State of the propagation
0073   /// @param [in] stepper Stepper of the propagation
0074   /// @param [in] navigator Navigator of the propagation
0075   ///
0076   /// @return Boolean statement if the particle is still in the detector
0077   template <typename propagator_state_t, typename stepper_t,
0078             typename navigator_t>
0079   bool operator()(propagator_state_t& state, const stepper_t& stepper,
0080                   const navigator_t& /*navigator*/,
0081                   const Logger& /*logger*/) const {
0082     const double tolerance = state.options.surfaceTolerance;
0083     if (maxX - std::abs(stepper.position(state.stepping).x()) <= tolerance ||
0084         std::abs(stepper.position(state.stepping).y()) >= 0.5_m ||
0085         std::abs(stepper.position(state.stepping).z()) >= 0.5_m) {
0086       return true;
0087     }
0088     return false;
0089   }
0090 };
0091
0092 ///
0093 /// @brief Data collector while propagation
0094 ///
0095 struct StepCollector {
0096   ///
0097   /// @brief Data container for result analysis
0098   ///
0099   struct this_result {
0100     // Position of the propagator after each step
0101     std::vector<Vector3> position;
0102     // Momentum of the propagator after each step
0103     std::vector<Vector3> momentum;
0104   };
0105
0106   using result_type = this_result;
0107
0108   /// @brief Main call operator for the action list. It stores the data for
0109   /// analysis afterwards
0110   ///
0111   /// @tparam propagator_state_t Type of the propagator state
0112   /// @tparam stepper_t Type of the stepper
0113   /// @tparam navigator_t Type of the navigator
0114   ///
0115   /// @param [in] state State of the propagator
0116   /// @param [in] stepper Stepper of the propagation
0117   /// @param [in] navigator Navigator of the propagation
0118   /// @param [out] result Struct which is filled with the data
0119   template <typename propagator_state_t, typename stepper_t,
0120             typename navigator_t>
0121   void operator()(propagator_state_t& state, const stepper_t& stepper,
0122                   const navigator_t& /*navigator*/, result_type& result,
0123                   const Logger& /*logger*/) const {
0124     result.position.push_back(stepper.position(state.stepping));
0125     result.momentum.push_back(stepper.momentum(state.stepping));
0126   }
0127 };
0128
0129 BOOST_AUTO_TEST_SUITE(PropagatorSuite)
0130
0131 /// These tests are aiming to test whether the state setup is working properly
0132 BOOST_AUTO_TEST_CASE(sympy_stepper_state_test) {
0133   // Set up some variables
0134   auto bField = std::make_shared<ConstantBField>(Vector3(1., 2.5, 33.33));
0135
0136   Vector3 pos(1., 2., 3.);
0137   Vector3 dir(4., 5., 6.);
0138   double time = 7.;
0139   double absMom = 8.;
0140   double charge = -1.;
0141
0142   SympyStepper::Options esOptions(tgContext, mfContext);
0143
0144   SympyStepper es(bField);
0145
0146   // Test charged parameters without covariance matrix
0147   BoundTrackParameters cp = BoundTrackParameters::createCurvilinear(
0148       makeVector4(pos, time), dir, charge / absMom, std::nullopt,
0149       ParticleHypothesis::pion());
0150   SympyStepper::State esState = es.makeState(esOptions);
0151   es.initialize(esState, cp);
0152
0153   // Test the result & compare with the input/test for reasonable members
0154   BOOST_CHECK_EQUAL(esState.jacToGlobal, BoundToFreeMatrix::Zero());
0155   BOOST_CHECK_EQUAL(esState.jacTransport, FreeMatrix::Identity());
0156   BOOST_CHECK_EQUAL(esState.derivative, FreeVector::Zero());
0157   BOOST_CHECK(!esState.covTransport);
0158   BOOST_CHECK_EQUAL(esState.cov, Covariance::Zero());
0159   BOOST_CHECK_EQUAL(esState.pathAccumulated, 0.);
0160   BOOST_CHECK_EQUAL(esState.previousStepSize, 0.);
0161
0162   // Test without charge and covariance matrix
0163   BoundTrackParameters ncp = BoundTrackParameters::createCurvilinear(
0164       makeVector4(pos, time), dir, 1 / absMom, std::nullopt,
0165       ParticleHypothesis::pion0());
0166   es.initialize(esState, ncp);
0167   BOOST_CHECK_EQUAL(es.charge(esState), 0.);
0168
0169   // Test with covariance matrix
0170   Covariance cov = 8. * Covariance::Identity();
0171   ncp = BoundTrackParameters::createCurvilinear(makeVector4(pos, time), dir,
0172                                                 1 / absMom, cov,
0173                                                 ParticleHypothesis::pion0());
0174   es.initialize(esState, ncp);
0175   BOOST_CHECK_NE(esState.jacToGlobal, BoundToFreeMatrix::Zero());
0176   BOOST_CHECK(esState.covTransport);
0177   BOOST_CHECK_EQUAL(esState.cov, cov);
0178 }
0179
0180 /// These tests are aiming to test the functions of the SympyStepper
0181 /// The numerical correctness of the stepper is tested in the integration tests
0182 BOOST_AUTO_TEST_CASE(sympy_stepper_test) {
0183   // Set up some variables for the state
0184   Direction navDir = Direction::Backward();
0185   double stepSize = 123.;
0186   auto bField = std::make_shared<ConstantBField>(Vector3(1., 2.5, 33.33));
0187   auto bCache = bField->makeCache(mfContext);
0188
0189   // Construct the parameters
0190   Vector3 pos(1., 2., 3.);
0191   Vector3 dir = Vector3(4., 5., 6.).normalized();
0192   double time = 7.;
0193   double absMom = 8.;
0194   double charge = -1.;
0195   Covariance cov = 8. * Covariance::Identity();
0196   BoundTrackParameters cp = BoundTrackParameters::createCurvilinear(
0197       makeVector4(pos, time), dir, charge / absMom, cov,
0198       ParticleHypothesis::pion());
0199
0200   SympyStepper::Options esOptions(tgContext, mfContext);
0201   esOptions.maxStepSize = stepSize;
0202   esOptions.initialStepSize = 10_m;
0203
0204   // Build the stepper and the state
0205   SympyStepper es(bField);
0206   SympyStepper::State esState = es.makeState(esOptions);
0207   es.initialize(esState, cp);
0208
0209   // Test the getters
0210   CHECK_CLOSE_ABS(es.position(esState), pos, eps);
0211   CHECK_CLOSE_ABS(es.direction(esState), dir, eps);
0212   CHECK_CLOSE_ABS(es.absoluteMomentum(esState), absMom, eps);
0213   CHECK_CLOSE_ABS(es.charge(esState), charge, eps);
0214   CHECK_CLOSE_ABS(es.time(esState), time, eps);
0215   BOOST_CHECK_EQUAL(es.getField(esState, pos).value(),
0216                     bField->getField(pos, bCache).value());
0217
0218   // Step size modifies
0219   const std::string originalStepSize = esState.stepSize.toString();
0220
0221   es.updateStepSize(esState, -1337., ConstrainedStep::Type::Navigator);
0222   BOOST_CHECK_EQUAL(esState.previousStepSize, stepSize);
0223   BOOST_CHECK_EQUAL(esState.stepSize.value(), -1337.);
0224
0225   es.releaseStepSize(esState, ConstrainedStep::Type::Navigator);
0226   BOOST_CHECK_EQUAL(esState.stepSize.value(), stepSize);
0227   BOOST_CHECK_EQUAL(es.outputStepSize(esState), originalStepSize);
0228
0229   // Test the curvilinear state construction
0230   auto curvState = es.curvilinearState(esState);
0231   auto curvPars = std::get<0>(curvState);
0232   CHECK_CLOSE_ABS(curvPars.position(tgContext), cp.position(tgContext), eps);
0233   CHECK_CLOSE_ABS(curvPars.momentum(), cp.momentum(), 10e-6);
0234   CHECK_CLOSE_ABS(curvPars.charge(), cp.charge(), eps);
0235   CHECK_CLOSE_ABS(curvPars.time(), cp.time(), eps);
0236   BOOST_CHECK(curvPars.covariance().has_value());
0237   BOOST_CHECK_NE(*curvPars.covariance(), cov);
0238   CHECK_CLOSE_COVARIANCE(std::get<1>(curvState),
0239                          BoundMatrix(BoundMatrix::Identity()), eps);
0240   CHECK_CLOSE_ABS(std::get<2>(curvState), 0., eps);
0241
0242   // Test the update method
0243   Vector3 newPos(2., 4., 8.);
0244   Vector3 newMom(3., 9., 27.);
0245   double newTime(321.);
0246   es.update(esState, newPos, newMom.normalized(), charge / newMom.norm(),
0247             newTime);
0248   BOOST_CHECK_EQUAL(es.position(esState), newPos);
0249   BOOST_CHECK_EQUAL(es.direction(esState), newMom.normalized());
0250   BOOST_CHECK_EQUAL(es.absoluteMomentum(esState), newMom.norm());
0251   BOOST_CHECK_EQUAL(es.charge(esState), charge);
0252   BOOST_CHECK_EQUAL(es.time(esState), newTime);
0253
0254   // The covariance transport
0255   esState.cov = cov;
0256   es.transportCovarianceToCurvilinear(esState);
0257   BOOST_CHECK_NE(esState.cov, cov);
0258   BOOST_CHECK_NE(esState.jacToGlobal, BoundToFreeMatrix::Zero());
0259   BOOST_CHECK_EQUAL(esState.jacTransport, FreeMatrix::Identity());
0260   BOOST_CHECK_EQUAL(esState.derivative, FreeVector::Zero());
0261
0262   // Perform a step without and with covariance transport
0263   esState.cov = cov;
0264
0265   esState.covTransport = false;
0266   es.step(esState, navDir, nullptr).value();
0267   CHECK_CLOSE_COVARIANCE(esState.cov, cov, eps);
0268   BOOST_CHECK_NE(es.position(esState).norm(), newPos.norm());
0269   BOOST_CHECK_NE(es.direction(esState), newMom.normalized());
0270   BOOST_CHECK_EQUAL(es.charge(esState), charge);
0271   BOOST_CHECK_LT(es.time(esState), newTime);
0272   BOOST_CHECK_EQUAL(esState.derivative, FreeVector::Zero());
0273   BOOST_CHECK_EQUAL(esState.jacTransport, FreeMatrix::Identity());
0274
0275   esState.covTransport = true;
0276   es.step(esState, navDir, nullptr).value();
0277   CHECK_CLOSE_COVARIANCE(esState.cov, cov, eps);
0278   BOOST_CHECK_NE(es.position(esState).norm(), newPos.norm());
0279   BOOST_CHECK_NE(es.direction(esState), newMom.normalized());
0280   BOOST_CHECK_EQUAL(es.charge(esState), charge);
0281   BOOST_CHECK_LT(es.time(esState), newTime);
0282   BOOST_CHECK_NE(esState.derivative, FreeVector::Zero());
0283   BOOST_CHECK_NE(esState.jacTransport, FreeMatrix::Identity());
0284
0285   /// Test the state reset
0286   // Construct the parameters
0287   Vector3 pos2(1.5, -2.5, 3.5);
0288   Vector3 dir2 = Vector3(4.5, -5.5, 6.5).normalized();
0289   double time2 = 7.5;
0290   double absMom2 = 8.5;
0291   double charge2 = 1.;
0292   BoundSquareMatrix cov2 = 8.5 * Covariance::Identity();
0293   BoundTrackParameters cp2 = BoundTrackParameters::createCurvilinear(
0294       makeVector4(pos2, time2), dir2, charge2 / absMom2, cov2,
0295       ParticleHypothesis::pion());
0296   FreeVector freeParams = transformBoundToFreeParameters(
0297       cp2.referenceSurface(), tgContext, cp2.parameters());
0298   navDir = Direction::Forward();
0299
0300   auto copyState = [&](auto& field, const auto& state) {
0301     using field_t = std::decay_t<decltype(field)>;
0302     std::decay_t<decltype(state)> copy = es.makeState(esOptions);
0303     es.initialize(esState, cp);
0304     copy.pars = state.pars;
0305     copy.covTransport = state.covTransport;
0306     copy.cov = state.cov;
0307     copy.jacobian = state.jacobian;
0308     copy.jacToGlobal = state.jacToGlobal;
0309     copy.jacTransport = state.jacTransport;
0310     copy.derivative = state.derivative;
0311     copy.pathAccumulated = state.pathAccumulated;
0312     copy.stepSize = state.stepSize;
0313     copy.previousStepSize = state.previousStepSize;
0314
0315     copy.fieldCache = MagneticFieldProvider::Cache(
0316         std::in_place_type<typename field_t::Cache>,
0317         state.fieldCache.template as<typename field_t::Cache>());
0318
0319     return copy;
0320   };
0321
0322   // Reset all possible parameters
0323   SympyStepper::State esStateCopy = copyState(*bField, esState);
0324   BOOST_CHECK(cp2.covariance().has_value());
0325   es.initialize(esStateCopy, cp2.parameters(), *cp2.covariance(),
0326                 cp2.particleHypothesis(), cp2.referenceSurface());
0327   // Test all components
0328   BOOST_CHECK_NE(esStateCopy.jacToGlobal, BoundToFreeMatrix::Zero());
0329   BOOST_CHECK_NE(esStateCopy.jacToGlobal, esState.jacToGlobal);
0330   BOOST_CHECK_EQUAL(esStateCopy.jacTransport, FreeMatrix::Identity());
0331   BOOST_CHECK_EQUAL(esStateCopy.derivative, FreeVector::Zero());
0332   BOOST_CHECK(esStateCopy.covTransport);
0333   BOOST_CHECK_EQUAL(esStateCopy.cov, cov2);
0334   BOOST_CHECK_EQUAL(es.position(esStateCopy),
0335                     freeParams.template segment<3>(eFreePos0));
0336   BOOST_CHECK_EQUAL(es.direction(esStateCopy),
0337                     freeParams.template segment<3>(eFreeDir0).normalized());
0338   BOOST_CHECK_EQUAL(es.absoluteMomentum(esStateCopy),
0339                     std::abs(1. / freeParams[eFreeQOverP]));
0340   BOOST_CHECK_EQUAL(es.charge(esStateCopy), -es.charge(esState));
0341   BOOST_CHECK_EQUAL(es.time(esStateCopy), freeParams[eFreeTime]);
0342   BOOST_CHECK_EQUAL(esStateCopy.pathAccumulated, 0.);
0343   BOOST_CHECK_EQUAL(esStateCopy.stepSize.value(), stepSize);
0344   BOOST_CHECK_EQUAL(esStateCopy.previousStepSize, 0.);
0345
0346   /// Repeat with surface related methods
0347   std::shared_ptr<PlaneSurface> plane =
0348       CurvilinearSurface(pos, dir.normalized()).planeSurface();
0349   auto bp = BoundTrackParameters::create(
0350                 tgContext, plane, makeVector4(pos, time), dir, charge / absMom,
0351                 cov, ParticleHypothesis::pion())
0352                 .value();
0353   esOptions = SympyStepper::Options(tgContext, mfContext);
0354   esState = es.makeState(esOptions);
0355   es.initialize(esState, bp);
0356
0357   // Test the intersection in the context of a surface
0358   std::shared_ptr<PlaneSurface> targetSurface =
0359       CurvilinearSurface(pos + navDir * 2. * dir, dir).planeSurface();
0360   es.updateSurfaceStatus(esState, *targetSurface, 0, navDir,
0361                          BoundaryTolerance::None(), s_onSurfaceTolerance,
0362                          ConstrainedStep::Type::Navigator);
0363   CHECK_CLOSE_ABS(esState.stepSize.value(ConstrainedStep::Type::Navigator),
0364                   navDir * 2., eps);
0365
0366   const auto getNavigationTarget = [&](const Surface& s,
0367                                        const BoundaryTolerance& bt) {
0368     auto [intersection, intersectionIndex] =
0369         s.intersect(tgContext, es.position(esState),
0370                     navDir * es.direction(esState), bt)
0371             .closestWithIndex();
0372     return NavigationTarget(intersection, intersectionIndex, s, bt);
0373   };
0374
0375   // Test the step size modification in the context of a surface
0376   es.updateStepSize(
0377       esState, getNavigationTarget(*targetSurface, BoundaryTolerance::None()),
0378       navDir, ConstrainedStep::Type::Navigator);
0379   CHECK_CLOSE_ABS(esState.stepSize.value(), 2., eps);
0380   esState.stepSize.setUser(navDir * stepSize);
0381   es.releaseStepSize(esState, ConstrainedStep::Type::Navigator);
0382   es.updateStepSize(
0383       esState, getNavigationTarget(*targetSurface, BoundaryTolerance::None()),
0384       navDir, ConstrainedStep::Type::Navigator);
0385   CHECK_CLOSE_ABS(esState.stepSize.value(), 2., eps);
0386
0387   // Test the bound state construction
0388   auto boundState = es.boundState(esState, *plane).value();
0389   auto boundPars = std::get<0>(boundState);
0390   CHECK_CLOSE_ABS(boundPars.position(tgContext), bp.position(tgContext), eps);
0391   CHECK_CLOSE_ABS(boundPars.momentum(), bp.momentum(), 1e-7);
0392   CHECK_CLOSE_ABS(boundPars.charge(), bp.charge(), eps);
0393   CHECK_CLOSE_ABS(boundPars.time(), bp.time(), eps);
0394   BOOST_CHECK(boundPars.covariance().has_value());
0395   BOOST_CHECK_NE(*boundPars.covariance(), cov);
0396   CHECK_CLOSE_COVARIANCE(std::get<1>(boundState),
0397                          BoundMatrix(BoundMatrix::Identity()), eps);
0398   CHECK_CLOSE_ABS(std::get<2>(boundState), 0., eps);
0399
0400   // Transport the covariance in the context of a surface
0401   es.transportCovarianceToBound(esState, *plane);
0402   BOOST_CHECK_NE(esState.cov, cov);
0403   BOOST_CHECK_NE(esState.jacToGlobal, BoundToFreeMatrix::Zero());
0404   BOOST_CHECK_EQUAL(esState.jacTransport, FreeMatrix::Identity());
0405   BOOST_CHECK_EQUAL(esState.derivative, FreeVector::Zero());
0406
0407   // Update in context of a surface
0408   freeParams = transformBoundToFreeParameters(bp.referenceSurface(), tgContext,
0409                                               bp.parameters());
0410
0411   es.update(esState, freeParams, bp.parameters(), 2 * (*bp.covariance()),
0412             *plane);
0413   CHECK_CLOSE_OR_SMALL(es.position(esState), pos, eps, eps);
0414   CHECK_CLOSE_OR_SMALL(es.direction(esState), dir, eps, eps);
0415   CHECK_CLOSE_REL(es.absoluteMomentum(esState), absMom, eps);
0416   BOOST_CHECK_EQUAL(es.charge(esState), charge);
0417   CHECK_CLOSE_OR_SMALL(es.time(esState), time, eps, eps);
0418   CHECK_CLOSE_COVARIANCE(esState.cov, Covariance(2. * cov), eps);
0419
0420   // Test a case where no step size adjustment is required
0421   esState.options.stepTolerance = 2. * 4.4258e+09;
0422   double h0 = esState.stepSize.value();
0423   es.step(esState, Direction::Forward(), nullptr);
0424   CHECK_CLOSE_ABS(h0, esState.stepSize.value(), eps);
0425 }
0426
0427 BOOST_AUTO_TEST_SUITE_END()
0428
0429 }  // namespace ActsTests