Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2026-04-08 07:48:07 

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

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!