refactor!: Move and Grid Density finders to cpp (#2973)

This PR moves 

- `GaussianGridTrackDensity`
- `GridDensityVertexFinder`

to cpp files. It has to refactor the interface a bit to make this possible, mostly, the grid sizes become runtime arguments instead of compile-time parameters. This then also changes the internal grid eigen objects to be dynamic, which I think should be fine because we essentially only do lookups on them.

Part of:
- #2842

Blocked by:
- #2971

Core/include/Acts/Vertexing/GaussianGridTrackDensity.hpp

@@ -21,32 +21,43 @@ namespace Acts {
 /// The position of the highest track density (of either a single
 /// bin or the sum of a certain region) can be determined.
 /// Single tracks can be cached and removed from the overall density.
-///
-/// @tparam mainGridSize The size of the z-axis 1-dim main density grid
-/// @tparam trkGridSize The 2(!)-dim grid size of a single track, i.e.
-/// a single track is modelled as a (trkGridSize x trkGridSize) grid
-/// in the d0-z0 plane. Note: trkGridSize has to be an odd value.
-template <int mainGridSize = 2000, int trkGridSize = 15>
 class GaussianGridTrackDensity {
-  // Assert odd trkGridSize
-  static_assert(trkGridSize % 2);
-  // Assert bigger main grid than track grid
-  static_assert(mainGridSize > trkGridSize);
-
  public:
-  using MainGridVector = Eigen::Matrix<float, mainGridSize, 1>;
-  using TrackGridVector = Eigen::Matrix<float, trkGridSize, 1>;
+  using MainGridVector = Eigen::Matrix<float, Eigen::Dynamic, 1>;
+  using TrackGridVector = Eigen::Matrix<float, Eigen::Dynamic, 1>;
 
   /// The configuration struct
   struct Config {
     /// @param zMinMax_ The minimum and maximum z-values (in mm) that
     ///                 should be covered by the main 1-dim density grid along
     ///                 the z-axis
+    /// @tparam mainGridSize The size of the z-axis 1-dim main density grid
+    /// @tparam trkGridSize The 2(!)-dim grid size of a single track, i.e.
+    /// a single track is modelled as a (trkGridSize x trkGridSize) grid
+    /// in the d0-z0 plane. Note: trkGridSize has to be an odd value.
     /// @note The value of @p zMinMax_ together with @p mainGridSize determines the
     /// overall bin size to be used as seen below
-    Config(float zMinMax_ = 100) : zMinMax(zMinMax_) {
+    Config(float zMinMax_ = 100, int mainGridSize_ = 2000,
+           int trkGridSize_ = 15)
+        : mainGridSize(mainGridSize_),
+          trkGridSize(trkGridSize_),
+          zMinMax(zMinMax_) {
       binSize = 2. * zMinMax / mainGridSize;
+
+      if (trkGridSize % 2 == 0) {
+        throw std::runtime_error(
+            "GaussianGridTrackDensity: trkGridSize has to be an odd value!");
+      }
+      if (mainGridSize < trkGridSize) {
+        throw std::runtime_error(
+            "GaussianGridTrackDensity: mainGridSize has to be bigger than "
+            "trkGridSize!");
+      }
     }
+
+    int mainGridSize;
+    int trkGridSize;
+
     // Min and max z value of big grid
     float zMinMax;  // mm
 
@@ -105,6 +116,8 @@ class GaussianGridTrackDensity {
   void removeTrackGridFromMainGrid(int zBin, const TrackGridVector& trkGrid,
                                    MainGridVector& mainGrid) const;
 
+  const Config& config() const { return m_cfg; }
+
  private:
   /// @brief Helper function that actually adds the track to the
   /// main density grid
@@ -176,5 +189,3 @@ class GaussianGridTrackDensity {
 };
 
 }  // namespace Acts
-
-#include "Acts/Vertexing/GaussianGridTrackDensity.ipp"

Core/include/Acts/Vertexing/GridDensityVertexFinder.hpp

@@ -29,33 +29,18 @@ namespace Acts {
 /// density vector) needs to be calculated. All track contributions along the
 /// beam axis (main density grid) a superimposed and the z-value of the bin
 /// with the highest track density is returned as a vertex candidate.
-/// @tparam mainGridSize The size of the z-axis 1-dim main density grid
-/// @tparam trkGridSize The 2(!)-dim grid size of a single track, i.e.
-/// a single track is modelled as a (trkGridSize x trkGridSize) grid
-/// in the d0-z0 plane. Note: trkGridSize has to be an odd value.
-template <int mainGridSize = 2000, int trkGridSize = 15>
 class GridDensityVertexFinder final : public IVertexFinder {
-  // Assert odd trkGridSize
-  static_assert(trkGridSize % 2);
-  // Assert bigger main grid than track grid
-  static_assert(mainGridSize > trkGridSize);
-
-  using GridDensity = GaussianGridTrackDensity<mainGridSize, trkGridSize>;
-
  public:
-  using MainGridVector = typename GridDensity::MainGridVector;
-  using TrackGridVector = typename GridDensity::TrackGridVector;
+  using MainGridVector = GaussianGridTrackDensity::MainGridVector;
+  using TrackGridVector = GaussianGridTrackDensity::TrackGridVector;
 
   /// @brief The Config struct
   struct Config {
-    ///@param zMinMax min and max z value of big z-axis grid
-    Config(float zMinMax = 100)
-        : gridDensity(typename GridDensity::Config(zMinMax)) {}
     ///@param gDensity The grid density
-    Config(const GridDensity& gDensity) : gridDensity(gDensity) {}
+    Config(GaussianGridTrackDensity gDensity) : gridDensity(gDensity) {}
 
     // The grid density object
-    GridDensity gridDensity;
+    GaussianGridTrackDensity gridDensity;
 
     // Cache the main grid and the density contributions (trackGrid and z-bin)
     // for every single track.
@@ -82,8 +67,10 @@ class GridDensityVertexFinder final : public IVertexFinder {
   ///
   /// Only needed if cacheGridStateForTrackRemoval == true
   struct State {
+    State(MainGridVector mainGrid_) : mainGrid(std::move(mainGrid_)) {}
+
     // The main density grid
-    MainGridVector mainGrid = MainGridVector::Zero();
+    MainGridVector mainGrid;
     // Map to store z-bin and track grid (i.e. the density contribution of
     // a single track to the main grid) for every single track
     std::map<InputTrack, std::pair<int, TrackGridVector>> binAndTrackGridMap;
@@ -115,7 +102,9 @@ class GridDensityVertexFinder final : public IVertexFinder {
 
   IVertexFinder::State makeState(
       const Acts::MagneticFieldContext& /*mctx*/) const override {
-    return IVertexFinder::State{State{}};
+    return IVertexFinder::State{
+        std::in_place_type<State>,
+        MainGridVector{m_cfg.gridDensity.config().mainGridSize}};
   }
 
   void setTracksToRemove(
@@ -151,5 +140,3 @@ class GridDensityVertexFinder final : public IVertexFinder {
 };
 
 }  // namespace Acts
-
-#include "GridDensityVertexFinder.ipp"

Core/src/Vertexing/CMakeLists.txt

@@ -17,4 +17,6 @@ target_sources(
     TrackDensityVertexFinder.cpp
     GaussianTrackDensity.cpp
     ImpactPointEstimator.cpp
+    GaussianGridTrackDensity.cpp
+    GridDensityVertexFinder.cpp
 )

Core/include/Acts/Vertexing/GaussianGridTrackDensity.ipp → Core/src/Vertexing/GaussianGridTrackDensity.cpp

@@ -1,20 +1,23 @@
 // This file is part of the Acts project.
 //
-// Copyright (C) 2020 CERN for the benefit of the Acts project
+// Copyright (C) 2020-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/.
+
+#include "Acts/Vertexing/GaussianGridTrackDensity.hpp"
+
 #include "Acts/Utilities/AlgebraHelpers.hpp"
 #include "Acts/Vertexing/VertexingError.hpp"
 
 #include <algorithm>
 
-template <int mainGridSize, int trkGridSize>
-Acts::Result<float>
-Acts::GaussianGridTrackDensity<mainGridSize, trkGridSize>::getMaxZPosition(
+namespace Acts {
+
+Result<float> GaussianGridTrackDensity::getMaxZPosition(
     MainGridVector& mainGrid) const {
-  if (mainGrid == MainGridVector::Zero()) {
+  if (mainGrid.isZero()) {
     return VertexingError::EmptyInput;
   }
 
@@ -29,13 +32,12 @@ Acts::GaussianGridTrackDensity<mainGridSize, trkGridSize>::getMaxZPosition(
   }
 
   // Derive corresponding z value
-  return (zbin - mainGridSize / 2.0f + 0.5f) * m_cfg.binSize;
+  return (zbin - m_cfg.mainGridSize / 2.0f + 0.5f) * m_cfg.binSize;
 }
 
-template <int mainGridSize, int trkGridSize>
-Acts::Result<std::pair<float, float>> Acts::GaussianGridTrackDensity<
-    mainGridSize, trkGridSize>::getMaxZPositionAndWidth(MainGridVector&
-                                                            mainGrid) const {
+Result<std::pair<float, float>>
+GaussianGridTrackDensity::getMaxZPositionAndWidth(
+    MainGridVector& mainGrid) const {
   // Get z maximum value
   auto maxZRes = getMaxZPosition(mainGrid);
   if (!maxZRes.ok()) {
@@ -53,11 +55,9 @@ Acts::Result<std::pair<float, float>> Acts::GaussianGridTrackDensity<
   return returnPair;
 }
 
-template <int mainGridSize, int trkGridSize>
-std::pair<int, typename Acts::GaussianGridTrackDensity<
-                   mainGridSize, trkGridSize>::TrackGridVector>
-Acts::GaussianGridTrackDensity<mainGridSize, trkGridSize>::addTrack(
-    const Acts::BoundTrackParameters& trk, MainGridVector& mainGrid) const {
+std::pair<int, GaussianGridTrackDensity::TrackGridVector>
+GaussianGridTrackDensity::addTrack(const BoundTrackParameters& trk,
+                                   MainGridVector& mainGrid) const {
   SquareMatrix2 cov = trk.spatialImpactParameterCovariance().value();
   float d0 = trk.parameters()[0];
   float z0 = trk.parameters()[1];
@@ -66,17 +66,17 @@ Acts::GaussianGridTrackDensity<mainGridSize, trkGridSize>::addTrack(
   int dOffset = static_cast<int>(std::floor(d0 / m_cfg.binSize - 0.5) + 1);
   // Check if current track does affect grid density
   // in central bins at z-axis
-  if (std::abs(dOffset) > (trkGridSize - 1) / 2.) {
+  if (std::abs(dOffset) > (m_cfg.trkGridSize - 1) / 2.) {
     // Current track is too far away to contribute
     // to track density at z-axis bins
-    return {-1, TrackGridVector::Zero()};
+    return {-1, TrackGridVector::Zero(m_cfg.trkGridSize)};
   }
 
   // Calculate bin in z
-  int zBin = int(z0 / m_cfg.binSize + mainGridSize / 2.);
+  int zBin = int(z0 / m_cfg.binSize + m_cfg.mainGridSize / 2.);
 
-  if (zBin < 0 || zBin >= mainGridSize) {
-    return {-1, TrackGridVector::Zero()};
+  if (zBin < 0 || zBin >= m_cfg.mainGridSize) {
+    return {-1, TrackGridVector::Zero(m_cfg.trkGridSize)};
   }
   // Calculate the positions of the bin centers
   float binCtrZ = (zBin + 0.5f) * m_cfg.binSize - m_cfg.zMinMax;
@@ -93,71 +93,61 @@ Acts::GaussianGridTrackDensity<mainGridSize, trkGridSize>::addTrack(
   return {zBin, trackGrid};
 }
 
-template <int mainGridSize, int trkGridSize>
-void Acts::GaussianGridTrackDensity<mainGridSize, trkGridSize>::
-    addTrackGridToMainGrid(int zBin, const TrackGridVector& trkGrid,
-                           MainGridVector& mainGrid) const {
+void GaussianGridTrackDensity::addTrackGridToMainGrid(
+    int zBin, const TrackGridVector& trkGrid, MainGridVector& mainGrid) const {
   modifyMainGridWithTrackGrid(zBin, trkGrid, mainGrid, +1);
 }
 
-template <int mainGridSize, int trkGridSize>
-void Acts::GaussianGridTrackDensity<mainGridSize, trkGridSize>::
-    removeTrackGridFromMainGrid(int zBin, const TrackGridVector& trkGrid,
-                                MainGridVector& mainGrid) const {
+void GaussianGridTrackDensity::removeTrackGridFromMainGrid(
+    int zBin, const TrackGridVector& trkGrid, MainGridVector& mainGrid) const {
   modifyMainGridWithTrackGrid(zBin, trkGrid, mainGrid, -1);
 }
 
-template <int mainGridSize, int trkGridSize>
-void Acts::GaussianGridTrackDensity<mainGridSize, trkGridSize>::
-    modifyMainGridWithTrackGrid(int zBin, const TrackGridVector& trkGrid,
-                                MainGridVector& mainGrid,
-                                int modifyModeSign) const {
-  int width = (trkGridSize - 1) / 2;
+void GaussianGridTrackDensity::modifyMainGridWithTrackGrid(
+    int zBin, const TrackGridVector& trkGrid, MainGridVector& mainGrid,
+    int modifyModeSign) const {
+  int width = (m_cfg.trkGridSize - 1) / 2;
   // Overlap left
   int leftOL = zBin - width;
   // Overlap right
-  int rightOL = zBin + width - mainGridSize + 1;
+  int rightOL = zBin + width - m_cfg.mainGridSize + 1;
   if (leftOL < 0) {
-    int totalTrkSize = trkGridSize + leftOL;
+    int totalTrkSize = m_cfg.trkGridSize + leftOL;
     mainGrid.segment(0, totalTrkSize) +=
         modifyModeSign * trkGrid.segment(-leftOL, totalTrkSize);
     return;
   }
   if (rightOL > 0) {
-    int totalTrkSize = trkGridSize - rightOL;
-    mainGrid.segment(mainGridSize - totalTrkSize, totalTrkSize) +=
+    int totalTrkSize = m_cfg.trkGridSize - rightOL;
+    mainGrid.segment(m_cfg.mainGridSize - totalTrkSize, totalTrkSize) +=
         modifyModeSign * trkGrid.segment(0, totalTrkSize);
     return;
   }
 
-  mainGrid.segment(zBin - width, trkGridSize) += modifyModeSign * trkGrid;
+  mainGrid.segment(zBin - width, m_cfg.trkGridSize) += modifyModeSign * trkGrid;
 }
 
-template <int mainGridSize, int trkGridSize>
-typename Acts::GaussianGridTrackDensity<mainGridSize,
-                                        trkGridSize>::TrackGridVector
-Acts::GaussianGridTrackDensity<mainGridSize, trkGridSize>::createTrackGrid(
-    float d0, float distCtrZ, const Acts::SquareMatrix2& cov) const {
-  TrackGridVector trackGrid(TrackGridVector::Zero());
-  float floorHalfTrkGridSize = static_cast<float>(trkGridSize) / 2 - 0.5f;
+GaussianGridTrackDensity::TrackGridVector
+GaussianGridTrackDensity::createTrackGrid(float d0, float distCtrZ,
+                                          const SquareMatrix2& cov) const {
+  TrackGridVector trackGrid(TrackGridVector::Zero(m_cfg.trkGridSize));
+  float floorHalfTrkGridSize = static_cast<float>(m_cfg.trkGridSize) / 2 - 0.5f;
 
   // Loop over columns
-  for (int j = 0; j < trkGridSize; j++) {
+  for (int j = 0; j < m_cfg.trkGridSize; j++) {
     float z = (j - floorHalfTrkGridSize) * m_cfg.binSize;
     trackGrid(j) = normal2D(-d0, z - distCtrZ, cov);
   }
   return trackGrid;
 }
 
-template <int mainGridSize, int trkGridSize>
-Acts::Result<float>
-Acts::GaussianGridTrackDensity<mainGridSize, trkGridSize>::estimateSeedWidth(
+Result<float> GaussianGridTrackDensity::estimateSeedWidth(
     MainGridVector& mainGrid, float maxZ) const {
-  if (mainGrid == MainGridVector::Zero()) {
+  if (mainGrid.isZero()) {
     return VertexingError::EmptyInput;
   }
   // Get z bin of max density z value
-  int zBin = int(maxZ / m_cfg.binSize + mainGridSize / 2.);
+  int zBin = int(maxZ / m_cfg.binSize + m_cfg.mainGridSize / 2.);
 
   const float maxValue = mainGrid(zBin);
   float gridValue = mainGrid(zBin);
@@ -195,9 +185,8 @@ Acts::GaussianGridTrackDensity<mainGridSize, trkGridSize>::estimateSeedWidth(
   return std::isnormal(width) ? width : 0.0f;
 }
 
-template <int mainGridSize, int trkGridSize>
-float Acts::GaussianGridTrackDensity<mainGridSize, trkGridSize>::normal2D(
-    float d, float z, const Acts::SquareMatrix2& cov) const {
+float GaussianGridTrackDensity::normal2D(float d, float z,
+                                         const SquareMatrix2& cov) const {
   float det = cov.determinant();
   float coef = 1 / std::sqrt(det);
   float expo =
@@ -206,9 +195,8 @@ float Acts::GaussianGridTrackDensity<mainGridSize, trkGridSize>::normal2D(
   return coef * safeExp(expo);
 }
 
-template <int mainGridSize, int trkGridSize>
-int Acts::GaussianGridTrackDensity<mainGridSize, trkGridSize>::
-    getHighestSumZPosition(MainGridVector& mainGrid) const {
+int GaussianGridTrackDensity::getHighestSumZPosition(
+    MainGridVector& mainGrid) const {
   // Checks the first (up to) 3 density maxima, if they are close, checks which
   // one has the highest surrounding density sum (the two neighboring bins)
 
@@ -255,9 +243,8 @@ int Acts::GaussianGridTrackDensity<mainGridSize, trkGridSize>::
   return zFirstMax;
 }
 
-template <int mainGridSize, int trkGridSize>
-double Acts::GaussianGridTrackDensity<mainGridSize, trkGridSize>::getDensitySum(
-    const MainGridVector& mainGrid, int pos) const {
+double GaussianGridTrackDensity::getDensitySum(const MainGridVector& mainGrid,
+                                               int pos) const {
   double sum = mainGrid(pos);
   // Sum up only the density contributions from the
   // neighboring bins if they are still within bounds
@@ -269,3 +256,5 @@ double Acts::GaussianGridTrackDensity<mainGridSize, trkGridSize>::getDensitySum(
   }
   return sum;
 }
+
+}  // namespace Acts