src/Surfaces/DiscSurface.cpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 #include "Acts/Surfaces/DiscSurface.hpp"
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Definitions/Units.hpp"
0013 #include "Acts/Geometry/GeometryObject.hpp"
0014 #include "Acts/Surfaces/BoundaryTolerance.hpp"
0015 #include "Acts/Surfaces/DiscBounds.hpp"
0016 #include "Acts/Surfaces/DiscTrapezoidBounds.hpp"
0017 #include "Acts/Surfaces/InfiniteBounds.hpp"
0018 #include "Acts/Surfaces/RadialBounds.hpp"
0019 #include "Acts/Surfaces/SurfaceBounds.hpp"
0020 #include "Acts/Surfaces/SurfaceError.hpp"
0021 #include "Acts/Surfaces/SurfaceMergingException.hpp"
0022 #include "Acts/Surfaces/detail/FacesHelper.hpp"
0023 #include "Acts/Surfaces/detail/MergeHelper.hpp"
0024 #include "Acts/Surfaces/detail/PlanarHelper.hpp"
0025 #include "Acts/Utilities/BinningType.hpp"
0026 #include "Acts/Utilities/Intersection.hpp"
0027 #include "Acts/Utilities/JacobianHelpers.hpp"
0028 #include "Acts/Utilities/ThrowAssert.hpp"
0029 #include "Acts/Utilities/detail/periodic.hpp"
0030
0031 #include <cmath>
0032 #include <stdexcept>
0033 #include <utility>
0034 #include <vector>
0035
0036 namespace Acts {
0037
0038 using VectorHelpers::perp;
0039 using VectorHelpers::phi;
0040
0041 DiscSurface::DiscSurface(const DiscSurface& other)
0042     : GeometryObject(), RegularSurface(other), m_bounds(other.m_bounds) {}
0043
0044 DiscSurface::DiscSurface(const GeometryContext& gctx, const DiscSurface& other,
0045                          const Transform3& shift)
0046     : GeometryObject(),
0047       RegularSurface(gctx, other, shift),
0048       m_bounds(other.m_bounds) {}
0049
0050 DiscSurface::DiscSurface(const Transform3& transform, double rmin, double rmax,
0051                          double hphisec)
0052     : GeometryObject(),
0053       RegularSurface(transform),
0054       m_bounds(std::make_shared<const RadialBounds>(rmin, rmax, hphisec)) {}
0055
0056 DiscSurface::DiscSurface(const Transform3& transform, double minhalfx,
0057                          double maxhalfx, double minR, double maxR,
0058                          double avephi, double stereo)
0059     : GeometryObject(),
0060       RegularSurface(transform),
0061       m_bounds(std::make_shared<const DiscTrapezoidBounds>(
0062           minhalfx, maxhalfx, minR, maxR, avephi, stereo)) {}
0063
0064 DiscSurface::DiscSurface(const Transform3& transform,
0065                          std::shared_ptr<const DiscBounds> dbounds)
0066     : GeometryObject(),
0067       RegularSurface(transform),
0068       m_bounds(std::move(dbounds)) {}
0069
0070 DiscSurface::DiscSurface(std::shared_ptr<const DiscBounds> dbounds,
0071                          const DetectorElementBase& detelement)
0072     : GeometryObject(),
0073       RegularSurface(detelement),
0074       m_bounds(std::move(dbounds)) {
0075   throw_assert(m_bounds, "nullptr as DiscBounds");
0076 }
0077
0078 DiscSurface& DiscSurface::operator=(const DiscSurface& other) {
0079   if (this != &other) {
0080     Surface::operator=(other);
0081     m_bounds = other.m_bounds;
0082   }
0083   return *this;
0084 }
0085
0086 Surface::SurfaceType DiscSurface::type() const {
0087   return Surface::Disc;
0088 }
0089
0090 Vector3 DiscSurface::localToGlobal(const GeometryContext& gctx,
0091                                    const Vector2& lposition) const {
0092   // create the position in the local 3d frame
0093   Vector3 loc3Dframe(lposition[0] * cos(lposition[1]),
0094                      lposition[0] * sin(lposition[1]), 0.);
0095   // transform to globalframe
0096   return transform(gctx) * loc3Dframe;
0097 }
0098
0099 Result<Vector2> DiscSurface::globalToLocal(const GeometryContext& gctx,
0100                                            const Vector3& position,
0101                                            double tolerance) const {
0102   // transport it to the globalframe
0103   Vector3 loc3Dframe = (transform(gctx).inverse()) * position;
0104   if (std::abs(loc3Dframe.z()) > std::abs(tolerance)) {
0105     return Result<Vector2>::failure(SurfaceError::GlobalPositionNotOnSurface);
0106   }
0107   return Result<Vector2>::success({perp(loc3Dframe), phi(loc3Dframe)});
0108 }
0109
0110 Vector2 DiscSurface::localPolarToLocalCartesian(const Vector2& locpol) const {
0111   const DiscTrapezoidBounds* dtbo =
0112       dynamic_cast<const DiscTrapezoidBounds*>(&(bounds()));
0113   if (dtbo != nullptr) {
0114     double rMedium = dtbo->rCenter();
0115     double phi = dtbo->get(DiscTrapezoidBounds::eAveragePhi);
0116
0117     Vector2 polarCenter(rMedium, phi);
0118     Vector2 cartCenter = localPolarToCartesian(polarCenter);
0119     Vector2 cartPos = localPolarToCartesian(locpol);
0120     Vector2 pos = cartPos - cartCenter;
0121
0122     Vector2 locPos(pos[0] * sin(phi) - pos[1] * cos(phi),
0123                    pos[1] * sin(phi) + pos[0] * cos(phi));
0124     return Vector2(locPos[0], locPos[1]);
0125   }
0126   return Vector2(locpol[0] * cos(locpol[1]), locpol[0] * sin(locpol[1]));
0127 }
0128
0129 Vector3 DiscSurface::localCartesianToGlobal(const GeometryContext& gctx,
0130                                             const Vector2& lposition) const {
0131   Vector3 loc3Dframe(lposition[0], lposition[1], 0.);
0132   return transform(gctx) * loc3Dframe;
0133 }
0134
0135 Vector2 DiscSurface::globalToLocalCartesian(const GeometryContext& gctx,
0136                                             const Vector3& position,
0137                                             double /*direction*/) const {
0138   Vector3 loc3Dframe = (transform(gctx).inverse()) * position;
0139   return Vector2(loc3Dframe.x(), loc3Dframe.y());
0140 }
0141
0142 std::string DiscSurface::name() const {
0143   return "Acts::DiscSurface";
0144 }
0145
0146 const SurfaceBounds& DiscSurface::bounds() const {
0147   if (m_bounds) {
0148     return (*(m_bounds.get()));
0149   }
0150   return s_noBounds;
0151 }
0152
0153 Polyhedron DiscSurface::polyhedronRepresentation(
0154     const GeometryContext& gctx, unsigned int quarterSegments) const {
0155   // Prepare vertices and faces
0156   std::vector<Vector3> vertices;
0157   // Understand the disc
0158   bool fullDisc = m_bounds->coversFullAzimuth();
0159   bool toCenter = m_bounds->rMin() < s_onSurfaceTolerance;
0160   // If you have bounds you can create a polyhedron representation
0161   bool exactPolyhedron = (m_bounds->type() == SurfaceBounds::eDiscTrapezoid);
0162   bool addCentreFromConvexFace = (m_bounds->type() != SurfaceBounds::eAnnulus);
0163   if (m_bounds) {
0164     auto vertices2D = m_bounds->vertices(quarterSegments);
0165     vertices.reserve(vertices2D.size() + 1);
0166     Vector3 wCenter(0., 0., 0);
0167     for (const auto& v2D : vertices2D) {
0168       vertices.push_back(transform(gctx) * Vector3(v2D.x(), v2D.y(), 0.));
0169       wCenter += (*vertices.rbegin());
0170     }
0171     // These are convex shapes, use the helper method
0172     // For rings there's a sweet spot when this stops working
0173     if (exactPolyhedron || toCenter || !fullDisc) {
0174       // Transform them into the vertex frame
0175       wCenter *= 1. / vertices.size();
0176       if (addCentreFromConvexFace) {
0177         vertices.push_back(wCenter);
0178       }
0179       auto [faces, triangularMesh] =
0180           detail::FacesHelper::convexFaceMesh(vertices, true);
0181       return Polyhedron(vertices, faces, triangularMesh, exactPolyhedron);
0182     } else {
0183       // Two concentric rings, we use the pure concentric method momentarily,
0184       // but that creates too  many unneccesarry faces, when only two
0185       // are needed to describe the mesh, @todo investigate merging flag
0186       auto [faces, triangularMesh] =
0187           detail::FacesHelper::cylindricalFaceMesh(vertices);
0188       return Polyhedron(vertices, faces, triangularMesh, exactPolyhedron);
0189     }
0190   }
0191   throw std::domain_error("Polyhedron repr of boundless surface not possible.");
0192 }
0193
0194 Vector2 DiscSurface::localPolarToCartesian(const Vector2& lpolar) const {
0195   return Vector2(lpolar[0] * cos(lpolar[1]), lpolar[0] * sin(lpolar[1]));
0196 }
0197
0198 Vector2 DiscSurface::localCartesianToPolar(const Vector2& lcart) const {
0199   return Vector2(lcart.norm(), std::atan2(lcart[1], lcart[0]));
0200 }
0201
0202 BoundToFreeMatrix DiscSurface::boundToFreeJacobian(
0203     const GeometryContext& gctx, const Vector3& position,
0204     const Vector3& direction) const {
0205   assert(isOnSurface(gctx, position, direction, BoundaryTolerance::Infinite()));
0206
0207   // The measurement frame of the surface
0208   RotationMatrix3 rframeT =
0209       referenceFrame(gctx, position, direction).transpose();
0210
0211   // calculate the transformation to local coordinates
0212   const Vector3 posLoc = transform(gctx).inverse() * position;
0213   const double lr = perp(posLoc);
0214   const double lphi = phi(posLoc);
0215   const double lcphi = std::cos(lphi);
0216   const double lsphi = std::sin(lphi);
0217   // rotate into the polar coorindates
0218   auto lx = rframeT.block<1, 3>(0, 0);
0219   auto ly = rframeT.block<1, 3>(1, 0);
0220   // Initialize the jacobian from local to global
0221   BoundToFreeMatrix jacToGlobal = BoundToFreeMatrix::Zero();
0222   // the local error components - rotated from reference frame
0223   jacToGlobal.block<3, 1>(eFreePos0, eBoundLoc0) = lcphi * lx + lsphi * ly;
0224   jacToGlobal.block<3, 1>(eFreePos0, eBoundLoc1) =
0225       lr * (lcphi * ly - lsphi * lx);
0226   // the time component
0227   jacToGlobal(eFreeTime, eBoundTime) = 1;
0228   // the momentum components
0229   jacToGlobal.block<3, 2>(eFreeDir0, eBoundPhi) =
0230       sphericalToFreeDirectionJacobian(direction);
0231   jacToGlobal(eFreeQOverP, eBoundQOverP) = 1;
0232   return jacToGlobal;
0233 }
0234
0235 FreeToBoundMatrix DiscSurface::freeToBoundJacobian(
0236     const GeometryContext& gctx, const Vector3& position,
0237     const Vector3& direction) const {
0238   using VectorHelpers::perp;
0239   using VectorHelpers::phi;
0240
0241   assert(isOnSurface(gctx, position, direction, BoundaryTolerance::Infinite()));
0242
0243   // The measurement frame of the surface
0244   RotationMatrix3 rframeT =
0245       referenceFrame(gctx, position, direction).transpose();
0246
0247   // calculate the transformation to local coordinates
0248   const Vector3 posLoc = transform(gctx).inverse() * position;
0249   const double lr = perp(posLoc);
0250   const double lphi = phi(posLoc);
0251   const double lcphi = std::cos(lphi);
0252   const double lsphi = std::sin(lphi);
0253   // rotate into the polar coorindates
0254   auto lx = rframeT.block<1, 3>(0, 0);
0255   auto ly = rframeT.block<1, 3>(1, 0);
0256   // Initialize the jacobian from global to local
0257   FreeToBoundMatrix jacToLocal = FreeToBoundMatrix::Zero();
0258   // Local position component
0259   jacToLocal.block<1, 3>(eBoundLoc0, eFreePos0) = lcphi * lx + lsphi * ly;
0260   jacToLocal.block<1, 3>(eBoundLoc1, eFreePos0) =
0261       (lcphi * ly - lsphi * lx) / lr;
0262   // Time element
0263   jacToLocal(eBoundTime, eFreeTime) = 1;
0264   // Directional and momentum elements for reference frame surface
0265   jacToLocal.block<2, 3>(eBoundPhi, eFreeDir0) =
0266       freeToSphericalDirectionJacobian(direction);
0267   jacToLocal(eBoundQOverP, eFreeQOverP) = 1;
0268   return jacToLocal;
0269 }
0270
0271 SurfaceMultiIntersection DiscSurface::intersect(
0272     const GeometryContext& gctx, const Vector3& position,
0273     const Vector3& direction, const BoundaryTolerance& boundaryTolerance,
0274     double tolerance) const {
0275   // Get the contextual transform
0276   auto gctxTransform = transform(gctx);
0277   // Use the intersection helper for planar surfaces
0278   auto intersection =
0279       PlanarHelper::intersect(gctxTransform, position, direction, tolerance);
0280   auto status = intersection.status();
0281   // Evaluate boundary check if requested (and reachable)
0282   if (intersection.status() != IntersectionStatus::unreachable &&
0283       m_bounds != nullptr && !boundaryTolerance.isInfinite()) {
0284     // Built-in local to global for speed reasons
0285     const auto& tMatrix = gctxTransform.matrix();
0286     const Vector3 vecLocal(intersection.position() - tMatrix.block<3, 1>(0, 3));
0287     const Vector2 lcartesian = tMatrix.block<3, 2>(0, 0).transpose() * vecLocal;
0288     if (auto absoluteBound = boundaryTolerance.asAbsoluteBoundOpt();
0289         absoluteBound.has_value() && m_bounds->coversFullAzimuth()) {
0290       double modifiedTolerance = tolerance + absoluteBound->tolerance0;
0291       if (!m_bounds->insideRadialBounds(VectorHelpers::perp(lcartesian),
0292                                         modifiedTolerance)) {
0293         status = IntersectionStatus::unreachable;
0294       }
0295     } else if (!insideBounds(localCartesianToPolar(lcartesian),
0296                              boundaryTolerance)) {
0297       status = IntersectionStatus::unreachable;
0298     }
0299   }
0300   return {{Intersection3D(intersection.position(), intersection.pathLength(),
0301                           status),
0302            Intersection3D::invalid()},
0303           this};
0304 }
0305
0306 ActsMatrix<2, 3> DiscSurface::localCartesianToBoundLocalDerivative(
0307     const GeometryContext& gctx, const Vector3& position) const {
0308   using VectorHelpers::perp;
0309   using VectorHelpers::phi;
0310   // The local frame transform
0311   const auto& sTransform = transform(gctx);
0312   // calculate the transformation to local coordinates
0313   const Vector3 localPos = sTransform.inverse() * position;
0314   const double lr = perp(localPos);
0315   const double lphi = phi(localPos);
0316   const double lcphi = std::cos(lphi);
0317   const double lsphi = std::sin(lphi);
0318   ActsMatrix<2, 3> loc3DToLocBound = ActsMatrix<2, 3>::Zero();
0319   loc3DToLocBound << lcphi, lsphi, 0, -lsphi / lr, lcphi / lr, 0;
0320
0321   return loc3DToLocBound;
0322 }
0323
0324 Vector3 DiscSurface::normal(const GeometryContext& gctx,
0325                             const Vector2& /*lposition*/) const {
0326   return normal(gctx);
0327 }
0328
0329 Vector3 DiscSurface::normal(const GeometryContext& gctx,
0330                             const Vector3& /*position*/) const {
0331   return normal(gctx);
0332 }
0333
0334 Vector3 DiscSurface::normal(const GeometryContext& gctx) const {
0335   // fast access via transform matrix (and not rotation())
0336   const auto& tMatrix = transform(gctx).matrix();
0337   return Vector3(tMatrix(0, 2), tMatrix(1, 2), tMatrix(2, 2));
0338 }
0339
0340 Vector3 DiscSurface::referencePosition(const GeometryContext& gctx,
0341                                        AxisDirection aDir) const {
0342   if (aDir == AxisDirection::AxisR || aDir == AxisDirection::AxisPhi) {
0343     double r = m_bounds->binningValueR();
0344     double phi = m_bounds->binningValuePhi();
0345     return localToGlobal(gctx, Vector2{r, phi}, Vector3{});
0346   }
0347   return center(gctx);
0348 }
0349
0350 double DiscSurface::referencePositionValue(const GeometryContext& gctx,
0351                                            AxisDirection aDir) const {
0352   if (aDir == AxisDirection::AxisR) {
0353     return VectorHelpers::perp(referencePosition(gctx, aDir));
0354   }
0355   if (aDir == AxisDirection::AxisPhi) {
0356     return VectorHelpers::phi(referencePosition(gctx, aDir));
0357   }
0358
0359   return GeometryObject::referencePositionValue(gctx, aDir);
0360 }
0361
0362 double DiscSurface::pathCorrection(const GeometryContext& gctx,
0363                                    const Vector3& /*position*/,
0364                                    const Vector3& direction) const {
0365   // we can ignore the global position here
0366   return 1. / std::abs(normal(gctx).dot(direction));
0367 }
0368
0369 std::pair<std::shared_ptr<DiscSurface>, bool> DiscSurface::mergedWith(
0370     const DiscSurface& other, AxisDirection direction, bool externalRotation,
0371     const Logger& logger) const {
0372   using namespace UnitLiterals;
0373
0374   ACTS_VERBOSE("Merging disc surfaces in " << direction << " direction");
0375
0376   if (m_associatedDetElement != nullptr ||
0377       other.m_associatedDetElement != nullptr) {
0378     throw SurfaceMergingException(getSharedPtr(), other.getSharedPtr(),
0379                                   "CylinderSurface::merge: surfaces are "
0380                                   "associated with a detector element");
0381   }
0382   assert(m_transform != nullptr && other.m_transform != nullptr);
0383
0384   Transform3 otherLocal = m_transform->inverse() * *other.m_transform;
0385
0386   constexpr auto tolerance = s_onSurfaceTolerance;
0387
0388   // surface cannot have any relative rotation
0389   if (std::abs(otherLocal.linear().col(eX)[eZ]) >= tolerance ||
0390       std::abs(otherLocal.linear().col(eY)[eZ]) >= tolerance) {
0391     ACTS_ERROR("DiscSurface::merge: surfaces have relative rotation");
0392     throw SurfaceMergingException(
0393         getSharedPtr(), other.getSharedPtr(),
0394         "DiscSurface::merge: surfaces have relative rotation");
0395   }
0396
0397   Vector3 translation = otherLocal.translation();
0398
0399   if (std::abs(translation[0]) > tolerance ||
0400       std::abs(translation[1]) > tolerance ||
0401       std::abs(translation[2]) > tolerance) {
0402     ACTS_ERROR(
0403         "DiscSurface::merge: surfaces have relative translation in x/y/z");
0404     throw SurfaceMergingException(
0405         getSharedPtr(), other.getSharedPtr(),
0406         "DiscSurface::merge: surfaces have relative translation in x/y/z");
0407   }
0408
0409   const auto* bounds = dynamic_cast<const RadialBounds*>(m_bounds.get());
0410   const auto* otherBounds =
0411       dynamic_cast<const RadialBounds*>(other.m_bounds.get());
0412
0413   if (bounds == nullptr || otherBounds == nullptr) {
0414     ACTS_ERROR("DiscSurface::merge: surfaces have bounds other than radial");
0415     throw SurfaceMergingException(
0416         getSharedPtr(), other.getSharedPtr(),
0417         "DiscSurface::merge: surfaces have bounds other than radial");
0418   }
0419
0420   double minR = bounds->get(RadialBounds::eMinR);
0421   double maxR = bounds->get(RadialBounds::eMaxR);
0422
0423   double hlPhi = bounds->get(RadialBounds::eHalfPhiSector);
0424   double avgPhi = bounds->get(RadialBounds::eAveragePhi);
0425   double minPhi = detail::radian_sym(-hlPhi + avgPhi);
0426   double maxPhi = detail::radian_sym(hlPhi + avgPhi);
0427
0428   ACTS_VERBOSE(" this: r =   [" << minR << ", " << maxR << "]");
0429   ACTS_VERBOSE("       phi = ["
0430                << minPhi / 1_degree << ", " << maxPhi / 1_degree << "] ~> "
0431                << avgPhi / 1_degree << " +- " << hlPhi / 1_degree);
0432
0433   double otherMinR = otherBounds->get(RadialBounds::eMinR);
0434   double otherMaxR = otherBounds->get(RadialBounds::eMaxR);
0435   double otherAvgPhi = otherBounds->get(RadialBounds::eAveragePhi);
0436   double otherHlPhi = otherBounds->get(RadialBounds::eHalfPhiSector);
0437   double otherMinPhi = detail::radian_sym(-otherHlPhi + otherAvgPhi);
0438   double otherMaxPhi = detail::radian_sym(otherHlPhi + otherAvgPhi);
0439
0440   ACTS_VERBOSE("other: r =   [" << otherMinR << ", " << otherMaxR << "]");
0441   ACTS_VERBOSE("       phi = [" << otherMinPhi / 1_degree << ", "
0442                                 << otherMaxPhi / 1_degree << "] ~> "
0443                                 << otherAvgPhi / 1_degree << " +- "
0444                                 << otherHlPhi / 1_degree);
0445
0446   if (direction == AxisDirection::AxisR) {
0447     if (std::abs(otherLocal.linear().col(eY)[eX]) >= tolerance &&
0448         (!bounds->coversFullAzimuth() || !otherBounds->coversFullAzimuth())) {
0449       throw SurfaceMergingException(getSharedPtr(), other.getSharedPtr(),
0450                                     "DiscSurface::merge: surfaces have "
0451                                     "relative rotation in z and phi sector");
0452     }
0453
0454     if (std::abs(minR - otherMaxR) > tolerance &&
0455         std::abs(maxR - otherMinR) > tolerance) {
0456       ACTS_ERROR("DiscSurface::merge: surfaces are not touching r");
0457       throw SurfaceMergingException(
0458           getSharedPtr(), other.getSharedPtr(),
0459           "DiscSurface::merge: surfaces are not touching in r");
0460     }
0461
0462     if (std::abs(avgPhi - otherAvgPhi) > tolerance) {
0463       ACTS_ERROR("DiscSurface::merge: surfaces have different average phi");
0464       throw SurfaceMergingException(
0465           getSharedPtr(), other.getSharedPtr(),
0466           "DiscSurface::merge: surfaces have different average phi");
0467     }
0468
0469     if (std::abs(hlPhi - otherHlPhi) > tolerance) {
0470       ACTS_ERROR("DiscSurface::merge: surfaces have different half phi sector");
0471       throw SurfaceMergingException(
0472           getSharedPtr(), other.getSharedPtr(),
0473           "DiscSurface::merge: surfaces have different half phi sector");
0474     }
0475
0476     double newMinR = std::min(minR, otherMinR);
0477     double newMaxR = std::max(maxR, otherMaxR);
0478     ACTS_VERBOSE("  new: r =   [" << newMinR << ", " << newMaxR << "]");
0479
0480     auto newBounds =
0481         std::make_shared<RadialBounds>(newMinR, newMaxR, hlPhi, avgPhi);
0482
0483     return {Surface::makeShared<DiscSurface>(*m_transform, newBounds),
0484             minR > otherMinR};
0485
0486   } else if (direction == AxisDirection::AxisPhi) {
0487     if (std::abs(maxR - otherMaxR) > tolerance ||
0488         std::abs(minR - otherMinR) > tolerance) {
0489       ACTS_ERROR("DiscSurface::merge: surfaces don't have same r bounds");
0490       throw SurfaceMergingException(
0491           getSharedPtr(), other.getSharedPtr(),
0492           "DiscSurface::merge: surfaces don't have same r bounds");
0493     }
0494
0495     // Figure out signed relative rotation
0496     Vector2 rotatedX = otherLocal.linear().col(eX).head<2>();
0497     double zrotation = std::atan2(rotatedX[1], rotatedX[0]);
0498
0499     ACTS_VERBOSE("this:  [" << avgPhi / 1_degree << " +- " << hlPhi / 1_degree
0500                             << "]");
0501     ACTS_VERBOSE("other: [" << otherAvgPhi / 1_degree << " +- "
0502                             << otherHlPhi / 1_degree << "]");
0503
0504     ACTS_VERBOSE("Relative rotation around local z: " << zrotation / 1_degree);
0505
0506     double prevOtherAvgPhi = otherAvgPhi;
0507     otherAvgPhi = detail::radian_sym(otherAvgPhi + zrotation);
0508     ACTS_VERBOSE("~> local other average phi: "
0509                  << otherAvgPhi / 1_degree
0510                  << " (was: " << prevOtherAvgPhi / 1_degree << ")");
0511
0512     try {
0513       auto [newHlPhi, newAvgPhi, reversed] = detail::mergedPhiSector(
0514           hlPhi, avgPhi, otherHlPhi, otherAvgPhi, logger, tolerance);
0515
0516       Transform3 newTransform = *m_transform;
0517
0518       if (externalRotation) {
0519         ACTS_VERBOSE("Modifying transform for external rotation of "
0520                      << newAvgPhi / 1_degree);
0521         newTransform = newTransform * AngleAxis3(newAvgPhi, Vector3::UnitZ());
0522         newAvgPhi = 0.;
0523       }
0524
0525       auto newBounds =
0526           std::make_shared<RadialBounds>(minR, maxR, newHlPhi, newAvgPhi);
0527
0528       return {Surface::makeShared<DiscSurface>(newTransform, newBounds),
0529               reversed};
0530     } catch (const std::invalid_argument& e) {
0531       throw SurfaceMergingException(getSharedPtr(), other.getSharedPtr(),
0532                                     e.what());
0533     }
0534
0535   } else {
0536     ACTS_ERROR("DiscSurface::merge: invalid direction " << direction);
0537
0538     throw SurfaceMergingException(getSharedPtr(), other.getSharedPtr(),
0539                                   "DiscSurface::merge: invalid direction " +
0540                                       axisDirectionName(direction));
0541   }
0542 }
0543
0544 }  // namespace Acts