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