Core/Propagator/JacobianEngineTests.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/unit_test.hpp>
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Definitions/TrackParametrization.hpp"
0013 #include "Acts/Geometry/GeometryContext.hpp"
0014 #include "Acts/Propagator/detail/JacobianEngine.hpp"
0015 #include "Acts/Surfaces/CurvilinearSurface.hpp"
0016 #include "Acts/Surfaces/PlaneSurface.hpp"
0017 #include "Acts/Surfaces/Surface.hpp"
0018
0019 namespace Acts::Test {
0020
0021 /// These tests do not test for a correct covariance transport but only for the
0022 /// correct conservation or modification of certain variables. A test suite for
0023 /// the numerical correctness is performed in the integration tests.
0024 ///
0025 BOOST_AUTO_TEST_CASE(jacobian_engine_to_bound) {
0026   // Create a test context
0027   GeometryContext tgContext = GeometryContext();
0028
0029   // Build a start vector
0030   Vector3 position{1., 2., 3.};
0031   double time = 4.;
0032   Vector3 direction = Vector3(5., 2., 7.).normalized();
0033   double qop = 0.125;
0034
0035   // Build a surface
0036   std::shared_ptr<PlaneSurface> pSurface =
0037       CurvilinearSurface(position, direction).planeSurface();
0038
0039   // Other rotated surface
0040   Vector3 odirection = Vector3(6., 2., 8.).normalized();
0041   std::shared_ptr<PlaneSurface> oSurface =
0042       CurvilinearSurface(position, odirection).planeSurface();
0043
0044   // The free parameter vector
0045   FreeVector freeParameters;
0046   freeParameters << position[0], position[1], position[2], time, direction[0],
0047       direction[1], direction[2], qop;
0048   // And its associated bound vector
0049   BoundVector boundParameters;
0050   boundParameters << 0., 0., VectorHelpers::phi(direction),
0051       VectorHelpers::theta(direction), qop, time;
0052
0053   // Build covariance matrices for bound and free case
0054   BoundSquareMatrix boundCovariance = 0.025 * BoundSquareMatrix::Identity();
0055   FreeSquareMatrix freeCovariance = 0.025 * FreeSquareMatrix::Identity();
0056
0057   FreeMatrix noTransportJacobian = FreeMatrix::Identity();
0058   FreeMatrix realTransportJacobian = 2 * FreeMatrix::Identity();
0059
0060   FreeToPathMatrix freeToPathDerivatives =
0061       pSurface->freeToPathDerivative(tgContext, position, direction);
0062   BoundToFreeMatrix boundToFreeJacobian =
0063       pSurface->boundToFreeJacobian(tgContext, position, direction);
0064
0065   // (1) curvilinear/bound to bound transport jacobian
0066   // a) test without actual transport into the same surface
0067   BoundMatrix b2bTransportJacobian;
0068   detail::boundToBoundTransportJacobian(
0069       tgContext, *pSurface, freeParameters, boundToFreeJacobian,
0070       noTransportJacobian, freeToPathDerivatives, b2bTransportJacobian);
0071   BoundSquareMatrix newBoundCovariance =
0072       b2bTransportJacobian * boundCovariance * b2bTransportJacobian.transpose();
0073   BOOST_CHECK(boundCovariance.isApprox(newBoundCovariance));
0074   // b) test without actual transport but to a new surface
0075   detail::boundToBoundTransportJacobian(
0076       tgContext, *oSurface, freeParameters, boundToFreeJacobian,
0077       noTransportJacobian, freeToPathDerivatives, b2bTransportJacobian);
0078   newBoundCovariance =
0079       b2bTransportJacobian * boundCovariance * b2bTransportJacobian.transpose();
0080   BOOST_CHECK(!boundCovariance.isApprox(newBoundCovariance));
0081   // c) test to "the same" surface with transport
0082   // (not really senseful, but should give a different result)
0083   detail::boundToBoundTransportJacobian(
0084       tgContext, *pSurface, freeParameters, boundToFreeJacobian,
0085       realTransportJacobian, freeToPathDerivatives, b2bTransportJacobian);
0086   newBoundCovariance =
0087       b2bTransportJacobian * boundCovariance * b2bTransportJacobian.transpose();
0088   BOOST_CHECK(!boundCovariance.isApprox(newBoundCovariance));
0089
0090   FreeToBoundMatrix f2bTransportJacobian;
0091   detail::freeToBoundTransportJacobian(
0092       tgContext, *pSurface, freeParameters, noTransportJacobian,
0093       freeToPathDerivatives, f2bTransportJacobian);
0094
0095   newBoundCovariance =
0096       f2bTransportJacobian * freeCovariance * f2bTransportJacobian.transpose();
0097   BOOST_CHECK(!boundCovariance.isApprox(newBoundCovariance));
0098 }
0099
0100 /// These tests do not test for a correct covariance transport but only for the
0101 /// correct conservation or modification of certain variables. A test suite for
0102 /// the numerical correctness is performed in the integration tests.
0103 ///
0104 BOOST_AUTO_TEST_CASE(jacobian_engine_to_curvilinear) {
0105   // Create a test context
0106   GeometryContext tgContext = GeometryContext();
0107
0108   // Build a start vector
0109   Vector3 position{1., 2., 3.};
0110   Vector3 direction = Vector3(5., 2., 7.).normalized();
0111
0112   // Build a surface, starting surface for curvilinear
0113   std::shared_ptr<PlaneSurface> pSurface =
0114       CurvilinearSurface(position, direction).planeSurface();
0115
0116   // Build covariance matrices for bound and free case
0117   BoundSquareMatrix boundCovariance = 0.025 * BoundSquareMatrix::Identity();
0118   FreeSquareMatrix freeCovariance = 0.025 * FreeSquareMatrix::Identity();
0119
0120   FreeMatrix noTransportJacobian = FreeMatrix::Identity();
0121
0122   FreeToPathMatrix freeToPathDerivatives =
0123       pSurface->freeToPathDerivative(tgContext, position, direction);
0124   BoundToFreeMatrix boundToFreeJacobian =
0125       CurvilinearSurface(direction).boundToFreeJacobian();
0126
0127   // (1) curvilinear/bound to curvilinear transport jacobian
0128   // a) test without actual transport into the same surface
0129   BoundMatrix b2cTransportJacobian;
0130   detail::boundToCurvilinearTransportJacobian(
0131       direction, boundToFreeJacobian, noTransportJacobian,
0132       freeToPathDerivatives, b2cTransportJacobian);
0133   BoundSquareMatrix newBoundCovariance =
0134       b2cTransportJacobian * boundCovariance * b2cTransportJacobian.transpose();
0135   BOOST_CHECK(boundCovariance.isApprox(newBoundCovariance));
0136   // b) test to another curvilinear frame at the same point (no transport)
0137   detail::boundToCurvilinearTransportJacobian(
0138       Vector3(4., 5., 6.).normalized(), boundToFreeJacobian,
0139       noTransportJacobian, freeToPathDerivatives, b2cTransportJacobian);
0140   newBoundCovariance =
0141       b2cTransportJacobian * boundCovariance * b2cTransportJacobian.transpose();
0142   BOOST_CHECK(!boundCovariance.isApprox(newBoundCovariance));
0143
0144   FreeToBoundMatrix f2cTransportJacobian =
0145       detail::freeToCurvilinearTransportJacobian(direction, noTransportJacobian,
0146                                                  freeToPathDerivatives);
0147
0148   newBoundCovariance =
0149       f2cTransportJacobian * freeCovariance * f2cTransportJacobian.transpose();
0150   BOOST_CHECK(!boundCovariance.isApprox(newBoundCovariance));
0151 }
0152
0153 /// These tests do not test for a correct covariance transport but only for the
0154 /// correct conservation or modification of certain variables. A test suite for
0155 /// the numerical correctness is performed in the integration tests.
0156 ///
0157 BOOST_AUTO_TEST_CASE(jacobian_engine_to_free) {
0158   // Create a test context
0159   GeometryContext tgContext = GeometryContext();
0160
0161   // Build a start vector
0162   Vector3 position{1., 2., 3.};
0163   Vector3 direction = Vector3(5., 2., 7.).normalized();
0164
0165   // Build a surface, starting surface for curvilinear
0166   std::shared_ptr<PlaneSurface> pSurface =
0167       CurvilinearSurface(position, direction).planeSurface();
0168
0169   // Build covariance matrices for bound and free case
0170   BoundSquareMatrix boundCovariance = 0.025 * BoundSquareMatrix::Identity();
0171   FreeSquareMatrix freeCovariance = 0.025 * FreeSquareMatrix::Identity();
0172
0173   FreeMatrix noTransportJacobian = FreeMatrix::Identity();
0174
0175   BoundToFreeMatrix boundToFreeJacobian =
0176       pSurface->boundToFreeJacobian(tgContext, position, direction);
0177
0178   // (1) bound to free
0179   BoundToFreeMatrix b2fTransportJacobian = detail::boundToFreeTransportJacobian(
0180       boundToFreeJacobian, noTransportJacobian);
0181
0182   FreeMatrix newFreeCovariance1 =
0183       b2fTransportJacobian * boundCovariance * b2fTransportJacobian.transpose();
0184   BOOST_CHECK(!newFreeCovariance1.isApprox(freeCovariance));
0185
0186   // (2) curvilinear to free
0187   boundToFreeJacobian = CurvilinearSurface(direction).boundToFreeJacobian();
0188   BoundToFreeMatrix c2fTransportJacobian = detail::boundToFreeTransportJacobian(
0189       boundToFreeJacobian, noTransportJacobian);
0190
0191   FreeMatrix newFreeCovariance2 =
0192       c2fTransportJacobian * boundCovariance * c2fTransportJacobian.transpose();
0193   BOOST_CHECK(!newFreeCovariance2.isApprox(freeCovariance));
0194   // But those should be similar/equal
0195   BOOST_CHECK(newFreeCovariance1.isApprox(newFreeCovariance2));
0196 }
0197
0198 }  // namespace Acts::Test