fix: minor bugs in grid vertex seeder (#2376)

Some small improvements in `AdaptiveGridTrackDensity` and `GaussianGridTrackDensity`:

- `1/(2 * M_PI)` is removed from the track density since constant factors have no influence on the position of the maximum z value.
- The central bin was twice as large as the others in `AdaptiveGridTrackDensity`. I repositioned the grid such that the central bin corresponds to the interval `[-binSize/2, binSize/2)`. As a consequence, I had to modify the unit tests.
- There was a sign mistake in the computation of the exponent. I removed it and check if the new sign is correct in the unit test `compare_to_analytical_solution_for_single_track`.
- There was a sign mistake and a typo in the computation of the seed width. I removed them and check that the new computation is correct in the unit tests `check_seed_width_estimation` and `compare_to_analytical_solution_for_single_track`.

Core/include/Acts/Vertexing/AdaptiveGridTrackDensity.hpp

@@ -1,6 +1,6 @@
 // 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
@@ -62,6 +62,16 @@ class AdaptiveGridTrackDensity {
 
   AdaptiveGridTrackDensity(const Config& cfg) : m_cfg(cfg) {}
 
+  /// @brief Calculates the bin center from the bin number
+  /// @param bin Bin number
+  /// @return Bin center
+  float getBinCenter(int bin) const;
+
+  /// @brief Calculates the bin number corresponding to a d or z value
+  /// @param value d or z value
+  /// @return Bin number
+  int getBin(float value) const;
+
   /// @brief Finds the maximum density of a DensityMap
   /// @param densityMap Map between z bins and corresponding density value
   /// @return Iterator of the map entry with the highest density value
@@ -95,28 +105,29 @@ class AdaptiveGridTrackDensity {
 
   /// @brief Removes a track from the overall grid density
   ///
-  /// @param trackDensityMap Map from z bins to corresponding track density. The track density comes from a single track.
-  /// @param mainDensityMap Map from z bins to corresponding track density. The track density comes an arbitrary number of tracks.
+  /// @param trackDensityMap Map from z bins to corresponding track density.
+  /// @note The track density comes from a single track.
+  /// @param mainDensityMap Map from z bins to corresponding track density.
+  /// @note The track density comes an arbitrary number of tracks.
   void subtractTrack(const DensityMap& trackDensityMap,
                      DensityMap& mainDensityMap) const;
 
  private:
-  /// @brief Function that creates a track density map, i.e., a map of z bins to corresponding density values coming from a single track.
+  /// @brief Function that creates a track density map, i.e., a map of z bins
+  /// to corresponding density values coming from a single track.
   ///
-  /// @param offset Offset in d0 direction, to account for the 2-dim part
-  /// of the Gaussian track distribution
-  /// @param cov The track covariance matrix
-  /// @param distCtrD The distance in d0 from the track position to its
-  /// bin center in the 2-dim grid
+  /// @param d0 Transverse impact parameter
+  /// @param z0 Longitudinal impact parameter
   /// @param centralZBin Central z bin of the track (where its density is the highest)
-  /// @param distCtrZ The distance in z0 from the track position to its
-  /// bin center in the 2-dim grid
-  DensityMap createTrackGrid(int offset, const SquareMatrix2& cov,
-                             float distCtrD, int centralZBin,
-                             float distCtrZ) const;
+  /// @param cov 2x2 impact parameter covariance matrix
+  DensityMap createTrackGrid(float d0, float z0, int centralZBin,
+                             const Acts::SquareMatrix2& cov) const;
 
   /// @brief Function that estimates the seed width based on the full width
   /// at half maximum (FWHM) of the maximum density peak
+  /// @note This only works if the maximum is sufficiently isolated since
+  /// overlapping neighboring peaks might lead to an overestimation of the
+  /// seed width.
   ///
   /// @param densityMap Map from z bins to corresponding track density
   /// @param maxZ z-position of the maximum density value
@@ -126,6 +137,7 @@ class AdaptiveGridTrackDensity {
                                   float maxZ) const;
 
   /// @brief Helper to retrieve values according to a 2-dim normal distribution
+  /// @note This function is defined in coordinate system centered around d0 and z0
   float normal2D(float d, float z, const SquareMatrix2& cov) const;
 
   /// @brief Checks the (up to) first three density maxima (only those that have

Core/include/Acts/Vertexing/AdaptiveGridTrackDensity.ipp

@@ -1,6 +1,6 @@
 // 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
@@ -10,6 +10,16 @@
 
 #include <algorithm>
 
+template <int trkGridSize>
+float Acts::AdaptiveGridTrackDensity<trkGridSize>::getBinCenter(int bin) const {
+  return bin * m_cfg.binSize;
+}
+
+template <int trkGridSize>
+int Acts::AdaptiveGridTrackDensity<trkGridSize>::getBin(float value) const {
+  return static_cast<int>(std::floor(value / m_cfg.binSize - 0.5) + 1);
+}
+
 template <int trkGridSize>
 typename Acts::AdaptiveGridTrackDensity<trkGridSize>::DensityMap::const_iterator
 Acts::AdaptiveGridTrackDensity<trkGridSize>::highestDensityEntry(
@@ -21,6 +31,7 @@ Acts::AdaptiveGridTrackDensity<trkGridSize>::highestDensityEntry(
       });
   return maxEntry;
 }
+
 template <int trkGridSize>
 Acts::Result<float>
 Acts::AdaptiveGridTrackDensity<trkGridSize>::getMaxZPosition(
@@ -39,8 +50,7 @@ Acts::AdaptiveGridTrackDensity<trkGridSize>::getMaxZPosition(
   }
 
   // Derive corresponding z value
-  int sign = (zBin > 0) ? +1 : -1;
-  return (zBin + sign * 0.5f) * m_cfg.binSize;
+  return getBinCenter(zBin);
 }
 
 template <int trkGridSize>
@@ -72,30 +82,17 @@ Acts::AdaptiveGridTrackDensity<trkGridSize>::addTrack(
   float d0 = trk.parameters()[0];
   float z0 = trk.parameters()[1];
 
-  // Calculate offset in d direction to central bin at z-axis
-  int dOffset = static_cast<int>(std::floor(d0 / m_cfg.binSize - 0.5) + 1);
+  // Calculate bin in d direction
+  int centralDBin = getBin(d0);
   // Check if current track affects grid density
-  // in central bins at z-axis
-  if (std::abs(dOffset) > (trkGridSize - 1) / 2.) {
+  if (std::abs(centralDBin) > (trkGridSize - 1) / 2.) {
     DensityMap emptyTrackDensityMap;
     return emptyTrackDensityMap;
   }
-  // Calculate bin in z
-  int centralZBin = int(z0 / m_cfg.binSize);
-
-  // Calculate the positions of the bin centers
-  float binCtrD = dOffset * m_cfg.binSize;
+  // Calculate bin in z direction
+  int centralZBin = getBin(z0);
 
-  int sign = (z0 > 0) ? +1 : -1;
-  float binCtrZ = (centralZBin + sign * 0.5f) * m_cfg.binSize;
-
-  // Calculate the distance between IP values and their
-  // corresponding bin centers
-  float distCtrD = d0 - binCtrD;
-  float distCtrZ = z0 - binCtrZ;
-
-  DensityMap trackDensityMap =
-      createTrackGrid(dOffset, cov, distCtrD, centralZBin, distCtrZ);
+  DensityMap trackDensityMap = createTrackGrid(d0, z0, centralZBin, cov);
 
   for (const auto& densityEntry : trackDensityMap) {
     int zBin = densityEntry.first;
@@ -122,19 +119,16 @@ void Acts::AdaptiveGridTrackDensity<trkGridSize>::subtractTrack(
 template <int trkGridSize>
 typename Acts::AdaptiveGridTrackDensity<trkGridSize>::DensityMap
 Acts::AdaptiveGridTrackDensity<trkGridSize>::createTrackGrid(
-    int offset, const Acts::SquareMatrix2& cov, float distCtrD, int centralZBin,
-    float distCtrZ) const {
+    float d0, float z0, int centralZBin, const Acts::SquareMatrix2& cov) const {
   DensityMap trackDensityMap;
 
   int halfTrkGridSize = (trkGridSize - 1) / 2;
-  float i = halfTrkGridSize + offset;
-  float d = (i - static_cast<float>(trkGridSize) / 2 + 0.5f) * m_cfg.binSize;
-
   int firstZBin = centralZBin - halfTrkGridSize;
   // Loop over columns
   for (int j = 0; j < trkGridSize; j++) {
-    float z = (j - static_cast<float>(trkGridSize) / 2 + 0.5f) * m_cfg.binSize;
-    trackDensityMap[firstZBin + j] = normal2D(d + distCtrD, z + distCtrZ, cov);
+    int zBin = firstZBin + j;
+    float z = getBinCenter(zBin);
+    trackDensityMap[zBin] = normal2D(-d0, z - z0, cov);
   }
   return trackDensityMap;
 }
@@ -148,40 +142,56 @@ Acts::AdaptiveGridTrackDensity<trkGridSize>::estimateSeedWidth(
   }
 
   // Get z bin of max density and max density value
-  int sign = (maxZ > 0) ? +1 : -1;
-  int zMaxBin = int(maxZ / m_cfg.binSize - sign * 0.5f);
+  int zMaxBin = getBin(maxZ);
   const float maxValue = densityMap.at(zMaxBin);
 
   int rhmBin = zMaxBin;
   float gridValue = maxValue;
+  // Boolean indicating whether we find a filled bin that has a densityValue <=
+  // maxValue/2
+  bool binFilled = true;
   while (gridValue > maxValue / 2) {
     // Check if we are still operating on continuous z values
     if (densityMap.count(rhmBin + 1) == 0) {
+      binFilled = false;
       break;
     }
     rhmBin += 1;
     gridValue = densityMap.at(rhmBin);
   }
 
   // Use linear approximation to find better z value for FWHM between bins
-  float deltaZ1 =
-      (maxValue / 2 - densityMap.at(rhmBin - 1)) *
-      (m_cfg.binSize / (densityMap.at(rhmBin - 1) - densityMap.at(rhmBin)));
+  float rightDensity = 0;
+  if (binFilled) {
+    rightDensity = densityMap.at(rhmBin);
+  }
+  float leftDensity = densityMap.at(rhmBin - 1);
+  float deltaZ1 = m_cfg.binSize * (maxValue / 2 - leftDensity) /
+                  (rightDensity - leftDensity);
+
   int lhmBin = zMaxBin;
   gridValue = maxValue;
+  binFilled = true;
   while (gridValue > maxValue / 2) {
     // Check if we are still operating on continuous z values
     if (densityMap.count(lhmBin - 1) == 0) {
+      binFilled = false;
       break;
     }
     lhmBin -= 1;
     gridValue = densityMap.at(lhmBin);
   }
 
   // Use linear approximation to find better z value for FWHM between bins
-  float deltaZ2 =
-      (maxValue / 2 - densityMap.at(lhmBin + 1)) *
-      (m_cfg.binSize / (densityMap.at(rhmBin + 1) - densityMap.at(rhmBin)));
+  rightDensity = densityMap.at(lhmBin + 1);
+  if (binFilled) {
+    leftDensity = densityMap.at(lhmBin);
+  } else {
+    leftDensity = 0;
+  }
+  float deltaZ2 = m_cfg.binSize * (rightDensity - maxValue / 2) /
+                  (rightDensity - leftDensity);
+
   float fwhm = (rhmBin - lhmBin) * m_cfg.binSize - deltaZ1 - deltaZ2;
 
   // FWHM = 2.355 * sigma
@@ -194,10 +204,10 @@ template <int trkGridSize>
 float Acts::AdaptiveGridTrackDensity<trkGridSize>::normal2D(
     float d, float z, const Acts::SquareMatrix2& cov) const {
   float det = cov.determinant();
-  float coef = 1 / (2 * M_PI * std::sqrt(det));
+  float coef = 1 / std::sqrt(det);
   float expo =
       -1 / (2 * det) *
-      (cov(1, 1) * d * d - d * z * (cov(0, 1) + cov(1, 0)) + cov(0, 0) * z * z);
+      (cov(1, 1) * d * d - (cov(0, 1) + cov(1, 0)) * d * z + cov(0, 0) * z * z);
   return coef * safeExp(expo);
 }
 
@@ -209,43 +219,43 @@ int Acts::AdaptiveGridTrackDensity<trkGridSize>::highestDensitySumBin(
 
   // The global maximum
   auto firstMax = highestDensityEntry(densityMap);
-  int zFirstMax = firstMax->first;
+  int zBinFirstMax = firstMax->first;
   float valueFirstMax = firstMax->second;
-  float firstSum = getDensitySum(densityMap, zFirstMax);
+  float firstSum = getDensitySum(densityMap, zBinFirstMax);
   float densityDeviation = valueFirstMax * m_cfg.maxRelativeDensityDev;
 
   // Get the second highest maximum
-  densityMap.at(zFirstMax) = 0;
+  densityMap.at(zBinFirstMax) = 0;
   auto secondMax = highestDensityEntry(densityMap);
-  int zSecondMax = secondMax->first;
+  int zBinSecondMax = secondMax->first;
   float valueSecondMax = secondMax->second;
   float secondSum = 0;
   if (valueFirstMax - valueSecondMax < densityDeviation) {
-    secondSum = getDensitySum(densityMap, zSecondMax);
+    secondSum = getDensitySum(densityMap, zBinSecondMax);
   }
 
   // Get the third highest maximum
-  densityMap.at(zSecondMax) = 0;
+  densityMap.at(zBinSecondMax) = 0;
   auto thirdMax = highestDensityEntry(densityMap);
-  int zThirdMax = thirdMax->first;
+  int zBinThirdMax = thirdMax->first;
   float valueThirdMax = thirdMax->second;
   float thirdSum = 0;
   if (valueFirstMax - valueThirdMax < densityDeviation) {
-    thirdSum = getDensitySum(densityMap, zThirdMax);
+    thirdSum = getDensitySum(densityMap, zBinThirdMax);
   }
 
   // Revert back to original values
-  densityMap.at(zFirstMax) = valueFirstMax;
-  densityMap.at(zSecondMax) = valueSecondMax;
+  densityMap.at(zBinFirstMax) = valueFirstMax;
+  densityMap.at(zBinSecondMax) = valueSecondMax;
 
   // Return the z-bin position of the highest density sum
   if (secondSum > firstSum && secondSum > thirdSum) {
-    return zSecondMax;
+    return zBinSecondMax;
   }
   if (thirdSum > secondSum && thirdSum > firstSum) {
-    return zThirdMax;
+    return zBinThirdMax;
   }
-  return zFirstMax;
+  return zBinFirstMax;
 }
 
 template <int trkGridSize>

Core/include/Acts/Vertexing/GaussianGridTrackDensity.hpp

@@ -130,18 +130,18 @@ class GaussianGridTrackDensity {
   /// @brief Function that creates a 1-dim track grid (i.e. a vector)
   /// with the correct density contribution of a track along the z-axis
   ///
-  /// @param offset Offset in d0 direction, to account for the 2-dim part
-  /// of the Gaussian track distribution
-  /// @param cov The track covariance matrix
-  /// @param distCtrD The distance in d0 from the track position to its
-  /// bin center in the 2-dim grid
+  /// @param d0 Transverse impact parameter
   /// @param distCtrZ The distance in z0 from the track position to its
   /// bin center in the 2-dim grid
-  TrackGridVector createTrackGrid(int offset, const SquareMatrix2& cov,
-                                  float distCtrD, float distCtrZ) const;
+  /// @param cov The track covariance matrix
+  TrackGridVector createTrackGrid(float d0, float distCtrZ,
+                                  const Acts::SquareMatrix2& cov) const;
 
   /// @brief Function that estimates the seed width based on the FWHM of
   /// the maximum density peak
+  /// @note This only works if the maximum is sufficiently isolated since
+  /// overlapping neighboring peaks might lead to an overestimation of the
+  /// seed width.
   ///
   /// @param mainGrid The main 1-dim density grid along the z-axis
   /// @param maxZ z-position of the maximum density value
@@ -150,6 +150,7 @@ class GaussianGridTrackDensity {
   Result<float> estimateSeedWidth(MainGridVector& mainGrid, float maxZ) const;
 
   /// @brief Helper to retrieve values according to a 2-dim normal distribution
+  /// @note This function is defined in coordinate system centered around d0 and z0
   float normal2D(float d, float z, const SquareMatrix2& cov) const;
 
   /// @brief Checks the (up to) first three density maxima (only those that have