Core/Propagator/JacobianTests.cpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2018-2019 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 http://mozilla.org/MPL/2.0/.
0008
0009 #include <boost/test/data/test_case.hpp>
0010 #include <boost/test/tools/output_test_stream.hpp>
0011 #include <boost/test/unit_test.hpp>
0012
0013 #include "Acts/Definitions/Algebra.hpp"
0014 #include "Acts/Definitions/TrackParametrization.hpp"
0015 #include "Acts/Definitions/Units.hpp"
0016 #include "Acts/EventData/GenericCurvilinearTrackParameters.hpp"
0017 #include "Acts/EventData/TrackParameters.hpp"
0018 #include "Acts/Geometry/GeometryContext.hpp"
0019 #include "Acts/MagneticField/ConstantBField.hpp"
0020 #include "Acts/MagneticField/MagneticFieldContext.hpp"
0021 #include "Acts/Propagator/AtlasStepper.hpp"
0022 #include "Acts/Propagator/EigenStepper.hpp"
0023 #include "Acts/Surfaces/CylinderSurface.hpp"
0024 #include "Acts/Surfaces/DiscSurface.hpp"
0025 #include "Acts/Surfaces/PerigeeSurface.hpp"
0026 #include "Acts/Surfaces/PlaneSurface.hpp"
0027 #include "Acts/Surfaces/StrawSurface.hpp"
0028 #include "Acts/Surfaces/Surface.hpp"
0029 #include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
0030
0031 #include <algorithm>
0032 #include <array>
0033 #include <cstddef>
0034 #include <memory>
0035 #include <optional>
0036 #include <type_traits>
0037 #include <utility>
0038
0039 using namespace Acts::UnitLiterals;
0040
0041 namespace Acts::Test {
0042
0043 using BFieldType = ConstantBField;
0044 using EigenStepperType = EigenStepper<>;
0045 using AtlasStepperType = AtlasStepper;
0046 using Covariance = BoundSquareMatrix;
0047
0048 // Create a test context
0049 GeometryContext tgContext = GeometryContext();
0050 MagneticFieldContext mfContext = MagneticFieldContext();
0051
0052 static auto bField = std::make_shared<BFieldType>(Vector3{0, 0, 1_T});
0053
0054 /// Helper method to create a transform for a plane
0055 /// to mimic detector situations, the plane is roughly
0056 /// perpendicular to the track
0057 ///
0058 /// @param nnomal The nominal normal direction
0059 /// @param angleT Rotation around the norminal normal
0060 /// @param angleU Rotation around the original U axis
0061 Transform3 createCylindricTransform(const Vector3& nposition, double angleX,
0062                                     double angleY) {
0063   Transform3 ctransform;
0064   ctransform.setIdentity();
0065   ctransform.pretranslate(nposition);
0066   ctransform.prerotate(AngleAxis3(angleX, Vector3::UnitX()));
0067   ctransform.prerotate(AngleAxis3(angleY, Vector3::UnitY()));
0068   return ctransform;
0069 }
0070
0071 /// Helper method to create a transform for a plane
0072 /// to mimic detector situations, the plane is roughly
0073 /// perpendicular to the track
0074 ///
0075 /// @param nnomal The nominal normal direction
0076 /// @param angleT Rotation around the norminal normal
0077 /// @param angleU Rotation around the original U axis
0078 Transform3 createPlanarTransform(const Vector3& nposition,
0079                                  const Vector3& nnormal, double angleT,
0080                                  double angleU) {
0081   // the rotation of the destination surface
0082   Vector3 T = nnormal.normalized();
0083   Vector3 U = std::abs(T.dot(Vector3::UnitZ())) < 0.99
0084                   ? Vector3::UnitZ().cross(T).normalized()
0085                   : Vector3::UnitX().cross(T).normalized();
0086   Vector3 V = T.cross(U);
0087   // that's the plane curvilinear Rotation
0088   RotationMatrix3 curvilinearRotation;
0089   curvilinearRotation.col(0) = U;
0090   curvilinearRotation.col(1) = V;
0091   curvilinearRotation.col(2) = T;
0092   // curvilinear surfaces are boundless
0093   Transform3 ctransform{curvilinearRotation};
0094   ctransform.pretranslate(nposition);
0095   ctransform.prerotate(AngleAxis3(angleT, T));
0096   ctransform.prerotate(AngleAxis3(angleU, U));
0097   //
0098   return ctransform;
0099 }
0100
0101 /// Helper method : convert into Acts matrix
0102 /// It takes the double array from AtlasStepper
0103 /// and transforms it into an ActsMatrixD
0104 ///
0105 /// @param P is the pointer to the array
0106 ///
0107 /// Translation is (for lookup)
0108 ///                   /dL0    /dL1    /dPhi   /dThe   /dCM   /dT
0109 /// X  ->P[0]  dX /   P[ 8]   P[16]   P[24]   P[32]   P[40]  P[48]
0110 /// Y  ->P[1]  dY /   P[ 9]   P[17]   P[25]   P[33]   P[41]  P[49]
0111 /// Z  ->P[2]  dZ /   P[10]   P[18]   P[26]   P[34]   P[42]  P[50]
0112 /// T  ->P[3]  dT/    P[11]   P[19]   P[27]   P[35]   P[43]  P[51]
0113 /// Ax ->P[4]  dAx/   P[12]   P[20]   P[28]   P[36]   P[44]  P[52]
0114 /// Ay ->P[5]  dAy/   P[13]   P[21]   P[29]   P[37]   P[45]  P[53]
0115 /// Az ->P[6]  dAz/   P[14]   P[22]   P[30]   P[38]   P[46]  P[54]
0116 /// CM ->P[7]  dCM/   P[15]   P[23]   P[31]   P[39]   P[47]  P[55]
0117
0118 BoundToFreeMatrix convertToMatrix(const std::array<double, 60> P) {
0119   // initialize to zero
0120   BoundToFreeMatrix jMatrix = BoundToFreeMatrix::Zero();
0121   for (std::size_t j = 0; j < eBoundSize; ++j) {
0122     for (std::size_t i = 0; i < eFreeSize; ++i) {
0123       std::size_t ijc = eFreeSize + j * eFreeSize + i;
0124       jMatrix(i, j) = P[ijc];
0125     }
0126   }
0127   return jMatrix;
0128 }
0129
0130 /// Helper method : tests the jacobian to Global
0131 /// for a templated Parameters object
0132 ///
0133 /// @tparam Parameters the parameter type
0134 /// @param pars the parameter object
0135 template <typename Parameters>
0136 void testJacobianToGlobal(const Parameters& pars) {
0137   // Jacobian creation for Propagator/Steppers
0138   // a) ATLAS stepper
0139   AtlasStepperType::State astepState(tgContext, bField->makeCache(mfContext),
0140                                      pars);
0141   // b) Eigen stepper
0142   EigenStepperType::State estepState(tgContext, bField->makeCache(mfContext),
0143                                      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   CurvilinearTrackParameters curvilinear(Vector4(341., 412., 93., 0.),
0159                                          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   auto pSurface = Surface::makeShared<PlaneSurface>(sPosition, sNormal);
0215
0216   Covariance cov;
0217   cov << 10_mm, 0, 0, 0, 0, 0, 0, 10_mm, 0, 0, 0, 0, 0, 0, 0.1, 0, 0, 0, 0, 0,
0218       0, 0.1, 0, 0, 0, 0, 0, 0, 1. / (10_GeV), 0, 0, 0, 0, 0, 0, 0;
0219
0220   BoundVector pars;
0221   pars << 12.34, -8722., 2.134, 0.85, 1. / (100_GeV), 0;
0222
0223   BoundTrackParameters atPlane(pSurface, pars, std::move(cov),
0224                                ParticleHypothesis::pion());
0225
0226   // run the test
0227   testJacobianToGlobal(atPlane);
0228 }
0229
0230 /// This tests the jacobian of local perigee -> global
0231 BOOST_AUTO_TEST_CASE(JacobianPerigeeToGlobalTest) {
0232   // Create a surface & parameters with covariance on the surface
0233   auto pSurface = Surface::makeShared<PerigeeSurface>(Vector3({0., 0., 0.}));
0234
0235   Covariance cov;
0236   cov << 10_mm, 0, 0, 0, 0, 0, 0, 10_mm, 0, 0, 0, 0, 0, 0, 0.1, 0, 0, 0, 0, 0,
0237       0, 0.1, 0, 0, 0, 0, 0, 0, 1. / (10_GeV), 0, 0, 0, 0, 0, 0, 0;
0238   BoundVector pars;
0239   pars << -3.34, -822., -0.734, 0.85, 1. / (100_GeV), 0;
0240
0241   BoundTrackParameters perigee(pSurface, pars, std::move(cov),
0242                                ParticleHypothesis::pion());
0243
0244   // run the test
0245   testJacobianToGlobal(perigee);
0246 }
0247
0248 /// This tests the jacobian of local straw -> global
0249 BOOST_AUTO_TEST_CASE(JacobianStrawToGlobalTest) {
0250   // Create a surface & parameters with covariance on the surface
0251   auto sTransform = createCylindricTransform({1019., -52., 382.}, 0.4, -0.3);
0252   auto sSurface = Surface::makeShared<StrawSurface>(sTransform, 10., 1000.);
0253
0254   Covariance cov;
0255   cov << 10_mm, 0, 0, 0, 0, 0, 0, 10_mm, 0, 0, 0, 0, 0, 0, 0.1, 0, 0, 0, 0, 0,
0256       0, 0.1, 0, 0, 0, 0, 0, 0, 1. / (10_GeV), 0, 0, 0, 0, 0, 0, 0;
0257
0258   BoundVector pars;
0259   pars << -8.34, 812., 0.734, 0.25, 1. / (100_GeV), 0;
0260
0261   BoundTrackParameters atStraw(sSurface, pars, std::move(cov),
0262                                ParticleHypothesis::pion());
0263
0264   // run the test
0265   testJacobianToGlobal(atStraw);
0266 }
0267
0268 }  // namespace Acts::Test