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