refactor(geometry): Surface merging returns ordering

Core/include/Acts/Surfaces/CylinderSurface.hpp

@@ -252,14 +252,14 @@ class CylinderSurface : public RegularSurface {
   /// @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 externalRotation If true, any phi rotation is done in the transform
   /// @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;
+  /// @return The merged cylinder surface and the ordering of input surfaces
+  std::pair<std::shared_ptr<CylinderSurface>, bool> mergedWith(
+      const CylinderSurface& other, BinningValue direction,
+      bool externalRotation, const Logger& logger = getDummyLogger()) const;
 
  protected:
   std::shared_ptr<const CylinderBounds> m_bounds;  //!< bounds (shared)

Core/include/Acts/Surfaces/DiscSurface.hpp

@@ -335,14 +335,14 @@ class DiscSurface : public RegularSurface {
   /// @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 externalRotation If true, any phi rotation is done in the transform
   /// @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;
+  /// @return The merged disc surface and the ordering of input surfaces
+  std::pair<std::shared_ptr<DiscSurface>, bool> mergedWith(
+      const DiscSurface& other, BinningValue direction, bool externalRotation,
+      const Logger& logger = getDummyLogger()) const;
 
  protected:
   std::shared_ptr<const DiscBounds> m_bounds;  ///< bounds (shared)

Core/include/Acts/Surfaces/Surface.hpp

@@ -528,4 +528,6 @@ class Surface : public virtual GeometryObject,
       const Vector3& direction) const;
 };
 
+std::ostream& operator<<(std::ostream& os, const Surface& srf);
+
 }  // namespace Acts

Core/include/Acts/Surfaces/detail/MergeHelper.hpp

