Core/Propagator/StraightLineStepperTests.cpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 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/EventData/Charge.hpp"
0015 #include "Acts/EventData/GenericBoundTrackParameters.hpp"
0016 #include "Acts/EventData/GenericCurvilinearTrackParameters.hpp"
0017 #include "Acts/EventData/TrackParameters.hpp"
0018 #include "Acts/EventData/TransformationHelpers.hpp"
0019 #include "Acts/EventData/detail/CorrectedTransformationFreeToBound.hpp"
0020 #include "Acts/Geometry/GeometryContext.hpp"
0021 #include "Acts/MagneticField/MagneticFieldContext.hpp"
0022 #include "Acts/Propagator/ConstrainedStep.hpp"
0023 #include "Acts/Propagator/StraightLineStepper.hpp"
0024 #include "Acts/Surfaces/BoundaryCheck.hpp"
0025 #include "Acts/Surfaces/PlaneSurface.hpp"
0026 #include "Acts/Surfaces/Surface.hpp"
0027 #include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
0028 #include "Acts/Utilities/Helpers.hpp"
0029 #include "Acts/Utilities/Result.hpp"
0030
0031 #include <functional>
0032 #include <limits>
0033 #include <memory>
0034 #include <optional>
0035 #include <string>
0036 #include <tuple>
0037 #include <utility>
0038
0039 using Acts::VectorHelpers::makeVector4;
0040
0041 namespace Acts::Test {
0042
0043 using Covariance = BoundSquareMatrix;
0044
0045 /// @brief Simplified propagator state
0046 struct PropState {
0047   /// @brief Constructor
0048   PropState(Direction direction, StraightLineStepper::State sState)
0049       : stepping(std::move(sState)) {
0050     options.direction = direction;
0051   }
0052   /// State of the straight line stepper
0053   StraightLineStepper::State stepping;
0054   /// Propagator options which only carry the particle's mass
0055   struct {
0056     Direction direction = Direction::Forward;
0057   } options;
0058 };
0059
0060 struct MockNavigator {};
0061
0062 static constexpr MockNavigator mockNavigator;
0063
0064 static constexpr auto eps = 2 * std::numeric_limits<double>::epsilon();
0065
0066 /// These tests are aiming to test whether the state setup is working properly
0067 BOOST_AUTO_TEST_CASE(straight_line_stepper_state_test) {
0068   // Set up some variables
0069   GeometryContext tgContext = GeometryContext();
0070   MagneticFieldContext mfContext = MagneticFieldContext();
0071   double stepSize = 123.;
0072   double tolerance = 234.;
0073
0074   Vector3 pos(1., 2., 3.);
0075   Vector3 dir(4., 5., 6.);
0076   double time = 7.;
0077   double absMom = 8.;
0078   double charge = -1.;
0079
0080   // Test charged parameters without covariance matrix
0081   CurvilinearTrackParameters cp(makeVector4(pos, time), dir, charge / absMom,
0082                                 std::nullopt, ParticleHypothesis::pion());
0083   StraightLineStepper::State slsState(tgContext, mfContext, cp, stepSize,
0084                                       tolerance);
0085
0086   StraightLineStepper sls;
0087
0088   // Test the result & compare with the input/test for reasonable members
0089   BOOST_CHECK_EQUAL(slsState.jacToGlobal, BoundToFreeMatrix::Zero());
0090   BOOST_CHECK_EQUAL(slsState.jacTransport, FreeMatrix::Identity());
0091   BOOST_CHECK_EQUAL(slsState.derivative, FreeVector::Zero());
0092   BOOST_CHECK(!slsState.covTransport);
0093   BOOST_CHECK_EQUAL(slsState.cov, Covariance::Zero());
0094   CHECK_CLOSE_OR_SMALL(sls.position(slsState), pos, eps, eps);
0095   CHECK_CLOSE_OR_SMALL(sls.direction(slsState), dir.normalized(), eps, eps);
0096   CHECK_CLOSE_REL(sls.absoluteMomentum(slsState), absMom, eps);
0097   BOOST_CHECK_EQUAL(sls.charge(slsState), charge);
0098   CHECK_CLOSE_OR_SMALL(sls.time(slsState), time, eps, eps);
0099   BOOST_CHECK_EQUAL(slsState.pathAccumulated, 0.);
0100   BOOST_CHECK_EQUAL(slsState.stepSize.value(), stepSize);
0101   BOOST_CHECK_EQUAL(slsState.previousStepSize, 0.);
0102   BOOST_CHECK_EQUAL(slsState.tolerance, tolerance);
0103
0104   // Test without charge and covariance matrix
0105   CurvilinearTrackParameters ncp(makeVector4(pos, time), dir, 1 / absMom,
0106                                  std::nullopt, ParticleHypothesis::pion0());
0107   slsState = StraightLineStepper::State(tgContext, mfContext, ncp, stepSize,
0108                                         tolerance);
0109
0110   // Test with covariance matrix
0111   Covariance cov = 8. * Covariance::Identity();
0112   ncp = CurvilinearTrackParameters(makeVector4(pos, time), dir, 1 / absMom, cov,
0113                                    ParticleHypothesis::pion0());
0114   slsState = StraightLineStepper::State(tgContext, mfContext, ncp, stepSize,
0115                                         tolerance);
0116   BOOST_CHECK_NE(slsState.jacToGlobal, BoundToFreeMatrix::Zero());
0117   BOOST_CHECK(slsState.covTransport);
0118   BOOST_CHECK_EQUAL(slsState.cov, cov);
0119 }
0120
0121 /// These tests are aiming to test the functions of the StraightLineStepper
0122 /// The numerical correctness of the stepper is tested in the integration tests
0123 BOOST_AUTO_TEST_CASE(straight_line_stepper_test) {
0124   // Set up some variables for the state
0125   GeometryContext tgContext = GeometryContext();
0126   MagneticFieldContext mfContext = MagneticFieldContext();
0127   Direction navDir = Direction::Backward;
0128   double stepSize = 123.;
0129   double tolerance = 234.;
0130
0131   // Construct the parameters
0132   Vector3 pos(1., 2., 3.);
0133   Vector3 dir = Vector3(4., 5., 6.).normalized();
0134   double time = 7.;
0135   double absMom = 8.;
0136   double charge = -1.;
0137   Covariance cov = 8. * Covariance::Identity();
0138   CurvilinearTrackParameters cp(makeVector4(pos, time), dir, charge / absMom,
0139                                 cov, ParticleHypothesis::pion());
0140
0141   // Build the state and the stepper
0142   StraightLineStepper::State slsState(tgContext, mfContext, cp, stepSize,
0143                                       tolerance);
0144   StraightLineStepper sls;
0145
0146   // Test the getters
0147   CHECK_CLOSE_ABS(sls.position(slsState), pos, 1e-6);
0148   CHECK_CLOSE_ABS(sls.direction(slsState), dir, 1e-6);
0149   CHECK_CLOSE_ABS(sls.absoluteMomentum(slsState), absMom, 1e-6);
0150   BOOST_CHECK_EQUAL(sls.charge(slsState), charge);
0151   BOOST_CHECK_EQUAL(sls.time(slsState), time);
0152
0153   //~ BOOST_CHECK_EQUAL(sls.overstepLimit(slsState), tolerance);
0154
0155   // Step size modifies
0156   const std::string originalStepSize = slsState.stepSize.toString();
0157
0158   sls.updateStepSize(slsState, -1337., ConstrainedStep::actor);
0159   BOOST_CHECK_EQUAL(slsState.previousStepSize, stepSize);
0160   BOOST_CHECK_EQUAL(slsState.stepSize.value(), -1337.);
0161
0162   sls.releaseStepSize(slsState, ConstrainedStep::actor);
0163   BOOST_CHECK_EQUAL(slsState.stepSize.value(), stepSize);
0164   BOOST_CHECK_EQUAL(sls.outputStepSize(slsState), originalStepSize);
0165
0166   // Test the curvilinear state construction
0167   auto curvState = sls.curvilinearState(slsState);
0168   auto curvPars = std::get<0>(curvState);
0169   CHECK_CLOSE_ABS(curvPars.position(tgContext), cp.position(tgContext), 1e-6);
0170   CHECK_CLOSE_ABS(curvPars.absoluteMomentum(), cp.absoluteMomentum(), 1e-6);
0171   CHECK_CLOSE_ABS(curvPars.charge(), cp.charge(), 1e-6);
0172   CHECK_CLOSE_ABS(curvPars.time(), cp.time(), 1e-6);
0173   BOOST_CHECK(curvPars.covariance().has_value());
0174   BOOST_CHECK_NE(*curvPars.covariance(), cov);
0175   CHECK_CLOSE_COVARIANCE(std::get<1>(curvState),
0176                          BoundMatrix(BoundMatrix::Identity()), 1e-6);
0177   CHECK_CLOSE_ABS(std::get<2>(curvState), 0., 1e-6);
0178
0179   // Test the update method
0180   Vector3 newPos(2., 4., 8.);
0181   Vector3 newMom(3., 9., 27.);
0182   double newTime(321.);
0183   sls.update(slsState, newPos, newMom.normalized(), charge / newMom.norm(),
0184              newTime);
0185   CHECK_CLOSE_ABS(sls.position(slsState), newPos, 1e-6);
0186   CHECK_CLOSE_ABS(sls.direction(slsState), newMom.normalized(), 1e-6);
0187   CHECK_CLOSE_ABS(sls.absoluteMomentum(slsState), newMom.norm(), 1e-6);
0188   BOOST_CHECK_EQUAL(sls.charge(slsState), charge);
0189   BOOST_CHECK_EQUAL(sls.time(slsState), newTime);
0190
0191   // The covariance transport
0192   slsState.cov = cov;
0193   sls.transportCovarianceToCurvilinear(slsState);
0194   BOOST_CHECK_NE(slsState.cov, cov);
0195   BOOST_CHECK_NE(slsState.jacToGlobal, BoundToFreeMatrix::Zero());
0196   BOOST_CHECK_EQUAL(slsState.jacTransport, FreeMatrix::Identity());
0197   BOOST_CHECK_EQUAL(slsState.derivative, FreeVector::Zero());
0198
0199   // Perform a step without and with covariance transport
0200   slsState.cov = cov;
0201   PropState ps(navDir, slsState);
0202
0203   ps.stepping.covTransport = false;
0204   double h = sls.step(ps, mockNavigator).value();
0205   BOOST_CHECK_EQUAL(ps.stepping.stepSize.value(), stepSize);
0206   BOOST_CHECK_EQUAL(ps.stepping.stepSize.value(), h * navDir);
0207   CHECK_CLOSE_COVARIANCE(ps.stepping.cov, cov, 1e-6);
0208   BOOST_CHECK_GT(sls.position(ps.stepping).norm(), newPos.norm());
0209   CHECK_CLOSE_ABS(sls.direction(ps.stepping), newMom.normalized(), 1e-6);
0210   CHECK_CLOSE_ABS(sls.absoluteMomentum(ps.stepping), newMom.norm(), 1e-6);
0211   CHECK_CLOSE_ABS(sls.charge(ps.stepping), charge, 1e-6);
0212   BOOST_CHECK_LT(sls.time(ps.stepping), newTime);
0213   BOOST_CHECK_EQUAL(ps.stepping.derivative, FreeVector::Zero());
0214   BOOST_CHECK_EQUAL(ps.stepping.jacTransport, FreeMatrix::Identity());
0215
0216   ps.stepping.covTransport = true;
0217   double h2 = sls.step(ps, mockNavigator).value();
0218   BOOST_CHECK_EQUAL(ps.stepping.stepSize.value(), stepSize);
0219   BOOST_CHECK_EQUAL(h2, h);
0220   CHECK_CLOSE_COVARIANCE(ps.stepping.cov, cov, 1e-6);
0221   BOOST_CHECK_GT(sls.position(ps.stepping).norm(), newPos.norm());
0222   CHECK_CLOSE_ABS(sls.direction(ps.stepping), newMom.normalized(), 1e-6);
0223   CHECK_CLOSE_ABS(sls.absoluteMomentum(ps.stepping), newMom.norm(), 1e-6);
0224   CHECK_CLOSE_ABS(sls.charge(ps.stepping), charge, 1e-6);
0225   BOOST_CHECK_LT(sls.time(ps.stepping), newTime);
0226   BOOST_CHECK_NE(ps.stepping.derivative, FreeVector::Zero());
0227   BOOST_CHECK_NE(ps.stepping.jacTransport, FreeMatrix::Identity());
0228
0229   /// Test the state reset
0230   // Construct the parameters
0231   Vector3 pos2(1.5, -2.5, 3.5);
0232   Vector3 dir2 = Vector3(4.5, -5.5, 6.5).normalized();
0233   double time2 = 7.5;
0234   double absMom2 = 8.5;
0235   double charge2 = 1.;
0236   BoundSquareMatrix cov2 = 8.5 * Covariance::Identity();
0237   CurvilinearTrackParameters cp2(makeVector4(pos2, time2), dir2,
0238                                  charge2 / absMom2, cov2,
0239                                  ParticleHypothesis::pion());
0240   BOOST_CHECK(cp2.covariance().has_value());
0241   FreeVector freeParams = transformBoundToFreeParameters(
0242       cp2.referenceSurface(), tgContext, cp2.parameters());
0243   navDir = Direction::Forward;
0244   double stepSize2 = -2. * stepSize;
0245
0246   // Reset all possible parameters
0247   StraightLineStepper::State slsStateCopy(ps.stepping);
0248   sls.resetState(slsStateCopy, cp2.parameters(), *cp2.covariance(),
0249                  cp2.referenceSurface(), stepSize2);
0250   // Test all components
0251   BOOST_CHECK_NE(slsStateCopy.jacToGlobal, BoundToFreeMatrix::Zero());
0252   BOOST_CHECK_NE(slsStateCopy.jacToGlobal, ps.stepping.jacToGlobal);
0253   BOOST_CHECK_EQUAL(slsStateCopy.jacTransport, FreeMatrix::Identity());
0254   BOOST_CHECK_EQUAL(slsStateCopy.derivative, FreeVector::Zero());
0255   BOOST_CHECK(slsStateCopy.covTransport);
0256   BOOST_CHECK_EQUAL(slsStateCopy.cov, cov2);
0257   CHECK_CLOSE_ABS(sls.position(slsStateCopy),
0258                   freeParams.template segment<3>(eFreePos0), 1e-6);
0259   CHECK_CLOSE_ABS(sls.direction(slsStateCopy),
0260                   freeParams.template segment<3>(eFreeDir0).normalized(), 1e-6);
0261   CHECK_CLOSE_ABS(sls.absoluteMomentum(slsStateCopy),
0262                   std::abs(1. / freeParams[eFreeQOverP]), 1e-6);
0263   CHECK_CLOSE_ABS(sls.charge(slsStateCopy), -sls.charge(ps.stepping), 1e-6);
0264   CHECK_CLOSE_ABS(sls.time(slsStateCopy), freeParams[eFreeTime], 1e-6);
0265   BOOST_CHECK_EQUAL(slsStateCopy.pathAccumulated, 0.);
0266   BOOST_CHECK_EQUAL(slsStateCopy.stepSize.value(), navDir * stepSize2);
0267   BOOST_CHECK_EQUAL(slsStateCopy.previousStepSize,
0268                     ps.stepping.previousStepSize);
0269   BOOST_CHECK_EQUAL(slsStateCopy.tolerance, ps.stepping.tolerance);
0270
0271   // Reset all possible parameters except the step size
0272   slsStateCopy = ps.stepping;
0273   sls.resetState(slsStateCopy, cp2.parameters(), *cp2.covariance(),
0274                  cp2.referenceSurface());
0275   // Test all components
0276   BOOST_CHECK_NE(slsStateCopy.jacToGlobal, BoundToFreeMatrix::Zero());
0277   BOOST_CHECK_NE(slsStateCopy.jacToGlobal, ps.stepping.jacToGlobal);
0278   BOOST_CHECK_EQUAL(slsStateCopy.jacTransport, FreeMatrix::Identity());
0279   BOOST_CHECK_EQUAL(slsStateCopy.derivative, FreeVector::Zero());
0280   BOOST_CHECK(slsStateCopy.covTransport);
0281   BOOST_CHECK_EQUAL(slsStateCopy.cov, cov2);
0282   CHECK_CLOSE_ABS(sls.position(slsStateCopy),
0283                   freeParams.template segment<3>(eFreePos0), 1e-6);
0284   CHECK_CLOSE_ABS(sls.direction(slsStateCopy),
0285                   freeParams.template segment<3>(eFreeDir0), 1e-6);
0286   CHECK_CLOSE_ABS(sls.absoluteMomentum(slsStateCopy),
0287                   std::abs(1. / freeParams[eFreeQOverP]), 1e-6);
0288   CHECK_CLOSE_ABS(sls.charge(slsStateCopy), -sls.charge(ps.stepping), 1e-6);
0289   CHECK_CLOSE_ABS(sls.time(slsStateCopy), freeParams[eFreeTime], 1e-6);
0290   BOOST_CHECK_EQUAL(slsStateCopy.pathAccumulated, 0.);
0291   BOOST_CHECK_EQUAL(slsStateCopy.stepSize.value(),
0292                     std::numeric_limits<double>::max());
0293   BOOST_CHECK_EQUAL(slsStateCopy.previousStepSize,
0294                     ps.stepping.previousStepSize);
0295   BOOST_CHECK_EQUAL(slsStateCopy.tolerance, ps.stepping.tolerance);
0296
0297   // Reset the least amount of parameters
0298   slsStateCopy = ps.stepping;
0299   sls.resetState(slsStateCopy, cp2.parameters(), *cp2.covariance(),
0300                  cp2.referenceSurface());
0301   // Test all components
0302   BOOST_CHECK_NE(slsStateCopy.jacToGlobal, BoundToFreeMatrix::Zero());
0303   BOOST_CHECK_NE(slsStateCopy.jacToGlobal, ps.stepping.jacToGlobal);
0304   BOOST_CHECK_EQUAL(slsStateCopy.jacTransport, FreeMatrix::Identity());
0305   BOOST_CHECK_EQUAL(slsStateCopy.derivative, FreeVector::Zero());
0306   BOOST_CHECK(slsStateCopy.covTransport);
0307   BOOST_CHECK_EQUAL(slsStateCopy.cov, cov2);
0308   CHECK_CLOSE_ABS(sls.position(slsStateCopy),
0309                   freeParams.template segment<3>(eFreePos0), 1e-6);
0310   CHECK_CLOSE_ABS(sls.direction(slsStateCopy),
0311                   freeParams.template segment<3>(eFreeDir0).normalized(), 1e-6);
0312   CHECK_CLOSE_ABS(sls.absoluteMomentum(slsStateCopy),
0313                   std::abs(1. / freeParams[eFreeQOverP]), 1e-6);
0314   CHECK_CLOSE_ABS(sls.charge(slsStateCopy), -sls.charge(ps.stepping), 1e-6);
0315   CHECK_CLOSE_ABS(sls.time(slsStateCopy), freeParams[eFreeTime], 1e-6);
0316   BOOST_CHECK_EQUAL(slsStateCopy.pathAccumulated, 0.);
0317   BOOST_CHECK_EQUAL(slsStateCopy.stepSize.value(),
0318                     std::numeric_limits<double>::max());
0319   BOOST_CHECK_EQUAL(slsStateCopy.previousStepSize,
0320                     ps.stepping.previousStepSize);
0321   BOOST_CHECK_EQUAL(slsStateCopy.tolerance, ps.stepping.tolerance);
0322
0323   /// Repeat with surface related methods
0324   auto plane = Surface::makeShared<PlaneSurface>(pos, dir);
0325   auto bp = BoundTrackParameters::create(
0326                 plane, tgContext, makeVector4(pos, time), dir, charge / absMom,
0327                 cov, ParticleHypothesis::pion())
0328                 .value();
0329   slsState =
0330       StraightLineStepper::State(tgContext, mfContext, cp, stepSize, tolerance);
0331
0332   // Test the intersection in the context of a surface
0333   auto targetSurface =
0334       Surface::makeShared<PlaneSurface>(pos + navDir * 2. * dir, dir);
0335   sls.updateSurfaceStatus(slsState, *targetSurface, 0, navDir,
0336                           BoundaryCheck(false));
0337   CHECK_CLOSE_ABS(slsState.stepSize.value(ConstrainedStep::actor), navDir * 2.,
0338                   1e-6);
0339
0340   // Test the step size modification in the context of a surface
0341   sls.updateStepSize(
0342       slsState,
0343       targetSurface
0344           ->intersect(slsState.geoContext, sls.position(slsState),
0345                       navDir * sls.direction(slsState), BoundaryCheck(false))
0346           .closest(),
0347       navDir, false);
0348   CHECK_CLOSE_ABS(slsState.stepSize.value(), 2, 1e-6);
0349   slsState.stepSize.setUser(navDir * stepSize);
0350   sls.updateStepSize(
0351       slsState,
0352       targetSurface
0353           ->intersect(slsState.geoContext, sls.position(slsState),
0354                       navDir * sls.direction(slsState), BoundaryCheck(false))
0355           .closest(),
0356       navDir, true);
0357   CHECK_CLOSE_ABS(slsState.stepSize.value(), 2, 1e-6);
0358
0359   // Test the bound state construction
0360   FreeToBoundCorrection freeToBoundCorrection(false);
0361   auto boundState =
0362       sls.boundState(slsState, *plane, true, freeToBoundCorrection).value();
0363   auto boundPars = std::get<0>(boundState);
0364   CHECK_CLOSE_ABS(boundPars.position(tgContext), bp.position(tgContext), 1e-6);
0365   CHECK_CLOSE_ABS(boundPars.momentum(), bp.momentum(), 1e-6);
0366   CHECK_CLOSE_ABS(boundPars.charge(), bp.charge(), 1e-6);
0367   CHECK_CLOSE_ABS(boundPars.time(), bp.time(), 1e-6);
0368   BOOST_CHECK(boundPars.covariance().has_value());
0369   BOOST_CHECK_NE(*boundPars.covariance(), cov);
0370   CHECK_CLOSE_COVARIANCE(std::get<1>(boundState),
0371                          BoundMatrix(BoundMatrix::Identity()), 1e-6);
0372   CHECK_CLOSE_ABS(std::get<2>(boundState), 0., 1e-6);
0373
0374   // Transport the covariance in the context of a surface
0375   sls.transportCovarianceToBound(slsState, *plane, freeToBoundCorrection);
0376   BOOST_CHECK_NE(slsState.cov, cov);
0377   BOOST_CHECK_NE(slsState.jacToGlobal, BoundToFreeMatrix::Zero());
0378   BOOST_CHECK_EQUAL(slsState.jacTransport, FreeMatrix::Identity());
0379   BOOST_CHECK_EQUAL(slsState.derivative, FreeVector::Zero());
0380
0381   // Update in context of a surface
0382   freeParams = transformBoundToFreeParameters(bp.referenceSurface(), tgContext,
0383                                               bp.parameters());
0384
0385   BOOST_CHECK(bp.covariance().has_value());
0386   sls.update(slsState, freeParams, bp.parameters(), 2 * (*bp.covariance()),
0387              *plane);
0388   CHECK_CLOSE_OR_SMALL(sls.position(slsState), pos, eps, eps);
0389   BOOST_CHECK_EQUAL(sls.charge(slsState), charge);
0390   CHECK_CLOSE_OR_SMALL(sls.time(slsState), time, eps, eps);
0391   CHECK_CLOSE_COVARIANCE(slsState.cov, Covariance(2. * cov), 1e-6);
0392 }
0393
0394 }  // namespace Acts::Test