Core/Utilities/JacobianHelpersTests.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/Utilities/JacobianHelpers.hpp"
0013 #include "Acts/Utilities/UnitVectors.hpp"
0014 #include "ActsTests/CommonHelpers/FloatComparisons.hpp"
0015
0016 #include <cmath>
0017 #include <limits>
0018 #include <numbers>
0019
0020 using namespace Acts;
0021
0022 namespace {
0023 constexpr double eps = 1e-6;
0024
0025 // d(p_x, p_y, p_z) / d(phi, theta, qOverP) evaluated by central finite
0026 // differences of the explicit spherical-to-Cartesian momentum map.
0027 Matrix<3, 3> numericSphericalToFreeMomentumJacobian(double phi, double theta,
0028                                                     double qOverP,
0029                                                     double charge) {
0030   auto momentum = [&](double phi_, double theta_, double qOverP_) -> Vector3 {
0031     return (charge / qOverP_) * makeDirectionFromPhiTheta(phi_, theta_);
0032   };
0033
0034   const double h = 1e-6;
0035   Matrix<3, 3> jacobian;
0036   jacobian.col(0) =
0037       (momentum(phi + h, theta, qOverP) - momentum(phi - h, theta, qOverP)) /
0038       (2 * h);
0039   jacobian.col(1) =
0040       (momentum(phi, theta + h, qOverP) - momentum(phi, theta - h, qOverP)) /
0041       (2 * h);
0042   jacobian.col(2) =
0043       (momentum(phi, theta, qOverP + h) - momentum(phi, theta, qOverP - h)) /
0044       (2 * h);
0045   return jacobian;
0046 }
0047 }  // namespace
0048
0049 namespace ActsTests {
0050
0051 BOOST_AUTO_TEST_SUITE(JacobianHelpersSuite)
0052
0053 BOOST_AUTO_TEST_CASE(SphericalToFreeMomentumJacobianMatchesFiniteDifference) {
0054   const double phi = 0.6;
0055   const double theta = 1.1;
0056   const double charge = 1.;
0057   const double qOverP = 0.5;
0058   const double momentum = std::abs(charge / qOverP);
0059
0060   const Vector3 direction = makeDirectionFromPhiTheta(phi, theta);
0061
0062   const Matrix<3, 3> analytic =
0063       sphericalToFreeMomentumJacobian(direction, qOverP, momentum);
0064   const Matrix<3, 3> numeric =
0065       numericSphericalToFreeMomentumJacobian(phi, theta, qOverP, charge);
0066
0067   CHECK_CLOSE_ABS(analytic, numeric, eps);
0068 }
0069
0070 BOOST_AUTO_TEST_CASE(FreeToSphericalMomentumJacobianIsInverse) {
0071   const double phi = -1.2;
0072   const double theta = 0.8;
0073   const double charge = -1.;
0074   const double qOverP = -0.25;
0075   const double momentum = std::abs(charge / qOverP);
0076
0077   const Vector3 direction = makeDirectionFromPhiTheta(phi, theta);
0078   const Vector3 momentumVector = momentum * direction;
0079
0080   const Matrix<3, 3> forward =
0081       sphericalToFreeMomentumJacobian(direction, qOverP, momentum);
0082   const Matrix<3, 3> inverse =
0083       freeToSphericalMomentumJacobian(momentumVector, charge);
0084
0085   const Matrix<3, 3> identity = Matrix<3, 3>::Identity();
0086
0087   CHECK_CLOSE_ABS(inverse * forward, identity, eps);
0088   CHECK_CLOSE_ABS(forward * inverse, identity, eps);
0089 }
0090
0091 BOOST_AUTO_TEST_SUITE_END()
0092
0093 }  // namespace ActsTests