refactor: Introduce CurvilinearSurface (acts-project#2814)

I noticed that the concept of a curvilinear surface is a quite spread over the codebase and I wanted to improve this. Here I introduce a surface like class which does not actually inherit from `Acts::Surface` but shares some of the functionality. This way the jacobians are put in a single place and can be used from somewhere else in an expressive fashion.

Afterwards it might make sense to let the `CurvilinearTrackParameters` depend on this special surface instead of a `PlaneSurface` or to merge them with `BoundTrackParameters` completely.

related issues
- acts-project#2812 

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

Core/include/Acts/EventData/GenericCurvilinearTrackParameters.hpp

@@ -10,7 +10,7 @@
 
 #include "Acts/EventData/GenericBoundTrackParameters.hpp"
 #include "Acts/EventData/TrackParametersConcept.hpp"
-#include "Acts/Surfaces/PlaneSurface.hpp"
+#include "Acts/Surfaces/CurvilinearSurface.hpp"
 
 namespace Acts {
 
@@ -46,7 +46,7 @@ class GenericCurvilinearTrackParameters
                                     Scalar qOverP,
                                     std::optional<CovarianceMatrix> cov,
                                     ParticleHypothesis particleHypothesis)
-      : Base(Surface::makeShared<PlaneSurface>(pos4.segment<3>(ePos0), dir),
+      : Base(CurvilinearSurface(pos4.segment<3>(ePos0), dir).surface(),
              transformFreeToCurvilinearParameters(pos4[eTime], dir, qOverP),
              std::move(cov), std::move(particleHypothesis)) {}
 
@@ -62,8 +62,9 @@ class GenericCurvilinearTrackParameters
                                     Scalar theta, Scalar qOverP,
                                     std::optional<CovarianceMatrix> cov,
                                     ParticleHypothesis particleHypothesis)
-      : Base(Surface::makeShared<PlaneSurface>(
-                 pos4.segment<3>(ePos0), makeDirectionFromPhiTheta(phi, theta)),
+      : Base(CurvilinearSurface(pos4.segment<3>(ePos0),
+                                makeDirectionFromPhiTheta(phi, theta))
+                 .surface(),
              transformFreeToCurvilinearParameters(pos4[eTime], phi, theta,
                                                   qOverP),
              std::move(cov), std::move(particleHypothesis)) {}

Core/include/Acts/EventData/detail/TestTrackState.hpp

@@ -12,6 +12,7 @@
 #include "Acts/EventData/VectorMultiTrajectory.hpp"
 #include "Acts/EventData/detail/GenerateParameters.hpp"
 #include "Acts/EventData/detail/TestSourceLink.hpp"
+#include "Acts/Surfaces/CurvilinearSurface.hpp"
 #include "Acts/Utilities/CalibrationContext.hpp"
 
 #include <random>
@@ -34,8 +35,8 @@ struct TestTrackState {
   // @param std::size_t nMeasurement either 1 or 2
   template <typename rng_t>
   TestTrackState(rng_t& rng, std::size_t measdim)
-      : surface(Surface::makeShared<PlaneSurface>(Vector3::Zero(),
-                                                  Vector3::UnitZ())),
+      : surface(
+            CurvilinearSurface(Vector3::Zero(), Vector3::UnitZ()).surface()),
         // set bogus parameters first since they are not default-constructible
         predicted(surface, BoundVector::Zero(), std::nullopt,
                   ParticleHypothesis::pion()),

Core/include/Acts/Propagator/detail/JacobianEngine.hpp

@@ -65,8 +65,6 @@ 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
 ///

Core/include/Acts/Surfaces/CurvilinearSurface.hpp

@@ -0,0 +1,120 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2023-2024 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 <iosfwd>
+#include <memory>
+#include <string>
+
+namespace Acts {
+
+class Surface;
+class PlaneSurface;
+
+/// @brief Utility class for curvilinear surfaces
+///
+class CurvilinearSurface final {
+ public:
+  explicit CurvilinearSurface(const Vector3& direction)
+      : m_direction{direction} {}
+
+  /// Constructor with direction vector
+  CurvilinearSurface(const Vector3& position, const Vector3& direction)
+      : m_position{position}, m_direction{direction} {}
+
+  /// Return method for the surface center by reference
+  /// @note the center is always recalculated in order to not keep a cache
+  /// @return center position by value
+  Vector3 center() const { return m_position; }
+
+  /// Return the surface normal at a given @p position and @p direction.
+  /// This method is fully generic, and valid for all surface types.
+  /// @note For some surface types, the @p direction is ignored, but
+  ///       it is **not safe** to pass in a zero vector!
+  /// @return The normal vector at the given position and direction
+  Vector3 normal() const { return m_direction; }
+
+  bool isStandardRepresentation() const;
+
+  /// Return method for the reference frame
+  /// This is the frame in which the covariance matrix is defined (specialized
+  /// by all surfaces)
+  ///
+  /// @return RotationMatrix3 which defines the three axes of the measurement
+  /// frame
+  RotationMatrix3 referenceFrame() const;
+
+  /// Return method for the surface Transform3 by reference
+  /// In case a detector element is associated the surface transform
+  /// is just forwarded to the detector element in order to keep the
+  /// (mis-)alignment cache cetrally handled
+  ///
+  /// @return the contextual transform
+  Transform3 transform() const;
+
+  /// Calculate the jacobian from local to global which the surface knows best,
+  /// hence the calculation is done here.
+  ///
+  /// @return Jacobian from local to global
+  BoundToFreeMatrix boundToFreeJacobian() const;
+
+  /// Calculate the jacobian from global to local which the surface knows best,
+  /// hence the calculation is done here.
+  ///
+  /// @return Jacobian from global to local
+  FreeToBoundMatrix freeToBoundJacobian() const;
+
+  /// Calculate the derivative of path length at the geometry constraint or
+  /// point-of-closest-approach w.r.t. free parameters. The calculation is
+  /// identical for all surfaces where the reference frame does not depend on
+  /// the direction
+  ///
+  /// @return Derivative of path length w.r.t. free parameters
+  FreeToPathMatrix freeToPathDerivative() const;
+
+  /// Output Method for std::ostream
+  ///
+  /// @param sl is the ostream to be dumped into
+  std::ostream& toStream(std::ostream& sl) const;
+
+  /// Output into a std::string
+  ///
+  /// @return the string representation of the curvilinear surface
+  std::string toString() const;
+
+  /// Return the plane surface representation of the curvilinear surface
+  ///
+  /// @return the plane surface representation of the curvilinear surface
+  std::shared_ptr<PlaneSurface> planeSurface() const;
+
+  /// Return the surface representation of the curvilinear surface
+  ///
+  /// @note same as planeSurface() but returns the base class.
+  ///   This is useful if the type of the surface is not relevant.
+  ///
+  /// @return the surface representation of the curvilinear surface
+  std::shared_ptr<Surface> surface() const;
+
+  /// Output operator
+  ///
+  /// @param os is the ostream to be dumped into
+  /// @param surface is the CurvilinearSurface to be dumped
+  /// @return the ostream
+  friend std::ostream& operator<<(std::ostream& os,
+                                  const CurvilinearSurface& surface);
+
+ private:
+  Vector3 m_position = Vector3::Zero();
+  Vector3 m_direction = Vector3::UnitZ();
+};
+
+}  // namespace Acts

Core/include/Acts/Surfaces/PlaneSurface.hpp

@@ -59,6 +59,8 @@ class PlaneSurface : public RegularSurface {
   PlaneSurface(const GeometryContext& gctx, const PlaneSurface& other,
                const Transform3& transform);
 
+  /// @deprecated Use `CurvilinearSurface` instead
+  ///
   /// Dedicated Constructor with normal vector
   /// This is for curvilinear surfaces which are by definition boundless
   ///

Core/src/Propagator/detail/JacobianEngine.cpp

@@ -9,83 +9,14 @@
 #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/CurvilinearSurface.hpp"
 #include "Acts/Surfaces/Surface.hpp"
 #include "Acts/Utilities/AlgebraHelpers.hpp"
-#include "Acts/Utilities/JacobianHelpers.hpp"
-#include "Acts/Utilities/VectorHelpers.hpp"
-
-#include <cmath>
 
 namespace Acts {
 
-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) {
-    // We normally operate in curvilinear coordinates defined as follows
-    freeToCurvJacobian(eBoundLoc0, eFreePos0) = -sinPhi;
-    freeToCurvJacobian(eBoundLoc0, eFreePos1) = cosPhi;
-    freeToCurvJacobian(eBoundLoc1, eFreePos0) = -cosPhi * cosTheta;
-    freeToCurvJacobian(eBoundLoc1, eFreePos1) = -sinPhi * cosTheta;
-    freeToCurvJacobian(eBoundLoc1, eFreePos2) = sinTheta;
-  } else {
-    // Under grazing incidence to z, the above coordinate system definition
-    // becomes numerically unstable, and we need to switch to another one
-    const ActsScalar x = direction(0);  // == cos(phi) * sin(theta)
-    const ActsScalar y = direction(1);  // == sin(phi) * sin(theta)
-    const ActsScalar z = direction(2);  // == cos(theta)
-    const ActsScalar c = std::hypot(y, z);
-    const ActsScalar invC = 1. / c;
-    freeToCurvJacobian(eBoundLoc0, eFreePos1) = -z * invC;
-    freeToCurvJacobian(eBoundLoc0, eFreePos2) = y * invC;
-    freeToCurvJacobian(eBoundLoc1, eFreePos0) = c;
-    freeToCurvJacobian(eBoundLoc1, eFreePos1) = -x * y * invC;
-    freeToCurvJacobian(eBoundLoc1, eFreePos2) = -x * z * invC;
-  }
-  // Time parameter
-  freeToCurvJacobian(eBoundTime, eFreeTime) = 1.;
-  // Directional and momentum parameters for curvilinear
-  freeToCurvJacobian(eBoundPhi, eFreeDir0) = -sinPhi * invSinTheta;
-  freeToCurvJacobian(eBoundPhi, eFreeDir1) = cosPhi * invSinTheta;
-  freeToCurvJacobian(eBoundTheta, eFreeDir0) = cosPhi * cosTheta;
-  freeToCurvJacobian(eBoundTheta, eFreeDir1) = sinPhi * cosTheta;
-  freeToCurvJacobian(eBoundTheta, eFreeDir2) = -sinTheta;
-  freeToCurvJacobian(eBoundQOverP, eFreeQOverP) = 1.;
-
-  return freeToCurvJacobian;
-}
-
-BoundToFreeMatrix detail::curvilinearToFreeJacobian(const Vector3& direction) {
-  auto [cosPhi, sinPhi, cosTheta, sinTheta] =
-      VectorHelpers::evaluateTrigonomics(direction);
-
-  // Prepare the jacobian to free
-  BoundToFreeMatrix curvToFreeJacobian = BoundToFreeMatrix::Zero();
-
-  curvToFreeJacobian(eFreePos0, eBoundLoc0) = -sinPhi;
-  curvToFreeJacobian(eFreePos0, eBoundLoc1) = -cosPhi * cosTheta;
-  curvToFreeJacobian(eFreePos1, eBoundLoc0) = cosPhi;
-  curvToFreeJacobian(eFreePos1, eBoundLoc1) = -sinPhi * cosTheta;
-  curvToFreeJacobian(eFreePos2, eBoundLoc1) = sinTheta;
-  // Time parameter: stays as is
-  curvToFreeJacobian(eFreeTime, eBoundTime) = 1;
-  curvToFreeJacobian(eFreeDir0, eBoundPhi) = -sinTheta * sinPhi;
-  curvToFreeJacobian(eFreeDir0, eBoundTheta) = cosTheta * cosPhi;
-  curvToFreeJacobian(eFreeDir1, eBoundPhi) = sinTheta * cosPhi;
-  curvToFreeJacobian(eFreeDir1, eBoundTheta) = cosTheta * sinPhi;
-  curvToFreeJacobian(eFreeDir2, eBoundTheta) = -sinTheta;
-  // Q/P parameter: stays as is
-  curvToFreeJacobian(eFreeQOverP, eBoundQOverP) = 1;
-
-  return curvToFreeJacobian;
-}
-
 void detail::boundToBoundTransportJacobian(
     const GeometryContext& geoContext, const Surface& surface,
     const FreeVector& freeParameters,
@@ -131,7 +62,8 @@ void detail::boundToCurvilinearTransportJacobian(
     const FreeVector& freeToPathDerivatives,
     BoundMatrix& fullTransportJacobian) {
   // Calculate the jacobian from global to local at the curvilinear surface
-  FreeToBoundMatrix freeToBoundJacobian = freeToCurvilinearJacobian(direction);
+  FreeToBoundMatrix freeToBoundJacobian =
+      CurvilinearSurface(direction).freeToBoundJacobian();
 
   // Update the jacobian to include the derivative of the path length at the
   // curvilinear surface w.r.t. the free parameters
@@ -182,7 +114,7 @@ FreeToBoundMatrix detail::freeToCurvilinearTransportJacobian(
 
   // Since the jacobian to local needs to calculated for the bound parameters
   // here, it is convenient to do the same here
-  return freeToCurvilinearJacobian(direction) *
+  return CurvilinearSurface(direction).freeToBoundJacobian() *
          (freeTransportJacobian - freeToPathDerivatives * sfactors);
 }
 
@@ -214,20 +146,7 @@ void detail::reinitializeJacobians(FreeMatrix& freeTransportJacobian,
   // 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;
+  boundToFreeJacobian = CurvilinearSurface(direction).boundToFreeJacobian();
 }
 
 }  // namespace Acts

Core/src/Surfaces/CMakeLists.txt

@@ -26,4 +26,5 @@ target_sources(
     detail/AlignmentHelper.cpp
     VerticesHelper.cpp
     RegularSurface.cpp
+    CurvilinearSurface.cpp
 )