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