Tests/IntegrationTests/PropagationTests.hpp

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 #pragma once
0010
0011 #include "Acts/Definitions/Direction.hpp"
0012 #include "Acts/Definitions/Units.hpp"
0013 #include "Acts/EventData/TrackParameters.hpp"
0014 #include "Acts/Geometry/CuboidVolumeBuilder.hpp"
0015 #include "Acts/Geometry/GeometryContext.hpp"
0016 #include "Acts/Geometry/TrackingGeometry.hpp"
0017 #include "Acts/Geometry/TrackingGeometryBuilder.hpp"
0018 #include "Acts/MagneticField/MagneticFieldContext.hpp"
0019 #include "Acts/Material/HomogeneousVolumeMaterial.hpp"
0020 #include "Acts/Material/Interactions.hpp"
0021 #include "Acts/Material/Material.hpp"
0022 #include "Acts/Surfaces/CylinderSurface.hpp"
0023 #include "Acts/Surfaces/DiscSurface.hpp"
0024 #include "Acts/Surfaces/PlaneSurface.hpp"
0025 #include "Acts/Surfaces/RectangleBounds.hpp"
0026 #include "Acts/Surfaces/StrawSurface.hpp"
0027 #include "Acts/Utilities/Logger.hpp"
0028 #include "Acts/Utilities/MathHelpers.hpp"
0029 #include "Acts/Utilities/UnitVectors.hpp"
0030 #include "Acts/Utilities/detail/periodic.hpp"
0031 #include "ActsTests/CommonHelpers/FloatComparisons.hpp"
0032 #include "ActsTests/CommonHelpers/PredefinedMaterials.hpp"
0033
0034 #include <utility>
0035
0036 inline std::shared_ptr<const Acts::TrackingGeometry> createDenseBlock(
0037     const Acts::GeometryContext& geoCtx) {
0038   using namespace Acts;
0039   using namespace ActsTests;
0040   using namespace UnitLiterals;
0041
0042   CuboidVolumeBuilder::VolumeConfig vConf;
0043   vConf.position = {0., 0., 0.};
0044   vConf.length = {4_m, 4_m, 4_m};
0045   vConf.volumeMaterial =
0046       std::make_shared<const HomogeneousVolumeMaterial>(makeBeryllium());
0047   CuboidVolumeBuilder::Config conf;
0048   conf.volumeCfg.push_back(vConf);
0049   conf.position = {0., 0., 0.};
0050   conf.length = {4_m, 4_m, 4_m};
0051   CuboidVolumeBuilder cvb(conf);
0052   TrackingGeometryBuilder::Config tgbCfg;
0053   tgbCfg.trackingVolumeBuilders.push_back(
0054       [=](const auto& context, const auto& inner, const auto&) {
0055         return cvb.trackingVolume(context, inner, nullptr);
0056       });
0057
0058   return TrackingGeometryBuilder(tgbCfg).trackingGeometry(geoCtx);
0059 }
0060
0061 inline std::tuple<std::shared_ptr<const Acts::TrackingGeometry>,
0062                   std::vector<const Acts::Surface*>>
0063 createDenseTelescope(const Acts::GeometryContext& geoCtx,
0064                      Acts::Material material, double thickness) {
0065   using namespace Acts;
0066   using namespace UnitLiterals;
0067
0068   CuboidVolumeBuilder::Config conf;
0069   conf.position = {0., 0., 0.};
0070   conf.length = {4_m, 2_m, 2_m};
0071
0072   {
0073     CuboidVolumeBuilder::VolumeConfig start;
0074     start.position = {-1_m, 0, 0};
0075     start.length = {1_m, 1_m, 1_m};
0076     start.name = "start";
0077
0078     conf.volumeCfg.push_back(start);
0079   }
0080
0081   if (thickness < 1_m) {
0082     CuboidVolumeBuilder::VolumeConfig gap;
0083     gap.position = {-0.25 * (1_m + thickness), 0, 0};
0084     gap.length = {0.5 * (1_m - thickness), 1_m, 1_m};
0085     gap.name = "gap1";
0086
0087     conf.volumeCfg.push_back(gap);
0088   }
0089
0090   {
0091     CuboidVolumeBuilder::VolumeConfig dense;
0092     dense.position = {0, 0, 0};
0093     dense.length = {thickness, 1_m, 1_m};
0094     dense.volumeMaterial =
0095         std::make_shared<const HomogeneousVolumeMaterial>(material);
0096     dense.name = "dense";
0097
0098     conf.volumeCfg.push_back(dense);
0099   }
0100
0101   if (thickness < 1_m) {
0102     CuboidVolumeBuilder::VolumeConfig gap;
0103     gap.position = {0.25 * (1_m + thickness), 0, 0};
0104     gap.length = {0.5 * (1_m - thickness), 1_m, 1_m};
0105     gap.name = "gap2";
0106
0107     conf.volumeCfg.push_back(gap);
0108   }
0109
0110   {
0111     CuboidVolumeBuilder::SurfaceConfig surface;
0112     surface.position = {1.5_m, 0, 0};
0113     surface.rotation =
0114         Eigen::AngleAxisd(std::numbers::pi / 2, Eigen::Vector3d::UnitY())
0115             .matrix();
0116     surface.rBounds = std::make_shared<RectangleBounds>(1_m, 1_m);
0117
0118     CuboidVolumeBuilder::LayerConfig layer;
0119     layer.surfaceCfg.push_back(surface);
0120
0121     CuboidVolumeBuilder::VolumeConfig end;
0122     end.position = {1_m, 0, 0};
0123     end.length = {1_m, 1_m, 1_m};
0124     end.layerCfg.push_back(layer);
0125     end.name = "end";
0126
0127     conf.volumeCfg.push_back(end);
0128   }
0129
0130   CuboidVolumeBuilder cvb(conf);
0131
0132   TrackingGeometryBuilder::Config tgbCfg;
0133   tgbCfg.trackingVolumeBuilders.push_back(
0134       [=](const auto& context, const auto& inner, const auto&) {
0135         return cvb.trackingVolume(context, inner, nullptr);
0136       });
0137   auto detector = TrackingGeometryBuilder(tgbCfg).trackingGeometry(geoCtx);
0138
0139   std::vector<const Surface*> surfaces;
0140   detector->visitSurfaces(
0141       [&](const Surface* surface) { surfaces.push_back(surface); });
0142
0143   return {std::move(detector), std::move(surfaces)};
0144 }
0145
0146 // parameter construction helpers
0147
0148 /// Construct (initial) curvilinear parameters.
0149 inline Acts::BoundTrackParameters makeParametersCurvilinear(double phi,
0150                                                             double theta,
0151                                                             double absMom,
0152                                                             double charge) {
0153   using namespace Acts;
0154   using namespace UnitLiterals;
0155
0156   Vector4 pos4 = Vector4::Zero();
0157   auto particleHypothesis = ParticleHypothesis::pionLike(std::abs(charge));
0158   return BoundTrackParameters::createCurvilinear(
0159       pos4, phi, theta, particleHypothesis.qOverP(absMom, charge), std::nullopt,
0160       particleHypothesis);
0161 }
0162
0163 /// Construct (initial) curvilinear parameters with covariance.
0164 inline Acts::BoundTrackParameters makeParametersCurvilinearWithCovariance(
0165     double phi, double theta, double absMom, double charge) {
0166   using namespace Acts;
0167   using namespace UnitLiterals;
0168
0169   BoundVector stddev = BoundVector::Zero();
0170   // TODO use momentum-dependent resolutions
0171   stddev[eBoundLoc0] = 15_um;
0172   stddev[eBoundLoc1] = 80_um;
0173   stddev[eBoundTime] = 20_ps;
0174   stddev[eBoundPhi] = 20_mrad;
0175   stddev[eBoundTheta] = 30_mrad;
0176   stddev[eBoundQOverP] = 1_e / 10_GeV;
0177   BoundSquareMatrix corr = BoundSquareMatrix::Identity();
0178   corr(eBoundLoc0, eBoundLoc1) = corr(eBoundLoc1, eBoundLoc0) = 0.125;
0179   corr(eBoundLoc0, eBoundPhi) = corr(eBoundPhi, eBoundLoc0) = 0.25;
0180   corr(eBoundLoc1, eBoundTheta) = corr(eBoundTheta, eBoundLoc1) = -0.25;
0181   corr(eBoundTime, eBoundQOverP) = corr(eBoundQOverP, eBoundTime) = 0.125;
0182   corr(eBoundPhi, eBoundTheta) = corr(eBoundTheta, eBoundPhi) = -0.25;
0183   corr(eBoundPhi, eBoundQOverP) = corr(eBoundPhi, eBoundQOverP) = -0.125;
0184   corr(eBoundTheta, eBoundQOverP) = corr(eBoundTheta, eBoundQOverP) = 0.5;
0185   BoundSquareMatrix cov = stddev.asDiagonal() * corr * stddev.asDiagonal();
0186
0187   Vector4 pos4 = Vector4::Zero();
0188   auto particleHypothesis = ParticleHypothesis::pionLike(std::abs(charge));
0189   return BoundTrackParameters::createCurvilinear(
0190       pos4, phi, theta, particleHypothesis.qOverP(absMom, charge), cov,
0191       particleHypothesis);
0192 }
0193
0194 /// Construct (initial) neutral curvilinear parameters.
0195 inline Acts::BoundTrackParameters makeParametersCurvilinearNeutral(
0196     double phi, double theta, double absMom) {
0197   using namespace Acts;
0198   using namespace UnitLiterals;
0199
0200   Vector4 pos4 = Vector4::Zero();
0201   return BoundTrackParameters::createCurvilinear(
0202       pos4, phi, theta, 1 / absMom, std::nullopt, ParticleHypothesis::pion0());
0203 }
0204
0205 // helpers to compare track parameters
0206
0207 /// Check that two parameters object are consistent within the tolerances.
0208 ///
0209 /// \warning Does not check that they are defined on the same surface.
0210 inline void checkParametersConsistency(const Acts::BoundTrackParameters& cmp,
0211                                        const Acts::BoundTrackParameters& ref,
0212                                        const Acts::GeometryContext& geoCtx,
0213                                        double epsPos, double epsDir,
0214                                        double epsMom) {
0215   using namespace Acts;
0216
0217   // check stored parameters
0218   CHECK_CLOSE_ABS(cmp.template get<eBoundLoc0>(),
0219                   ref.template get<eBoundLoc0>(), epsPos);
0220   CHECK_CLOSE_ABS(cmp.template get<eBoundLoc1>(),
0221                   ref.template get<eBoundLoc1>(), epsPos);
0222   CHECK_CLOSE_ABS(cmp.template get<eBoundTime>(),
0223                   ref.template get<eBoundTime>(), epsPos);
0224   // check phi equivalence with circularity
0225   CHECK_SMALL(detail::radian_sym(cmp.template get<eBoundPhi>() -
0226                                  ref.template get<eBoundPhi>()),
0227               epsDir);
0228   CHECK_CLOSE_ABS(cmp.template get<eBoundTheta>(),
0229                   ref.template get<eBoundTheta>(), epsDir);
0230   CHECK_CLOSE_ABS(cmp.template get<eBoundQOverP>(),
0231                   ref.template get<eBoundQOverP>(), epsMom);
0232   // check derived parameters
0233   CHECK_CLOSE_ABS(cmp.position(geoCtx), ref.position(geoCtx), epsPos);
0234   CHECK_CLOSE_ABS(cmp.time(), ref.time(), epsPos);
0235   CHECK_CLOSE_ABS(cmp.direction(), ref.direction(), epsDir);
0236   CHECK_CLOSE_ABS(cmp.absoluteMomentum(), ref.absoluteMomentum(), epsMom);
0237   // charge should be identical not just similar
0238   BOOST_CHECK_EQUAL(cmp.charge(), ref.charge());
0239 }
0240
0241 /// Check that two parameters covariances are consistent within the tolerances.
0242 ///
0243 /// \warning Does not check that the parameters value itself are consistent.
0244 inline void checkCovarianceConsistency(const Acts::BoundTrackParameters& cmp,
0245                                        const Acts::BoundTrackParameters& ref,
0246                                        double relativeTolerance) {
0247   // either both or none have covariance set
0248   if (cmp.covariance().has_value()) {
0249     // comparison parameters have covariance but the reference does not
0250     BOOST_CHECK(ref.covariance().has_value());
0251   }
0252   if (ref.covariance().has_value()) {
0253     // reference parameters have covariance but the comparison does not
0254     BOOST_CHECK(cmp.covariance().has_value());
0255   }
0256   if (cmp.covariance().has_value() && ref.covariance().has_value()) {
0257     CHECK_CLOSE_COVARIANCE(cmp.covariance().value(), ref.covariance().value(),
0258                            relativeTolerance);
0259   }
0260 }
0261
0262 // helpers to construct target surfaces from track states
0263
0264 /// Construct the transformation from the curvilinear to the global coordinates.
0265 inline Acts::Transform3 createCurvilinearTransform(
0266     const Acts::BoundTrackParameters& params,
0267     const Acts::GeometryContext& geoCtx) {
0268   using namespace Acts;
0269
0270   Vector3 unitW = params.direction();
0271   auto [unitU, unitV] = createCurvilinearUnitVectors(unitW);
0272
0273   RotationMatrix3 rotation = RotationMatrix3::Zero();
0274   rotation.col(0) = unitU;
0275   rotation.col(1) = unitV;
0276   rotation.col(2) = unitW;
0277   Translation3 offset(params.position(geoCtx));
0278   Transform3 toGlobal = offset * rotation;
0279
0280   return toGlobal;
0281 }
0282
0283 /// Construct a z-cylinder centered at zero with the track on its surface.
0284 struct ZCylinderSurfaceBuilder {
0285   std::shared_ptr<Acts::CylinderSurface> operator()(
0286       const Acts::BoundTrackParameters& params,
0287       const Acts::GeometryContext& geoCtx) {
0288     using namespace Acts;
0289
0290     auto radius = params.position(geoCtx).template head<2>().norm();
0291     auto halfz = std::numeric_limits<double>::max();
0292     return Surface::makeShared<CylinderSurface>(Transform3::Identity(), radius,
0293                                                 halfz);
0294   }
0295 };
0296
0297 /// Construct a disc at track position with plane normal along track tangent.
0298 struct DiscSurfaceBuilder {
0299   std::shared_ptr<Acts::DiscSurface> operator()(
0300       const Acts::BoundTrackParameters& params,
0301       const Acts::GeometryContext& geoCtx) {
0302     using namespace Acts;
0303     using namespace UnitLiterals;
0304
0305     auto cl = createCurvilinearTransform(params, geoCtx);
0306     // shift the origin of the plane so the local particle position does not
0307     // sit directly at the rho=0,phi=undefined singularity
0308     // TODO this is a hack do avoid issues with the numerical covariance
0309     //      transport that does not work well at rho=0,
0310     Vector3 localOffset = Vector3::Zero();
0311     localOffset[ePos0] = 1_cm;
0312     localOffset[ePos1] = -1_cm;
0313     Vector3 globalOriginDelta = cl.linear() * localOffset;
0314     cl.pretranslate(globalOriginDelta);
0315
0316     return Surface::makeShared<DiscSurface>(cl);
0317   }
0318 };
0319
0320 /// Construct a plane at track position with plane normal along track tangent.
0321 struct PlaneSurfaceBuilder {
0322   std::shared_ptr<Acts::PlaneSurface> operator()(
0323       const Acts::BoundTrackParameters& params,
0324       const Acts::GeometryContext& geoCtx) {
0325     using namespace Acts;
0326
0327     return Surface::makeShared<PlaneSurface>(
0328         createCurvilinearTransform(params, geoCtx));
0329   }
0330 };
0331
0332 /// Construct a z-straw at the track position.
0333 struct ZStrawSurfaceBuilder {
0334   std::shared_ptr<Acts::StrawSurface> operator()(
0335       const Acts::BoundTrackParameters& params,
0336       const Acts::GeometryContext& geoCtx) {
0337     using namespace Acts;
0338
0339     return Surface::makeShared<StrawSurface>(
0340         Transform3(Translation3(params.position(geoCtx))));
0341   }
0342 };
0343
0344 // helper functions to run the propagation with additional checks
0345
0346 /// Propagate the initial parameters for the given pathlength in space.
0347 ///
0348 /// Use a negative path length to indicate backward propagation.
0349 template <typename propagator_t,
0350           typename options_t = typename propagator_t::template Options<>>
0351 inline std::pair<Acts::BoundTrackParameters, double> transportFreely(
0352     const propagator_t& propagator, const Acts::GeometryContext& geoCtx,
0353     const Acts::MagneticFieldContext& magCtx,
0354     const Acts::BoundTrackParameters& initialParams, double pathLength) {
0355   using namespace Acts;
0356   using namespace UnitLiterals;
0357
0358   // setup propagation options
0359   options_t options(geoCtx, magCtx);
0360   options.direction = Direction::fromScalar(pathLength);
0361   options.pathLimit = pathLength;
0362   options.surfaceTolerance = 1_nm;
0363   options.stepping.stepTolerance = 1_nm;
0364
0365   auto result = propagator.propagate(initialParams, options);
0366   BOOST_CHECK(result.ok());
0367   BOOST_CHECK(result.value().endParameters);
0368
0369   return {*result.value().endParameters, result.value().pathLength};
0370 }
0371
0372 /// Propagate the initial parameters to the target surface.
0373 template <typename propagator_t,
0374           typename options_t = typename propagator_t::template Options<>>
0375 inline std::pair<Acts::BoundTrackParameters, double> transportToSurface(
0376     const propagator_t& propagator, const Acts::GeometryContext& geoCtx,
0377     const Acts::MagneticFieldContext& magCtx,
0378     const Acts::BoundTrackParameters& initialParams,
0379     const Acts::Surface& targetSurface, double pathLimit) {
0380   using namespace Acts;
0381   using namespace UnitLiterals;
0382
0383   // setup propagation options
0384   options_t options(geoCtx, magCtx);
0385   options.direction = Direction::Forward();
0386   options.pathLimit = pathLimit;
0387   options.surfaceTolerance = 1_nm;
0388   options.stepping.stepTolerance = 1_nm;
0389
0390   auto result = propagator.propagate(initialParams, targetSurface, options);
0391   BOOST_CHECK(result.ok());
0392   BOOST_CHECK(result.value().endParameters);
0393
0394   return {*result.value().endParameters, result.value().pathLength};
0395 }
0396
0397 // self-consistency tests for a single propagator
0398
0399 /// Propagate the initial parameters the given path length along its
0400 /// trajectory and then propagate the final parameters back. Verify that the
0401 /// propagated parameters match the initial ones.
0402 template <typename propagator_t,
0403           typename options_t = typename propagator_t::template Options<>>
0404 inline void runForwardBackwardTest(
0405     const propagator_t& propagator, const Acts::GeometryContext& geoCtx,
0406     const Acts::MagneticFieldContext& magCtx,
0407     const Acts::BoundTrackParameters& initialParams, double pathLength,
0408     double epsPos, double epsDir, double epsMom) {
0409   // propagate parameters Acts::Direction::Forward()
0410   auto [fwdParams, fwdPathLength] = transportFreely<propagator_t, options_t>(
0411       propagator, geoCtx, magCtx, initialParams, pathLength);
0412   CHECK_CLOSE_ABS(fwdPathLength, pathLength, epsPos);
0413   // propagate propagated parameters back again
0414   auto [bwdParams, bwdPathLength] = transportFreely<propagator_t, options_t>(
0415       propagator, geoCtx, magCtx, fwdParams, -pathLength);
0416   CHECK_CLOSE_ABS(bwdPathLength, -pathLength, epsPos);
0417   // check that initial and back-propagated parameters match
0418   checkParametersConsistency(initialParams, bwdParams, geoCtx, epsPos, epsDir,
0419                              epsMom);
0420 }
0421
0422 /// Propagate the initial parameters once for the given path length and
0423 /// use the propagated parameters to define a target surface. Propagate the
0424 /// initial parameters again to the target surface. Verify that the surface has
0425 /// been found and the parameters are consistent.
0426 template <typename propagator_t, typename surface_builder_t,
0427           typename options_t = typename propagator_t::template Options<>>
0428 inline void runToSurfaceTest(const propagator_t& propagator,
0429                              const Acts::GeometryContext& geoCtx,
0430                              const Acts::MagneticFieldContext& magCtx,
0431                              const Acts::BoundTrackParameters& initialParams,
0432                              double pathLength,
0433                              surface_builder_t&& buildTargetSurface,
0434                              double epsPos, double epsDir, double epsMom) {
0435   // free propagation for the given path length
0436   auto [freeParams, freePathLength] = transportFreely<propagator_t, options_t>(
0437       propagator, geoCtx, magCtx, initialParams, pathLength);
0438   CHECK_CLOSE_ABS(freePathLength, pathLength, epsPos);
0439   // build a target surface at the propagated position
0440   auto surface = buildTargetSurface(freeParams, geoCtx);
0441   BOOST_CHECK(surface);
0442
0443   // bound propagation onto the target surface
0444   // increase path length limit to ensure the surface can be reached
0445   auto [surfParams, surfPathLength] =
0446       transportToSurface<propagator_t, options_t>(propagator, geoCtx, magCtx,
0447                                                   initialParams, *surface,
0448                                                   1.5 * pathLength);
0449   CHECK_CLOSE_ABS(surfPathLength, pathLength, epsPos);
0450
0451   // check that the to-surface propagation matches the defining free parameters
0452   CHECK_CLOSE_ABS(surfParams.position(geoCtx), freeParams.position(geoCtx),
0453                   epsPos);
0454   CHECK_CLOSE_ABS(surfParams.time(), freeParams.time(), epsPos);
0455   CHECK_CLOSE_ABS(surfParams.direction(), freeParams.direction(), epsDir);
0456   CHECK_CLOSE_ABS(surfParams.absoluteMomentum(), freeParams.absoluteMomentum(),
0457                   epsMom);
0458   CHECK_CLOSE_ABS(surfPathLength, freePathLength, epsPos);
0459 }
0460
0461 // consistency tests between two propagators
0462
0463 /// Propagate the initial parameters along their trajectory for the given path
0464 /// length using two different propagators and verify consistent output.
0465 template <typename cmp_propagator_t, typename ref_propagator_t>
0466 inline void runForwardComparisonTest(
0467     const cmp_propagator_t& cmpPropagator,
0468     const ref_propagator_t& refPropagator, const Acts::GeometryContext& geoCtx,
0469     const Acts::MagneticFieldContext& magCtx,
0470     const Acts::BoundTrackParameters& initialParams, double pathLength,
0471     double epsPos, double epsDir, double epsMom, double tolCov) {
0472   // propagate twice using the two different propagators
0473   auto [cmpParams, cmpPath] = transportFreely<cmp_propagator_t>(
0474       cmpPropagator, geoCtx, magCtx, initialParams, pathLength);
0475   auto [refParams, refPath] = transportFreely<ref_propagator_t>(
0476       refPropagator, geoCtx, magCtx, initialParams, pathLength);
0477   // check parameter comparison
0478   checkParametersConsistency(cmpParams, refParams, geoCtx, epsPos, epsDir,
0479                              epsMom);
0480   checkCovarianceConsistency(cmpParams, refParams, tolCov);
0481   CHECK_CLOSE_ABS(cmpPath, pathLength, epsPos);
0482   CHECK_CLOSE_ABS(refPath, pathLength, epsPos);
0483   CHECK_CLOSE_ABS(cmpPath, refPath, epsPos);
0484 }
0485
0486 /// Propagate the initial parameters along their trajectory for the given path
0487 /// length using the reference propagator. Use the propagated track parameters
0488 /// to define a target plane. Propagate the initial parameters using two
0489 /// different propagators and verify consistent output.
0490 template <typename cmp_propagator_t, typename ref_propagator_t,
0491           typename surface_builder_t>
0492 inline void runToSurfaceComparisonTest(
0493     const cmp_propagator_t& cmpPropagator,
0494     const ref_propagator_t& refPropagator, const Acts::GeometryContext& geoCtx,
0495     const Acts::MagneticFieldContext& magCtx,
0496     const Acts::BoundTrackParameters& initialParams, double pathLength,
0497     surface_builder_t&& buildTargetSurface, double epsPos, double epsDir,
0498     double epsMom, double tolCov) {
0499   // free propagation with the reference propagator for the given path length
0500   auto [freeParams, freePathLength] = transportFreely<ref_propagator_t>(
0501       refPropagator, geoCtx, magCtx, initialParams, pathLength);
0502   CHECK_CLOSE_ABS(freePathLength, pathLength, epsPos);
0503
0504   // build a target surface at the propagated position
0505   auto surface = buildTargetSurface(freeParams, geoCtx);
0506   BOOST_CHECK(surface);
0507
0508   // propagate twice to the surface using the two different propagators
0509   // increase path length limit to ensure the surface can be reached
0510   auto [cmpParams, cmpPath] = transportToSurface<cmp_propagator_t>(
0511       cmpPropagator, geoCtx, magCtx, initialParams, *surface, 1.5 * pathLength);
0512   auto [refParams, refPath] = transportToSurface<ref_propagator_t>(
0513       refPropagator, geoCtx, magCtx, initialParams, *surface, 1.5 * pathLength);
0514   // check parameter comparison
0515   checkParametersConsistency(cmpParams, refParams, geoCtx, epsPos, epsDir,
0516                              epsMom);
0517   checkCovarianceConsistency(cmpParams, refParams, tolCov);
0518   CHECK_CLOSE_ABS(cmpPath, pathLength, epsPos);
0519   CHECK_CLOSE_ABS(refPath, pathLength, epsPos);
0520   CHECK_CLOSE_ABS(cmpPath, refPath, epsPos);
0521 }
0522
0523 template <typename propagator_t>
0524 inline void runDenseForwardTest(
0525     const propagator_t& propagator, const Acts::GeometryContext& geoCtx,
0526     const Acts::MagneticFieldContext& magCtx, double p, double q,
0527     const Acts::Surface& targetSurface, const Acts::Material& material,
0528     double thickness, const Acts::Logger& logger = Acts::getDummyLogger()) {
0529   using namespace Acts;
0530   using namespace UnitLiterals;
0531
0532   const auto particleHypothesis = ParticleHypothesis::muon();
0533   const float mass = particleHypothesis.mass();
0534   const float absQ = std::abs(q);
0535   const double initialQOverP = particleHypothesis.qOverP(p, q);
0536   const auto initialParams = BoundTrackParameters::createCurvilinear(
0537       Vector4(-1.5_m, 0, 0, 0), Vector3::UnitX(), initialQOverP,
0538       BoundVector::Constant(1e-16).asDiagonal(), particleHypothesis);
0539
0540   typename propagator_t::template Options<> options(geoCtx, magCtx);
0541   options.maxSteps = 10000;
0542   options.stepping.maxStepSize = 100_mm;
0543   options.stepping.dense.meanEnergyLoss = true;
0544   options.stepping.dense.momentumCutOff = 1_MeV;
0545
0546   const auto result =
0547       propagator.propagate(initialParams, targetSurface, options);
0548
0549   BOOST_CHECK(result.ok());
0550   CHECK_CLOSE_REL(3_m, result->pathLength, 1e-6);
0551   BOOST_CHECK(result->endParameters);
0552
0553   BoundTrackParameters endParams = result->endParameters.value();
0554
0555   BOOST_CHECK(endParams.covariance());
0556   CHECK_CLOSE_ABS(initialParams.position(geoCtx) + Acts::Vector3(3_m, 0, 0),
0557                   endParams.position(geoCtx), 1e-6);
0558   CHECK_CLOSE_ABS(initialParams.direction(), endParams.direction(), 1e-6);
0559
0560   const auto& cov = endParams.covariance().value();
0561
0562   const double endQOverP = endParams.qOverP();
0563   const double endP = endParams.absoluteMomentum();
0564   const double endEloss =
0565       Acts::fastHypot(p, mass) - Acts::fastHypot(endP, mass);
0566   const double endErrX = std::sqrt(cov(eBoundLoc0, eBoundLoc0));
0567   const double endErrY = std::sqrt(cov(eBoundLoc1, eBoundLoc1));
0568   const double endErrQOverP = std::sqrt(cov(eBoundQOverP, eBoundQOverP));
0569   const double endErrP =
0570       (absQ / std::pow(endParams.qOverP(), 2)) * endErrQOverP;
0571   const double endErrE = (p / Acts::fastHypot(p, mass)) * endErrP;
0572   const double endErrTheta = std::sqrt(cov(eBoundTheta, eBoundTheta));
0573   const double endErrPhi = std::sqrt(cov(eBoundPhi, eBoundPhi));
0574
0575   ACTS_VERBOSE("input q/p = " << initialQOverP);
0576   ACTS_VERBOSE("input p = " << p);
0577   ACTS_VERBOSE("output q/p = " << endQOverP);
0578   ACTS_VERBOSE("output p = " << endP);
0579   ACTS_VERBOSE("output t = " << endParams.time());
0580   ACTS_VERBOSE("output eloss = " << endEloss);
0581   ACTS_VERBOSE("output err x = " << endErrX);
0582   ACTS_VERBOSE("output err y = " << endErrY);
0583   ACTS_VERBOSE("output err q/p = " << endErrQOverP);
0584   ACTS_VERBOSE("output err p = " << endErrP);
0585   ACTS_VERBOSE("output err E = " << endErrE);
0586   ACTS_VERBOSE("output err theta = " << endErrTheta);
0587   ACTS_VERBOSE("output err phi = " << endErrPhi);
0588
0589   BOOST_CHECK_CLOSE(endErrX, endErrY, 1e-6);
0590   BOOST_CHECK_CLOSE(endErrTheta, endErrPhi, 1e-6);
0591
0592   const float elossInitial = computeEnergyLossMean(
0593       MaterialSlab(material, thickness), particleHypothesis.absolutePdg(), mass,
0594       initialQOverP, absQ);
0595   const float elossFinal = computeEnergyLossMean(
0596       MaterialSlab(material, thickness), particleHypothesis.absolutePdg(), mass,
0597       endQOverP, absQ);
0598
0599   ACTS_VERBOSE("ref eloss initial = " << elossInitial);
0600   ACTS_VERBOSE("ref eloss final = " << elossFinal);
0601
0602   // energy loss will be smaller than elossInitial because of the increasing
0603   // radiation losses
0604   BOOST_CHECK_LT(endEloss, elossInitial);
0605   BOOST_CHECK_GT(endEloss, elossFinal);
0606
0607   const float theta0Initial = computeMultipleScatteringTheta0(
0608       MaterialSlab(material, thickness), particleHypothesis.absolutePdg(), mass,
0609       initialQOverP, absQ);
0610   const float theta0Final = computeMultipleScatteringTheta0(
0611       MaterialSlab(material, thickness), particleHypothesis.absolutePdg(), mass,
0612       endQOverP, absQ);
0613
0614   ACTS_VERBOSE("ref theta0 initial = " << theta0Initial);
0615   ACTS_VERBOSE("ref theta0 final = " << theta0Final);
0616
0617   // angle errors will be larger than theta0Initial because of the energy loss
0618   BOOST_CHECK_GE(endErrTheta, theta0Initial);
0619   BOOST_CHECK_LT(endErrTheta, theta0Final);
0620   CHECK_CLOSE_REL(endErrTheta, 0.5 * (theta0Initial + theta0Final), 0.1);
0621
0622   // estimates the positional uncertainty after crossing the material
0623   // block by integrating the angle variance over substeps
0624   const auto estimateVarX = [&](double qOverP, std::size_t substeps) {
0625     double previousTheta0 = 0;
0626     double varAngle = 0;
0627     double varPosition = 0;
0628     double covAnglePosition = 0;
0629
0630     const auto step = [&](double substep, double theta0) {
0631       double deltaVarTheta = square(theta0) - square(previousTheta0);
0632       double deltaVarPos = varAngle * square(substep) +
0633                            2 * covAnglePosition * substep +
0634                            deltaVarTheta * (square(substep) / 3);
0635       double deltaCovAnglePosition =
0636           varAngle * substep + deltaVarTheta * substep / 2;
0637       previousTheta0 = theta0;
0638       varAngle += deltaVarTheta;
0639       varPosition += deltaVarPos;
0640       covAnglePosition += deltaCovAnglePosition;
0641     };
0642
0643     // step through the material block
0644     for (std::size_t i = 0; i < substeps; ++i) {
0645       const float theta0 = computeMultipleScatteringTheta0(
0646           MaterialSlab(material, (i + 1) * thickness / substeps),
0647           particleHypothesis.absolutePdg(), mass, qOverP, absQ);
0648       step(thickness / substeps, theta0);
0649     }
0650
0651     // step to the target surface
0652     step(1_m, previousTheta0);
0653
0654     return varPosition;
0655   };
0656
0657   // as a lower bound for sigma_x we use the initial theta0
0658   const double lowerBoundVarX = estimateVarX(initialQOverP, 10);
0659   // as an upper bound for sigma_x we use the final theta0
0660   const double upperBoundVarX = estimateVarX(endQOverP, 10);
0661
0662   const double lowerBoundSigmaX = std::sqrt(lowerBoundVarX);
0663   const double upperBoundSigmaX = std::sqrt(upperBoundVarX);
0664
0665   ACTS_VERBOSE("ref lower bound sigma x = " << lowerBoundSigmaX);
0666   ACTS_VERBOSE("ref upper bound sigma x = " << upperBoundSigmaX);
0667
0668   BOOST_CHECK_GT(endErrX, lowerBoundSigmaX);
0669   BOOST_CHECK_LT(endErrX, upperBoundSigmaX);
0670   CHECK_CLOSE_REL(endErrX, 0.5 * (lowerBoundSigmaX + upperBoundSigmaX), 0.2);
0671
0672   const float qOverPSigma = computeEnergyLossLandauSigmaQOverP(
0673       MaterialSlab(material, thickness), mass, initialQOverP, absQ);
0674
0675   ACTS_VERBOSE("ref sigma q/p = " << qOverPSigma);
0676
0677   CHECK_CLOSE_REL(endErrQOverP, qOverPSigma, 1e-4);
0678 }