refactor!: Refactor Surface::boundToFreeJacobian to use free instea…

…d of bound vector (acts-project#2811)

This avoid free->bound->free roundtrips and aligns the interface with `freeToBoundJacobian`

related issues
- acts-project#2810

blocked by
- acts-project#2789
- acts-project#2782
- acts-project#2781

Core/include/Acts/Propagator/EigenStepper.hpp

@@ -92,7 +92,7 @@ class EigenStepper {
         // set the covariance transport flag to true and copy
         covTransport = true;
         cov = BoundSquareMatrix(*par.covariance());
-        jacToGlobal = surface.boundToFreeJacobian(gctx, par.parameters());
+        jacToGlobal = surface.boundToFreeJacobian(gctx, pars);
       }
     }
 

Core/include/Acts/Propagator/EigenStepper.ipp

@@ -29,17 +29,17 @@ void Acts::EigenStepper<E, A>::resetState(State& state,
                                           const BoundSquareMatrix& cov,
                                           const Surface& surface,
                                           const double stepSize) const {
+  FreeVector freeParams = detail::transformBoundToFreeParameters(
+      surface, state.geoContext, boundParams);
+
   // Update the stepping state
-  update(state,
-         detail::transformBoundToFreeParameters(surface, state.geoContext,
-                                                boundParams),
-         boundParams, cov, surface);
+  state.pars = freeParams;
+  state.cov = cov;
   state.stepSize = ConstrainedStep(stepSize);
   state.pathAccumulated = 0.;
 
   // Reinitialize the stepping jacobian
-  state.jacToGlobal =
-      surface.boundToFreeJacobian(state.geoContext, boundParams);
+  state.jacToGlobal = surface.boundToFreeJacobian(state.geoContext, freeParams);
   state.jacobian = BoundMatrix::Identity();
   state.jacTransport = FreeMatrix::Identity();
   state.derivative = FreeVector::Zero();
@@ -107,13 +107,12 @@ auto Acts::EigenStepper<E, A>::curvilinearState(State& state,
 template <typename E, typename A>
 void Acts::EigenStepper<E, A>::update(State& state,
                                       const FreeVector& freeParams,
-                                      const BoundVector& boundParams,
+                                      const BoundVector& /*boundParams*/,
                                       const Covariance& covariance,
                                       const Surface& surface) const {
   state.pars = freeParams;
   state.cov = covariance;
-  state.jacToGlobal =
-      surface.boundToFreeJacobian(state.geoContext, boundParams);
+  state.jacToGlobal = surface.boundToFreeJacobian(state.geoContext, freeParams);
 }
 
 template <typename E, typename A>

Core/include/Acts/Propagator/StraightLineStepper.hpp

@@ -86,7 +86,7 @@ class StraightLineStepper {
         // set the covariance transport flag to true and copy
         covTransport = true;
         cov = BoundSquareMatrix(*par.covariance());
-        jacToGlobal = surface.boundToFreeJacobian(gctx, par.parameters());
+        jacToGlobal = surface.boundToFreeJacobian(gctx, pars);
       }
     }
 

Core/include/Acts/Surfaces/DiscSurface.hpp

@@ -234,17 +234,17 @@ class DiscSurface : public RegularSurface {
   /// hence the calculation is done here.
   ///
   /// @param gctx The current geometry context object, e.g. alignment
-  /// @param boundParams is the bound parameters vector
+  /// @param parameters is the free parameters vector
   ///
   /// @return Jacobian from local to global
   BoundToFreeMatrix boundToFreeJacobian(
-      const GeometryContext& gctx, const BoundVector& boundParams) const final;
+      const GeometryContext& gctx, const FreeVector& parameters) const final;
 
   /// Calculate the jacobian from global to local which the surface knows best,
   /// hence the calculation is done here.
   ///
   /// @param gctx The current geometry context object, e.g. alignment
-  /// @param parameters is the free parameters
+  /// @param parameters is the free parameters vector
   ///
   /// @return Jacobian from global to local
   FreeToBoundMatrix freeToBoundJacobian(

Core/include/Acts/Surfaces/LineSurface.hpp

@@ -121,17 +121,17 @@ class LineSurface : public Surface {
   /// hence the calculation is done here.
   ///
   /// @param gctx The current geometry context object, e.g. alignment
-  /// @param boundParams is the bound parameters vector
+  /// @param parameters is the free parameters vector
   ///
   /// @return Jacobian from local to global
   BoundToFreeMatrix boundToFreeJacobian(
-      const GeometryContext& gctx, const BoundVector& boundParams) const final;
+      const GeometryContext& gctx, const FreeVector& parameters) const final;
 
   /// Calculate the derivative of path length at the geometry constraint or
   /// point-of-closest-approach w.r.t. free parameters
   ///
   /// @param gctx The current geometry context object, e.g. alignment
-  /// @param parameters is the free parameters
+  /// @param parameters is the free parameters vector
   ///
   /// @return Derivative of path length w.r.t. free parameters
   FreeToPathMatrix freeToPathDerivative(

Core/include/Acts/Surfaces/Surface.hpp

@@ -326,11 +326,11 @@ class Surface : public virtual GeometryObject,
   /// "Acts/EventData/detail/coordinate_transformations.hpp"
   ///
   /// @param gctx The current geometry context object, e.g. alignment
-  /// @param boundParams is the bound parameters vector
+  /// @param parameters is the free parameters vector
   ///
   /// @return Jacobian from local to global
   virtual BoundToFreeMatrix boundToFreeJacobian(
-      const GeometryContext& gctx, const BoundVector& boundParams) const;
+      const GeometryContext& gctx, const FreeVector& parameters) const;
 
   /// Calculate the jacobian from global to local which the surface knows best,
   /// hence the calculation is done here.
@@ -343,7 +343,7 @@ class Surface : public virtual GeometryObject,
   /// "Acts/EventData/detail/coordinate_transformations.hpp"
   ///
   /// @param gctx The current geometry context object, e.g. alignment
-  /// @param parameters is the free parameters
+  /// @param parameters is the free parameters vector
   ///
   /// @return Jacobian from global to local
   virtual FreeToBoundMatrix freeToBoundJacobian(

Core/include/Acts/Surfaces/SurfaceConcept.hpp

@@ -68,11 +68,7 @@ concept SurfaceConcept = requires(S s, const S cs, S s2, const S cs2,
     } -> std::same_as<RotationMatrix3>;
 
   {
-    cs.boundToFreeJacobian(gctx, BoundVector{})
-    } -> std::same_as<BoundToFreeMatrix>;
-
-  {
-    cs.boundToFreeJacobian(gctx, BoundVector{})
+    cs.boundToFreeJacobian(gctx, FreeVector{})
     } -> std::same_as<BoundToFreeMatrix>;
 
   {

Core/include/Acts/TrackFitting/detail/GsfActor.hpp

@@ -529,7 +529,8 @@ struct GsfActor {
       }
 
       auto& cmp = *res;
-      cmp.jacToGlobal() = surface.boundToFreeJacobian(state.geoContext, pars);
+      cmp.jacToGlobal() =
+          surface.boundToFreeJacobian(state.geoContext, cmp.pars());
       cmp.pathAccumulated() = state.stepping.pathAccumulated;
       cmp.jacobian() = Acts::BoundMatrix::Identity();
       cmp.derivative() = Acts::FreeVector::Zero();

Core/src/Propagator/StraightLineStepper.cpp

@@ -33,13 +33,12 @@ StraightLineStepper::curvilinearState(State& state, bool transportCov) const {
 }
 
 void StraightLineStepper::update(State& state, const FreeVector& freeParams,
-                                 const BoundVector& boundParams,
+                                 const BoundVector& /*boundParams*/,
                                  const Covariance& covariance,
                                  const Surface& surface) const {
   state.pars = freeParams;
   state.cov = covariance;
-  state.jacToGlobal =
-      surface.boundToFreeJacobian(state.geoContext, boundParams);
+  state.jacToGlobal = surface.boundToFreeJacobian(state.geoContext, freeParams);
 }
 
 void StraightLineStepper::update(State& state, const Vector3& uposition,
@@ -70,17 +69,17 @@ void StraightLineStepper::resetState(State& state,
                                      const BoundSquareMatrix& cov,
                                      const Surface& surface,
                                      const double stepSize) const {
+  FreeVector freeParams = detail::transformBoundToFreeParameters(
+      surface, state.geoContext, boundParams);
+
   // Update the stepping state
-  update(state,
-         detail::transformBoundToFreeParameters(surface, state.geoContext,
-                                                boundParams),
-         boundParams, cov, surface);
+  state.pars = freeParams;
+  state.cov = cov;
   state.stepSize = ConstrainedStep(stepSize);
   state.pathAccumulated = 0.;
 
   // Reinitialize the stepping jacobian
-  state.jacToGlobal =
-      surface.boundToFreeJacobian(state.geoContext, boundParams);
+  state.jacToGlobal = surface.boundToFreeJacobian(state.geoContext, freeParams);
   state.jacobian = BoundMatrix::Identity();
   state.jacTransport = FreeMatrix::Identity();
   state.derivative = FreeVector::Zero();

Core/src/Propagator/detail/CovarianceEngine.cpp

@@ -209,7 +209,7 @@ Acts::Result<Acts::BoundTrackParameters> detail::boundToBoundConversion(
 
   if (boundParameters.covariance().has_value()) {
     Acts::BoundToFreeMatrix boundToFreeJacobian =
-        sourceSurface.boundToFreeJacobian(gctx, boundParameters.parameters());
+        sourceSurface.boundToFreeJacobian(gctx, freePars);
 
     Acts::FreeMatrix freeTransportJacobian = FreeMatrix::Identity();
 

Core/src/Propagator/detail/JacobianEngine.cpp

@@ -198,15 +198,8 @@ Result<void> detail::reinitializeJacobians(
   if (!lpResult.ok()) {
     return lpResult.error();
   }
-  // Transform from free to bound parameters
-  Result<BoundVector> boundParameters = detail::transformFreeToBoundParameters(
-      freeParameters, surface, geoContext);
-  if (!boundParameters.ok()) {
-    return boundParameters.error();
-  }
   // Reset the jacobian from local to global
-  boundToFreeJacobian =
-      surface.boundToFreeJacobian(geoContext, *boundParameters);
+  boundToFreeJacobian = surface.boundToFreeJacobian(geoContext, freeParameters);
   return Result<void>::success();
 }
 

Core/src/Surfaces/DiscSurface.cpp

@@ -19,12 +19,10 @@
 #include "Acts/Surfaces/SurfaceError.hpp"
 #include "Acts/Surfaces/detail/FacesHelper.hpp"
 #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 <algorithm>
 #include <cmath>
 #include <stdexcept>
 #include <utility>
@@ -211,32 +209,32 @@ Acts::Vector2 Acts::DiscSurface::localCartesianToPolar(
 }
 
 Acts::BoundToFreeMatrix Acts::DiscSurface::boundToFreeJacobian(
-    const GeometryContext& gctx, const BoundVector& boundParams) const {
-  // Transform from bound to free parameters
-  FreeVector freeParams =
-      detail::transformBoundToFreeParameters(*this, gctx, boundParams);
+    const GeometryContext& gctx, const FreeVector& parameters) const {
   // The global position
-  const Vector3 position = freeParams.segment<3>(eFreePos0);
+  const Vector3 position = parameters.segment<3>(eFreePos0);
   // The direction
-  const Vector3 direction = freeParams.segment<3>(eFreeDir0);
+  const Vector3 direction = parameters.segment<3>(eFreeDir0);
 
   assert(isOnSurface(gctx, position, direction, BoundaryCheck(false)));
 
-  // special polar coordinates for the Disc
-  double lrad = boundParams[eBoundLoc0];
-  double lphi = boundParams[eBoundLoc1];
-  double lcphi = std::cos(lphi);
-  double lsphi = std::sin(lphi);
-  // retrieve the reference frame
-  const auto rframe = referenceFrame(gctx, position, direction);
+  // The measurement frame of the surface
+  RotationMatrix3 rframeT =
+      referenceFrame(gctx, position, direction).transpose();
+  // calculate the transformation to local coordinates
+  const Vector3 posLoc = transform(gctx).inverse() * position;
+  const double lr = perp(posLoc);
+  const double lphi = phi(posLoc);
+  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);
   // 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) =
-      lcphi * rframe.block<3, 1>(0, 0) + lsphi * rframe.block<3, 1>(0, 1);
+  jacToGlobal.block<3, 1>(eFreePos0, eBoundLoc0) = lcphi * lx + lsphi * ly;
   jacToGlobal.block<3, 1>(eFreePos0, eBoundLoc1) =
-      lrad *
-      (lcphi * rframe.block<3, 1>(0, 1) - lsphi * rframe.block<3, 1>(0, 0));
+      lr * (lcphi * ly - lsphi * lx);
   // the time component
   jacToGlobal(eFreeTime, eBoundTime) = 1;
   // the momentum components
@@ -250,6 +248,7 @@ Acts::FreeToBoundMatrix Acts::DiscSurface::freeToBoundJacobian(
     const GeometryContext& gctx, const FreeVector& parameters) const {
   using VectorHelpers::perp;
   using VectorHelpers::phi;
+
   // The global position
   const auto position = parameters.segment<3>(eFreePos0);
   // The direction
@@ -258,9 +257,9 @@ Acts::FreeToBoundMatrix Acts::DiscSurface::freeToBoundJacobian(
   RotationMatrix3 rframeT =
       referenceFrame(gctx, position, direction).transpose();
   // calculate the transformation to local coordinates
-  const Vector3 pos_loc = transform(gctx).inverse() * position;
-  const double lr = perp(pos_loc);
-  const double lphi = phi(pos_loc);
+  const Vector3 posLoc = transform(gctx).inverse() * position;
+  const double lr = perp(posLoc);
+  const double lphi = phi(posLoc);
   const double lcphi = std::cos(lphi);
   const double lsphi = std::sin(lphi);
   // rotate into the polar coorindates

Core/src/Surfaces/LineSurface.cpp

@@ -197,30 +197,22 @@ Acts::SurfaceMultiIntersection Acts::LineSurface::intersect(
 }
 
 Acts::BoundToFreeMatrix Acts::LineSurface::boundToFreeJacobian(
-    const GeometryContext& gctx, const BoundVector& boundParams) const {
-  // Transform from bound to free parameters
-  FreeVector freeParams =
-      detail::transformBoundToFreeParameters(*this, gctx, boundParams);
+    const GeometryContext& gctx, const FreeVector& parameters) const {
   // The global position
-  Vector3 position = freeParams.segment<3>(eFreePos0);
+  Vector3 position = parameters.segment<3>(eFreePos0);
   // The direction
-  Vector3 direction = freeParams.segment<3>(eFreeDir0);
-  // retrieve the reference frame
-  auto rframe = referenceFrame(gctx, position, direction);
+  Vector3 direction = parameters.segment<3>(eFreeDir0);
 
   assert(isOnSurface(gctx, position, direction, BoundaryCheck(false)));
 
-  // Initialize the jacobian from local to global
-  BoundToFreeMatrix jacToGlobal = BoundToFreeMatrix::Zero();
+  // retrieve the reference frame
+  auto rframe = referenceFrame(gctx, position, direction);
+
+  Vector2 local = *globalToLocal(gctx, position, direction,
+                                 std::numeric_limits<double>::max());
 
-  // the local error components - given by the reference frame
-  jacToGlobal.topLeftCorner<3, 2>() = rframe.topLeftCorner<3, 2>();
-  // the time component
-  jacToGlobal(eFreeTime, eBoundTime) = 1;
-  // the momentum components
-  jacToGlobal.block<3, 2>(eFreeDir0, eBoundPhi) =
-      sphericalToFreeDirectionJacobian(direction);
-  jacToGlobal(eFreeQOverP, eBoundQOverP) = 1;
+  BoundToFreeMatrix jacToGlobal =
+      Surface::boundToFreeJacobian(gctx, parameters);
 
   // For the derivative of global position with bound angles, refer the
   // following white paper:
@@ -239,10 +231,8 @@ Acts::BoundToFreeMatrix Acts::LineSurface::boundToFreeJacobian(
   dDThetaY -=
       rframe.block<3, 1>(0, 0) * (rframe.block<3, 1>(0, 0).dot(dDThetaY));
   // set the jacobian components for global d/ phi/Theta
-  jacToGlobal.block<3, 1>(eFreePos0, eBoundPhi) =
-      dDPhiY * boundParams[eBoundLoc0] * ipdn;
-  jacToGlobal.block<3, 1>(eFreePos0, eBoundTheta) =
-      dDThetaY * boundParams[eBoundLoc0] * ipdn;
+  jacToGlobal.block<3, 1>(eFreePos0, eBoundPhi) = dDPhiY * local.x() * ipdn;
+  jacToGlobal.block<3, 1>(eFreePos0, eBoundTheta) = dDThetaY * local.x() * ipdn;
 
   return jacToGlobal;
 }

Core/src/Surfaces/Surface.cpp

@@ -272,14 +272,10 @@ Acts::RotationMatrix3 Acts::Surface::referenceFrame(
 }
 
 Acts::BoundToFreeMatrix Acts::Surface::boundToFreeJacobian(
-    const GeometryContext& gctx, const BoundVector& boundParams) const {
-  // Transform from bound to free parameters
-  FreeVector freeParams =
-      detail::transformBoundToFreeParameters(*this, gctx, boundParams);
-  // The global position
-  const Vector3 position = freeParams.segment<3>(eFreePos0);
+    const GeometryContext& gctx, const FreeVector& parameters) const {
+  const Vector3 position = parameters.segment<3>(eFreePos0);
   // The direction
-  const Vector3 direction = freeParams.segment<3>(eFreeDir0);
+  const Vector3 direction = parameters.segment<3>(eFreeDir0);
   // retrieve the reference frame
   const auto rframe = referenceFrame(gctx, position, direction);