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