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 "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
0028
0029 #include <array>
0030 #include <cstddef>
0031 #include <memory>
0032 #include <optional>
0033 #include <utility>
0034
0035 using namespace Acts::UnitLiterals;
0036
0037 namespace Acts::Test {
0038
0039 using BFieldType = ConstantBField;
0040 using EigenStepperType = EigenStepper<>;
0041 using AtlasStepperType = AtlasStepper;
0042 using Covariance = BoundSquareMatrix;
0043
0044 // Create a test context
0045 GeometryContext tgContext = GeometryContext();
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 /// This tests the jacobian of local curvilinear -> global
0153 BOOST_AUTO_TEST_CASE(JacobianCurvilinearToGlobalTest) {
0154   // Create curvilinear parameters
0155   Covariance cov;
0156   cov << 10_mm, 0, 0, 0, 0, 0, 0, 10_mm, 0, 0, 0, 0, 0, 0, 0.1, 0, 0, 0, 0, 0,
0157       0, 0.1, 0, 0, 0, 0, 0, 0, 1. / (10_GeV), 0, 0, 0, 0, 0, 0, 0;
0158   BoundTrackParameters curvilinear = BoundTrackParameters::createCurvilinear(
0159       Vector4(341., 412., 93., 0.), Vector3(1.2, 8.3, 0.45), 1 / 10.0, cov,
0160       ParticleHypothesis::pion());
0161
0162   // run the test
0163   testJacobianToGlobal(curvilinear);
0164 }
0165
0166 /// This tests the jacobian of local cylinder -> global
0167 BOOST_AUTO_TEST_CASE(JacobianCylinderToGlobalTest) {
0168   // the cylinder transform and surface
0169   auto cTransform = createCylindricTransform({10., -5., 0.}, 0.004, 0.03);
0170   auto cSurface = Surface::makeShared<CylinderSurface>(cTransform, 200., 1000.);
0171
0172   Covariance cov;
0173   cov << 10_mm, 0, 0, 0, 0, 0, 0, 10_mm, 0, 0, 0, 0, 0, 0, 0.1, 0, 0, 0, 0, 0,
0174       0, 0.1, 0, 0, 0, 0, 0, 0, 1. / (10_GeV), 0, 0, 0, 0, 0, 0, 0;
0175
0176   BoundVector pars;
0177   pars << 182.34, -82., 0.134, 0.85, 1. / (100_GeV), 0;
0178
0179   BoundTrackParameters atCylinder(cSurface, pars, std::move(cov),
0180                                   ParticleHypothesis::pion());
0181
0182   // run the test
0183   testJacobianToGlobal(atCylinder);
0184 }
0185
0186 /// This tests the jacobian of local disc -> global
0187 BOOST_AUTO_TEST_CASE(JacobianDiscToGlobalTest) {
0188   // the disc transform and surface
0189   auto dTransform = createPlanarTransform(
0190       {10., -5., 0.}, Vector3(0.23, 0.07, 1.).normalized(), 0.004, 0.03);
0191   auto dSurface = Surface::makeShared<DiscSurface>(dTransform, 200., 1000.);
0192
0193   Covariance cov;
0194   cov << 10_mm, 0, 0, 0, 0, 0, 0, 10_mm, 0, 0, 0, 0, 0, 0, 0.1, 0, 0, 0, 0, 0,
0195       0, 0.1, 0, 0, 0, 0, 0, 0, 1. / (10_GeV), 0, 0, 0, 0, 0, 0, 0;
0196
0197   BoundVector pars;
0198   pars << 192.34, 1.823, 0.734, 0.235, 1. / (100_GeV), 0;
0199
0200   BoundTrackParameters atDisc(dSurface, pars, std::move(cov),
0201                               ParticleHypothesis::pion());
0202
0203   // run the test
0204   testJacobianToGlobal(atDisc);
0205 }
0206
0207 /// This tests the jacobian of local plane -> global
0208 BOOST_AUTO_TEST_CASE(JacobianPlaneToGlobalTest) {
0209   // Let's create a surface somewhere in space
0210   Vector3 sPosition(3421., 112., 893.);
0211   Vector3 sNormal = Vector3(1.2, -0.3, 0.05).normalized();
0212
0213   // Create a surface & parameters with covariance on the surface
0214   std::shared_ptr<PlaneSurface> pSurface =
0215       CurvilinearSurface(sPosition, sNormal).planeSurface();
0216
0217   Covariance cov;
0218   cov << 10_mm, 0, 0, 0, 0, 0, 0, 10_mm, 0, 0, 0, 0, 0, 0, 0.1, 0, 0, 0, 0, 0,
0219       0, 0.1, 0, 0, 0, 0, 0, 0, 1. / (10_GeV), 0, 0, 0, 0, 0, 0, 0;
0220
0221   BoundVector pars;
0222   pars << 12.34, -8722., 2.134, 0.85, 1. / (100_GeV), 0;
0223
0224   BoundTrackParameters atPlane(pSurface, pars, std::move(cov),
0225                                ParticleHypothesis::pion());
0226
0227   // run the test
0228   testJacobianToGlobal(atPlane);
0229 }
0230
0231 /// This tests the jacobian of local perigee -> global
0232 BOOST_AUTO_TEST_CASE(JacobianPerigeeToGlobalTest) {
0233   // Create a surface & parameters with covariance on the surface
0234   auto pSurface = Surface::makeShared<PerigeeSurface>(Vector3({0., 0., 0.}));
0235
0236   Covariance cov;
0237   cov << 10_mm, 0, 0, 0, 0, 0, 0, 10_mm, 0, 0, 0, 0, 0, 0, 0.1, 0, 0, 0, 0, 0,
0238       0, 0.1, 0, 0, 0, 0, 0, 0, 1. / (10_GeV), 0, 0, 0, 0, 0, 0, 0;
0239   BoundVector pars;
0240   pars << -3.34, -822., -0.734, 0.85, 1. / (100_GeV), 0;
0241
0242   BoundTrackParameters perigee(pSurface, pars, std::move(cov),
0243                                ParticleHypothesis::pion());
0244
0245   // run the test
0246   testJacobianToGlobal(perigee);
0247 }
0248
0249 /// This tests the jacobian of local straw -> global
0250 BOOST_AUTO_TEST_CASE(JacobianStrawToGlobalTest) {
0251   // Create a surface & parameters with covariance on the surface
0252   auto sTransform = createCylindricTransform({1019., -52., 382.}, 0.4, -0.3);
0253   auto sSurface = Surface::makeShared<StrawSurface>(sTransform, 10., 1000.);
0254
0255   Covariance cov;
0256   cov << 10_mm, 0, 0, 0, 0, 0, 0, 10_mm, 0, 0, 0, 0, 0, 0, 0.1, 0, 0, 0, 0, 0,
0257       0, 0.1, 0, 0, 0, 0, 0, 0, 1. / (10_GeV), 0, 0, 0, 0, 0, 0, 0;
0258
0259   BoundVector pars;
0260   pars << -8.34, 812., 0.734, 0.25, 1. / (100_GeV), 0;
0261
0262   BoundTrackParameters atStraw(sSurface, pars, std::move(cov),
0263                                ParticleHypothesis::pion());
0264
0265   // run the test
0266   testJacobianToGlobal(atStraw);
0267 }
0268
0269 }  // namespace Acts::Test