Core/Propagator/LoopProtectionTests.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/data/test_case.hpp>
0010 #include <boost/test/unit_test.hpp>
0011
0012 #include "Acts/Definitions/Algebra.hpp"
0013 #include "Acts/Definitions/Direction.hpp"
0014 #include "Acts/Definitions/Units.hpp"
0015 #include "Acts/EventData/GenericCurvilinearTrackParameters.hpp"
0016 #include "Acts/EventData/ParticleHypothesis.hpp"
0017 #include "Acts/EventData/TrackParameters.hpp"
0018 #include "Acts/Geometry/GeometryContext.hpp"
0019 #include "Acts/MagneticField/ConstantBField.hpp"
0020 #include "Acts/MagneticField/MagneticFieldContext.hpp"
0021 #include "Acts/Propagator/ActorList.hpp"
0022 #include "Acts/Propagator/EigenStepper.hpp"
0023 #include "Acts/Propagator/Propagator.hpp"
0024 #include "Acts/Propagator/StandardAborters.hpp"
0025 #include "Acts/Propagator/detail/LoopProtection.hpp"
0026 #include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
0027 #include "Acts/Utilities/Logger.hpp"
0028 #include "Acts/Utilities/Result.hpp"
0029
0030 #include <algorithm>
0031 #include <array>
0032 #include <cmath>
0033 #include <limits>
0034 #include <memory>
0035 #include <numbers>
0036 #include <optional>
0037 #include <random>
0038 #include <string>
0039 #include <tuple>
0040 #include <utility>
0041
0042 namespace bdata = boost::unit_test::data;
0043 using namespace Acts::UnitLiterals;
0044 using namespace Acts::detail;
0045
0046 namespace Acts::Test {
0047
0048 // Create a test context
0049 GeometryContext tgContext = GeometryContext();
0050 MagneticFieldContext mfContext = MagneticFieldContext();
0051
0052 /// @brief mockup of stepping state
0053 struct SteppingState {
0054   /// Parameters
0055   Vector3 pos = Vector3(0., 0., 0.);
0056   Vector3 dir = Vector3(0., 0., 1);
0057   double p = 100_MeV;
0058 };
0059
0060 /// @brief mockup of stepping state
0061 struct Stepper {
0062   Vector3 field = Vector3(0., 0., 2_T);
0063
0064   /// Get the field for the stepping, it checks first if the access is still
0065   /// within the Cell, and updates the cell if necessary.
0066   ///
0067   /// @param [in,out] state is the propagation state associated with the track
0068   ///                 the magnetic field cell is used (and potentially
0069   ///                 updated)
0070   /// @param [in] pos is the field position
0071   Result<Vector3> getField(SteppingState& /*state*/,
0072                            const Vector3& /*pos*/) const {
0073     // get the field from the cell
0074     return Result<Vector3>::success(field);
0075   }
0076
0077   /// Access method - position
0078   Vector3 position(const SteppingState& state) const { return state.pos; }
0079
0080   /// Access method - direction
0081   Vector3 direction(const SteppingState& state) const { return state.dir; }
0082
0083   /// Access method - momentum
0084   double absoluteMomentum(const SteppingState& state) const { return state.p; }
0085 };
0086
0087 /// @brief mockup of navigation state
0088 struct NavigationState {
0089   bool navigationBreak = false;
0090 };
0091
0092 /// @brief mockup of the Propagator Options
0093 struct Options {
0094   /// Absolute maximum path length
0095   double pathLimit = std::numeric_limits<double>::max();
0096   bool loopProtection = true;
0097   double loopFraction = 0.5;
0098   Direction direction = Direction::Forward();
0099
0100   bool debug = false;
0101   std::string debugString;
0102   int debugMsgWidth = 60;
0103   int debugPfxWidth = 30;
0104
0105   /// Contains: target aborters
0106   ActorList<PathLimitReached> abortList;
0107
0108   const Acts::Logger& logger = Acts::getDummyLogger();
0109 };
0110
0111 /// @brief mockup of propagtor state
0112 struct PropagatorState {
0113   /// Contains: stepping state
0114   SteppingState stepping;
0115   /// Contains: navigation state
0116   NavigationState navigation;
0117   /// Contains: options
0118   Options options;
0119 };
0120
0121 // This test case checks that no segmentation fault appears
0122 // - this tests the collection of surfaces
0123 BOOST_DATA_TEST_CASE(
0124     loop_aborter_test,
0125     bdata::random((bdata::engine = std::mt19937(), bdata::seed = 21,
0126                    bdata::distribution = std::uniform_real_distribution<double>(
0127                        -std::numbers::pi, std::numbers::pi))) ^
0128         bdata::random(
0129             (bdata::engine = std::mt19937(), bdata::seed = 22,
0130              bdata::distribution = std::uniform_real_distribution<double>(
0131                  -std::numbers::pi, std::numbers::pi))) ^
0132         bdata::xrange(1),
0133     phi, deltaPhi, index) {
0134   (void)index;
0135   (void)deltaPhi;
0136
0137   PropagatorState pState;
0138   pState.stepping.dir = Vector3(cos(phi), sin(phi), 0.);
0139   pState.stepping.p = 100_MeV;
0140
0141   Stepper pStepper;
0142
0143   auto& pathLimit = pState.options.abortList.get<PathLimitReached>();
0144   auto initialLimit = pathLimit.internalLimit;
0145
0146   detail::setupLoopProtection(
0147       pState, pStepper, pathLimit, false,
0148       *Acts::getDefaultLogger("LoopProt", Logging::INFO));
0149
0150   auto updatedLimit =
0151       pState.options.abortList.get<PathLimitReached>().internalLimit;
0152   BOOST_CHECK_LT(updatedLimit, initialLimit);
0153 }
0154
0155 using BField = ConstantBField;
0156 using EigenStepper = Acts::EigenStepper<>;
0157 using EigenPropagator = Propagator<EigenStepper>;
0158
0159 const int ntests = 100;
0160 const int skip = 0;
0161
0162 // This test case checks that the propagator with loop LoopProtection
0163 // stops where expected
0164 BOOST_DATA_TEST_CASE(
0165     propagator_loop_protection_test,
0166     bdata::random((bdata::engine = std::mt19937(), bdata::seed = 20,
0167                    bdata::distribution = std::uniform_real_distribution<double>(
0168                        0.5_GeV, 10_GeV))) ^
0169         bdata::random(
0170             (bdata::engine = std::mt19937(), bdata::seed = 21,
0171              bdata::distribution = std::uniform_real_distribution<double>(
0172                  -std::numbers::pi, std::numbers::pi))) ^
0173         bdata::random(
0174             (bdata::engine = std::mt19937(), bdata::seed = 22,
0175              bdata::distribution = std::uniform_real_distribution<double>(
0176                  1., std::numbers::pi - 1.))) ^
0177         bdata::random((bdata::engine = std::mt19937(), bdata::seed = 23,
0178                        bdata::distribution =
0179                            std::uniform_int_distribution<std::uint8_t>(0, 1))) ^
0180         bdata::xrange(ntests),
0181     pT, phi, theta, charge, index) {
0182   if (index < skip) {
0183     return;
0184   }
0185
0186   double px = pT * cos(phi);
0187   double py = pT * sin(phi);
0188   double pz = pT / tan(theta);
0189   double p = pT / sin(theta);
0190   double q = -1 + 2 * charge;
0191
0192   const double Bz = 2_T;
0193   auto bField = std::make_shared<BField>(Vector3{0, 0, Bz});
0194   EigenStepper estepper(bField);
0195   EigenPropagator epropagator(std::move(estepper));
0196
0197   // define start parameters
0198   CurvilinearTrackParameters start(Vector4(0, 0, 0, 42), phi, theta, q / p,
0199                                    std::nullopt, ParticleHypothesis::pion());
0200
0201   using PropagatorOptions = EigenPropagator::Options<ActorList<>>;
0202   PropagatorOptions options(tgContext, mfContext);
0203   options.maxSteps = 1e6;
0204   const auto& result = epropagator.propagate(start, options).value();
0205
0206   // this test assumes state.options.loopFraction = 0.5
0207   CHECK_CLOSE_REL(px, -result.endParameters->momentum().x(), 1e-2);
0208   CHECK_CLOSE_REL(py, -result.endParameters->momentum().y(), 1e-2);
0209   CHECK_CLOSE_REL(pz, result.endParameters->momentum().z(), 1e-2);
0210 }
0211
0212 }  // namespace Acts::Test