refactor: Reduce number of surface allocations in vertexing (#2443)

Core/include/Acts/Vertexing/AdaptiveMultiVertexFitter.ipp

@@ -257,6 +257,11 @@ Acts::Result<void> Acts::
         const VertexingOptions<input_track_t>& vertexingOptions) const {
   for (auto vtx : state.vertexCollection) {
     VertexInfo<input_track_t>& currentVtxInfo = state.vtxInfoMap[vtx];
+
+    const std::shared_ptr<PerigeeSurface> vtxPerigeeSurface =
+        Surface::makeShared<PerigeeSurface>(
+            VectorHelpers::position(state.vtxInfoMap[vtx].oldPosition));
+
     for (const auto& trk : currentVtxInfo.trackLinks) {
       auto& trkAtVtx = state.tracksAtVerticesMap.at(std::make_pair(trk, vtx));
 
@@ -270,9 +275,9 @@ Acts::Result<void> Acts::
         // Check if linearization state exists or need to be relinearized
         if (not trkAtVtx.isLinearized || state.vtxInfoMap[vtx].relinearize) {
           auto result = linearizer.linearizeTrack(
-              m_extractParameters(*trk), state.vtxInfoMap[vtx].oldPosition,
-              vertexingOptions.geoContext, vertexingOptions.magFieldContext,
-              state.linearizerState);
+              m_extractParameters(*trk), state.vtxInfoMap[vtx].oldPosition[3],
+              *vtxPerigeeSurface, vertexingOptions.geoContext,
+              vertexingOptions.magFieldContext, state.linearizerState);
           if (!result.ok()) {
             return result.error();
           }

Core/include/Acts/Vertexing/FullBilloirVertexFitter.ipp

@@ -101,15 +101,22 @@ Acts::FullBilloirVertexFitter<input_track_t, linearizer_t>::fit(
     double newChi2 = 0;
     BilloirVertex billoirVertex;
 
+    Vector3 linPointPos = VectorHelpers::position(linPoint);
+    // Make Perigee surface at linPointPos, transverse plane of Perigee
+    // corresponds the global x-y plane
+    const std::shared_ptr<PerigeeSurface> perigeeSurface =
+        Surface::makeShared<PerigeeSurface>(linPointPos);
+
     // iterate over all tracks
     for (std::size_t iTrack = 0; iTrack < nTracks; ++iTrack) {
       const input_track_t* trackContainer = paramVector[iTrack];
 
       const auto& trackParams = extractParameters(*trackContainer);
 
       auto result = linearizer.linearizeTrack(
-          trackParams, linPoint, vertexingOptions.geoContext,
-          vertexingOptions.magFieldContext, state.linearizerState);
+          trackParams, linPoint[3], *perigeeSurface,
+          vertexingOptions.geoContext, vertexingOptions.magFieldContext,
+          state.linearizerState);
       if (!result.ok()) {
         return result.error();
       }

Core/include/Acts/Vertexing/HelicalTrackLinearizer.hpp

@@ -100,16 +100,18 @@ class HelicalTrackLinearizer {
   /// the PCA to a given Perigee surface
   ///
   /// @param params Parameters to linearize
-  /// @param linPoint Point which defines the Perigee.
+  /// @param linPointTime Time associated to the linearization point
   /// @note Transverse plane of the Perigee corresponding to @p linPoint is
   /// parallel to the global x-y plane
+  /// @param perigeeSurface Perigee surface belonging to @p linPoint
   /// @param gctx Geometry context
   /// @param mctx Magnetic field context
   /// @param state Linearizer state object
   ///
   /// @return Linearized track
   Result<LinearizedTrack> linearizeTrack(const BoundTrackParameters& params,
-                                         const Vector4& linPoint,
+                                         double linPointTime,
+                                         const Surface& perigeeSurface,
                                          const Acts::GeometryContext& gctx,
                                          const Acts::MagneticFieldContext& mctx,
                                          State& state) const;

Core/include/Acts/Vertexing/HelicalTrackLinearizer.ipp

@@ -12,15 +12,9 @@
 template <typename propagator_t, typename propagator_options_t>
 Acts::Result<Acts::LinearizedTrack> Acts::
     HelicalTrackLinearizer<propagator_t, propagator_options_t>::linearizeTrack(
-        const BoundTrackParameters& params, const Vector4& linPoint,
-        const Acts::GeometryContext& gctx,
+        const BoundTrackParameters& params, double linPointTime,
+        const Surface& perigeeSurface, const Acts::GeometryContext& gctx,
         const Acts::MagneticFieldContext& mctx, State& state) const {
-  // Make Perigee surface at linPointPos, transverse plane of Perigee
-  // corresponds the global x-y plane
-  Vector3 linPointPos = VectorHelpers::position(linPoint);
-  const std::shared_ptr<PerigeeSurface> perigeeSurface =
-      Surface::makeShared<PerigeeSurface>(linPointPos);
-
   // Create propagator options
   propagator_options_t pOptions(gctx, mctx);
 
@@ -34,7 +28,7 @@ Acts::Result<Acts::LinearizedTrack> Acts::
   // forward or backward to arrive at the PCA.
   auto intersection =
       perigeeSurface
-          ->intersect(gctx, params.position(gctx), params.direction(), false)
+          .intersect(gctx, params.position(gctx), params.direction(), false)
           .closest();
 
   // Setting the propagation direction using the intersection length from
@@ -44,14 +38,14 @@ Acts::Result<Acts::LinearizedTrack> Acts::
   pOptions.direction =
       Direction::fromScalarZeroAsPositive(intersection.pathLength());
 
-  // Propagate to the PCA of linPointPos
-  auto result = m_cfg.propagator->propagate(params, *perigeeSurface, pOptions);
+  // Propagate to the PCA of the reference point
+  auto result = m_cfg.propagator->propagate(params, perigeeSurface, pOptions);
   if (not result.ok()) {
     return result.error();
   }
 
   // Extracting the track parameters at said PCA - this corresponds to the
-  // Perigee representation of the track wrt linPointPos
+  // Perigee representation of the track wrt the reference point
   const auto& endParams = *result->endParameters;
   BoundVector paramsAtPCA = endParams.parameters();
 
@@ -143,10 +137,10 @@ Acts::Result<Acts::LinearizedTrack> Acts::
     ActsScalar h = (rho < 0.) ? -1 : 1;
 
     // Quantities from Eq. 5.34 in Ref. (1) (see .hpp)
-    ActsScalar X = pca(0) - linPointPos.x() + rho * sinPhi;
-    ActsScalar Y = pca(1) - linPointPos.y() - rho * cosPhi;
+    ActsScalar X = pca(0) - perigeeSurface.center(gctx).x() + rho * sinPhi;
+    ActsScalar Y = pca(1) - perigeeSurface.center(gctx).y() - rho * cosPhi;
     ActsScalar S2 = (X * X + Y * Y);
-    // S is the 2D distance from the helix center to linPointPos
+    // S is the 2D distance from the helix center to the reference point
     // in the x-y plane
     ActsScalar S = std::sqrt(S2);
 
@@ -204,6 +198,10 @@ Acts::Result<Acts::LinearizedTrack> Acts::
   // The parameter weight
   BoundSquareMatrix weightAtPCA = parCovarianceAtPCA.inverse();
 
+  Vector4 linPoint;
+  linPoint.head<3>() = perigeeSurface.center(gctx);
+  linPoint[3] = linPointTime;
+
   return LinearizedTrack(paramsAtPCA, parCovarianceAtPCA, weightAtPCA, linPoint,
                          positionJacobian, momentumJacobian, pca, momentumAtPCA,
                          constTerm);

Core/include/Acts/Vertexing/IterativeVertexFinder.hpp

@@ -200,10 +200,12 @@ class IterativeVertexFinder {
   ///
   /// @param params Track parameters
   /// @param vertex The vertex
+  /// @param perigeeSurface The perigee surface at vertex position
   /// @param vertexingOptions Vertexing options
   /// @param state The state object
   Result<double> getCompatibility(
       const BoundTrackParameters& params, const Vertex<InputTrack_t>& vertex,
+      const Surface& perigeeSurface,
       const VertexingOptions<InputTrack_t>& vertexingOptions,
       State& state) const;
 

Core/include/Acts/Vertexing/IterativeVertexFinder.ipp

@@ -220,12 +220,14 @@ template <typename vfitter_t, typename sfinder_t>
 Acts::Result<double>
 Acts::IterativeVertexFinder<vfitter_t, sfinder_t>::getCompatibility(
     const BoundTrackParameters& params, const Vertex<InputTrack_t>& vertex,
+    const Surface& perigeeSurface,
     const VertexingOptions<InputTrack_t>& vertexingOptions,
     State& state) const {
   // Linearize track
   auto result = m_cfg.linearizer.linearizeTrack(
-      params, vertex.fullPosition(), vertexingOptions.geoContext,
-      vertexingOptions.magFieldContext, state.linearizerState);
+      params, vertex.fullPosition()[3], perigeeSurface,
+      vertexingOptions.geoContext, vertexingOptions.magFieldContext,
+      state.linearizerState);
   if (!result.ok()) {
     return result.error();
   }
@@ -301,10 +303,15 @@ Acts::IterativeVertexFinder<vfitter_t, sfinder_t>::removeUsedCompatibleTracks(
   // m_cfg.cutOffTrackWeight threshold and are hence outliers
   ACTS_DEBUG("Number of outliers: " << perigeesToFit.size());
 
+  const std::shared_ptr<PerigeeSurface> myVertexPerigeeSurface =
+      Surface::makeShared<PerigeeSurface>(
+          VectorHelpers::position(myVertex.fullPosition()));
+
   for (const auto& myPerigeeToFit : perigeesToFit) {
     // calculate chi2 w.r.t. last fitted vertex
-    auto result = getCompatibility(m_extractParameters(*myPerigeeToFit),
-                                   myVertex, vertexingOptions, state);
+    auto result =
+        getCompatibility(m_extractParameters(*myPerigeeToFit), myVertex,
+                         *myVertexPerigeeSurface, vertexingOptions, state);
 
     if (!result.ok()) {
       return result.error();
@@ -428,6 +435,10 @@ Acts::IterativeVertexFinder<vfitter_t, sfinder_t>::reassignTracksToNewVertex(
     State& state) const {
   int numberOfAddedTracks = 0;
 
+  const std::shared_ptr<PerigeeSurface> currentVertexPerigeeSurface =
+      Surface::makeShared<PerigeeSurface>(
+          VectorHelpers::position(currentVertex.fullPosition()));
+
   // iterate over all vertices and check if tracks need to be reassigned
   // to new (current) vertex
   for (auto& vertexIt : vertexCollection) {
@@ -436,6 +447,10 @@ Acts::IterativeVertexFinder<vfitter_t, sfinder_t>::reassignTracksToNewVertex(
     auto tracksBegin = tracksAtVertex.begin();
     auto tracksEnd = tracksAtVertex.end();
 
+    const std::shared_ptr<PerigeeSurface> vertexItPerigeeSurface =
+        Surface::makeShared<PerigeeSurface>(
+            VectorHelpers::position(vertexIt.fullPosition()));
+
     for (auto tracksIter = tracksBegin; tracksIter != tracksEnd;) {
       // consider only tracks that are not too tightly assigned to other
       // vertex
@@ -449,14 +464,16 @@ Acts::IterativeVertexFinder<vfitter_t, sfinder_t>::reassignTracksToNewVertex(
 
       // compute compatibility
       auto resultNew = getCompatibility(trackPerigee, currentVertex,
+                                        *currentVertexPerigeeSurface,
                                         vertexingOptions, state);
       if (!resultNew.ok()) {
         return Result<bool>::failure(resultNew.error());
       }
       double chi2NewVtx = *resultNew;
 
       auto resultOld =
-          getCompatibility(trackPerigee, vertexIt, vertexingOptions, state);
+          getCompatibility(trackPerigee, vertexIt, *vertexItPerigeeSurface,
+                           vertexingOptions, state);
       if (!resultOld.ok()) {
         return Result<bool>::failure(resultOld.error());
       }

Core/include/Acts/Vertexing/LinearizerConcept.hpp

@@ -35,11 +35,12 @@ METHOD_TRAIT(linTrack_t, linearizeTrack);
 
          constexpr static bool linTrack_exists = has_method<const S, Result<LinearizedTrack>,
          linTrack_t, const BoundTrackParameters&,
-                     const Vector4&,
+                     double,
+                     const Surface&,
                      const Acts::GeometryContext&,
                      const Acts::MagneticFieldContext&,
                      typename S::State&>;
-  
+
         static_assert(linTrack_exists, "linearizeTrack method not found");
 
         constexpr static bool propagator_exists = exists<propagator_t, S>;

Core/include/Acts/Vertexing/NumericalTrackLinearizer.hpp

@@ -116,15 +116,17 @@ class NumericalTrackLinearizer {
   /// the PCA to a given Perigee surface
   ///
   /// @param params Parameters to linearize
-  /// @param linPoint Point which defines the Perigee.
+  /// @param linPointTime Time associated to the linearization point
   /// @note Transverse plane of the Perigee corresponding to @p linPoint is
   /// parallel to the global x-y plane
+  /// @param perigeeSurface Perigee surface belonging to @p linPoint
   /// @param gctx Geometry context
   /// @param mctx Magnetic field context
   ///
   /// @return Linearized track
   Result<LinearizedTrack> linearizeTrack(const BoundTrackParameters& params,
-                                         const Vector4& linPoint,
+                                         double linPointTime,
+                                         const Surface& perigeeSurface,
                                          const Acts::GeometryContext& gctx,
                                          const Acts::MagneticFieldContext& mctx,
                                          State& /*state*/) const;

Core/include/Acts/Vertexing/NumericalTrackLinearizer.ipp

@@ -13,16 +13,11 @@
 template <typename propagator_t, typename propagator_options_t>
 Acts::Result<Acts::LinearizedTrack>
 Acts::NumericalTrackLinearizer<propagator_t, propagator_options_t>::
-    linearizeTrack(const BoundTrackParameters& params, const Vector4& linPoint,
+    linearizeTrack(const BoundTrackParameters& params, double linPointTime,
+                   const Surface& perigeeSurface,
                    const Acts::GeometryContext& gctx,
                    const Acts::MagneticFieldContext& mctx,
                    State& /*state*/) const {
-  // Make Perigee surface at linPointPos, transverse plane of Perigee
-  // corresponds the global x-y plane
-  Vector3 linPointPos = VectorHelpers::position(linPoint);
-  std::shared_ptr<PerigeeSurface> perigeeSurface =
-      Surface::makeShared<PerigeeSurface>(linPointPos);
-
   // Create propagator options
   propagator_options_t pOptions(gctx, mctx);
 
@@ -36,7 +31,7 @@ Acts::NumericalTrackLinearizer<propagator_t, propagator_options_t>::
   // forward or backward to arrive at the PCA.
   auto intersection =
       perigeeSurface
-          ->intersect(gctx, params.position(gctx), params.direction(), false)
+          .intersect(gctx, params.position(gctx), params.direction(), false)
           .closest();
 
   // Setting the propagation direction using the intersection length from
@@ -46,17 +41,17 @@ Acts::NumericalTrackLinearizer<propagator_t, propagator_options_t>::
   pOptions.direction =
       Direction::fromScalarZeroAsPositive(intersection.pathLength());
 
-  // Propagate to the PCA of linPointPos
-  auto result = m_cfg.propagator->propagate(params, *perigeeSurface, pOptions);
+  // Propagate to the PCA of the reference point
+  auto result = m_cfg.propagator->propagate(params, perigeeSurface, pOptions);
   if (not result.ok()) {
     return result.error();
   }
 
-  // Extracting the Perigee representation of the track wrt linPointPos
+  // Extracting the Perigee representation of the track wrt the reference point
   auto endParams = *result->endParameters;
   BoundVector perigeeParams = endParams.parameters();
 
-  // Covariance and weight matrix at the PCA to "linPoint"
+  // Covariance and weight matrix at the PCA to the reference point
   BoundSquareMatrix parCovarianceAtPCA = endParams.covariance().value();
   BoundSquareMatrix weightAtPCA = parCovarianceAtPCA.inverse();
 
@@ -92,7 +87,7 @@ Acts::NumericalTrackLinearizer<propagator_t, propagator_options_t>::
   ActsMatrix<eBoundSize, eLinSize> completeJacobian =
       ActsMatrix<eBoundSize, eLinSize>::Zero(eBoundSize, eLinSize);
 
-  // Perigee parameters wrt linPoint after wiggling
+  // Perigee parameters wrt the reference point after wiggling
   BoundVector newPerigeeParams;
 
   // Check if wiggled angle theta are within definition range [0, pi]
@@ -121,16 +116,16 @@ Acts::NumericalTrackLinearizer<propagator_t, propagator_options_t>::
         paramVecCopy(eLinQOverP), std::nullopt, ParticleHypothesis::pion());
 
     // Obtain propagation direction
-    intersection = perigeeSurface
-                       ->intersect(gctx, paramVecCopy.template head<3>(),
-                                   wiggledDir, false)
-                       .closest();
+    intersection =
+        perigeeSurface
+            .intersect(gctx, paramVecCopy.template head<3>(), wiggledDir, false)
+            .closest();
     pOptions.direction =
         Direction::fromScalarZeroAsPositive(intersection.pathLength());
 
     // Propagate to the new PCA and extract Perigee parameters
     auto newResult = m_cfg.propagator->propagate(wiggledCurvilinearParams,
-                                                 *perigeeSurface, pOptions);
+                                                 perigeeSurface, pOptions);
     if (not newResult.ok()) {
       return newResult.error();
     }
@@ -157,6 +152,10 @@ Acts::NumericalTrackLinearizer<propagator_t, propagator_options_t>::
   BoundVector constTerm =
       perigeeParams - positionJacobian * pca - momentumJacobian * momentumAtPCA;
 
+  Vector4 linPoint;
+  linPoint.head<3>() = perigeeSurface.center(gctx);
+  linPoint[3] = linPointTime;
+
   return LinearizedTrack(perigeeParams, parCovarianceAtPCA, weightAtPCA,
                          linPoint, positionJacobian, momentumJacobian, pca,
                          momentumAtPCA, constTerm);

Tests/UnitTests/Core/Vertexing/KalmanVertexTrackUpdaterTests.cpp

@@ -153,11 +153,14 @@ BOOST_AUTO_TEST_CASE(Kalman_Vertex_TrackUpdater) {
     BoundTrackParameters params(perigeeSurface, paramVec, std::move(covMat),
                                 ParticleHypothesis::pion());
 
+    std::shared_ptr<PerigeeSurface> perigee =
+        Surface::makeShared<PerigeeSurface>(Vector3::Zero());
+
     // Linearized state of the track
     LinearizedTrack linTrack =
         linearizer
-            .linearizeTrack(params, Vector4::Zero(), geoContext,
-                            magFieldContext, linState)
+            .linearizeTrack(params, 0, *perigee, geoContext, magFieldContext,
+                            linState)
             .value();
 
     // Create TrackAtVertex

Tests/UnitTests/Core/Vertexing/KalmanVertexUpdaterTests.cpp

@@ -149,11 +149,14 @@ BOOST_AUTO_TEST_CASE(Kalman_Vertex_Updater) {
     BoundTrackParameters params(perigeeSurface, paramVec, std::move(covMat),
                                 ParticleHypothesis::pion());
 
+    std::shared_ptr<PerigeeSurface> perigee =
+        Surface::makeShared<PerigeeSurface>(Vector3::Zero());
+
     // Linearized state of the track
     LinearizedTrack linTrack =
         linearizer
-            .linearizeTrack(params, Vector4::Zero(), geoContext,
-                            magFieldContext, state)
+            .linearizeTrack(params, 0, *perigee, geoContext, magFieldContext,
+                            state)
             .value();
 
     // Create TrackAtVertex