Core/Propagator/PropagatorTests.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/ConstrainedStep.hpp"
0022 #include "Acts/Propagator/EigenStepper.hpp"
0023 #include "Acts/Propagator/EigenStepperDenseExtension.hpp"
0024 #include "Acts/Propagator/Propagator.hpp"
0025 #include "Acts/Propagator/StraightLineStepper.hpp"
0026 #include "Acts/Propagator/VoidNavigator.hpp"
0027 #include "Acts/Surfaces/CurvilinearSurface.hpp"
0028 #include "Acts/Surfaces/CylinderBounds.hpp"
0029 #include "Acts/Surfaces/CylinderSurface.hpp"
0030 #include "Acts/Surfaces/PlaneSurface.hpp"
0031 #include "Acts/Surfaces/Surface.hpp"
0032 #include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
0033 #include "Acts/Utilities/Logger.hpp"
0034 #include "Acts/Utilities/Result.hpp"
0035
0036 #include <cmath>
0037 #include <cstddef>
0038 #include <limits>
0039 #include <memory>
0040 #include <numbers>
0041 #include <optional>
0042 #include <random>
0043 #include <type_traits>
0044 #include <utility>
0045
0046 namespace bdata = boost::unit_test::data;
0047 using namespace Acts::UnitLiterals;
0048 using Acts::VectorHelpers::makeVector4;
0049 using Acts::VectorHelpers::perp;
0050
0051 namespace Acts::Test {
0052
0053 // Create a test context
0054 GeometryContext tgContext = GeometryContext();
0055 MagneticFieldContext mfContext = MagneticFieldContext();
0056
0057 using Covariance = BoundSquareMatrix;
0058
0059 /// An observer that measures the perpendicular distance
0060 struct PerpendicularMeasure {
0061   /// Simple result struct to be returned
0062   struct this_result {
0063     double distance = std::numeric_limits<double>::max();
0064   };
0065
0066   using result_type = this_result;
0067
0068   PerpendicularMeasure() = default;
0069
0070   template <typename propagator_state_t, typename stepper_t,
0071             typename navigator_t>
0072   void operator()(propagator_state_t& state, const stepper_t& stepper,
0073                   const navigator_t& /*navigator*/, result_type& result) const {
0074     result.distance = perp(stepper.position(state.stepping));
0075   }
0076 };
0077
0078 /// An observer that measures the perpendicular distance
0079 template <typename Surface>
0080 struct SurfaceObserver {
0081   // the surface to be intersected
0082   const Surface* surface = nullptr;
0083   // the tolerance for intersection
0084   double tolerance = 1e-5;
0085
0086   /// Simple result struct to be returned
0087   struct this_result {
0088     std::size_t surfaces_passed = 0;
0089     double surface_passed_r = std::numeric_limits<double>::max();
0090   };
0091
0092   using result_type = this_result;
0093
0094   template <typename propagator_state_t, typename stepper_t,
0095             typename navigator_t>
0096   void act(propagator_state_t& state, const stepper_t& stepper,
0097            const navigator_t& /*navigator*/, result_type& result,
0098            const Logger& /*logger*/) const {
0099     if (surface == nullptr || result.surfaces_passed != 0) {
0100       return;
0101     }
0102
0103     // calculate the distance to the surface
0104     const double distance =
0105         surface
0106             ->intersect(state.geoContext, stepper.position(state.stepping),
0107                         stepper.direction(state.stepping),
0108                         BoundaryTolerance::None())
0109             .closest()
0110             .pathLength();
0111
0112     // Adjust the step size so that we cannot cross the target surface
0113     state.stepping.stepSize.release(ConstrainedStep::Type::Actor);
0114     state.stepping.stepSize.update(distance * state.options.direction,
0115                                    ConstrainedStep::Type::Actor);
0116
0117     // return true if you fall below tolerance
0118     if (std::abs(distance) <= tolerance) {
0119       ++result.surfaces_passed;
0120       result.surface_passed_r = perp(stepper.position(state.stepping));
0121       state.stepping.stepSize.release(ConstrainedStep::Type::Actor);
0122     }
0123   }
0124 };
0125
0126 // Global definitions
0127 using BFieldType = ConstantBField;
0128 using EigenStepperType = EigenStepper<>;
0129 using EigenPropagatorType = Propagator<EigenStepperType>;
0130
0131 const double Bz = 2_T;
0132 auto bField = std::make_shared<BFieldType>(Vector3{0, 0, Bz});
0133 EigenStepperType estepper(bField);
0134 EigenPropagatorType epropagator(std::move(estepper), VoidNavigator(),
0135                                 getDefaultLogger("prop", Logging::VERBOSE));
0136
0137 auto mCylinder = std::make_shared<CylinderBounds>(10_mm, 1000_mm);
0138 auto mSurface =
0139     Surface::makeShared<CylinderSurface>(Transform3::Identity(), mCylinder);
0140 auto cCylinder = std::make_shared<CylinderBounds>(150_mm, 1000_mm);
0141 auto cSurface =
0142     Surface::makeShared<CylinderSurface>(Transform3::Identity(), cCylinder);
0143
0144 const int ntests = 5;
0145
0146 // This tests the Options
0147 BOOST_AUTO_TEST_CASE(PropagatorOptions_) {
0148   using NullOptionsType = EigenPropagatorType::Options<>;
0149   NullOptionsType null_options(tgContext, mfContext);
0150
0151   using ActorList = ActorList<PerpendicularMeasure>;
0152   using OptionsType = EigenPropagatorType::Options<ActorList>;
0153   OptionsType options(tgContext, mfContext);
0154 }
0155
0156 BOOST_DATA_TEST_CASE(
0157     cylinder_passage_observer_,
0158     bdata::random((bdata::engine = std::mt19937(), bdata::seed = 0,
0159                    bdata::distribution = std::uniform_real_distribution<double>(
0160                        0.4_GeV, 10_GeV))) ^
0161         bdata::random(
0162             (bdata::engine = std::mt19937(), bdata::seed = 1,
0163              bdata::distribution = std::uniform_real_distribution<double>(
0164                  -std::numbers::pi, std::numbers::pi))) ^
0165         bdata::random(
0166             (bdata::engine = std::mt19937(), bdata::seed = 2,
0167              bdata::distribution = std::uniform_real_distribution<double>(
0168                  1., std::numbers::pi - 1.))) ^
0169         bdata::random((bdata::engine = std::mt19937(), bdata::seed = 3,
0170                        bdata::distribution =
0171                            std::uniform_int_distribution<std::uint8_t>(0, 1))) ^
0172         bdata::random((bdata::engine = std::mt19937(), bdata::seed = 4,
0173                        bdata::distribution =
0174                            std::uniform_real_distribution<double>(-1_ns,
0175                                                                   1_ns))) ^
0176         bdata::xrange(ntests),
0177     pT, phi, theta, charge, time, index) {
0178   double dcharge = -1 + 2 * charge;
0179   (void)index;
0180
0181   using CylinderObserver = SurfaceObserver<CylinderSurface>;
0182   using ActorList = ActorList<CylinderObserver>;
0183
0184   // setup propagation options
0185   EigenPropagatorType::Options<ActorList> options(tgContext, mfContext);
0186
0187   options.pathLimit = 20_m;
0188   options.stepping.maxStepSize = 1_cm;
0189
0190   // set the surface to be passed by
0191   options.actorList.get<CylinderObserver>().surface = mSurface.get();
0192
0193   using so_result = typename CylinderObserver::result_type;
0194
0195   // define start parameters
0196   double x = 0;
0197   double y = 0;
0198   double z = 0;
0199   double px = pT * cos(phi);
0200   double py = pT * sin(phi);
0201   double pz = pT / tan(theta);
0202   double q = dcharge;
0203   Vector3 pos(x, y, z);
0204   Vector3 mom(px, py, pz);
0205   BoundTrackParameters start = BoundTrackParameters::createCurvilinear(
0206       makeVector4(pos, time), mom.normalized(), q / mom.norm(), std::nullopt,
0207       ParticleHypothesis::pion());
0208   // propagate to the cylinder surface
0209   const auto& result = epropagator.propagate(start, *cSurface, options).value();
0210   auto& sor = result.get<so_result>();
0211
0212   BOOST_CHECK_EQUAL(sor.surfaces_passed, 1u);
0213   CHECK_CLOSE_ABS(sor.surface_passed_r, 10., 1e-5);
0214 }
0215
0216 BOOST_DATA_TEST_CASE(
0217     curvilinear_additive_,
0218     bdata::random((bdata::engine = std::mt19937(), bdata::seed = 0,
0219                    bdata::distribution = std::uniform_real_distribution<double>(
0220                        0.4_GeV, 10_GeV))) ^
0221         bdata::random(
0222             (bdata::engine = std::mt19937(), bdata::seed = 1,
0223              bdata::distribution = std::uniform_real_distribution<double>(
0224                  -std::numbers::pi, std::numbers::pi))) ^
0225         bdata::random(
0226             (bdata::engine = std::mt19937(), bdata::seed = 2,
0227              bdata::distribution = std::uniform_real_distribution<double>(
0228                  1., std::numbers::pi - 1.))) ^
0229         bdata::random((bdata::engine = std::mt19937(), bdata::seed = 3,
0230                        bdata::distribution =
0231                            std::uniform_int_distribution<std::uint8_t>(0, 1))) ^
0232         bdata::random((bdata::engine = std::mt19937(), bdata::seed = 4,
0233                        bdata::distribution =
0234                            std::uniform_real_distribution<double>(-1_ns,
0235                                                                   1_ns))) ^
0236         bdata::xrange(ntests),
0237     pT, phi, theta, charge, time, index) {
0238   double dcharge = -1 + 2 * charge;
0239   (void)index;
0240
0241   // setup propagation options - the tow step options
0242   EigenPropagatorType::Options<> options_2s(tgContext, mfContext);
0243   options_2s.pathLimit = 50_cm;
0244   options_2s.stepping.maxStepSize = 1_cm;
0245
0246   // define start parameters
0247   double x = 0;
0248   double y = 0;
0249   double z = 0;
0250   double px = pT * cos(phi);
0251   double py = pT * sin(phi);
0252   double pz = pT / tan(theta);
0253   double q = dcharge;
0254   Vector3 pos(x, y, z);
0255   Vector3 mom(px, py, pz);
0256   /// a covariance matrix to transport
0257   Covariance cov;
0258   // take some major correlations (off-diagonals)
0259   cov << 10_mm, 0, 0.123, 0, 0.5, 0, 0, 10_mm, 0, 0.162, 0, 0, 0.123, 0, 0.1, 0,
0260       0, 0, 0, 0.162, 0, 0.1, 0, 0, 0.5, 0, 0, 0, 1. / (10_GeV), 0, 0, 0, 0, 0,
0261       0, 0;
0262   BoundTrackParameters start = BoundTrackParameters::createCurvilinear(
0263       makeVector4(pos, time), mom.normalized(), q / mom.norm(), cov,
0264       ParticleHypothesis::pion());
0265   // propagate to a path length of 100 with two steps of 50
0266   const auto& mid_parameters =
0267       epropagator.propagate(start, options_2s).value().endParameters;
0268   const auto& end_parameters_2s =
0269       epropagator.propagate(*mid_parameters, options_2s).value().endParameters;
0270
0271   // setup propagation options - the one step options
0272   EigenPropagatorType::Options<> options_1s(tgContext, mfContext);
0273   options_1s.pathLimit = 100_cm;
0274   options_1s.stepping.maxStepSize = 1_cm;
0275   // propagate to a path length of 100 in one step
0276   const auto& end_parameters_1s =
0277       epropagator.propagate(start, options_1s).value().endParameters;
0278
0279   // test that the propagation is additive
0280   CHECK_CLOSE_REL(end_parameters_1s->position(tgContext),
0281                   end_parameters_2s->position(tgContext), 0.001);
0282
0283   BOOST_CHECK(end_parameters_1s->covariance().has_value());
0284   const auto& cov_1s = *(end_parameters_1s->covariance());
0285   BOOST_CHECK(end_parameters_2s->covariance().has_value());
0286   const auto& cov_2s = *(end_parameters_2s->covariance());
0287
0288   // CHECK_CLOSE_COVARIANCE(cov_1s, cov_2s, 0.001);
0289   for (unsigned int i = 0; i < cov_1s.rows(); i++) {
0290     for (unsigned int j = 0; j < cov_1s.cols(); j++) {
0291       CHECK_CLOSE_OR_SMALL(cov_1s(i, j), cov_2s(i, j), 0.001, 1e-6);
0292     }
0293   }
0294 }
0295
0296 BOOST_DATA_TEST_CASE(
0297     cylinder_additive_,
0298     bdata::random((bdata::engine = std::mt19937(), bdata::seed = 0,
0299                    bdata::distribution = std::uniform_real_distribution<double>(
0300                        0.4_GeV, 10_GeV))) ^
0301         bdata::random(
0302             (bdata::engine = std::mt19937(), bdata::seed = 1,
0303              bdata::distribution = std::uniform_real_distribution<double>(
0304                  -std::numbers::pi, std::numbers::pi))) ^
0305         bdata::random(
0306             (bdata::engine = std::mt19937(), bdata::seed = 2,
0307              bdata::distribution = std::uniform_real_distribution<double>(
0308                  1., std::numbers::pi - 1.))) ^
0309         bdata::random((bdata::engine = std::mt19937(), bdata::seed = 3,
0310                        bdata::distribution =
0311                            std::uniform_int_distribution<std::uint8_t>(0, 1))) ^
0312         bdata::random((bdata::engine = std::mt19937(), bdata::seed = 4,
0313                        bdata::distribution =
0314                            std::uniform_real_distribution<double>(-1_ns,
0315                                                                   1_ns))) ^
0316         bdata::xrange(ntests),
0317     pT, phi, theta, charge, time, index) {
0318   double dcharge = -1 + 2 * charge;
0319   (void)index;
0320
0321   // setup propagation options - 2 setp options
0322   EigenPropagatorType::Options<> options_2s(tgContext, mfContext);
0323   options_2s.pathLimit = 10_m;
0324   options_2s.stepping.maxStepSize = 1_cm;
0325
0326   // define start parameters
0327   double x = 0;
0328   double y = 0;
0329   double z = 0;
0330   double px = pT * cos(phi);
0331   double py = pT * sin(phi);
0332   double pz = pT / tan(theta);
0333   double q = dcharge;
0334   Vector3 pos(x, y, z);
0335   Vector3 mom(px, py, pz);
0336   /// a covariance matrix to transport
0337   Covariance cov;
0338   // take some major correlations (off-diagonals)
0339   cov << 10_mm, 0, 0.123, 0, 0.5, 0, 0, 10_mm, 0, 0.162, 0, 0, 0.123, 0, 0.1, 0,
0340       0, 0, 0, 0.162, 0, 0.1, 0, 0, 0.5, 0, 0, 0, 1. / (10_GeV), 0, 0, 0, 0, 0,
0341       0, 0;
0342   BoundTrackParameters start = BoundTrackParameters::createCurvilinear(
0343       makeVector4(pos, time), mom.normalized(), q / mom.norm(), cov,
0344       ParticleHypothesis::pion());
0345   // propagate to a final surface with one stop in between
0346   const auto& mid_parameters =
0347       epropagator.propagate(start, *mSurface, options_2s).value().endParameters;
0348
0349   const auto& end_parameters_2s =
0350       epropagator.propagate(*mid_parameters, *cSurface, options_2s)
0351           .value()
0352           .endParameters;
0353
0354   // setup propagation options - one step options
0355   EigenPropagatorType::Options<> options_1s(tgContext, mfContext);
0356   options_1s.pathLimit = 10_m;
0357   options_1s.stepping.maxStepSize = 1_cm;
0358   // propagate to a final surface in one stop
0359   const auto& end_parameters_1s =
0360       epropagator.propagate(start, *cSurface, options_1s).value().endParameters;
0361
0362   // test that the propagation is additive
0363   CHECK_CLOSE_REL(end_parameters_1s->position(tgContext),
0364                   end_parameters_2s->position(tgContext), 0.001);
0365
0366   BOOST_CHECK(end_parameters_1s->covariance().has_value());
0367   const auto& cov_1s = (*(end_parameters_1s->covariance()));
0368   BOOST_CHECK(end_parameters_2s->covariance().has_value());
0369   const auto& cov_2s = (*(end_parameters_2s->covariance()));
0370
0371   // CHECK_CLOSE_COVARIANCE(cov_1s, cov_2s, 0.001);
0372   for (unsigned int i = 0; i < cov_1s.rows(); i++) {
0373     for (unsigned int j = 0; j < cov_1s.cols(); j++) {
0374       CHECK_CLOSE_OR_SMALL(cov_1s(i, j), cov_2s(i, j), 0.001, 1e-6);
0375     }
0376   }
0377 }
0378
0379 BOOST_AUTO_TEST_CASE(BasicPropagatorInterface) {
0380   auto field = std::make_shared<ConstantBField>(Vector3{0, 0, 2_T});
0381   EigenStepper<> eigenStepper{field};
0382   VoidNavigator navigator{};
0383
0384   std::shared_ptr<PlaneSurface> startSurface =
0385       CurvilinearSurface(Vector3::Zero(), Vector3::UnitX()).planeSurface();
0386   std::shared_ptr<PlaneSurface> targetSurface =
0387       CurvilinearSurface(Vector3::UnitX() * 20_mm, Vector3::UnitX())
0388           .planeSurface();
0389
0390   BoundVector startPars;
0391   startPars << 0, 0, 0, std::numbers::pi / 2., 1 / 1_GeV, 0;
0392
0393   BoundTrackParameters startParameters{startSurface, startPars, std::nullopt,
0394                                        ParticleHypothesis::pion()};
0395
0396   BoundTrackParameters startCurv = BoundTrackParameters::createCurvilinear(
0397       Vector4::Zero(), Vector3::UnitX(), 1. / 1_GeV, std::nullopt,
0398       ParticleHypothesis::pion());
0399
0400   GeometryContext gctx;
0401   MagneticFieldContext mctx;
0402
0403   {
0404     EigenPropagatorType::Options<> options{gctx, mctx};
0405
0406     Propagator propagator{eigenStepper, navigator};
0407     static_assert(std::is_base_of_v<BasePropagator, decltype(propagator)>,
0408                   "Propagator does not inherit from BasePropagator");
0409     const BasePropagator* base =
0410         static_cast<const BasePropagator*>(&propagator);
0411
0412     // Ensure the propagation does the same thing
0413     auto result =
0414         propagator.propagate(startParameters, *targetSurface, options);
0415     BOOST_REQUIRE(result.ok());
0416     BOOST_CHECK_EQUAL(&result.value().endParameters.value().referenceSurface(),
0417                       targetSurface.get());
0418
0419     auto resultBase =
0420         base->propagateToSurface(startParameters, *targetSurface,
0421                                  static_cast<PropagatorPlainOptions>(options));
0422
0423     BOOST_REQUIRE(resultBase.ok());
0424     BOOST_CHECK_EQUAL(&resultBase.value().referenceSurface(),
0425                       targetSurface.get());
0426
0427     BOOST_CHECK_EQUAL(result.value().endParameters.value().parameters(),
0428                       resultBase.value().parameters());
0429
0430     // Propagation call with curvilinear also works
0431     auto resultCurv =
0432         base->propagateToSurface(startCurv, *targetSurface,
0433                                  static_cast<PropagatorPlainOptions>(options));
0434     BOOST_CHECK(resultCurv.ok());
0435   }
0436
0437   StraightLineStepper slStepper{};
0438   {
0439     Propagator<StraightLineStepper>::Options<> options{gctx, mctx};
0440
0441     Propagator propagator{slStepper, navigator};
0442     static_assert(std::is_base_of_v<BasePropagator, decltype(propagator)>,
0443                   "Propagator does not inherit from BasePropagator");
0444     const BasePropagator* base =
0445         static_cast<const BasePropagator*>(&propagator);
0446
0447     // Ensure the propagation does the same thing
0448     auto result =
0449         propagator.propagate(startParameters, *targetSurface, options);
0450     BOOST_REQUIRE(result.ok());
0451     BOOST_CHECK_EQUAL(&result.value().endParameters.value().referenceSurface(),
0452                       targetSurface.get());
0453
0454     auto resultBase =
0455         base->propagateToSurface(startParameters, *targetSurface,
0456                                  static_cast<PropagatorPlainOptions>(options));
0457
0458     BOOST_REQUIRE(resultBase.ok());
0459     BOOST_CHECK_EQUAL(&resultBase.value().referenceSurface(),
0460                       targetSurface.get());
0461
0462     BOOST_CHECK_EQUAL(result.value().endParameters.value().parameters(),
0463                       resultBase.value().parameters());
0464
0465     // Propagation call with curvilinear also works
0466     auto resultCurv =
0467         base->propagateToSurface(startCurv, *targetSurface,
0468                                  static_cast<PropagatorPlainOptions>(options));
0469     BOOST_CHECK(resultCurv.ok());
0470   }
0471
0472   EigenStepper<EigenStepperDenseExtension> denseEigenStepper{field};
0473
0474   {
0475     Propagator propagator{denseEigenStepper, navigator};
0476     static_assert(!std::is_base_of_v<BasePropagator, decltype(propagator)>,
0477                   "Propagator unexpectedly inherits from BasePropagator");
0478   }
0479 }
0480 }  // namespace Acts::Test