Core/Propagator/JacobianTests.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/Definitions/Units.hpp"
0014 #include "Acts/EventData/TrackParameters.hpp"
0015 #include "Acts/Geometry/GeometryContext.hpp"
0016 #include "Acts/MagneticField/ConstantBField.hpp"
0017 #include "Acts/MagneticField/MagneticFieldContext.hpp"
0018 #include "Acts/Propagator/AtlasStepper.hpp"
0019 #include "Acts/Propagator/EigenStepper.hpp"
0020 #include "Acts/Surfaces/CurvilinearSurface.hpp"
0021 #include "Acts/Surfaces/CylinderSurface.hpp"
0022 #include "Acts/Surfaces/DiscSurface.hpp"
0023 #include "Acts/Surfaces/PerigeeSurface.hpp"
0024 #include "Acts/Surfaces/PlaneSurface.hpp"
0025 #include "Acts/Surfaces/StrawSurface.hpp"
0026 #include "Acts/Surfaces/Surface.hpp"
0027 #include "ActsTests/CommonHelpers/FloatComparisons.hpp"
0028
0029 #include <array>
0030 #include <cstddef>
0031 #include <memory>
0032 #include <optional>
0033 #include <utility>
0034
0035 using namespace Acts;
0036 using namespace Acts::UnitLiterals;
0037
0038 namespace ActsTests {
0039
0040 using BFieldType = ConstantBField;
0041 using EigenStepperType = EigenStepper<>;
0042 using AtlasStepperType = AtlasStepper;
0043 using Covariance = BoundSquareMatrix;
0044
0045 // Create a test context
0046 GeometryContext tgContext = GeometryContext();
0047 MagneticFieldContext mfContext = MagneticFieldContext();
0048
0049 static auto bField = std::make_shared<BFieldType>(Vector3{0, 0, 1_T});
0050
0051 /// Helper method to create a transform for a plane
0052 /// to mimic detector situations, the plane is roughly
0053 /// perpendicular to the track
0054 ///
0055 /// @param nnomal The nominal normal direction
0056 /// @param angleT Rotation around the norminal normal
0057 /// @param angleU Rotation around the original U axis
0058 Transform3 createCylindricTransform(const Vector3& nposition, double angleX,
0059                                     double angleY) {
0060   Transform3 ctransform;
0061   ctransform.setIdentity();
0062   ctransform.pretranslate(nposition);
0063   ctransform.prerotate(AngleAxis3(angleX, Vector3::UnitX()));
0064   ctransform.prerotate(AngleAxis3(angleY, Vector3::UnitY()));
0065   return ctransform;
0066 }
0067
0068 /// Helper method to create a transform for a plane
0069 /// to mimic detector situations, the plane is roughly
0070 /// perpendicular to the track
0071 ///
0072 /// @param nnomal The nominal normal direction
0073 /// @param angleT Rotation around the norminal normal
0074 /// @param angleU Rotation around the original U axis
0075 Transform3 createPlanarTransform(const Vector3& nposition,
0076                                  const Vector3& nnormal, double angleT,
0077                                  double angleU) {
0078   // the rotation of the destination surface
0079   Vector3 T = nnormal.normalized();
0080   Vector3 U = std::abs(T.dot(Vector3::UnitZ())) < 0.99
0081                   ? Vector3::UnitZ().cross(T).normalized()
0082                   : Vector3::UnitX().cross(T).normalized();
0083   Vector3 V = T.cross(U);
0084   // that's the plane curvilinear Rotation
0085   RotationMatrix3 curvilinearRotation;
0086   curvilinearRotation.col(0) = U;
0087   curvilinearRotation.col(1) = V;
0088   curvilinearRotation.col(2) = T;
0089   // curvilinear surfaces are boundless
0090   Transform3 ctransform{curvilinearRotation};
0091   ctransform.pretranslate(nposition);
0092   ctransform.prerotate(AngleAxis3(angleT, T));
0093   ctransform.prerotate(AngleAxis3(angleU, U));
0094   //
0095   return ctransform;
0096 }
0097
0098 /// Helper method : convert into Acts matrix
0099 /// It takes the double array from AtlasStepper
0100 /// and transforms it into an ActsMatrixD
0101 ///
0102 /// @param P is the pointer to the array
0103 ///
0104 /// Translation is (for lookup)
0105 ///                   /dL0    /dL1    /dPhi   /dThe   /dCM   /dT
0106 /// X  ->P[0]  dX /   P[ 8]   P[16]   P[24]   P[32]   P[40]  P[48]
0107 /// Y  ->P[1]  dY /   P[ 9]   P[17]   P[25]   P[33]   P[41]  P[49]
0108 /// Z  ->P[2]  dZ /   P[10]   P[18]   P[26]   P[34]   P[42]  P[50]
0109 /// T  ->P[3]  dT/    P[11]   P[19]   P[27]   P[35]   P[43]  P[51]
0110 /// Ax ->P[4]  dAx/   P[12]   P[20]   P[28]   P[36]   P[44]  P[52]
0111 /// Ay ->P[5]  dAy/   P[13]   P[21]   P[29]   P[37]   P[45]  P[53]
0112 /// Az ->P[6]  dAz/   P[14]   P[22]   P[30]   P[38]   P[46]  P[54]
0113 /// CM ->P[7]  dCM/   P[15]   P[23]   P[31]   P[39]   P[47]  P[55]
0114
0115 BoundToFreeMatrix convertToMatrix(const std::array<double, 60> P) {
0116   // initialize to zero
0117   BoundToFreeMatrix jMatrix = BoundToFreeMatrix::Zero();
0118   for (std::size_t j = 0; j < eBoundSize; ++j) {
0119     for (std::size_t i = 0; i < eFreeSize; ++i) {
0120       std::size_t ijc = eFreeSize + j * eFreeSize + i;
0121       jMatrix(i, j) = P[ijc];
0122     }
0123   }
0124   return jMatrix;
0125 }
0126
0127 /// Helper method : tests the jacobian to Global
0128 /// for a templated Parameters object
0129 ///
0130 /// @tparam Parameters the parameter type
0131 /// @param pars the parameter object
0132 template <typename Parameters>
0133 void testJacobianToGlobal(const Parameters& pars) {
0134   // Jacobian creation for Propagator/Steppers
0135   // a) ATLAS stepper
0136   AtlasStepperType astep(bField);
0137   AtlasStepperType::State astepState =
0138       astep.makeState(AtlasStepperType::Options(tgContext, mfContext));
0139   astep.initialize(astepState, pars);
0140   // b) Eigen stepper
0141   EigenStepperType estep(bField);
0142   EigenStepperType::State estepState =
0143       estep.makeState(EigenStepperType::Options(tgContext, mfContext));
0144   estep.initialize(estepState, pars);
0145
0146   // create the matrices
0147   auto asMatrix = convertToMatrix(astepState.pVector);
0148
0149   // cross comparison checks
0150   CHECK_CLOSE_OR_SMALL(asMatrix, estepState.jacToGlobal, 1e-6, 1e-9);
0151 }
0152
0153 BOOST_AUTO_TEST_SUITE(PropagatorSuite)
0154
0155 /// This tests the jacobian of local curvilinear -> global
0156 BOOST_AUTO_TEST_CASE(JacobianCurvilinearToGlobalTest) {
0157   // Create curvilinear parameters
0158   Covariance cov;
0159   cov << 10_mm, 0, 0, 0, 0, 0, 0, 10_mm, 0, 0, 0, 0, 0, 0, 0.1, 0, 0, 0, 0, 0,
0160       0, 0.1, 0, 0, 0, 0, 0, 0, 1. / (10_GeV), 0, 0, 0, 0, 0, 0, 0;
0161   BoundTrackParameters curvilinear = BoundTrackParameters::createCurvilinear(
0162       Vector4(341., 412., 93., 0.), Vector3(1.2, 8.3, 0.45), 1 / 10.0, cov,
0163       ParticleHypothesis::pion());
0164
0165   // run the test
0166   testJacobianToGlobal(curvilinear);
0167 }
0168
0169 /// This tests the jacobian of local cylinder -> global
0170 BOOST_AUTO_TEST_CASE(JacobianCylinderToGlobalTest) {
0171   // the cylinder transform and surface
0172   auto cTransform = createCylindricTransform({10., -5., 0.}, 0.004, 0.03);
0173   auto cSurface = Surface::makeShared<CylinderSurface>(cTransform, 200., 1000.);
0174
0175   Covariance cov;
0176   cov << 10_mm, 0, 0, 0, 0, 0, 0, 10_mm, 0, 0, 0, 0, 0, 0, 0.1, 0, 0, 0, 0, 0,
0177       0, 0.1, 0, 0, 0, 0, 0, 0, 1. / (10_GeV), 0, 0, 0, 0, 0, 0, 0;
0178
0179   BoundVector pars;
0180   pars << 182.34, -82., 0.134, 0.85, 1. / (100_GeV), 0;
0181
0182   BoundTrackParameters atCylinder(cSurface, pars, std::move(cov),
0183                                   ParticleHypothesis::pion());
0184
0185   // run the test
0186   testJacobianToGlobal(atCylinder);
0187 }
0188
0189 /// This tests the jacobian of local disc -> global
0190 BOOST_AUTO_TEST_CASE(JacobianDiscToGlobalTest) {
0191   // the disc transform and surface
0192   auto dTransform = createPlanarTransform(
0193       {10., -5., 0.}, Vector3(0.23, 0.07, 1.).normalized(), 0.004, 0.03);
0194   auto dSurface = Surface::makeShared<DiscSurface>(dTransform, 200., 1000.);
0195
0196   Covariance cov;
0197   cov << 10_mm, 0, 0, 0, 0, 0, 0, 10_mm, 0, 0, 0, 0, 0, 0, 0.1, 0, 0, 0, 0, 0,
0198       0, 0.1, 0, 0, 0, 0, 0, 0, 1. / (10_GeV), 0, 0, 0, 0, 0, 0, 0;
0199
0200   BoundVector pars;
0201   pars << 192.34, 1.823, 0.734, 0.235, 1. / (100_GeV), 0;
0202
0203   BoundTrackParameters atDisc(dSurface, pars, std::move(cov),
0204                               ParticleHypothesis::pion());
0205
0206   // run the test
0207   testJacobianToGlobal(atDisc);
0208 }
0209
0210 /// This tests the jacobian of local plane -> global
0211 BOOST_AUTO_TEST_CASE(JacobianPlaneToGlobalTest) {
0212   // Let's create a surface somewhere in space
0213   Vector3 sPosition(3421., 112., 893.);
0214   Vector3 sNormal = Vector3(1.2, -0.3, 0.05).normalized();
0215
0216   // Create a surface & parameters with covariance on the surface
0217   std::shared_ptr<PlaneSurface> pSurface =
0218       CurvilinearSurface(sPosition, sNormal).planeSurface();
0219
0220   Covariance cov;
0221   cov << 10_mm, 0, 0, 0, 0, 0, 0, 10_mm, 0, 0, 0, 0, 0, 0, 0.1, 0, 0, 0, 0, 0,
0222       0, 0.1, 0, 0, 0, 0, 0, 0, 1. / (10_GeV), 0, 0, 0, 0, 0, 0, 0;
0223
0224   BoundVector pars;
0225   pars << 12.34, -8722., 2.134, 0.85, 1. / (100_GeV), 0;
0226
0227   BoundTrackParameters atPlane(pSurface, pars, std::move(cov),
0228                                ParticleHypothesis::pion());
0229
0230   // run the test
0231   testJacobianToGlobal(atPlane);
0232 }
0233
0234 /// This tests the jacobian of local perigee -> global
0235 BOOST_AUTO_TEST_CASE(JacobianPerigeeToGlobalTest) {
0236   // Create a surface & parameters with covariance on the surface
0237   auto pSurface = Surface::makeShared<PerigeeSurface>(Vector3({0., 0., 0.}));
0238
0239   Covariance cov;
0240   cov << 10_mm, 0, 0, 0, 0, 0, 0, 10_mm, 0, 0, 0, 0, 0, 0, 0.1, 0, 0, 0, 0, 0,
0241       0, 0.1, 0, 0, 0, 0, 0, 0, 1. / (10_GeV), 0, 0, 0, 0, 0, 0, 0;
0242   BoundVector pars;
0243   pars << -3.34, -822., -0.734, 0.85, 1. / (100_GeV), 0;
0244
0245   BoundTrackParameters perigee(pSurface, pars, std::move(cov),
0246                                ParticleHypothesis::pion());
0247
0248   // run the test
0249   testJacobianToGlobal(perigee);
0250 }
0251
0252 /// This tests the jacobian of local straw -> global
0253 BOOST_AUTO_TEST_CASE(JacobianStrawToGlobalTest) {
0254   // Create a surface & parameters with covariance on the surface
0255   auto sTransform = createCylindricTransform({1019., -52., 382.}, 0.4, -0.3);
0256   auto sSurface = Surface::makeShared<StrawSurface>(sTransform, 10., 1000.);
0257
0258   Covariance cov;
0259   cov << 10_mm, 0, 0, 0, 0, 0, 0, 10_mm, 0, 0, 0, 0, 0, 0, 0.1, 0, 0, 0, 0, 0,
0260       0, 0.1, 0, 0, 0, 0, 0, 0, 1. / (10_GeV), 0, 0, 0, 0, 0, 0, 0;
0261
0262   BoundVector pars;
0263   pars << -8.34, 812., 0.734, 0.25, 1. / (100_GeV), 0;
0264
0265   BoundTrackParameters atStraw(sSurface, pars, std::move(cov),
0266                                ParticleHypothesis::pion());
0267
0268   // run the test
0269   testJacobianToGlobal(atStraw);
0270 }
0271
0272 BOOST_AUTO_TEST_SUITE_END()
0273
0274 }  // namespace ActsTests