From 2c080a5acedefcd8aed42f17f0357f7e48d1224f Mon Sep 17 00:00:00 2001 From: andiwand Date: Mon, 11 Dec 2023 09:28:21 +0100 Subject: [PATCH 1/9] draft --- .../Acts/Propagator/CovarianceTransport.hpp | 120 -------- .../Propagator/detail/CovarianceEngine.hpp | 28 +- .../Acts/Propagator/detail/JacobianEngine.hpp | 82 ++--- .../Acts/Utilities/JacobianHelpers.hpp | 68 +++++ Core/include/Acts/Utilities/VectorHelpers.hpp | 6 +- Core/src/Propagator/CMakeLists.txt | 1 - Core/src/Propagator/CovarianceTransport.cpp | 120 -------- .../Propagator/detail/CovarianceEngine.cpp | 281 ++---------------- Core/src/Propagator/detail/JacobianEngine.cpp | 146 +++++---- Core/src/Surfaces/DiscSurface.cpp | 44 +-- Core/src/Surfaces/LineSurface.cpp | 13 +- Core/src/Surfaces/Surface.cpp | 21 +- .../Acts/Plugins/EDM4hep/EDM4hepUtil.hpp | 1 - Plugins/EDM4hep/src/EDM4hepUtil.cpp | 19 +- Tests/Benchmarks/CMakeLists.txt | 1 - .../CovarianceTransportBenchmark.cpp | 92 ------ .../UnitTests/Core/Propagator/CMakeLists.txt | 1 - .../Propagator/CovarianceTransportTests.cpp | 237 --------------- 18 files changed, 257 insertions(+), 1024 deletions(-) delete mode 100644 Core/include/Acts/Propagator/CovarianceTransport.hpp create mode 100644 Core/include/Acts/Utilities/JacobianHelpers.hpp delete mode 100644 Core/src/Propagator/CovarianceTransport.cpp delete mode 100644 Tests/Benchmarks/CovarianceTransportBenchmark.cpp delete mode 100644 Tests/UnitTests/Core/Propagator/CovarianceTransportTests.cpp diff --git a/Core/include/Acts/Propagator/CovarianceTransport.hpp b/Core/include/Acts/Propagator/CovarianceTransport.hpp deleted file mode 100644 index adb22669a68..00000000000 --- a/Core/include/Acts/Propagator/CovarianceTransport.hpp +++ /dev/null @@ -1,120 +0,0 @@ -// This file is part of the Acts project. -// -// Copyright (C) 2021 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/Definitions/TrackParametrization.hpp" -#include "Acts/Geometry/GeometryContext.hpp" -#include "Acts/Surfaces/Surface.hpp" - -#include -#include -#include -#include - -namespace Acts { - -using VariantCovariance = std::variant; - -using VariantTransportJacobian = - std::variant; - -/// Helper struct holding the necessary cache for the Covariance -/// transport in various parametrisations. -struct CovarianceCache { - /// Internal cache state, indicates correct setup - bool applyTransport = false; - - /// Internal cache state, at surface indication - const Surface* atSurface = nullptr; - - /// Internal cache state, at position - Vector3 atPosition = Vector3::Zero(); - - /// Variant: the currently held covariance - VariantCovariance covariance; - - /// Optional for starting from bound or curvilinear - std::optional boundToFreeJacobian = std::nullopt; - /// Options for starting from free - std::optional> anglesToDirectionJacobian = std::nullopt; - std::optional> directionToAnglesJacobian = std::nullopt; - - /// Non-variant: the free transport jacobian - FreeMatrix freeTransportJacobian = FreeMatrix::Identity(); - /// Non-variant: the free derivatives - FreeVector freeToPathDerivatives = FreeVector::Zero(); - - /// Defaulted constructor, gives invalid cache - CovarianceCache() = default; - - /// Constructor from bound & surface - /// - /// @param gctx The current geometry context - /// @param surface The surface of the bound representation - /// @param position The position of the representation - /// @param boundParameters The bound parameters at the surface - /// @param boundCovariance The bound covariance to be propagated - /// - /// This constructor will set the variant covariance type to - /// a bound matrix, remember the surface & establish the - /// jacobian between bound and free parametrisation. - CovarianceCache(const GeometryContext& gctx, const Surface& surface, - const Vector3& position, const BoundVector& boundParameters, - const BoundSquareMatrix& boundCovariance); - - /// Constructor from curvilinear - /// - /// @param position The position of the representation - /// @param direction The direction of at the representation - /// @param boundCovariance The bound covariance to be propagated - /// - /// This constructor will set the variant covariance type to - /// a bound matrix, remember the surface & establish the - /// jacobian between bound and free parametrisation. - CovarianceCache(const Vector3& position, const Vector3& direction, - const BoundSquareMatrix& boundCovariance); - - /// Construction from free - /// - /// @param freeParameters The free parameters - /// @param freeCovariance The free covariance to be propagated - /// - CovarianceCache(const FreeVector& freeParameters, - const FreeSquareMatrix& freeCovariance); -}; - -/// Transport the covariance to a new (surface) bound definition -/// -/// @param gctx [in] The current geometry context -/// @param surface [in] The surface of the bound representation -/// @param parameters [in] The free parametrisation on the surface -/// @param cCache [in,out] the covariance cache (to be modified) -/// -/// @return a variant transport jacobian -std::tuple -transportCovarianceToBound(const GeometryContext& gctx, const Surface& surface, - const FreeVector& parameters, - CovarianceCache& cCache); - -/// Transport the covariance to a new (curvilinear) bound definition -/// -/// @param direction [in] The track direction at the new curvilinear definition -/// @param cCache [in,out] the covariance cache (to be modified) -std::tuple -transportCovarianceToCurvilinear(const Vector3& direction, - CovarianceCache& cCache); - -/// Transport the covariance to a new free definition -/// -/// @param cCache [in,out] the covariance cache (to be modified) -std::tuple -transportCovarianceToFree(CovarianceCache& cCache); - -} // namespace Acts diff --git a/Core/include/Acts/Propagator/detail/CovarianceEngine.hpp b/Core/include/Acts/Propagator/detail/CovarianceEngine.hpp index e5b61043844..eb530fd5e5c 100644 --- a/Core/include/Acts/Propagator/detail/CovarianceEngine.hpp +++ b/Core/include/Acts/Propagator/detail/CovarianceEngine.hpp @@ -25,13 +25,13 @@ #include namespace Acts { +namespace detail { /// @brief These functions perform the transport of a covariance matrix using /// given Jacobians. The required data is provided by the stepper object /// with some additional data. Since this is a purely algebraic problem the /// calculations are identical for @c StraightLineStepper and @c EigenStepper. /// As a consequence the methods can be located in a separate file. -namespace detail { /// Create and return the bound state at the current position /// @@ -43,13 +43,13 @@ namespace detail { /// @param [in, out] jacobian Full jacobian since the last reset /// @param [in, out] transportJacobian Global jacobian since the last reset /// @param [in, out] derivatives Path length derivatives of the free, nominal -/// parameters +/// parameters /// @param [in, out] jacToGlobal Projection jacobian of the last bound -/// parametrisation to free parameters +/// parametrisation to free parameters /// @param [in, out] parameters Free, nominal parametrisation /// @param [in] particleHypothesis Particle hypothesis /// @param [in] covTransport Decision whether the covariance transport should be -/// performed +/// performed /// @param [in] accumulatedPath Propagated distance /// @param [in] surface Target surface on which the state is represented /// @param [in] freeToBoundCorrection Correction for non-linearity effect during transform from free to bound @@ -64,8 +64,7 @@ Result> boundState( FreeVector& derivatives, BoundToFreeMatrix& jacToGlobal, FreeVector& parameters, const ParticleHypothesis& particleHypothesis, bool covTransport, double accumulatedPath, const Surface& surface, - const FreeToBoundCorrection& freeToBoundCorrection = - FreeToBoundCorrection(false)); + const FreeToBoundCorrection& freeToBoundCorrection); /// Create and return a curvilinear state at the current position /// @@ -75,13 +74,13 @@ Result> boundState( /// @param [in, out] jacobian Full jacobian since the last reset /// @param [in, out] transportJacobian Global jacobian since the last reset /// @param [in, out] derivatives Path length derivatives of the free, nominal -/// parameters +/// parameters /// @param [in, out] jacToGlobal Projection jacobian of the last bound -/// parametrisation to free parameters +/// parametrisation to free parameters /// @param [in] parameters Free, nominal parametrisation /// @param [in] particleHypothesis Particle hypothesis /// @param [in] covTransport Decision whether the covariance transport should be -/// performed +/// performed /// @param [in] accumulatedPath Propagated distance /// /// @return A curvilinear state: @@ -96,7 +95,7 @@ std::tuple curvilinearState( double accumulatedPath); /// @brief Method for on-demand covariance transport of a bound/curvilinear to -/// another bound representation. +/// another bound representation. /// /// @param [in] geoContext The geometry context /// @param [in, out] boundCovariance The covariance matrix of the state @@ -104,7 +103,7 @@ std::tuple curvilinearState( /// @param [in, out] freeTransportJacobian Global jacobian since the last reset /// @param [in, out] freeToPathDerivatives Path length derivatives /// @param [in, out] boundToFreeJacobian Projection jacobian of the last bound -/// parametrisation to free parameters +/// parametrisation to free parameters /// @param [in, out] freeParameters Free, nominal parametrisation /// @param [in] surface is the surface to which the covariance is /// forwarded to @@ -117,18 +116,17 @@ void transportCovarianceToBound( BoundMatrix& fullTransportJacobian, FreeMatrix& freeTransportJacobian, FreeVector& freeToPathDerivatives, BoundToFreeMatrix& boundToFreeJacobian, FreeVector& freeParameters, const Surface& surface, - const FreeToBoundCorrection& freeToBoundCorrection = - FreeToBoundCorrection(false)); + const FreeToBoundCorrection& freeToBoundCorrection); /// @brief Method for on-demand covariance transport of a bound/curvilinear -/// to a new curvilinear representation. +/// to a new curvilinear representation. /// /// @param [in, out] boundCovariance The covariance matrix of the state /// @param [in, out] fullTransportJacobian Full jacobian since the last reset /// @param [in, out] freeTransportJacobian Global jacobian since the last reset /// @param [in, out] freeToPathDerivatives Path length derivatives /// @param [in, out] boundToFreeJacobian Projection jacobian of the last bound -/// parametrisation to free parameters +/// parametrisation to free parameters /// @param [in] direction Normalised direction vector /// void transportCovarianceToCurvilinear(BoundSquareMatrix& boundCovariance, diff --git a/Core/include/Acts/Propagator/detail/JacobianEngine.hpp b/Core/include/Acts/Propagator/detail/JacobianEngine.hpp index a3d38240690..1f5bc8da6a6 100644 --- a/Core/include/Acts/Propagator/detail/JacobianEngine.hpp +++ b/Core/include/Acts/Propagator/detail/JacobianEngine.hpp @@ -12,20 +12,19 @@ #include "Acts/Utilities/detail/ReferenceWrapperAnyCompat.hpp" #include "Acts/Definitions/Algebra.hpp" -#include "Acts/EventData/TrackParameters.hpp" #include "Acts/Geometry/GeometryContext.hpp" #include "Acts/Surfaces/Surface.hpp" namespace Acts { +namespace detail { /// @brief These functions perform the calculation of the Jacobians for the /// the covariance transport. This is a purely algebraic problem the /// calculations are identical for @c StraightLineStepper and @c EigenStepper. /// As a consequence the methods can be located in a separate file. -namespace detail { /// @brief Evaluate the projection Jacobian from free to curvilinear parameters -/// without transport jacobian. +/// without transport jacobian. /// /// @param [in] direction Normalised direction vector /// @@ -33,31 +32,15 @@ namespace detail { FreeToBoundMatrix freeToCurvilinearJacobian(const Vector3& direction); /// @brief Evaluate the projection Jacobian from curvilinear to free parameters -/// without transport jacobian. +/// without transport jacobian. /// /// @param [in] direction Normalised direction vector /// /// @return Projection Jacobian BoundToFreeMatrix curvilinearToFreeJacobian(const Vector3& direction); -/// @brief This function calculates the jacobian from a free parameterisation -/// to a mixed one with angular representation. -/// -/// @param direction The direction of the track -/// -/// @return the 7x8 jacobian for direction to angles relation -ActsMatrix<7, 8> directionToAnglesJacobian(const Vector3& direction); - -/// @brief This function calculates the jacobian from a free parameterisation -/// to a mixed one with angular representation. -/// -/// @param direction The direction of the track -/// -/// @return the 8x7 jacobian for angles to direction relation -ActsMatrix<8, 7> anglesToDirectionJacobian(const Vector3& direction); - /// @brief This function calculates the full transport jacobian from a bound -/// curvilinear representation to a new bound representation +/// curvilinear representation to a new bound representation /// /// @note Modifications of the jacobian related to the /// projection onto a surface is considered. Since a variation of the start @@ -97,7 +80,7 @@ BoundMatrix boundToBoundTransportJacobian( /// @param [in] boundToFreeJacobian Jacobian from bound to free at start /// @param [in] freeTransportJacobian Transport jacobian free to free /// @param [in] freeToPathDerivatives Path length derivatives for free -/// parameters +/// parameters /// /// @note The parameter @p surface is only required if projected to bound /// parameters. In the case of curvilinear parameters the geometry and the @@ -110,8 +93,7 @@ BoundMatrix boundToCurvilinearTransportJacobian( const FreeVector& freeToPathDerivatives); /// @brief This function calculates the full jacobian from a given -/// bound/curvilinear parameterisation from a new free -/// parameterisation. +/// bound/curvilinear parameterisation from a new free parameterisation. /// /// @param [in] boundToFreeJacobian Jacobian from bound to free at start /// @param [in] freeTransportJacobian Transport jacobian free to free @@ -122,8 +104,7 @@ BoundToFreeMatrix boundToFreeTransportJacobian( const FreeMatrix& freeTransportJacobian); /// @brief This function calculates the full jacobian from a given -/// free parameterisation to a new curvilinear bound -/// parameterisation. +/// free parameterisation to a new curvilinear bound parameterisation. /// /// @note Modifications of the jacobian related to the /// projection onto the target surface is considered. Since a variation of @@ -134,12 +115,12 @@ BoundToFreeMatrix boundToFreeTransportJacobian( /// @param [in] geoContext The geometry Context /// @param [in] freeParameters Free, nominal parametrisation /// @param [in] directionToAnglesJacobian The relation jacobian from dir to -/// angle +/// angle /// @param [in] anglesToDirectionJacobian The relation jacobian from angle to -/// dir +/// dir /// @param [in] freeTransportJacobian Transport jacobian free to free /// @param [in] freeToPathDerivatives Path length derivatives for free -/// parameters +/// parameters /// @param [in] surface The target surface /// /// @return the 6x8 transport jacobian from bound to free @@ -158,9 +139,9 @@ FreeToBoundMatrix freeToBoundTransportJacobian( /// /// @param [in] direction Normalised direction vector /// @param [in] directionToAnglesJacobian The relation jacobian from dir to -/// angle +/// angle /// @param [in] anglesToDirectionJacobian The relation jacobian from angle to -/// dir +/// dir /// @param [in] freeTransportJacobian Transport jacobian free to free /// @param [in] freeToPathDerivatives Path length derivatives for free /// @@ -174,9 +155,9 @@ FreeToBoundMatrix freeToCurvilinearTransportJacobian( /// @brief This function calculates the free transfport jacobian from a free /// parameterisation. /// @param [in] directionToAnglesJacobian The relation jacobian from dir to -/// angle +/// angle /// @param [in] anglesToDirectionJacobian The relation jacobian from angle to -/// dir +/// dir /// @param [in] freeTransportJacobian Transport jacobian free to free /// /// @return a 8x8 transport jacobian from free to free @@ -185,6 +166,39 @@ FreeMatrix freeToFreeTransportJacobian( const ActsMatrix<8, 7>& anglesToDirectionJacobian, const FreeMatrix& freeTransportJacobian); -} // namespace detail +/// @brief This function reinitialises the state members required for the +/// covariance transport +/// +/// @param [in] geoContext The geometry context +/// @param [in, out] freeTransportJacobian The transport jacobian from start +/// free to final free parameters +/// @param [in, out] freeToPathDerivatives Path length derivatives of the free, +/// nominal parameters +/// @param [in, out] boundToFreeJacobian Projection jacobian of the last bound +/// parametrisation to free parameters +/// @param [in] freeParameters Free, nominal parametrisation +/// @param [in] surface The reference surface of the local parametrisation +Result reinitializeJacobians(const GeometryContext& geoContext, + FreeMatrix& freeTransportJacobian, + FreeVector& freeToPathDerivatives, + BoundToFreeMatrix& boundToFreeJacobian, + const FreeVector& freeParameters, + const Surface& surface); + +/// @brief This function reinitialises the state members required for the +/// covariance transport +/// +/// @param [in, out] freeTransportJacobian The transport jacobian from start +/// free to final free parameters +/// @param [in, out] derivatives Path length derivatives of the free, nominal +/// parameters +/// @param [in, out] boundToFreeJacobian Projection jacobian of the last bound +/// parametrisation to free parameters +/// @param [in] direction Normalised direction vector +void reinitializeJacobians(FreeMatrix& freeTransportJacobian, + FreeVector& freeToPathDerivatives, + BoundToFreeMatrix& boundToFreeJacobian, + const Vector3& direction); +} // namespace detail } // namespace Acts 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/CMakeLists.txt b/Core/src/Propagator/CMakeLists.txt index df3d22a2251..4e1f1d8aca0 100644 --- a/Core/src/Propagator/CMakeLists.txt +++ b/Core/src/Propagator/CMakeLists.txt @@ -1,7 +1,6 @@ target_sources( ActsCore PRIVATE - CovarianceTransport.cpp EigenStepperError.cpp MultiStepperError.cpp PropagatorError.cpp diff --git a/Core/src/Propagator/CovarianceTransport.cpp b/Core/src/Propagator/CovarianceTransport.cpp deleted file mode 100644 index 92cd2d6054a..00000000000 --- a/Core/src/Propagator/CovarianceTransport.cpp +++ /dev/null @@ -1,120 +0,0 @@ -// This file is part of the Acts project. -// -// Copyright (C) 2021 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/. - -#include "Acts/Propagator/CovarianceTransport.hpp" - -#include "Acts/Propagator/detail/JacobianEngine.hpp" - -Acts::CovarianceCache::CovarianceCache(const GeometryContext& gctx, - const Surface& surface, - const Vector3& position, - const BoundVector& boundParameters, - const BoundSquareMatrix& boundCovariance) - : applyTransport(true), atSurface(&surface), atPosition(position) { - covariance.template emplace(boundCovariance); - boundToFreeJacobian = surface.boundToFreeJacobian(gctx, boundParameters); -} - -Acts::CovarianceCache::CovarianceCache(const Vector3& position, - const Vector3& direction, - const BoundSquareMatrix& boundCovariance) - : applyTransport(true), atPosition(position) { - covariance.emplace(boundCovariance); - boundToFreeJacobian = detail::curvilinearToFreeJacobian(direction); -} - -Acts::CovarianceCache::CovarianceCache(const FreeVector& freeParameters, - const FreeSquareMatrix& freeCovariance) - : applyTransport(true), atPosition(freeParameters.segment<3>(eFreePos0)) { - covariance.emplace(freeCovariance); - anglesToDirectionJacobian = - detail::anglesToDirectionJacobian(freeParameters.segment<3>(eFreeDir0)); - directionToAnglesJacobian = - detail::directionToAnglesJacobian(freeParameters.segment<3>(eFreeDir0)); -} - -std::tuple -Acts::transportCovarianceToBound(const GeometryContext& gctx, - const Surface& surface, - const FreeVector& freeParameters, - CovarianceCache& cCache) { - if (cCache.boundToFreeJacobian.has_value()) { - // Create the full transport jacobian: bound/curvilinear to bound - const auto ftJacobian = detail::boundToBoundTransportJacobian( - gctx, freeParameters, cCache.boundToFreeJacobian.value(), - cCache.freeTransportJacobian, cCache.freeToPathDerivatives, surface); - // Perform the transport - const auto& covariance = std::get(cCache.covariance); - BoundSquareMatrix newCovariance = - ftJacobian * covariance * ftJacobian.transpose(); - return {newCovariance, ftJacobian}; - } else { - // Create the full transport jacobian: free to bound - const auto ftJacobian = detail::freeToBoundTransportJacobian( - gctx, freeParameters, cCache.directionToAnglesJacobian.value(), - cCache.anglesToDirectionJacobian.value(), cCache.freeTransportJacobian, - cCache.freeToPathDerivatives, surface); - // Perform the transport - const auto& covariance = std::get(cCache.covariance); - BoundSquareMatrix newCovariance = - ftJacobian * covariance * ftJacobian.transpose(); - return {newCovariance, ftJacobian}; - } -} - -std::tuple -Acts::transportCovarianceToCurvilinear(const Vector3& direction, - CovarianceCache& cCache) { - if (cCache.boundToFreeJacobian.has_value()) { - // Create the full transport jacobian: bound/curvilinear to - // curvilinear - auto ftJacobian = detail::boundToCurvilinearTransportJacobian( - direction, cCache.boundToFreeJacobian.value(), - cCache.freeTransportJacobian, cCache.freeToPathDerivatives); - // Perform the transport - const auto& covariance = std::get(cCache.covariance); - BoundSquareMatrix newCovariance = - ftJacobian * covariance * ftJacobian.transpose(); - return {newCovariance, ftJacobian}; - } else { - // Create the full transport jacobian: free to curvilinear - const auto ftJacobian = detail::freeToCurvilinearTransportJacobian( - direction, cCache.directionToAnglesJacobian.value(), - cCache.anglesToDirectionJacobian.value(), cCache.freeTransportJacobian, - cCache.freeToPathDerivatives); - // Perform the transport - const auto& covariance = std::get(cCache.covariance); - BoundSquareMatrix newCovariance = - ftJacobian * covariance * ftJacobian.transpose(); - return {newCovariance, ftJacobian}; - } -} - -std::tuple -Acts::transportCovarianceToFree(CovarianceCache& cCache) { - if (cCache.boundToFreeJacobian.has_value()) { - // Create the full transport jacobian: bound/curvilinear to free - const auto ftJacobian = detail::boundToFreeTransportJacobian( - cCache.boundToFreeJacobian.value(), cCache.freeTransportJacobian); - // Perform the transport - const auto& covariance = std::get(cCache.covariance); - FreeSquareMatrix newCovariance = - ftJacobian * covariance * ftJacobian.transpose(); - return {newCovariance, ftJacobian}; - } else { - // Create the full transport jacobian: free to free - const auto ftJacobian = detail::freeToFreeTransportJacobian( - cCache.directionToAnglesJacobian.value(), - cCache.anglesToDirectionJacobian.value(), cCache.freeTransportJacobian); - // Perform the transport - const auto& covariance = std::get(cCache.covariance); - FreeSquareMatrix newCovariance = - ftJacobian * covariance * ftJacobian.transpose(); - return {newCovariance, ftJacobian}; - } -} diff --git a/Core/src/Propagator/detail/CovarianceEngine.cpp b/Core/src/Propagator/detail/CovarianceEngine.cpp index 2d65323d0bc..575b093cc2e 100644 --- a/Core/src/Propagator/detail/CovarianceEngine.cpp +++ b/Core/src/Propagator/detail/CovarianceEngine.cpp @@ -9,262 +9,28 @@ #include "Acts/Propagator/detail/CovarianceEngine.hpp" #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" #include "Acts/EventData/detail/TransformationBoundToFree.hpp" #include "Acts/EventData/detail/TransformationFreeToBound.hpp" -#include "Acts/Utilities/AlgebraHelpers.hpp" +#include "Acts/Propagator/detail/JacobianEngine.hpp" #include "Acts/Utilities/Result.hpp" -#include -#include #include #include -#include #include namespace Acts { -namespace { + /// Some type defs using Jacobian = BoundMatrix; - using BoundState = std::tuple; using CurvilinearState = std::tuple; -/// @brief Evaluate the projection Jacobian from free to curvilinear parameters -/// -/// @param [in] direction Normalised direction vector -/// -/// @return Projection Jacobian -FreeToBoundMatrix freeToCurvilinearJacobian(const Vector3& direction) { - // 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; - // prepare the jacobian to curvilinear - FreeToBoundMatrix jacToCurv = FreeToBoundMatrix::Zero(); - if (std::abs(cosTheta) < s_curvilinearProjTolerance) { - // 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; - } 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; - } - // Time parameter - jacToCurv(5, 3) = 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.; - - return jacToCurv; -} - -/// @brief This function calculates the full jacobian from local parameters at -/// the start surface to bound parameters at the final surface -/// -/// @note Modifications of the jacobian related to the -/// projection onto a surface is considered. Since a variation of the start -/// parameters within a given uncertainty would lead to a variation of the end -/// parameters, these need to be propagated onto the target surface. This an -/// approximated approach to treat the (assumed) small change. -/// -/// @param [in] geoContext The geometry Context -/// @param [in] freeParameters Free, nominal parametrisation -/// @param [in] boundToFreeJacobian The projection jacobian from start local -/// to start free parameters -/// @param [in] freeTransportJacobian The transport jacobian from start free to -/// final free parameters -/// @param [in] freeToPathDerivatives Path length derivatives of the final free -/// parameters -/// @param [in, out] fullTransportJacobian The full jacobian from start local to -/// bound parameters at the final surface -/// @param [in] surface The final surface onto which the projection should be -/// performed -void boundToBoundJacobian(const GeometryContext& geoContext, - const FreeVector& freeParameters, - const BoundToFreeMatrix& boundToFreeJacobian, - const FreeMatrix& freeTransportJacobian, - const FreeVector& freeToPathDerivatives, - BoundMatrix& fullTransportJacobian, - const Surface& surface) { - // Calculate the derivative of path length at the final surface or the - // point-of-closest approach w.r.t. free parameters - const FreeToPathMatrix freeToPath = - surface.freeToPathDerivative(geoContext, freeParameters); - // Calculate the jacobian from free to bound at the final surface - FreeToBoundMatrix freeToBoundJacobian = - surface.freeToBoundJacobian(geoContext, freeParameters); - // Calculate the full jacobian from the local/bound parameters at the start - // surface to local/bound parameters at the final surface - // @note jac(locA->locB) = jac(gloB->locB)*(1+ - // pathCorrectionFactor(gloB))*jacTransport(gloA->gloB) *jac(locA->gloA) - fullTransportJacobian = - freeToBoundJacobian * - (FreeMatrix::Identity() + freeToPathDerivatives * freeToPath) * - freeTransportJacobian * boundToFreeJacobian; -} - -/// @brief This function calculates the full jacobian from local parameters at -/// the start surface to final curvilinear parameters -/// -/// @note Modifications of the jacobian related to the -/// projection onto a curvilinear surface is considered. Since a variation of -/// the start parameters within a given uncertainty would lead to a variation of -/// the end parameters, these need to be propagated onto the target surface. -/// This is an approximated approach to treat the (assumed) small change. -/// -/// @param [in] direction Normalised direction vector -/// @param [in] boundToFreeJacobian The projection jacobian from local start -/// to global final parameters -/// @param [in] freeTransportJacobian The transport jacobian from start free to -/// final free parameters -/// @param [in] freeToPathDerivatives Path length derivatives of the final free -/// parameters -/// @param [in, out] jacFull The full jacobian from start local to curvilinear -/// parameters -/// -/// @note The parameter @p surface is only required if projected to bound -/// parameters. In the case of curvilinear parameters the geometry and the -/// position is known and the calculation can be simplified -void boundToCurvilinearJacobian(const Vector3& direction, - const BoundToFreeMatrix& boundToFreeJacobian, - const FreeMatrix& freeTransportJacobian, - const FreeVector& freeToPathDerivatives, - BoundMatrix& fullTransportJacobian) { - // Calculate the jacobian from global to local at the curvilinear surface - FreeToBoundMatrix freeToBoundJacobian = freeToCurvilinearJacobian(direction); - - // Update the jacobian to include the derivative of the path length at the - // curvilinear surface w.r.t. the free parameters - freeToBoundJacobian.topLeftCorner<6, 3>() += - (freeToBoundJacobian * freeToPathDerivatives) * - (-1.0 * direction).transpose(); - - // Calculate the full jocobian from the local parameters at the start surface - // to curvilinear parameters - // @note jac(locA->locB) = jac(gloB->locB)*(1+ - // pathCorrectionFactor(gloB))*jacTransport(gloA->gloB) *jac(locA->gloA) - fullTransportJacobian = - blockedMult(freeToBoundJacobian, - blockedMult(freeTransportJacobian, boundToFreeJacobian)); -} - -/// @brief This function reinitialises the state members required for the -/// covariance transport -/// -/// @param [in] geoContext The geometry context -/// @param [in, out] freeTransportJacobian The transport jacobian from start -/// free to final free parameters -/// @param [in, out] freeToPathDerivatives Path length derivatives of the free, -/// nominal parameters -/// @param [in, out] boundToFreeJacobian Projection jacobian of the last bound -/// parametrisation to free parameters -/// @param [in] freeParameters Free, nominal parametrisation -/// @param [in] surface The reference surface of the local parametrisation -Result reinitializeJacobians(const GeometryContext& geoContext, - FreeMatrix& freeTransportJacobian, - FreeVector& freeToPathDerivatives, - BoundToFreeMatrix& boundToFreeJacobian, - const FreeVector& freeParameters, - const Surface& surface) { - using VectorHelpers::phi; - using VectorHelpers::theta; - - // Reset the jacobians - freeTransportJacobian = FreeMatrix::Identity(); - freeToPathDerivatives = FreeVector::Zero(); - - // Get the local position - const Vector3 position = freeParameters.segment<3>(eFreePos0); - const Vector3 direction = freeParameters.segment<3>(eFreeDir0); - auto lpResult = surface.globalToLocal(geoContext, position, direction); - if (!lpResult.ok()) { - return lpResult.error(); - } - // Transform from free to bound parameters - Result boundParameters = detail::transformFreeToBoundParameters( - freeParameters, surface, geoContext); - if (!boundParameters.ok()) { - return boundParameters.error(); - } - // Reset the jacobian from local to global - boundToFreeJacobian = - surface.boundToFreeJacobian(geoContext, *boundParameters); - return Result::success(); -} - -/// @brief This function reinitialises the state members required for the -/// covariance transport -/// -/// @param [in, out] freeTransportJacobian The transport jacobian from start -/// free to final free parameters -/// @param [in, out] derivatives Path length derivatives of the free, nominal -/// parameters -/// @param [in, out] boundToFreeJacobian Projection jacobian of the last bound -/// parametrisation to free parameters -/// @param [in] direction Normalised direction vector -void reinitializeJacobians(FreeMatrix& freeTransportJacobian, - FreeVector& freeToPathDerivatives, - BoundToFreeMatrix& boundToFreeJacobian, - const Vector3& direction) { - // Reset the jacobians - freeTransportJacobian = FreeMatrix::Identity(); - 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; - - boundToFreeJacobian(eFreePos0, eBoundLoc0) = -sinPhi; - boundToFreeJacobian(eFreePos0, eBoundLoc1) = -cosPhi * cosTheta; - boundToFreeJacobian(eFreePos1, eBoundLoc0) = cosPhi; - 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(eFreeQOverP, eBoundQOverP) = 1; -} -} // namespace - -namespace detail { - -Result boundState( +Result detail::boundState( const GeometryContext& geoContext, BoundSquareMatrix& covarianceMatrix, BoundMatrix& jacobian, FreeMatrix& transportJacobian, FreeVector& derivatives, BoundToFreeMatrix& jacToGlobal, @@ -300,14 +66,12 @@ Result boundState( jacobian, accumulatedPath); } -CurvilinearState curvilinearState(BoundSquareMatrix& covarianceMatrix, - BoundMatrix& jacobian, - FreeMatrix& transportJacobian, - FreeVector& derivatives, - BoundToFreeMatrix& jacToGlobal, - const FreeVector& parameters, - const ParticleHypothesis& particleHypothesis, - bool covTransport, double accumulatedPath) { +CurvilinearState detail::curvilinearState( + BoundSquareMatrix& covarianceMatrix, BoundMatrix& jacobian, + FreeMatrix& transportJacobian, FreeVector& derivatives, + BoundToFreeMatrix& jacToGlobal, const FreeVector& parameters, + const ParticleHypothesis& particleHypothesis, bool covTransport, + double accumulatedPath) { const Vector3& direction = parameters.segment<3>(eFreeDir0); // Covariance transport @@ -340,7 +104,7 @@ CurvilinearState curvilinearState(BoundSquareMatrix& covarianceMatrix, accumulatedPath); } -void transportCovarianceToBound( +void detail::transportCovarianceToBound( const GeometryContext& geoContext, BoundSquareMatrix& boundCovariance, BoundMatrix& fullTransportJacobian, FreeMatrix& freeTransportJacobian, FreeVector& freeToPathDerivatives, BoundToFreeMatrix& boundToFreeJacobian, @@ -348,9 +112,9 @@ void transportCovarianceToBound( const FreeToBoundCorrection& freeToBoundCorrection) { // Calculate the full jacobian from local parameters at the start surface to // current bound parameters - boundToBoundJacobian(geoContext, freeParameters, boundToFreeJacobian, - freeTransportJacobian, freeToPathDerivatives, - fullTransportJacobian, surface); + fullTransportJacobian = boundToBoundTransportJacobian( + geoContext, freeParameters, boundToFreeJacobian, freeTransportJacobian, + freeToPathDerivatives, surface); bool correction = false; if (freeToBoundCorrection) { @@ -395,17 +159,15 @@ void transportCovarianceToBound( freeParameters, surface); } -void transportCovarianceToCurvilinear(BoundSquareMatrix& boundCovariance, - BoundMatrix& fullTransportJacobian, - FreeMatrix& freeTransportJacobian, - FreeVector& freeToPathDerivatives, - BoundToFreeMatrix& boundToFreeJacobian, - const Vector3& direction) { +void detail::transportCovarianceToCurvilinear( + BoundSquareMatrix& boundCovariance, BoundMatrix& fullTransportJacobian, + FreeMatrix& freeTransportJacobian, FreeVector& freeToPathDerivatives, + BoundToFreeMatrix& boundToFreeJacobian, const Vector3& direction) { // Calculate the full jacobian from local parameters at the start surface to // current curvilinear parameters - boundToCurvilinearJacobian(direction, boundToFreeJacobian, - freeTransportJacobian, freeToPathDerivatives, - fullTransportJacobian); + fullTransportJacobian = boundToCurvilinearTransportJacobian( + direction, boundToFreeJacobian, freeTransportJacobian, + freeToPathDerivatives); // Apply the actual covariance transport to get covariance of the current // curvilinear parameters @@ -421,5 +183,4 @@ void transportCovarianceToCurvilinear(BoundSquareMatrix& boundCovariance, boundToFreeJacobian, direction); } -} // namespace detail } // namespace Acts diff --git a/Core/src/Propagator/detail/JacobianEngine.cpp b/Core/src/Propagator/detail/JacobianEngine.cpp index e215f58fe4f..a84b7488eec 100644 --- a/Core/src/Propagator/detail/JacobianEngine.cpp +++ b/Core/src/Propagator/detail/JacobianEngine.cpp @@ -6,24 +6,25 @@ // 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/. +#include "Acts/Propagator/detail/JacobianEngine.hpp" + #include "Acts/Definitions/Algebra.hpp" #include "Acts/Definitions/Tolerance.hpp" #include "Acts/Definitions/TrackParametrization.hpp" #include "Acts/Geometry/GeometryContext.hpp" #include "Acts/Surfaces/Surface.hpp" +#include "Acts/Utilities/AlgebraHelpers.hpp" +#include "Acts/Utilities/JacobianHelpers.hpp" #include "Acts/Utilities/VectorHelpers.hpp" -#include #include -#include namespace Acts { -namespace detail { - -FreeToBoundMatrix freeToCurvilinearJacobian(const Vector3& direction) { - auto [cosPhi, sinPhi, cosTheta, sinTheta, invSinTheta] = +FreeToBoundMatrix detail::freeToCurvilinearJacobian(const Vector3& direction) { + 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) { @@ -60,8 +61,8 @@ FreeToBoundMatrix freeToCurvilinearJacobian(const Vector3& direction) { return freeToCurvJacobian; } -BoundToFreeMatrix curvilinearToFreeJacobian(const Vector3& direction) { - auto [cosPhi, sinPhi, cosTheta, sinTheta, invSinTheta] = +BoundToFreeMatrix detail::curvilinearToFreeJacobian(const Vector3& direction) { + auto [cosPhi, sinPhi, cosTheta, sinTheta] = VectorHelpers::evaluateTrigonomics(direction); // Prepare the jacobian to free @@ -85,48 +86,7 @@ BoundToFreeMatrix curvilinearToFreeJacobian(const Vector3& direction) { return curvToFreeJacobian; } -ActsMatrix<8, 7> anglesToDirectionJacobian(const Vector3& direction) { - ActsMatrix<8, 7> jacobian = ActsMatrix<8, 7>::Zero(); - - auto [cosPhi, sinPhi, cosTheta, sinTheta, invSinTheta] = - VectorHelpers::evaluateTrigonomics(direction); - - jacobian(0, 0) = 1.; - jacobian(1, 1) = 1.; - jacobian(2, 2) = 1.; - jacobian(3, 3) = 1.; - jacobian(7, 6) = 1.; - - jacobian(4, 4) = -sinTheta * sinPhi; - jacobian(4, 5) = cosTheta * cosPhi; - jacobian(5, 4) = sinTheta * cosPhi; - jacobian(5, 5) = cosTheta * sinPhi; - jacobian(6, 5) = -sinTheta; - return jacobian; -} - -ActsMatrix<7, 8> directionToAnglesJacobian(const Vector3& direction) { - ActsMatrix<7, 8> jacobian = ActsMatrix<7, 8>::Zero(); - - auto [cosPhi, sinPhi, cosTheta, sinTheta, invSinTheta] = - VectorHelpers::evaluateTrigonomics(direction); - - jacobian(0, 0) = 1.; - jacobian(1, 1) = 1.; - jacobian(2, 2) = 1.; - jacobian(3, 3) = 1.; - jacobian(6, 7) = 1.; - - jacobian(4, 4) = -sinPhi * invSinTheta; - jacobian(4, 5) = cosPhi * invSinTheta; - jacobian(5, 4) = cosPhi * cosTheta; - jacobian(5, 5) = sinPhi * cosTheta; - jacobian(5, 6) = -sinTheta; - - return jacobian; -} - -BoundMatrix boundToBoundTransportJacobian( +BoundMatrix detail::boundToBoundTransportJacobian( const GeometryContext& geoContext, const FreeVector& freeParameters, const BoundToFreeMatrix& boundToFreeJacobian, const FreeMatrix& freeTransportJacobian, @@ -141,29 +101,35 @@ BoundMatrix boundToBoundTransportJacobian( // https://acts.readthedocs.io/en/latest/white_papers/correction-for-transport-jacobian.html // Calculate the full jacobian from the local/bound parameters at the start // surface to local/bound parameters at the final surface + // @note jac(locA->locB) = jac(gloB->locB)*(1+ + // pathCorrectionFactor(gloB))*jacTransport(gloA->gloB) *jac(locA->gloA) return freeToBoundJacobian * (FreeMatrix::Identity() + freeToPathDerivatives * freeToPath) * freeTransportJacobian * boundToFreeJacobian; } -BoundMatrix boundToCurvilinearTransportJacobian( +BoundMatrix detail::boundToCurvilinearTransportJacobian( const Vector3& direction, const BoundToFreeMatrix& boundToFreeJacobian, const FreeMatrix& freeTransportJacobian, const FreeVector& freeToPathDerivatives) { - // Calculate the derivative of path length at the curvilinear surface - // w.r.t. free parameters - FreeToPathMatrix freeToPath = FreeToPathMatrix::Zero(); - freeToPath.segment<3>(eFreePos0) = -1.0 * direction; // Calculate the jacobian from global to local at the curvilinear surface FreeToBoundMatrix freeToBoundJacobian = freeToCurvilinearJacobian(direction); - // Calculate the full jacobian from the local parameters at the start surface + + // Update the jacobian to include the derivative of the path length at the + // curvilinear surface w.r.t. the free parameters + freeToBoundJacobian.topLeftCorner<6, 3>() += + (freeToBoundJacobian * freeToPathDerivatives) * + (-1.0 * direction).transpose(); + + // Calculate the full jocobian from the local parameters at the start surface // to curvilinear parameters - return freeToBoundJacobian * - (FreeMatrix::Identity() + freeToPathDerivatives * freeToPath) * - freeTransportJacobian * boundToFreeJacobian; + // @note jac(locA->locB) = jac(gloB->locB)*(1+ + // pathCorrectionFactor(gloB))*jacTransport(gloA->gloB) *jac(locA->gloA) + return blockedMult(freeToBoundJacobian, + blockedMult(freeTransportJacobian, boundToFreeJacobian)); } -BoundToFreeMatrix boundToFreeTransportJacobian( +BoundToFreeMatrix detail::boundToFreeTransportJacobian( const BoundToFreeMatrix& boundToFreeJacobian, const FreeMatrix& freeTransportJacobian) { // Calculate the full jacobian, in this case simple a product of @@ -171,7 +137,7 @@ BoundToFreeMatrix boundToFreeTransportJacobian( return (freeTransportJacobian * boundToFreeJacobian); } -FreeToBoundMatrix freeToBoundTransportJacobian( +FreeToBoundMatrix detail::freeToBoundTransportJacobian( const GeometryContext& geoContext, const FreeVector& freeParameters, const ActsMatrix<7, 8>& directionToAnglesJacobian, const ActsMatrix<8, 7>& anglesToDirectionJacobian, @@ -190,7 +156,7 @@ FreeToBoundMatrix freeToBoundTransportJacobian( directionToAnglesJacobian; } -FreeToBoundMatrix freeToCurvilinearTransportJacobian( +FreeToBoundMatrix detail::freeToCurvilinearTransportJacobian( const Vector3& direction, const ActsMatrix<7, 8>& directionToAnglesJacobian, const ActsMatrix<8, 7>& anglesToDirectionJacobian, const FreeMatrix& freeTransportJacobian, @@ -207,7 +173,7 @@ FreeToBoundMatrix freeToCurvilinearTransportJacobian( directionToAnglesJacobian; } -FreeMatrix freeToFreeTransportJacobian( +FreeMatrix detail::freeToFreeTransportJacobian( const ActsMatrix<7, 8>& directionToAnglesJacobian, const ActsMatrix<8, 7>& anglesToDirectionJacobian, const FreeMatrix& freeTransportJacobian) { @@ -215,5 +181,57 @@ FreeMatrix freeToFreeTransportJacobian( directionToAnglesJacobian; } -} // namespace detail +Result detail::reinitializeJacobians( + const GeometryContext& geoContext, FreeMatrix& freeTransportJacobian, + FreeVector& freeToPathDerivatives, BoundToFreeMatrix& boundToFreeJacobian, + const FreeVector& freeParameters, const Surface& surface) { + using VectorHelpers::phi; + using VectorHelpers::theta; + + // Reset the jacobians + freeTransportJacobian = FreeMatrix::Identity(); + freeToPathDerivatives = FreeVector::Zero(); + + // Get the local position + const Vector3 position = freeParameters.segment<3>(eFreePos0); + const Vector3 direction = freeParameters.segment<3>(eFreeDir0); + auto lpResult = surface.globalToLocal(geoContext, position, direction); + if (!lpResult.ok()) { + return lpResult.error(); + } + // Transform from free to bound parameters + Result boundParameters = detail::transformFreeToBoundParameters( + freeParameters, surface, geoContext); + if (!boundParameters.ok()) { + return boundParameters.error(); + } + // Reset the jacobian from local to global + boundToFreeJacobian = + surface.boundToFreeJacobian(geoContext, *boundParameters); + return Result::success(); +} + +void detail::reinitializeJacobians(FreeMatrix& freeTransportJacobian, + FreeVector& freeToPathDerivatives, + BoundToFreeMatrix& boundToFreeJacobian, + const Vector3& direction) { + // Reset the jacobians + freeTransportJacobian = FreeMatrix::Identity(); + freeToPathDerivatives = FreeVector::Zero(); + boundToFreeJacobian = BoundToFreeMatrix::Zero(); + + auto [cosPhi, sinPhi, cosTheta, sinTheta] = + VectorHelpers::evaluateTrigonomics(direction); + + boundToFreeJacobian(eFreePos0, eBoundLoc0) = -sinPhi; + boundToFreeJacobian(eFreePos0, eBoundLoc1) = -cosPhi * cosTheta; + boundToFreeJacobian(eFreePos1, eBoundLoc0) = cosPhi; + boundToFreeJacobian(eFreePos1, eBoundLoc1) = -sinPhi * cosTheta; + boundToFreeJacobian(eFreePos2, eBoundLoc1) = sinTheta; + boundToFreeJacobian(eFreeTime, eBoundTime) = 1; + boundToFreeJacobian.block<3, 2>(eFreeDir0, eBoundPhi) = + sphericalToFreeDirectionJacobian(direction); + boundToFreeJacobian(eFreeQOverP, eBoundQOverP) = 1; +} + } // namespace Acts 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 5b0b165004e..c011498af6f 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; // the projection of direction onto ref frame normal 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/Plugins/EDM4hep/include/Acts/Plugins/EDM4hep/EDM4hepUtil.hpp b/Plugins/EDM4hep/include/Acts/Plugins/EDM4hep/EDM4hepUtil.hpp index d3a8db38d47..0ee893415b7 100644 --- a/Plugins/EDM4hep/include/Acts/Plugins/EDM4hep/EDM4hepUtil.hpp +++ b/Plugins/EDM4hep/include/Acts/Plugins/EDM4hep/EDM4hepUtil.hpp @@ -18,7 +18,6 @@ #include "Acts/EventData/detail/TransformationBoundToFree.hpp" #include "Acts/EventData/detail/TransformationFreeToBound.hpp" #include "Acts/Geometry/GeometryContext.hpp" -#include "Acts/Propagator/CovarianceTransport.hpp" #include "Acts/Surfaces/PerigeeSurface.hpp" #include "Acts/Surfaces/Surface.hpp" #include "Acts/Utilities/Logger.hpp" diff --git a/Plugins/EDM4hep/src/EDM4hepUtil.cpp b/Plugins/EDM4hep/src/EDM4hepUtil.cpp index 13c7f01d1ba..3e0361020ab 100644 --- a/Plugins/EDM4hep/src/EDM4hepUtil.cpp +++ b/Plugins/EDM4hep/src/EDM4hepUtil.cpp @@ -8,6 +8,7 @@ #include "Acts/Plugins/EDM4hep/EDM4hepUtil.hpp" +#include "Acts/Definitions/Algebra.hpp" #include "Acts/Definitions/Units.hpp" #include "Acts/EventData/MultiTrajectory.hpp" #include "Acts/EventData/MultiTrajectoryHelpers.hpp" @@ -144,25 +145,9 @@ Parameters convertTrackParametersToEdm4hep(const Acts::GeometryContext& gctx, // Only run covariance conversion if we have a covariance input if (params.covariance()) { - auto boundToFree = - refSurface->boundToFreeJacobian(gctx, params.parameters()); - Acts::FreeMatrix freeCov = - boundToFree * params.covariance().value() * boundToFree.transpose(); - - // ensure we have free pars - if (!freePars.has_value()) { - freePars = makeFreePars(); - } - - Acts::CovarianceCache covCache{freePars.value(), freeCov}; - auto [varNewCov, varNewJac] = Acts::transportCovarianceToBound( - gctx, *refSurface, freePars.value(), covCache); - auto targetCov = std::get(varNewCov); - Acts::ActsSquareMatrix<6> J = jacobianToEdm4hep( targetPars[eBoundTheta], targetPars[eBoundQOverP], Bz); - Acts::ActsSquareMatrix<6> cIn = targetCov.template topLeftCorner<6, 6>(); - result.covariance = J * cIn * J.transpose(); + result.covariance = J * params.covariance().value() * J.transpose(); } result.values[0] = targetPars[Acts::eBoundLoc0]; diff --git a/Tests/Benchmarks/CMakeLists.txt b/Tests/Benchmarks/CMakeLists.txt index a08cb4d0a9b..e452add9ed2 100644 --- a/Tests/Benchmarks/CMakeLists.txt +++ b/Tests/Benchmarks/CMakeLists.txt @@ -23,7 +23,6 @@ endmacro() add_benchmark(AtlasStepper AtlasStepperBenchmark.cpp) add_benchmark(BoundaryCheck BoundaryCheckBenchmark.cpp) add_benchmark(BinUtility BinUtilityBenchmark.cpp) -add_benchmark(CovarianceTransport CovarianceTransportBenchmark.cpp) add_benchmark(EigenStepper EigenStepperBenchmark.cpp) add_benchmark(SolenoidField SolenoidFieldBenchmark.cpp) add_benchmark(SurfaceIntersection SurfaceIntersectionBenchmark.cpp) diff --git a/Tests/Benchmarks/CovarianceTransportBenchmark.cpp b/Tests/Benchmarks/CovarianceTransportBenchmark.cpp deleted file mode 100644 index ce6152468bb..00000000000 --- a/Tests/Benchmarks/CovarianceTransportBenchmark.cpp +++ /dev/null @@ -1,92 +0,0 @@ -// This file is part of the Acts project. -// -// Copyright (C) 2021 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/. - -#include "Acts/Definitions/Algebra.hpp" -#include "Acts/Definitions/TrackParametrization.hpp" -#include "Acts/Geometry/GeometryContext.hpp" -#include "Acts/Propagator/CovarianceTransport.hpp" -#include "Acts/Surfaces/PlaneSurface.hpp" -#include "Acts/Tests/CommonHelpers/BenchmarkTools.hpp" -#include "Acts/Utilities/Logger.hpp" - -#include - -int main(int argc, char* argv[]) { - using namespace Acts; - - std::size_t iterations = 3; - std::size_t runs = 1000; - if (argc >= 2) { - iterations = std::stoi(argv[1]); - } - if (argc >= 3) { - runs = std::stoi(argv[2]); - } - - // Create a test context - GeometryContext tgContext = GeometryContext(); - - Vector3 position{1., 2., 3.}; - ActsScalar time = 4.; - ActsScalar qop = 0.125; - Vector3 direction = Vector3(5., 6., 7.).normalized(); - auto planeSurface = Surface::makeShared(position, direction); - auto otherSurface = Surface::makeShared( - position, Vector3(6., 7., 8.).normalized()); - - // Free & bound parameters - FreeVector freeParameters; - freeParameters << position[0], position[1], position[2], time, direction[0], - direction[1], direction[2], qop; - - BoundVector boundParameters; - boundParameters << 0., 0., VectorHelpers::phi(direction), - VectorHelpers::theta(direction), qop, time; - - std::minstd_rand rng; - std::uniform_real_distribution uniform(0.5, 0.95); - - unsigned int sillyCounter = 0; - - ACTS_LOCAL_LOGGER( - getDefaultLogger("CovarianceTransport", Acts::Logging::Level(0))); - - const auto cov_transport_bound_bound = Acts::Test::microBenchmark( - [&] { - BoundSquareMatrix boundCovariance = - uniform(rng) * BoundSquareMatrix::Identity(); - boundCovariance(eBoundLoc0, eBoundPhi) = 0.076; - boundCovariance(eBoundPhi, eBoundLoc0) = 0.076; - boundCovariance(eBoundLoc0, eBoundQOverP) = -0.022; - boundCovariance(eBoundQOverP, eBoundLoc0) = -0.022; - boundCovariance(eBoundLoc1, eBoundTheta) = -0.007; - boundCovariance(eBoundTheta, eBoundLoc1) = -0.007; - - CovarianceCache covCache(tgContext, *planeSurface, position, - boundParameters, boundCovariance); - - // Transport bound to another bound surface - auto covJacAtBound = transportCovarianceToBound( - tgContext, *otherSurface, freeParameters, covCache); - - // Mimic access to the result - const auto& variantCovariance = std::get<0>(covJacAtBound); - const auto& variantJacobian = std::get<1>(covJacAtBound); - - const auto& covariance = std::get(variantCovariance); - const auto& jacobian = std::get(variantJacobian); - - if (covariance(eBoundLoc0, eBoundLoc0) > 0 && - jacobian(eBoundLoc0, eBoundLoc1) > 0.) { - ++sillyCounter; - } - }, - iterations, runs); - - ACTS_INFO("Execution stats: " << cov_transport_bound_bound); -} diff --git a/Tests/UnitTests/Core/Propagator/CMakeLists.txt b/Tests/UnitTests/Core/Propagator/CMakeLists.txt index e3ba11aff87..4a9925f1b8a 100644 --- a/Tests/UnitTests/Core/Propagator/CMakeLists.txt +++ b/Tests/UnitTests/Core/Propagator/CMakeLists.txt @@ -2,7 +2,6 @@ add_unittest(AtlasStepper AtlasStepperTests.cpp) add_unittest(Auctioneer AuctioneerTests.cpp) add_unittest(ConstrainedStep ConstrainedStepTests.cpp) add_unittest(CovarianceEngine CovarianceEngineTests.cpp) -add_unittest(CovarianceTransport CovarianceTransportTests.cpp) add_unittest(DirectNavigator DirectNavigatorTests.cpp) add_unittest(Extrapolator ExtrapolatorTests.cpp) add_unittest(Jacobian JacobianTests.cpp) diff --git a/Tests/UnitTests/Core/Propagator/CovarianceTransportTests.cpp b/Tests/UnitTests/Core/Propagator/CovarianceTransportTests.cpp deleted file mode 100644 index 018934849d0..00000000000 --- a/Tests/UnitTests/Core/Propagator/CovarianceTransportTests.cpp +++ /dev/null @@ -1,237 +0,0 @@ -// This file is part of the Acts project. -// -// Copyright (C) 2021 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/. - -#include - -#include "Acts/Definitions/Algebra.hpp" -#include "Acts/Definitions/TrackParametrization.hpp" -#include "Acts/Geometry/GeometryContext.hpp" -#include "Acts/Propagator/CovarianceTransport.hpp" -#include "Acts/Surfaces/PlaneSurface.hpp" -#include "Acts/Surfaces/Surface.hpp" -#include "Acts/Utilities/Helpers.hpp" - -#include -#include -#include -#include - -namespace Acts { -namespace Test { - -BOOST_AUTO_TEST_CASE(covariance_transport_invalid) { - CovarianceCache covCache; - BOOST_CHECK(!covCache.applyTransport); -} - -BOOST_AUTO_TEST_CASE(covariance_transport_bound_start) { - // Create a test context - GeometryContext tgContext = GeometryContext(); - - Vector3 position{1., 2., 3.}; - ActsScalar time = 4.; - ActsScalar qop = 0.125; - Vector3 direction = Vector3(5., 6., 7.).normalized(); - auto planeSurface = Surface::makeShared(position, direction); - auto otherSurface = Surface::makeShared( - position, Vector3(6., 7., 8.).normalized()); - - // Free & bound parameters - FreeVector freeParameters; - freeParameters << position[0], position[1], position[2], time, direction[0], - direction[1], direction[2], qop; - - BoundVector boundParameters; - boundParameters << 0., 0., VectorHelpers::phi(direction), - VectorHelpers::theta(direction), qop, time; - - BoundSquareMatrix boundCovariance = 8. * BoundSquareMatrix::Identity(); - - CovarianceCache covCache(tgContext, *planeSurface, position, boundParameters, - boundCovariance); - // Test that the transport should be applied now - BOOST_CHECK(covCache.applyTransport); - // Test that the set covariance is 5x5 and what it has been set - BOOST_CHECK_EQUAL(std::get(covCache.covariance), - boundCovariance); - BOOST_CHECK_THROW(std::get(covCache.covariance), - std::bad_variant_access); - // Test that the bound to free jacobian has a value - BOOST_CHECK(covCache.boundToFreeJacobian.has_value()); - BOOST_CHECK(!covCache.boundToFreeJacobian.value().isApprox( - BoundToFreeMatrix::Zero())); - // Test that the free transport jacobian, derivative is properly set up - BOOST_CHECK_EQUAL(covCache.freeTransportJacobian, FreeMatrix::Identity()); - BOOST_CHECK_EQUAL(covCache.freeToPathDerivatives, FreeVector::Zero()); - // Check that the surface is set - BOOST_CHECK_EQUAL(covCache.atSurface, planeSurface.get()); - BOOST_CHECK_EQUAL(covCache.atPosition, position); - - // Transport bound to the same bound surface - auto covJacAtBound = transportCovarianceToBound(tgContext, *planeSurface, - freeParameters, covCache); - BOOST_CHECK(boundCovariance.isApprox( - std::get(covCache.covariance))); - BOOST_CHECK_THROW(std::get(covCache.covariance), - std::bad_variant_access); - - auto variantCovariance = std::get<0>(covJacAtBound); - auto variantJacobian = std::get<1>(covJacAtBound); - BOOST_CHECK_THROW(std::get(variantCovariance), - std::bad_variant_access); - - // Transport bound to curvilinear on the same place - auto covJacAtCurv = transportCovarianceToCurvilinear(direction, covCache); - BOOST_CHECK(boundCovariance.isApprox( - std::get(covCache.covariance))); - BOOST_CHECK_THROW(std::get(covCache.covariance), - std::bad_variant_access); - - variantCovariance = std::get<0>(covJacAtCurv); - variantJacobian = std::get<1>(covJacAtCurv); - BOOST_CHECK_THROW(std::get(variantCovariance), - std::bad_variant_access); - - // Transport bound to Free on the same place - auto covJacAtFree = transportCovarianceToFree(covCache); - variantCovariance = std::get<0>(covJacAtFree); - variantJacobian = std::get<1>(covJacAtFree); - BOOST_CHECK_THROW(std::get(variantCovariance), - std::bad_variant_access); - - // Transport bound to another bound surface - covJacAtBound = transportCovarianceToBound(tgContext, *otherSurface, - freeParameters, covCache); - - variantCovariance = std::get<0>(covJacAtBound); - variantJacobian = std::get<1>(covJacAtBound); - BOOST_CHECK_THROW(std::get(variantCovariance), - std::bad_variant_access); -} - -BOOST_AUTO_TEST_CASE(covariance_transport_curvilinear_start) { - // Create a test context - GeometryContext tgContext = GeometryContext(); - - Vector3 position{1., 2., 3.}; - Vector3 direction = Vector3(5., 6., 7.).normalized(); - ActsScalar time = 4.; - ActsScalar qop = 0.125; - - // Free & bound parameters - FreeVector freeParameters; - freeParameters << position[0], position[1], position[2], time, direction[0], - direction[1], direction[2], qop; - - auto planeSurface = Surface::makeShared(position, direction); - - BoundSquareMatrix boundCovariance = 8. * BoundSquareMatrix::Identity(); - - CovarianceCache covCache(position, direction, boundCovariance); - // Test that the transport should be applied now - BOOST_CHECK(covCache.applyTransport); - // Test that the set covariance is 5x5 and what it has been set - BOOST_CHECK_EQUAL(std::get(covCache.covariance), - boundCovariance); - BOOST_CHECK_THROW(std::get(covCache.covariance), - std::bad_variant_access); - - // Test that the bound to free jacobian has a value - BOOST_CHECK(covCache.boundToFreeJacobian.has_value()); - BOOST_CHECK(!covCache.boundToFreeJacobian.value().isApprox( - BoundToFreeMatrix::Zero())); - // Test that the free transport jacobian, derivative is properly set up - BOOST_CHECK_EQUAL(covCache.freeTransportJacobian, FreeMatrix::Identity()); - BOOST_CHECK_EQUAL(covCache.freeToPathDerivatives, FreeVector::Zero()); - // Check that the surface is set - BOOST_CHECK_EQUAL(covCache.atSurface, nullptr); - BOOST_CHECK_EQUAL(covCache.atPosition, position); - - // Transport bound to the same bound surface - auto covJacAtBound = transportCovarianceToBound(tgContext, *planeSurface, - freeParameters, covCache); - auto variantCovariance = std::get<0>(covJacAtBound); - auto variantJacobian = std::get<1>(covJacAtBound); - BOOST_CHECK_THROW(std::get(variantCovariance), - std::bad_variant_access); - - // Transport bound to curvilinear on the same place - auto covJacAtCurv = transportCovarianceToCurvilinear(direction, covCache); - variantCovariance = std::get<0>(covJacAtCurv); - variantJacobian = std::get<1>(covJacAtCurv); - BOOST_CHECK_THROW(std::get(variantCovariance), - std::bad_variant_access); - - // Transport bound to Free on the same place - auto covJacAtFree = transportCovarianceToFree(covCache); - variantCovariance = std::get<0>(covJacAtFree); - variantJacobian = std::get<1>(covJacAtFree); - BOOST_CHECK_THROW(std::get(variantCovariance), - std::bad_variant_access); -} - -BOOST_AUTO_TEST_CASE(covariance_transport_free_start) { - // Create a test context - GeometryContext tgContext = GeometryContext(); - - // Some start parameters - Vector3 position{1., 2., 3.}; - ActsScalar time = 4.; - Vector3 direction = Vector3(5., 6., 7.).normalized(); - - auto planeSurface = Surface::makeShared(position, direction); - - ActsScalar qop = 0.125; - FreeVector freeParameters; - freeParameters << position[0], position[1], position[2], time, direction[0], - direction[1], direction[2], qop; - - FreeSquareMatrix freeCovariance = 8. * FreeSquareMatrix::Identity(); - - CovarianceCache covCache(freeParameters, freeCovariance); - BOOST_CHECK(covCache.applyTransport); - // Test that the set covariance is 5x5 and what it has been set - BOOST_CHECK_THROW(std::get(covCache.covariance), - std::bad_variant_access); - BOOST_CHECK_EQUAL(std::get(covCache.covariance), - freeCovariance); - // Test that the bound to free jacobian has NO value - BOOST_CHECK(!covCache.boundToFreeJacobian.has_value()); - // Test that the free transport jacobian, derivative is properly set up - BOOST_CHECK_EQUAL(covCache.freeTransportJacobian, FreeMatrix::Identity()); - BOOST_CHECK_EQUAL(covCache.freeToPathDerivatives, FreeVector::Zero()); - // Check that the surface is NOT set - BOOST_CHECK_EQUAL(covCache.atSurface, nullptr); - BOOST_CHECK_EQUAL(covCache.atPosition, position); - - // Transport bound to the same bound surface - auto covJacAtBound = transportCovarianceToBound(tgContext, *planeSurface, - freeParameters, covCache); - - auto variantCovariance = std::get<0>(covJacAtBound); - auto variantJacobian = std::get<1>(covJacAtBound); - BOOST_CHECK_THROW(std::get(variantCovariance), - std::bad_variant_access); - - // Transport bound to curvilinear on the same place - auto covJacAtCurv = transportCovarianceToCurvilinear(direction, covCache); - variantCovariance = std::get<0>(covJacAtCurv); - variantJacobian = std::get<1>(covJacAtCurv); - BOOST_CHECK_THROW(std::get(variantCovariance), - std::bad_variant_access); - - // Transport bound to Free on the same place - auto covJacAtFree = transportCovarianceToFree(covCache); - variantCovariance = std::get<0>(covJacAtFree); - variantJacobian = std::get<1>(covJacAtFree); - BOOST_CHECK_THROW(std::get(variantCovariance), - std::bad_variant_access); -} - -} // namespace Test -} // namespace Acts From 88fd155f9d7cf13826dd922ad798344bb04f1d1d Mon Sep 17 00:00:00 2001 From: andiwand Date: Mon, 11 Dec 2023 09:55:52 +0100 Subject: [PATCH 2/9] draft --- .../Acts/Propagator/detail/JacobianEngine.hpp | 31 ++----------------- Core/src/Propagator/detail/JacobianEngine.cpp | 30 +++++------------- .../Core/Propagator/JacobianEngineTests.cpp | 30 ++---------------- 3 files changed, 12 insertions(+), 79 deletions(-) diff --git a/Core/include/Acts/Propagator/detail/JacobianEngine.hpp b/Core/include/Acts/Propagator/detail/JacobianEngine.hpp index 1f5bc8da6a6..83b7e2a10df 100644 --- a/Core/include/Acts/Propagator/detail/JacobianEngine.hpp +++ b/Core/include/Acts/Propagator/detail/JacobianEngine.hpp @@ -114,20 +114,13 @@ BoundToFreeMatrix boundToFreeTransportJacobian( /// /// @param [in] geoContext The geometry Context /// @param [in] freeParameters Free, nominal parametrisation -/// @param [in] directionToAnglesJacobian The relation jacobian from dir to -/// angle -/// @param [in] anglesToDirectionJacobian The relation jacobian from angle to -/// dir /// @param [in] freeTransportJacobian Transport jacobian free to free -/// @param [in] freeToPathDerivatives Path length derivatives for free -/// parameters +/// @param [in] freeToPathDerivatives Path length derivatives for free parameters /// @param [in] surface The target surface /// /// @return the 6x8 transport jacobian from bound to free FreeToBoundMatrix freeToBoundTransportJacobian( const GeometryContext& geoContext, const FreeVector& freeParameters, - const ActsMatrix<7, 8>& directionToAnglesJacobian, - const ActsMatrix<8, 7>& anglesToDirectionJacobian, const FreeMatrix& freeTransportJacobian, const FreeVector& freeToPathDerivatives, const Surface& surface); @@ -138,34 +131,14 @@ FreeToBoundMatrix freeToBoundTransportJacobian( /// approximated approach to treat the (assumed) small change. /// /// @param [in] direction Normalised direction vector -/// @param [in] directionToAnglesJacobian The relation jacobian from dir to -/// angle -/// @param [in] anglesToDirectionJacobian The relation jacobian from angle to -/// dir /// @param [in] freeTransportJacobian Transport jacobian free to free /// @param [in] freeToPathDerivatives Path length derivatives for free /// /// @return a 6x8 transport jacobian from curvilinear to free FreeToBoundMatrix freeToCurvilinearTransportJacobian( - const Vector3& direction, const ActsMatrix<7, 8>& directionToAnglesJacobian, - const ActsMatrix<8, 7>& anglesToDirectionJacobian, - const FreeMatrix& freeTransportJacobian, + const Vector3& direction, const FreeMatrix& freeTransportJacobian, const FreeVector& freeToPathDerivatives); -/// @brief This function calculates the free transfport jacobian from a free -/// parameterisation. -/// @param [in] directionToAnglesJacobian The relation jacobian from dir to -/// angle -/// @param [in] anglesToDirectionJacobian The relation jacobian from angle to -/// dir -/// @param [in] freeTransportJacobian Transport jacobian free to free -/// -/// @return a 8x8 transport jacobian from free to free -FreeMatrix freeToFreeTransportJacobian( - const ActsMatrix<7, 8>& directionToAnglesJacobian, - const ActsMatrix<8, 7>& anglesToDirectionJacobian, - const FreeMatrix& freeTransportJacobian); - /// @brief This function reinitialises the state members required for the /// covariance transport /// diff --git a/Core/src/Propagator/detail/JacobianEngine.cpp b/Core/src/Propagator/detail/JacobianEngine.cpp index a84b7488eec..0a261bce5b0 100644 --- a/Core/src/Propagator/detail/JacobianEngine.cpp +++ b/Core/src/Propagator/detail/JacobianEngine.cpp @@ -11,6 +11,7 @@ #include "Acts/Definitions/Algebra.hpp" #include "Acts/Definitions/Tolerance.hpp" #include "Acts/Definitions/TrackParametrization.hpp" +#include "Acts/EventData/detail/TransformationFreeToBound.hpp" #include "Acts/Geometry/GeometryContext.hpp" #include "Acts/Surfaces/Surface.hpp" #include "Acts/Utilities/AlgebraHelpers.hpp" @@ -139,46 +140,29 @@ BoundToFreeMatrix detail::boundToFreeTransportJacobian( FreeToBoundMatrix detail::freeToBoundTransportJacobian( const GeometryContext& geoContext, const FreeVector& freeParameters, - const ActsMatrix<7, 8>& directionToAnglesJacobian, - const ActsMatrix<8, 7>& anglesToDirectionJacobian, const FreeMatrix& freeTransportJacobian, const FreeVector& freeToPathDerivatives, const Surface& surface) { // Calculate the jacobian from free to bound at the final surface FreeToBoundMatrix freeToBoundJacobian = surface.freeToBoundJacobian(geoContext, freeParameters); - // Calculate the form factors for the derivatives - auto transport = freeTransportJacobian * anglesToDirectionJacobian; FreeToPathMatrix sVec = surface.freeToPathDerivative(geoContext, freeParameters); // Return the jacobian to local return freeToBoundJacobian * - (transport + freeToPathDerivatives * sVec * transport) * - directionToAnglesJacobian; + (freeTransportJacobian + + freeToPathDerivatives * sVec * freeTransportJacobian); } FreeToBoundMatrix detail::freeToCurvilinearTransportJacobian( - const Vector3& direction, const ActsMatrix<7, 8>& directionToAnglesJacobian, - const ActsMatrix<8, 7>& anglesToDirectionJacobian, - const FreeMatrix& freeTransportJacobian, + const Vector3& direction, const FreeMatrix& freeTransportJacobian, const FreeVector& freeToPathDerivatives) { - const ActsMatrix<8, 7> transport = - freeTransportJacobian * anglesToDirectionJacobian; - auto sfactors = - direction.transpose() * transport.template topLeftCorner<3, 7>(); + auto sfactors = direction.transpose() * + freeTransportJacobian.template topLeftCorner<3, 8>(); // Since the jacobian to local needs to calculated for the bound parameters // here, it is convenient to do the same here return freeToCurvilinearJacobian(direction) * - (transport - freeToPathDerivatives * sfactors) * - directionToAnglesJacobian; -} - -FreeMatrix detail::freeToFreeTransportJacobian( - const ActsMatrix<7, 8>& directionToAnglesJacobian, - const ActsMatrix<8, 7>& anglesToDirectionJacobian, - const FreeMatrix& freeTransportJacobian) { - return freeTransportJacobian * anglesToDirectionJacobian * - directionToAnglesJacobian; + (freeTransportJacobian - freeToPathDerivatives * sfactors); } Result detail::reinitializeJacobians( diff --git a/Tests/UnitTests/Core/Propagator/JacobianEngineTests.cpp b/Tests/UnitTests/Core/Propagator/JacobianEngineTests.cpp index 07562814971..5183df23c1a 100644 --- a/Tests/UnitTests/Core/Propagator/JacobianEngineTests.cpp +++ b/Tests/UnitTests/Core/Propagator/JacobianEngineTests.cpp @@ -90,16 +90,6 @@ BOOST_AUTO_TEST_CASE(jacobian_engine_helper) { CHECK_CLOSE_REL(c2fJacobian(eFreeDir2, eBoundTheta), -sinTheta, 1e-5); // Q/P parameter: stays as is CHECK_CLOSE_REL(c2fJacobian(eFreeQOverP, eBoundQOverP), 1, 1e-5); - - // (3) Test angle - directio relational jacobians - direction = Vector3(2., 3., 4.).normalized(); - - ActsMatrix<7, 8> d2aJacobian = detail::directionToAnglesJacobian(direction); - ActsMatrix<8, 7> a2dJacobian = detail::anglesToDirectionJacobian(direction); - - auto roundTrip = a2dJacobian * d2aJacobian; - BOOST_CHECK_EQUAL(roundTrip.cols(), 8); - BOOST_CHECK_EQUAL(roundTrip.rows(), 8); } /// These tests do not test for a correct covariance transport but only for the @@ -168,15 +158,8 @@ BOOST_AUTO_TEST_CASE(jacobian_engine_to_bound) { b2bTransportJacobian * boundCovariance * b2bTransportJacobian.transpose(); BOOST_CHECK(!boundCovariance.isApprox(newBoundCovariance)); - // (2) free to bound transport jacobian - const ActsMatrix<7, 8>& directionToAnglesJacobian = - detail::directionToAnglesJacobian(direction); - const ActsMatrix<8, 7>& anglesToDirectionJacobian = - detail::anglesToDirectionJacobian(direction); - FreeToBoundMatrix f2bTransportJacobian = detail::freeToBoundTransportJacobian( - tgContext, freeParameters, directionToAnglesJacobian, - anglesToDirectionJacobian, noTransportJacobian, freeToPathDerivatives, + tgContext, freeParameters, noTransportJacobian, freeToPathDerivatives, *pSurface); newBoundCovariance = @@ -238,16 +221,9 @@ BOOST_AUTO_TEST_CASE(jacobian_engine_to_curvilinear) { b2cTransportJacobian * boundCovariance * b2cTransportJacobian.transpose(); BOOST_CHECK(!boundCovariance.isApprox(newBoundCovariance)); - // (2) free to bound transport jacobian - const ActsMatrix<7, 8>& directionToAnglesJacobian = - detail::directionToAnglesJacobian(direction); - const ActsMatrix<8, 7>& anglesToDirectionJacobian = - detail::anglesToDirectionJacobian(direction); - FreeToBoundMatrix f2cTransportJacobian = - detail::freeToCurvilinearTransportJacobian( - direction, directionToAnglesJacobian, anglesToDirectionJacobian, - noTransportJacobian, freeToPathDerivatives); + detail::freeToCurvilinearTransportJacobian(direction, noTransportJacobian, + freeToPathDerivatives); newBoundCovariance = f2cTransportJacobian * freeCovariance * f2cTransportJacobian.transpose(); From 0527e57cb48435d58e3ba2de73eb0cc62e22557c Mon Sep 17 00:00:00 2001 From: andiwand Date: Mon, 11 Dec 2023 14:40:52 +0100 Subject: [PATCH 3/9] refactor params --- Core/include/Acts/Propagator/EigenStepper.ipp | 12 +-- .../Propagator/detail/CovarianceEngine.hpp | 68 +++++++------- .../Acts/Propagator/detail/JacobianEngine.hpp | 43 ++++----- Core/src/Propagator/StraightLineStepper.cpp | 9 +- .../Propagator/detail/CovarianceEngine.cpp | 88 ++++++++++--------- Core/src/Propagator/detail/JacobianEngine.cpp | 37 ++++---- 6 files changed, 131 insertions(+), 126 deletions(-) diff --git a/Core/include/Acts/Propagator/EigenStepper.ipp b/Core/include/Acts/Propagator/EigenStepper.ipp index bf0c29dde72..f4dee1c9158 100644 --- a/Core/include/Acts/Propagator/EigenStepper.ipp +++ b/Core/include/Acts/Propagator/EigenStepper.ipp @@ -51,9 +51,9 @@ auto Acts::EigenStepper::boundState( const FreeToBoundCorrection& freeToBoundCorrection) const -> Result { return detail::boundState( - state.geoContext, state.cov, state.jacobian, state.jacTransport, + state.geoContext, surface, state.cov, state.jacobian, state.jacTransport, state.derivative, state.jacToGlobal, state.pars, state.particleHypothesis, - state.covTransport && transportCov, state.pathAccumulated, surface, + state.covTransport && transportCov, state.pathAccumulated, freeToBoundCorrection); } @@ -101,10 +101,10 @@ template void Acts::EigenStepper::transportCovarianceToBound( State& state, const Surface& surface, const FreeToBoundCorrection& freeToBoundCorrection) const { - detail::transportCovarianceToBound( - state.geoContext.get(), state.cov, state.jacobian, state.jacTransport, - state.derivative, state.jacToGlobal, state.pars, surface, - freeToBoundCorrection); + detail::transportCovarianceToBound(state.geoContext.get(), surface, state.cov, + state.jacobian, state.jacTransport, + state.derivative, state.jacToGlobal, + state.pars, freeToBoundCorrection); } template diff --git a/Core/include/Acts/Propagator/detail/CovarianceEngine.hpp b/Core/include/Acts/Propagator/detail/CovarianceEngine.hpp index eb530fd5e5c..c169f682c7d 100644 --- a/Core/include/Acts/Propagator/detail/CovarianceEngine.hpp +++ b/Core/include/Acts/Propagator/detail/CovarianceEngine.hpp @@ -19,10 +19,7 @@ #include "Acts/Surfaces/Surface.hpp" #include "Acts/Utilities/Result.hpp" -#include -#include #include -#include namespace Acts { namespace detail { @@ -36,48 +33,49 @@ namespace detail { /// Create and return the bound state at the current position /// /// @brief It does not check if the transported state is at the surface, this -/// needs to be guaranteed by the propagator +/// needs to be guaranteed by the propagator /// /// @param [in] geoContext The geometry context -/// @param [in, out] covarianceMatrix The covariance matrix of the state -/// @param [in, out] jacobian Full jacobian since the last reset -/// @param [in, out] transportJacobian Global jacobian since the last reset -/// @param [in, out] derivatives Path length derivatives of the free, nominal -/// parameters -/// @param [in, out] jacToGlobal Projection jacobian of the last bound +/// @param [in, out] boundCovariance The covariance matrix of the state +/// @param [in, out] fullTransportJacobian Full jacobian since the last reset +/// @param [in, out] freeTransportJacobian Global jacobian since the last reset +/// @param [in, out] freeToPathDerivatives Path length derivatives of the free, +/// nominal parameters +/// @param [in, out] boundToFreeJacobian Projection jacobian of the last bound /// parametrisation to free parameters -/// @param [in, out] parameters Free, nominal parametrisation +/// @param [in, out] freeParameters Free, nominal parametrisation /// @param [in] particleHypothesis Particle hypothesis /// @param [in] covTransport Decision whether the covariance transport should be /// performed /// @param [in] accumulatedPath Propagated distance /// @param [in] surface Target surface on which the state is represented -/// @param [in] freeToBoundCorrection Correction for non-linearity effect during transform from free to bound +/// @param [in] freeToBoundCorrection Correction for non-linearity effect during +/// transform from free to bound /// /// @return A bound state: /// - the parameters at the surface /// - the stepwise jacobian towards it (from last bound) /// - and the path length (from start - for ordering) Result> boundState( - const GeometryContext& geoContext, BoundSquareMatrix& covarianceMatrix, - BoundMatrix& jacobian, FreeMatrix& transportJacobian, - FreeVector& derivatives, BoundToFreeMatrix& jacToGlobal, - FreeVector& parameters, const ParticleHypothesis& particleHypothesis, - bool covTransport, double accumulatedPath, const Surface& surface, - const FreeToBoundCorrection& freeToBoundCorrection); + const GeometryContext& geoContext, const Surface& surface, + BoundSquareMatrix& boundCovariance, BoundMatrix& fullTransportJacobian, + FreeMatrix& freeTransportJacobian, FreeVector& freeToPathDerivatives, + BoundToFreeMatrix& boundToFreeJacobian, FreeVector& freeParameters, + const ParticleHypothesis& particleHypothesis, bool covTransport, + double accumulatedPath, const FreeToBoundCorrection& freeToBoundCorrection); /// Create and return a curvilinear state at the current position /// /// @brief This creates a curvilinear state. /// -/// @param [in, out] covarianceMatrix The covariance matrix of the state -/// @param [in, out] jacobian Full jacobian since the last reset -/// @param [in, out] transportJacobian Global jacobian since the last reset -/// @param [in, out] derivatives Path length derivatives of the free, nominal -/// parameters -/// @param [in, out] jacToGlobal Projection jacobian of the last bound +/// @param [in, out] boundCovariance The covariance matrix of the state +/// @param [in, out] fullTransportJacobian Full jacobian since the last reset +/// @param [in, out] freeTransportJacobian Global jacobian since the last reset +/// @param [in, out] freeToPathDerivatives Path length derivatives of the free, +/// nominal parameters +/// @param [in, out] boundToFreeJacobian Projection jacobian of the last bound /// parametrisation to free parameters -/// @param [in] parameters Free, nominal parametrisation +/// @param [in] freeParameters Free, nominal parametrisation /// @param [in] particleHypothesis Particle hypothesis /// @param [in] covTransport Decision whether the covariance transport should be /// performed @@ -88,9 +86,9 @@ Result> boundState( /// - the stepweise jacobian towards it (from last bound) /// - and the path length (from start - for ordering) std::tuple curvilinearState( - BoundSquareMatrix& covarianceMatrix, BoundMatrix& jacobian, - FreeMatrix& transportJacobian, FreeVector& derivatives, - BoundToFreeMatrix& jacToGlobal, const FreeVector& parameters, + BoundSquareMatrix& boundCovariance, BoundMatrix& fullTransportJacobian, + FreeMatrix& transportJacobian, FreeVector& freeToPathDerivatives, + BoundToFreeMatrix& boundToFreeJacobian, const FreeVector& freeParameters, const ParticleHypothesis& particleHypothesis, bool covTransport, double accumulatedPath); @@ -98,6 +96,7 @@ std::tuple curvilinearState( /// another bound representation. /// /// @param [in] geoContext The geometry context +/// @param [in] surface is the surface to which the covariance is forwarded to /// @param [in, out] boundCovariance The covariance matrix of the state /// @param [in, out] fullTransportJacobian Full jacobian since the last reset /// @param [in, out] freeTransportJacobian Global jacobian since the last reset @@ -105,17 +104,16 @@ std::tuple curvilinearState( /// @param [in, out] boundToFreeJacobian Projection jacobian of the last bound /// parametrisation to free parameters /// @param [in, out] freeParameters Free, nominal parametrisation -/// @param [in] surface is the surface to which the covariance is -/// forwarded to -/// @param [in] freeToBoundCorrection Correction for non-linearity effect during transform from free to bound +/// @param [in] freeToBoundCorrection Correction for non-linearity effect during +/// transform from free to bound /// /// @note No check is done if the position is actually on the surface /// void transportCovarianceToBound( - const GeometryContext& geoContext, BoundSquareMatrix& boundCovariance, - BoundMatrix& fullTransportJacobian, FreeMatrix& freeTransportJacobian, - FreeVector& freeToPathDerivatives, BoundToFreeMatrix& boundToFreeJacobian, - FreeVector& freeParameters, const Surface& surface, + const GeometryContext& geoContext, const Surface& surface, + BoundSquareMatrix& boundCovariance, BoundMatrix& fullTransportJacobian, + FreeMatrix& freeTransportJacobian, FreeVector& freeToPathDerivatives, + BoundToFreeMatrix& boundToFreeJacobian, FreeVector& freeParameters, const FreeToBoundCorrection& freeToBoundCorrection); /// @brief Method for on-demand covariance transport of a bound/curvilinear diff --git a/Core/include/Acts/Propagator/detail/JacobianEngine.hpp b/Core/include/Acts/Propagator/detail/JacobianEngine.hpp index 83b7e2a10df..859bb60cb37 100644 --- a/Core/include/Acts/Propagator/detail/JacobianEngine.hpp +++ b/Core/include/Acts/Propagator/detail/JacobianEngine.hpp @@ -49,22 +49,23 @@ BoundToFreeMatrix curvilinearToFreeJacobian(const Vector3& direction); /// approximated approach to treat the (assumed) small change. /// /// @param [in] geoContext The geometry Context +/// @param [in] surface Target surface /// @param [in] freeParameters Free, nominal parametrisation /// @param [in] boundToFreeJacobian Jacobian from bound to free at start /// @param [in] freeTransportJacobian Transport jacobian free to free /// @param [in] freeToPathDerivatives Path length derivatives for free -/// parameters -/// @param [in] surface Target surface +/// parameters +/// @param [out] fullTransportJacobian A 6x6 transport jacobian from bound to bound /// /// @note jac(locA->locB) = jac(gloB->locB)*(1+ -/// pathCorrectionFactor(gloB))*jacTransport(gloA->gloB) *jac(locA->gloA) -/// -/// @return a 6x6 transport jacobian from bound to bound -BoundMatrix boundToBoundTransportJacobian( - const GeometryContext& geoContext, const FreeVector& freeParameters, - const BoundToFreeMatrix& boundToFreeJacobian, - const FreeMatrix& freeTransportJacobian, - const FreeVector& freeToPathDerivatives, const Surface& surface); +/// pathCorrectionFactor(gloB))*jacTransport(gloA->gloB) *jac(locA->gloA) +void boundToBoundTransportJacobian(const GeometryContext& geoContext, + const Surface& surface, + const FreeVector& freeParameters, + const BoundToFreeMatrix& boundToFreeJacobian, + const FreeMatrix& freeTransportJacobian, + const FreeVector& freeToPathDerivatives, + BoundMatrix& fullTransportJacobian); /// @brief This function calculates the full jacobian from a given /// bound/curvilinear parameterisation from a surface to new curvilinear @@ -81,16 +82,16 @@ BoundMatrix boundToBoundTransportJacobian( /// @param [in] freeTransportJacobian Transport jacobian free to free /// @param [in] freeToPathDerivatives Path length derivatives for free /// parameters +/// @param [out] fullTransportJacobian A 6x6 transport jacobian from curilinear to bound /// /// @note The parameter @p surface is only required if projected to bound /// parameters. In the case of curvilinear parameters the geometry and the /// position is known and the calculation can be simplified -/// -/// @return a 6x6 transport jacobian from curilinear to bound -BoundMatrix boundToCurvilinearTransportJacobian( +void boundToCurvilinearTransportJacobian( const Vector3& direction, const BoundToFreeMatrix& boundToFreeJacobian, const FreeMatrix& freeTransportJacobian, - const FreeVector& freeToPathDerivatives); + const FreeVector& freeToPathDerivatives, + BoundMatrix& fullTransportJacobian); /// @brief This function calculates the full jacobian from a given /// bound/curvilinear parameterisation from a new free parameterisation. @@ -113,16 +114,16 @@ BoundToFreeMatrix boundToFreeTransportJacobian( /// This is an approximated approach to treat the (assumed) small change. /// /// @param [in] geoContext The geometry Context +/// @param [in] surface The target surface /// @param [in] freeParameters Free, nominal parametrisation /// @param [in] freeTransportJacobian Transport jacobian free to free /// @param [in] freeToPathDerivatives Path length derivatives for free parameters -/// @param [in] surface The target surface /// /// @return the 6x8 transport jacobian from bound to free FreeToBoundMatrix freeToBoundTransportJacobian( - const GeometryContext& geoContext, const FreeVector& freeParameters, - const FreeMatrix& freeTransportJacobian, - const FreeVector& freeToPathDerivatives, const Surface& surface); + const GeometryContext& geoContext, const Surface& surface, + const FreeVector& freeParameters, const FreeMatrix& freeTransportJacobian, + const FreeVector& freeToPathDerivatives); /// @brief This function calculates the full transport jacobian from a free /// parameterisation to a bound one. Since a variation of the start @@ -143,6 +144,7 @@ FreeToBoundMatrix freeToCurvilinearTransportJacobian( /// covariance transport /// /// @param [in] geoContext The geometry context +/// @param [in] surface The reference surface of the local parametrisation /// @param [in, out] freeTransportJacobian The transport jacobian from start /// free to final free parameters /// @param [in, out] freeToPathDerivatives Path length derivatives of the free, @@ -150,13 +152,12 @@ FreeToBoundMatrix freeToCurvilinearTransportJacobian( /// @param [in, out] boundToFreeJacobian Projection jacobian of the last bound /// parametrisation to free parameters /// @param [in] freeParameters Free, nominal parametrisation -/// @param [in] surface The reference surface of the local parametrisation Result reinitializeJacobians(const GeometryContext& geoContext, + const Surface& surface, FreeMatrix& freeTransportJacobian, FreeVector& freeToPathDerivatives, BoundToFreeMatrix& boundToFreeJacobian, - const FreeVector& freeParameters, - const Surface& surface); + const FreeVector& freeParameters); /// @brief This function reinitialises the state members required for the /// covariance transport diff --git a/Core/src/Propagator/StraightLineStepper.cpp b/Core/src/Propagator/StraightLineStepper.cpp index e82a4fe39c5..ec66b8b864f 100644 --- a/Core/src/Propagator/StraightLineStepper.cpp +++ b/Core/src/Propagator/StraightLineStepper.cpp @@ -18,9 +18,9 @@ StraightLineStepper::boundState( State& state, const Surface& surface, bool transportCov, const FreeToBoundCorrection& freeToBoundCorrection) const { return detail::boundState( - state.geoContext, state.cov, state.jacobian, state.jacTransport, + state.geoContext, surface, state.cov, state.jacobian, state.jacTransport, state.derivative, state.jacToGlobal, state.pars, state.particleHypothesis, - state.covTransport && transportCov, state.pathAccumulated, surface, + state.covTransport && transportCov, state.pathAccumulated, freeToBoundCorrection); } @@ -61,9 +61,8 @@ void StraightLineStepper::transportCovarianceToBound( State& state, const Surface& surface, const FreeToBoundCorrection& freeToBoundCorrection) const { detail::transportCovarianceToBound( - state.geoContext, state.cov, state.jacobian, state.jacTransport, - state.derivative, state.jacToGlobal, state.pars, surface, - freeToBoundCorrection); + state.geoContext, surface, state.cov, state.jacobian, state.jacTransport, + state.derivative, state.jacToGlobal, state.pars, freeToBoundCorrection); } void StraightLineStepper::resetState(State& state, diff --git a/Core/src/Propagator/detail/CovarianceEngine.cpp b/Core/src/Propagator/detail/CovarianceEngine.cpp index 575b093cc2e..01198857ab9 100644 --- a/Core/src/Propagator/detail/CovarianceEngine.cpp +++ b/Core/src/Propagator/detail/CovarianceEngine.cpp @@ -31,31 +31,33 @@ using CurvilinearState = std::tuple; Result detail::boundState( - const GeometryContext& geoContext, BoundSquareMatrix& covarianceMatrix, - BoundMatrix& jacobian, FreeMatrix& transportJacobian, - FreeVector& derivatives, BoundToFreeMatrix& jacToGlobal, - FreeVector& parameters, const ParticleHypothesis& particleHypothesis, - bool covTransport, double accumulatedPath, const Surface& surface, + const GeometryContext& geoContext, const Surface& surface, + BoundSquareMatrix& boundCovariance, BoundMatrix& fullTransportJacobian, + FreeMatrix& freeTransportJacobian, FreeVector& freeToPathDerivatives, + BoundToFreeMatrix& boundToFreeJacobian, FreeVector& freeParameters, + const ParticleHypothesis& particleHypothesis, bool covTransport, + double accumulatedPath, const FreeToBoundCorrection& freeToBoundCorrection) { // Covariance transport std::optional cov = std::nullopt; if (covTransport) { // Initialize the jacobian from start local to final local - jacobian = BoundMatrix::Identity(); + fullTransportJacobian = BoundMatrix::Identity(); // Calculate the jacobian and transport the covarianceMatrix to final local. // Then reinitialize the transportJacobian, derivatives and the - // jacToGlobal - transportCovarianceToBound(geoContext, covarianceMatrix, jacobian, - transportJacobian, derivatives, jacToGlobal, - parameters, surface, freeToBoundCorrection); + // boundToFreeJacobian + transportCovarianceToBound(geoContext, surface, boundCovariance, + fullTransportJacobian, freeTransportJacobian, + freeToPathDerivatives, boundToFreeJacobian, + freeParameters, freeToBoundCorrection); } - if (covarianceMatrix != BoundSquareMatrix::Zero()) { - cov = covarianceMatrix; + if (boundCovariance != BoundSquareMatrix::Zero()) { + cov = boundCovariance; } // Create the bound parameters - Result bv = - detail::transformFreeToBoundParameters(parameters, surface, geoContext); + Result bv = detail::transformFreeToBoundParameters( + freeParameters, surface, geoContext); if (!bv.ok()) { return bv.error(); } @@ -63,58 +65,58 @@ Result detail::boundState( return std::make_tuple( BoundTrackParameters(surface.getSharedPtr(), *bv, std::move(cov), particleHypothesis), - jacobian, accumulatedPath); + fullTransportJacobian, accumulatedPath); } CurvilinearState detail::curvilinearState( - BoundSquareMatrix& covarianceMatrix, BoundMatrix& jacobian, - FreeMatrix& transportJacobian, FreeVector& derivatives, - BoundToFreeMatrix& jacToGlobal, const FreeVector& parameters, + BoundSquareMatrix& boundCovariance, BoundMatrix& fullTransportJacobian, + FreeMatrix& freeTransportJacobian, FreeVector& freeToPathDerivatives, + BoundToFreeMatrix& boundToFreeJacobian, const FreeVector& freeParameters, const ParticleHypothesis& particleHypothesis, bool covTransport, double accumulatedPath) { - const Vector3& direction = parameters.segment<3>(eFreeDir0); + const Vector3& direction = freeParameters.segment<3>(eFreeDir0); // Covariance transport std::optional cov = std::nullopt; if (covTransport) { // Initialize the jacobian from start local to final local - jacobian = BoundMatrix::Identity(); + fullTransportJacobian = BoundMatrix::Identity(); // Calculate the jacobian and transport the covarianceMatrix to final local. // Then reinitialize the transportJacobian, derivatives and the - // jacToGlobal - transportCovarianceToCurvilinear(covarianceMatrix, jacobian, - transportJacobian, derivatives, - jacToGlobal, direction); + // boundToFreeJacobian + transportCovarianceToCurvilinear( + boundCovariance, fullTransportJacobian, freeTransportJacobian, + freeToPathDerivatives, boundToFreeJacobian, direction); } - if (covarianceMatrix != BoundSquareMatrix::Zero()) { - cov = covarianceMatrix; + if (boundCovariance != BoundSquareMatrix::Zero()) { + cov = boundCovariance; } // Create the curvilinear parameters Vector4 pos4 = Vector4::Zero(); - pos4[ePos0] = parameters[eFreePos0]; - pos4[ePos1] = parameters[eFreePos1]; - pos4[ePos2] = parameters[eFreePos2]; - pos4[eTime] = parameters[eFreeTime]; + pos4[ePos0] = freeParameters[eFreePos0]; + pos4[ePos1] = freeParameters[eFreePos1]; + pos4[ePos2] = freeParameters[eFreePos2]; + pos4[eTime] = freeParameters[eFreeTime]; CurvilinearTrackParameters curvilinearParams( - pos4, direction, parameters[eFreeQOverP], std::move(cov), + pos4, direction, freeParameters[eFreeQOverP], std::move(cov), particleHypothesis); // Create the curvilinear state - return std::make_tuple(std::move(curvilinearParams), jacobian, + return std::make_tuple(std::move(curvilinearParams), fullTransportJacobian, accumulatedPath); } void detail::transportCovarianceToBound( - const GeometryContext& geoContext, BoundSquareMatrix& boundCovariance, - BoundMatrix& fullTransportJacobian, FreeMatrix& freeTransportJacobian, - FreeVector& freeToPathDerivatives, BoundToFreeMatrix& boundToFreeJacobian, - FreeVector& freeParameters, const Surface& surface, + const GeometryContext& geoContext, const Surface& surface, + BoundSquareMatrix& boundCovariance, BoundMatrix& fullTransportJacobian, + FreeMatrix& freeTransportJacobian, FreeVector& freeToPathDerivatives, + BoundToFreeMatrix& boundToFreeJacobian, FreeVector& freeParameters, const FreeToBoundCorrection& freeToBoundCorrection) { // Calculate the full jacobian from local parameters at the start surface to // current bound parameters - fullTransportJacobian = boundToBoundTransportJacobian( - geoContext, freeParameters, boundToFreeJacobian, freeTransportJacobian, - freeToPathDerivatives, surface); + boundToBoundTransportJacobian(geoContext, surface, freeParameters, + boundToFreeJacobian, freeTransportJacobian, + freeToPathDerivatives, fullTransportJacobian); bool correction = false; if (freeToBoundCorrection) { @@ -154,9 +156,9 @@ void detail::transportCovarianceToBound( // ->The transportJacobian is reinitialized to Identity // ->The derivatives is reinitialized to Zero // ->The boundToFreeJacobian is initialized to that at the current surface - reinitializeJacobians(geoContext, freeTransportJacobian, + reinitializeJacobians(geoContext, surface, freeTransportJacobian, freeToPathDerivatives, boundToFreeJacobian, - freeParameters, surface); + freeParameters); } void detail::transportCovarianceToCurvilinear( @@ -165,9 +167,9 @@ void detail::transportCovarianceToCurvilinear( BoundToFreeMatrix& boundToFreeJacobian, const Vector3& direction) { // Calculate the full jacobian from local parameters at the start surface to // current curvilinear parameters - fullTransportJacobian = boundToCurvilinearTransportJacobian( + boundToCurvilinearTransportJacobian( direction, boundToFreeJacobian, freeTransportJacobian, - freeToPathDerivatives); + freeToPathDerivatives, fullTransportJacobian); // Apply the actual covariance transport to get covariance of the current // curvilinear parameters diff --git a/Core/src/Propagator/detail/JacobianEngine.cpp b/Core/src/Propagator/detail/JacobianEngine.cpp index 0a261bce5b0..8c4512ea4a7 100644 --- a/Core/src/Propagator/detail/JacobianEngine.cpp +++ b/Core/src/Propagator/detail/JacobianEngine.cpp @@ -87,11 +87,13 @@ BoundToFreeMatrix detail::curvilinearToFreeJacobian(const Vector3& direction) { return curvToFreeJacobian; } -BoundMatrix detail::boundToBoundTransportJacobian( - const GeometryContext& geoContext, const FreeVector& freeParameters, +void detail::boundToBoundTransportJacobian( + const GeometryContext& geoContext, const Surface& surface, + const FreeVector& freeParameters, const BoundToFreeMatrix& boundToFreeJacobian, const FreeMatrix& freeTransportJacobian, - const FreeVector& freeToPathDerivatives, const Surface& surface) { + const FreeVector& freeToPathDerivatives, + BoundMatrix& fullTransportJacobian) { // Calculate the derivative of path length at the final surface or the // point-of-closest approach w.r.t. free parameters const FreeToPathMatrix freeToPath = @@ -104,15 +106,17 @@ BoundMatrix detail::boundToBoundTransportJacobian( // surface to local/bound parameters at the final surface // @note jac(locA->locB) = jac(gloB->locB)*(1+ // pathCorrectionFactor(gloB))*jacTransport(gloA->gloB) *jac(locA->gloA) - return freeToBoundJacobian * - (FreeMatrix::Identity() + freeToPathDerivatives * freeToPath) * - freeTransportJacobian * boundToFreeJacobian; + fullTransportJacobian = + freeToBoundJacobian * + (FreeMatrix::Identity() + freeToPathDerivatives * freeToPath) * + freeTransportJacobian * boundToFreeJacobian; } -BoundMatrix detail::boundToCurvilinearTransportJacobian( +void detail::boundToCurvilinearTransportJacobian( const Vector3& direction, const BoundToFreeMatrix& boundToFreeJacobian, const FreeMatrix& freeTransportJacobian, - const FreeVector& freeToPathDerivatives) { + const FreeVector& freeToPathDerivatives, + BoundMatrix& fullTransportJacobian) { // Calculate the jacobian from global to local at the curvilinear surface FreeToBoundMatrix freeToBoundJacobian = freeToCurvilinearJacobian(direction); @@ -126,8 +130,9 @@ BoundMatrix detail::boundToCurvilinearTransportJacobian( // to curvilinear parameters // @note jac(locA->locB) = jac(gloB->locB)*(1+ // pathCorrectionFactor(gloB))*jacTransport(gloA->gloB) *jac(locA->gloA) - return blockedMult(freeToBoundJacobian, - blockedMult(freeTransportJacobian, boundToFreeJacobian)); + fullTransportJacobian = + blockedMult(freeToBoundJacobian, + blockedMult(freeTransportJacobian, boundToFreeJacobian)); } BoundToFreeMatrix detail::boundToFreeTransportJacobian( @@ -139,9 +144,9 @@ BoundToFreeMatrix detail::boundToFreeTransportJacobian( } FreeToBoundMatrix detail::freeToBoundTransportJacobian( - const GeometryContext& geoContext, const FreeVector& freeParameters, - const FreeMatrix& freeTransportJacobian, - const FreeVector& freeToPathDerivatives, const Surface& surface) { + const GeometryContext& geoContext, const Surface& surface, + const FreeVector& freeParameters, const FreeMatrix& freeTransportJacobian, + const FreeVector& freeToPathDerivatives) { // Calculate the jacobian from free to bound at the final surface FreeToBoundMatrix freeToBoundJacobian = surface.freeToBoundJacobian(geoContext, freeParameters); @@ -166,9 +171,9 @@ FreeToBoundMatrix detail::freeToCurvilinearTransportJacobian( } Result detail::reinitializeJacobians( - const GeometryContext& geoContext, FreeMatrix& freeTransportJacobian, - FreeVector& freeToPathDerivatives, BoundToFreeMatrix& boundToFreeJacobian, - const FreeVector& freeParameters, const Surface& surface) { + const GeometryContext& geoContext, const Surface& surface, + FreeMatrix& freeTransportJacobian, FreeVector& freeToPathDerivatives, + BoundToFreeMatrix& boundToFreeJacobian, const FreeVector& freeParameters) { using VectorHelpers::phi; using VectorHelpers::theta; From e3fbd1b2e03e940d2d5b60e1f98dd2dd7b56c73c Mon Sep 17 00:00:00 2001 From: andiwand Date: Mon, 11 Dec 2023 14:47:35 +0100 Subject: [PATCH 4/9] fix --- .../Acts/Propagator/detail/JacobianEngine.hpp | 14 ++++--- .../Propagator/detail/LoopStepperUtils.hpp | 10 ++--- Core/src/Propagator/detail/JacobianEngine.cpp | 11 +++--- .../Core/Propagator/CovarianceEngineTests.cpp | 22 +++++------ .../Core/Propagator/JacobianEngineTests.cpp | 38 ++++++++++--------- 5 files changed, 50 insertions(+), 45 deletions(-) diff --git a/Core/include/Acts/Propagator/detail/JacobianEngine.hpp b/Core/include/Acts/Propagator/detail/JacobianEngine.hpp index 859bb60cb37..d2642bc7dc2 100644 --- a/Core/include/Acts/Propagator/detail/JacobianEngine.hpp +++ b/Core/include/Acts/Propagator/detail/JacobianEngine.hpp @@ -118,12 +118,14 @@ BoundToFreeMatrix boundToFreeTransportJacobian( /// @param [in] freeParameters Free, nominal parametrisation /// @param [in] freeTransportJacobian Transport jacobian free to free /// @param [in] freeToPathDerivatives Path length derivatives for free parameters -/// -/// @return the 6x8 transport jacobian from bound to free -FreeToBoundMatrix freeToBoundTransportJacobian( - const GeometryContext& geoContext, const Surface& surface, - const FreeVector& freeParameters, const FreeMatrix& freeTransportJacobian, - const FreeVector& freeToPathDerivatives); +/// @param [out] fullTransportJacobian The 6x8 transport jacobian from bound to free +/// +void freeToBoundTransportJacobian(const GeometryContext& geoContext, + const Surface& surface, + const FreeVector& freeParameters, + const FreeMatrix& freeTransportJacobian, + const FreeVector& freeToPathDerivatives, + FreeToBoundMatrix& fullTransportJacobian); /// @brief This function calculates the full transport jacobian from a free /// parameterisation to a bound one. Since a variation of the start diff --git a/Core/include/Acts/Propagator/detail/LoopStepperUtils.hpp b/Core/include/Acts/Propagator/detail/LoopStepperUtils.hpp index d683b9d76dc..55925db5e0b 100644 --- a/Core/include/Acts/Propagator/detail/LoopStepperUtils.hpp +++ b/Core/include/Acts/Propagator/detail/LoopStepperUtils.hpp @@ -126,11 +126,11 @@ struct LoopComponentProxy Result boundState( const Surface& surface, bool transportCov, const FreeToBoundCorrection& freeToBoundCorrection) { - return detail::boundState( - all_state.geoContext, cov(), jacobian(), jacTransport(), derivative(), - jacToGlobal(), pars(), all_state.particleHypothesis, - all_state.covTransport && transportCov, cmp.state.pathAccumulated, - surface, freeToBoundCorrection); + return detail::boundState(all_state.geoContext, surface, cov(), jacobian(), + jacTransport(), derivative(), jacToGlobal(), + pars(), all_state.particleHypothesis, + all_state.covTransport && transportCov, + cmp.state.pathAccumulated, freeToBoundCorrection); } void update(const FreeVector& freeParams, const BoundVector& boundParams, diff --git a/Core/src/Propagator/detail/JacobianEngine.cpp b/Core/src/Propagator/detail/JacobianEngine.cpp index 8c4512ea4a7..39fab57c4fd 100644 --- a/Core/src/Propagator/detail/JacobianEngine.cpp +++ b/Core/src/Propagator/detail/JacobianEngine.cpp @@ -143,19 +143,20 @@ BoundToFreeMatrix detail::boundToFreeTransportJacobian( return (freeTransportJacobian * boundToFreeJacobian); } -FreeToBoundMatrix detail::freeToBoundTransportJacobian( +void detail::freeToBoundTransportJacobian( const GeometryContext& geoContext, const Surface& surface, const FreeVector& freeParameters, const FreeMatrix& freeTransportJacobian, - const FreeVector& freeToPathDerivatives) { + const FreeVector& freeToPathDerivatives, + FreeToBoundMatrix& fullTransportJacobian) { // Calculate the jacobian from free to bound at the final surface FreeToBoundMatrix freeToBoundJacobian = surface.freeToBoundJacobian(geoContext, freeParameters); FreeToPathMatrix sVec = surface.freeToPathDerivative(geoContext, freeParameters); // Return the jacobian to local - return freeToBoundJacobian * - (freeTransportJacobian + - freeToPathDerivatives * sVec * freeTransportJacobian); + fullTransportJacobian = freeToBoundJacobian * (freeTransportJacobian + + freeToPathDerivatives * sVec * + freeTransportJacobian); } FreeToBoundMatrix detail::freeToCurvilinearTransportJacobian( diff --git a/Tests/UnitTests/Core/Propagator/CovarianceEngineTests.cpp b/Tests/UnitTests/Core/Propagator/CovarianceEngineTests.cpp index 8a2d1fdc3cf..d7a106fa59e 100644 --- a/Tests/UnitTests/Core/Propagator/CovarianceEngineTests.cpp +++ b/Tests/UnitTests/Core/Propagator/CovarianceEngineTests.cpp @@ -82,8 +82,8 @@ BOOST_AUTO_TEST_CASE(covariance_engine_test) { FreeToBoundCorrection freeToBoundCorrection(false); auto surface = Surface::makeShared(position, direction); detail::transportCovarianceToBound( - tgContext, covariance, jacobian, transportJacobian, derivatives, - boundToFreeJacobian, parameters, *surface, freeToBoundCorrection); + tgContext, *surface, covariance, jacobian, transportJacobian, derivatives, + boundToFreeJacobian, parameters, freeToBoundCorrection); BOOST_CHECK_NE(covariance, Covariance::Identity()); BOOST_CHECK_NE(jacobian, 2. * Jacobian::Identity()); @@ -121,9 +121,9 @@ BOOST_AUTO_TEST_CASE(covariance_engine_test) { // Produce a bound state without covariance matrix covarianceBefore = covariance; auto boundResult = - detail::boundState(tgContext, covariance, jacobian, transportJacobian, - derivatives, boundToFreeJacobian, parameters, - particleHypothesis, false, 1337., *surface, + detail::boundState(tgContext, *surface, covariance, jacobian, + transportJacobian, derivatives, boundToFreeJacobian, + parameters, particleHypothesis, false, 1337., freeToBoundCorrection) .value(); BOOST_CHECK(std::get<0>(curvResult).covariance().has_value()); @@ -139,9 +139,9 @@ BOOST_AUTO_TEST_CASE(covariance_engine_test) { // Produce a bound state with covariance matrix boundResult = - detail::boundState(tgContext, covariance, jacobian, transportJacobian, - derivatives, boundToFreeJacobian, parameters, - ParticleHypothesis::pion(), true, 1337., *surface, + detail::boundState(tgContext, *surface, covariance, jacobian, + transportJacobian, derivatives, boundToFreeJacobian, + parameters, ParticleHypothesis::pion(), true, 1337., freeToBoundCorrection) .value(); BOOST_CHECK(std::get<0>(boundResult).covariance().has_value()); @@ -155,9 +155,9 @@ BOOST_AUTO_TEST_CASE(covariance_engine_test) { // Produce a bound state with free to bound correction boundResult = - detail::boundState(tgContext, covariance, jacobian, transportJacobian, - derivatives, boundToFreeJacobian, parameters, - ParticleHypothesis::pion(), true, 1337., *surface, + detail::boundState(tgContext, *surface, covariance, jacobian, + transportJacobian, derivatives, boundToFreeJacobian, + parameters, ParticleHypothesis::pion(), true, 1337., freeToBoundCorrection) .value(); BOOST_CHECK(std::get<0>(boundResult).covariance().has_value()); diff --git a/Tests/UnitTests/Core/Propagator/JacobianEngineTests.cpp b/Tests/UnitTests/Core/Propagator/JacobianEngineTests.cpp index 5183df23c1a..3062658cadb 100644 --- a/Tests/UnitTests/Core/Propagator/JacobianEngineTests.cpp +++ b/Tests/UnitTests/Core/Propagator/JacobianEngineTests.cpp @@ -136,31 +136,33 @@ BOOST_AUTO_TEST_CASE(jacobian_engine_to_bound) { // (1) curvilinear/bound to bound transport jacobian // a) test without actual transport into the same surface - BoundMatrix b2bTransportJacobian = detail::boundToBoundTransportJacobian( - tgContext, freeParameters, boundToFreeJacobian, noTransportJacobian, - freeToPathDerivatives, *pSurface); + BoundMatrix b2bTransportJacobian; + detail::boundToBoundTransportJacobian( + tgContext, *pSurface, freeParameters, boundToFreeJacobian, + noTransportJacobian, freeToPathDerivatives, b2bTransportJacobian); BoundSquareMatrix newBoundCovariance = b2bTransportJacobian * boundCovariance * b2bTransportJacobian.transpose(); BOOST_CHECK(boundCovariance.isApprox(newBoundCovariance)); // b) test without actual transport but to a new surface - b2bTransportJacobian = detail::boundToBoundTransportJacobian( - tgContext, freeParameters, boundToFreeJacobian, noTransportJacobian, - freeToPathDerivatives, *oSurface); + detail::boundToBoundTransportJacobian( + tgContext, *oSurface, freeParameters, boundToFreeJacobian, + noTransportJacobian, freeToPathDerivatives, b2bTransportJacobian); newBoundCovariance = b2bTransportJacobian * boundCovariance * b2bTransportJacobian.transpose(); BOOST_CHECK(!boundCovariance.isApprox(newBoundCovariance)); // c) test to "the same" surface with transport // (not really senseful, but should give a different result) - b2bTransportJacobian = detail::boundToBoundTransportJacobian( - tgContext, freeParameters, boundToFreeJacobian, realTransportJacobian, - freeToPathDerivatives, *pSurface); + detail::boundToBoundTransportJacobian( + tgContext, *pSurface, freeParameters, boundToFreeJacobian, + realTransportJacobian, freeToPathDerivatives, b2bTransportJacobian); newBoundCovariance = b2bTransportJacobian * boundCovariance * b2bTransportJacobian.transpose(); BOOST_CHECK(!boundCovariance.isApprox(newBoundCovariance)); - FreeToBoundMatrix f2bTransportJacobian = detail::freeToBoundTransportJacobian( - tgContext, freeParameters, noTransportJacobian, freeToPathDerivatives, - *pSurface); + FreeToBoundMatrix f2bTransportJacobian; + detail::freeToBoundTransportJacobian( + tgContext, *pSurface, freeParameters, noTransportJacobian, + freeToPathDerivatives, f2bTransportJacobian); newBoundCovariance = f2bTransportJacobian * freeCovariance * f2bTransportJacobian.transpose(); @@ -206,17 +208,17 @@ BOOST_AUTO_TEST_CASE(jacobian_engine_to_curvilinear) { // (1) curvilinear/bound to curvilinear transport jacobian // a) test without actual transport into the same surface - BoundMatrix b2cTransportJacobian = - detail::boundToCurvilinearTransportJacobian( - direction, boundToFreeJacobian, noTransportJacobian, - freeToPathDerivatives); + BoundMatrix b2cTransportJacobian; + detail::boundToCurvilinearTransportJacobian( + direction, boundToFreeJacobian, noTransportJacobian, + freeToPathDerivatives, b2cTransportJacobian); BoundSquareMatrix newBoundCovariance = b2cTransportJacobian * boundCovariance * b2cTransportJacobian.transpose(); BOOST_CHECK(boundCovariance.isApprox(newBoundCovariance)); // b) test to another curvilinear frame at the same point (no transport) - b2cTransportJacobian = detail::boundToCurvilinearTransportJacobian( + detail::boundToCurvilinearTransportJacobian( Vector3(4., 5., 6.).normalized(), boundToFreeJacobian, - noTransportJacobian, freeToPathDerivatives); + noTransportJacobian, freeToPathDerivatives, b2cTransportJacobian); newBoundCovariance = b2cTransportJacobian * boundCovariance * b2cTransportJacobian.transpose(); BOOST_CHECK(!boundCovariance.isApprox(newBoundCovariance)); From 1cf045ddb01919b47850a433c755e97e38fdf4b9 Mon Sep 17 00:00:00 2001 From: andiwand Date: Mon, 11 Dec 2023 16:31:07 +0100 Subject: [PATCH 5/9] minor --- Core/src/Propagator/detail/JacobianEngine.cpp | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/Core/src/Propagator/detail/JacobianEngine.cpp b/Core/src/Propagator/detail/JacobianEngine.cpp index 39fab57c4fd..5c071b4fc41 100644 --- a/Core/src/Propagator/detail/JacobianEngine.cpp +++ b/Core/src/Propagator/detail/JacobianEngine.cpp @@ -140,7 +140,7 @@ BoundToFreeMatrix detail::boundToFreeTransportJacobian( const FreeMatrix& freeTransportJacobian) { // Calculate the full jacobian, in this case simple a product of // jacobian(transport in free) * jacobian(bound to free) - return (freeTransportJacobian * boundToFreeJacobian); + return freeTransportJacobian * boundToFreeJacobian; } void detail::freeToBoundTransportJacobian( @@ -175,9 +175,6 @@ Result detail::reinitializeJacobians( const GeometryContext& geoContext, const Surface& surface, FreeMatrix& freeTransportJacobian, FreeVector& freeToPathDerivatives, BoundToFreeMatrix& boundToFreeJacobian, const FreeVector& freeParameters) { - using VectorHelpers::phi; - using VectorHelpers::theta; - // Reset the jacobians freeTransportJacobian = FreeMatrix::Identity(); freeToPathDerivatives = FreeVector::Zero(); From 8408067d916b1c17a89194984e8ebaad35d4230d Mon Sep 17 00:00:00 2001 From: andiwand Date: Fri, 22 Dec 2023 16:07:16 +0100 Subject: [PATCH 6/9] fix after merging main --- Core/src/Propagator/detail/CovarianceEngine.cpp | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/Core/src/Propagator/detail/CovarianceEngine.cpp b/Core/src/Propagator/detail/CovarianceEngine.cpp index 9f59fad3830..321d706bb76 100644 --- a/Core/src/Propagator/detail/CovarianceEngine.cpp +++ b/Core/src/Propagator/detail/CovarianceEngine.cpp @@ -218,9 +218,10 @@ Acts::Result detail::boundToBoundConversion( freeToPathDerivatives.segment<3>(eFreeDir0) = bField.cross(freePars.segment<3>(eFreeDir0)); - BoundMatrix boundToBoundJac = detail::boundToBoundTransportJacobian( - gctx, freePars, boundToFreeJacobian, freeTransportJacobian, - freeToPathDerivatives, targetSurface); + BoundMatrix boundToBoundJac; + detail::boundToBoundTransportJacobian( + gctx, targetSurface, freePars, boundToFreeJacobian, + freeTransportJacobian, freeToPathDerivatives, boundToBoundJac); covOut = boundToBoundJac * (*boundParameters.covariance()) * boundToBoundJac.transpose(); From 2a2b0d1677ac5b81bdcfa26d2e6fb91485156737 Mon Sep 17 00:00:00 2001 From: Andreas Stefl Date: Tue, 16 Jan 2024 16:08:49 +0100 Subject: [PATCH 7/9] Update Core/include/Acts/Propagator/detail/JacobianEngine.hpp Co-authored-by: Paul Gessinger --- Core/include/Acts/Propagator/detail/JacobianEngine.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Core/include/Acts/Propagator/detail/JacobianEngine.hpp b/Core/include/Acts/Propagator/detail/JacobianEngine.hpp index d2642bc7dc2..5fe5448514c 100644 --- a/Core/include/Acts/Propagator/detail/JacobianEngine.hpp +++ b/Core/include/Acts/Propagator/detail/JacobianEngine.hpp @@ -105,7 +105,7 @@ BoundToFreeMatrix boundToFreeTransportJacobian( const FreeMatrix& freeTransportJacobian); /// @brief This function calculates the full jacobian from a given -/// free parameterisation to a new curvilinear bound parameterisation. +/// free parametrisation to a new curvilinear bound parametrisation. /// /// @note Modifications of the jacobian related to the /// projection onto the target surface is considered. Since a variation of From 2bd95b00c245ff7f08b350c168b270eaceb7c04f Mon Sep 17 00:00:00 2001 From: andiwand Date: Tue, 16 Jan 2024 16:37:36 +0100 Subject: [PATCH 8/9] improve doc --- .../Acts/Propagator/detail/JacobianEngine.hpp | 15 +++++++++++++-- 1 file changed, 13 insertions(+), 2 deletions(-) diff --git a/Core/include/Acts/Propagator/detail/JacobianEngine.hpp b/Core/include/Acts/Propagator/detail/JacobianEngine.hpp index 5fe5448514c..7beab2461cd 100644 --- a/Core/include/Acts/Propagator/detail/JacobianEngine.hpp +++ b/Core/include/Acts/Propagator/detail/JacobianEngine.hpp @@ -143,7 +143,12 @@ FreeToBoundMatrix freeToCurvilinearTransportJacobian( const FreeVector& freeToPathDerivatives); /// @brief This function reinitialises the state members required for the -/// covariance transport +/// covariance transport for usual surfaces +/// +/// Reinitialize jacobian components: +/// ->The transportJacobian is reinitialized to Identity +/// ->The derivatives is reinitialized to Zero +/// ->The boundToFreeJacobian is initialized to that at the current surface /// /// @param [in] geoContext The geometry context /// @param [in] surface The reference surface of the local parametrisation @@ -162,7 +167,13 @@ Result reinitializeJacobians(const GeometryContext& geoContext, const FreeVector& freeParameters); /// @brief This function reinitialises the state members required for the -/// covariance transport +/// covariance transport for curvilinear surfaces +/// +/// Reinitialize jacobian components: +/// ->The free transportJacobian is reinitialized to Identity +/// ->The path derivatives is reinitialized to Zero +/// ->The boundToFreeJacobian is reinitialized to that at the current +/// curvilinear surface /// /// @param [in, out] freeTransportJacobian The transport jacobian from start /// free to final free parameters From 8dd76754a47e44835cc08b67312ac3299f11749f Mon Sep 17 00:00:00 2001 From: andiwand Date: Tue, 16 Jan 2024 16:52:50 +0100 Subject: [PATCH 9/9] revive function which is used by Athena --- .../Acts/Propagator/detail/JacobianEngine.hpp | 44 +++++++++++++++---- Core/src/Propagator/detail/JacobianEngine.cpp | 12 +++++ 2 files changed, 47 insertions(+), 9 deletions(-) diff --git a/Core/include/Acts/Propagator/detail/JacobianEngine.hpp b/Core/include/Acts/Propagator/detail/JacobianEngine.hpp index 7beab2461cd..e6798262256 100644 --- a/Core/include/Acts/Propagator/detail/JacobianEngine.hpp +++ b/Core/include/Acts/Propagator/detail/JacobianEngine.hpp @@ -42,19 +42,18 @@ BoundToFreeMatrix curvilinearToFreeJacobian(const Vector3& direction); /// @brief This function calculates the full transport jacobian from a bound /// curvilinear representation to a new bound representation /// -/// @note Modifications of the jacobian related to the -/// projection onto a surface is considered. Since a variation of the start -/// parameters within a given uncertainty would lead to a variation of the end -/// parameters, these need to be propagated onto the target surface. This an -/// approximated approach to treat the (assumed) small change. +/// @note Modifications of the jacobian related to the projection onto a surface is +/// considered. Since a variation of the start parameters within a given +/// uncertainty would lead to a variation of the end parameters, these need to +/// be propagated onto the target surface. This an approximated approach to +/// treat the (assumed) small change. /// /// @param [in] geoContext The geometry Context /// @param [in] surface Target surface /// @param [in] freeParameters Free, nominal parametrisation /// @param [in] boundToFreeJacobian Jacobian from bound to free at start /// @param [in] freeTransportJacobian Transport jacobian free to free -/// @param [in] freeToPathDerivatives Path length derivatives for free -/// parameters +/// @param [in] freeToPathDerivatives Path length derivatives for free parameters /// @param [out] fullTransportJacobian A 6x6 transport jacobian from bound to bound /// /// @note jac(locA->locB) = jac(gloB->locB)*(1+ @@ -67,6 +66,34 @@ void boundToBoundTransportJacobian(const GeometryContext& geoContext, const FreeVector& freeToPathDerivatives, BoundMatrix& fullTransportJacobian); +/// TODO this should be removed but Athena depends on it +/// +/// @brief This function calculates the full transport jacobian from a bound +/// curvilinear representation to a new bound representation +/// +/// @note Modifications of the jacobian related to the projection onto a surface is +/// considered. Since a variation of the start parameters within a given +/// uncertainty would lead to a variation of the end parameters, these need to +/// be propagated onto the target surface. This an approximated approach to +/// treat the (assumed) small change. +/// +/// @param [in] geoContext The geometry Context +/// @param [in] freeParameters Free, nominal parametrisation +/// @param [in] boundToFreeJacobian Jacobian from bound to free at start +/// @param [in] freeTransportJacobian Transport jacobian free to free +/// @param [in] freeToPathDerivatives Path length derivatives for free parameters +/// @param [in] surface Target surface +/// +/// @note jac(locA->locB) = jac(gloB->locB)*(1+ +/// pathCorrectionFactor(gloB))*jacTransport(gloA->gloB) *jac(locA->gloA) +/// +/// @return a 6x6 transport jacobian from bound to bound +BoundMatrix boundToBoundTransportJacobian( + const GeometryContext& geoContext, const FreeVector& freeParameters, + const BoundToFreeMatrix& boundToFreeJacobian, + const FreeMatrix& freeTransportJacobian, + const FreeVector& freeToPathDerivatives, const Surface& surface); + /// @brief This function calculates the full jacobian from a given /// bound/curvilinear parameterisation from a surface to new curvilinear /// parameterisation. @@ -80,8 +107,7 @@ void boundToBoundTransportJacobian(const GeometryContext& geoContext, /// @param [in] direction Normalised direction vector /// @param [in] boundToFreeJacobian Jacobian from bound to free at start /// @param [in] freeTransportJacobian Transport jacobian free to free -/// @param [in] freeToPathDerivatives Path length derivatives for free -/// parameters +/// @param [in] freeToPathDerivatives Path length derivatives for free parameters /// @param [out] fullTransportJacobian A 6x6 transport jacobian from curilinear to bound /// /// @note The parameter @p surface is only required if projected to bound diff --git a/Core/src/Propagator/detail/JacobianEngine.cpp b/Core/src/Propagator/detail/JacobianEngine.cpp index 5c071b4fc41..1b98b070b70 100644 --- a/Core/src/Propagator/detail/JacobianEngine.cpp +++ b/Core/src/Propagator/detail/JacobianEngine.cpp @@ -112,6 +112,18 @@ void detail::boundToBoundTransportJacobian( freeTransportJacobian * boundToFreeJacobian; } +BoundMatrix detail::boundToBoundTransportJacobian( + const GeometryContext& geoContext, const FreeVector& freeParameters, + const BoundToFreeMatrix& boundToFreeJacobian, + const FreeMatrix& freeTransportJacobian, + const FreeVector& freeToPathDerivatives, const Surface& surface) { + BoundMatrix result; + detail::boundToBoundTransportJacobian( + geoContext, surface, freeParameters, boundToFreeJacobian, + freeTransportJacobian, freeToPathDerivatives, result); + return result; +} + void detail::boundToCurvilinearTransportJacobian( const Vector3& direction, const BoundToFreeMatrix& boundToFreeJacobian, const FreeMatrix& freeTransportJacobian,