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