feat: Add check for SurfaceArray grid + transform handling improvement

Author: paulgessinger

Core/include/Acts/Surfaces/SurfaceArray.hpp

@@ -558,6 +558,8 @@ class SurfaceArray {
   /// @return Reference to the surface grid lookup interface
   const ISurfaceGridLookup& gridLookup() const { return *p_gridLookup; }
 
+  void checkGrid(AnyGridConstView<SurfaceVector> grid);
+
  private:
   std::unique_ptr<ISurfaceGridLookup> p_gridLookup;
   // this vector makes sure we have shared ownership over the surfaces

Core/src/Geometry/SurfaceArrayCreator.cpp

@@ -57,6 +57,7 @@ std::unique_ptr<SurfaceArray> SurfaceArrayCreator::surfaceArrayOnCylinder(
   const double layerTolerance = protoLayer.range(AxisDirection::AxisR) * 0.5;
 
   auto surface = Surface::makeShared<CylinderSurface>(fullTransform, R, halfZ);
+  ACTS_VERBOSE("- projection surface is: " << surface->toString(gctx));
   std::unique_ptr<SurfaceArray::ISurfaceGridLookup> sl =
       makeSurfaceGridLookup2D<AxisBoundaryType::Closed,
                               AxisBoundaryType::Bound>(
@@ -150,8 +151,18 @@ std::unique_ptr<SurfaceArray> SurfaceArrayCreator::surfaceArrayOnDisc(
   const double Rmax = protoLayer.max(AxisDirection::AxisR, true);
   const double layerThickness = protoLayer.range(AxisDirection::AxisZ) * 0.5;
   ACTS_VERBOSE("- z-position of disc estimated as " << Z);
+  ACTS_VERBOSE("- full transform is \n" << fullTransform.matrix());
+
+  if (fullTransform.translation().norm() < s_transformEquivalentTolerance) {
+    ACTS_VERBOSE(
+        "input transform does not have translation: putting projection surface "
+        "at center of gravity in z");
+    fullTransform.translate(Vector3::UnitZ() * Z);
+  }
 
   auto surface = Surface::makeShared<DiscSurface>(fullTransform, Rmin, Rmax);
+  ACTS_VERBOSE("- projection surface is: " << surface->toString(gctx));
+
   std::unique_ptr<SurfaceArray::ISurfaceGridLookup> sl =
       makeSurfaceGridLookup2D<AxisBoundaryType::Bound,
                               AxisBoundaryType::Closed>(
@@ -253,6 +264,13 @@ std::unique_ptr<SurfaceArray> SurfaceArrayCreator::surfaceArrayOnDisc(
   const double layerThickness = protoLayer.range(AxisDirection::AxisZ) * 0.5;
   ACTS_VERBOSE("- z-position of disc estimated as " << Z);
 
+  if (fullTransform.translation().norm() < s_transformEquivalentTolerance) {
+    ACTS_VERBOSE(
+        "input transform does not have translation: putting projection surface "
+        "at center of gravity in z");
+    fullTransform.translate(Vector3::UnitZ() * Z);
+  }
+
   auto surface = Surface::makeShared<DiscSurface>(fullTransform, Rmin, Rmax);
   std::unique_ptr<SurfaceArray::ISurfaceGridLookup> sl =
       makeSurfaceGridLookup2D<AxisBoundaryType::Bound,

Core/src/Navigation/SurfaceArrayNavigationPolicy.cpp

@@ -9,6 +9,7 @@
 #include "Acts/Navigation/SurfaceArrayNavigationPolicy.hpp"
 
 #include "Acts/Geometry/GeometryContext.hpp"
+#include "Acts/Geometry/ProtoLayer.hpp"
 #include "Acts/Geometry/SurfaceArrayCreator.hpp"
 #include "Acts/Geometry/TrackingVolume.hpp"
 #include "Acts/Navigation/NavigationStream.hpp"
@@ -118,7 +119,8 @@ SurfaceArrayNavigationPolicy::SurfaceArrayNavigationPolicy(
     throw std::runtime_error("Cannot create surface array with zero surfaces");
   }
 
-  ProtoLayer protoLayer(gctx, surfaces);
+  ProtoLayer protoLayer(
+      gctx, surfaces, Transform3{volume.localToGlobalTransform(gctx).linear()});
 
   if (config.layerType == LayerType::Disc) {
     auto [binsR, binsPhi] = config.bins;

Core/src/Surfaces/SurfaceArray.cpp

@@ -11,9 +11,12 @@
 #include "Acts/Definitions/Algebra.hpp"
 #include "Acts/Geometry/SurfaceArrayCreator.hpp"
 #include "Acts/Utilities/Helpers.hpp"
+#include "Acts/Utilities/Ranges.hpp"
 
+#include <ranges>
 #include <utility>
 
+using namespace Acts::Ranges;
 namespace Acts {
 
 // implementation for pure virtual destructor of ISurfaceGridLookup
@@ -25,7 +28,13 @@ SurfaceArray::SurfaceArray(std::unique_ptr<ISurfaceGridLookup> gridLookup,
     : p_gridLookup(std::move(gridLookup)),
       m_surfaces(std::move(surfaces)),
       m_surfacesRawPointers(unpackSmartPointers(m_surfaces)),
-      m_transform(transform) {}
+      m_transform(transform) {
+  if (p_gridLookup != nullptr) {
+    if (const auto& grid = p_gridLookup->getGridView()) {
+      checkGrid(grid.value());
+    }
+  }
+}
 
 SurfaceArray::SurfaceArray(std::shared_ptr<const Surface> srf)
     : p_gridLookup(std::make_unique<SingleElementLookup>(srf.get())),
@@ -74,4 +83,32 @@ std::ostream& SurfaceArray::toStream(const GeometryContext& /*gctx*/,
   return sl;
 }
 
+void SurfaceArray::checkGrid(AnyGridConstView<SurfaceVector> grid) {
+  std::set allSurfaces =
+      m_surfaces |
+      std::views::transform([](const auto& sp) { return sp.get(); }) |
+      to<std::set>;
+  std::set<const Surface*> seenSurface;
+  auto bins = grid.numLocalBins();
+  for (std::size_t i = 0; i <= bins.at(0); ++i) {
+    for (std::size_t j = 0; j <= bins.at(1); ++j) {
+      const auto& surfaces = grid.atLocalBins({i, j});
+      for (const auto& srf : surfaces) {
+        seenSurface.insert(srf);
+      }
+    }
+  }
+
+  if (allSurfaces != seenSurface) {
+    std::set<const Surface*> diff;
+    std::ranges::set_difference(allSurfaces, seenSurface,
+                                std::inserter(diff, diff.begin()));
+
+    throw std::logic_error(
+        std::format("SurfaceArray grid does not contain all surfaces provided! "
+                    "{} surfaces not seen",
+                    diff.size()));
+  }
+}
+
 }  // namespace Acts