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