Core/Propagator/PropagatorTests.cpp

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