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