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