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/Intersection.hpp"
0026 #include "Acts/Utilities/JacobianHelpers.hpp"
0027 #include "Acts/Utilities/MathHelpers.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] * std::cos(lposition[1]),
0094                      lposition[0] * std::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] * std::sin(phi) - pos[1] * std::cos(phi),
0123                    pos[1] * std::sin(phi) + pos[0] * std::cos(phi));
0124     return Vector2(locPos[0], locPos[1]);
0125   }
0126   return Vector2(locpol[0] * std::cos(locpol[1]),
0127                  locpol[0] * std::sin(locpol[1]));
0128 }
0129
0130 Vector3 DiscSurface::localCartesianToGlobal(const GeometryContext& gctx,
0131                                             const Vector2& lposition) const {
0132   Vector3 loc3Dframe(lposition[0], lposition[1], 0.);
0133   return transform(gctx) * loc3Dframe;
0134 }
0135
0136 Vector2 DiscSurface::globalToLocalCartesian(const GeometryContext& gctx,
0137                                             const Vector3& position,
0138                                             double /*direction*/) const {
0139   Vector3 loc3Dframe = (transform(gctx).inverse()) * position;
0140   return Vector2(loc3Dframe.x(), loc3Dframe.y());
0141 }
0142
0143 std::string DiscSurface::name() const {
0144   return "Acts::DiscSurface";
0145 }
0146
0147 const SurfaceBounds& DiscSurface::bounds() const {
0148   if (m_bounds) {
0149     return *m_bounds;
0150   }
0151   return s_noBounds;
0152 }
0153
0154 Polyhedron DiscSurface::polyhedronRepresentation(
0155     const GeometryContext& gctx, unsigned int quarterSegments) const {
0156   // Prepare vertices and faces
0157   std::vector<Vector3> vertices;
0158   // Understand the disc
0159   bool fullDisc = m_bounds->coversFullAzimuth();
0160   bool toCenter = m_bounds->rMin() < s_onSurfaceTolerance;
0161   // If you have bounds you can create a polyhedron representation
0162   bool exactPolyhedron = (m_bounds->type() == SurfaceBounds::eDiscTrapezoid);
0163   bool addCentreFromConvexFace = (m_bounds->type() != SurfaceBounds::eAnnulus);
0164   if (m_bounds) {
0165     auto vertices2D = m_bounds->vertices(quarterSegments);
0166     vertices.reserve(vertices2D.size() + 1);
0167     Vector3 wCenter(0., 0., 0);
0168     for (const auto& v2D : vertices2D) {
0169       vertices.push_back(transform(gctx) * Vector3(v2D.x(), v2D.y(), 0.));
0170       wCenter += (*vertices.rbegin());
0171     }
0172     // These are convex shapes, use the helper method
0173     // For rings there's a sweet spot when this stops working
0174     if (exactPolyhedron || toCenter || !fullDisc) {
0175       // Transform them into the vertex frame
0176       wCenter *= 1. / vertices.size();
0177       if (addCentreFromConvexFace) {
0178         vertices.push_back(wCenter);
0179       }
0180       auto [faces, triangularMesh] = detail::FacesHelper::convexFaceMesh(
0181           vertices, addCentreFromConvexFace);
0182       return Polyhedron(vertices, faces, triangularMesh, exactPolyhedron);
0183     } else {
0184       // Two concentric rings, we use the pure concentric method momentarily,
0185       // but that creates too  many unneccesarry faces, when only two
0186       // are needed to describe the mesh, @todo investigate merging flag
0187       auto [faces, triangularMesh] =
0188           detail::FacesHelper::cylindricalFaceMesh(vertices);
0189       return Polyhedron(vertices, faces, triangularMesh, exactPolyhedron);
0190     }
0191   }
0192   throw std::domain_error("Polyhedron repr of boundless surface not possible.");
0193 }
0194
0195 Vector2 DiscSurface::localPolarToCartesian(const Vector2& lpolar) const {
0196   return Vector2(lpolar[0] * std::cos(lpolar[1]),
0197                  lpolar[0] * std::sin(lpolar[1]));
0198 }
0199
0200 Vector2 DiscSurface::localCartesianToPolar(const Vector2& lcart) const {
0201   return Vector2(lcart.norm(), std::atan2(lcart[1], lcart[0]));
0202 }
0203
0204 BoundToFreeMatrix DiscSurface::boundToFreeJacobian(
0205     const GeometryContext& gctx, const Vector3& position,
0206     const Vector3& direction) const {
0207   assert(isOnSurface(gctx, position, direction, BoundaryTolerance::Infinite()));
0208
0209   // The measurement frame of the surface
0210   RotationMatrix3 rframeT =
0211       referenceFrame(gctx, position, direction).transpose();
0212
0213   // calculate the transformation to local coordinates
0214   const Vector3 posLoc = transform(gctx).inverse() * position;
0215   const double lr = perp(posLoc);
0216   const double lphi = phi(posLoc);
0217   const double lcphi = std::cos(lphi);
0218   const double lsphi = std::sin(lphi);
0219   // rotate into the polar coorindates
0220   auto lx = rframeT.block<1, 3>(0, 0);
0221   auto ly = rframeT.block<1, 3>(1, 0);
0222   // Initialize the jacobian from local to global
0223   BoundToFreeMatrix jacToGlobal = BoundToFreeMatrix::Zero();
0224   // the local error components - rotated from reference frame
0225   jacToGlobal.block<3, 1>(eFreePos0, eBoundLoc0) = lcphi * lx + lsphi * ly;
0226   jacToGlobal.block<3, 1>(eFreePos0, eBoundLoc1) =
0227       lr * (lcphi * ly - lsphi * lx);
0228   // the time component
0229   jacToGlobal(eFreeTime, eBoundTime) = 1;
0230   // the momentum components
0231   jacToGlobal.block<3, 2>(eFreeDir0, eBoundPhi) =
0232       sphericalToFreeDirectionJacobian(direction);
0233   jacToGlobal(eFreeQOverP, eBoundQOverP) = 1;
0234   return jacToGlobal;
0235 }
0236
0237 FreeToBoundMatrix DiscSurface::freeToBoundJacobian(
0238     const GeometryContext& gctx, const Vector3& position,
0239     const Vector3& direction) const {
0240   using VectorHelpers::perp;
0241   using VectorHelpers::phi;
0242
0243   assert(isOnSurface(gctx, position, direction, BoundaryTolerance::Infinite()));
0244
0245   // The measurement frame of the surface
0246   RotationMatrix3 rframeT =
0247       referenceFrame(gctx, position, direction).transpose();
0248
0249   // calculate the transformation to local coordinates
0250   const Vector3 posLoc = transform(gctx).inverse() * position;
0251   const double lr = perp(posLoc);
0252   const double lphi = phi(posLoc);
0253   const double lcphi = std::cos(lphi);
0254   const double lsphi = std::sin(lphi);
0255   // rotate into the polar coorindates
0256   auto lx = rframeT.block<1, 3>(0, 0);
0257   auto ly = rframeT.block<1, 3>(1, 0);
0258   // Initialize the jacobian from global to local
0259   FreeToBoundMatrix jacToLocal = FreeToBoundMatrix::Zero();
0260   // Local position component
0261   jacToLocal.block<1, 3>(eBoundLoc0, eFreePos0) = lcphi * lx + lsphi * ly;
0262   jacToLocal.block<1, 3>(eBoundLoc1, eFreePos0) =
0263       (lcphi * ly - lsphi * lx) / lr;
0264   // Time element
0265   jacToLocal(eBoundTime, eFreeTime) = 1;
0266   // Directional and momentum elements for reference frame surface
0267   jacToLocal.block<2, 3>(eBoundPhi, eFreeDir0) =
0268       freeToSphericalDirectionJacobian(direction);
0269   jacToLocal(eBoundQOverP, eFreeQOverP) = 1;
0270   return jacToLocal;
0271 }
0272
0273 MultiIntersection3D DiscSurface::intersect(
0274     const GeometryContext& gctx, const Vector3& position,
0275     const Vector3& direction, const BoundaryTolerance& boundaryTolerance,
0276     double tolerance) const {
0277   // Get the contextual transform
0278   const Transform3& gctxTransform = transform(gctx);
0279   // Use the intersection helper for planar surfaces
0280   const Intersection3D intersection =
0281       PlanarHelper::intersect(gctxTransform, position, direction, tolerance);
0282   IntersectionStatus status = intersection.status();
0283   if (status == IntersectionStatus::unreachable) {
0284     return MultiIntersection3D(Intersection3D::Invalid());
0285   }
0286   if (m_bounds == nullptr || boundaryTolerance.isInfinite()) {
0287     return MultiIntersection3D(intersection);
0288   }
0289   // Built-in local to global for speed reasons
0290   const auto& tMatrix = gctxTransform.matrix();
0291   const Vector3 fromCenter =
0292       intersection.position() - tMatrix.block<3, 1>(0, 3);
0293   if (m_bounds->coversFullAzimuth() && boundaryTolerance.isNone()) {
0294     // avoids `atan2` in case of full phi coverage
0295     const double r2 = fromCenter.squaredNorm();
0296     const bool isInside =
0297         (r2 >= square(m_bounds->rMin())) && (r2 <= square(m_bounds->rMax()));
0298     if (!isInside) {
0299       status = IntersectionStatus::unreachable;
0300     }
0301   } else {
0302     const Vector2 localCartesian =
0303         tMatrix.block<3, 2>(0, 0).transpose() * fromCenter;
0304     const bool isInside =
0305         insideBounds(localCartesianToPolar(localCartesian), boundaryTolerance);
0306     if (!isInside) {
0307       status = IntersectionStatus::unreachable;
0308     }
0309   }
0310   return MultiIntersection3D(Intersection3D(intersection.position(),
0311                                             intersection.pathLength(), status));
0312 }
0313
0314 ActsMatrix<2, 3> DiscSurface::localCartesianToBoundLocalDerivative(
0315     const GeometryContext& gctx, const Vector3& position) const {
0316   using VectorHelpers::perp;
0317   using VectorHelpers::phi;
0318   // The local frame transform
0319   const auto& sTransform = transform(gctx);
0320   // calculate the transformation to local coordinates
0321   const Vector3 localPos = sTransform.inverse() * position;
0322   const double lr = perp(localPos);
0323   const double lphi = phi(localPos);
0324   const double lcphi = std::cos(lphi);
0325   const double lsphi = std::sin(lphi);
0326   ActsMatrix<2, 3> loc3DToLocBound = ActsMatrix<2, 3>::Zero();
0327   loc3DToLocBound << lcphi, lsphi, 0, -lsphi / lr, lcphi / lr, 0;
0328
0329   return loc3DToLocBound;
0330 }
0331
0332 Vector3 DiscSurface::normal(const GeometryContext& gctx,
0333                             const Vector2& /*lposition*/) const {
0334   return normal(gctx);
0335 }
0336
0337 Vector3 DiscSurface::normal(const GeometryContext& gctx,
0338                             const Vector3& /*position*/) const {
0339   return normal(gctx);
0340 }
0341
0342 Vector3 DiscSurface::normal(const GeometryContext& gctx) const {
0343   // fast access via transform matrix (and not rotation())
0344   const auto& tMatrix = transform(gctx).matrix();
0345   return Vector3(tMatrix(0, 2), tMatrix(1, 2), tMatrix(2, 2));
0346 }
0347
0348 Vector3 DiscSurface::referencePosition(const GeometryContext& gctx,
0349                                        AxisDirection aDir) const {
0350   if (aDir == AxisDirection::AxisR || aDir == AxisDirection::AxisPhi) {
0351     double r = m_bounds->binningValueR();
0352     double phi = m_bounds->binningValuePhi();
0353     return localToGlobal(gctx, Vector2{r, phi}, Vector3{});
0354   }
0355   return center(gctx);
0356 }
0357
0358 double DiscSurface::referencePositionValue(const GeometryContext& gctx,
0359                                            AxisDirection aDir) const {
0360   if (aDir == AxisDirection::AxisR) {
0361     return VectorHelpers::perp(referencePosition(gctx, aDir));
0362   }
0363   if (aDir == AxisDirection::AxisPhi) {
0364     return VectorHelpers::phi(referencePosition(gctx, aDir));
0365   }
0366
0367   return GeometryObject::referencePositionValue(gctx, aDir);
0368 }
0369
0370 double DiscSurface::pathCorrection(const GeometryContext& gctx,
0371                                    const Vector3& /*position*/,
0372                                    const Vector3& direction) const {
0373   // we can ignore the global position here
0374   return 1. / std::abs(normal(gctx).dot(direction));
0375 }
0376
0377 std::pair<std::shared_ptr<DiscSurface>, bool> DiscSurface::mergedWith(
0378     const DiscSurface& other, AxisDirection direction, bool externalRotation,
0379     const Logger& logger) const {
0380   using namespace UnitLiterals;
0381
0382   ACTS_VERBOSE("Merging disc surfaces in " << direction << " direction");
0383
0384   if (m_associatedDetElement != nullptr ||
0385       other.m_associatedDetElement != nullptr) {
0386     throw SurfaceMergingException(getSharedPtr(), other.getSharedPtr(),
0387                                   "CylinderSurface::merge: surfaces are "
0388                                   "associated with a detector element");
0389   }
0390   assert(m_transform != nullptr && other.m_transform != nullptr);
0391
0392   Transform3 otherLocal = m_transform->inverse() * *other.m_transform;
0393
0394   constexpr auto tolerance = s_onSurfaceTolerance;
0395
0396   // surface cannot have any relative rotation
0397   if (std::abs(otherLocal.linear().col(eX)[eZ]) >= tolerance ||
0398       std::abs(otherLocal.linear().col(eY)[eZ]) >= tolerance) {
0399     ACTS_ERROR("DiscSurface::merge: surfaces have relative rotation");
0400     throw SurfaceMergingException(
0401         getSharedPtr(), other.getSharedPtr(),
0402         "DiscSurface::merge: surfaces have relative rotation");
0403   }
0404
0405   Vector3 translation = otherLocal.translation();
0406
0407   if (std::abs(translation[0]) > tolerance ||
0408       std::abs(translation[1]) > tolerance ||
0409       std::abs(translation[2]) > tolerance) {
0410     ACTS_ERROR(
0411         "DiscSurface::merge: surfaces have relative translation in x/y/z");
0412     throw SurfaceMergingException(
0413         getSharedPtr(), other.getSharedPtr(),
0414         "DiscSurface::merge: surfaces have relative translation in x/y/z");
0415   }
0416
0417   const auto* bounds = dynamic_cast<const RadialBounds*>(m_bounds.get());
0418   const auto* otherBounds =
0419       dynamic_cast<const RadialBounds*>(other.m_bounds.get());
0420
0421   if (bounds == nullptr || otherBounds == nullptr) {
0422     ACTS_ERROR("DiscSurface::merge: surfaces have bounds other than radial");
0423     throw SurfaceMergingException(
0424         getSharedPtr(), other.getSharedPtr(),
0425         "DiscSurface::merge: surfaces have bounds other than radial");
0426   }
0427
0428   double minR = bounds->get(RadialBounds::eMinR);
0429   double maxR = bounds->get(RadialBounds::eMaxR);
0430
0431   double hlPhi = bounds->get(RadialBounds::eHalfPhiSector);
0432   double avgPhi = bounds->get(RadialBounds::eAveragePhi);
0433   double minPhi = detail::radian_sym(-hlPhi + avgPhi);
0434   double maxPhi = detail::radian_sym(hlPhi + avgPhi);
0435
0436   ACTS_VERBOSE(" this: r =   [" << minR << ", " << maxR << "]");
0437   ACTS_VERBOSE("       phi = ["
0438                << minPhi / 1_degree << ", " << maxPhi / 1_degree << "] ~> "
0439                << avgPhi / 1_degree << " +- " << hlPhi / 1_degree);
0440
0441   double otherMinR = otherBounds->get(RadialBounds::eMinR);
0442   double otherMaxR = otherBounds->get(RadialBounds::eMaxR);
0443   double otherAvgPhi = otherBounds->get(RadialBounds::eAveragePhi);
0444   double otherHlPhi = otherBounds->get(RadialBounds::eHalfPhiSector);
0445   double otherMinPhi = detail::radian_sym(-otherHlPhi + otherAvgPhi);
0446   double otherMaxPhi = detail::radian_sym(otherHlPhi + otherAvgPhi);
0447
0448   ACTS_VERBOSE("other: r =   [" << otherMinR << ", " << otherMaxR << "]");
0449   ACTS_VERBOSE("       phi = [" << otherMinPhi / 1_degree << ", "
0450                                 << otherMaxPhi / 1_degree << "] ~> "
0451                                 << otherAvgPhi / 1_degree << " +- "
0452                                 << otherHlPhi / 1_degree);
0453
0454   if (direction == AxisDirection::AxisR) {
0455     if (std::abs(otherLocal.linear().col(eY)[eX]) >= tolerance &&
0456         (!bounds->coversFullAzimuth() || !otherBounds->coversFullAzimuth())) {
0457       throw SurfaceMergingException(getSharedPtr(), other.getSharedPtr(),
0458                                     "DiscSurface::merge: surfaces have "
0459                                     "relative rotation in z and phi sector");
0460     }
0461
0462     if (std::abs(minR - otherMaxR) > tolerance &&
0463         std::abs(maxR - otherMinR) > tolerance) {
0464       ACTS_ERROR("DiscSurface::merge: surfaces are not touching r");
0465       throw SurfaceMergingException(
0466           getSharedPtr(), other.getSharedPtr(),
0467           "DiscSurface::merge: surfaces are not touching in r");
0468     }
0469
0470     if (std::abs(avgPhi - otherAvgPhi) > tolerance) {
0471       ACTS_ERROR("DiscSurface::merge: surfaces have different average phi");
0472       throw SurfaceMergingException(
0473           getSharedPtr(), other.getSharedPtr(),
0474           "DiscSurface::merge: surfaces have different average phi");
0475     }
0476
0477     if (std::abs(hlPhi - otherHlPhi) > tolerance) {
0478       ACTS_ERROR("DiscSurface::merge: surfaces have different half phi sector");
0479       throw SurfaceMergingException(
0480           getSharedPtr(), other.getSharedPtr(),
0481           "DiscSurface::merge: surfaces have different half phi sector");
0482     }
0483
0484     double newMinR = std::min(minR, otherMinR);
0485     double newMaxR = std::max(maxR, otherMaxR);
0486     ACTS_VERBOSE("  new: r =   [" << newMinR << ", " << newMaxR << "]");
0487
0488     auto newBounds =
0489         std::make_shared<RadialBounds>(newMinR, newMaxR, hlPhi, avgPhi);
0490
0491     return {Surface::makeShared<DiscSurface>(*m_transform, newBounds),
0492             minR > otherMinR};
0493
0494   } else if (direction == AxisDirection::AxisPhi) {
0495     if (std::abs(maxR - otherMaxR) > tolerance ||
0496         std::abs(minR - otherMinR) > tolerance) {
0497       ACTS_ERROR("DiscSurface::merge: surfaces don't have same r bounds");
0498       throw SurfaceMergingException(
0499           getSharedPtr(), other.getSharedPtr(),
0500           "DiscSurface::merge: surfaces don't have same r bounds");
0501     }
0502
0503     // Figure out signed relative rotation
0504     Vector2 rotatedX = otherLocal.linear().col(eX).head<2>();
0505     double zrotation = std::atan2(rotatedX[1], rotatedX[0]);
0506
0507     ACTS_VERBOSE("this:  [" << avgPhi / 1_degree << " +- " << hlPhi / 1_degree
0508                             << "]");
0509     ACTS_VERBOSE("other: [" << otherAvgPhi / 1_degree << " +- "
0510                             << otherHlPhi / 1_degree << "]");
0511
0512     ACTS_VERBOSE("Relative rotation around local z: " << zrotation / 1_degree);
0513
0514     double prevOtherAvgPhi = otherAvgPhi;
0515     otherAvgPhi = detail::radian_sym(otherAvgPhi + zrotation);
0516     ACTS_VERBOSE("~> local other average phi: "
0517                  << otherAvgPhi / 1_degree
0518                  << " (was: " << prevOtherAvgPhi / 1_degree << ")");
0519
0520     try {
0521       auto [newHlPhi, newAvgPhi, reversed] = detail::mergedPhiSector(
0522           hlPhi, avgPhi, otherHlPhi, otherAvgPhi, logger, tolerance);
0523
0524       Transform3 newTransform = *m_transform;
0525
0526       if (externalRotation) {
0527         ACTS_VERBOSE("Modifying transform for external rotation of "
0528                      << newAvgPhi / 1_degree);
0529         newTransform = newTransform * AngleAxis3(newAvgPhi, Vector3::UnitZ());
0530         newAvgPhi = 0.;
0531       }
0532
0533       auto newBounds =
0534           std::make_shared<RadialBounds>(minR, maxR, newHlPhi, newAvgPhi);
0535
0536       return {Surface::makeShared<DiscSurface>(newTransform, newBounds),
0537               reversed};
0538     } catch (const std::invalid_argument& e) {
0539       throw SurfaceMergingException(getSharedPtr(), other.getSharedPtr(),
0540                                     e.what());
0541     }
0542
0543   } else {
0544     ACTS_ERROR("DiscSurface::merge: invalid direction " << direction);
0545
0546     throw SurfaceMergingException(getSharedPtr(), other.getSharedPtr(),
0547                                   "DiscSurface::merge: invalid direction " +
0548                                       axisDirectionName(direction));
0549   }
0550 }
0551 const std::shared_ptr<const DiscBounds>& DiscSurface::boundsPtr() const {
0552   return m_bounds;
0553 }
0554
0555 void DiscSurface::assignSurfaceBounds(
0556     std::shared_ptr<const DiscBounds> newBounds) {
0557   m_bounds = std::move(newBounds);
0558 }
0559
0560 }  // namespace Acts