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/GeometryContext.hpp"
0015 #include "Acts/MagneticField/MagneticFieldContext.hpp"
0016 #include "Acts/Surfaces/CylinderSurface.hpp"
0017 #include "Acts/Surfaces/DiscSurface.hpp"
0018 #include "Acts/Surfaces/PlaneSurface.hpp"
0019 #include "Acts/Surfaces/StrawSurface.hpp"
0020 #include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
0021 #include "Acts/Utilities/UnitVectors.hpp"
0022 #include "Acts/Utilities/detail/periodic.hpp"
0023
0024 #include <numbers>
0025 #include <utility>
0026
0027 // parameter construction helpers
0028
0029 /// Construct (initial) curvilinear parameters.
0030 inline Acts::CurvilinearTrackParameters makeParametersCurvilinear(
0031     double phi, double theta, double absMom, double charge) {
0032   using namespace Acts;
0033   using namespace Acts::UnitLiterals;
0034
0035   // phi is ill-defined in forward/backward tracks. normalize the value to
0036   // ensure parameter comparisons give correct answers.
0037   if (!((0 < theta) && (theta < std::numbers::pi))) {
0038     phi = 0;
0039   }
0040
0041   Vector4 pos4 = Vector4::Zero();
0042   auto particleHypothesis = ParticleHypothesis::pionLike(std::abs(charge));
0043   return CurvilinearTrackParameters(pos4, phi, theta,
0044                                     particleHypothesis.qOverP(absMom, charge),
0045                                     std::nullopt, particleHypothesis);
0046 }
0047
0048 /// Construct (initial) curvilinear parameters with covariance.
0049 inline Acts::CurvilinearTrackParameters makeParametersCurvilinearWithCovariance(
0050     double phi, double theta, double absMom, double charge) {
0051   using namespace Acts;
0052   using namespace Acts::UnitLiterals;
0053
0054   // phi is ill-defined in forward/backward tracks. normalize the value to
0055   // ensure parameter comparisons give correct answers.
0056   if (!((0 < theta) && (theta < std::numbers::pi))) {
0057     phi = 0;
0058   }
0059
0060   BoundVector stddev = BoundVector::Zero();
0061   // TODO use momentum-dependent resolutions
0062   stddev[eBoundLoc0] = 15_um;
0063   stddev[eBoundLoc1] = 80_um;
0064   stddev[eBoundTime] = 25_ns;
0065   stddev[eBoundPhi] = 1_degree;
0066   stddev[eBoundTheta] = 1.5_degree;
0067   stddev[eBoundQOverP] = 1_e / 10_GeV;
0068   BoundSquareMatrix corr = BoundSquareMatrix::Identity();
0069   corr(eBoundLoc0, eBoundLoc1) = corr(eBoundLoc1, eBoundLoc0) = 0.125;
0070   corr(eBoundLoc0, eBoundPhi) = corr(eBoundPhi, eBoundLoc0) = 0.25;
0071   corr(eBoundLoc1, eBoundTheta) = corr(eBoundTheta, eBoundLoc1) = -0.25;
0072   corr(eBoundTime, eBoundQOverP) = corr(eBoundQOverP, eBoundTime) = 0.125;
0073   corr(eBoundPhi, eBoundTheta) = corr(eBoundTheta, eBoundPhi) = -0.25;
0074   corr(eBoundPhi, eBoundQOverP) = corr(eBoundPhi, eBoundQOverP) = -0.125;
0075   corr(eBoundTheta, eBoundQOverP) = corr(eBoundTheta, eBoundQOverP) = 0.5;
0076   BoundSquareMatrix cov = stddev.asDiagonal() * corr * stddev.asDiagonal();
0077
0078   Vector4 pos4 = Vector4::Zero();
0079   auto particleHypothesis = ParticleHypothesis::pionLike(std::abs(charge));
0080   return CurvilinearTrackParameters(pos4, phi, theta,
0081                                     particleHypothesis.qOverP(absMom, charge),
0082                                     cov, particleHypothesis);
0083 }
0084
0085 /// Construct (initial) neutral curvilinear parameters.
0086 inline Acts::CurvilinearTrackParameters makeParametersCurvilinearNeutral(
0087     double phi, double theta, double absMom) {
0088   using namespace Acts;
0089   using namespace Acts::UnitLiterals;
0090
0091   // phi is ill-defined in forward/backward tracks. normalize the value to
0092   // ensure parameter comparisons give correct answers.
0093   if (!((0 < theta) && (theta < std::numbers::pi))) {
0094     phi = 0;
0095   }
0096
0097   Vector4 pos4 = Vector4::Zero();
0098   return CurvilinearTrackParameters(pos4, phi, theta, 1 / absMom, std::nullopt,
0099                                     ParticleHypothesis::pion0());
0100 }
0101
0102 // helpers to compare track parameters
0103
0104 /// Check that two parameters object are consistent within the tolerances.
0105 ///
0106 /// \warning Does not check that they are defined on the same surface.
0107 inline void checkParametersConsistency(const Acts::BoundTrackParameters& cmp,
0108                                        const Acts::BoundTrackParameters& ref,
0109                                        const Acts::GeometryContext& geoCtx,
0110                                        double epsPos, double epsDir,
0111                                        double epsMom) {
0112   using namespace Acts;
0113
0114   // check stored parameters
0115   CHECK_CLOSE_ABS(cmp.template get<eBoundLoc0>(),
0116                   ref.template get<eBoundLoc0>(), epsPos);
0117   CHECK_CLOSE_ABS(cmp.template get<eBoundLoc1>(),
0118                   ref.template get<eBoundLoc1>(), epsPos);
0119   CHECK_CLOSE_ABS(cmp.template get<eBoundTime>(),
0120                   ref.template get<eBoundTime>(), epsPos);
0121   // check phi equivalence with circularity
0122   CHECK_SMALL(detail::radian_sym(cmp.template get<eBoundPhi>() -
0123                                  ref.template get<eBoundPhi>()),
0124               epsDir);
0125   CHECK_CLOSE_ABS(cmp.template get<eBoundTheta>(),
0126                   ref.template get<eBoundTheta>(), epsDir);
0127   CHECK_CLOSE_ABS(cmp.template get<eBoundQOverP>(),
0128                   ref.template get<eBoundQOverP>(), epsMom);
0129   // check derived parameters
0130   CHECK_CLOSE_ABS(cmp.position(geoCtx), ref.position(geoCtx), epsPos);
0131   CHECK_CLOSE_ABS(cmp.time(), ref.time(), epsPos);
0132   CHECK_CLOSE_ABS(cmp.direction(), ref.direction(), epsDir);
0133   CHECK_CLOSE_ABS(cmp.absoluteMomentum(), ref.absoluteMomentum(), epsMom);
0134   // charge should be identical not just similar
0135   BOOST_CHECK_EQUAL(cmp.charge(), ref.charge());
0136 }
0137
0138 /// Check that two parameters covariances are consistent within the tolerances.
0139 ///
0140 /// \warning Does not check that the parameters value itself are consistent.
0141 inline void checkCovarianceConsistency(const Acts::BoundTrackParameters& cmp,
0142                                        const Acts::BoundTrackParameters& ref,
0143                                        double relativeTolerance) {
0144   // either both or none have covariance set
0145   if (cmp.covariance().has_value()) {
0146     // comparison parameters have covariance but the reference does not
0147     BOOST_CHECK(ref.covariance().has_value());
0148   }
0149   if (ref.covariance().has_value()) {
0150     // reference parameters have covariance but the comparison does not
0151     BOOST_CHECK(cmp.covariance().has_value());
0152   }
0153   if (cmp.covariance().has_value() && ref.covariance().has_value()) {
0154     CHECK_CLOSE_COVARIANCE(cmp.covariance().value(), ref.covariance().value(),
0155                            relativeTolerance);
0156   }
0157 }
0158
0159 // helpers to construct target surfaces from track states
0160
0161 /// Construct the transformation from the curvilinear to the global coordinates.
0162 inline Acts::Transform3 makeCurvilinearTransform(
0163     const Acts::BoundTrackParameters& params,
0164     const Acts::GeometryContext& geoCtx) {
0165   Acts::Vector3 unitW = params.direction();
0166   auto [unitU, unitV] = Acts::makeCurvilinearUnitVectors(unitW);
0167
0168   Acts::RotationMatrix3 rotation = Acts::RotationMatrix3::Zero();
0169   rotation.col(0) = unitU;
0170   rotation.col(1) = unitV;
0171   rotation.col(2) = unitW;
0172   Acts::Translation3 offset(params.position(geoCtx));
0173   Acts::Transform3 toGlobal = offset * rotation;
0174
0175   return toGlobal;
0176 }
0177
0178 /// Construct a z-cylinder centered at zero with the track on its surface.
0179 struct ZCylinderSurfaceBuilder {
0180   std::shared_ptr<Acts::CylinderSurface> operator()(
0181       const Acts::BoundTrackParameters& params,
0182       const Acts::GeometryContext& geoCtx) {
0183     auto radius = params.position(geoCtx).template head<2>().norm();
0184     auto halfz = std::numeric_limits<double>::max();
0185     return Acts::Surface::makeShared<Acts::CylinderSurface>(
0186         Acts::Transform3::Identity(), radius, halfz);
0187   }
0188 };
0189
0190 /// Construct a disc at track position with plane normal along track tangent.
0191 struct DiscSurfaceBuilder {
0192   std::shared_ptr<Acts::DiscSurface> operator()(
0193       const Acts::BoundTrackParameters& params,
0194       const Acts::GeometryContext& geoCtx) {
0195     using namespace Acts;
0196     using namespace Acts::UnitLiterals;
0197
0198     auto cl = makeCurvilinearTransform(params, geoCtx);
0199     // shift the origin of the plane so the local particle position does not
0200     // sit directly at the rho=0,phi=undefined singularity
0201     // TODO this is a hack do avoid issues with the numerical covariance
0202     //      transport that does not work well at rho=0,
0203     Acts::Vector3 localOffset = Acts::Vector3::Zero();
0204     localOffset[Acts::ePos0] = 1_cm;
0205     localOffset[Acts::ePos1] = -1_cm;
0206     Acts::Vector3 globalOriginDelta = cl.linear() * localOffset;
0207     cl.pretranslate(globalOriginDelta);
0208
0209     return Acts::Surface::makeShared<Acts::DiscSurface>(cl);
0210   }
0211 };
0212
0213 /// Construct a plane at track position with plane normal along track tangent.
0214 struct PlaneSurfaceBuilder {
0215   std::shared_ptr<Acts::PlaneSurface> operator()(
0216       const Acts::BoundTrackParameters& params,
0217       const Acts::GeometryContext& geoCtx) {
0218     return Acts::Surface::makeShared<Acts::PlaneSurface>(
0219         makeCurvilinearTransform(params, geoCtx));
0220   }
0221 };
0222
0223 /// Construct a z-straw at the track position.
0224 struct ZStrawSurfaceBuilder {
0225   std::shared_ptr<Acts::StrawSurface> operator()(
0226       const Acts::BoundTrackParameters& params,
0227       const Acts::GeometryContext& geoCtx) {
0228     return Acts::Surface::makeShared<Acts::StrawSurface>(
0229         Acts::Transform3(Acts::Translation3(params.position(geoCtx))));
0230   }
0231 };
0232
0233 // helper functions to run the propagation with additional checks
0234
0235 /// Propagate the initial parameters for the given pathlength in space.
0236 ///
0237 /// Use a negative path length to indicate backward propagation.
0238 template <typename propagator_t,
0239           typename options_t = typename propagator_t::template Options<>>
0240 inline std::pair<Acts::CurvilinearTrackParameters, double> transportFreely(
0241     const propagator_t& propagator, const Acts::GeometryContext& geoCtx,
0242     const Acts::MagneticFieldContext& magCtx,
0243     const Acts::CurvilinearTrackParameters& initialParams, double pathLength) {
0244   using namespace Acts::UnitLiterals;
0245
0246   // setup propagation options
0247   options_t options(geoCtx, magCtx);
0248   options.direction = Acts::Direction::fromScalar(pathLength);
0249   options.pathLimit = pathLength;
0250   options.surfaceTolerance = 1_nm;
0251   options.stepping.stepTolerance = 1_nm;
0252
0253   auto result = propagator.propagate(initialParams, options);
0254   BOOST_CHECK(result.ok());
0255   BOOST_CHECK(result.value().endParameters);
0256
0257   return {*result.value().endParameters, result.value().pathLength};
0258 }
0259
0260 /// Propagate the initial parameters to the target surface.
0261 template <typename propagator_t,
0262           typename options_t = typename propagator_t::template Options<>>
0263 inline std::pair<Acts::BoundTrackParameters, double> transportToSurface(
0264     const propagator_t& propagator, const Acts::GeometryContext& geoCtx,
0265     const Acts::MagneticFieldContext& magCtx,
0266     const Acts::CurvilinearTrackParameters& initialParams,
0267     const Acts::Surface& targetSurface, double pathLimit) {
0268   using namespace Acts::UnitLiterals;
0269
0270   // setup propagation options
0271   options_t options(geoCtx, magCtx);
0272   options.direction = Acts::Direction::Forward();
0273   options.pathLimit = pathLimit;
0274   options.surfaceTolerance = 1_nm;
0275   options.stepping.stepTolerance = 1_nm;
0276
0277   auto result = propagator.propagate(initialParams, targetSurface, options);
0278   BOOST_CHECK(result.ok());
0279   BOOST_CHECK(result.value().endParameters);
0280
0281   return {*result.value().endParameters, result.value().pathLength};
0282 }
0283
0284 // self-consistency tests for a single propagator
0285
0286 /// Propagate the initial parameters the given path length along its
0287 /// trajectory and then propagate the final parameters back. Verify that the
0288 /// propagated parameters match the initial ones.
0289 template <typename propagator_t,
0290           typename options_t = typename propagator_t::template Options<>>
0291 inline void runForwardBackwardTest(
0292     const propagator_t& propagator, const Acts::GeometryContext& geoCtx,
0293     const Acts::MagneticFieldContext& magCtx,
0294     const Acts::CurvilinearTrackParameters& initialParams, double pathLength,
0295     double epsPos, double epsDir, double epsMom) {
0296   // propagate parameters Acts::Direction::Forward()
0297   auto [fwdParams, fwdPathLength] = transportFreely<propagator_t, options_t>(
0298       propagator, geoCtx, magCtx, initialParams, pathLength);
0299   CHECK_CLOSE_ABS(fwdPathLength, pathLength, epsPos);
0300   // propagate propagated parameters back again
0301   auto [bwdParams, bwdPathLength] = transportFreely<propagator_t, options_t>(
0302       propagator, geoCtx, magCtx, fwdParams, -pathLength);
0303   CHECK_CLOSE_ABS(bwdPathLength, -pathLength, epsPos);
0304   // check that initial and back-propagated parameters match
0305   checkParametersConsistency(initialParams, bwdParams, geoCtx, epsPos, epsDir,
0306                              epsMom);
0307 }
0308
0309 /// Propagate the initial parameters once for the given path length and
0310 /// use the propagated parameters to define a target surface. Propagate the
0311 /// initial parameters again to the target surface. Verify that the surface has
0312 /// been found and the parameters are consistent.
0313 template <typename propagator_t, typename surface_builder_t,
0314           typename options_t = typename propagator_t::template Options<>>
0315 inline void runToSurfaceTest(
0316     const propagator_t& propagator, const Acts::GeometryContext& geoCtx,
0317     const Acts::MagneticFieldContext& magCtx,
0318     const Acts::CurvilinearTrackParameters& initialParams, double pathLength,
0319     surface_builder_t&& buildTargetSurface, double epsPos, double epsDir,
0320     double epsMom) {
0321   // free propagation for the given path length
0322   auto [freeParams, freePathLength] = transportFreely<propagator_t, options_t>(
0323       propagator, geoCtx, magCtx, initialParams, pathLength);
0324   CHECK_CLOSE_ABS(freePathLength, pathLength, epsPos);
0325   // build a target surface at the propagated position
0326   auto surface = buildTargetSurface(freeParams, geoCtx);
0327   BOOST_CHECK(surface);
0328
0329   // bound propagation onto the target surface
0330   // increase path length limit to ensure the surface can be reached
0331   auto [surfParams, surfPathLength] =
0332       transportToSurface<propagator_t, options_t>(propagator, geoCtx, magCtx,
0333                                                   initialParams, *surface,
0334                                                   1.5 * pathLength);
0335   CHECK_CLOSE_ABS(surfPathLength, pathLength, epsPos);
0336
0337   // check that the to-surface propagation matches the defining free parameters
0338   CHECK_CLOSE_ABS(surfParams.position(geoCtx), freeParams.position(geoCtx),
0339                   epsPos);
0340   CHECK_CLOSE_ABS(surfParams.time(), freeParams.time(), epsPos);
0341   CHECK_CLOSE_ABS(surfParams.direction(), freeParams.direction(), epsDir);
0342   CHECK_CLOSE_ABS(surfParams.absoluteMomentum(), freeParams.absoluteMomentum(),
0343                   epsMom);
0344   CHECK_CLOSE_ABS(surfPathLength, freePathLength, epsPos);
0345 }
0346
0347 // consistency tests between two propagators
0348
0349 /// Propagate the initial parameters along their trajectory for the given path
0350 /// length using two different propagators and verify consistent output.
0351 template <typename cmp_propagator_t, typename ref_propagator_t>
0352 inline void runForwardComparisonTest(
0353     const cmp_propagator_t& cmpPropagator,
0354     const ref_propagator_t& refPropagator, const Acts::GeometryContext& geoCtx,
0355     const Acts::MagneticFieldContext& magCtx,
0356     const Acts::CurvilinearTrackParameters& initialParams, double pathLength,
0357     double epsPos, double epsDir, double epsMom, double tolCov) {
0358   // propagate twice using the two different propagators
0359   auto [cmpParams, cmpPath] = transportFreely<cmp_propagator_t>(
0360       cmpPropagator, geoCtx, magCtx, initialParams, pathLength);
0361   auto [refParams, refPath] = transportFreely<ref_propagator_t>(
0362       refPropagator, geoCtx, magCtx, initialParams, pathLength);
0363   // check parameter comparison
0364   checkParametersConsistency(cmpParams, refParams, geoCtx, epsPos, epsDir,
0365                              epsMom);
0366   checkCovarianceConsistency(cmpParams, refParams, tolCov);
0367   CHECK_CLOSE_ABS(cmpPath, pathLength, epsPos);
0368   CHECK_CLOSE_ABS(refPath, pathLength, epsPos);
0369   CHECK_CLOSE_ABS(cmpPath, refPath, epsPos);
0370 }
0371
0372 /// Propagate the initial parameters along their trajectory for the given path
0373 /// length using the reference propagator. Use the propagated track parameters
0374 /// to define a target plane. Propagate the initial parameters using two
0375 /// different propagators and verify consistent output.
0376 template <typename cmp_propagator_t, typename ref_propagator_t,
0377           typename surface_builder_t>
0378 inline void runToSurfaceComparisonTest(
0379     const cmp_propagator_t& cmpPropagator,
0380     const ref_propagator_t& refPropagator, const Acts::GeometryContext& geoCtx,
0381     const Acts::MagneticFieldContext& magCtx,
0382     const Acts::CurvilinearTrackParameters& initialParams, double pathLength,
0383     surface_builder_t&& buildTargetSurface, double epsPos, double epsDir,
0384     double epsMom, double tolCov) {
0385   // free propagation with the reference propagator for the given path length
0386   auto [freeParams, freePathLength] = transportFreely<ref_propagator_t>(
0387       refPropagator, geoCtx, magCtx, initialParams, pathLength);
0388   CHECK_CLOSE_ABS(freePathLength, pathLength, epsPos);
0389
0390   // build a target surface at the propagated position
0391   auto surface = buildTargetSurface(freeParams, geoCtx);
0392   BOOST_CHECK(surface);
0393
0394   // propagate twice to the surface using the two different propagators
0395   // increase path length limit to ensure the surface can be reached
0396   auto [cmpParams, cmpPath] = transportToSurface<cmp_propagator_t>(
0397       cmpPropagator, geoCtx, magCtx, initialParams, *surface, 1.5 * pathLength);
0398   auto [refParams, refPath] = transportToSurface<ref_propagator_t>(
0399       refPropagator, geoCtx, magCtx, initialParams, *surface, 1.5 * pathLength);
0400   // check parameter comparison
0401   checkParametersConsistency(cmpParams, refParams, geoCtx, epsPos, epsDir,
0402                              epsMom);
0403   checkCovarianceConsistency(cmpParams, refParams, tolCov);
0404   CHECK_CLOSE_ABS(cmpPath, pathLength, epsPos);
0405   CHECK_CLOSE_ABS(refPath, pathLength, epsPos);
0406   CHECK_CLOSE_ABS(cmpPath, refPath, epsPos);
0407 }