@@ -23,7 +23,7 @@ namespace Acts::detail {
 /// 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(
+std::tuple<ActsScalar, ActsScalar, bool> mergedPhiSector(
     ActsScalar hlPhi1, ActsScalar avgPhi1, ActsScalar hlPhi2,
     ActsScalar avgPhi2, const Logger& logger = getDummyLogger(),
     ActsScalar tolerance = s_onSurfaceTolerance);

Core/src/Surfaces/CylinderSurface.cpp

@@ -368,20 +368,32 @@ 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 {
+std::pair<std::shared_ptr<Acts::CylinderSurface>, bool>
+Acts::CylinderSurface::mergedWith(const CylinderSurface& other,
+                                  BinningValue direction, bool externalRotation,
+                                  const Logger& logger) const {
   using namespace Acts::UnitLiterals;
 
   ACTS_DEBUG("Merging cylinder surfaces in " << binningValueName(direction)
                                              << " direction");
 
-  Transform3 otherLocal = transform(gctx).inverse() * other.transform(gctx);
+  if (m_associatedDetElement != nullptr ||
+      other.m_associatedDetElement != nullptr) {
+    throw SurfaceMergingException(getSharedPtr(), other.getSharedPtr(),
+                                  "CylinderSurface::merge: surfaces are "
+                                  "associated with a detector element");
+  }
+
+  assert(m_transform != nullptr && other.m_transform != nullptr);
+
+  Transform3 otherLocal = m_transform->inverse() * *other.m_transform;
 
   constexpr auto tolerance = s_onSurfaceTolerance;
 
   // surface cannot have any relative rotation
-  if (!otherLocal.linear().isApprox(RotationMatrix3::Identity())) {
+
+  if (std::abs(otherLocal.linear().col(eX)[eZ]) >= tolerance ||
+      std::abs(otherLocal.linear().col(eY)[eZ]) >= tolerance) {
     ACTS_ERROR("CylinderSurface::merge: surfaces have relative rotation");
     throw SurfaceMergingException(
         getSharedPtr(), other.getSharedPtr(),
@@ -436,17 +448,31 @@ std::shared_ptr<Acts::CylinderSurface> Acts::CylinderSurface::mergedWith(
         "CylinderSurface::merge: surfaces have relative translation in x/y");
   }
 
+  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;
+
+  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);
+
   if (direction == Acts::BinningValue::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;
+    if (std::abs(otherLocal.linear().col(eY)[eX]) >= tolerance &&
+        (!bounds().coversFullAzimuth() ||
+         !other.bounds().coversFullAzimuth())) {
+      throw SurfaceMergingException(getSharedPtr(), other.getSharedPtr(),
+                                    "CylinderSurface::merge: surfaces have "
+                                    "relative rotation in z and phi sector");
+    }
 
     ACTS_VERBOSE("this: [" << minZ << ", " << maxZ << "] ~> "
                            << (minZ + maxZ) / 2.0 << " +- " << hlZ);
@@ -455,7 +481,6 @@ std::shared_ptr<Acts::CylinderSurface> Acts::CylinderSurface::mergedWith(
     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");
@@ -464,19 +489,26 @@ std::shared_ptr<Acts::CylinderSurface> Acts::CylinderSurface::mergedWith(
           "CylinderSurface::merge: surfaces have incompatible z bounds");
     }
 
+    if (hlPhi != otherHlPhi || avgPhi != otherAvgPhi) {
+      throw SurfaceMergingException(getSharedPtr(), other.getSharedPtr(),
+                                    "CylinderSurface::merge: surfaces have "
+                                    "different phi sectors");
+    }
+
     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);
+    auto newBounds = std::make_shared<CylinderBounds>(r, newHlZ, hlPhi, avgPhi);
 
     Transform3 newTransform =
-        transform(gctx) * Translation3{Vector3::UnitZ() * newMidZ};
+        *m_transform * Translation3{Vector3::UnitZ() * newMidZ};
 
-    return Surface::makeShared<CylinderSurface>(newTransform, newBounds);
+    return {Surface::makeShared<CylinderSurface>(newTransform, newBounds),
+            zShift < 0};
 
   } else if (direction == Acts::BinningValue::binRPhi) {
     // no z shift is allowed
@@ -490,20 +522,47 @@ std::shared_ptr<Acts::CylinderSurface> Acts::CylinderSurface::mergedWith(
           "rPhi merging");
     }
 
-    ActsScalar hlPhi = bounds().get(CylinderBounds::eHalfPhiSector);
-    ActsScalar avgPhi = bounds().get(CylinderBounds::eAveragePhi);
+    if (hlZ != otherHlZ) {
+      throw SurfaceMergingException(getSharedPtr(), other.getSharedPtr(),
+                                    "CylinderSurface::merge: surfaces have "
+                                    "different z bounds");
+    }
 
-    ActsScalar otherHlPhi = other.bounds().get(CylinderBounds::eHalfPhiSector);
-    ActsScalar otherAvgPhi = other.bounds().get(CylinderBounds::eAveragePhi);
+    // Figure out signed relative rotation
+    Vector2 rotatedX = otherLocal.linear().col(eX).head<2>();
+    ActsScalar zrotation = std::atan2(rotatedX[1], rotatedX[0]);
+
+    ACTS_VERBOSE("this:  [" << avgPhi / 1_degree << " +- " << hlPhi / 1_degree
+                            << "]");
+    ACTS_VERBOSE("other: [" << otherAvgPhi / 1_degree << " +- "
+                            << otherHlPhi / 1_degree << "]");
+
+    ACTS_VERBOSE("Relative rotation around local z: " << zrotation / 1_degree);
+
+    ActsScalar prevOtherAvgPhi = otherAvgPhi;
+    otherAvgPhi = detail::radian_sym(otherAvgPhi + zrotation);
+    ACTS_VERBOSE("~> local other average phi: "
+                 << otherAvgPhi / 1_degree
+                 << " (was: " << prevOtherAvgPhi / 1_degree << ")");
 
     try {
-      auto [newHlPhi, newAvgPhi] = detail::mergedPhiSector(
+      auto [newHlPhi, newAvgPhi, reversed] = detail::mergedPhiSector(
           hlPhi, avgPhi, otherHlPhi, otherAvgPhi, logger, tolerance);
 
+      Transform3 newTransform = *m_transform;
+
+      if (externalRotation) {
+        ACTS_VERBOSE("Modifying transform for external rotation of "
+                     << newAvgPhi / 1_degree);
+        newTransform = newTransform * AngleAxis3(newAvgPhi, Vector3::UnitZ());
+        newAvgPhi = 0.;
+      }
+
       auto newBounds = std::make_shared<CylinderBounds>(
           r, bounds().get(CylinderBounds::eHalfLengthZ), newHlPhi, newAvgPhi);
 
-      return Surface::makeShared<CylinderSurface>(transform(gctx), newBounds);
+      return {Surface::makeShared<CylinderSurface>(newTransform, newBounds),
+              reversed};
     } catch (const std::invalid_argument& e) {
       throw SurfaceMergingException(getSharedPtr(), other.getSharedPtr(),
                                     e.what());

Core/src/Surfaces/DiscSurface.cpp

@@ -379,20 +379,29 @@ double Acts::DiscSurface::pathCorrection(const GeometryContext& gctx,
   return 1. / std::abs(normal(gctx).dot(direction));
 }
 
-std::shared_ptr<Acts::DiscSurface> Acts::DiscSurface::mergedWith(
-    const GeometryContext& gctx, const DiscSurface& other,
-    BinningValue direction, const Logger& logger) const {
+std::pair<std::shared_ptr<Acts::DiscSurface>, bool>
+Acts::DiscSurface::mergedWith(const DiscSurface& other, BinningValue direction,
+                              bool externalRotation,
+                              const Logger& logger) const {
   using namespace Acts::UnitLiterals;
 
-  ACTS_DEBUG("Merging disc surfaces in " << binningValueName(direction)
-                                         << " direction");
+  ACTS_DEBUG("Merging disc surfaces in " << direction << " direction");
 
-  Transform3 otherLocal = transform(gctx).inverse() * other.transform(gctx);
+  if (m_associatedDetElement != nullptr ||
+      other.m_associatedDetElement != nullptr) {
+    throw SurfaceMergingException(getSharedPtr(), other.getSharedPtr(),
+                                  "CylinderSurface::merge: surfaces are "
+                                  "associated with a detector element");
+  }
+  assert(m_transform != nullptr && other.m_transform != nullptr);
+
+  Transform3 otherLocal = m_transform->inverse() * *other.m_transform;
 
   constexpr auto tolerance = s_onSurfaceTolerance;
 
   // surface cannot have any relative rotation
-  if (!otherLocal.linear().isApprox(RotationMatrix3::Identity())) {
+  if (std::abs(otherLocal.linear().col(eX)[eZ]) >= tolerance ||
+      std::abs(otherLocal.linear().col(eY)[eZ]) >= tolerance) {
     ACTS_ERROR("DiscSurface::merge: surfaces have relative rotation");
     throw SurfaceMergingException(
         getSharedPtr(), other.getSharedPtr(),
@@ -449,6 +458,13 @@ std::shared_ptr<Acts::DiscSurface> Acts::DiscSurface::mergedWith(
                                 << otherHlPhi / 1_degree);
 
   if (direction == Acts::BinningValue::binR) {
+    if (std::abs(otherLocal.linear().col(eY)[eX]) >= tolerance &&
+        (!bounds->coversFullAzimuth() || !otherBounds->coversFullAzimuth())) {
+      throw SurfaceMergingException(getSharedPtr(), other.getSharedPtr(),
+                                    "DiscSurface::merge: surfaces have "
+                                    "relative rotation in z and phi sector");
+    }
+
     if (std::abs(minR - otherMaxR) > tolerance &&
         std::abs(maxR - otherMinR) > tolerance) {
       ACTS_ERROR("DiscSurface::merge: surfaces are not touching r");
@@ -478,7 +494,8 @@ std::shared_ptr<Acts::DiscSurface> Acts::DiscSurface::mergedWith(
     auto newBounds =
         std::make_shared<RadialBounds>(newMinR, newMaxR, hlPhi, avgPhi);
 
-    return Surface::makeShared<DiscSurface>(transform(gctx), newBounds);
+    return {Surface::makeShared<DiscSurface>(*m_transform, newBounds),
+            minR > otherMinR};
 
   } else if (direction == Acts::BinningValue::binPhi) {
     if (std::abs(maxR - otherMaxR) > tolerance ||
@@ -489,22 +506,49 @@ std::shared_ptr<Acts::DiscSurface> Acts::DiscSurface::mergedWith(
           "DiscSurface::merge: surfaces don't have same r bounds");
     }
 
+    // Figure out signed relative rotation
+    Vector2 rotatedX = otherLocal.linear().col(eX).head<2>();
+    ActsScalar zrotation = std::atan2(rotatedX[1], rotatedX[0]);
+
+    ACTS_VERBOSE("this:  [" << avgPhi / 1_degree << " +- " << hlPhi / 1_degree
+                            << "]");
+    ACTS_VERBOSE("other: [" << otherAvgPhi / 1_degree << " +- "
+                            << otherHlPhi / 1_degree << "]");
+
+    ACTS_VERBOSE("Relative rotation around local z: " << zrotation / 1_degree);
+
+    ActsScalar prevOtherAvgPhi = otherAvgPhi;
+    otherAvgPhi = detail::radian_sym(otherAvgPhi + zrotation);
+    ACTS_VERBOSE("~> local other average phi: "
+                 << otherAvgPhi / 1_degree
+                 << " (was: " << prevOtherAvgPhi / 1_degree << ")");
+
     try {
-      auto [newHlPhi, newAvgPhi] = detail::mergedPhiSector(
+      auto [newHlPhi, newAvgPhi, reversed] = detail::mergedPhiSector(
           hlPhi, avgPhi, otherHlPhi, otherAvgPhi, logger, tolerance);
 
+      Transform3 newTransform = *m_transform;
+
+      if (externalRotation) {
+        ACTS_VERBOSE("Modifying transform for external rotation of "
+                     << newAvgPhi / 1_degree);
+        newTransform = newTransform * AngleAxis3(newAvgPhi, Vector3::UnitZ());
+        newAvgPhi = 0.;
+      }
+
       auto newBounds =
           std::make_shared<RadialBounds>(minR, maxR, newHlPhi, newAvgPhi);
 
-      return Surface::makeShared<DiscSurface>(transform(gctx), newBounds);
+      return {Surface::makeShared<DiscSurface>(newTransform, newBounds),
+              reversed};
     } catch (const std::invalid_argument& e) {
       throw SurfaceMergingException(getSharedPtr(), other.getSharedPtr(),
                                     e.what());
     }
 
   } else {
-    ACTS_ERROR("DiscSurface::merge: invalid direction "
-               << binningValueName(direction));
+    ACTS_ERROR("DiscSurface::merge: invalid direction " << direction);
+
     throw SurfaceMergingException(
         getSharedPtr(), other.getSharedPtr(),
         "DiscSurface::merge: invalid direction " + binningValueName(direction));