Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2026-05-30 07:56:28

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