File indexing completed on 2026-05-30 07:56:28
0001
0002
0003
0004
0005
0006
0007
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
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
0077 Vector3 CylinderSurface::referencePosition(const GeometryContext& gctx,
0078 AxisDirection aDir) const {
0079
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
0086
0087
0088
0089 return center(gctx);
0090 }
0091
0092
0093 RotationMatrix3 CylinderSurface::referenceFrame(
0094 const GeometryContext& gctx, const Vector3& position,
0095 const Vector3& ) const {
0096 RotationMatrix3 mFrame;
0097
0098
0099 Vector3 measY = rotSymmetryAxis(gctx);
0100
0101 Vector3 measDepth = normal(gctx, position);
0102
0103 Vector3 measX(measY.cross(measDepth).normalized());
0104
0105 mFrame.col(0) = measX;
0106 mFrame.col(1) = measY;
0107 mFrame.col(2) = measDepth;
0108
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
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
0131
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
0162 Vector3 pos3D = sfTransform.inverse() * position;
0163
0164 pos3D.z() = 0.;
0165
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
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
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
0202 double R = bounds().get(CylinderBounds::eR);
0203
0204
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
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
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
0227 auto qe = intersectionSolver(gctxTransform, position, direction);
0228
0229
0230 if (qe.solutions == 0) {
0231 return MultiIntersection3D(Intersection3D::Invalid(),
0232 Intersection3D::Invalid());
0233 }
0234
0235
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
0242 auto boundaryCheck = [&](const Vector3& solution,
0243 IntersectionStatus status) -> IntersectionStatus {
0244
0245 if (boundaryTolerance.isInfinite()) {
0246 return status;
0247 }
0248 if (boundaryTolerance.isNone() && bounds().coversFullAzimuth()) {
0249
0250
0251 const auto& tMatrix = gctxTransform.matrix();
0252
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
0264 status1 = boundaryCheck(solution1, status1);
0265
0266 Intersection3D first(solution1, qe.first, status1);
0267 if (qe.solutions == 1) {
0268 return MultiIntersection3D(first, first);
0269 }
0270
0271 Vector3 solution2 = position + qe.second * direction;
0272 IntersectionStatus status2 = std::abs(qe.second) < std::abs(tolerance)
0273 ? IntersectionStatus::onSurface
0274 : IntersectionStatus::reachable;
0275
0276 status2 = boundaryCheck(solution2, status2);
0277 Intersection3D second(solution2, qe.second, status2);
0278
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
0291 const auto pcRowVec = (position - center(gctx)).transpose().eval();
0292
0293 const auto& rotation = localToGlobalTransform(gctx).rotation();
0294
0295 const auto& localXAxis = rotation.col(0);
0296 const auto& localYAxis = rotation.col(1);
0297 const auto& localZAxis = rotation.col(2);
0298
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
0304 const auto norm = 1 / (1 - dz * dz);
0305
0306 const auto& dirRowVec = direction.transpose();
0307
0308
0309
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
0319 const auto [rotToLocalXAxis, rotToLocalYAxis, rotToLocalZAxis] =
0320 detail::rotationToLocalAxesDerivative(rotation);
0321
0322
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
0338 const auto& sTransform = localToGlobalTransform(gctx);
0339
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
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
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
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
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
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
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
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 }