src/Surfaces/CylinderBounds.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/CylinderBounds.hpp"
0010
0011 #include "Acts/Surfaces/detail/VerticesHelper.hpp"
0012 #include "Acts/Utilities/VectorHelpers.hpp"
0013 #include "Acts/Utilities/detail/periodic.hpp"
0014
0015 #include <algorithm>
0016 #include <cmath>
0017 #include <iomanip>
0018 #include <iostream>
0019 #include <numbers>
0020 #include <utility>
0021
0022 namespace Acts {
0023
0024 using VectorHelpers::perp;
0025 using VectorHelpers::phi;
0026
0027 std::vector<double> CylinderBounds::values() const {
0028   return {m_values.begin(), m_values.end()};
0029 }
0030
0031 Vector2 CylinderBounds::shifted(const Vector2& lposition) const {
0032   return {detail::radian_sym((lposition[0] / get(eR)) - get(eAveragePhi)),
0033           lposition[1]};
0034 }
0035
0036 bool CylinderBounds::inside(const Vector2& lposition) const {
0037   double bevelMinZ = get(eBevelMinZ);
0038   double bevelMaxZ = get(eBevelMaxZ);
0039
0040   double halfLengthZ = get(eHalfLengthZ);
0041   double halfPhi = get(eHalfPhiSector);
0042
0043   if (bevelMinZ == 0. || bevelMaxZ == 0.) {
0044     return detail::VerticesHelper::isInsideRectangle(
0045         shifted(lposition), Vector2(-halfPhi, -halfLengthZ),
0046         Vector2(halfPhi, halfLengthZ));
0047   }
0048
0049   double radius = get(eR);
0050   // Beleved sides will unwrap to a trapezoid
0051   ///////////////////////////////////
0052   //  ________
0053   // /| .  . |\ r/phi
0054   // \|______|/ r/phi
0055   // -Z   0  Z
0056   ///////////////////////////////////
0057   double localx =
0058       lposition[0] > radius ? 2 * radius - lposition[0] : lposition[0];
0059   Vector2 shiftedlposition = shifted(lposition);
0060   if ((std::abs(shiftedlposition[0]) <= halfPhi &&
0061        std::abs(shiftedlposition[1]) <= halfLengthZ)) {
0062     return true;
0063   }
0064
0065   if ((lposition[1] >= -(localx * std::tan(bevelMinZ) + halfLengthZ)) &&
0066       (lposition[1] <= (localx * std::tan(bevelMaxZ) + halfLengthZ))) {
0067     return true;
0068   }
0069
0070   Vector2 lowerLeft = {-radius, -halfLengthZ};
0071   Vector2 middleLeft = {0., -(halfLengthZ + radius * std::tan(bevelMinZ))};
0072   Vector2 upperLeft = {radius, -halfLengthZ};
0073   Vector2 upperRight = {radius, halfLengthZ};
0074   Vector2 middleRight = {0., (halfLengthZ + radius * std::tan(bevelMaxZ))};
0075   Vector2 lowerRight = {-radius, halfLengthZ};
0076   std::array<Vector2, 6> vertices{lowerLeft,  middleLeft,  upperLeft,
0077                                   upperRight, middleRight, lowerRight};
0078
0079   return detail::VerticesHelper::isInsidePolygon(lposition, vertices);
0080 }
0081
0082 Vector2 CylinderBounds::closestPoint(const Vector2& lposition,
0083                                      const SquareMatrix2& metric) const {
0084   double bevelMinZ = get(eBevelMinZ);
0085   double bevelMaxZ = get(eBevelMaxZ);
0086   double halfLengthZ = get(eHalfLengthZ);
0087   double radius = get(eR);
0088
0089   Vector2 lowerLeft = {-radius, -halfLengthZ};
0090   Vector2 middleLeft = {0., -(halfLengthZ + radius * std::tan(bevelMinZ))};
0091   Vector2 upperLeft = {radius, -halfLengthZ};
0092   Vector2 upperRight = {radius, halfLengthZ};
0093   Vector2 middleRight = {0., (halfLengthZ + radius * std::tan(bevelMaxZ))};
0094   Vector2 lowerRight = {-radius, halfLengthZ};
0095   Vector2 vertices[] = {lowerLeft,  middleLeft,  upperLeft,
0096                         upperRight, middleRight, lowerRight};
0097
0098   return detail::VerticesHelper::computeClosestPointOnPolygon(lposition,
0099                                                               vertices, metric);
0100 }
0101
0102 Vector2 CylinderBounds::center() const {
0103   // For cylinder bounds in local coordinates (rphi, z),
0104   // centroid is at (averagePhi, 0) since z extends symmetrically
0105   return Vector2(get(eAveragePhi), 0.0);
0106 }
0107
0108 std::ostream& CylinderBounds::toStream(std::ostream& sl) const {
0109   sl << std::setiosflags(std::ios::fixed);
0110   sl << std::setprecision(7);
0111   sl << "Acts::CylinderBounds: (radius, halfLengthZ, halfPhiSector, "
0112         "averagePhi, bevelMinZ, bevelMaxZ) = ";
0113   sl << "(" << get(eR) << ", " << get(eHalfLengthZ) << ", ";
0114   sl << get(eHalfPhiSector) << ", " << get(eAveragePhi) << ", ";
0115   sl << get(eBevelMinZ) << ", " << get(eBevelMaxZ) << ")";
0116   sl << std::setprecision(-1);
0117   return sl;
0118 }
0119
0120 std::vector<Vector3> CylinderBounds::circleVertices(
0121     const Transform3 transform, unsigned int quarterSegments) const {
0122   std::vector<Vector3> vertices;
0123
0124   double avgPhi = get(eAveragePhi);
0125   double halfPhi = get(eHalfPhiSector);
0126
0127   std::vector<double> phiRef = {};
0128   if (bool fullCylinder = coversFullAzimuth(); fullCylinder) {
0129     phiRef = {avgPhi};
0130   }
0131
0132   // Write the two bows/circles on either side
0133   std::vector<int> sides = {-1, 1};
0134   for (auto& side : sides) {
0135     // Helper method to create the segment
0136     auto svertices = detail::VerticesHelper::segmentVertices(
0137         {get(eR), get(eR)}, avgPhi - halfPhi, avgPhi + halfPhi, phiRef,
0138         quarterSegments, Vector3(0., 0., side * get(eHalfLengthZ)), transform);
0139     vertices.insert(vertices.end(), svertices.begin(), svertices.end());
0140   }
0141
0142   double bevelMinZ = get(eBevelMinZ);
0143   double bevelMaxZ = get(eBevelMaxZ);
0144
0145   // Modify the vertices position if bevel is defined
0146   if ((bevelMinZ != 0. || bevelMaxZ != 0.) && vertices.size() % 2 == 0) {
0147     auto halfWay = vertices.end() - vertices.size() / 2;
0148     double mult{1};
0149     auto invTransform = transform.inverse();
0150     auto func = [&mult, &transform, &invTransform](Vector3& v) {
0151       v = invTransform * v;
0152       v(2) += v(1) * mult;
0153       v = transform * v;
0154     };
0155     if (bevelMinZ != 0.) {
0156       mult = std::tan(-bevelMinZ);
0157       std::for_each(vertices.begin(), halfWay, func);
0158     }
0159     if (bevelMaxZ != 0.) {
0160       mult = std::tan(bevelMaxZ);
0161       std::for_each(halfWay, vertices.end(), func);
0162     }
0163   }
0164   return vertices;
0165 }
0166
0167 void CylinderBounds::checkConsistency() noexcept(false) {
0168   if (get(eR) <= 0.) {
0169     throw std::invalid_argument(
0170         "CylinderBounds: invalid radial setup: radius is negative");
0171   }
0172   if (get(eHalfLengthZ) <= 0.) {
0173     throw std::invalid_argument(
0174         "CylinderBounds: invalid length setup: half length is negative");
0175   }
0176   if (get(eHalfPhiSector) <= 0. || get(eHalfPhiSector) > std::numbers::pi) {
0177     throw std::invalid_argument("CylinderBounds: invalid phi sector setup.");
0178   }
0179   if (get(eAveragePhi) != detail::radian_sym(get(eAveragePhi)) &&
0180       std::abs(std::abs(get(eAveragePhi)) - std::numbers::pi) > s_epsilon) {
0181     throw std::invalid_argument("CylinderBounds: invalid phi positioning.");
0182   }
0183   if (get(eBevelMinZ) != detail::radian_sym(get(eBevelMinZ))) {
0184     throw std::invalid_argument("CylinderBounds: invalid bevel at min Z.");
0185   }
0186   if (get(eBevelMaxZ) != detail::radian_sym(get(eBevelMaxZ))) {
0187     throw std::invalid_argument("CylinderBounds: invalid bevel at max Z.");
0188   }
0189 }
0190
0191 }  // namespace Acts