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), pars);
0138   // b) Eigen stepper
0139   EigenStepperType estep(bField);
0140   EigenStepperType::State estepState =
0141       estep.makeState(EigenStepperType::Options(tgContext, mfContext), pars);
0142
0143   // create the matrices
0144   auto asMatrix = convertToMatrix(astepState.pVector);
0145
0146   // cross comparison checks
0147   CHECK_CLOSE_OR_SMALL(asMatrix, estepState.jacToGlobal, 1e-6, 1e-9);
0148 }
0149
0150 /// This tests the jacobian of local curvilinear -> global
0151 BOOST_AUTO_TEST_CASE(JacobianCurvilinearToGlobalTest) {
0152   // Create curvilinear parameters
0153   Covariance cov;
0154   cov << 10_mm, 0, 0, 0, 0, 0, 0, 10_mm, 0, 0, 0, 0, 0, 0, 0.1, 0, 0, 0, 0, 0,
0155       0, 0.1, 0, 0, 0, 0, 0, 0, 1. / (10_GeV), 0, 0, 0, 0, 0, 0, 0;
0156   CurvilinearTrackParameters curvilinear(Vector4(341., 412., 93., 0.),
0157                                          Vector3(1.2, 8.3, 0.45), 1 / 10.0, cov,
0158                                          ParticleHypothesis::pion());
0159
0160   // run the test
0161   testJacobianToGlobal(curvilinear);
0162 }
0163
0164 /// This tests the jacobian of local cylinder -> global
0165 BOOST_AUTO_TEST_CASE(JacobianCylinderToGlobalTest) {
0166   // the cylinder transform and surface
0167   auto cTransform = createCylindricTransform({10., -5., 0.}, 0.004, 0.03);
0168   auto cSurface = Surface::makeShared<CylinderSurface>(cTransform, 200., 1000.);
0169
0170   Covariance cov;
0171   cov << 10_mm, 0, 0, 0, 0, 0, 0, 10_mm, 0, 0, 0, 0, 0, 0, 0.1, 0, 0, 0, 0, 0,
0172       0, 0.1, 0, 0, 0, 0, 0, 0, 1. / (10_GeV), 0, 0, 0, 0, 0, 0, 0;
0173
0174   BoundVector pars;
0175   pars << 182.34, -82., 0.134, 0.85, 1. / (100_GeV), 0;
0176
0177   BoundTrackParameters atCylinder(cSurface, pars, std::move(cov),
0178                                   ParticleHypothesis::pion());
0179
0180   // run the test
0181   testJacobianToGlobal(atCylinder);
0182 }
0183
0184 /// This tests the jacobian of local disc -> global
0185 BOOST_AUTO_TEST_CASE(JacobianDiscToGlobalTest) {
0186   // the disc transform and surface
0187   auto dTransform = createPlanarTransform(
0188       {10., -5., 0.}, Vector3(0.23, 0.07, 1.).normalized(), 0.004, 0.03);
0189   auto dSurface = Surface::makeShared<DiscSurface>(dTransform, 200., 1000.);
0190
0191   Covariance cov;
0192   cov << 10_mm, 0, 0, 0, 0, 0, 0, 10_mm, 0, 0, 0, 0, 0, 0, 0.1, 0, 0, 0, 0, 0,
0193       0, 0.1, 0, 0, 0, 0, 0, 0, 1. / (10_GeV), 0, 0, 0, 0, 0, 0, 0;
0194
0195   BoundVector pars;
0196   pars << 192.34, 1.823, 0.734, 0.235, 1. / (100_GeV), 0;
0197
0198   BoundTrackParameters atDisc(dSurface, pars, std::move(cov),
0199                               ParticleHypothesis::pion());
0200
0201   // run the test
0202   testJacobianToGlobal(atDisc);
0203 }
0204
0205 /// This tests the jacobian of local plane -> global
0206 BOOST_AUTO_TEST_CASE(JacobianPlaneToGlobalTest) {
0207   // Let's create a surface somewhere in space
0208   Vector3 sPosition(3421., 112., 893.);
0209   Vector3 sNormal = Vector3(1.2, -0.3, 0.05).normalized();
0210
0211   // Create a surface & parameters with covariance on the surface
0212   auto pSurface = CurvilinearSurface(sPosition, sNormal).planeSurface();
0213
0214   Covariance cov;
0215   cov << 10_mm, 0, 0, 0, 0, 0, 0, 10_mm, 0, 0, 0, 0, 0, 0, 0.1, 0, 0, 0, 0, 0,
0216       0, 0.1, 0, 0, 0, 0, 0, 0, 1. / (10_GeV), 0, 0, 0, 0, 0, 0, 0;
0217
0218   BoundVector pars;
0219   pars << 12.34, -8722., 2.134, 0.85, 1. / (100_GeV), 0;
0220
0221   BoundTrackParameters atPlane(pSurface, pars, std::move(cov),
0222                                ParticleHypothesis::pion());
0223
0224   // run the test
0225   testJacobianToGlobal(atPlane);
0226 }
0227
0228 /// This tests the jacobian of local perigee -> global
0229 BOOST_AUTO_TEST_CASE(JacobianPerigeeToGlobalTest) {
0230   // Create a surface & parameters with covariance on the surface
0231   auto pSurface = Surface::makeShared<PerigeeSurface>(Vector3({0., 0., 0.}));
0232
0233   Covariance cov;
0234   cov << 10_mm, 0, 0, 0, 0, 0, 0, 10_mm, 0, 0, 0, 0, 0, 0, 0.1, 0, 0, 0, 0, 0,
0235       0, 0.1, 0, 0, 0, 0, 0, 0, 1. / (10_GeV), 0, 0, 0, 0, 0, 0, 0;
0236   BoundVector pars;
0237   pars << -3.34, -822., -0.734, 0.85, 1. / (100_GeV), 0;
0238
0239   BoundTrackParameters perigee(pSurface, pars, std::move(cov),
0240                                ParticleHypothesis::pion());
0241
0242   // run the test
0243   testJacobianToGlobal(perigee);
0244 }
0245
0246 /// This tests the jacobian of local straw -> global
0247 BOOST_AUTO_TEST_CASE(JacobianStrawToGlobalTest) {
0248   // Create a surface & parameters with covariance on the surface
0249   auto sTransform = createCylindricTransform({1019., -52., 382.}, 0.4, -0.3);
0250   auto sSurface = Surface::makeShared<StrawSurface>(sTransform, 10., 1000.);
0251
0252   Covariance cov;
0253   cov << 10_mm, 0, 0, 0, 0, 0, 0, 10_mm, 0, 0, 0, 0, 0, 0, 0.1, 0, 0, 0, 0, 0,
0254       0, 0.1, 0, 0, 0, 0, 0, 0, 1. / (10_GeV), 0, 0, 0, 0, 0, 0, 0;
0255
0256   BoundVector pars;
0257   pars << -8.34, 812., 0.734, 0.25, 1. / (100_GeV), 0;
0258
0259   BoundTrackParameters atStraw(sSurface, pars, std::move(cov),
0260                                ParticleHypothesis::pion());
0261
0262   // run the test
0263   testJacobianToGlobal(atStraw);
0264 }
0265
0266 }  // namespace Acts::Test