Core/Propagator/LoopProtectionTests.cpp

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