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