Core/Propagator/HelixStepperTests.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/TrackParametrization.hpp"
0014 #include "Acts/Definitions/Units.hpp"
0015 #include "Acts/EventData/BoundTrackParameters.hpp"
0016 #include "Acts/Geometry/GeometryContext.hpp"
0017 #include "Acts/MagneticField/ConstantBField.hpp"
0018 #include "Acts/MagneticField/MagneticFieldContext.hpp"
0019 #include "Acts/MagneticField/MagneticFieldProvider.hpp"
0020 #include "Acts/Propagator/HelixStepper.hpp"
0021 #include "Acts/Propagator/StraightLineStepper.hpp"
0022 #include "ActsTests/CommonHelpers/FloatComparisons.hpp"
0023
0024 #include <array>
0025 #include <cmath>
0026 #include <memory>
0027 #include <numbers>
0028 #include <optional>
0029
0030 using namespace Acts;
0031 using namespace Acts::UnitLiterals;
0032 using Acts::VectorHelpers::makeVector4;
0033
0034 namespace ActsTests {
0035
0036 namespace {
0037
0038 const auto tgContext = GeometryContext::dangerouslyDefaultConstruct();
0039 const MagneticFieldContext mfContext;
0040
0041 /// Field with a constant gradient of Bz along x, to force step rejections.
0042 class GradientField final : public MagneticFieldProvider {
0043  public:
0044   struct Cache {
0045     explicit Cache(const MagneticFieldContext& /*mctx*/) {}
0046   };
0047
0048   GradientField(double bz, double dBzdx) : m_bz(bz), m_dBzdx(dBzdx) {}
0049
0050   MagneticFieldProvider::Cache makeCache(
0051       const MagneticFieldContext& mctx) const override {
0052     return MagneticFieldProvider::Cache(std::in_place_type<Cache>, mctx);
0053   }
0054
0055   Result<Vector3> getField(
0056       const Vector3& position,
0057       MagneticFieldProvider::Cache& /*cache*/) const override {
0058     return Result<Vector3>::success(
0059         Vector3(0., 0., m_bz + m_dBzdx * position.x()));
0060   }
0061
0062  private:
0063   double m_bz;
0064   double m_dBzdx;
0065 };
0066
0067 BoundTrackParameters makeParameters(const Vector3& pos, const Vector3& dir,
0068                                     double qop,
0069                                     std::optional<BoundMatrix> cov) {
0070   return BoundTrackParameters::createCurvilinear(
0071       makeVector4(pos, 0.), dir, qop, cov, ParticleHypothesis::pion());
0072 }
0073
0074 /// Run one step of fixed size from the free parameters and return them.
0075 FreeVector stepFrom(const HelixStepper& stepper, const FreeVector& start,
0076                     double h, FreeMatrix* jacTransport = nullptr) {
0077   HelixStepper::Options options(tgContext, mfContext);
0078   options.initialStepSize = std::abs(h);
0079   options.maxStepSize = std::abs(h);
0080   HelixStepper::State state = stepper.makeState(options);
0081   stepper.initialize(state, makeParameters(Vector3::Zero(), Vector3::UnitX(),
0082                                            1., BoundMatrix::Identity()));
0083   state.pars = start;
0084   state.field.reset();
0085   auto res = stepper.step(state, Direction::fromScalar(h), nullptr);
0086   BOOST_REQUIRE(res.ok());
0087   CHECK_CLOSE_REL(*res, h, 1e-15);
0088   if (jacTransport != nullptr) {
0089     *jacTransport = state.jacTransport;
0090   }
0091   return state.pars;
0092 }
0093
0094 }  // namespace
0095
0096 BOOST_AUTO_TEST_SUITE(PropagatorSuite)
0097
0098 BOOST_AUTO_TEST_CASE(helix_stepper_state_test) {
0099   auto bField = std::make_shared<ConstantBField>(Vector3(0., 0., 2_T));
0100   HelixStepper stepper(bField);
0101
0102   HelixStepper::Options options(tgContext, mfContext);
0103   options.maxStepSize = 123.;
0104   HelixStepper::State state = stepper.makeState(options);
0105
0106   Vector3 pos(1., 2., 3.);
0107   Vector3 dir = Vector3(4., 5., 6.).normalized();
0108   stepper.initialize(state,
0109                      makeParameters(pos, dir, -1. / 8_GeV, std::nullopt));
0110
0111   BOOST_CHECK(!stepper.hasCovariance(state));
0112   BOOST_CHECK(!state.field.has_value());
0113   CHECK_CLOSE_ABS(stepper.position(state), pos, 1e-12);
0114   CHECK_CLOSE_ABS(stepper.direction(state), dir, 1e-12);
0115   CHECK_CLOSE_REL(stepper.absoluteMomentum(state), 8_GeV, 1e-12);
0116   BOOST_CHECK_EQUAL(stepper.charge(state), -1.);
0117   BOOST_CHECK_EQUAL(state.stepSize.value(), 123.);
0118
0119   stepper.initialize(
0120       state, makeParameters(pos, dir, -1. / 8_GeV, BoundMatrix::Identity()));
0121   BOOST_CHECK(stepper.hasCovariance(state));
0122   BOOST_CHECK_NE(state.jacToGlobal, BoundToFreeMatrix::Zero());
0123   BOOST_CHECK_EQUAL(state.jacTransport, FreeMatrix::Identity());
0124 }
0125
0126 /// A full turn of the helix returns to the start, shifted along the field.
0127 BOOST_AUTO_TEST_CASE(helix_stepper_full_turn) {
0128   const double bz = 2_T;
0129   const double p = 1_GeV;
0130   const double qop = 1. / p;
0131   auto bField = std::make_shared<ConstantBField>(Vector3(0., 0., bz));
0132   HelixStepper stepper(bField);
0133
0134   const Vector3 pos(10., -20., 30.);
0135   const Vector3 dir = Vector3(1., 1., 0.5).normalized();
0136   // The direction turns about the field once per this path length
0137   const double turn = 2. * std::numbers::pi / (std::abs(qop) * bz);
0138
0139   for (int nSteps : {1, 7}) {
0140     HelixStepper::Options options(tgContext, mfContext);
0141     options.initialStepSize = turn / nSteps;
0142     options.maxStepSize = turn / nSteps;
0143     HelixStepper::State state = stepper.makeState(options);
0144     stepper.initialize(state, makeParameters(pos, dir, qop, std::nullopt));
0145     for (int i = 0; i < nSteps; ++i) {
0146       BOOST_REQUIRE(stepper.step(state, Direction::Forward(), nullptr).ok());
0147     }
0148     CHECK_CLOSE_ABS(stepper.position(state),
0149                     Vector3(pos + Vector3(0., 0., turn * dir.z())), 1e-9);
0150     CHECK_CLOSE_ABS(stepper.direction(state), dir, 1e-12);
0151     BOOST_CHECK_EQUAL(state.statistics.nRejectedSteps, 0u);
0152   }
0153 }
0154
0155 /// The helix with a vanishing field is a straight line.
0156 BOOST_AUTO_TEST_CASE(helix_stepper_straight_limit) {
0157   auto bField = std::make_shared<ConstantBField>(Vector3::Zero());
0158   HelixStepper helix(bField);
0159   StraightLineStepper line;
0160
0161   const double h = 123.;
0162   HelixStepper::Options helixOptions(tgContext, mfContext);
0163   helixOptions.maxStepSize = h;
0164   StraightLineStepper::Options lineOptions(tgContext, mfContext);
0165   lineOptions.maxStepSize = h;
0166
0167   auto params =
0168       makeParameters(Vector3(1., 2., 3.), Vector3(4., 5., 6.).normalized(),
0169                      -1. / 2_GeV, BoundMatrix::Identity());
0170   HelixStepper::State helixState = helix.makeState(helixOptions);
0171   helix.initialize(helixState, params);
0172   StraightLineStepper::State lineState = line.makeState(lineOptions);
0173   line.initialize(lineState, params);
0174
0175   BOOST_REQUIRE(helix.step(helixState, Direction::Forward(), nullptr).ok());
0176   BOOST_REQUIRE(line.step(lineState, Direction::Forward(), nullptr).ok());
0177
0178   CHECK_CLOSE_ABS(helixState.pars, lineState.pars, 1e-12);
0179   CHECK_CLOSE_ABS(helixState.jacTransport, lineState.jacTransport, 1e-12);
0180   CHECK_CLOSE_ABS(helixState.derivative, lineState.derivative, 1e-12);
0181 }
0182
0183 /// The analytic step jacobian matches central finite differences.
0184 BOOST_AUTO_TEST_CASE(helix_stepper_jacobian) {
0185   auto bField = std::make_shared<ConstantBField>(Vector3(0.1_T, -0.2_T, 2_T));
0186   HelixStepper stepper(bField);
0187
0188   FreeVector start = FreeVector::Zero();
0189   start.segment<3>(eFreePos0) = Vector3(1., 2., 3.);
0190   start[eFreeTime] = 4.;
0191   start.segment<3>(eFreeDir0) = Vector3(4., -5., 6.).normalized();
0192   start[eFreeQOverP] = -1. / 0.7_GeV;
0193
0194   // Small and large turning angles, forward and backward
0195   for (double h : {1e-3, 50., 1500., -400.}) {
0196     FreeMatrix jac;
0197     stepFrom(stepper, start, h, &jac);
0198
0199     std::array<double, eFreeSize> deltas{};
0200     deltas.fill(1e-6);
0201     deltas[eFreeQOverP] = 1e-6 * std::abs(start[eFreeQOverP]);
0202
0203     FreeMatrix numeric;
0204     for (std::size_t j = 0; j < eFreeSize; ++j) {
0205       FreeVector plus = start;
0206       FreeVector minus = start;
0207       plus[j] += deltas[j];
0208       minus[j] -= deltas[j];
0209       numeric.col(j) =
0210           (stepFrom(stepper, plus, h) - stepFrom(stepper, minus, h)) /
0211           (2. * deltas[j]);
0212     }
0213
0214     BOOST_TEST_CONTEXT("h = " << h) {
0215       CHECK_CLOSE_OR_SMALL(jac, numeric, 1e-5, 1e-6);
0216     }
0217   }
0218 }
0219
0220 /// A field gradient shortens the step.
0221 BOOST_AUTO_TEST_CASE(helix_stepper_step_size_control) {
0222   const auto params = makeParameters(Vector3::Zero(), Vector3::UnitX(),
0223                                      1. / 1_GeV, std::nullopt);
0224
0225   HelixStepper::Options options(tgContext, mfContext);
0226   options.stepTolerance = 1e-4;
0227   options.initialStepSize = 1_m;
0228   options.maxStepSize = 1_m;
0229
0230   HelixStepper gradient(std::make_shared<GradientField>(2_T, 1_T / 1_m));
0231   HelixStepper::State gradientState = gradient.makeState(options);
0232   gradient.initialize(gradientState, params);
0233   auto gradientStep =
0234       gradient.step(gradientState, Direction::Forward(), nullptr);
0235   BOOST_REQUIRE(gradientStep.ok());
0236   BOOST_CHECK_LT(*gradientStep, 1_m);
0237   BOOST_CHECK_GT(gradientState.statistics.nRejectedSteps, 0u);
0238   // The field at the end of the step is cached for the next step
0239   BOOST_REQUIRE(gradientState.field.has_value());
0240   CHECK_CLOSE_ABS(
0241       *gradientState.field,
0242       Vector3(0., 0., 2_T + 1_T / 1_m * gradient.position(gradientState).x()),
0243       1e-15);
0244 }
0245
0246 BOOST_AUTO_TEST_SUITE_END()
0247
0248 }  // namespace ActsTests