From c51f48a7a667cdb54fe912420303fb13455d81df Mon Sep 17 00:00:00 2001 From: Andreas Stefl Date: Thu, 21 Dec 2023 22:49:44 +0100 Subject: [PATCH 1/2] refactor!: Remove navigator layer bounds check option (#2851) This option is supposed to improve the navigation but can also cause regressions. Layers without bounds overlap with boundaries therefore we might end up outside the current tracking volume while targeting a layer. I propose that this should be rather fixed by the used with proper layer bounds that extend to the volume boundaries. --- Core/include/Acts/Propagator/Navigator.hpp | 5 --- .../PropagationBentTracks.cpp | 37 ------------------- 2 files changed, 42 deletions(-) diff --git a/Core/include/Acts/Propagator/Navigator.hpp b/Core/include/Acts/Propagator/Navigator.hpp index 2babb8d5e3c..7a6725fe930 100644 --- a/Core/include/Acts/Propagator/Navigator.hpp +++ b/Core/include/Acts/Propagator/Navigator.hpp @@ -113,10 +113,6 @@ class Navigator { bool resolveMaterial = true; /// stop at every surface regardless what it is bool resolvePassive = false; - - /// Whether to perform boundary checks for layer resolving (improves - /// navigation for bended tracks) - BoundaryCheck boundaryCheckLayerResolving = BoundaryCheck(true); }; /// @brief Nested State struct @@ -1187,7 +1183,6 @@ class Navigator { // Create the navigation options // - and get the compatible layers, start layer will be excluded NavigationOptions navOpts; - navOpts.boundaryCheck = m_cfg.boundaryCheckLayerResolving; navOpts.resolveSensitive = m_cfg.resolveSensitive; navOpts.resolveMaterial = m_cfg.resolveMaterial; navOpts.resolvePassive = m_cfg.resolvePassive; diff --git a/Tests/IntegrationTests/PropagationBentTracks.cpp b/Tests/IntegrationTests/PropagationBentTracks.cpp index 02edb9ac953..3aac989681e 100644 --- a/Tests/IntegrationTests/PropagationBentTracks.cpp +++ b/Tests/IntegrationTests/PropagationBentTracks.cpp @@ -89,25 +89,6 @@ BOOST_AUTO_TEST_CASE(with_boundary_check_no_bfield) { 0); } -BOOST_AUTO_TEST_CASE(without_boundary_check_no_bfield) { - auto navCfg = Acts::Navigator::Config{}; - navCfg.boundaryCheckLayerResolving = Acts::BoundaryCheck(false); - const auto xPositions = xPositionsOfPassedSurfaces(navCfg, 0.0_T); - - // without bField we exit at the side so we don't hit the surfaces at x ~ - // 2000 and also not the boundary surface at x = 3000, regardless of the - // boundary checking - BOOST_CHECK_EQUAL(std::count(xPositions.begin(), xPositions.end(), 999.0), 1); - BOOST_CHECK_EQUAL(std::count(xPositions.begin(), xPositions.end(), 1001.0), - 1); - BOOST_CHECK_EQUAL(std::count(xPositions.begin(), xPositions.end(), 1999.0), - 0); - BOOST_CHECK_EQUAL(std::count(xPositions.begin(), xPositions.end(), 2001.0), - 0); - BOOST_CHECK_EQUAL(std::count(xPositions.begin(), xPositions.end(), 3000.0), - 0); -} - BOOST_AUTO_TEST_CASE(with_boundary_check_with_bfield) { auto navCfg = Acts::Navigator::Config{}; const auto xPositions = xPositionsOfPassedSurfaces(navCfg, 0.5_T); @@ -124,21 +105,3 @@ BOOST_AUTO_TEST_CASE(with_boundary_check_with_bfield) { BOOST_CHECK_EQUAL(std::count(xPositions.begin(), xPositions.end(), 3000.0), 1); } - -BOOST_AUTO_TEST_CASE(no_boundary_check_with_bfield) { - auto navCfg = Acts::Navigator::Config{}; - navCfg.boundaryCheckLayerResolving = Acts::BoundaryCheck(false); - const auto xPositions = xPositionsOfPassedSurfaces(navCfg, 0.5_T); - - // Without boundary check at layer resolving, we also hit the surfaces at x ~ - // 2000 - BOOST_CHECK_EQUAL(std::count(xPositions.begin(), xPositions.end(), 999.0), 1); - BOOST_CHECK_EQUAL(std::count(xPositions.begin(), xPositions.end(), 1001.0), - 1); - BOOST_CHECK_EQUAL(std::count(xPositions.begin(), xPositions.end(), 1999.0), - 1); - BOOST_CHECK_EQUAL(std::count(xPositions.begin(), xPositions.end(), 2001.0), - 1); - BOOST_CHECK_EQUAL(std::count(xPositions.begin(), xPositions.end(), 3000.0), - 1); -} From 7b26be28bd15b1a27dddff1f493ae1f838088ceb Mon Sep 17 00:00:00 2001 From: Andreas Stefl Date: Fri, 22 Dec 2023 12:03:01 +0100 Subject: [PATCH 2/2] refactor: Use common direction transform Jacobian (#2782) I found the same code in a few places and there seems to be a bit of simplification opportunity in the math --- .../Acts/Utilities/JacobianHelpers.hpp | 68 +++++++++++++++++++ Core/include/Acts/Utilities/VectorHelpers.hpp | 6 +- .../Propagator/detail/CovarianceEngine.cpp | 65 +++++++----------- Core/src/Propagator/detail/JacobianEngine.cpp | 10 +-- Core/src/Surfaces/DiscSurface.cpp | 44 ++++-------- Core/src/Surfaces/LineSurface.cpp | 13 +--- Core/src/Surfaces/Surface.cpp | 21 ++---- Examples/Python/tests/root_file_hashes.txt | 28 ++++---- 8 files changed, 137 insertions(+), 118 deletions(-) create mode 100644 Core/include/Acts/Utilities/JacobianHelpers.hpp diff --git a/Core/include/Acts/Utilities/JacobianHelpers.hpp b/Core/include/Acts/Utilities/JacobianHelpers.hpp new file mode 100644 index 00000000000..a666fe2ecfe --- /dev/null +++ b/Core/include/Acts/Utilities/JacobianHelpers.hpp @@ -0,0 +1,68 @@ +// This file is part of the Acts project. +// +// Copyright (C) 2023 CERN for the benefit of the Acts project +// +// This Source Code Form is subject to the terms of the Mozilla Public +// License, v. 2.0. If a copy of the MPL was not distributed with this +// file, You can obtain one at http://mozilla.org/MPL/2.0/. + +#pragma once + +#include "Acts/Definitions/Algebra.hpp" +#include "Acts/Utilities/VectorHelpers.hpp" + +namespace Acts { + +/// @brief Calculates the Jacobian for spherical to free +/// direction vector transformation +/// +/// @note We use the direction vector as an input because +/// the trigonometric simplify that way +/// +/// @param direction The normalised direction vector +/// +/// @return The Jacobian d(dir_x, dir_y, dir_z) / d(phi, theta) +/// +inline ActsMatrix<3, 2> sphericalToFreeDirectionJacobian( + const Vector3& direction) { + auto [cosPhi, sinPhi, cosTheta, sinTheta] = + VectorHelpers::evaluateTrigonomics(direction); + + // clang-format off + ActsMatrix<3, 2> jacobian; + jacobian << + -direction.y(), cosTheta * cosPhi, + direction.x(), cosTheta * sinPhi, + 0, -sinTheta; + // clang-format on + + return jacobian; +} + +/// @brief Calculates the Jacobian for free to spherical +/// direction vector transformation +/// +/// @note We use the direction vector as an input because +/// the trigonometric simplify that way +/// +/// @param direction The normalised direction vector +/// +/// @return The Jacobian d(phi, theta) / d(dir_x, dir_y, dir_z) +/// +inline ActsMatrix<2, 3> freeToSphericalDirectionJacobian( + const Vector3& direction) { + auto [cosPhi, sinPhi, cosTheta, sinTheta] = + VectorHelpers::evaluateTrigonomics(direction); + ActsScalar invSinTheta = 1. / sinTheta; + + // clang-format off + ActsMatrix<2, 3> jacobian; + jacobian << + -sinPhi * invSinTheta, cosPhi * invSinTheta, 0, + cosPhi * cosTheta, sinPhi * cosTheta, -sinTheta; + // clang-format on + + return jacobian; +} + +} // namespace Acts diff --git a/Core/include/Acts/Utilities/VectorHelpers.hpp b/Core/include/Acts/Utilities/VectorHelpers.hpp index 42208ffd60f..9760343f71b 100644 --- a/Core/include/Acts/Utilities/VectorHelpers.hpp +++ b/Core/include/Acts/Utilities/VectorHelpers.hpp @@ -132,18 +132,18 @@ double eta(const Eigen::MatrixBase& v) noexcept { /// @param direction for this evaluatoin /// /// @return cos(phi), sin(phi), cos(theta), sin(theta), 1/sin(theta) -inline std::array evaluateTrigonomics(const Vector3& direction) { +inline std::array evaluateTrigonomics(const Vector3& direction) { const ActsScalar x = direction(0); // == cos(phi) * sin(theta) const ActsScalar y = direction(1); // == sin(phi) * sin(theta) const ActsScalar z = direction(2); // == cos(theta) // can be turned into cosine/sine const ActsScalar cosTheta = z; - const ActsScalar sinTheta = std::hypot(x, y); + const ActsScalar sinTheta = std::sqrt(1 - z * z); const ActsScalar invSinTheta = 1. / sinTheta; const ActsScalar cosPhi = x * invSinTheta; const ActsScalar sinPhi = y * invSinTheta; - return {cosPhi, sinPhi, cosTheta, sinTheta, invSinTheta}; + return {cosPhi, sinPhi, cosTheta, sinTheta}; } /// Helper method to extract the binning value from a 3D vector. diff --git a/Core/src/Propagator/detail/CovarianceEngine.cpp b/Core/src/Propagator/detail/CovarianceEngine.cpp index 6db4a7af236..e196f24a97f 100644 --- a/Core/src/Propagator/detail/CovarianceEngine.cpp +++ b/Core/src/Propagator/detail/CovarianceEngine.cpp @@ -10,6 +10,7 @@ #include "Acts/Definitions/Common.hpp" #include "Acts/Definitions/Tolerance.hpp" +#include "Acts/Definitions/TrackParametrization.hpp" #include "Acts/EventData/GenericBoundTrackParameters.hpp" #include "Acts/EventData/GenericCurvilinearTrackParameters.hpp" #include "Acts/EventData/detail/CorrectedTransformationFreeToBound.hpp" @@ -17,6 +18,7 @@ #include "Acts/EventData/detail/TransformationFreeToBound.hpp" #include "Acts/Propagator/detail/JacobianEngine.hpp" #include "Acts/Utilities/AlgebraHelpers.hpp" +#include "Acts/Utilities/JacobianHelpers.hpp" #include "Acts/Utilities/Result.hpp" #include @@ -36,6 +38,7 @@ using CurvilinearState = std::tuple; /// @brief Evaluate the projection Jacobian from free to curvilinear parameters +/// without transport jacobian. /// /// @param [in] direction Normalised direction vector /// @@ -45,41 +48,36 @@ FreeToBoundMatrix freeToCurvilinearJacobian(const Vector3& direction) { const double x = direction(0); // == cos(phi) * sin(theta) const double y = direction(1); // == sin(phi) * sin(theta) const double z = direction(2); // == cos(theta) - // can be turned into cosine/sine - const double cosTheta = z; - const double sinTheta = std::hypot(x, y); - const double invSinTheta = 1. / sinTheta; - const double cosPhi = x * invSinTheta; - const double sinPhi = y * invSinTheta; // prepare the jacobian to curvilinear FreeToBoundMatrix jacToCurv = FreeToBoundMatrix::Zero(); - if (std::abs(cosTheta) < s_curvilinearProjTolerance) { + if (std::abs(z) < s_curvilinearProjTolerance) { + auto [cosPhi, sinPhi, cosTheta, sinTheta] = + VectorHelpers::evaluateTrigonomics(direction); // We normally operate in curvilinear coordinates defined as follows - jacToCurv(0, 0) = -sinPhi; - jacToCurv(0, 1) = cosPhi; - jacToCurv(1, 0) = -cosPhi * cosTheta; - jacToCurv(1, 1) = -sinPhi * cosTheta; - jacToCurv(1, 2) = sinTheta; + jacToCurv(eBoundLoc0, eFreePos0) = -sinPhi; + jacToCurv(eBoundLoc0, eFreePos1) = cosPhi; + // jacToCurv(eBoundLoc0, eFreePos2) = 0; + jacToCurv(eBoundLoc1, eFreePos0) = -cosPhi * cosTheta; + jacToCurv(eBoundLoc1, eFreePos1) = -sinPhi * cosTheta; + jacToCurv(eBoundLoc1, eFreePos2) = sinTheta; } else { // Under grazing incidence to z, the above coordinate system definition // becomes numerically unstable, and we need to switch to another one const double c = std::hypot(y, z); const double invC = 1. / c; - jacToCurv(0, 1) = -z * invC; - jacToCurv(0, 2) = y * invC; - jacToCurv(1, 0) = c; - jacToCurv(1, 1) = -x * y * invC; - jacToCurv(1, 2) = -x * z * invC; + // jacToCurv(eBoundLoc0, eFreePos0) = 0; + jacToCurv(eBoundLoc0, eFreePos1) = -z * invC; + jacToCurv(eBoundLoc0, eFreePos2) = y * invC; + jacToCurv(eBoundLoc1, eFreePos0) = c; + jacToCurv(eBoundLoc1, eFreePos1) = -x * y * invC; + jacToCurv(eBoundLoc1, eFreePos2) = -x * z * invC; } // Time parameter - jacToCurv(5, 3) = 1.; + jacToCurv(eBoundTime, eFreeTime) = 1.; // Directional and momentum parameters for curvilinear - jacToCurv(2, 4) = -sinPhi * invSinTheta; - jacToCurv(2, 5) = cosPhi * invSinTheta; - jacToCurv(3, 4) = cosPhi * cosTheta; - jacToCurv(3, 5) = sinPhi * cosTheta; - jacToCurv(3, 6) = -sinTheta; - jacToCurv(4, 7) = 1.; + jacToCurv.block<2, 3>(eBoundPhi, eFreeDir0) = + freeToSphericalDirectionJacobian(direction); + jacToCurv(eBoundQOverP, eFreeQOverP) = 1.; return jacToCurv; } @@ -237,16 +235,8 @@ void reinitializeJacobians(FreeMatrix& freeTransportJacobian, freeToPathDerivatives = FreeVector::Zero(); boundToFreeJacobian = BoundToFreeMatrix::Zero(); - // Optimized trigonometry on the propagation direction - const double x = direction(0); // == cos(phi) * sin(theta) - const double y = direction(1); // == sin(phi) * sin(theta) - const double z = direction(2); // == cos(theta) - // can be turned into cosine/sine - const double cosTheta = z; - const double sinTheta = std::hypot(x, y); - const double invSinTheta = 1. / sinTheta; - const double cosPhi = x * invSinTheta; - const double sinPhi = y * invSinTheta; + auto [cosPhi, sinPhi, cosTheta, sinTheta] = + VectorHelpers::evaluateTrigonomics(direction); boundToFreeJacobian(eFreePos0, eBoundLoc0) = -sinPhi; boundToFreeJacobian(eFreePos0, eBoundLoc1) = -cosPhi * cosTheta; @@ -254,11 +244,8 @@ void reinitializeJacobians(FreeMatrix& freeTransportJacobian, boundToFreeJacobian(eFreePos1, eBoundLoc1) = -sinPhi * cosTheta; boundToFreeJacobian(eFreePos2, eBoundLoc1) = sinTheta; boundToFreeJacobian(eFreeTime, eBoundTime) = 1; - boundToFreeJacobian(eFreeDir0, eBoundPhi) = -sinTheta * sinPhi; - boundToFreeJacobian(eFreeDir0, eBoundTheta) = cosTheta * cosPhi; - boundToFreeJacobian(eFreeDir1, eBoundPhi) = sinTheta * cosPhi; - boundToFreeJacobian(eFreeDir1, eBoundTheta) = cosTheta * sinPhi; - boundToFreeJacobian(eFreeDir2, eBoundTheta) = -sinTheta; + boundToFreeJacobian.block<3, 2>(eFreeDir0, eBoundPhi) = + sphericalToFreeDirectionJacobian(direction); boundToFreeJacobian(eFreeQOverP, eBoundQOverP) = 1; } } // namespace diff --git a/Core/src/Propagator/detail/JacobianEngine.cpp b/Core/src/Propagator/detail/JacobianEngine.cpp index e215f58fe4f..12c939b005b 100644 --- a/Core/src/Propagator/detail/JacobianEngine.cpp +++ b/Core/src/Propagator/detail/JacobianEngine.cpp @@ -22,8 +22,9 @@ namespace Acts { namespace detail { FreeToBoundMatrix freeToCurvilinearJacobian(const Vector3& direction) { - auto [cosPhi, sinPhi, cosTheta, sinTheta, invSinTheta] = + auto [cosPhi, sinPhi, cosTheta, sinTheta] = VectorHelpers::evaluateTrigonomics(direction); + ActsScalar invSinTheta = 1. / sinTheta; // Prepare the jacobian to curvilinear FreeToBoundMatrix freeToCurvJacobian = FreeToBoundMatrix::Zero(); if (std::abs(cosTheta) < s_curvilinearProjTolerance) { @@ -61,7 +62,7 @@ FreeToBoundMatrix freeToCurvilinearJacobian(const Vector3& direction) { } BoundToFreeMatrix curvilinearToFreeJacobian(const Vector3& direction) { - auto [cosPhi, sinPhi, cosTheta, sinTheta, invSinTheta] = + auto [cosPhi, sinPhi, cosTheta, sinTheta] = VectorHelpers::evaluateTrigonomics(direction); // Prepare the jacobian to free @@ -88,7 +89,7 @@ BoundToFreeMatrix curvilinearToFreeJacobian(const Vector3& direction) { ActsMatrix<8, 7> anglesToDirectionJacobian(const Vector3& direction) { ActsMatrix<8, 7> jacobian = ActsMatrix<8, 7>::Zero(); - auto [cosPhi, sinPhi, cosTheta, sinTheta, invSinTheta] = + auto [cosPhi, sinPhi, cosTheta, sinTheta] = VectorHelpers::evaluateTrigonomics(direction); jacobian(0, 0) = 1.; @@ -108,8 +109,9 @@ ActsMatrix<8, 7> anglesToDirectionJacobian(const Vector3& direction) { ActsMatrix<7, 8> directionToAnglesJacobian(const Vector3& direction) { ActsMatrix<7, 8> jacobian = ActsMatrix<7, 8>::Zero(); - auto [cosPhi, sinPhi, cosTheta, sinTheta, invSinTheta] = + auto [cosPhi, sinPhi, cosTheta, sinTheta] = VectorHelpers::evaluateTrigonomics(direction); + ActsScalar invSinTheta = 1. / sinTheta; jacobian(0, 0) = 1.; jacobian(1, 1) = 1.; diff --git a/Core/src/Surfaces/DiscSurface.cpp b/Core/src/Surfaces/DiscSurface.cpp index a9140343663..78c5c773ed1 100644 --- a/Core/src/Surfaces/DiscSurface.cpp +++ b/Core/src/Surfaces/DiscSurface.cpp @@ -21,6 +21,7 @@ #include "Acts/Surfaces/detail/PlanarHelper.hpp" #include "Acts/Utilities/Helpers.hpp" #include "Acts/Utilities/Intersection.hpp" +#include "Acts/Utilities/JacobianHelpers.hpp" #include "Acts/Utilities/ThrowAssert.hpp" #include @@ -218,34 +219,26 @@ Acts::BoundToFreeMatrix Acts::DiscSurface::boundToFreeJacobian( const Vector3 position = freeParams.segment<3>(eFreePos0); // The direction const Vector3 direction = freeParams.segment<3>(eFreeDir0); - // Get the sines and cosines directly - const double cos_theta = std::cos(boundParams[eBoundTheta]); - const double sin_theta = std::sin(boundParams[eBoundTheta]); - const double cos_phi = std::cos(boundParams[eBoundPhi]); - const double sin_phi = std::sin(boundParams[eBoundPhi]); // special polar coordinates for the Disc double lrad = boundParams[eBoundLoc0]; double lphi = boundParams[eBoundLoc1]; - double lcos_phi = cos(lphi); - double lsin_phi = sin(lphi); + double lcphi = std::cos(lphi); + double lsphi = std::sin(lphi); // retrieve the reference frame const auto rframe = referenceFrame(gctx, position, direction); // Initialize the jacobian from local to global BoundToFreeMatrix jacToGlobal = BoundToFreeMatrix::Zero(); // the local error components - rotated from reference frame jacToGlobal.block<3, 1>(eFreePos0, eBoundLoc0) = - lcos_phi * rframe.block<3, 1>(0, 0) + lsin_phi * rframe.block<3, 1>(0, 1); + lcphi * rframe.block<3, 1>(0, 0) + lsphi * rframe.block<3, 1>(0, 1); jacToGlobal.block<3, 1>(eFreePos0, eBoundLoc1) = - lrad * (lcos_phi * rframe.block<3, 1>(0, 1) - - lsin_phi * rframe.block<3, 1>(0, 0)); + lrad * + (lcphi * rframe.block<3, 1>(0, 1) - lsphi * rframe.block<3, 1>(0, 0)); // the time component jacToGlobal(eFreeTime, eBoundTime) = 1; // the momentum components - jacToGlobal(eFreeDir0, eBoundPhi) = (-sin_theta) * sin_phi; - jacToGlobal(eFreeDir0, eBoundTheta) = cos_theta * cos_phi; - jacToGlobal(eFreeDir1, eBoundPhi) = sin_theta * cos_phi; - jacToGlobal(eFreeDir1, eBoundTheta) = cos_theta * sin_phi; - jacToGlobal(eFreeDir2, eBoundTheta) = (-sin_theta); + jacToGlobal.block<3, 2>(eFreeDir0, eBoundPhi) = + sphericalToFreeDirectionJacobian(direction); jacToGlobal(eFreeQOverP, eBoundQOverP) = 1; return jacToGlobal; } @@ -258,16 +251,6 @@ Acts::FreeToBoundMatrix Acts::DiscSurface::freeToBoundJacobian( const auto position = parameters.segment<3>(eFreePos0); // The direction const auto direction = parameters.segment<3>(eFreeDir0); - // Optimized trigonometry on the propagation direction - const double x = direction(0); // == cos(phi) * sin(theta) - const double y = direction(1); // == sin(phi) * sin(theta) - const double z = direction(2); // == cos(theta) - // can be turned into cosine/sine - const double cosTheta = z; - const double sinTheta = std::hypot(x, y); - const double invSinTheta = 1. / sinTheta; - const double cosPhi = x * invSinTheta; - const double sinPhi = y * invSinTheta; // The measurement frame of the surface RotationMatrix3 rframeT = referenceFrame(gctx, position, direction).transpose(); @@ -275,8 +258,8 @@ Acts::FreeToBoundMatrix Acts::DiscSurface::freeToBoundJacobian( const Vector3 pos_loc = transform(gctx).inverse() * position; const double lr = perp(pos_loc); const double lphi = phi(pos_loc); - const double lcphi = cos(lphi); - const double lsphi = sin(lphi); + const double lcphi = std::cos(lphi); + const double lsphi = std::sin(lphi); // rotate into the polar coorindates auto lx = rframeT.block<1, 3>(0, 0); auto ly = rframeT.block<1, 3>(1, 0); @@ -289,11 +272,8 @@ Acts::FreeToBoundMatrix Acts::DiscSurface::freeToBoundJacobian( // Time element jacToLocal(eBoundTime, eFreeTime) = 1; // Directional and momentum elements for reference frame surface - jacToLocal(eBoundPhi, eFreeDir0) = -sinPhi * invSinTheta; - jacToLocal(eBoundPhi, eFreeDir1) = cosPhi * invSinTheta; - jacToLocal(eBoundTheta, eFreeDir0) = cosPhi * cosTheta; - jacToLocal(eBoundTheta, eFreeDir1) = sinPhi * cosTheta; - jacToLocal(eBoundTheta, eFreeDir2) = -sinTheta; + jacToLocal.block<2, 3>(eBoundPhi, eFreeDir0) = + freeToSphericalDirectionJacobian(direction); jacToLocal(eBoundQOverP, eFreeQOverP) = 1; return jacToLocal; } diff --git a/Core/src/Surfaces/LineSurface.cpp b/Core/src/Surfaces/LineSurface.cpp index e6d0c045cd0..d7c60b34b14 100644 --- a/Core/src/Surfaces/LineSurface.cpp +++ b/Core/src/Surfaces/LineSurface.cpp @@ -18,6 +18,7 @@ #include "Acts/Surfaces/detail/AlignmentHelper.hpp" #include "Acts/Utilities/Helpers.hpp" #include "Acts/Utilities/Intersection.hpp" +#include "Acts/Utilities/JacobianHelpers.hpp" #include "Acts/Utilities/ThrowAssert.hpp" #include @@ -204,11 +205,6 @@ Acts::BoundToFreeMatrix Acts::LineSurface::boundToFreeJacobian( Vector3 position = freeParams.segment<3>(eFreePos0); // The direction Vector3 direction = freeParams.segment<3>(eFreeDir0); - // Get the sines and cosines directly - double cosTheta = std::cos(boundParams[eBoundTheta]); - double sinTheta = std::sin(boundParams[eBoundTheta]); - double cosPhi = std::cos(boundParams[eBoundPhi]); - double sinPhi = std::sin(boundParams[eBoundPhi]); // retrieve the reference frame auto rframe = referenceFrame(gctx, position, direction); @@ -220,11 +216,8 @@ Acts::BoundToFreeMatrix Acts::LineSurface::boundToFreeJacobian( // the time component jacToGlobal(eFreeTime, eBoundTime) = 1; // the momentum components - jacToGlobal(eFreeDir0, eBoundPhi) = -sinTheta * sinPhi; - jacToGlobal(eFreeDir0, eBoundTheta) = cosTheta * cosPhi; - jacToGlobal(eFreeDir1, eBoundPhi) = sinTheta * cosPhi; - jacToGlobal(eFreeDir1, eBoundTheta) = cosTheta * sinPhi; - jacToGlobal(eFreeDir2, eBoundTheta) = -sinTheta; + jacToGlobal.block<3, 2>(eFreeDir0, eBoundPhi) = + sphericalToFreeDirectionJacobian(direction); jacToGlobal(eFreeQOverP, eBoundQOverP) = 1; // For the derivative of global position with bound angles, refer the diff --git a/Core/src/Surfaces/Surface.cpp b/Core/src/Surfaces/Surface.cpp index 5fcee6b7fec..13064a506ca 100644 --- a/Core/src/Surfaces/Surface.cpp +++ b/Core/src/Surfaces/Surface.cpp @@ -12,6 +12,7 @@ #include "Acts/EventData/detail/TransformationBoundToFree.hpp" #include "Acts/Surfaces/SurfaceBounds.hpp" #include "Acts/Surfaces/detail/AlignmentHelper.hpp" +#include "Acts/Utilities/JacobianHelpers.hpp" #include "Acts/Utilities/VectorHelpers.hpp" #include @@ -260,9 +261,6 @@ Acts::BoundToFreeMatrix Acts::Surface::boundToFreeJacobian( const Vector3 position = freeParams.segment<3>(eFreePos0); // The direction const Vector3 direction = freeParams.segment<3>(eFreeDir0); - // Use fast evaluation function of sin/cos - auto [cosPhi, sinPhi, cosTheta, sinTheta, invSinTheta] = - VectorHelpers::evaluateTrigonomics(direction); // retrieve the reference frame const auto rframe = referenceFrame(gctx, position, direction); // Initialize the jacobian from local to global @@ -272,11 +270,8 @@ Acts::BoundToFreeMatrix Acts::Surface::boundToFreeJacobian( // the time component jacToGlobal(eFreeTime, eBoundTime) = 1; // the momentum components - jacToGlobal(eFreeDir0, eBoundPhi) = (-sinTheta) * sinPhi; - jacToGlobal(eFreeDir0, eBoundTheta) = cosTheta * cosPhi; - jacToGlobal(eFreeDir1, eBoundPhi) = sinTheta * cosPhi; - jacToGlobal(eFreeDir1, eBoundTheta) = cosTheta * sinPhi; - jacToGlobal(eFreeDir2, eBoundTheta) = (-sinTheta); + jacToGlobal.block<3, 2>(eFreeDir0, eBoundPhi) = + sphericalToFreeDirectionJacobian(direction); jacToGlobal(eFreeQOverP, eBoundQOverP) = 1; return jacToGlobal; } @@ -287,9 +282,6 @@ Acts::FreeToBoundMatrix Acts::Surface::freeToBoundJacobian( const auto position = parameters.segment<3>(eFreePos0); // The direction const auto direction = parameters.segment<3>(eFreeDir0); - // Use fast evaluation function of sin/cos - auto [cosPhi, sinPhi, cosTheta, sinTheta, invSinTheta] = - VectorHelpers::evaluateTrigonomics(direction); // The measurement frame of the surface RotationMatrix3 rframeT = referenceFrame(gctx, position, direction).transpose(); @@ -300,11 +292,8 @@ Acts::FreeToBoundMatrix Acts::Surface::freeToBoundJacobian( // Time component jacToLocal(eBoundTime, eFreeTime) = 1; // Directional and momentum elements for reference frame surface - jacToLocal(eBoundPhi, eFreeDir0) = -sinPhi * invSinTheta; - jacToLocal(eBoundPhi, eFreeDir1) = cosPhi * invSinTheta; - jacToLocal(eBoundTheta, eFreeDir0) = cosPhi * cosTheta; - jacToLocal(eBoundTheta, eFreeDir1) = sinPhi * cosTheta; - jacToLocal(eBoundTheta, eFreeDir2) = -sinTheta; + jacToLocal.block<2, 3>(eBoundPhi, eFreeDir0) = + freeToSphericalDirectionJacobian(direction); jacToLocal(eBoundQOverP, eFreeQOverP) = 1; return jacToLocal; } diff --git a/Examples/Python/tests/root_file_hashes.txt b/Examples/Python/tests/root_file_hashes.txt index 408cec7540e..a7038b2e630 100644 --- a/Examples/Python/tests/root_file_hashes.txt +++ b/Examples/Python/tests/root_file_hashes.txt @@ -17,21 +17,21 @@ test_itk_seeding__particles.root: cf270be1d1e2eb2e847a8206471fc8e90625453d2a3b5e test_itk_seeding__particles_simulation.root: 6bbbca85d0702e845cfdde71c4625600932c142a8e2442f06cf43d6104260723 test_propagation__propagation_steps.root: ba2e20d66f9f58850a9248bfaaafecbf0e8257020bb879b4572b783917f4d19b test_material_recording__geant4_material_tracks.root: de4e05a7e328eb237a111155382a3c3f8b17715b11c7c855e0fbc71c128c3fab -test_truth_tracking_kalman[generic-0.0]__trackstates_fitter.root: c4133e0b53334c75a6151f87a81e3e89fec35f8b5f2169f4b1c7550745bec780 +test_truth_tracking_kalman[generic-0.0]__trackstates_fitter.root: fee25bb8c40770d178ffedb77e914013e828031bc45e0c1760ff22174fe0b1e1 test_truth_tracking_kalman[generic-0.0]__tracksummary_fitter.root: 8ab6f20d72822bd5027781c46f982d6be8f8fdf7cfee5081a1e5139b805a2883 test_truth_tracking_kalman[generic-0.0]__performance_track_finder.root: ada9cda2ec3c648b144bdaa081d6eff923c80f3d484c4196006e847428cf67a8 test_truth_tracking_kalman[generic-1000.0]__trackstates_fitter.root: 15c07008eba724100448204a88d27372d44c1d6feb520e9ee07cd425765492a8 test_truth_tracking_kalman[generic-1000.0]__tracksummary_fitter.root: 381ffeb5e9b12fb8c12f0f510f4e637bca7a47d48e3b7a2d156dc8c8ab31957e test_truth_tracking_kalman[generic-1000.0]__performance_track_finder.root: ada9cda2ec3c648b144bdaa081d6eff923c80f3d484c4196006e847428cf67a8 -test_truth_tracking_kalman[odd-0.0]__trackstates_fitter.root: e41f58fc755e72ca17fa47577fd1b2ac9625126646f5c9449dc64c1a8deb9980 +test_truth_tracking_kalman[odd-0.0]__trackstates_fitter.root: 60ec798037cdbdeda32699ebc039be8b96dc0b358d299525e749c585c9e2ae0c test_truth_tracking_kalman[odd-0.0]__tracksummary_fitter.root: c4736f81d25f949f16062c4795eaaa3b1589e98c438272effef962948967e1ce test_truth_tracking_kalman[odd-0.0]__performance_track_finder.root: 39aec6316cceb90e314e16b02947faa691c18f57c3a851a25e547a8fc05a4593 test_truth_tracking_kalman[odd-1000.0]__trackstates_fitter.root: 93aa34d1346d87d9f53834ddc519077e47f15f0527cc56e9adf93382ba9d3c0d test_truth_tracking_kalman[odd-1000.0]__tracksummary_fitter.root: 99007dc5e8f596419ad428d0991dc7cdb99eed7dc81c05bd38c86606790cd45b test_truth_tracking_kalman[odd-1000.0]__performance_track_finder.root: 39aec6316cceb90e314e16b02947faa691c18f57c3a851a25e547a8fc05a4593 -test_truth_tracking_gsf[generic]__trackstates_gsf.root: 68981c5e5660f885b46e76e10fef07b20ccd2ef414c237e1f4a25a039886c610 +test_truth_tracking_gsf[generic]__trackstates_gsf.root: 477d806ed2e4a3367295165d5527980ff388273017cbda89706c619691358046 test_truth_tracking_gsf[generic]__tracksummary_gsf.root: 56f9a38033ef3d5c381074a245ea58c94c623ef11f9db60f7143b2cac3316b4a -test_truth_tracking_gsf[odd]__trackstates_gsf.root: c62a19d33a83ea0d87469ba1582e20f8c408700220d5df69092ce42557a437f5 +test_truth_tracking_gsf[odd]__trackstates_gsf.root: e900094f5889e933237d76c27e84bf216d7b771530ca4f836113c9bb32c55335 test_truth_tracking_gsf[odd]__tracksummary_gsf.root: bb40fd30b954fdf8fd1322dd5144b2708ee00aa2093a97d75e13d91c628f2401 test_particle_gun__particles.root: 7eec62018b6944fea565dad75aa41ef87d1f2737b2a814fbab189817ac8180fe test_material_mapping__material-map_tracks.root: 6e1441c418ff0b17983c2d0483248cc1dee6b77b09d0ca9d03c742c9d1373630 @@ -44,22 +44,22 @@ test_digitization_example_input[smeared]__particles.root: 7eec62018b6944fea565da test_digitization_example_input[smeared]__measurements.root: 0c168d371d0130c68d1ee44bd77eeeb3cf702a77c2afbf12bed8354b61a29262 test_digitization_example_input[geometric]__particles.root: 7eec62018b6944fea565dad75aa41ef87d1f2737b2a814fbab189817ac8180fe test_digitization_example_input[geometric]__measurements.root: 0c6d88b4de3ee7365103b8f0d6be6b4db3d7b7f2a59d3db58a1e5f89fa8130b3 -test_ckf_tracks_example[generic-full_seeding]__trackstates_ckf.root: f6dfb74b453adfbe079035da166d8c684ae292236d9837b959e59a5412de25ae -test_ckf_tracks_example[generic-full_seeding]__tracksummary_ckf.root: 627db7f14581706ff2dcc477f681cd925593a3f738b6a4235b0a5195631b3fa8 +test_ckf_tracks_example[generic-full_seeding]__trackstates_ckf.root: f7f98fd27b4513f531d8a3ef08d65737273cb7809fe95c5ec61548b434d5870d +test_ckf_tracks_example[generic-full_seeding]__tracksummary_ckf.root: cbcdec52dde2e30c76fb9f5b0623198b0fbf9b5d0cf67dbb829cfaa16936942e test_ckf_tracks_example[generic-full_seeding]__performance_seeding_trees.root: 0e0676ffafdb27112fbda50d1cf627859fa745760f98073261dcf6db3f2f991e -test_ckf_tracks_example[generic-truth_estimated]__trackstates_ckf.root: 2ea5d63221354d434eec9a49fa5a34523e45b98dc43bb76246e729dcc944d0ff +test_ckf_tracks_example[generic-truth_estimated]__trackstates_ckf.root: 3b1031c8420b0c9ecdff950636f6b246b8e12c1510974715403cac3cd5e3a0be test_ckf_tracks_example[generic-truth_estimated]__tracksummary_ckf.root: b9666af7fb37032a79b48d7ed87a21e3a797a12e7e55a62ff6c0e1402eef2d20 test_ckf_tracks_example[generic-truth_estimated]__performance_seeding.root: 1facb05c066221f6361b61f015cdf0918e94d9f3fce2269ec7b6a4dffeb2bc7e -test_ckf_tracks_example[generic-truth_smeared]__trackstates_ckf.root: c84e34a6a20381849e532412b1e6c30a9161edfecbd6c10531654a602b166471 +test_ckf_tracks_example[generic-truth_smeared]__trackstates_ckf.root: 24ce6e572c0682ae5ad3034816b58347d1578742a17fb2e5f7b7cec7a028fa0e test_ckf_tracks_example[generic-truth_smeared]__tracksummary_ckf.root: 3c7effc176ec1b10cd8f6fd3aceb6168343189a9f30f6c0c874b16c0f836d1d3 -test_ckf_tracks_example[odd-full_seeding]__trackstates_ckf.root: 15ac7c58345025e9123daae37b8ec966855641b74b11eb1475ca21a1399d4aae -test_ckf_tracks_example[odd-full_seeding]__tracksummary_ckf.root: 5fee89ae9b56f0816ce53a13bae1d63fa5044b56f90f125423b214b7389278a5 +test_ckf_tracks_example[odd-full_seeding]__trackstates_ckf.root: 3453df24d10b5cc1c85007456212b937ae1dc30292588e8a06794d9e27f24b6a +test_ckf_tracks_example[odd-full_seeding]__tracksummary_ckf.root: 8c0848b7bce23c76abb2ccb20d2a468c99f9106cf584fe75a58e2e29dca87930 test_ckf_tracks_example[odd-full_seeding]__performance_seeding_trees.root: 43c58577aafe07645e5660c4f43904efadf91d8cda45c5c04c248bbe0f59814f -test_ckf_tracks_example[odd-truth_estimated]__trackstates_ckf.root: 4c35dfb1ed48fc7ba23fb98901827009a07c5194da3967db2e5f1f19323e130f -test_ckf_tracks_example[odd-truth_estimated]__tracksummary_ckf.root: ab964bd60050eadc3ec2b6238ee5bb960cc6a0cd770830ffa065e7ccd21f92be +test_ckf_tracks_example[odd-truth_estimated]__trackstates_ckf.root: baa5b89bea7a26cac7f228f1f2c5b1bbd00410595be005b29737cbf2bf3a0fac +test_ckf_tracks_example[odd-truth_estimated]__tracksummary_ckf.root: 37ff338bc0c9a22eb5ee2afa3d0486c15ce1ab9d553b5704e2714ce505c43e73 test_ckf_tracks_example[odd-truth_estimated]__performance_seeding.root: 1a36b7017e59f1c08602ef3c2cb0483c51df248f112e3780c66594110719c575 -test_ckf_tracks_example[odd-truth_smeared]__trackstates_ckf.root: da11fb58148e81db3175253dab6657008ddee3b365e3b82fddcc90e78e358011 -test_ckf_tracks_example[odd-truth_smeared]__tracksummary_ckf.root: 9cc039f0aaba74f6474adcbefae1c276c3280a8bea142be26452ea07534d8e03 +test_ckf_tracks_example[odd-truth_smeared]__trackstates_ckf.root: 398fbbb721bc7b536158c8b83ca887e9f03d39f1421845827261616d9ab1eb94 +test_ckf_tracks_example[odd-truth_smeared]__tracksummary_ckf.root: da753244139c3b4111d271b4e80249d37dbc4c940b4cc0697b2739ebfbea308b test_vertex_fitting_reading[Truth-False-100]__performance_vertexing.root: 76ef6084d758dfdfc0151ddec2170e12d73394424e3dac4ffe46f0f339ec8293 test_vertex_fitting_reading[Iterative-False-100]__performance_vertexing.root: 60372210c830a04f95ceb78c6c68a9b0de217746ff59e8e73053750c837b57eb test_vertex_fitting_reading[Iterative-True-100]__performance_vertexing.root: e34f217d524a5051dbb04a811d3407df3ebe2cc4bb7f54f6bda0847dbd7b52c3