ci: Monitor grid seeding (#2355)

Since we will extend the adaptive grid seeder to accommodate time seeding, I propose to add it to the physmon.

CI/physmon/phys_perf_mon.sh

@@ -172,6 +172,26 @@ function full_chain() {
         "AMVF ${suffix}" \
         amvf_${suffix}
 
+    if [ $suffix == seeded ]; then
+	    Examples/Scripts/generic_plotter.py \
+            $outdir/performance_amvf_gridseeder_${suffix}.root \
+            vertexing \
+            $outdir/performance_amvf_gridseeder_${suffix}_hist.root \
+            --silent \
+            --config CI/physmon/vertexing_config.yml
+        ec=$(($ec | $?))
+
+        # remove ntuple file because it's large
+        rm $outdir/performance_amvf_gridseeder_${suffix}.root
+
+        run \
+            $outdir/performance_amvf_gridseeder_${suffix}_hist.root \
+            $refdir/performance_amvf_gridseeder_${suffix}_hist.root \
+            --title "AMVF (+grid seeder) ${suffix}" \
+            -o $outdir/amvf_gridseeder_${suffix}.html \
+            -p $outdir/amvf_gridseeder_${suffix}_plots
+    fi
+
     Examples/Scripts/generic_plotter.py \
         $outdir/tracksummary_ckf_${suffix}.root \
         tracksummary \
@@ -285,6 +305,23 @@ if [[ "$mode" == "all" || "$mode" == "fullchains" ]]; then
         $refdir/performance_amvf_ttbar_hist.root \
         "AMVF ttbar" \
         amvf_ttbar
+
+    Examples/Scripts/generic_plotter.py \
+        $outdir/performance_amvf_gridseeder_ttbar.root \
+        vertexing \
+        $outdir/performance_amvf_gridseeder_ttbar_hist.root \
+        --silent \
+        --config CI/physmon/vertexing_config.yml
+    ec=$(($ec | $?))
+
+    # remove ntuple file because it's large
+    rm $outdir/performance_amvf_gridseeder_ttbar.root
+
+    run \
+        $outdir/performance_amvf_gridseeder_ttbar_hist.root \
+        $refdir/performance_amvf_gridseeder_ttbar_hist.root \
+        "AMVF (+grid seeder) ttbar" \
+        amvf_gridseeder_ttbar
 fi
 
 if [[ "$mode" == "all" || "$mode" == "gsf" ]]; then

CI/physmon/workflows/physmon_ckf_tracking.py

@@ -142,6 +142,7 @@ def run_ckf_tracking(truthSmearedSeeded, truthEstimatedSeeded, label):
         addVertexFitting(
             s,
             setup.field,
+            seeder=acts.VertexSeedFinder.GaussianSeeder,
             associatedParticles=None
             if label in ["seeded", "orthogonal"]
             else "particles_input",
@@ -154,6 +155,7 @@ def run_ckf_tracking(truthSmearedSeeded, truthEstimatedSeeded, label):
         addVertexFitting(
             s,
             setup.field,
+            seeder=acts.VertexSeedFinder.GaussianSeeder,
             associatedParticles=None
             if label in ["seeded", "orthogonal"]
             else "particles_input",
@@ -163,6 +165,20 @@ def run_ckf_tracking(truthSmearedSeeded, truthEstimatedSeeded, label):
             outputDirRoot=tp / "amvf",
         )
 
+        # Use the adaptive grid vertex seeder in combination with the AMVF
+        # To avoid having too many physmon cases, we only do this for the label "seeded"
+        if label == "seeded":
+            addVertexFitting(
+                s,
+                setup.field,
+                seeder=acts.VertexSeedFinder.AdaptiveGridSeeder,
+                associatedParticles=None,
+                outputProtoVertices="amvf_gridseeder_protovertices",
+                outputVertices="amvf_gridseeder_fittedVertices",
+                vertexFinder=VertexFinder.AMVF,
+                outputDirRoot=tp / "amvf_gridseeder",
+            )
+
         s.run()
         del s
 
@@ -172,17 +188,28 @@ def run_ckf_tracking(truthSmearedSeeded, truthEstimatedSeeded, label):
                 tp / f"performance_{vertexing}.root",
             )
 
-        for stem in [
-            "performance_ckf",
-            "tracksummary_ckf",
-            "performance_ivf",
-            "performance_amvf",
-        ] + (
-            ["performance_seeding", "performance_ambi"]
-            if label in ["seeded", "orthogonal"]
-            else ["performance_seeding"]
-            if label == "truth_estimated"
-            else []
+        if label == "seeded":
+            vertexing = "amvf_gridseeder"
+            shutil.move(
+                tp / f"{vertexing}/performance_vertexing.root",
+                tp / f"performance_{vertexing}.root",
+            )
+
+        for stem in (
+            [
+                "performance_ckf",
+                "tracksummary_ckf",
+                "performance_ivf",
+                "performance_amvf",
+            ]
+            + (["performance_amvf_gridseeder"] if label == "seeded" else [])
+            + (
+                ["performance_seeding", "performance_ambi"]
+                if label in ["seeded", "orthogonal"]
+                else ["performance_seeding"]
+                if label == "truth_estimated"
+                else []
+            )
         ):
             perf_file = tp / f"{stem}.root"
             assert perf_file.exists(), "Performance file not found"

CI/physmon/workflows/physmon_track_finding_ttbar.py

@@ -154,26 +154,39 @@
     addVertexFitting(
         s,
         setup.field,
+        seeder=acts.VertexSeedFinder.GaussianSeeder,
         associatedParticles=None,
         outputProtoVertices="amvf_protovertices",
         outputVertices="amvf_fittedVertices",
         vertexFinder=VertexFinder.AMVF,
         outputDirRoot=tp / "amvf",
     )
 
+    addVertexFitting(
+        s,
+        setup.field,
+        seeder=acts.VertexSeedFinder.AdaptiveGridSeeder,
+        associatedParticles=None,
+        outputProtoVertices="amvf_gridseeder_protovertices",
+        outputVertices="amvf_gridseeder_fittedVertices",
+        vertexFinder=VertexFinder.AMVF,
+        outputDirRoot=tp / "amvf_gridseeder",
+    )
+
     s.run()
     del s
 
-    vertexing = "amvf"
-    shutil.move(
-        tp / f"{vertexing}/performance_vertexing.root",
-        tp / f"performance_{vertexing}.root",
-    )
+    for vertexing in ["amvf", "amvf_gridseeder"]:
+        shutil.move(
+            tp / f"{vertexing}/performance_vertexing.root",
+            tp / f"performance_{vertexing}.root",
+        )
 
     for stem in [
         "performance_ckf",
         "tracksummary_ckf",
         "performance_amvf",
+        "performance_amvf_gridseeder",
     ] + (["performance_seeding", "performance_ambi"]):
         perf_file = tp / f"{stem}.root"
         assert perf_file.exists(), "Performance file not found"

Core/include/Acts/Utilities/AlgebraHelpers.hpp

@@ -167,14 +167,16 @@ inline ActsMatrix<A::RowsAtCompileTime, B::ColsAtCompileTime> blockedMult(
   }
 }
 
-/// Calculate the inverse of an Eigen matrix after checking if it can be
-/// numerically inverted. This allows to catch potential FPEs before they occur.
+/// FPE "safe" functions
 ///
 /// Our main motivation for this is that users might have a strict FPE policy
 /// which would flag every single occurrence as a failure and then sombody has
 /// to investigate. Since we are processing a high number of events and floating
 /// point numbers sometimes work in mysterious ways the caller of this function
 /// might want to hide FPEs and handle them in a more controlled way.
+
+/// Calculate the inverse of an Eigen matrix after checking if it can be
+/// numerically inverted. This allows to catch potential FPEs before they occur.
 ///
 /// @tparam Derived Eigen derived concrete type
 /// @tparam Result Eigen result type defaulted to input type
@@ -196,4 +198,25 @@ std::optional<ResultType> safeInverse(const MatrixType& m) noexcept {
   return std::nullopt;
 }
 
+/// Calculate the exponential function while avoiding FPEs.
+/// @note The boundary values of -50.0 and 50.0 might need to be adapted when
+/// using this function with doubles
+///
+/// @param val argument for which the exponential function should be evaluated.
+///
+/// @return 0 in the case of underflow, std::numeric_limits<T>::infinity in the
+/// case of overflow, std::exp(val) else
+template <typename T>
+T safeExp(T val) noexcept {
+  if (val < -50.0) {
+    return 0.0;
+  }
+
+  if (val > 50.0) {
+    return std::numeric_limits<T>::infinity();
+  }
+
+  return std::exp(val);
+}
+
 }  // namespace Acts

Core/include/Acts/Vertexing/AdaptiveGridTrackDensity.ipp

@@ -5,6 +5,7 @@
 // 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/Utilities/AlgebraHelpers.hpp"
 #include "Acts/Vertexing/VertexingError.hpp"
 
 #include <algorithm>
@@ -197,7 +198,7 @@ float Acts::AdaptiveGridTrackDensity<trkGridSize>::normal2D(
   float expo =
       -1 / (2 * det) *
       (cov(1, 1) * d * d - d * z * (cov(0, 1) + cov(1, 0)) + cov(0, 0) * z * z);
-  return coef * std::exp(expo);
+  return coef * safeExp(expo);
 }
 
 template <int trkGridSize>

Core/include/Acts/Vertexing/GaussianGridTrackDensity.ipp

@@ -5,6 +5,7 @@
 // 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/Utilities/AlgebraHelpers.hpp"
 #include "Acts/Vertexing/VertexingError.hpp"
 
 #include <algorithm>
@@ -207,7 +208,7 @@ float Acts::GaussianGridTrackDensity<mainGridSize, trkGridSize>::normal2D(
   float expo =
       -1 / (2 * det) *
       (cov(1, 1) * d * d - d * z * (cov(0, 1) + cov(1, 0)) + cov(0, 0) * z * z);
-  return coef * std::exp(expo);
+  return coef * safeExp(expo);
 }
 
 template <int mainGridSize, int trkGridSize>

Examples/Algorithms/Vertexing/include/ActsExamples/Vertexing/AdaptiveMultiVertexFinderAlgorithm.hpp

@@ -17,6 +17,7 @@
 #include "Acts/Propagator/Propagator.hpp"
 #include "Acts/Utilities/Helpers.hpp"
 #include "Acts/Utilities/Logger.hpp"
+#include "Acts/Vertexing/AdaptiveGridDensityVertexFinder.hpp"
 #include "Acts/Vertexing/AdaptiveMultiVertexFinder.hpp"
 #include "Acts/Vertexing/AdaptiveMultiVertexFitter.hpp"
 #include "Acts/Vertexing/GaussianTrackDensity.hpp"
@@ -56,14 +57,13 @@ class AdaptiveMultiVertexFinderAlgorithm final : public IAlgorithm {
   using Linearizer = Acts::HelicalTrackLinearizer<Propagator>;
   using Fitter =
       Acts::AdaptiveMultiVertexFitter<Acts::BoundTrackParameters, Linearizer>;
-  using Seeder = Acts::TrackDensityVertexFinder<
-      Fitter, Acts::GaussianTrackDensity<Acts::BoundTrackParameters>>;
-  using Finder = Acts::AdaptiveMultiVertexFinder<Fitter, Seeder>;
   using Options = Acts::VertexingOptions<Acts::BoundTrackParameters>;
 
   using VertexCollection =
       std::vector<Acts::Vertex<Acts::BoundTrackParameters>>;
 
+  enum class SeedFinder { GaussianSeeder, AdaptiveGridSeeder };
+
   struct Config {
     /// Optional. Input track parameters collection
     std::string inputTrackParameters;
@@ -73,18 +73,31 @@ class AdaptiveMultiVertexFinderAlgorithm final : public IAlgorithm {
     std::string outputProtoVertices;
     /// Output vertex collection
     std::string outputVertices = "vertices";
+    /// Enum member determining the choice of the vertex seed finder
+    SeedFinder seedFinder;
     /// The magnetic field
     std::shared_ptr<Acts::MagneticFieldProvider> bField;
   };
 
   AdaptiveMultiVertexFinderAlgorithm(const Config& config,
                                      Acts::Logging::Level level);
-  /// Find vertices using the adaptive multi vertex finder algorithm.
+
+  /// Set up vertex seeder and call the function executeAfterSeederChoice.
   ///
   /// @param ctx is the algorithm context with event information
   /// @return a process code indication success or failure
   ProcessCode execute(const AlgorithmContext& ctx) const final;
 
+  /// Find vertices using the adaptive multi vertex finder algorithm.
+  ///
+  /// @param ctx is the algorithm context with event information
+  /// @param seedFinder is the vertex seed finder
+  /// @return a process code indication success or failure
+  template <typename vseeder_t, typename vfinder_t>
+  ProcessCode executeAfterSeederChoice(
+      const ActsExamples::AlgorithmContext& ctx,
+      const vseeder_t& seedFinder) const;
+
   /// Get readonly access to the config parameters
   const Config& config() const { return m_cfg; }
 

Examples/Algorithms/Vertexing/src/AdaptiveMultiVertexFinderAlgorithm.cpp

@@ -57,6 +57,31 @@ ActsExamples::AdaptiveMultiVertexFinderAlgorithm::
 ActsExamples::ProcessCode
 ActsExamples::AdaptiveMultiVertexFinderAlgorithm::execute(
     const ActsExamples::AlgorithmContext& ctx) const {
+  if (m_cfg.seedFinder == SeedFinder::GaussianSeeder) {
+    using Seeder = Acts::TrackDensityVertexFinder<
+        Fitter, Acts::GaussianTrackDensity<Acts::BoundTrackParameters>>;
+    using Finder = Acts::AdaptiveMultiVertexFinder<Fitter, Seeder>;
+    Seeder seedFinder;
+    return executeAfterSeederChoice<Seeder, Finder>(ctx, seedFinder);
+  } else if (m_cfg.seedFinder == SeedFinder::AdaptiveGridSeeder) {
+    using Seeder = Acts::AdaptiveGridDensityVertexFinder<109, Fitter>;
+    using Finder = Acts::AdaptiveMultiVertexFinder<Fitter, Seeder>;
+    // The seeder config argument corresponds to the bin size in mm
+    Seeder::Config seederConfig(0.05);
+    Seeder seedFinder(seederConfig);
+    return executeAfterSeederChoice<Seeder, Finder>(ctx, seedFinder);
+  } else {
+    return ActsExamples::ProcessCode::ABORT;
+  }
+}
+
+template <typename vseeder_t, typename vfinder_t>
+ActsExamples::ProcessCode
+ActsExamples::AdaptiveMultiVertexFinderAlgorithm::executeAfterSeederChoice(
+    const ActsExamples::AlgorithmContext& ctx,
+    const vseeder_t& seedFinder) const {
+  using Finder = vfinder_t;
+
   // Set up EigenStepper
   Acts::EigenStepper<> stepper(m_cfg.bField);
 
@@ -83,11 +108,9 @@ ActsExamples::AdaptiveMultiVertexFinderAlgorithm::execute(
   fitterCfg.doSmoothing = true;
   Fitter fitter(fitterCfg, logger().cloneWithSuffix("AMVFitter"));
 
-  // Set up the vertex seed finder
-  Seeder seedFinder;
-
-  Finder::Config finderConfig(std::move(fitter), seedFinder, ipEstimator,
-                              std::move(linearizer), m_cfg.bField);
+  typename Finder::Config finderConfig(std::move(fitter), seedFinder,
+                                       ipEstimator, std::move(linearizer),
+                                       m_cfg.bField);
   // We do not want to use a beamspot constraint here
   finderConfig.useBeamSpotConstraint = false;
   finderConfig.tracksMaxZinterval = 1. * Acts::UnitConstants::mm;
@@ -121,7 +144,7 @@ ActsExamples::AdaptiveMultiVertexFinderAlgorithm::execute(
   //////////////////////////////////////////////
 
   // The vertex finder state
-  Finder::State state;
+  typename Finder::State state;
 
   // Default vertexing options, this is where e.g. a constraint could be set
   Options finderOpts(ctx.geoContext, ctx.magFieldContext);

Examples/Python/python/acts/examples/reconstruction.py

@@ -1436,6 +1436,7 @@ def addAmbiguityResolutionMLDBScan(
 def addVertexFitting(
     s,
     field,
+    seeder: Optional[acts.VertexSeedFinder] = acts.VertexSeedFinder.GaussianSeeder,
     trajectories: Optional[str] = "trajectories",
     trackParameters: Optional[str] = None,
     associatedParticles: Optional[str] = None,
@@ -1453,6 +1454,8 @@ def addVertexFitting(
     s: Sequencer
         the sequencer module to which we add the Seeding steps (returned from addVertexFitting)
     field : magnetic field
+    seeder : enum member
+        determines vertex seeder, can be acts.seeder.GaussianSeeder or acts.seeder.AdaptiveGridSeeder
     outputDirRoot : Path|str, path, None
         the output folder for the Root output, None triggers no output
     associatedParticles : str, "associatedTruthParticles"
@@ -1525,6 +1528,7 @@ def addVertexFitting(
     elif vertexFinder == VertexFinder.AMVF:
         findVertices = AdaptiveMultiVertexFinderAlgorithm(
             level=customLogLevel(),
+            seedFinder=seeder,
             bField=field,
             inputTrajectories=trajectories,
             inputTrackParameters=trackParameters,
@@ -1571,7 +1575,6 @@ def addSingleSeedVertexFinding(
     inputSpacePoints: Optional[str] = "spacepoints",
     outputVertices: Optional[str] = "fittedSeedVertices",
 ) -> None:
-
     from acts.examples import (
         SingleSeedVertexFinderAlgorithm,
         VertexPerformanceWriter,