Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-10-30 07:55:26 

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/Utilities/Logger.hpp"
0033 #include "Acts/Utilities/Result.hpp"
0034 #include "ActsTests/CommonHelpers/FloatComparisons.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 
0048 using namespace Acts;
0049 using namespace Acts::UnitLiterals;
0050 using Acts::VectorHelpers::makeVector4;
0051 using Acts::VectorHelpers::perp;
0052 
0053 namespace ActsTests {
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 BOOST_AUTO_TEST_SUITE(PropagatorSuite)
0149 
0150 // This tests the Options
0151 BOOST_AUTO_TEST_CASE(PropagatorOptions_) {
0152   using NullOptionsType = EigenPropagatorType::Options<>;
0153   NullOptionsType null_options(tgContext, mfContext);
0154 
0155   using ActorList = ActorList<PerpendicularMeasure>;
0156   using OptionsType = EigenPropagatorType::Options<ActorList>;
0157   OptionsType options(tgContext, mfContext);
0158 }
0159 
0160 BOOST_DATA_TEST_CASE(
0161     cylinder_passage_observer_,
0162     bdata::random((bdata::engine = std::mt19937(), bdata::seed = 0,
0163                    bdata::distribution = std::uniform_real_distribution<double>(
0164                        0.4_GeV, 10_GeV))) ^
0165         bdata::random(
0166             (bdata::engine = std::mt19937(), bdata::seed = 1,
0167              bdata::distribution = std::uniform_real_distribution<double>(
0168                  -std::numbers::pi, std::numbers::pi))) ^
0169         bdata::random(
0170             (bdata::engine = std::mt19937(), bdata::seed = 2,
0171              bdata::distribution = std::uniform_real_distribution<double>(
0172                  1., std::numbers::pi - 1.))) ^
0173         bdata::random((bdata::engine = std::mt19937(), bdata::seed = 3,
0174                        bdata::distribution =
0175                            std::uniform_int_distribution<std::uint8_t>(0, 1))) ^
0176         bdata::random((bdata::engine = std::mt19937(), bdata::seed = 4,
0177                        bdata::distribution =
0178                            std::uniform_real_distribution<double>(-1_ns,
0179                                                                   1_ns))) ^
0180         bdata::xrange(ntests),
0181     pT, phi, theta, charge, time, index) {
0182   double dcharge = -1 + 2 * charge;
0183   (void)index;
0184 
0185   using CylinderObserver = SurfaceObserver<CylinderSurface>;
0186   using ActorList = ActorList<CylinderObserver>;
0187 
0188   // setup propagation options
0189   EigenPropagatorType::Options<ActorList> options(tgContext, mfContext);
0190 
0191   options.pathLimit = 20_m;
0192   options.stepping.maxStepSize = 1_cm;
0193 
0194   // set the surface to be passed by
0195   options.actorList.get<CylinderObserver>().surface = mSurface.get();
0196 
0197   using so_result = typename CylinderObserver::result_type;
0198 
0199   // define start parameters
0200   double x = 0;
0201   double y = 0;
0202   double z = 0;
0203   double px = pT * cos(phi);
0204   double py = pT * sin(phi);
0205   double pz = pT / tan(theta);
0206   double q = dcharge;
0207   Vector3 pos(x, y, z);
0208   Vector3 mom(px, py, pz);
0209   BoundTrackParameters start = BoundTrackParameters::createCurvilinear(
0210       makeVector4(pos, time), mom.normalized(), q / mom.norm(), std::nullopt,
0211       ParticleHypothesis::pion());
0212   // propagate to the cylinder surface
0213   const auto& result = epropagator.propagate(start, *cSurface, options).value();
0214   auto& sor = result.get<so_result>();
0215 
0216   BOOST_CHECK_EQUAL(sor.surfaces_passed, 1u);
0217   CHECK_CLOSE_ABS(sor.surface_passed_r, 10., 1e-5);
0218 }
0219 
0220 BOOST_DATA_TEST_CASE(
0221     curvilinear_additive_,
0222     bdata::random((bdata::engine = std::mt19937(), bdata::seed = 0,
0223                    bdata::distribution = std::uniform_real_distribution<double>(
0224                        0.4_GeV, 10_GeV))) ^
0225         bdata::random(
0226             (bdata::engine = std::mt19937(), bdata::seed = 1,
0227              bdata::distribution = std::uniform_real_distribution<double>(
0228                  -std::numbers::pi, std::numbers::pi))) ^
0229         bdata::random(
0230             (bdata::engine = std::mt19937(), bdata::seed = 2,
0231              bdata::distribution = std::uniform_real_distribution<double>(
0232                  1., std::numbers::pi - 1.))) ^
0233         bdata::random((bdata::engine = std::mt19937(), bdata::seed = 3,
0234                        bdata::distribution =
0235                            std::uniform_int_distribution<std::uint8_t>(0, 1))) ^
0236         bdata::random((bdata::engine = std::mt19937(), bdata::seed = 4,
0237                        bdata::distribution =
0238                            std::uniform_real_distribution<double>(-1_ns,
0239                                                                   1_ns))) ^
0240         bdata::xrange(ntests),
0241     pT, phi, theta, charge, time, index) {
0242   double dcharge = -1 + 2 * charge;
0243   (void)index;
0244 
0245   // setup propagation options - the tow step options
0246   EigenPropagatorType::Options<> options_2s(tgContext, mfContext);
0247   options_2s.pathLimit = 50_cm;
0248   options_2s.stepping.maxStepSize = 1_cm;
0249 
0250   // define start parameters
0251   double x = 0;
0252   double y = 0;
0253   double z = 0;
0254   double px = pT * cos(phi);
0255   double py = pT * sin(phi);
0256   double pz = pT / tan(theta);
0257   double q = dcharge;
0258   Vector3 pos(x, y, z);
0259   Vector3 mom(px, py, pz);
0260   /// a covariance matrix to transport
0261   Covariance cov;
0262   // take some major correlations (off-diagonals)
0263   cov << 10_mm, 0, 0.123, 0, 0.5, 0, 0, 10_mm, 0, 0.162, 0, 0, 0.123, 0, 0.1, 0,
0264       0, 0, 0, 0.162, 0, 0.1, 0, 0, 0.5, 0, 0, 0, 1. / (10_GeV), 0, 0, 0, 0, 0,
0265       0, 0;
0266   BoundTrackParameters start = BoundTrackParameters::createCurvilinear(
0267       makeVector4(pos, time), mom.normalized(), q / mom.norm(), cov,
0268       ParticleHypothesis::pion());
0269   // propagate to a path length of 100 with two steps of 50
0270   const auto& mid_parameters =
0271       epropagator.propagate(start, options_2s).value().endParameters;
0272   const auto& end_parameters_2s =
0273       epropagator.propagate(*mid_parameters, options_2s).value().endParameters;
0274 
0275   // setup propagation options - the one step options
0276   EigenPropagatorType::Options<> options_1s(tgContext, mfContext);
0277   options_1s.pathLimit = 100_cm;
0278   options_1s.stepping.maxStepSize = 1_cm;
0279   // propagate to a path length of 100 in one step
0280   const auto& end_parameters_1s =
0281       epropagator.propagate(start, options_1s).value().endParameters;
0282 
0283   // test that the propagation is additive
0284   CHECK_CLOSE_REL(end_parameters_1s->position(tgContext),
0285                   end_parameters_2s->position(tgContext), 0.001);
0286 
0287   BOOST_CHECK(end_parameters_1s->covariance().has_value());
0288   const auto& cov_1s = *(end_parameters_1s->covariance());
0289   BOOST_CHECK(end_parameters_2s->covariance().has_value());
0290   const auto& cov_2s = *(end_parameters_2s->covariance());
0291 
0292   // CHECK_CLOSE_COVARIANCE(cov_1s, cov_2s, 0.001);
0293   for (unsigned int i = 0; i < cov_1s.rows(); i++) {
0294     for (unsigned int j = 0; j < cov_1s.cols(); j++) {
0295       CHECK_CLOSE_OR_SMALL(cov_1s(i, j), cov_2s(i, j), 0.001, 1e-6);
0296     }
0297   }
0298 }
0299 
0300 BOOST_DATA_TEST_CASE(
0301     cylinder_additive_,
0302     bdata::random((bdata::engine = std::mt19937(), bdata::seed = 0,
0303                    bdata::distribution = std::uniform_real_distribution<double>(
0304                        0.4_GeV, 10_GeV))) ^
0305         bdata::random(
0306             (bdata::engine = std::mt19937(), bdata::seed = 1,
0307              bdata::distribution = std::uniform_real_distribution<double>(
0308                  -std::numbers::pi, std::numbers::pi))) ^
0309         bdata::random(
0310             (bdata::engine = std::mt19937(), bdata::seed = 2,
0311              bdata::distribution = std::uniform_real_distribution<double>(
0312                  1., std::numbers::pi - 1.))) ^
0313         bdata::random((bdata::engine = std::mt19937(), bdata::seed = 3,
0314                        bdata::distribution =
0315                            std::uniform_int_distribution<std::uint8_t>(0, 1))) ^
0316         bdata::random((bdata::engine = std::mt19937(), bdata::seed = 4,
0317                        bdata::distribution =
0318                            std::uniform_real_distribution<double>(-1_ns,
0319                                                                   1_ns))) ^
0320         bdata::xrange(ntests),
0321     pT, phi, theta, charge, time, index) {
0322   double dcharge = -1 + 2 * charge;
0323   (void)index;
0324 
0325   // setup propagation options - 2 setp options
0326   EigenPropagatorType::Options<> options_2s(tgContext, mfContext);
0327   options_2s.pathLimit = 10_m;
0328   options_2s.stepping.maxStepSize = 1_cm;
0329 
0330   // define start parameters
0331   double x = 0;
0332   double y = 0;
0333   double z = 0;
0334   double px = pT * cos(phi);
0335   double py = pT * sin(phi);
0336   double pz = pT / tan(theta);
0337   double q = dcharge;
0338   Vector3 pos(x, y, z);
0339   Vector3 mom(px, py, pz);
0340   /// a covariance matrix to transport
0341   Covariance cov;
0342   // take some major correlations (off-diagonals)
0343   cov << 10_mm, 0, 0.123, 0, 0.5, 0, 0, 10_mm, 0, 0.162, 0, 0, 0.123, 0, 0.1, 0,
0344       0, 0, 0, 0.162, 0, 0.1, 0, 0, 0.5, 0, 0, 0, 1. / (10_GeV), 0, 0, 0, 0, 0,
0345       0, 0;
0346   BoundTrackParameters start = BoundTrackParameters::createCurvilinear(
0347       makeVector4(pos, time), mom.normalized(), q / mom.norm(), cov,
0348       ParticleHypothesis::pion());
0349   // propagate to a final surface with one stop in between
0350   const auto& mid_parameters =
0351       epropagator.propagate(start, *mSurface, options_2s).value().endParameters;
0352 
0353   const auto& end_parameters_2s =
0354       epropagator.propagate(*mid_parameters, *cSurface, options_2s)
0355           .value()
0356           .endParameters;
0357 
0358   // setup propagation options - one step options
0359   EigenPropagatorType::Options<> options_1s(tgContext, mfContext);
0360   options_1s.pathLimit = 10_m;
0361   options_1s.stepping.maxStepSize = 1_cm;
0362   // propagate to a final surface in one stop
0363   const auto& end_parameters_1s =
0364       epropagator.propagate(start, *cSurface, options_1s).value().endParameters;
0365 
0366   // test that the propagation is additive
0367   CHECK_CLOSE_REL(end_parameters_1s->position(tgContext),
0368                   end_parameters_2s->position(tgContext), 0.001);
0369 
0370   BOOST_CHECK(end_parameters_1s->covariance().has_value());
0371   const auto& cov_1s = (*(end_parameters_1s->covariance()));
0372   BOOST_CHECK(end_parameters_2s->covariance().has_value());
0373   const auto& cov_2s = (*(end_parameters_2s->covariance()));
0374 
0375   // CHECK_CLOSE_COVARIANCE(cov_1s, cov_2s, 0.001);
0376   for (unsigned int i = 0; i < cov_1s.rows(); i++) {
0377     for (unsigned int j = 0; j < cov_1s.cols(); j++) {
0378       CHECK_CLOSE_OR_SMALL(cov_1s(i, j), cov_2s(i, j), 0.001, 1e-6);
0379     }
0380   }
0381 }
0382 
0383 BOOST_AUTO_TEST_CASE(BasicPropagatorInterface) {
0384   auto field = std::make_shared<ConstantBField>(Vector3{0, 0, 2_T});
0385   EigenStepper<> eigenStepper{field};
0386   VoidNavigator navigator{};
0387 
0388   std::shared_ptr<PlaneSurface> startSurface =
0389       CurvilinearSurface(Vector3::Zero(), Vector3::UnitX()).planeSurface();
0390   std::shared_ptr<PlaneSurface> targetSurface =
0391       CurvilinearSurface(Vector3::UnitX() * 20_mm, Vector3::UnitX())
0392           .planeSurface();
0393 
0394   BoundVector startPars;
0395   startPars << 0, 0, 0, std::numbers::pi / 2., 1 / 1_GeV, 0;
0396 
0397   BoundTrackParameters startParameters{startSurface, startPars, std::nullopt,
0398                                        ParticleHypothesis::pion()};
0399 
0400   BoundTrackParameters startCurv = BoundTrackParameters::createCurvilinear(
0401       Vector4::Zero(), Vector3::UnitX(), 1. / 1_GeV, std::nullopt,
0402       ParticleHypothesis::pion());
0403 
0404   GeometryContext gctx;
0405   MagneticFieldContext mctx;
0406 
0407   {
0408     EigenPropagatorType::Options<> options{gctx, mctx};
0409 
0410     Propagator propagator{eigenStepper, navigator};
0411     static_assert(std::is_base_of_v<BasePropagator, decltype(propagator)>,
0412                   "Propagator does not inherit from BasePropagator");
0413     const BasePropagator* base =
0414         static_cast<const BasePropagator*>(&propagator);
0415 
0416     // Ensure the propagation does the same thing
0417     auto result =
0418         propagator.propagate(startParameters, *targetSurface, options);
0419     BOOST_REQUIRE(result.ok());
0420     BOOST_CHECK_EQUAL(&result.value().endParameters.value().referenceSurface(),
0421                       targetSurface.get());
0422 
0423     auto resultBase =
0424         base->propagateToSurface(startParameters, *targetSurface,
0425                                  static_cast<PropagatorPlainOptions>(options));
0426 
0427     BOOST_REQUIRE(resultBase.ok());
0428     BOOST_CHECK_EQUAL(&resultBase.value().referenceSurface(),
0429                       targetSurface.get());
0430 
0431     BOOST_CHECK_EQUAL(result.value().endParameters.value().parameters(),
0432                       resultBase.value().parameters());
0433 
0434     // Propagation call with curvilinear also works
0435     auto resultCurv =
0436         base->propagateToSurface(startCurv, *targetSurface,
0437                                  static_cast<PropagatorPlainOptions>(options));
0438     BOOST_CHECK(resultCurv.ok());
0439   }
0440 
0441   StraightLineStepper slStepper{};
0442   {
0443     Propagator<StraightLineStepper>::Options<> options{gctx, mctx};
0444 
0445     Propagator propagator{slStepper, navigator};
0446     static_assert(std::is_base_of_v<BasePropagator, decltype(propagator)>,
0447                   "Propagator does not inherit from BasePropagator");
0448     const BasePropagator* base =
0449         static_cast<const BasePropagator*>(&propagator);
0450 
0451     // Ensure the propagation does the same thing
0452     auto result =
0453         propagator.propagate(startParameters, *targetSurface, options);
0454     BOOST_REQUIRE(result.ok());
0455     BOOST_CHECK_EQUAL(&result.value().endParameters.value().referenceSurface(),
0456                       targetSurface.get());
0457 
0458     auto resultBase =
0459         base->propagateToSurface(startParameters, *targetSurface,
0460                                  static_cast<PropagatorPlainOptions>(options));
0461 
0462     BOOST_REQUIRE(resultBase.ok());
0463     BOOST_CHECK_EQUAL(&resultBase.value().referenceSurface(),
0464                       targetSurface.get());
0465 
0466     BOOST_CHECK_EQUAL(result.value().endParameters.value().parameters(),
0467                       resultBase.value().parameters());
0468 
0469     // Propagation call with curvilinear also works
0470     auto resultCurv =
0471         base->propagateToSurface(startCurv, *targetSurface,
0472                                  static_cast<PropagatorPlainOptions>(options));
0473     BOOST_CHECK(resultCurv.ok());
0474   }
0475 
0476   EigenStepper<EigenStepperDenseExtension> denseEigenStepper{field};
0477 
0478   {
0479     Propagator propagator{denseEigenStepper, navigator};
0480     static_assert(!std::is_base_of_v<BasePropagator, decltype(propagator)>,
0481                   "Propagator unexpectedly inherits from BasePropagator");
0482   }
0483 }
0484 
0485 BOOST_AUTO_TEST_SUITE_END()
0486 
0487 }  // namespace ActsTests
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!