feat: Cylinder and Disc surface merging (acts-project#3288)

This PR allows merging of two cylinder and disc surfaces. Merging is allowed in the $z$ and $r\times\phi$ direction for cylinders, and in the $r$ and $\phi$ direction for discs. The surfaces need to be consistent, i.e. their bounds dimensions and transforms need to line up. - For **cylinders**: - $z$ merging: $z$ edges must be aligned, no relative $xy$ shift or any rotation is allowed. - $r\times\phi$ merging: $r\times\phi$ edges must be aligned, no $xyz$ shift or any rotation is allowed - For **discs**: - $r$ merging: $r$ edges must be aligned, no $xyz$ shift or rotation is allowed - $\phi$ merging: $\phi$ edges must be aligned, no $xyz$ shift or rotation is allowed. ![Cylinder_Merging](https://github.com/acts-project/acts/assets/1058585/ed7b2a9d-764b-473a-9c5a-675651232786) ![Disc_Merging](https://github.com/acts-project/acts/assets/1058585/3fe97e7e-94ee-4d1a-8277-8cd01104b06a) Blocked by: - acts-project#3290
Matthewharri · Jun 18, 2024 · 8e6eba4 · 8e6eba4
1 parent 7cf4883
commit 8e6eba4
Show file tree

Hide file tree

Showing 14 changed files with 1,007 additions and 28 deletions.
diff --git a/Core/include/Acts/Surfaces/CylinderBounds.hpp b/Core/include/Acts/Surfaces/CylinderBounds.hpp
@@ -159,25 +159,4 @@ inline bool CylinderBounds::coversFullAzimuth() const {
   return m_closed;
 }
 
-inline void CylinderBounds::checkConsistency() noexcept(false) {
-  if (get(eR) <= 0.) {
-    throw std::invalid_argument("CylinderBounds: invalid radial setup.");
-  }
-  if (get(eHalfLengthZ) <= 0.) {
-    throw std::invalid_argument("CylinderBounds: invalid length setup.");
-  }
-  if (get(eHalfPhiSector) <= 0. || get(eHalfPhiSector) > M_PI) {
-    throw std::invalid_argument("CylinderBounds: invalid phi sector setup.");
-  }
-  if (get(eAveragePhi) != detail::radian_sym(get(eAveragePhi))) {
-    throw std::invalid_argument("CylinderBounds: invalid phi positioning.");
-  }
-  if (get(eBevelMinZ) != detail::radian_sym(get(eBevelMinZ))) {
-    throw std::invalid_argument("CylinderBounds: invalid bevel at min Z.");
-  }
-  if (get(eBevelMaxZ) != detail::radian_sym(get(eBevelMaxZ))) {
-    throw std::invalid_argument("CylinderBounds: invalid bevel at max Z.");
-  }
-}
-
 }  // namespace Acts
diff --git a/Core/include/Acts/Surfaces/CylinderSurface.hpp b/Core/include/Acts/Surfaces/CylinderSurface.hpp
@@ -21,6 +21,7 @@
 #include "Acts/Surfaces/SurfaceConcept.hpp"
 #include "Acts/Utilities/BinningType.hpp"
 #include "Acts/Utilities/Concepts.hpp"
+#include "Acts/Utilities/Logger.hpp"
 #include "Acts/Utilities/Result.hpp"
 #include "Acts/Utilities/detail/RealQuadraticEquation.hpp"
 
@@ -247,6 +248,19 @@ class CylinderSurface : public RegularSurface {
   ActsMatrix<2, 3> localCartesianToBoundLocalDerivative(
       const GeometryContext& gctx, const Vector3& position) const final;
 
+  /// Merge two cylinder surfaces into a single one.
+  /// @image html Cylinder_Merging.svg
+  /// @note The surfaces need to be *compatible*, i.e. have cylinder bounds
+  ///       that align, and have the same radius
+  /// @param gctx The current geometry context object, e.g. alignment
+  /// @param other The other cylinder surface to merge with
+  /// @param direction The binning direction: either @c binZ or @c binRPhi
+  /// @param logger The logger to use
+  /// @return The merged cylinder surface
+  std::shared_ptr<CylinderSurface> mergedWith(
+      const GeometryContext& gctx, const CylinderSurface& other,
+      BinningValue direction, const Logger& logger = getDummyLogger()) const;
+
  protected:
   std::shared_ptr<const CylinderBounds> m_bounds;  //!< bounds (shared)
 

diff --git a/Core/include/Acts/Surfaces/DiscSurface.hpp b/Core/include/Acts/Surfaces/DiscSurface.hpp
@@ -331,6 +331,19 @@ class DiscSurface : public RegularSurface {
   ActsMatrix<2, 3> localCartesianToBoundLocalDerivative(
       const GeometryContext& gctx, const Vector3& position) const final;
 
+  /// Merge two disc surfaces into a single one.
+  /// @image html Disc_Merging.svg
+  /// @note The surfaces need to be *compatible*, i.e. have disc bounds
+  ///       that align
+  /// @param gctx The current geometry context object, e.g. alignment
+  /// @param other The other disc surface to merge with
+  /// @param direction The binning direction: either @c binR or @c binPhi
+  /// @param logger The logger to use
+  /// @return The merged disc surface
+  std::shared_ptr<DiscSurface> mergedWith(
+      const GeometryContext& gctx, const DiscSurface& other,
+      BinningValue direction, const Logger& logger = getDummyLogger()) const;
+
  protected:
   std::shared_ptr<const DiscBounds> m_bounds;  ///< bounds (shared)
 };

diff --git a/Core/include/Acts/Surfaces/detail/MergeHelper.hpp b/Core/include/Acts/Surfaces/detail/MergeHelper.hpp
@@ -0,0 +1,31 @@
+// This file is part of the Acts project.
+//
+// Copyright (C) 2024 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 http://mozilla.org/MPL/2.0/.
+
+#pragma once
+
+#include "Acts/Definitions/Algebra.hpp"
+#include "Acts/Definitions/Tolerance.hpp"
+#include "Acts/Definitions/Units.hpp"
+#include "Acts/Utilities/Logger.hpp"
+#include "Acts/Utilities/detail/periodic.hpp"
+
+#include <utility>
+
+namespace Acts::detail {
+
+/// This helper function calculates the combination
+/// of two phi sectors, defined by a phi half-length +
+/// a half phi sector in the range [0,pi). The two
+/// ranges need to line up, i.e. that one of the sector
+/// ends exactly where the other one starts.
+std::pair<ActsScalar, ActsScalar> mergedPhiSector(
+    ActsScalar hlPhi1, ActsScalar avgPhi1, ActsScalar hlPhi2,
+    ActsScalar avgPhi2, const Logger& logger = getDummyLogger(),
+    ActsScalar tolerance = s_onSurfaceTolerance);
+
+}  // namespace Acts::detail
diff --git a/Core/src/Surfaces/CMakeLists.txt b/Core/src/Surfaces/CMakeLists.txt
@@ -28,4 +28,5 @@ target_sources(
     CurvilinearSurface.cpp
     detail/AlignmentHelper.cpp
     detail/AnnulusBoundsHelper.cpp
+    detail/MergeHelper.cpp
 )
diff --git a/Core/src/Surfaces/CylinderBounds.cpp b/Core/src/Surfaces/CylinderBounds.cpp
@@ -191,3 +191,25 @@ std::vector<Acts::Vector3> Acts::CylinderBounds::createCircles(
   }
   return vertices;
 }
+
+void Acts::CylinderBounds::checkConsistency() noexcept(false) {
+  if (get(eR) <= 0.) {
+    throw std::invalid_argument("CylinderBounds: invalid radial setup.");
+  }
+  if (get(eHalfLengthZ) <= 0.) {
+    throw std::invalid_argument("CylinderBounds: invalid length setup.");
+  }
+  if (get(eHalfPhiSector) <= 0. || get(eHalfPhiSector) > M_PI) {
+    throw std::invalid_argument("CylinderBounds: invalid phi sector setup.");
+  }
+  if (get(eAveragePhi) != detail::radian_sym(get(eAveragePhi)) &&
+      std::abs(std::abs(get(eAveragePhi)) - M_PI) > s_epsilon) {
+    throw std::invalid_argument("CylinderBounds: invalid phi positioning.");
+  }
+  if (get(eBevelMinZ) != detail::radian_sym(get(eBevelMinZ))) {
+    throw std::invalid_argument("CylinderBounds: invalid bevel at min Z.");
+  }
+  if (get(eBevelMaxZ) != detail::radian_sym(get(eBevelMaxZ))) {
+    throw std::invalid_argument("CylinderBounds: invalid bevel at max Z.");
+  }
+}
diff --git a/Core/src/Surfaces/CylinderSurface.cpp b/Core/src/Surfaces/CylinderSurface.cpp
@@ -8,17 +8,26 @@
 
 #include "Acts/Surfaces/CylinderSurface.hpp"
 
+#include "Acts/Definitions/Algebra.hpp"
+#include "Acts/Definitions/Tolerance.hpp"
+#include "Acts/Definitions/Units.hpp"
 #include "Acts/Geometry/GeometryObject.hpp"
+#include "Acts/Surfaces/CylinderBounds.hpp"
 #include "Acts/Surfaces/SurfaceError.hpp"
 #include "Acts/Surfaces/detail/AlignmentHelper.hpp"
 #include "Acts/Surfaces/detail/FacesHelper.hpp"
+#include "Acts/Surfaces/detail/MergeHelper.hpp"
+#include "Acts/Utilities/BinningType.hpp"
 #include "Acts/Utilities/Helpers.hpp"
 #include "Acts/Utilities/Intersection.hpp"
 #include "Acts/Utilities/ThrowAssert.hpp"
+#include "Acts/Utilities/detail/periodic.hpp"
 
 #include <algorithm>
 #include <cassert>
 #include <cmath>
+#include <iostream>
+#include <memory>
 #include <utility>
 #include <vector>
 
@@ -353,3 +362,137 @@ Acts::CylinderSurface::localCartesianToBoundLocalDerivative(
 
   return loc3DToLocBound;
 }
+
+std::shared_ptr<Acts::CylinderSurface> Acts::CylinderSurface::mergedWith(
+    const GeometryContext& gctx, const CylinderSurface& other,
+    BinningValue direction, const Logger& logger) const {
+  using namespace Acts::UnitLiterals;
+
+  ACTS_DEBUG("Merging cylinder surfaces in " << binningValueNames()[direction]
+                                             << " direction")
+
+  Transform3 otherLocal = transform(gctx).inverse() * other.transform(gctx);
+
+  constexpr auto tolerance = s_onSurfaceTolerance;
+
+  // surface cannot have any relative rotation
+  if (!otherLocal.linear().isApprox(RotationMatrix3::Identity())) {
+    ACTS_ERROR("CylinderSurface::merge: surfaces have relative rotation");
+    throw std::invalid_argument(
+        "CylinderSurface::merge: surfaces have relative rotation");
+  }
+
+  auto checkNoBevel = [&logger](const auto& bounds) {
+    if (bounds.get(CylinderBounds::eBevelMinZ) != 0.0) {
+      ACTS_ERROR(
+          "CylinderVolumeStack requires all volumes to have a bevel angle of "
+          "0");
+      throw std::invalid_argument(
+          "CylinderVolumeStack requires all volumes to have a bevel angle of "
+          "0");
+    }
+
+    if (bounds.get(CylinderBounds::eBevelMaxZ) != 0.0) {
+      ACTS_ERROR(
+          "CylinderVolumeStack requires all volumes to have a bevel angle of "
+          "0");
+      throw std::invalid_argument(
+          "CylinderVolumeStack requires all volumes to have a bevel angle of "
+          "0");
+    }
+  };
+
+  checkNoBevel(bounds());
+  checkNoBevel(other.bounds());
+
+  // radii need to be identical
+  if (std::abs(bounds().get(CylinderBounds::eR) -
+               other.bounds().get(CylinderBounds::eR)) > tolerance) {
+    ACTS_ERROR("CylinderSurface::merge: surfaces have different radii");
+    throw std::invalid_argument(
+        "CylinderSurface::merge: surfaces have different radii");
+  }
+
+  ActsScalar r = bounds().get(CylinderBounds::eR);
+
+  // no translation in x/z is allowed
+  Vector3 translation = otherLocal.translation();
+
+  if (std::abs(translation[0]) > tolerance ||
+      std::abs(translation[1]) > tolerance) {
+    ACTS_ERROR(
+        "CylinderSurface::merge: surfaces have relative translation in x/y");
+    throw std::invalid_argument(
+        "CylinderSurface::merge: surfaces have relative translation in x/y");
+  }
+
+  if (direction == Acts::binZ) {
+    // z shift must match the bounds
+
+    ActsScalar hlZ = bounds().get(CylinderBounds::eHalfLengthZ);
+    ActsScalar minZ = -hlZ;
+    ActsScalar maxZ = hlZ;
+
+    ActsScalar zShift = translation[2];
+    ActsScalar otherHlZ = other.bounds().get(CylinderBounds::eHalfLengthZ);
+    ActsScalar otherMinZ = -otherHlZ + zShift;
+    ActsScalar otherMaxZ = otherHlZ + zShift;
+
+    ACTS_VERBOSE("this: [" << minZ << ", " << maxZ << "] ~> "
+                           << (minZ + maxZ) / 2.0 << " +- " << hlZ);
+    ACTS_VERBOSE("zShift: " << zShift);
+
+    ACTS_VERBOSE("other: [" << otherMinZ << ", " << otherMaxZ << "] ~> "
+                            << (otherMinZ + otherMaxZ) / 2.0 << " +- "
+                            << otherHlZ);
+
+    if (std::abs(maxZ - otherMinZ) > tolerance &&
+        std::abs(minZ - otherMaxZ) > tolerance) {
+      ACTS_ERROR("CylinderSurface::merge: surfaces have incompatible z bounds");
+      throw std::invalid_argument(
+          "CylinderSurface::merge: surfaces have incompatible z bounds");
+    }
+
+    ActsScalar newMaxZ = std::max(maxZ, otherMaxZ);
+    ActsScalar newMinZ = std::min(minZ, otherMinZ);
+    ActsScalar newHlZ = (newMaxZ - newMinZ) / 2.0;
+    ActsScalar newMidZ = (newMaxZ + newMinZ) / 2.0;
+    ACTS_VERBOSE("merged: [" << newMinZ << ", " << newMaxZ << "] ~> " << newMidZ
+                             << " +- " << newHlZ);
+
+    auto newBounds = std::make_shared<CylinderBounds>(r, newHlZ);
+
+    Transform3 newTransform =
+        transform(gctx) * Translation3{Vector3::UnitZ() * newMidZ};
+
+    return Surface::makeShared<CylinderSurface>(newTransform, newBounds);
+
+  } else if (direction == Acts::binRPhi) {
+    // no z shift is allowed
+    if (std::abs(translation[2]) > tolerance) {
+      ACTS_ERROR(
+          "CylinderSurface::merge: surfaces have relative translation in z for "
+          "rPhi merging");
+      throw std::invalid_argument(
+          "CylinderSurface::merge: surfaces have relative translation in z for "
+          "rPhi merging");
+    }
+
+    ActsScalar hlPhi = bounds().get(CylinderBounds::eHalfPhiSector);
+    ActsScalar avgPhi = bounds().get(CylinderBounds::eAveragePhi);
+
+    ActsScalar otherHlPhi = other.bounds().get(CylinderBounds::eHalfPhiSector);
+    ActsScalar otherAvgPhi = other.bounds().get(CylinderBounds::eAveragePhi);
+
+    auto [newHlPhi, newAvgPhi] = detail::mergedPhiSector(
+        hlPhi, avgPhi, otherHlPhi, otherAvgPhi, logger, tolerance);
+
+    auto newBounds = std::make_shared<CylinderBounds>(
+        r, bounds().get(CylinderBounds::eHalfLengthZ), newHlPhi, newAvgPhi);
+
+    return Surface::makeShared<CylinderSurface>(transform(gctx), newBounds);
+  } else {
+    throw std::invalid_argument("CylinderSurface::merge: invalid direction " +
+                                binningValueNames()[direction]);
+  }
+}