Merge branch 'main' into feat-index-grid-navigation-policy

Author: kodiakhq[bot]

CI/physmon/workflows/physmon_trackrefitting_kf.py

@@ -24,7 +24,7 @@
         field=setup.field,
         digiConfigFile=setup.digiConfig,
         outputDir=tp,
-        reverseFilteringMomThreshold=float("inf"),
+        reverseFilteringMomThreshold=float("inf"),  # use reverse kf smoothing
         reverseFilteringCovarianceScaling=100.0,
         s=s,
     )

Core/include/Acts/Propagator/DirectNavigator.hpp

@@ -296,26 +296,33 @@ class DirectNavigator {
     // Move the sequence to the next surface
     state.nextSurface();
 
-    if (!state.endOfSurfaces()) {
-      ACTS_VERBOSE("Next surface candidate is  "
+    while (!state.endOfSurfaces()) {
+      ACTS_VERBOSE("Next surface candidate is "
                    << state.navSurface().geometryId() << ". "
                    << state.remainingSurfaces() << " out of "
                    << state.options.surfaces.size()
                    << " surfaces remain to try.");
-    } else {
-      ACTS_VERBOSE("End of surfaces reached, navigation break.");
-      state.navigationBreak = true;
-      return NavigationTarget::None();
+
+      // Establish & update the surface status
+      // TODO we do not know the intersection index - passing the closer one
+      const Surface& surface = state.navSurface();
+      const double farLimit = std::numeric_limits<double>::max();
+      const NavigationTarget target = chooseIntersection(
+          state.options.geoContext, surface, position, direction,
+          BoundaryTolerance::Infinite(), state.options.nearLimit, farLimit,
+          state.options.surfaceTolerance);
+      if (target.isValid()) {
+        return target;
+      }
+
+      ACTS_VERBOSE("No valid intersection found with surface "
+                   << surface.geometryId() << ", trying next surface.");
+      state.nextSurface();
     }
 
-    // Establish & update the surface status
-    // TODO we do not know the intersection index - passing the closer one
-    const Surface& surface = state.navSurface();
-    const double farLimit = std::numeric_limits<double>::max();
-    return chooseIntersection(state.options.geoContext, surface, position,
-                              direction, BoundaryTolerance::Infinite(),
-                              state.options.nearLimit, farLimit,
-                              state.options.surfaceTolerance);
+    ACTS_VERBOSE("End of surfaces reached, navigation break.");
+    state.navigationBreak = true;
+    return NavigationTarget::None();
   }
 
   /// @brief Check if the current target is still valid
@@ -361,7 +368,7 @@ class DirectNavigator {
 
     // Set the current surface
     state.currentSurface = &state.navSurface();
-    ACTS_VERBOSE("Current surface set to  "
+    ACTS_VERBOSE("Current surface set to "
                  << state.currentSurface->geometryId());
   }
 
@@ -376,7 +383,8 @@ class DirectNavigator {
 
     for (auto [intersectionIndex, intersection] :
          Acts::enumerate(intersections)) {
-      if (detail::checkPathLength(intersection.pathLength(), nearLimit,
+      if (intersection.isValid() &&
+          detail::checkPathLength(intersection.pathLength(), nearLimit,
                                   farLimit, logger())) {
         return NavigationTarget(intersection, intersectionIndex, surface,
                                 boundaryTolerance);

Core/include/Acts/Utilities/ProtoAxisHelpers.hpp

@@ -0,0 +1,127 @@
+// This file is part of the ACTS project.
+//
+// Copyright (C) 2016 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 https://mozilla.org/MPL/2.0/.
+#pragma once
+
+#include "Acts/Definitions/Algebra.hpp"
+#include "Acts/Utilities/BinningData.hpp"
+#include "Acts/Utilities/ProtoAxis.hpp"
+
+#include <span>
+#include <string>
+#include <vector>
+
+namespace Acts::ProtoAxisHelpers {
+
+/// @brief Get the number of bins from a ProtoAxis
+/// @param axis DirectedProtoAxis object
+/// @return Number of bins in the axis
+inline std::size_t binsOfProtoAxis(const Acts::DirectedProtoAxis& axis) {
+  return axis.getAxis().getNBins();
+}
+
+/// @brief Get the total number of bins from multiple ProtoAxes
+/// @param axes Span of DirectedProtoAxis objects
+/// @return Total number of bins across all axes
+inline std::size_t totalBinsFromProtoAxes(
+    std::span<const Acts::DirectedProtoAxis> axes) {
+  if (axes.empty() || axes.size() > 3) {
+    throw std::runtime_error(
+        "Unsupported number of axes, must be 1-3, instead got " +
+        std::to_string(axes.size()) + ")");
+  }
+  return axes[0].getAxis().getNBins() *
+         (axes.size() > 1 ? axes[1].getAxis().getNBins() : 1) *
+         (axes.size() > 2 ? axes[2].getAxis().getNBins() : 1);
+}
+
+/// @brief Get the number of bins from a specific ProtoAxis in a collection
+/// @param axes DirectedProtoAxis span
+/// @param ba Bin axis index
+/// @return Number of bins in the specified axis
+inline std::size_t binsFromProtoAxes(
+    std::span<const Acts::DirectedProtoAxis> axes, std::size_t ba) {
+  if (axes.empty() || axes.size() > 3) {
+    throw std::runtime_error(
+        "Unsupported number of axes, must be 1-3, instead got " +
+        std::to_string(axes.size()) + ")");
+  }
+  Acts::BinningData bd(axes[ba]);
+  return bd.bins();
+}
+
+/// @brief Get the bin index from a ProtoAxis using local coordinates
+/// @param axis DirectedProtoAxis object
+/// @param lp Local position vector
+/// @return Bin index corresponding to the local position
+inline std::size_t binFromProtoAxis(const Acts::DirectedProtoAxis& axis,
+                                    const Acts::Vector2& lp) {
+  Acts::BinningData bd(axis);
+  return bd.searchLocal(lp);
+}
+
+/// @brief Get the bin index from a ProtoAxis using global coordinates
+/// @param axis DirectedProtoAxis object
+/// @param gp Global position vector
+/// @return Bin index corresponding to the global position
+inline std::size_t binFromProtoAxis(const Acts::DirectedProtoAxis& axis,
+                                    const Acts::Vector3& gp) {
+  Acts::BinningData bd(axis);
+  return bd.searchGlobal(gp);
+}
+
+/// @brief Get the bin triple from multiple ProtoAxes using global coordinates
+/// @param axes Span of DirectedProtoAxis objects
+/// @param gp Global position vector
+/// @return Array of bin indices corresponding to the global position for each axis
+inline std::array<std::size_t, 3> binTripleFromProtoAxes(
+    std::span<const Acts::DirectedProtoAxis> axes, const Acts::Vector3& gp) {
+  const Acts::Vector3& bPosition = gp;
+  std::array<std::size_t, 3> bTriple = {0, 0, 0};
+  if (axes.empty() || axes.size() > 3) {
+    throw std::runtime_error(
+        "Unsupported number of axes, must be 1-3, instead got " +
+        std::to_string(axes.size()) + ")");
+  }
+  if (axes.size() == 1) {
+    Acts::BinningData bd0(axes[0]);
+    bTriple[0] = bd0.searchGlobal(bPosition);
+  }
+  if (axes.size() == 2) {
+    Acts::BinningData bd1(axes[1]);
+    bTriple[1] = bd1.searchGlobal(bPosition);
+  }
+  if (axes.size() == 3) {
+    Acts::BinningData bd2(axes[2]);
+    bTriple[2] = bd2.searchGlobal(bPosition);
+  }
+  return bTriple;
+}
+
+/// @brief Get the maximum bin index from a specific ProtoAxis in a collection
+/// @param axes DirectedProtoAxis span
+/// @param ba Bin axis index
+/// @return Maximum bin index in the specified axis
+inline std::size_t maxBin(std::span<const Acts::DirectedProtoAxis> axes,
+                          std::size_t ba = 0) {
+  if (axes.empty() || axes.size() > 3) {
+    throw std::runtime_error(
+        "Unsupported number of axes, must be 1-3, instead got " +
+        std::to_string(axes.size()) + ")");
+  }
+  std::vector<Acts::BinningData> binningDataVec;
+  binningDataVec.reserve(axes.size());
+  for (const auto& axis : axes) {
+    binningDataVec.emplace_back(axis);
+  }
+  if (ba >= binningDataVec.size()) {
+    return 0;
+  }
+  return (binningDataVec.at(ba).bins() - 1);
+}
+
+}  // namespace Acts::ProtoAxisHelpers

Core/src/Material/AccumulatedSurfaceMaterial.cpp

@@ -10,6 +10,7 @@
 
 #include "Acts/Material/BinnedSurfaceMaterial.hpp"
 #include "Acts/Material/HomogeneousSurfaceMaterial.hpp"
+#include "Acts/Utilities/ProtoAxisHelpers.hpp"
 
 #include <utility>
 

Examples/Algorithms/TrackFitting/src/RefittingAlgorithm.cpp

@@ -124,6 +124,8 @@ ProcessCode RefittingAlgorithm::execute(const AlgorithmContext& ctx) const {
       if (refittedTrack.hasReferenceSurface()) {
         ACTS_VERBOSE("Refitted parameters for track " << itrack);
         ACTS_VERBOSE("  " << track.parameters().transpose());
+        ACTS_VERBOSE("Measurements: " << refittedTrack.nMeasurements());
+        ACTS_VERBOSE("Outliers: " << refittedTrack.nOutliers());
       } else {
         ACTS_DEBUG("No refitted parameters for track " << itrack);
       }
@@ -135,9 +137,13 @@ ProcessCode RefittingAlgorithm::execute(const AlgorithmContext& ctx) const {
     ++itrack;
   }
 
-  std::stringstream ss;
-  trackStateContainer->statistics().toStream(ss);
-  ACTS_DEBUG(ss.str());
+  ACTS_DEBUG("Fitted tracks: " << trackContainer->size());
+
+  if (logger().doPrint(Acts::Logging::DEBUG)) {
+    std::stringstream ss;
+    trackStateContainer->statistics().toStream(ss);
+    ACTS_DEBUG(ss.str());
+  }
 
   ConstTrackContainer constTracks{
       std::make_shared<Acts::ConstVectorTrackContainer>(

Examples/Algorithms/TrackFitting/src/TrackFittingAlgorithm.cpp

@@ -171,6 +171,8 @@ ProcessCode TrackFittingAlgorithm::execute(const AlgorithmContext& ctx) const {
     }
   }
 
+  ACTS_DEBUG("Fitted tracks: " << trackContainer->size());
+
   if (logger().doPrint(Acts::Logging::DEBUG)) {
     std::stringstream ss;
     trackStateContainer->statistics().toStream(ss);

Examples/Scripts/Python/truth_tracking_gsf_refitting.py

@@ -21,6 +21,7 @@ def runRefittingGsf(
     field: acts.MagneticFieldProvider,
     digiConfigFile: Path,
     outputDir: Path,
+    reverseFilteringCovarianceScaling=100.0,
     inputParticlePath: Optional[Path] = None,
     inputSimHitsPath: Optional[Path] = None,
     decorators=[],
@@ -38,8 +39,8 @@ def runRefittingGsf(
         inputHitsPath=inputSimHitsPath,
         decorators=decorators,
         generatedParticleType=acts.PdgParticle.eElectron,
-        reverseFilteringMomThreshold=0.0,
-        reverseFilteringCovarianceScaling=1.0,
+        reverseFilteringMomThreshold=0 * u.GeV,  # use direct smoothing
+        reverseFilteringCovarianceScaling=reverseFilteringCovarianceScaling,
         s=s,
     )
 
@@ -54,7 +55,7 @@ def runRefittingGsf(
         "componentMergeMethod": acts.examples.ComponentMergeMethod.maxWeight,
         "mixtureReductionAlgorithm": acts.examples.MixtureReductionAlgorithm.KLDistance,
         "weightCutoff": 1.0e-4,
-        "reverseFilteringCovarianceScaling": 100.0,
+        "reverseFilteringCovarianceScaling": reverseFilteringCovarianceScaling,
         "level": acts.logging.INFO,
     }
 

Examples/Scripts/Python/truth_tracking_kalman_refitting.py

@@ -22,16 +22,16 @@ def runRefittingKf(
     outputDir: Path,
     multipleScattering: bool = True,
     energyLoss: bool = True,
-    reverseFilteringMomThreshold=0 * u.GeV,
-    reverseFilteringCovarianceScaling=1.0,
+    reverseFilteringMomThreshold=float("inf"),
+    reverseFilteringCovarianceScaling=100.0,
     s: acts.examples.Sequencer = None,
 ):
     s = runTruthTrackingKalman(
         trackingGeometry,
         field,
         digiConfigFile=digiConfigFile,
         outputDir=outputDir,
-        reverseFilteringMomThreshold=reverseFilteringMomThreshold,
+        reverseFilteringMomThreshold=0 * u.GeV,  # use direct smoothing
         reverseFilteringCovarianceScaling=reverseFilteringCovarianceScaling,
         s=s,
     )