Core/Propagator/StraightLineStepperTests.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/EventData/GenericBoundTrackParameters.hpp"
0016 #include "Acts/EventData/TrackParameters.hpp"
0017 #include "Acts/EventData/TransformationHelpers.hpp"
0018 #include "Acts/EventData/detail/CorrectedTransformationFreeToBound.hpp"
0019 #include "Acts/Geometry/GeometryContext.hpp"
0020 #include "Acts/MagneticField/MagneticFieldContext.hpp"
0021 #include "Acts/Propagator/ConstrainedStep.hpp"
0022 #include "Acts/Propagator/StraightLineStepper.hpp"
0023 #include "Acts/Surfaces/BoundaryTolerance.hpp"
0024 #include "Acts/Surfaces/CurvilinearSurface.hpp"
0025 #include "Acts/Surfaces/PlaneSurface.hpp"
0026 #include "Acts/Utilities/Result.hpp"
0027 #include "ActsTests/CommonHelpers/FloatComparisons.hpp"
0028
0029 #include <limits>
0030 #include <memory>
0031 #include <optional>
0032 #include <string>
0033
0034 using namespace Acts;
0035 using namespace Acts::UnitLiterals;
0036 using Acts::VectorHelpers::makeVector4;
0037
0038 namespace ActsTests {
0039
0040 using Covariance = BoundSquareMatrix;
0041
0042 static constexpr auto eps = 2 * std::numeric_limits<double>::epsilon();
0043
0044 BOOST_AUTO_TEST_SUITE(PropagatorSuite)
0045
0046 /// These tests are aiming to test whether the state setup is working properly
0047 BOOST_AUTO_TEST_CASE(straight_line_stepper_state_test) {
0048   // Set up some variables
0049   GeometryContext tgContext = GeometryContext();
0050   MagneticFieldContext mfContext = MagneticFieldContext();
0051   double stepSize = 123.;
0052
0053   Vector3 pos(1., 2., 3.);
0054   Vector3 dir(4., 5., 6.);
0055   double time = 7.;
0056   double absMom = 8.;
0057   double charge = -1.;
0058
0059   StraightLineStepper::Options slsOptions(tgContext, mfContext);
0060   slsOptions.maxStepSize = stepSize;
0061
0062   // Test charged parameters without covariance matrix
0063   BoundTrackParameters cp = BoundTrackParameters::createCurvilinear(
0064       makeVector4(pos, time), dir, charge / absMom, std::nullopt,
0065       ParticleHypothesis::pion());
0066
0067   StraightLineStepper sls;
0068   StraightLineStepper::State slsState = sls.makeState(slsOptions);
0069   sls.initialize(slsState, cp);
0070
0071   // Test the result & compare with the input/test for reasonable members
0072   BOOST_CHECK_EQUAL(slsState.jacToGlobal, BoundToFreeMatrix::Zero());
0073   BOOST_CHECK_EQUAL(slsState.jacTransport, FreeMatrix::Identity());
0074   BOOST_CHECK_EQUAL(slsState.derivative, FreeVector::Zero());
0075   BOOST_CHECK(!slsState.covTransport);
0076   BOOST_CHECK_EQUAL(slsState.cov, Covariance::Zero());
0077   CHECK_CLOSE_OR_SMALL(sls.position(slsState), pos, eps, eps);
0078   CHECK_CLOSE_OR_SMALL(sls.direction(slsState), dir.normalized(), eps, eps);
0079   CHECK_CLOSE_REL(sls.absoluteMomentum(slsState), absMom, eps);
0080   BOOST_CHECK_EQUAL(sls.charge(slsState), charge);
0081   CHECK_CLOSE_OR_SMALL(sls.time(slsState), time, eps, eps);
0082   BOOST_CHECK_EQUAL(slsState.pathAccumulated, 0.);
0083   BOOST_CHECK_EQUAL(slsState.stepSize.value(), stepSize);
0084   BOOST_CHECK_EQUAL(slsState.previousStepSize, 0.);
0085
0086   // Test without charge and covariance matrix
0087   BoundTrackParameters ncp = BoundTrackParameters::createCurvilinear(
0088       makeVector4(pos, time), dir, 1 / absMom, std::nullopt,
0089       ParticleHypothesis::pion0());
0090   sls.initialize(slsState, ncp);
0091
0092   // Test with covariance matrix
0093   Covariance cov = 8. * Covariance::Identity();
0094   ncp = BoundTrackParameters::createCurvilinear(makeVector4(pos, time), dir,
0095                                                 1 / absMom, cov,
0096                                                 ParticleHypothesis::pion0());
0097   sls.initialize(slsState, ncp);
0098   BOOST_CHECK_NE(slsState.jacToGlobal, BoundToFreeMatrix::Zero());
0099   BOOST_CHECK(slsState.covTransport);
0100   BOOST_CHECK_EQUAL(slsState.cov, cov);
0101 }
0102
0103 /// These tests are aiming to test the functions of the StraightLineStepper
0104 /// The numerical correctness of the stepper is tested in the integration tests
0105 BOOST_AUTO_TEST_CASE(straight_line_stepper_test) {
0106   // Set up some variables for the state
0107   GeometryContext tgContext = GeometryContext();
0108   MagneticFieldContext mfContext = MagneticFieldContext();
0109   Direction navDir = Direction::Backward();
0110   double stepSize = 123.;
0111
0112   // Construct the parameters
0113   Vector3 pos(1., 2., 3.);
0114   Vector3 dir = Vector3(4., 5., 6.).normalized();
0115   double time = 7.;
0116   double absMom = 8.;
0117   double charge = -1.;
0118   Covariance cov = 8. * Covariance::Identity();
0119   BoundTrackParameters cp = BoundTrackParameters::createCurvilinear(
0120       makeVector4(pos, time), dir, charge / absMom, cov,
0121       ParticleHypothesis::pion());
0122
0123   StraightLineStepper::Options options(tgContext, mfContext);
0124   options.maxStepSize = stepSize;
0125   options.initialStepSize = 10_m;
0126
0127   // Build the stepper and the state
0128   StraightLineStepper sls;
0129   StraightLineStepper::State slsState = sls.makeState(options);
0130   sls.initialize(slsState, cp);
0131
0132   // Test the getters
0133   CHECK_CLOSE_ABS(sls.position(slsState), pos, 1e-6);
0134   CHECK_CLOSE_ABS(sls.direction(slsState), dir, 1e-6);
0135   CHECK_CLOSE_ABS(sls.absoluteMomentum(slsState), absMom, 1e-6);
0136   BOOST_CHECK_EQUAL(sls.charge(slsState), charge);
0137   BOOST_CHECK_EQUAL(sls.time(slsState), time);
0138
0139   // Step size modifies
0140   const std::string originalStepSize = slsState.stepSize.toString();
0141
0142   sls.updateStepSize(slsState, -1337., ConstrainedStep::Type::Navigator);
0143   BOOST_CHECK_EQUAL(slsState.previousStepSize, stepSize);
0144   BOOST_CHECK_EQUAL(slsState.stepSize.value(), -1337.);
0145
0146   sls.releaseStepSize(slsState, ConstrainedStep::Type::Navigator);
0147   BOOST_CHECK_EQUAL(slsState.stepSize.value(), stepSize);
0148   BOOST_CHECK_EQUAL(sls.outputStepSize(slsState), originalStepSize);
0149
0150   // Test the curvilinear state construction
0151   auto curvState = sls.curvilinearState(slsState);
0152   auto curvPars = std::get<0>(curvState);
0153   CHECK_CLOSE_ABS(curvPars.position(tgContext), cp.position(tgContext), 1e-6);
0154   CHECK_CLOSE_ABS(curvPars.absoluteMomentum(), cp.absoluteMomentum(), 1e-6);
0155   CHECK_CLOSE_ABS(curvPars.charge(), cp.charge(), 1e-6);
0156   CHECK_CLOSE_ABS(curvPars.time(), cp.time(), 1e-6);
0157   BOOST_CHECK(curvPars.covariance().has_value());
0158   BOOST_CHECK_NE(*curvPars.covariance(), cov);
0159   CHECK_CLOSE_COVARIANCE(std::get<1>(curvState),
0160                          BoundMatrix(BoundMatrix::Identity()), 1e-6);
0161   CHECK_CLOSE_ABS(std::get<2>(curvState), 0., 1e-6);
0162
0163   // Test the update method
0164   Vector3 newPos(2., 4., 8.);
0165   Vector3 newMom(3., 9., 27.);
0166   double newTime(321.);
0167   sls.update(slsState, newPos, newMom.normalized(), charge / newMom.norm(),
0168              newTime);
0169   CHECK_CLOSE_ABS(sls.position(slsState), newPos, 1e-6);
0170   CHECK_CLOSE_ABS(sls.direction(slsState), newMom.normalized(), 1e-6);
0171   CHECK_CLOSE_ABS(sls.absoluteMomentum(slsState), newMom.norm(), 1e-6);
0172   BOOST_CHECK_EQUAL(sls.charge(slsState), charge);
0173   BOOST_CHECK_EQUAL(sls.time(slsState), newTime);
0174
0175   // The covariance transport
0176   slsState.cov = cov;
0177   sls.transportCovarianceToCurvilinear(slsState);
0178   BOOST_CHECK_NE(slsState.cov, cov);
0179   BOOST_CHECK_NE(slsState.jacToGlobal, BoundToFreeMatrix::Zero());
0180   BOOST_CHECK_EQUAL(slsState.jacTransport, FreeMatrix::Identity());
0181   BOOST_CHECK_EQUAL(slsState.derivative, FreeVector::Zero());
0182
0183   // Perform a step without and with covariance transport
0184   slsState.cov = cov;
0185
0186   slsState.covTransport = false;
0187   double h = sls.step(slsState, navDir, nullptr).value();
0188   BOOST_CHECK_EQUAL(slsState.stepSize.value(), stepSize);
0189   BOOST_CHECK_EQUAL(slsState.stepSize.value(), h * navDir);
0190   CHECK_CLOSE_COVARIANCE(slsState.cov, cov, 1e-6);
0191   BOOST_CHECK_GT(sls.position(slsState).norm(), newPos.norm());
0192   CHECK_CLOSE_ABS(sls.direction(slsState), newMom.normalized(), 1e-6);
0193   CHECK_CLOSE_ABS(sls.absoluteMomentum(slsState), newMom.norm(), 1e-6);
0194   CHECK_CLOSE_ABS(sls.charge(slsState), charge, 1e-6);
0195   BOOST_CHECK_LT(sls.time(slsState), newTime);
0196   BOOST_CHECK_EQUAL(slsState.derivative, FreeVector::Zero());
0197   BOOST_CHECK_EQUAL(slsState.jacTransport, FreeMatrix::Identity());
0198
0199   slsState.covTransport = true;
0200   double h2 = sls.step(slsState, navDir, nullptr).value();
0201   BOOST_CHECK_EQUAL(slsState.stepSize.value(), stepSize);
0202   BOOST_CHECK_EQUAL(h2, h);
0203   CHECK_CLOSE_COVARIANCE(slsState.cov, cov, 1e-6);
0204   BOOST_CHECK_GT(sls.position(slsState).norm(), newPos.norm());
0205   CHECK_CLOSE_ABS(sls.direction(slsState), newMom.normalized(), 1e-6);
0206   CHECK_CLOSE_ABS(sls.absoluteMomentum(slsState), newMom.norm(), 1e-6);
0207   CHECK_CLOSE_ABS(sls.charge(slsState), charge, 1e-6);
0208   BOOST_CHECK_LT(sls.time(slsState), newTime);
0209   BOOST_CHECK_NE(slsState.derivative, FreeVector::Zero());
0210   BOOST_CHECK_NE(slsState.jacTransport, FreeMatrix::Identity());
0211
0212   /// Test the state reset
0213   // Construct the parameters
0214   Vector3 pos2(1.5, -2.5, 3.5);
0215   Vector3 dir2 = Vector3(4.5, -5.5, 6.5).normalized();
0216   double time2 = 7.5;
0217   double absMom2 = 8.5;
0218   double charge2 = 1.;
0219   BoundSquareMatrix cov2 = 8.5 * Covariance::Identity();
0220   BoundTrackParameters cp2 = BoundTrackParameters::createCurvilinear(
0221       makeVector4(pos2, time2), dir2, charge2 / absMom2, cov2,
0222       ParticleHypothesis::pion());
0223   BOOST_CHECK(cp2.covariance().has_value());
0224   FreeVector freeParams = transformBoundToFreeParameters(
0225       cp2.referenceSurface(), tgContext, cp2.parameters());
0226   navDir = Direction::Forward();
0227
0228   // Reset all possible parameters
0229   StraightLineStepper::State slsStateCopy = slsState;
0230   sls.initialize(slsStateCopy, cp2.parameters(), cp2.covariance(),
0231                  cp2.particleHypothesis(), cp2.referenceSurface());
0232   // Test all components
0233   BOOST_CHECK_NE(slsStateCopy.jacToGlobal, BoundToFreeMatrix::Zero());
0234   BOOST_CHECK_NE(slsStateCopy.jacToGlobal, slsState.jacToGlobal);
0235   BOOST_CHECK_EQUAL(slsStateCopy.jacTransport, FreeMatrix::Identity());
0236   BOOST_CHECK_EQUAL(slsStateCopy.derivative, FreeVector::Zero());
0237   BOOST_CHECK(slsStateCopy.covTransport);
0238   BOOST_CHECK_EQUAL(slsStateCopy.cov, cov2);
0239   CHECK_CLOSE_ABS(sls.position(slsStateCopy),
0240                   freeParams.template segment<3>(eFreePos0), 1e-6);
0241   CHECK_CLOSE_ABS(sls.direction(slsStateCopy),
0242                   freeParams.template segment<3>(eFreeDir0).normalized(), 1e-6);
0243   CHECK_CLOSE_ABS(sls.absoluteMomentum(slsStateCopy),
0244                   std::abs(1. / freeParams[eFreeQOverP]), 1e-6);
0245   CHECK_CLOSE_ABS(sls.charge(slsStateCopy), -sls.charge(slsState), 1e-6);
0246   CHECK_CLOSE_ABS(sls.time(slsStateCopy), freeParams[eFreeTime], 1e-6);
0247   BOOST_CHECK_EQUAL(slsStateCopy.pathAccumulated, 0.);
0248   BOOST_CHECK_EQUAL(slsStateCopy.stepSize.value(), stepSize);
0249   BOOST_CHECK_EQUAL(slsStateCopy.previousStepSize, 0.);
0250
0251   /// Repeat with surface related methods
0252   std::shared_ptr<PlaneSurface> plane =
0253       CurvilinearSurface(pos, dir).planeSurface();
0254   auto bp = BoundTrackParameters::create(
0255                 tgContext, plane, makeVector4(pos, time), dir, charge / absMom,
0256                 cov, ParticleHypothesis::pion())
0257                 .value();
0258   slsState = sls.makeState(options);
0259   sls.initialize(slsState, bp);
0260
0261   // Test the intersection in the context of a surface
0262   std::shared_ptr<PlaneSurface> targetSurface =
0263       CurvilinearSurface(pos + navDir * 2. * dir, dir).planeSurface();
0264   sls.updateSurfaceStatus(slsState, *targetSurface, 0, navDir,
0265                           BoundaryTolerance::Infinite(), s_onSurfaceTolerance,
0266                           ConstrainedStep::Type::Navigator);
0267   CHECK_CLOSE_ABS(slsState.stepSize.value(ConstrainedStep::Type::Navigator),
0268                   navDir * 2., 1e-6);
0269
0270   const auto getNavigationTarget = [&](const Surface& s,
0271                                        const BoundaryTolerance& bt) {
0272     auto [intersection, intersectionIndex] =
0273         s.intersect(tgContext, sls.position(slsState),
0274                     navDir * sls.direction(slsState), bt)
0275             .closestWithIndex();
0276     return NavigationTarget(intersection, intersectionIndex, s, bt);
0277   };
0278
0279   // Test the step size modification in the context of a surface
0280   sls.updateStepSize(
0281       slsState,
0282       getNavigationTarget(*targetSurface, BoundaryTolerance::Infinite()),
0283       navDir, ConstrainedStep::Type::Navigator);
0284   CHECK_CLOSE_ABS(slsState.stepSize.value(), 2, 1e-6);
0285   slsState.stepSize.setUser(navDir * stepSize);
0286   sls.releaseStepSize(slsState, ConstrainedStep::Type::Navigator);
0287   sls.updateStepSize(
0288       slsState,
0289       getNavigationTarget(*targetSurface, BoundaryTolerance::Infinite()),
0290       navDir, ConstrainedStep::Type::Navigator);
0291   CHECK_CLOSE_ABS(slsState.stepSize.value(), 2, 1e-6);
0292
0293   // Test the bound state construction
0294   FreeToBoundCorrection freeToBoundCorrection(false);
0295   auto boundState =
0296       sls.boundState(slsState, *plane, true, freeToBoundCorrection).value();
0297   auto boundPars = std::get<0>(boundState);
0298   CHECK_CLOSE_ABS(boundPars.position(tgContext), bp.position(tgContext), 1e-6);
0299   CHECK_CLOSE_ABS(boundPars.momentum(), bp.momentum(), 1e-6);
0300   CHECK_CLOSE_ABS(boundPars.charge(), bp.charge(), 1e-6);
0301   CHECK_CLOSE_ABS(boundPars.time(), bp.time(), 1e-6);
0302   BOOST_CHECK(boundPars.covariance().has_value());
0303   BOOST_CHECK_NE(*boundPars.covariance(), cov);
0304   CHECK_CLOSE_COVARIANCE(std::get<1>(boundState),
0305                          BoundMatrix(BoundMatrix::Identity()), 1e-6);
0306   CHECK_CLOSE_ABS(std::get<2>(boundState), 0., 1e-6);
0307
0308   // Transport the covariance in the context of a surface
0309   sls.transportCovarianceToBound(slsState, *plane, freeToBoundCorrection);
0310   BOOST_CHECK_NE(slsState.cov, cov);
0311   BOOST_CHECK_NE(slsState.jacToGlobal, BoundToFreeMatrix::Zero());
0312   BOOST_CHECK_EQUAL(slsState.jacTransport, FreeMatrix::Identity());
0313   BOOST_CHECK_EQUAL(slsState.derivative, FreeVector::Zero());
0314
0315   // Update in context of a surface
0316   freeParams = transformBoundToFreeParameters(bp.referenceSurface(), tgContext,
0317                                               bp.parameters());
0318
0319   BOOST_CHECK(bp.covariance().has_value());
0320   sls.update(slsState, freeParams, bp.parameters(), 2 * (*bp.covariance()),
0321              *plane);
0322   CHECK_CLOSE_OR_SMALL(sls.position(slsState), pos, eps, eps);
0323   BOOST_CHECK_EQUAL(sls.charge(slsState), charge);
0324   CHECK_CLOSE_OR_SMALL(sls.time(slsState), time, eps, eps);
0325   CHECK_CLOSE_COVARIANCE(slsState.cov, Covariance(2. * cov), 1e-6);
0326 }
0327
0328 BOOST_AUTO_TEST_SUITE_END()
0329
0330 }  // namespace ActsTests