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