Core/Propagator/ExtrapolatorTests.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/TrackParametrization.hpp"
0014 #include "Acts/Definitions/Units.hpp"
0015 #include "Acts/EventData/GenericCurvilinearTrackParameters.hpp"
0016 #include "Acts/EventData/TrackParameters.hpp"
0017 #include "Acts/Geometry/GeometryContext.hpp"
0018 #include "Acts/MagneticField/ConstantBField.hpp"
0019 #include "Acts/MagneticField/MagneticFieldContext.hpp"
0020 #include "Acts/Propagator/ActorList.hpp"
0021 #include "Acts/Propagator/EigenStepper.hpp"
0022 #include "Acts/Propagator/MaterialInteractor.hpp"
0023 #include "Acts/Propagator/Navigator.hpp"
0024 #include "Acts/Propagator/Propagator.hpp"
0025 #include "Acts/Propagator/SurfaceCollector.hpp"
0026 #include "Acts/Surfaces/Surface.hpp"
0027 #include "Acts/Tests/CommonHelpers/CylindricalTrackingGeometry.hpp"
0028 #include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
0029 #include "Acts/Utilities/Result.hpp"
0030
0031 #include <cmath>
0032 #include <cstdint>
0033 #include <memory>
0034 #include <numbers>
0035 #include <optional>
0036 #include <random>
0037 #include <utility>
0038
0039 namespace bdata = boost::unit_test::data;
0040 using namespace Acts::UnitLiterals;
0041
0042 namespace Acts::Test {
0043
0044 // Create a test context
0045 GeometryContext tgContext = GeometryContext();
0046 MagneticFieldContext mfContext = MagneticFieldContext();
0047
0048 // Global definitions
0049 CylindricalTrackingGeometry cGeometry(tgContext);
0050 auto tGeometry = cGeometry();
0051
0052 // Get the navigator and provide the TrackingGeometry
0053 Navigator navigator({tGeometry});
0054
0055 using BFieldType = ConstantBField;
0056 using EigenStepperType = EigenStepper<>;
0057 using EigenPropagatorType = Propagator<EigenStepperType, Navigator>;
0058 using Covariance = BoundSquareMatrix;
0059
0060 auto bField = std::make_shared<BFieldType>(Vector3{0, 0, 2_T});
0061 EigenStepperType estepper(bField);
0062 EigenPropagatorType epropagator(std::move(estepper), std::move(navigator));
0063
0064 const int ntests = 100;
0065
0066 // A plane selector for the SurfaceCollector
0067 struct PlaneSelector {
0068   /// Call operator
0069   /// @param sf The input surface to be checked
0070   bool operator()(const Surface& sf) const {
0071     return (sf.type() == Surface::Plane);
0072   }
0073 };
0074
0075 // This test case checks that no segmentation fault appears
0076 // - simple extrapolation test
0077 BOOST_DATA_TEST_CASE(
0078     test_extrapolation_,
0079     bdata::random((bdata::engine = std::mt19937(), bdata::seed = 0,
0080                    bdata::distribution = std::uniform_real_distribution<double>(
0081                        0.4_GeV, 10_GeV))) ^
0082         bdata::random(
0083             (bdata::engine = std::mt19937(), bdata::seed = 1,
0084              bdata::distribution = std::uniform_real_distribution<double>(
0085                  -std::numbers::pi, std::numbers::pi))) ^
0086         bdata::random(
0087             (bdata::engine = std::mt19937(), bdata::seed = 2,
0088              bdata::distribution = std::uniform_real_distribution<double>(
0089                  1., std::numbers::pi - 1.))) ^
0090         bdata::random((bdata::engine = std::mt19937(), bdata::seed = 3,
0091                        bdata::distribution =
0092                            std::uniform_int_distribution<std::uint8_t>(0, 1))) ^
0093         bdata::xrange(ntests),
0094     pT, phi, theta, charge, index) {
0095   double p = pT / sin(theta);
0096   double q = -1 + 2 * charge;
0097   (void)index;
0098
0099   // define start parameters
0100   /// a covariance matrix to transport
0101   Covariance cov;
0102   // take some major correlations (off-diagonals)
0103   cov << 10_mm, 0, 0.123, 0, 0.5, 0, 0, 10_mm, 0, 0.162, 0, 0, 0.123, 0, 0.1, 0,
0104       0, 0, 0, 0.162, 0, 0.1, 0, 0, 0.5, 0, 0, 0, 1. / (10_GeV), 0, 0, 0, 0, 0,
0105       0, 0;
0106   CurvilinearTrackParameters start(Vector4(0, 0, 0, 0), phi, theta, q / p, cov,
0107                                    ParticleHypothesis::pion());
0108
0109   EigenPropagatorType::Options<> options(tgContext, mfContext);
0110   options.stepping.maxStepSize = 10_cm;
0111   options.pathLimit = 25_cm;
0112
0113   BOOST_CHECK(
0114       epropagator.propagate(start, options).value().endParameters.has_value());
0115 }
0116
0117 // This test case checks that no segmentation fault appears
0118 // - this tests the collection of surfaces
0119 BOOST_DATA_TEST_CASE(
0120     test_surface_collection_,
0121     bdata::random((bdata::engine = std::mt19937(), bdata::seed = 10,
0122                    bdata::distribution = std::uniform_real_distribution<double>(
0123                        0.4_GeV, 10_GeV))) ^
0124         bdata::random(
0125             (bdata::engine = std::mt19937(), bdata::seed = 11,
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 = 12,
0130              bdata::distribution = std::uniform_real_distribution<double>(
0131                  1., std::numbers::pi - 1.))) ^
0132         bdata::random((bdata::engine = std::mt19937(), bdata::seed = 13,
0133                        bdata::distribution =
0134                            std::uniform_int_distribution<std::uint8_t>(0, 1))) ^
0135         bdata::xrange(ntests),
0136     pT, phi, theta, charge, index) {
0137   double p = pT / sin(theta);
0138   double q = -1 + 2 * charge;
0139   (void)index;
0140
0141   // define start parameters
0142   /// a covariance matrix to transport
0143   Covariance cov;
0144   // take some major correlations (off-diagonals)
0145   cov << 10_mm, 0, 0.123, 0, 0.5, 0, 0, 10_mm, 0, 0.162, 0, 0, 0.123, 0, 0.1, 0,
0146       0, 0, 0, 0.162, 0, 0.1, 0, 0, 0.5, 0, 0, 0, 1. / (10_GeV), 0, 0, 0, 0, 0,
0147       0, 0;
0148   CurvilinearTrackParameters start(Vector4(0, 0, 0, 0), phi, theta, q / p, cov,
0149                                    ParticleHypothesis::pion());
0150
0151   // A PlaneSelector for the SurfaceCollector
0152   using PlaneCollector = SurfaceCollector<PlaneSelector>;
0153
0154   EigenPropagatorType::Options<ActorList<PlaneCollector>> options(tgContext,
0155                                                                   mfContext);
0156
0157   options.stepping.maxStepSize = 10_cm;
0158   options.pathLimit = 25_cm;
0159
0160   const auto& result = epropagator.propagate(start, options).value();
0161   auto collector_result = result.get<PlaneCollector::result_type>();
0162
0163   // step through the surfaces and go step by step
0164   EigenPropagatorType::Options<> optionsEmpty(tgContext, mfContext);
0165
0166   optionsEmpty.stepping.maxStepSize = 25_cm;
0167   // Try propagation from start to each surface
0168   for (const auto& colsf : collector_result.collected) {
0169     const auto& csurface = colsf.surface;
0170     // Avoid going to the same surface
0171     // @todo: decide on strategy and write unit test for this
0172     if (csurface == &start.referenceSurface()) {
0173       continue;
0174     }
0175     // Extrapolate & check
0176     const auto& cresult = epropagator.propagate(start, *csurface, optionsEmpty)
0177                               .value()
0178                               .endParameters;
0179     BOOST_CHECK(cresult.has_value());
0180   }
0181 }
0182
0183 // This test case checks that no segmentation fault appears
0184 // - this tests the collection of surfaces
0185 BOOST_DATA_TEST_CASE(
0186     test_material_interactor_,
0187     bdata::random((bdata::engine = std::mt19937(), bdata::seed = 20,
0188                    bdata::distribution = std::uniform_real_distribution<double>(
0189                        0.4_GeV, 10_GeV))) ^
0190         bdata::random(
0191             (bdata::engine = std::mt19937(), bdata::seed = 21,
0192              bdata::distribution = std::uniform_real_distribution<double>(
0193                  -std::numbers::pi, std::numbers::pi))) ^
0194         bdata::random(
0195             (bdata::engine = std::mt19937(), bdata::seed = 22,
0196              bdata::distribution = std::uniform_real_distribution<double>(
0197                  1., std::numbers::pi - 1.))) ^
0198         bdata::random((bdata::engine = std::mt19937(), bdata::seed = 23,
0199                        bdata::distribution =
0200                            std::uniform_int_distribution<std::uint8_t>(0, 1))) ^
0201         bdata::xrange(ntests),
0202     pT, phi, theta, charge, index) {
0203   double p = pT / sin(theta);
0204   double q = -1 + 2 * charge;
0205   (void)index;
0206
0207   // define start parameters
0208   /// a covariance matrix to transport
0209   Covariance cov;
0210   // take some major correlations (off-diagonals)
0211   cov << 10_mm, 0, 0.123, 0, 0.5, 0, 0, 10_mm, 0, 0.162, 0, 0, 0.123, 0, 0.1, 0,
0212       0, 0, 0, 0.162, 0, 0.1, 0, 0, 0.5, 0, 0, 0, 1. / (10_GeV), 0, 0, 0, 0, 0,
0213       0, 0;
0214   CurvilinearTrackParameters start(Vector4(0, 0, 0, 0), phi, theta, q / p, cov,
0215                                    ParticleHypothesis::pion());
0216
0217   EigenPropagatorType::Options<ActorList<MaterialInteractor>> options(
0218       tgContext, mfContext);
0219   options.stepping.maxStepSize = 25_cm;
0220   options.pathLimit = 25_cm;
0221
0222   const auto& result = epropagator.propagate(start, options).value();
0223   if (result.endParameters) {
0224     // test that you actually lost some energy
0225     BOOST_CHECK_LT(result.endParameters->absoluteMomentum(),
0226                    start.absoluteMomentum());
0227   }
0228 }
0229
0230 // This test case checks that no segmentation fault appears
0231 // - this tests the loop protection
0232 BOOST_DATA_TEST_CASE(
0233     loop_protection_test,
0234     bdata::random((bdata::engine = std::mt19937(), bdata::seed = 20,
0235                    bdata::distribution = std::uniform_real_distribution<double>(
0236                        0.1_GeV, 0.5_GeV))) ^
0237         bdata::random(
0238             (bdata::engine = std::mt19937(), bdata::seed = 21,
0239              bdata::distribution = std::uniform_real_distribution<double>(
0240                  -std::numbers::pi, std::numbers::pi))) ^
0241         bdata::random(
0242             (bdata::engine = std::mt19937(), bdata::seed = 22,
0243              bdata::distribution = std::uniform_real_distribution<double>(
0244                  1., std::numbers::pi - 1.))) ^
0245         bdata::random((bdata::engine = std::mt19937(), bdata::seed = 23,
0246                        bdata::distribution =
0247                            std::uniform_int_distribution<std::uint8_t>(0, 1))) ^
0248         bdata::xrange(ntests),
0249     pT, phi, theta, charge, index) {
0250   double p = pT / sin(theta);
0251   double q = -1 + 2 * charge;
0252   (void)index;
0253
0254   // define start parameters
0255   /// a covariance matrix to transport
0256   Covariance cov;
0257   // take some major correlations (off-diagonals)
0258   cov << 10_mm, 0, 0.123, 0, 0.5, 0, 0, 10_mm, 0, 0.162, 0, 0, 0.123, 0, 0.1, 0,
0259       0, 0, 0, 0.162, 0, 0.1, 0, 0, 0.5, 0, 0, 0, 1. / (10_GeV), 0, 0, 0, 0, 0,
0260       0, 0;
0261   CurvilinearTrackParameters start(Vector4(0, 0, 0, 0), phi, theta, q / p, cov,
0262                                    ParticleHypothesis::pion());
0263
0264   // Action list and abort list
0265   EigenPropagatorType::Options<ActorList<MaterialInteractor>> options(
0266       tgContext, mfContext);
0267   options.stepping.maxStepSize = 25_cm;
0268   options.pathLimit = 1500_mm;
0269
0270   const auto& status = epropagator.propagate(start, options).value();
0271   // this test assumes state.options.loopFraction = 0.5
0272   // maximum momentum allowed
0273   auto bCache = bField->makeCache(mfContext);
0274   double pmax =
0275       options.pathLimit *
0276       bField->getField(start.position(tgContext), bCache).value().norm() /
0277       std::numbers::pi;
0278   if (p < pmax) {
0279     BOOST_CHECK_LT(status.pathLength, options.pathLimit);
0280   } else {
0281     CHECK_CLOSE_REL(status.pathLength, options.pathLimit, 1e-3);
0282   }
0283 }
0284
0285 }  // namespace Acts::Test