src/Surfaces/CylinderSurface.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/CylinderSurface.hpp"
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Definitions/Tolerance.hpp"
0013 #include "Acts/Definitions/Units.hpp"
0014 #include "Acts/Geometry/GeometryObject.hpp"
0015 #include "Acts/Surfaces/CylinderBounds.hpp"
0016 #include "Acts/Surfaces/SurfaceError.hpp"
0017 #include "Acts/Surfaces/SurfaceMergingException.hpp"
0018 #include "Acts/Surfaces/detail/AlignmentHelper.hpp"
0019 #include "Acts/Surfaces/detail/FacesHelper.hpp"
0020 #include "Acts/Surfaces/detail/MergeHelper.hpp"
0021 #include "Acts/Utilities/Intersection.hpp"
0022 #include "Acts/Utilities/ThrowAssert.hpp"
0023 #include "Acts/Utilities/detail/periodic.hpp"
0024
0025 #include <algorithm>
0026 #include <cassert>
0027 #include <cmath>
0028 #include <iostream>
0029 #include <memory>
0030 #include <stdexcept>
0031 #include <utility>
0032 #include <vector>
0033
0034 namespace Acts {
0035
0036 using VectorHelpers::perp;
0037 using VectorHelpers::phi;
0038
0039 CylinderSurface::CylinderSurface(const CylinderSurface& other)
0040     : GeometryObject(), RegularSurface(other), m_bounds(other.m_bounds) {}
0041
0042 CylinderSurface::CylinderSurface(const GeometryContext& gctx,
0043                                  const CylinderSurface& other,
0044                                  const Transform3& shift)
0045     : GeometryObject(),
0046       RegularSurface(gctx, other, shift),
0047       m_bounds(other.m_bounds) {}
0048
0049 CylinderSurface::CylinderSurface(const Transform3& transform, double radius,
0050                                  double halfz, double halfphi, double avphi,
0051                                  double bevelMinZ, double bevelMaxZ)
0052     : RegularSurface(transform),
0053       m_bounds(std::make_shared<const CylinderBounds>(
0054           radius, halfz, halfphi, avphi, bevelMinZ, bevelMaxZ)) {}
0055
0056 CylinderSurface::CylinderSurface(std::shared_ptr<const CylinderBounds> cbounds,
0057                                  const DetectorElementBase& detelement)
0058     : RegularSurface(detelement), m_bounds(std::move(cbounds)) {
0059   // surfaces representing a detector element must have bounds
0060   throw_assert(m_bounds, "CylinderBounds must not be nullptr");
0061 }
0062
0063 CylinderSurface::CylinderSurface(const Transform3& transform,
0064                                  std::shared_ptr<const CylinderBounds> cbounds)
0065     : RegularSurface(transform), m_bounds(std::move(cbounds)) {
0066   throw_assert(m_bounds, "CylinderBounds must not be nullptr");
0067 }
0068
0069 CylinderSurface& CylinderSurface::operator=(const CylinderSurface& other) {
0070   if (this != &other) {
0071     Surface::operator=(other);
0072     m_bounds = other.m_bounds;
0073   }
0074   return *this;
0075 }
0076
0077 // return the binning position for ordering in the BinnedArray
0078 Vector3 CylinderSurface::referencePosition(const GeometryContext& gctx,
0079                                            AxisDirection aDir) const {
0080   // special binning type for R-type methods
0081   if (aDir == AxisDirection::AxisR || aDir == AxisDirection::AxisRPhi) {
0082     double R = bounds().get(CylinderBounds::eR);
0083     double phi = bounds().get(CylinderBounds::eAveragePhi);
0084     return localToGlobal(gctx, Vector2{phi * R, 0});
0085   }
0086   // give the center as default for all of these binning types
0087   // AxisDirection::AxisX, AxisDirection::AxisY, AxisDirection::AxisZ,
0088   // AxisDirection::AxisR, AxisDirection::AxisPhi, AxisDirection::AxisRPhi,
0089   // AxisDirection::AxisTheta, AxisDirection::AxisEta
0090   return center(gctx);
0091 }
0092
0093 // return the measurement frame: it's the tangential plane
0094 RotationMatrix3 CylinderSurface::referenceFrame(
0095     const GeometryContext& gctx, const Vector3& position,
0096     const Vector3& /*direction*/) const {
0097   RotationMatrix3 mFrame;
0098   // construct the measurement frame
0099   // measured Y is the z axis
0100   Vector3 measY = rotSymmetryAxis(gctx);
0101   // measured z is the position normalized transverse (in local)
0102   Vector3 measDepth = normal(gctx, position);
0103   // measured X is what comoes out of it
0104   Vector3 measX(measY.cross(measDepth).normalized());
0105   // assign the columnes
0106   mFrame.col(0) = measX;
0107   mFrame.col(1) = measY;
0108   mFrame.col(2) = measDepth;
0109   // return the rotation matrix
0110   return mFrame;
0111 }
0112
0113 Surface::SurfaceType CylinderSurface::type() const {
0114   return Surface::Cylinder;
0115 }
0116
0117 Vector3 CylinderSurface::localToGlobal(const GeometryContext& gctx,
0118                                        const Vector2& lposition) const {
0119   // create the position in the local 3d frame
0120   double r = bounds().get(CylinderBounds::eR);
0121   double phi = lposition[0] / r;
0122   Vector3 position(r * cos(phi), r * sin(phi), lposition[1]);
0123   return transform(gctx) * position;
0124 }
0125
0126 Result<Vector2> CylinderSurface::globalToLocal(const GeometryContext& gctx,
0127                                                const Vector3& position,
0128                                                double tolerance) const {
0129   double inttol = tolerance;
0130   if (tolerance == s_onSurfaceTolerance) {
0131     // transform default value!
0132     // @TODO: check if s_onSurfaceTolerance would do here
0133     inttol = bounds().get(CylinderBounds::eR) * 0.0001;
0134   }
0135   if (inttol < 0.01) {
0136     inttol = 0.01;
0137   }
0138   const Transform3& sfTransform = transform(gctx);
0139   Transform3 inverseTrans(sfTransform.inverse());
0140   Vector3 loc3Dframe(inverseTrans * position);
0141   if (std::abs(perp(loc3Dframe) - bounds().get(CylinderBounds::eR)) > inttol) {
0142     return Result<Vector2>::failure(SurfaceError::GlobalPositionNotOnSurface);
0143   }
0144   return Result<Vector2>::success(
0145       {bounds().get(CylinderBounds::eR) * phi(loc3Dframe), loc3Dframe.z()});
0146 }
0147
0148 std::string CylinderSurface::name() const {
0149   return "Acts::CylinderSurface";
0150 }
0151
0152 Vector3 CylinderSurface::normal(const GeometryContext& gctx,
0153                                 const Vector2& lposition) const {
0154   double phi = lposition[0] / m_bounds->get(CylinderBounds::eR);
0155   Vector3 localNormal(cos(phi), sin(phi), 0.);
0156   return transform(gctx).linear() * localNormal;
0157 }
0158
0159 Vector3 CylinderSurface::normal(const GeometryContext& gctx,
0160                                 const Vector3& position) const {
0161   const Transform3& sfTransform = transform(gctx);
0162   // get it into the cylinder frame
0163   Vector3 pos3D = sfTransform.inverse() * position;
0164   // set the z coordinate to 0
0165   pos3D.z() = 0.;
0166   // normalize and rotate back into global
0167   return sfTransform.linear() * pos3D.normalized();
0168 }
0169
0170 double CylinderSurface::pathCorrection(const GeometryContext& gctx,
0171                                        const Vector3& position,
0172                                        const Vector3& direction) const {
0173   Vector3 normalT = normal(gctx, position);
0174   double cosAlpha = normalT.dot(direction);
0175   return std::abs(1. / cosAlpha);
0176 }
0177
0178 const CylinderBounds& CylinderSurface::bounds() const {
0179   return *m_bounds;
0180 }
0181
0182 Polyhedron CylinderSurface::polyhedronRepresentation(
0183     const GeometryContext& gctx, unsigned int quarterSegments) const {
0184   auto ctrans = transform(gctx);
0185
0186   // Prepare vertices and faces
0187   std::vector<Vector3> vertices =
0188       bounds().circleVertices(ctrans, quarterSegments);
0189   auto [faces, triangularMesh] =
0190       detail::FacesHelper::cylindricalFaceMesh(vertices);
0191   return Polyhedron(vertices, faces, triangularMesh, false);
0192 }
0193
0194 Vector3 CylinderSurface::rotSymmetryAxis(const GeometryContext& gctx) const {
0195   // fast access via transform matrix (and not rotation())
0196   return transform(gctx).matrix().block<3, 1>(0, 2);
0197 }
0198
0199 detail::RealQuadraticEquation CylinderSurface::intersectionSolver(
0200     const Transform3& transform, const Vector3& position,
0201     const Vector3& direction) const {
0202   // Solve for radius R
0203   double R = bounds().get(CylinderBounds::eR);
0204
0205   // Get the transformation matrtix
0206   const auto& tMatrix = transform.matrix();
0207   Vector3 caxis = tMatrix.block<3, 1>(0, 2).transpose();
0208   Vector3 ccenter = tMatrix.block<3, 1>(0, 3).transpose();
0209
0210   // Check documentation for explanation
0211   Vector3 pc = position - ccenter;
0212   Vector3 pcXcd = pc.cross(caxis);
0213   Vector3 ldXcd = direction.cross(caxis);
0214   double a = ldXcd.dot(ldXcd);
0215   double b = 2. * (ldXcd.dot(pcXcd));
0216   double c = pcXcd.dot(pcXcd) - (R * R);
0217   // And solve the qaudratic equation
0218   return detail::RealQuadraticEquation(a, b, c);
0219 }
0220
0221 SurfaceMultiIntersection CylinderSurface::intersect(
0222     const GeometryContext& gctx, const Vector3& position,
0223     const Vector3& direction, const BoundaryTolerance& boundaryTolerance,
0224     double tolerance) const {
0225   const auto& gctxTransform = transform(gctx);
0226
0227   // Solve the quadratic equation
0228   auto qe = intersectionSolver(gctxTransform, position, direction);
0229
0230   // If no valid solution return a non-valid surfaceIntersection
0231   if (qe.solutions == 0) {
0232     return {{Intersection3D::invalid(), Intersection3D::invalid()},
0233             *this,
0234             boundaryTolerance};
0235   }
0236
0237   // Check the validity of the first solution
0238   Vector3 solution1 = position + qe.first * direction;
0239   IntersectionStatus status1 = std::abs(qe.first) < std::abs(tolerance)
0240                                    ? IntersectionStatus::onSurface
0241                                    : IntersectionStatus::reachable;
0242
0243   // Helper method for boundary check
0244   auto boundaryCheck = [&](const Vector3& solution,
0245                            IntersectionStatus status) -> IntersectionStatus {
0246     // No check to be done, return current status
0247     if (boundaryTolerance.isInfinite()) {
0248       return status;
0249     }
0250     if (boundaryTolerance.isNone() && bounds().coversFullAzimuth()) {
0251       // Project out the current Z value via local z axis
0252       // Built-in local to global for speed reasons
0253       const auto& tMatrix = gctxTransform.matrix();
0254       // Create the reference vector in local
0255       const Vector3 vecLocal(solution - tMatrix.block<3, 1>(0, 3));
0256       double cZ = vecLocal.dot(tMatrix.block<3, 1>(0, 2));
0257       double hZ = bounds().get(CylinderBounds::eHalfLengthZ) + tolerance;
0258       return std::abs(cZ) < std::abs(hZ) ? status
0259                                          : IntersectionStatus::unreachable;
0260     }
0261     return isOnSurface(gctx, solution, direction, boundaryTolerance)
0262                ? status
0263                : IntersectionStatus::unreachable;
0264   };
0265   // Check first solution for boundary compatibility
0266   status1 = boundaryCheck(solution1, status1);
0267   // Set the intersection
0268   Intersection3D first(solution1, qe.first, status1);
0269   if (qe.solutions == 1) {
0270     return {{first, first}, *this, boundaryTolerance};
0271   }
0272   // Check the validity of the second solution
0273   Vector3 solution2 = position + qe.second * direction;
0274   IntersectionStatus status2 = std::abs(qe.second) < std::abs(tolerance)
0275                                    ? IntersectionStatus::onSurface
0276                                    : IntersectionStatus::reachable;
0277   // Check first solution for boundary compatibility
0278   status2 = boundaryCheck(solution2, status2);
0279   Intersection3D second(solution2, qe.second, status2);
0280   // Order based on path length
0281   if (first.pathLength() <= second.pathLength()) {
0282     return {{first, second}, *this, boundaryTolerance};
0283   }
0284   return {{second, first}, *this, boundaryTolerance};
0285 }
0286
0287 AlignmentToPathMatrix CylinderSurface::alignmentToPathDerivative(
0288     const GeometryContext& gctx, const Vector3& position,
0289     const Vector3& direction) const {
0290   assert(isOnSurface(gctx, position, direction, BoundaryTolerance::Infinite()));
0291
0292   // The vector between position and center
0293   const auto pcRowVec = (position - center(gctx)).transpose().eval();
0294   // The rotation
0295   const auto& rotation = transform(gctx).rotation();
0296   // The local frame x/y/z axis
0297   const auto& localXAxis = rotation.col(0);
0298   const auto& localYAxis = rotation.col(1);
0299   const auto& localZAxis = rotation.col(2);
0300   // The local coordinates
0301   const auto localPos = (rotation.transpose() * position).eval();
0302   const auto dx = direction.dot(localXAxis);
0303   const auto dy = direction.dot(localYAxis);
0304   const auto dz = direction.dot(localZAxis);
0305   // The normalization factor
0306   const auto norm = 1 / (1 - dz * dz);
0307   // The direction transpose
0308   const auto& dirRowVec = direction.transpose();
0309   // The derivative of path w.r.t. the local axes
0310   // @note The following calculations assume that the intersection of the track
0311   // with the cylinder always satisfy: perp(localPos) = R
0312   const auto localXAxisToPath =
0313       (-2 * norm * (dx * pcRowVec + localPos.x() * dirRowVec)).eval();
0314   const auto localYAxisToPath =
0315       (-2 * norm * (dy * pcRowVec + localPos.y() * dirRowVec)).eval();
0316   const auto localZAxisToPath =
0317       (-4 * norm * norm * (dx * localPos.x() + dy * localPos.y()) * dz *
0318        dirRowVec)
0319           .eval();
0320   // Calculate the derivative of local frame axes w.r.t its rotation
0321   const auto [rotToLocalXAxis, rotToLocalYAxis, rotToLocalZAxis] =
0322       detail::rotationToLocalAxesDerivative(rotation);
0323   // Initialize the derivative of propagation path w.r.t. local frame
0324   // translation (origin) and rotation
0325   AlignmentToPathMatrix alignToPath = AlignmentToPathMatrix::Zero();
0326   alignToPath.segment<3>(eAlignmentCenter0) =
0327       2 * norm * (dx * localXAxis.transpose() + dy * localYAxis.transpose());
0328   alignToPath.segment<3>(eAlignmentRotation0) =
0329       localXAxisToPath * rotToLocalXAxis + localYAxisToPath * rotToLocalYAxis +
0330       localZAxisToPath * rotToLocalZAxis;
0331
0332   return alignToPath;
0333 }
0334
0335 ActsMatrix<2, 3> CylinderSurface::localCartesianToBoundLocalDerivative(
0336     const GeometryContext& gctx, const Vector3& position) const {
0337   using VectorHelpers::perp;
0338   using VectorHelpers::phi;
0339   // The local frame transform
0340   const auto& sTransform = transform(gctx);
0341   // calculate the transformation to local coordinates
0342   const Vector3 localPos = sTransform.inverse() * position;
0343   const double lr = perp(localPos);
0344   const double lphi = phi(localPos);
0345   const double lcphi = std::cos(lphi);
0346   const double lsphi = std::sin(lphi);
0347   // Solve for radius R
0348   double R = bounds().get(CylinderBounds::eR);
0349   ActsMatrix<2, 3> loc3DToLocBound = ActsMatrix<2, 3>::Zero();
0350   loc3DToLocBound << -R * lsphi / lr, R * lcphi / lr, 0, 0, 0, 1;
0351
0352   return loc3DToLocBound;
0353 }
0354
0355 std::pair<std::shared_ptr<CylinderSurface>, bool> CylinderSurface::mergedWith(
0356     const CylinderSurface& other, AxisDirection direction,
0357     bool externalRotation, const Logger& logger) const {
0358   using namespace UnitLiterals;
0359
0360   ACTS_VERBOSE("Merging cylinder surfaces in " << axisDirectionName(direction)
0361                                                << " direction");
0362
0363   if (m_associatedDetElement != nullptr ||
0364       other.m_associatedDetElement != nullptr) {
0365     throw SurfaceMergingException(getSharedPtr(), other.getSharedPtr(),
0366                                   "CylinderSurface::merge: surfaces are "
0367                                   "associated with a detector element");
0368   }
0369
0370   assert(m_transform != nullptr && other.m_transform != nullptr);
0371
0372   Transform3 otherLocal = m_transform->inverse() * *other.m_transform;
0373
0374   constexpr auto tolerance = s_onSurfaceTolerance;
0375
0376   // surface cannot have any relative rotation
0377
0378   if (std::abs(otherLocal.linear().col(eX)[eZ]) >= tolerance ||
0379       std::abs(otherLocal.linear().col(eY)[eZ]) >= tolerance) {
0380     ACTS_ERROR("CylinderSurface::merge: surfaces have relative rotation");
0381     throw SurfaceMergingException(
0382         getSharedPtr(), other.getSharedPtr(),
0383         "CylinderSurface::merge: surfaces have relative rotation");
0384   }
0385
0386   auto checkNoBevel = [this, &logger, &other](const auto& bounds) {
0387     if (bounds.get(CylinderBounds::eBevelMinZ) != 0.0) {
0388       ACTS_ERROR(
0389           "CylinderVolumeStack requires all volumes to have a bevel angle of "
0390           "0");
0391       throw SurfaceMergingException(
0392           getSharedPtr(), other.getSharedPtr(),
0393           "CylinderVolumeStack requires all volumes to have a bevel angle of "
0394           "0");
0395     }
0396
0397     if (bounds.get(CylinderBounds::eBevelMaxZ) != 0.0) {
0398       ACTS_ERROR(
0399           "CylinderVolumeStack requires all volumes to have a bevel angle of "
0400           "0");
0401       throw SurfaceMergingException(
0402           getSharedPtr(), other.getSharedPtr(),
0403           "CylinderVolumeStack requires all volumes to have a bevel angle of "
0404           "0");
0405     }
0406   };
0407
0408   checkNoBevel(bounds());
0409   checkNoBevel(other.bounds());
0410
0411   // radii need to be identical
0412   if (std::abs(bounds().get(CylinderBounds::eR) -
0413                other.bounds().get(CylinderBounds::eR)) > tolerance) {
0414     ACTS_ERROR("CylinderSurface::merge: surfaces have different radii");
0415     throw SurfaceMergingException(
0416         getSharedPtr(), other.getSharedPtr(),
0417         "CylinderSurface::merge: surfaces have different radii");
0418   }
0419
0420   double r = bounds().get(CylinderBounds::eR);
0421
0422   // no translation in x/z is allowed
0423   Vector3 translation = otherLocal.translation();
0424
0425   if (std::abs(translation[0]) > tolerance ||
0426       std::abs(translation[1]) > tolerance) {
0427     ACTS_ERROR(
0428         "CylinderSurface::merge: surfaces have relative translation in x/y");
0429     throw SurfaceMergingException(
0430         getSharedPtr(), other.getSharedPtr(),
0431         "CylinderSurface::merge: surfaces have relative translation in x/y");
0432   }
0433
0434   double hlZ = bounds().get(CylinderBounds::eHalfLengthZ);
0435   double minZ = -hlZ;
0436   double maxZ = hlZ;
0437
0438   double zShift = translation[2];
0439   double otherHlZ = other.bounds().get(CylinderBounds::eHalfLengthZ);
0440   double otherMinZ = -otherHlZ + zShift;
0441   double otherMaxZ = otherHlZ + zShift;
0442
0443   double hlPhi = bounds().get(CylinderBounds::eHalfPhiSector);
0444   double avgPhi = bounds().get(CylinderBounds::eAveragePhi);
0445
0446   double otherHlPhi = other.bounds().get(CylinderBounds::eHalfPhiSector);
0447   double otherAvgPhi = other.bounds().get(CylinderBounds::eAveragePhi);
0448
0449   if (direction == AxisDirection::AxisZ) {
0450     // z shift must match the bounds
0451
0452     if (std::abs(otherLocal.linear().col(eY)[eX]) >= tolerance &&
0453         (!bounds().coversFullAzimuth() ||
0454          !other.bounds().coversFullAzimuth())) {
0455       throw SurfaceMergingException(getSharedPtr(), other.getSharedPtr(),
0456                                     "CylinderSurface::merge: surfaces have "
0457                                     "relative rotation in z and phi sector");
0458     }
0459
0460     ACTS_VERBOSE("this: [" << minZ << ", " << maxZ << "] ~> "
0461                            << (minZ + maxZ) / 2.0 << " +- " << hlZ);
0462     ACTS_VERBOSE("zShift: " << zShift);
0463
0464     ACTS_VERBOSE("other: [" << otherMinZ << ", " << otherMaxZ << "] ~> "
0465                             << (otherMinZ + otherMaxZ) / 2.0 << " +- "
0466                             << otherHlZ);
0467     if (std::abs(maxZ - otherMinZ) > tolerance &&
0468         std::abs(minZ - otherMaxZ) > tolerance) {
0469       ACTS_ERROR("CylinderSurface::merge: surfaces have incompatible z bounds");
0470       throw SurfaceMergingException(
0471           getSharedPtr(), other.getSharedPtr(),
0472           "CylinderSurface::merge: surfaces have incompatible z bounds");
0473     }
0474
0475     if (hlPhi != otherHlPhi || avgPhi != otherAvgPhi) {
0476       throw SurfaceMergingException(getSharedPtr(), other.getSharedPtr(),
0477                                     "CylinderSurface::merge: surfaces have "
0478                                     "different phi sectors");
0479     }
0480
0481     double newMaxZ = std::max(maxZ, otherMaxZ);
0482     double newMinZ = std::min(minZ, otherMinZ);
0483     double newHlZ = (newMaxZ - newMinZ) / 2.0;
0484     double newMidZ = (newMaxZ + newMinZ) / 2.0;
0485     ACTS_VERBOSE("merged: [" << newMinZ << ", " << newMaxZ << "] ~> " << newMidZ
0486                              << " +- " << newHlZ);
0487
0488     auto newBounds = std::make_shared<CylinderBounds>(r, newHlZ, hlPhi, avgPhi);
0489
0490     Transform3 newTransform =
0491         *m_transform * Translation3{Vector3::UnitZ() * newMidZ};
0492
0493     return {Surface::makeShared<CylinderSurface>(newTransform, newBounds),
0494             zShift < 0};
0495
0496   } else if (direction == AxisDirection::AxisRPhi) {
0497     // no z shift is allowed
0498     if (std::abs(translation[2]) > tolerance) {
0499       ACTS_ERROR(
0500           "CylinderSurface::merge: surfaces have relative translation in z for "
0501           "rPhi merging");
0502       throw SurfaceMergingException(
0503           getSharedPtr(), other.getSharedPtr(),
0504           "CylinderSurface::merge: surfaces have relative translation in z for "
0505           "rPhi merging");
0506     }
0507
0508     if (hlZ != otherHlZ) {
0509       throw SurfaceMergingException(getSharedPtr(), other.getSharedPtr(),
0510                                     "CylinderSurface::merge: surfaces have "
0511                                     "different z bounds");
0512     }
0513
0514     // Figure out signed relative rotation
0515     Vector2 rotatedX = otherLocal.linear().col(eX).head<2>();
0516     double zrotation = std::atan2(rotatedX[1], rotatedX[0]);
0517
0518     ACTS_VERBOSE("this:  [" << avgPhi / 1_degree << " +- " << hlPhi / 1_degree
0519                             << "]");
0520     ACTS_VERBOSE("other: [" << otherAvgPhi / 1_degree << " +- "
0521                             << otherHlPhi / 1_degree << "]");
0522
0523     ACTS_VERBOSE("Relative rotation around local z: " << zrotation / 1_degree);
0524
0525     double prevOtherAvgPhi = otherAvgPhi;
0526     otherAvgPhi = detail::radian_sym(otherAvgPhi + zrotation);
0527     ACTS_VERBOSE("~> local other average phi: "
0528                  << otherAvgPhi / 1_degree
0529                  << " (was: " << prevOtherAvgPhi / 1_degree << ")");
0530
0531     try {
0532       auto [newHlPhi, newAvgPhi, reversed] = detail::mergedPhiSector(
0533           hlPhi, avgPhi, otherHlPhi, otherAvgPhi, logger, tolerance);
0534
0535       Transform3 newTransform = *m_transform;
0536
0537       if (externalRotation) {
0538         ACTS_VERBOSE("Modifying transform for external rotation of "
0539                      << newAvgPhi / 1_degree);
0540         newTransform = newTransform * AngleAxis3(newAvgPhi, Vector3::UnitZ());
0541         newAvgPhi = 0.;
0542       }
0543
0544       auto newBounds = std::make_shared<CylinderBounds>(
0545           r, bounds().get(CylinderBounds::eHalfLengthZ), newHlPhi, newAvgPhi);
0546
0547       return {Surface::makeShared<CylinderSurface>(newTransform, newBounds),
0548               reversed};
0549     } catch (const std::invalid_argument& e) {
0550       throw SurfaceMergingException(getSharedPtr(), other.getSharedPtr(),
0551                                     e.what());
0552     }
0553   } else {
0554     throw SurfaceMergingException(getSharedPtr(), other.getSharedPtr(),
0555                                   "CylinderSurface::merge: invalid direction " +
0556                                       axisDirectionName(direction));
0557   }
0558 }
0559
0560 }  // namespace Acts