Core/Propagator/CovarianceEngineTests.cpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 #include <boost/test/data/test_case.hpp>
0010 #include <boost/test/unit_test.hpp>
0011
0012 #include "Acts/Definitions/Algebra.hpp"
0013 #include "Acts/Definitions/TrackParametrization.hpp"
0014 #include "Acts/EventData/ParticleHypothesis.hpp"
0015 #include "Acts/EventData/TrackParameters.hpp"
0016 #include "Acts/EventData/TransformationHelpers.hpp"
0017 #include "Acts/EventData/detail/CorrectedTransformationFreeToBound.hpp"
0018 #include "Acts/Geometry/GeometryContext.hpp"
0019 #include "Acts/MagneticField/ConstantBField.hpp"
0020 #include "Acts/MagneticField/MagneticFieldContext.hpp"
0021 #include "Acts/Propagator/EigenStepper.hpp"
0022 #include "Acts/Propagator/Propagator.hpp"
0023 #include "Acts/Propagator/VoidNavigator.hpp"
0024 #include "Acts/Propagator/detail/CovarianceEngine.hpp"
0025 #include "Acts/Surfaces/PerigeeSurface.hpp"
0026 #include "Acts/Surfaces/PlaneSurface.hpp"
0027 #include "Acts/Surfaces/Surface.hpp"
0028 #include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
0029 #include "Acts/Utilities/Result.hpp"
0030
0031 #include <cmath>
0032 #include <memory>
0033 #include <numbers>
0034 #include <optional>
0035 #include <random>
0036 #include <tuple>
0037 #include <utility>
0038
0039 namespace bdata = boost::unit_test::data;
0040
0041 namespace Acts::Test {
0042
0043 Acts::GeometryContext gctx;
0044 Acts::MagneticFieldContext mctx;
0045
0046 using namespace Acts::UnitLiterals;
0047
0048 using Covariance = BoundSquareMatrix;
0049 using Jacobian = BoundMatrix;
0050
0051 /// These tests do not test for a correct covariance transport but only for the
0052 /// correct conservation or modification of certain variables. A test suite for
0053 /// the numerical correctness is performed in the integration tests.
0054 BOOST_AUTO_TEST_CASE(covariance_engine_test) {
0055   // Create a test context
0056   GeometryContext tgContext = GeometryContext();
0057
0058   auto particleHypothesis = ParticleHypothesis::pion();
0059
0060   // Build a start vector
0061   Vector3 position{1., 2., 3.};
0062   double time = 4.;
0063   Vector3 direction{sqrt(5. / 22.), 3. * sqrt(2. / 55.), 7. / sqrt(110.)};
0064   double qop = 0.125;
0065   FreeVector parameters, startParameters;
0066   parameters << position[0], position[1], position[2], time, direction[0],
0067       direction[1], direction[2], qop;
0068   startParameters = parameters;
0069
0070   // Build covariance matrix, jacobians and related components
0071   Covariance covariance = Covariance::Identity();
0072   Jacobian jacobian = 2. * Jacobian::Identity();
0073   FreeMatrix transportJacobian = 3. * FreeMatrix::Identity();
0074   FreeVector derivatives;
0075   derivatives << 9., 10., 11., 12., 13., 14., 15., 16.;
0076   BoundToFreeMatrix boundToFreeJacobian = 4. * BoundToFreeMatrix::Identity();
0077
0078   // Covariance transport to curvilinear coordinates
0079   detail::transportCovarianceToCurvilinear(covariance, jacobian,
0080                                            transportJacobian, derivatives,
0081                                            boundToFreeJacobian, direction);
0082
0083   // Tests to see that the right components are (un-)changed
0084   BOOST_CHECK_NE(covariance, Covariance::Identity());
0085   BOOST_CHECK_NE(jacobian, 2. * Jacobian::Identity());
0086   BOOST_CHECK_EQUAL(transportJacobian, FreeMatrix::Identity());
0087   BOOST_CHECK_EQUAL(derivatives, FreeVector::Zero());
0088   BOOST_CHECK_NE(boundToFreeJacobian, 4. * BoundToFreeMatrix::Identity());
0089   BOOST_CHECK_EQUAL(
0090       direction, Vector3(sqrt(5. / 22.), 3. * sqrt(2. / 55.), 7. / sqrt(110.)));
0091
0092   // Reset
0093   covariance = Covariance::Identity();
0094   jacobian = 2. * Jacobian::Identity();
0095   transportJacobian = 3. * FreeMatrix::Identity();
0096   derivatives << 9., 10., 11., 12., 13., 14., 15., 16.;
0097   boundToFreeJacobian = 4. * BoundToFreeMatrix::Identity();
0098
0099   // Repeat transport to surface
0100   FreeToBoundCorrection freeToBoundCorrection(false);
0101   auto surface = CurvilinearSurface(position, direction).planeSurface();
0102   detail::transportCovarianceToBound(
0103       tgContext, *surface, covariance, jacobian, transportJacobian, derivatives,
0104       boundToFreeJacobian, parameters, freeToBoundCorrection);
0105
0106   BOOST_CHECK_NE(covariance, Covariance::Identity());
0107   BOOST_CHECK_NE(jacobian, 2. * Jacobian::Identity());
0108   BOOST_CHECK_EQUAL(transportJacobian, FreeMatrix::Identity());
0109   BOOST_CHECK_EQUAL(derivatives, FreeVector::Zero());
0110   BOOST_CHECK_NE(boundToFreeJacobian, 4. * BoundToFreeMatrix::Identity());
0111   BOOST_CHECK_EQUAL(parameters, startParameters);
0112
0113   // Produce a curvilinear state without covariance matrix
0114   auto curvResult = detail::curvilinearState(
0115       covariance, jacobian, transportJacobian, derivatives, boundToFreeJacobian,
0116       parameters, particleHypothesis, false, 1337.);
0117   BOOST_CHECK(!std::get<0>(curvResult).covariance().has_value());
0118   BOOST_CHECK_EQUAL(std::get<2>(curvResult), 1337.);
0119
0120   // Reset
0121   covariance = Covariance::Identity();
0122   jacobian = 2. * Jacobian::Identity();
0123   transportJacobian = 3. * FreeMatrix::Identity();
0124   derivatives << 9., 10., 11., 12., 13., 14., 15., 16.;
0125   boundToFreeJacobian = 4. * BoundToFreeMatrix::Identity();
0126
0127   // Produce a curvilinear state with covariance matrix
0128   curvResult = detail::curvilinearState(
0129       covariance, jacobian, transportJacobian, derivatives, boundToFreeJacobian,
0130       parameters, particleHypothesis, true, 1337.);
0131   BOOST_CHECK(std::get<0>(curvResult).covariance().has_value());
0132   BOOST_CHECK_NE(*(std::get<0>(curvResult).covariance()),
0133                  Covariance::Identity());
0134   BOOST_CHECK_NE(std::get<1>(curvResult), 2. * Jacobian::Identity());
0135   BOOST_CHECK_EQUAL(std::get<2>(curvResult), 1337.);
0136
0137   // Produce a bound state without covariance matrix
0138   auto covarianceBefore = covariance;
0139   auto boundResult =
0140       detail::boundState(tgContext, *surface, covariance, jacobian,
0141                          transportJacobian, derivatives, boundToFreeJacobian,
0142                          parameters, particleHypothesis, false, 1337.,
0143                          freeToBoundCorrection)
0144           .value();
0145   BOOST_CHECK(std::get<0>(curvResult).covariance().has_value());
0146   BOOST_CHECK_EQUAL(*(std::get<0>(curvResult).covariance()), covarianceBefore);
0147   BOOST_CHECK_EQUAL(std::get<2>(boundResult), 1337.);
0148
0149   // Reset
0150   covariance = Covariance::Identity();
0151   jacobian = 2. * Jacobian::Identity();
0152   transportJacobian = 3. * FreeMatrix::Identity();
0153   derivatives << 9., 10., 11., 12., 13., 14., 15., 16.;
0154   boundToFreeJacobian = 4. * BoundToFreeMatrix::Identity();
0155
0156   // Produce a bound state with covariance matrix
0157   boundResult =
0158       detail::boundState(tgContext, *surface, covariance, jacobian,
0159                          transportJacobian, derivatives, boundToFreeJacobian,
0160                          parameters, ParticleHypothesis::pion(), true, 1337.,
0161                          freeToBoundCorrection)
0162           .value();
0163   BOOST_CHECK(std::get<0>(boundResult).covariance().has_value());
0164   BOOST_CHECK_NE(*(std::get<0>(boundResult).covariance()),
0165                  Covariance::Identity());
0166   BOOST_CHECK_NE(std::get<1>(boundResult), 2. * Jacobian::Identity());
0167   BOOST_CHECK_EQUAL(std::get<2>(boundResult), 1337.);
0168
0169   // Reset
0170   freeToBoundCorrection.apply = true;
0171
0172   // Produce a bound state with free to bound correction
0173   boundResult =
0174       detail::boundState(tgContext, *surface, covariance, jacobian,
0175                          transportJacobian, derivatives, boundToFreeJacobian,
0176                          parameters, ParticleHypothesis::pion(), true, 1337.,
0177                          freeToBoundCorrection)
0178           .value();
0179   BOOST_CHECK(std::get<0>(boundResult).covariance().has_value());
0180   BOOST_CHECK_NE(*(std::get<0>(boundResult).covariance()),
0181                  Covariance::Identity());
0182 }
0183
0184 std::pair<BoundVector, BoundMatrix> boundToBound(const BoundVector& parIn,
0185                                                  const BoundMatrix& covIn,
0186                                                  const Surface& srfA,
0187                                                  const Surface& srfB,
0188                                                  const Vector3& bField) {
0189   Acts::BoundTrackParameters boundParamIn{srfA.getSharedPtr(), parIn, covIn,
0190                                           ParticleHypothesis::pion()};
0191
0192   auto converted =
0193       detail::boundToBoundConversion(gctx, boundParamIn, srfB, bField).value();
0194
0195   return {converted.parameters(), converted.covariance().value()};
0196 }
0197
0198 using propagator_t = Propagator<EigenStepper<>, VoidNavigator>;
0199
0200 BoundVector localToLocal(const propagator_t& prop, const BoundVector& local,
0201                          const Surface& src, const Surface& dst) {
0202   using PropagatorOptions = typename propagator_t::template Options<>;
0203   PropagatorOptions options{gctx, mctx};
0204   options.stepping.stepTolerance = 1e-10;
0205   options.surfaceTolerance = 1e-10;
0206
0207   BoundTrackParameters start{src.getSharedPtr(), local, std::nullopt,
0208                              ParticleHypothesis::pion()};
0209
0210   auto res = prop.propagate(start, dst, options).value();
0211   auto endParameters = res.endParameters.value();
0212
0213   BOOST_CHECK_EQUAL(&endParameters.referenceSurface(), &dst);
0214
0215   BoundVector out = endParameters.parameters();
0216   out[eBoundTime] = local[eBoundTime];
0217   return out;
0218 }
0219
0220 propagator_t makePropagator(const Vector3& bField) {
0221   return propagator_t{EigenStepper<>{std::make_shared<ConstantBField>(bField)},
0222                       VoidNavigator{}};
0223 }
0224
0225 BoundMatrix numericalBoundToBoundJacobian(const propagator_t& prop,
0226                                           const BoundVector& parA,
0227                                           const Surface& srfA,
0228                                           const Surface& srfB) {
0229   double h = 1e-4;
0230   BoundMatrix J;
0231   for (std::size_t i = 0; i < 6; i++) {
0232     for (std::size_t j = 0; j < 6; j++) {
0233       BoundVector parInitial1 = parA;
0234       BoundVector parInitial2 = parA;
0235       parInitial1[i] -= h;
0236       parInitial2[i] += h;
0237       BoundVector parFinal1 = localToLocal(prop, parInitial1, srfA, srfB);
0238       BoundVector parFinal2 = localToLocal(prop, parInitial2, srfA, srfB);
0239
0240       J(j, i) = (parFinal2[j] - parFinal1[j]) / (2 * h);
0241     }
0242   }
0243   return J;
0244 }
0245
0246 unsigned int getNextSeed() {
0247   static unsigned int seed = 10;
0248   return ++seed;
0249 }
0250
0251 auto makeDist(double a, double b) {
0252   return bdata::random(
0253       (bdata::engine = std::mt19937{}, bdata::seed = getNextSeed(),
0254        bdata::distribution = std::uniform_real_distribution<double>(a, b)));
0255 }
0256
0257 const auto locDist = makeDist(-5_mm, 5_mm);
0258 const auto bFieldDist = makeDist(0, 3_T);
0259 const auto angleDist = makeDist(-2 * std::numbers::pi, 2 * std::numbers::pi);
0260 const auto posDist = makeDist(-50_mm, 50_mm);
0261
0262 #define MAKE_SURFACE()                                    \
0263   [&]() {                                                 \
0264     Transform3 transformA = Transform3::Identity();       \
0265     transformA = AngleAxis3(Rx, Vector3::UnitX());        \
0266     transformA = AngleAxis3(Ry, Vector3::UnitY());        \
0267     transformA = AngleAxis3(Rz, Vector3::UnitZ());        \
0268     transformA.translation() << gx, gy, gz;               \
0269     return Surface::makeShared<PlaneSurface>(transformA); \
0270   }()
0271
0272 BOOST_DATA_TEST_CASE(CovarianceConversionSamePlane,
0273                      (bFieldDist ^ bFieldDist ^ bFieldDist ^ angleDist ^
0274                       angleDist ^ angleDist ^ posDist ^ posDist ^ posDist ^
0275                       locDist ^ locDist) ^
0276                          bdata::xrange(100),
0277                      Bx, By, Bz, Rx, Ry, Rz, gx, gy, gz, l0, l1, index) {
0278   (void)index;
0279   const Vector3 bField{Bx, By, Bz};
0280
0281   auto planeSurfaceA = MAKE_SURFACE();
0282   auto planeSurfaceB =
0283       Surface::makeShared<PlaneSurface>(planeSurfaceA->transform(gctx));
0284
0285   BoundMatrix covA;
0286   covA.setZero();
0287   covA.diagonal() << 1, 2, 3, 4, 5, 6;
0288
0289   BoundVector parA;
0290   parA << l0, l1, std::numbers::pi / 4., std::numbers::pi / 2. * 0.9,
0291       -1 / 1_GeV, 5_ns;
0292
0293   // identical surface, this should work
0294   auto [parB, covB] =
0295       boundToBound(parA, covA, *planeSurfaceA, *planeSurfaceB, bField);
0296
0297   // these should be the same because the plane surface are the same
0298   CHECK_CLOSE_ABS(parA, parB, 1e-9);
0299   CHECK_CLOSE_COVARIANCE(covA, covB, 1e-9);
0300
0301   // now go back
0302   auto [parA2, covA2] =
0303       boundToBound(parB, covB, *planeSurfaceB, *planeSurfaceA, bField);
0304   CHECK_CLOSE_ABS(parA, parA2, 1e-9);
0305   CHECK_CLOSE_COVARIANCE(covA, covA2, 1e-9);
0306
0307   auto prop = makePropagator(bField);
0308
0309   BoundMatrix J =
0310       numericalBoundToBoundJacobian(prop, parA, *planeSurfaceA, *planeSurfaceB);
0311   BoundMatrix covC = J * covA * J.transpose();
0312   CHECK_CLOSE_COVARIANCE(covB, covC, 1e-6);
0313 }
0314
0315 BOOST_DATA_TEST_CASE(CovarianceConversionRotatedPlane,
0316                      (bFieldDist ^ bFieldDist ^ bFieldDist ^ angleDist ^
0317                       angleDist ^ angleDist ^ posDist ^ posDist ^ posDist ^
0318                       locDist ^ locDist ^ angleDist) ^
0319                          bdata::xrange(100),
0320                      Bx, By, Bz, Rx, Ry, Rz, gx, gy, gz, l0, l1, angle, index) {
0321   (void)index;
0322   const Vector3 bField{Bx, By, Bz};
0323
0324   auto planeSurfaceA = MAKE_SURFACE();
0325
0326   Transform3 transform;
0327   transform = planeSurfaceA->transform(gctx).rotation();
0328   transform = AngleAxis3(angle, planeSurfaceA->normal(gctx)) * transform;
0329   transform.translation() = planeSurfaceA->transform(gctx).translation();
0330   auto planeSurfaceB = Surface::makeShared<PlaneSurface>(transform);
0331
0332   // sanity check that the normal didn't change
0333   CHECK_CLOSE_ABS(planeSurfaceA->normal(gctx), planeSurfaceB->normal(gctx),
0334                   1e-9);
0335
0336   BoundMatrix covA;
0337   covA.setZero();
0338   covA.diagonal() << 1, 2, 3, 4, 5, 6;
0339
0340   BoundVector parA;
0341   parA << l0, l1, std::numbers::pi / 4., std::numbers::pi / 2. * 0.9,
0342       -1 / 1_GeV, 5_ns;
0343
0344   auto [parB, covB] =
0345       boundToBound(parA, covA, *planeSurfaceA, *planeSurfaceB, bField);
0346   BoundVector exp = parA;
0347   // loc0 and loc1 are rotated
0348   exp.head<2>() = Eigen::Rotation2D<double>(-angle) * parA.head<2>();
0349
0350   CHECK_CLOSE_ABS(exp, parB, 1e-9);
0351
0352   // now go back
0353   auto [parA2, covA2] =
0354       boundToBound(parB, covB, *planeSurfaceB, *planeSurfaceA, bField);
0355   CHECK_CLOSE_ABS(parA, parA2, 1e-9);
0356   CHECK_CLOSE_COVARIANCE(covA, covA2, 1e-9);
0357
0358   auto prop = makePropagator(bField);
0359   BoundMatrix J =
0360       numericalBoundToBoundJacobian(prop, parA, *planeSurfaceA, *planeSurfaceB);
0361   BoundMatrix covC = J * covA * J.transpose();
0362   CHECK_CLOSE_COVARIANCE(covB, covC, 1e-6);
0363 }
0364
0365 BOOST_DATA_TEST_CASE(CovarianceConversionL0TiltedPlane,
0366                      (bFieldDist ^ bFieldDist ^ bFieldDist ^ angleDist ^
0367                       angleDist ^ angleDist ^ posDist ^ posDist ^ posDist ^
0368                       locDist ^ angleDist) ^
0369                          bdata::xrange(100),
0370                      Bx, By, Bz, Rx, Ry, Rz, gx, gy, gz, l1, angle, index) {
0371   (void)index;
0372   const Vector3 bField{Bx, By, Bz};
0373
0374   auto planeSurfaceA = MAKE_SURFACE();
0375
0376   // make plane that is slightly rotated
0377   Transform3 transform;
0378   transform = planeSurfaceA->transform(gctx).rotation();
0379
0380   // figure out rotation axis along local x
0381   Vector3 axis = planeSurfaceA->transform(gctx).rotation() * Vector3::UnitY();
0382   transform = AngleAxis3(angle, axis) * transform;
0383
0384   transform.translation() = planeSurfaceA->transform(gctx).translation();
0385
0386   auto planeSurfaceB = Surface::makeShared<PlaneSurface>(transform);
0387
0388   BoundVector parA;
0389   // loc 0 must be zero so we're on the intersection of both surfaces.
0390   parA << 0, l1, std::numbers::pi / 4., std::numbers::pi / 2. * 0.9, -1 / 1_GeV,
0391       5_ns;
0392
0393   BoundMatrix covA;
0394   covA.setZero();
0395   covA.diagonal() << 1, 2, 3, 4, 5, 6;
0396
0397   auto [parB, covB] =
0398       boundToBound(parA, covA, *planeSurfaceA, *planeSurfaceB, bField);
0399
0400   // now go back
0401   auto [parA2, covA2] =
0402       boundToBound(parB, covB, *planeSurfaceB, *planeSurfaceA, bField);
0403   CHECK_CLOSE_ABS(parA, parA2, 1e-9);
0404   CHECK_CLOSE_COVARIANCE(covA, covA2, 1e-7);
0405
0406   auto prop = makePropagator(bField);
0407   BoundMatrix J =
0408       numericalBoundToBoundJacobian(prop, parA, *planeSurfaceA, *planeSurfaceB);
0409   BoundMatrix covC = J * covA * J.transpose();
0410   CHECK_CLOSE_OR_SMALL((covB.template topLeftCorner<2, 2>()),
0411                        (covC.template topLeftCorner<2, 2>()), 1e-7, 1e-9);
0412   CHECK_CLOSE_OR_SMALL(covB.diagonal(), covC.diagonal(), 1e-7, 1e-9);
0413 }
0414
0415 BOOST_DATA_TEST_CASE(CovarianceConversionL1TiltedPlane,
0416                      (bFieldDist ^ bFieldDist ^ bFieldDist ^ angleDist ^
0417                       angleDist ^ angleDist ^ posDist ^ posDist ^ posDist ^
0418                       locDist ^ angleDist) ^
0419                          bdata::xrange(100),
0420                      Bx, By, Bz, Rx, Ry, Rz, gx, gy, gz, l0, angle, index) {
0421   (void)index;
0422   const Vector3 bField{Bx, By, Bz};
0423
0424   auto planeSurfaceA = MAKE_SURFACE();
0425
0426   // make plane that is slightly rotated
0427   Transform3 transform;
0428   transform = planeSurfaceA->transform(gctx).rotation();
0429
0430   Vector3 axis = planeSurfaceA->transform(gctx).rotation() * Vector3::UnitX();
0431   transform = AngleAxis3(angle, axis) * transform;
0432
0433   transform.translation() = planeSurfaceA->transform(gctx).translation();
0434
0435   auto planeSurfaceB = Surface::makeShared<PlaneSurface>(transform);
0436
0437   BoundVector parA;
0438   // loc 1 must be zero so we're on the intersection of both surfaces.
0439   parA << l0, 0, std::numbers::pi / 4., std::numbers::pi / 2. * 0.9, -1 / 1_GeV,
0440       5_ns;
0441
0442   BoundMatrix covA;
0443   covA.setZero();
0444   covA.diagonal() << 1, 2, 3, 4, 5, 6;
0445
0446   auto [parB, covB] =
0447       boundToBound(parA, covA, *planeSurfaceA, *planeSurfaceB, bField);
0448
0449   // now go back
0450   auto [parA2, covA2] =
0451       boundToBound(parB, covB, *planeSurfaceB, *planeSurfaceA, bField);
0452   CHECK_CLOSE_ABS(parA, parA2, 1e-9);
0453   // tolerance is a bit higher here
0454   CHECK_CLOSE_COVARIANCE(covA, covA2, 1e-6);
0455
0456   auto prop = makePropagator(bField);
0457   BoundMatrix J =
0458       numericalBoundToBoundJacobian(prop, parA, *planeSurfaceA, *planeSurfaceB);
0459   BoundMatrix covC = J * covA * J.transpose();
0460   CHECK_CLOSE_OR_SMALL((covB.template topLeftCorner<2, 2>()),
0461                        (covC.template topLeftCorner<2, 2>()), 1e-6, 1e-9);
0462   CHECK_CLOSE_OR_SMALL(covB.diagonal(), covC.diagonal(), 1e-6, 1e-9);
0463 }
0464
0465 BOOST_DATA_TEST_CASE(CovarianceConversionPerigee,
0466                      (bFieldDist ^ bFieldDist ^ bFieldDist ^ angleDist ^
0467                       angleDist ^ angleDist ^ posDist ^ posDist ^ posDist ^
0468                       locDist ^ locDist ^ angleDist ^ angleDist ^ angleDist) ^
0469                          bdata::xrange(100),
0470                      Bx, By, Bz, Rx, Ry, Rz, gx, gy, gz, l0, l1, pRx, pRy, pRz,
0471                      index) {
0472   (void)index;
0473   const Vector3 bField{Bx, By, Bz};
0474
0475   auto planeSurfaceA = MAKE_SURFACE();
0476
0477   BoundVector parA;
0478   parA << l0, l1, std::numbers::pi / 4., std::numbers::pi / 2. * 0.9,
0479       -1 / 1_GeV, 5_ns;
0480
0481   BoundMatrix covA;
0482   covA.setZero();
0483   covA.diagonal() << 1, 2, 3, 4, 5, 6;
0484
0485   Vector3 global = planeSurfaceA->localToGlobal(gctx, parA.head<2>());
0486
0487   Transform3 transform;
0488   transform.setIdentity();
0489   transform.rotate(AngleAxis3(pRx, Vector3::UnitX()));
0490   transform.rotate(AngleAxis3(pRy, Vector3::UnitY()));
0491   transform.rotate(AngleAxis3(pRz, Vector3::UnitZ()));
0492   transform.translation() = global;
0493
0494   auto perigee = Surface::makeShared<PerigeeSurface>(transform);
0495
0496   auto [parB, covB] =
0497       boundToBound(parA, covA, *planeSurfaceA, *perigee, bField);
0498
0499   // now go back
0500   auto [parA2, covA2] =
0501       boundToBound(parB, covB, *perigee, *planeSurfaceA, bField);
0502   CHECK_CLOSE_ABS(parA, parA2, 1e-9);
0503   CHECK_CLOSE_COVARIANCE(covA, covA2, 1e-9);
0504
0505   auto prop = makePropagator(bField);
0506   BoundMatrix J =
0507       numericalBoundToBoundJacobian(prop, parA, *planeSurfaceA, *perigee);
0508   BoundMatrix covC = J * covA * J.transpose();
0509   CHECK_CLOSE_ABS(covB, covC, 1e-7);
0510 }
0511
0512 }  // namespace Acts::Test