feat: Delayed Grid construction for Portals (#3718)

This PR does three things:

- `PortalLinkBase::merge` **no longer** does merging of grids or trivials into a combined grid. This has been observed to lead to accumulating floating point imprecision.
- `CompositePortalLink` gets a method to make a grid **if** (and only if) it is composed of a set of `TrivialPortalLink`s.
- `Portal::fuse` will attempt to convert `CompositePortalLink`s composed of only `TrivialPortalLink`s to a single `GridPortalLink` before fusing it with the other portal.

In combination, this avoids the floating point issues and produces correctly sized grids.

Part of #3502.

---

This pull request introduces several enhancements and bug fixes to the `Core/include/Acts/Geometry` module, focusing on improving the `CompositePortalLink` and `GridPortalLink` classes. The most important changes include the addition of new constructors, the introduction of the `makeGrid` method, and various adjustments to ensure compatibility and correctness.

### Enhancements to `CompositePortalLink`:

* Added new constructors to allow the creation of composite portals from multiple links and to handle nested composites. (`Core/include/Acts/Geometry/CompositePortalLink.hpp`, [[1]](diffhunk://#diff-248145d68399a17324b82d81d6e661a3ab739eb5b6f8d67f40145195ca465c36R55-R63) [[2]](diffhunk://#diff-5663ec0ea1d9723e610725aa0d0964e5c833bb90f431281c3802a9b9c5fa4314R101-R129)
* Introduced the `makeGrid` method to convert composite portal links into grid portal links under specific conditions. (`Core/include/Acts/Geometry/CompositePortalLink.hpp`, [Core/src/Geometry/CompositePortalLink.cppR180-R297](diffhunk://#diff-5663ec0ea1d9723e610725aa0d0964e5c833bb90f431281c3802a9b9c5fa4314R180-R297))

### Adjustments to `GridPortalLink`:

* Changed the type of `atLocalBins` methods to return `const TrackingVolume*` instead of `TrackingVolume*`. (`Core/include/Acts/Geometry/GridPortalLink.hpp`, [[1]](diffhunk://#diff-5cc33f33e4da7753510a3c7bf2481d12c34dd9c1344bfd624c0b5db1d70f214fL384-R389) [[2]](diffhunk://#diff-5cc33f33e4da7753510a3c7bf2481d12c34dd9c1344bfd624c0b5db1d70f214fL402-R402) [[3]](diffhunk://#diff-5cc33f33e4da7753510a3c7bf2481d12c34dd9c1344bfd624c0b5db1d70f214fL548-R547) [[4]](diffhunk://#diff-5cc33f33e4da7753510a3c7bf2481d12c34dd9c1344bfd624c0b5db1d70f214fL560-R559)
* Updated internal grid type to use `const TrackingVolume*`. (`Core/include/Acts/Geometry/GridPortalLink.hpp`, [Core/include/Acts/Geometry/GridPortalLink.hppL402-R402](diffhunk://#diff-5cc33f33e4da7753510a3c7bf2481d12c34dd9c1344bfd624c0b5db1d70f214fL402-R402))

### Bug Fixes and Code Improvements:

* Moved `mergedSurface` function to an anonymous namespace in the implementation file and refactored it for better error handling and type safety. (`Core/src/Geometry/CompositePortalLink.cpp`, [[1]](diffhunk://#diff-5663ec0ea1d9723e610725aa0d0964e5c833bb90f431281c3802a9b9c5fa4314R11-R78) [[2]](diffhunk://#diff-5663ec0ea1d9723e610725aa0d0964e5c833bb90f431281c3802a9b9c5fa4314L77-L106)
* Improved the `isSameSurface` function to include detailed checks for surface bounds and transformations. (`Core/src/Geometry/Portal.cpp`, [Core/src/Geometry/Portal.cppL308-R357](diffhunk://#diff-e32791625fda93fd367fc971619ea03be19128e91bbca7e8a09b5af399beb461L308-R357))
* Enhanced logging and error messages for better debugging and clarity. (`Core/src/Geometry/Portal.cpp`, [[1]](diffhunk://#diff-e32791625fda93fd367fc971619ea03be19128e91bbca7e8a09b5af399beb461R274-R277) [[2]](diffhunk://#diff-e32791625fda93fd367fc971619ea03be19128e91bbca7e8a09b5af399beb461R293-R313)

These changes collectively improve the functionality, safety, and maintainability of the `Acts` geometry module, particularly in handling complex portal link structures.

Core/include/Acts/Geometry/CompositePortalLink.hpp

@@ -17,6 +17,9 @@
 
 namespace Acts {
 
+class GridPortalLink;
+class Surface;
+
 /// Composite portal links can graft together other portal link instances, for
 /// example grids that could not be merged due to invalid binnings.
 ///
@@ -49,6 +52,15 @@ class CompositePortalLink final : public PortalLinkBase {
                       std::unique_ptr<PortalLinkBase> b, BinningValue direction,
                       bool flatten = true);
 
+  /// Construct a composite portal from any number of arbitrary other portal
+  /// links. The only requirement is that the portal link surfaces are
+  /// mergeable.
+  /// @param links The portal links
+  /// @param direction The binning direction
+  /// @param flatten If true, the composite will flatten any nested composite
+  CompositePortalLink(std::vector<std::unique_ptr<PortalLinkBase>> links,
+                      BinningValue direction, bool flatten = true);
+
   /// Print the composite portal link
   /// @param os The output stream
   void toStream(std::ostream& os) const override;
@@ -81,19 +93,22 @@ class CompositePortalLink final : public PortalLinkBase {
   /// @return The number of children
   std::size_t size() const;
 
- private:
-  /// Helper function to construct a merged surface from two portal links along
-  /// a given direction
-  /// @param a The first portal link
-  /// @param b The second portal link
-  /// @param direction The merging direction
-  /// @return The merged surface
-  static std::shared_ptr<RegularSurface> mergedSurface(const PortalLinkBase* a,
-                                                       const PortalLinkBase* b,
-                                                       BinningValue direction);
+  /// (Potentially) create a grid portal link that represents this composite
+  /// portal link.
+  /// @note This only works, if the composite is **flat** and only contains
+  ///       **trivial portal links**. If these preconditions are not met, this
+  ///       function returns a nullptr.
+  /// @param gctx The geometry context
+  /// @param logger The logger
+  /// @return The grid portal link
+  std::unique_ptr<GridPortalLink> makeGrid(const GeometryContext& gctx,
+                                           const Logger& logger) const;
 
+ private:
   boost::container::small_vector<std::unique_ptr<PortalLinkBase>, 4>
       m_children{};
+
+  BinningValue m_direction;
 };
 
 }  // namespace Acts

Core/include/Acts/Geometry/GridPortalLink.hpp

@@ -381,12 +381,12 @@ class GridPortalLink : public PortalLinkBase {
   /// Helper function to get grid bin content in type-eraased way.
   /// @param indices The bin indices
   /// @return The tracking volume at the bin
-  virtual TrackingVolume*& atLocalBins(IndexType indices) = 0;
+  virtual const TrackingVolume*& atLocalBins(IndexType indices) = 0;
 
   /// Helper function to get grid bin content in type-eraased way.
   /// @param indices The bin indices
   /// @return The tracking volume at the bin
-  virtual TrackingVolume* atLocalBins(IndexType indices) const = 0;
+  virtual const TrackingVolume* atLocalBins(IndexType indices) const = 0;
 
  private:
   BinningValue m_direction;
@@ -399,7 +399,7 @@ template <typename... Axes>
 class GridPortalLinkT final : public GridPortalLink {
  public:
   /// The internal grid type
-  using GridType = Grid<TrackingVolume*, Axes...>;
+  using GridType = Grid<const TrackingVolume*, Axes...>;
 
   /// The dimension of the grid
   static constexpr std::size_t DIM = sizeof...(Axes);
@@ -530,7 +530,6 @@ class GridPortalLinkT final : public GridPortalLink {
     return m_grid.atPosition(m_projection(position));
   }
 
- protected:
   /// Type erased access to the number of bins
   /// @return The number of bins in each direction
   IndexType numLocalBins() const override {
@@ -545,7 +544,7 @@ class GridPortalLinkT final : public GridPortalLink {
   /// Type erased local bin access
   /// @param indices The bin indices
   /// @return The tracking volume at the bin
-  TrackingVolume*& atLocalBins(IndexType indices) override {
+  const TrackingVolume*& atLocalBins(IndexType indices) override {
     throw_assert(indices.size() == DIM, "Invalid number of indices");
     typename GridType::index_t idx;
     for (std::size_t i = 0; i < DIM; i++) {
@@ -557,7 +556,7 @@ class GridPortalLinkT final : public GridPortalLink {
   /// Type erased local bin access
   /// @param indices The bin indices
   /// @return The tracking volume at the bin
-  TrackingVolume* atLocalBins(IndexType indices) const override {
+  const TrackingVolume* atLocalBins(IndexType indices) const override {
     throw_assert(indices.size() == DIM, "Invalid number of indices");
     typename GridType::index_t idx;
     for (std::size_t i = 0; i < DIM; i++) {

Core/include/Acts/Geometry/TrivialPortalLink.hpp

@@ -58,6 +58,10 @@ class TrivialPortalLink final : public PortalLinkBase {
       const GeometryContext& gctx, const Vector3& position,
       double tolerance = s_onSurfaceTolerance) const override;
 
+  /// Get the single volume that this trivial portal link is associated with
+  /// @return The target volume
+  const TrackingVolume& volume() const;
+
  private:
   TrackingVolume* m_volume;
 };

Core/src/Geometry/CompositePortalLink.cpp

@@ -8,21 +8,74 @@
 
 #include "Acts/Geometry/CompositePortalLink.hpp"
 
+#include "Acts/Geometry/GridPortalLink.hpp"
 #include "Acts/Geometry/PortalError.hpp"
+#include "Acts/Geometry/TrivialPortalLink.hpp"
 #include "Acts/Surfaces/CylinderSurface.hpp"
 #include "Acts/Surfaces/DiscSurface.hpp"
 #include "Acts/Surfaces/PlaneSurface.hpp"
+#include "Acts/Surfaces/RadialBounds.hpp"
 #include "Acts/Surfaces/RegularSurface.hpp"
+#include "Acts/Utilities/Axis.hpp"
 
+#include <algorithm>
 #include <iostream>
+#include <iterator>
 #include <stdexcept>
 
+#include <boost/algorithm/string/join.hpp>
+
 namespace Acts {
 
+namespace {
+std::shared_ptr<RegularSurface> mergedSurface(const Surface& a,
+                                              const Surface& b,
+                                              BinningValue direction) {
+  if (a.type() != b.type()) {
+    throw std::invalid_argument{"Cannot merge surfaces of different types"};
+  }
+
+  if (const auto* cylinderA = dynamic_cast<const CylinderSurface*>(&a);
+      cylinderA != nullptr) {
+    const auto& cylinderB = dynamic_cast<const CylinderSurface&>(b);
+
+    auto [merged, reversed] = cylinderA->mergedWith(cylinderB, direction, true);
+    return merged;
+  } else if (const auto* discA = dynamic_cast<const DiscSurface*>(&a);
+             discA != nullptr) {
+    const auto& discB = dynamic_cast<const DiscSurface&>(b);
+    auto [merged, reversed] = discA->mergedWith(discB, direction, true);
+    return merged;
+  } else if (const auto* planeA = dynamic_cast<const PlaneSurface*>(&a);
+             planeA != nullptr) {
+    throw std::logic_error{"Plane surfaces not implemented yet"};
+  } else {
+    throw std::invalid_argument{"Unsupported surface type"};
+  }
+}
+
+std::shared_ptr<RegularSurface> mergePortalLinks(
+    const std::vector<std::unique_ptr<PortalLinkBase>>& links,
+    BinningValue direction) {
+  assert(std::ranges::all_of(links,
+                             [](const auto& link) { return link != nullptr; }));
+  assert(!links.empty());
+
+  std::shared_ptr<RegularSurface> result = links.front()->surfacePtr();
+  for (auto it = std::next(links.begin()); it != links.end(); ++it) {
+    assert(result != nullptr);
+    result = mergedSurface(*result, it->get()->surface(), direction);
+  }
+
+  return result;
+}
+}  // namespace
+
 CompositePortalLink::CompositePortalLink(std::unique_ptr<PortalLinkBase> a,
                                          std::unique_ptr<PortalLinkBase> b,
                                          BinningValue direction, bool flatten)
-    : PortalLinkBase(mergedSurface(a.get(), b.get(), direction)) {
+    : PortalLinkBase(mergedSurface(a->surface(), b->surface(), direction)),
+      m_direction{direction} {
   if (!flatten) {
     m_children.push_back(std::move(a));
     m_children.push_back(std::move(b));
@@ -45,6 +98,35 @@ CompositePortalLink::CompositePortalLink(std::unique_ptr<PortalLinkBase> a,
   }
 }
 
+CompositePortalLink::CompositePortalLink(
+    std::vector<std::unique_ptr<PortalLinkBase>> links, BinningValue direction,
+    bool flatten)
+    : PortalLinkBase(mergePortalLinks(links, direction)),
+      m_direction(direction) {
+  if (!flatten) {
+    for (auto& child : links) {
+      m_children.push_back(std::move(child));
+    }
+
+  } else {
+    auto handle = [&](std::unique_ptr<PortalLinkBase> link) {
+      if (auto* composite = dynamic_cast<CompositePortalLink*>(link.get());
+          composite != nullptr) {
+        m_children.insert(
+            m_children.end(),
+            std::make_move_iterator(composite->m_children.begin()),
+            std::make_move_iterator(composite->m_children.end()));
+      } else {
+        m_children.push_back(std::move(link));
+      }
+    };
+
+    for (auto& child : links) {
+      handle(std::move(child));
+    }
+  }
+}
+
 Result<const TrackingVolume*> CompositePortalLink::resolveVolume(
     const GeometryContext& gctx, const Vector2& position,
     double tolerance) const {
@@ -74,36 +156,6 @@ Result<const TrackingVolume*> CompositePortalLink::resolveVolume(
   return PortalError::PositionNotOnAnyChildPortalLink;
 }
 
-std::shared_ptr<RegularSurface> CompositePortalLink::mergedSurface(
-    const PortalLinkBase* a, const PortalLinkBase* b, BinningValue direction) {
-  assert(a != nullptr);
-  assert(b != nullptr);
-  if (a->surface().type() != b->surface().type()) {
-    throw std::invalid_argument{"Cannot merge surfaces of different types"};
-  }
-
-  if (const auto* cylinderA =
-          dynamic_cast<const CylinderSurface*>(&a->surface());
-      cylinderA != nullptr) {
-    const auto& cylinderB = dynamic_cast<const CylinderSurface&>(b->surface());
-
-    auto [merged, reversed] = cylinderA->mergedWith(cylinderB, direction, true);
-    return merged;
-  } else if (const auto* discA =
-                 dynamic_cast<const DiscSurface*>(&a->surface());
-             discA != nullptr) {
-    const auto& discB = dynamic_cast<const DiscSurface&>(b->surface());
-    auto [merged, reversed] = discA->mergedWith(discB, direction, true);
-    return merged;
-  } else if (const auto* planeA =
-                 dynamic_cast<const PlaneSurface*>(&a->surface());
-             planeA != nullptr) {
-    throw std::logic_error{"Plane surfaces not implemented yet"};
-  } else {
-    throw std::invalid_argument{"Unsupported surface type"};
-  }
-}
-
 std::size_t CompositePortalLink::depth() const {
   std::size_t result = 1;
 
@@ -125,4 +177,122 @@ void CompositePortalLink::toStream(std::ostream& os) const {
   os << "CompositePortalLink";
 }
 
+std::unique_ptr<GridPortalLink> CompositePortalLink::makeGrid(
+    const GeometryContext& gctx, const Logger& logger) const {
+  ACTS_VERBOSE("Attempting to make a grid from a composite portal link");
+  std::vector<TrivialPortalLink*> trivialLinks;
+  std::ranges::transform(m_children, std::back_inserter(trivialLinks),
+                         [](const auto& child) {
+                           return dynamic_cast<TrivialPortalLink*>(child.get());
+                         });
+
+  if (std::ranges::any_of(trivialLinks,
+                          [](auto link) { return link == nullptr; })) {
+    ACTS_VERBOSE(
+        "Failed to make a grid from a composite portal link -> returning "
+        "nullptr");
+    return nullptr;
+  }
+
+  // we already know all children have surfaces that are compatible, we produced
+  // a merged surface for the overall dimensions.
+
+  auto printEdges = [](const auto& edges) -> std::string {
+    std::vector<std::string> strings;
+    std::ranges::transform(edges, std::back_inserter(strings),
+                           [](const auto& v) { return std::to_string(v); });
+    return boost::algorithm::join(strings, ", ");
+  };
+
+  if (surface().type() == Surface::SurfaceType::Cylinder) {
+    ACTS_VERBOSE("Combining composite into cylinder grid");
+
+    if (m_direction != BinningValue::binZ) {
+      ACTS_ERROR("Cylinder grid only supports binning in Z direction");
+      throw std::runtime_error{"Unsupported binning direction"};
+    }
+
+    std::vector<double> edges;
+    edges.reserve(m_children.size() + 1);
+
+    const Transform3& groupTransform = m_surface->transform(gctx);
+    Transform3 itransform = groupTransform.inverse();
+
+    std::ranges::sort(
+        trivialLinks, [&itransform, &gctx](const auto& a, const auto& b) {
+          return (itransform * a->surface().transform(gctx)).translation()[eZ] <
+                 (itransform * b->surface().transform(gctx)).translation()[eZ];
+        });
+
+    for (const auto& [i, child] : enumerate(trivialLinks)) {
+      const auto& bounds =
+          dynamic_cast<const CylinderBounds&>(child->surface().bounds());
+      Transform3 ltransform = itransform * child->surface().transform(gctx);
+      ActsScalar hlZ = bounds.get(CylinderBounds::eHalfLengthZ);
+      ActsScalar minZ = ltransform.translation()[eZ] - hlZ;
+      ActsScalar maxZ = ltransform.translation()[eZ] + hlZ;
+      if (i == 0) {
+        edges.push_back(minZ);
+      }
+      edges.push_back(maxZ);
+    }
+
+    ACTS_VERBOSE("~> Determined bin edges to be " << printEdges(edges));
+
+    Axis axis{AxisBound, edges};
+
+    auto grid = GridPortalLink::make(m_surface, m_direction, std::move(axis));
+    for (const auto& [i, child] : enumerate(trivialLinks)) {
+      grid->atLocalBins({i + 1}) = &child->volume();
+    }
+
+    return grid;
+
+  } else if (surface().type() == Surface::SurfaceType::Disc) {
+    ACTS_VERBOSE("Combining composite into disc grid");
+
+    if (m_direction != BinningValue::binR) {
+      ACTS_ERROR("Disc grid only supports binning in R direction");
+      throw std::runtime_error{"Unsupported binning direction"};
+    }
+
+    std::vector<double> edges;
+    edges.reserve(m_children.size() + 1);
+
+    std::ranges::sort(trivialLinks, [](const auto& a, const auto& b) {
+      const auto& boundsA =
+          dynamic_cast<const RadialBounds&>(a->surface().bounds());
+      const auto& boundsB =
+          dynamic_cast<const RadialBounds&>(b->surface().bounds());
+      return boundsA.get(RadialBounds::eMinR) <
+             boundsB.get(RadialBounds::eMinR);
+    });
+
+    for (const auto& [i, child] : enumerate(trivialLinks)) {
+      const auto& bounds =
+          dynamic_cast<const RadialBounds&>(child->surface().bounds());
+
+      if (i == 0) {
+        edges.push_back(bounds.get(RadialBounds::eMinR));
+      }
+      edges.push_back(bounds.get(RadialBounds::eMaxR));
+    }
+
+    ACTS_VERBOSE("~> Determined bin edges to be " << printEdges(edges));
+
+    Axis axis{AxisBound, edges};
+
+    auto grid = GridPortalLink::make(m_surface, m_direction, std::move(axis));
+    for (const auto& [i, child] : enumerate(trivialLinks)) {
+      grid->atLocalBins({i + 1}) = &child->volume();
+    }
+
+    return grid;
+  } else if (surface().type() == Surface::SurfaceType::Plane) {
+    throw std::runtime_error{"Plane surfaces not implemented yet"};
+  } else {
+    throw std::invalid_argument{"Unsupported surface type"};
+  }
+}
+
 }  // namespace Acts