refactor!: ImpactPointEstimator moves to cpp file (acts-project#2971)

This PR moves the code of the `ImpactPointEstimator` from the header into a compiled file. To allow this, I have to make some changes to the interface to remove templates.

In the process, I introduced `Eigen::Map<ActsDynamicVector>` as an argument to `getDistanceAndMomentum` and `getVertexCompatibility`.
Part of:
- acts-project#2842 

Blocked by:
- acts-project#2953

Core/include/Acts/Vertexing/AdaptiveMultiVertexFitter.hpp

@@ -40,7 +40,7 @@ class AdaptiveMultiVertexFitter {
   struct State {
     State(const MagneticFieldProvider& field,
           const Acts::MagneticFieldContext& magContext)
-        : ipState(field.makeCache(magContext)),
+        : ipState{field.makeCache(magContext)},
           fieldCache(field.makeCache(magContext)) {}
     // Vertex collection to be fitted
     std::vector<Vertex*> vertexCollection;

Core/include/Acts/Vertexing/ImpactPointEstimator.hpp

@@ -8,6 +8,7 @@
 
 #pragma once
 
+#include "Acts/Definitions/Algebra.hpp"
 #include "Acts/EventData/TrackParameters.hpp"
 #include "Acts/Geometry/GeometryContext.hpp"
 #include "Acts/MagneticField/MagneticFieldContext.hpp"
@@ -47,11 +48,6 @@ class ImpactPointEstimator {
  public:
   /// State struct
   struct State {
-    /// @brief The state constructor
-    ///
-    /// @param fieldCacheIn The magnetic field cache
-    State(MagneticFieldProvider::Cache fieldCacheIn)
-        : fieldCache(std::move(fieldCacheIn)) {}
     /// Magnetic field cache
     MagneticFieldProvider::Cache fieldCache;
   };
@@ -139,18 +135,22 @@ class ImpactPointEstimator {
   /// compatible.
   ///
   /// @tparam nDim Number of dimensions used to compute compatibility
-  /// @note If nDim = 3 we only consider spatial dimensions; if nDim = 4, we
-  /// also consider time. Other values are not allowed.
+  /// @note If @p nDim = 3 we only consider spatial dimensions; if nDim = 4, we
+  ///       also consider time. Other values are not allowed.
   /// @param gctx The Geometry context
   /// @param trkParams Track parameters at point of closest
   /// approach in 3D as retrieved by estimate3DImpactParameters
   /// @param vertexPos The vertex position
   ///
   /// @return The compatibility value
-  template <unsigned int nDim>
+  template <int nDim>
   Result<double> getVertexCompatibility(
       const GeometryContext& gctx, const BoundTrackParameters* trkParams,
-      const ActsVector<nDim>& vertexPos) const;
+      const ActsVector<nDim>& vertexPos) const {
+    static_assert(nDim == 3 || nDim == 4,
+                  "Only 3D and 4D vertex positions allowed");
+    return getVertexCompatibility(gctx, trkParams, {vertexPos.data(), nDim});
+  }
 
   /// @brief Calculate the distance between a track and a vertex by finding the
   /// corresponding 3D PCA. Returns also the momentum direction at the 3D PCA.
@@ -160,17 +160,27 @@ class ImpactPointEstimator {
   /// tracks we use the Newton method.
   ///
   /// @tparam nDim Number of dimensions used to compute compatibility
-  /// @note If nDim = 3 we only consider spatial dimensions; if nDim = 4, we
-  /// also consider time. Other values are not allowed.
+  /// @note If @p nDim = 3 we only consider spatial dimensions; if nDim = 4, we
+  ///       also consider time. Other values are not allowed.
   /// @param gctx Geometry context
   /// @param trkParams Track parameters
   /// @param vtxPos Vertex position
   /// @param state The state object
-  template <unsigned int nDim>
+  template <int nDim>
   Result<std::pair<Acts::ActsVector<nDim>, Acts::Vector3>>
   getDistanceAndMomentum(const GeometryContext& gctx,
                          const BoundTrackParameters& trkParams,
-                         const ActsVector<nDim>& vtxPos, State& state) const;
+                         const ActsVector<nDim>& vtxPos, State& state) const {
+    static_assert(nDim == 3 || nDim == 4,
+                  "Only 3D and 4D vertex positions allowed");
+    auto res =
+        getDistanceAndMomentum(gctx, trkParams, {vtxPos.data(), nDim}, state);
+    if (!res.ok()) {
+      return res.error();
+    }
+    auto& [distance, momentum] = *res;
+    return std::pair{distance.template head<nDim>(), momentum};
+  }
 
   /// @brief Calculates the impact parameters of a track w.r.t. a vertex. The
   /// corresponding errors are approximated by summing the variances of the
@@ -219,6 +229,14 @@ class ImpactPointEstimator {
       const MagneticFieldContext& mctx) const;
 
  private:
+  Result<std::pair<Acts::Vector4, Acts::Vector3>> getDistanceAndMomentum(
+      const GeometryContext& gctx, const BoundTrackParameters& trkParams,
+      Eigen::Map<const ActsDynamicVector> vtxPos, State& state) const;
+
+  Result<double> getVertexCompatibility(
+      const GeometryContext& gctx, const BoundTrackParameters* trkParams,
+      Eigen::Map<const ActsDynamicVector> vertexPos) const;
+
   /// Configuration object
   const Config m_cfg;
 
@@ -227,24 +245,6 @@ class ImpactPointEstimator {
 
   /// Private access to logging instance
   const Logger& logger() const { return *m_logger; }
-
-  /// @brief Performs a Newton approximation to retrieve a point
-  /// of closest approach in 3D to a reference position
-  ///
-  /// @param helixCenter Position of the helix center
-  /// @param vtxPos Vertex position
-  /// @param phi Azimuthal momentum angle
-  /// @note Modifying phi corresponds to moving along the track. This function
-  /// optimizes phi until we reach a 3D PCA.
-  /// @param theta Polar momentum angle (constant along the track)
-  /// @param rho Signed helix radius
-  ///
-  /// @return Phi value at 3D PCA
-  Result<double> performNewtonOptimization(const Vector3& helixCenter,
-                                           const Vector3& vtxPos, double phi,
-                                           double theta, double rho) const;
 };
 
 }  // namespace Acts
-
-#include "Acts/Vertexing/ImpactPointEstimator.ipp"

Core/include/Acts/Vertexing/IterativeVertexFinder.hpp

@@ -121,7 +121,7 @@ class IterativeVertexFinder final : public IVertexFinder {
     State(const MagneticFieldProvider& field,
           const Acts::MagneticFieldContext& _magContext)
         : magContext(_magContext),
-          ipState(field.makeCache(magContext)),
+          ipState{field.makeCache(magContext)},
           fieldCache(field.makeCache(magContext)) {}
 
     std::reference_wrapper<const Acts::MagneticFieldContext> magContext;

Core/include/Acts/Vertexing/VertexingError.hpp

@@ -21,6 +21,7 @@ enum class VertexingError {
   NotConverged,
   ElementNotFound,
   NoCovariance,
+  InvalidInput,
 };
 
 std::error_code make_error_code(Acts::VertexingError e);

Core/src/Vertexing/AdaptiveMultiVertexFitter.cpp

@@ -8,6 +8,7 @@
 
 #include "Acts/Vertexing/AdaptiveMultiVertexFitter.hpp"
 
+#include "Acts/Surfaces/PerigeeSurface.hpp"
 #include "Acts/Vertexing/KalmanVertexUpdater.hpp"
 #include "Acts/Vertexing/VertexingError.hpp"
 
@@ -236,14 +237,17 @@ Acts::Result<void> Acts::AdaptiveMultiVertexFitter::setAllVertexCompatibilities(
     // Set compatibility with current vertex
     Acts::Result<double> compatibilityResult(0.);
     if (m_cfg.useTime) {
-      compatibilityResult = m_cfg.ipEst.template getVertexCompatibility<4>(
+      compatibilityResult = m_cfg.ipEst.getVertexCompatibility(
           vertexingOptions.geoContext, &(vtxInfo.impactParams3D.at(trk)),
           vtxInfo.oldPosition);
     } else {
-      compatibilityResult = m_cfg.ipEst.template getVertexCompatibility<3>(
+      Acts::Vector3 vertexPosOnly =
+          VectorHelpers::position(vtxInfo.oldPosition);
+      compatibilityResult = m_cfg.ipEst.getVertexCompatibility(
           vertexingOptions.geoContext, &(vtxInfo.impactParams3D.at(trk)),
-          VectorHelpers::position(vtxInfo.oldPosition));
+          vertexPosOnly);
     }
+
     if (!compatibilityResult.ok()) {
       return compatibilityResult.error();
     }
@@ -288,9 +292,10 @@ Acts::Result<void> Acts::AdaptiveMultiVertexFitter::setWeightsAndUpdate(
           trkAtVtx.linearizedState = *result;
           trkAtVtx.isLinearized = true;
         }
-        // Update the vertex with the new track. The second template argument
-        // corresponds to the number of fitted vertex dimensions (i.e., 3 if we
-        // only fit spatial coordinates and 4 if we also fit time).
+        // Update the vertex with the new track. The second template
+        // argument corresponds to the number of fitted vertex dimensions
+        // (i.e., 3 if we only fit spatial coordinates and 4 if we also fit
+        // time).
         KalmanVertexUpdater::updateVertexWithTrack(*vtx, trkAtVtx,
                                                    m_cfg.useTime ? 4 : 3);
       } else {

Core/src/Vertexing/CMakeLists.txt

@@ -16,4 +16,5 @@ target_sources(
     NumericalTrackLinearizer.cpp
     TrackDensityVertexFinder.cpp
     GaussianTrackDensity.cpp
+    ImpactPointEstimator.cpp
 )