src/Geometry/CylinderVolumeBounds.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/Geometry/CylinderVolumeBounds.hpp"
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Definitions/Direction.hpp"
0013 #include "Acts/Definitions/Tolerance.hpp"
0014 #include "Acts/Surfaces/CylinderBounds.hpp"
0015 #include "Acts/Surfaces/CylinderSurface.hpp"
0016 #include "Acts/Surfaces/DiscSurface.hpp"
0017 #include "Acts/Surfaces/PlaneSurface.hpp"
0018 #include "Acts/Surfaces/RadialBounds.hpp"
0019 #include "Acts/Surfaces/RectangleBounds.hpp"
0020 #include "Acts/Surfaces/Surface.hpp"
0021 #include "Acts/Utilities/BoundingBox.hpp"
0022 #include "Acts/Utilities/detail/periodic.hpp"
0023
0024 #include <cmath>
0025 #include <numbers>
0026 #include <utility>
0027
0028 namespace Acts {
0029
0030 CylinderVolumeBounds::CylinderVolumeBounds(double rmin, double rmax,
0031                                            double halfz, double halfphi,
0032                                            double avgphi, double bevelMinZ,
0033                                            double bevelMaxZ)
0034     : m_values() {
0035   m_values[eMinR] = rmin;
0036   m_values[eMaxR] = rmax;
0037   m_values[eHalfLengthZ] = halfz;
0038   m_values[eHalfPhiSector] = halfphi;
0039   m_values[eAveragePhi] = avgphi;
0040   m_values[eBevelMinZ] = bevelMinZ;
0041   m_values[eBevelMaxZ] = bevelMaxZ;
0042   checkConsistency();
0043   buildSurfaceBounds();
0044 }
0045
0046 CylinderVolumeBounds::CylinderVolumeBounds(
0047     const std::array<double, eSize>& values)
0048     : m_values(values) {
0049   checkConsistency();
0050   buildSurfaceBounds();
0051 }
0052
0053 CylinderVolumeBounds::CylinderVolumeBounds(const CylinderBounds& cBounds,
0054                                            double thickness)
0055     : VolumeBounds() {
0056   double cR = cBounds.get(CylinderBounds::eR);
0057   if (thickness <= 0. || (cR - 0.5 * thickness) < 0.) {
0058     throw(std::invalid_argument(
0059         "CylinderVolumeBounds: invalid extrusion thickness."));
0060   }
0061   m_values[eMinR] = cR - 0.5 * thickness;
0062   m_values[eMaxR] = cR + 0.5 * thickness;
0063   m_values[eHalfLengthZ] = cBounds.get(CylinderBounds::eHalfLengthZ);
0064   m_values[eHalfPhiSector] = cBounds.get(CylinderBounds::eHalfPhiSector);
0065   m_values[eAveragePhi] = cBounds.get(CylinderBounds::eAveragePhi);
0066   m_values[eBevelMinZ] = cBounds.get(CylinderBounds::eBevelMinZ);
0067   m_values[eBevelMaxZ] = cBounds.get(CylinderBounds::eBevelMaxZ);
0068   buildSurfaceBounds();
0069 }
0070
0071 CylinderVolumeBounds::CylinderVolumeBounds(const RadialBounds& rBounds,
0072                                            double thickness)
0073     : VolumeBounds() {
0074   if (thickness <= 0.) {
0075     throw(std::invalid_argument(
0076         "CylinderVolumeBounds: invalid extrusion thickness."));
0077   }
0078   m_values[eMinR] = rBounds.get(RadialBounds::eMinR);
0079   m_values[eMaxR] = rBounds.get(RadialBounds::eMaxR);
0080   m_values[eHalfLengthZ] = 0.5 * thickness;
0081   m_values[eHalfPhiSector] = rBounds.get(RadialBounds::eHalfPhiSector);
0082   m_values[eAveragePhi] = rBounds.get(RadialBounds::eAveragePhi);
0083   m_values[eBevelMinZ] = 0.;
0084   m_values[eBevelMaxZ] = 0.;
0085   buildSurfaceBounds();
0086 }
0087
0088 std::vector<OrientedSurface> CylinderVolumeBounds::orientedSurfaces(
0089     const Transform3& transform) const {
0090   std::vector<OrientedSurface> oSurfaces;
0091   oSurfaces.reserve(6);
0092
0093   Translation3 vMinZ(0., 0., -get(eHalfLengthZ));
0094   Translation3 vMaxZ(0., 0., get(eHalfLengthZ));
0095   // Set up transform for beveled edges if they are defined
0096   double bevelMinZ = get(eBevelMinZ);
0097   double bevelMaxZ = get(eBevelMaxZ);
0098   Transform3 transMinZ, transMaxZ;
0099   if (bevelMinZ != 0.) {
0100     double sy = 1 - 1 / std::cos(bevelMinZ);
0101     transMinZ = transform * vMinZ *
0102                 Eigen::AngleAxisd(-bevelMinZ, Eigen::Vector3d(1., 0., 0.)) *
0103                 Eigen::Scaling(1., 1. + sy, 1.);
0104   } else {
0105     transMinZ = transform * vMinZ;
0106   }
0107   if (bevelMaxZ != 0.) {
0108     double sy = 1 - 1 / std::cos(bevelMaxZ);
0109     transMaxZ = transform * vMaxZ *
0110                 Eigen::AngleAxisd(bevelMaxZ, Eigen::Vector3d(1., 0., 0.)) *
0111                 Eigen::Scaling(1., 1. + sy, 1.);
0112   } else {
0113     transMaxZ = transform * vMaxZ;
0114   }
0115   // [0] Bottom Disc (negative z)
0116   auto dSurface = Surface::makeShared<DiscSurface>(transMinZ, m_discBounds);
0117   oSurfaces.push_back(
0118       OrientedSurface{std::move(dSurface), Direction::AlongNormal()});
0119   // [1] Top Disc (positive z)
0120   dSurface = Surface::makeShared<DiscSurface>(transMaxZ, m_discBounds);
0121   oSurfaces.push_back(
0122       OrientedSurface{std::move(dSurface), Direction::OppositeNormal()});
0123
0124   // [2] Outer Cylinder
0125   auto cSurface =
0126       Surface::makeShared<CylinderSurface>(transform, m_outerCylinderBounds);
0127   oSurfaces.push_back(
0128       OrientedSurface{std::move(cSurface), Direction::OppositeNormal()});
0129
0130   // [3] Inner Cylinder (optional)
0131   if (m_innerCylinderBounds != nullptr) {
0132     cSurface =
0133         Surface::makeShared<CylinderSurface>(transform, m_innerCylinderBounds);
0134     oSurfaces.push_back(
0135         OrientedSurface{std::move(cSurface), Direction::AlongNormal()});
0136   }
0137
0138   // [4] & [5] - Sectoral planes (optional)
0139   if (m_sectorPlaneBounds != nullptr) {
0140     // sectorPlane 1 (negative phi)
0141     const Transform3 sp1Transform =
0142         Transform3(transform *
0143                    AngleAxis3(get(eAveragePhi) - get(eHalfPhiSector),
0144                               Vector3(0., 0., 1.)) *
0145                    Translation3(0.5 * (get(eMinR) + get(eMaxR)), 0., 0.) *
0146                    AngleAxis3(std::numbers::pi / 2, Vector3(1., 0., 0.)));
0147     auto pSurface =
0148         Surface::makeShared<PlaneSurface>(sp1Transform, m_sectorPlaneBounds);
0149     oSurfaces.push_back(
0150         OrientedSurface{std::move(pSurface), Direction::AlongNormal()});
0151     // sectorPlane 2 (positive phi)
0152     const Transform3 sp2Transform =
0153         Transform3(transform *
0154                    AngleAxis3(get(eAveragePhi) + get(eHalfPhiSector),
0155                               Vector3(0., 0., 1.)) *
0156                    Translation3(0.5 * (get(eMinR) + get(eMaxR)), 0., 0.) *
0157                    AngleAxis3(-std::numbers::pi / 2, Vector3(1., 0., 0.)));
0158     pSurface =
0159         Surface::makeShared<PlaneSurface>(sp2Transform, m_sectorPlaneBounds);
0160     oSurfaces.push_back(
0161         OrientedSurface{std::move(pSurface), Direction::OppositeNormal()});
0162   }
0163   return oSurfaces;
0164 }
0165
0166 void CylinderVolumeBounds::buildSurfaceBounds() {
0167   if (get(eMinR) > s_epsilon) {
0168     m_innerCylinderBounds = std::make_shared<const CylinderBounds>(
0169         get(eMinR), get(eHalfLengthZ), get(eHalfPhiSector), get(eAveragePhi),
0170         get(eBevelMinZ), get(eBevelMaxZ));
0171   }
0172   m_outerCylinderBounds = std::make_shared<const CylinderBounds>(
0173       get(eMaxR), get(eHalfLengthZ), get(eHalfPhiSector), get(eAveragePhi),
0174       get(eBevelMinZ), get(eBevelMaxZ));
0175   m_discBounds = std::make_shared<const RadialBounds>(
0176       get(eMinR), get(eMaxR), get(eHalfPhiSector), get(eAveragePhi));
0177
0178   if (std::abs(get(eHalfPhiSector) - std::numbers::pi) > s_epsilon) {
0179     m_sectorPlaneBounds = std::make_shared<const RectangleBounds>(
0180         0.5 * (get(eMaxR) - get(eMinR)), get(eHalfLengthZ));
0181   }
0182 }
0183
0184 std::ostream& CylinderVolumeBounds::toStream(std::ostream& os) const {
0185   os << std::setiosflags(std::ios::fixed);
0186   os << std::setprecision(5);
0187   os << "CylinderVolumeBounds: (rMin, rMax, halfZ, halfPhi, "
0188         "averagePhi, minBevelZ, maxBevelZ) = ";
0189   os << get(eMinR) << ", " << get(eMaxR) << ", " << get(eHalfLengthZ) << ", "
0190      << get(eHalfPhiSector) << ", " << get(eAveragePhi) << ", "
0191      << get(eBevelMinZ) << ", " << get(eBevelMaxZ);
0192   return os;
0193 }
0194
0195 Volume::BoundingBox CylinderVolumeBounds::boundingBox(
0196     const Transform3* trf, const Vector3& envelope,
0197     const Volume* entity) const {
0198   double xmax = 0, xmin = 0, ymax = 0, ymin = 0;
0199   xmax = get(eMaxR);
0200
0201   if (get(eHalfPhiSector) > std::numbers::pi / 2.) {
0202     // more than half
0203     ymax = xmax;
0204     ymin = -xmax;
0205     xmin = xmax * std::cos(get(eHalfPhiSector));
0206   } else {
0207     // less than half
0208     ymax = get(eMaxR) * std::sin(get(eHalfPhiSector));
0209     ymin = -ymax;
0210     // in this case, xmin is given by the inner radius
0211     xmin = get(eMinR) * std::cos(get(eHalfPhiSector));
0212   }
0213
0214   Vector3 vmin(xmin, ymin, -get(eHalfLengthZ));
0215   Vector3 vmax(xmax, ymax, get(eHalfLengthZ));
0216
0217   // this is probably not perfect, but at least conservative
0218   Volume::BoundingBox box{entity, vmin - envelope, vmax + envelope};
0219   return trf == nullptr ? box : box.transformed(*trf);
0220 }
0221
0222 bool CylinderVolumeBounds::inside(const Vector3& pos, double tol) const {
0223   using VectorHelpers::perp;
0224   using VectorHelpers::phi;
0225   double ros = perp(pos);
0226   bool insidePhi = cos(phi(pos)) >= cos(get(eHalfPhiSector)) - tol;
0227   bool insideR = insidePhi
0228                      ? ((ros >= get(eMinR) - tol) && (ros <= get(eMaxR) + tol))
0229                      : false;
0230   bool insideZ =
0231       insideR ? (std::abs(pos.z()) <= get(eHalfLengthZ) + tol) : false;
0232   return (insideZ && insideR && insidePhi);
0233 }
0234
0235 Vector3 CylinderVolumeBounds::referenceOffset(AxisDirection aDir)
0236     const {  // the medium radius is taken for r-type binning
0237   if (aDir == Acts::AxisDirection::AxisR ||
0238       aDir == Acts::AxisDirection::AxisRPhi) {
0239     return Vector3(0.5 * (get(eMinR) + get(eMaxR)), 0., 0.);
0240   }
0241   return VolumeBounds::referenceOffset(aDir);
0242 }
0243
0244 double CylinderVolumeBounds::referenceBorder(AxisDirection aDir) const {
0245   if (aDir == Acts::AxisDirection::AxisR) {
0246     return 0.5 * (get(eMaxR) - get(eMinR));
0247   }
0248   if (aDir == Acts::AxisDirection::AxisZ) {
0249     return get(eHalfLengthZ);
0250   }
0251   return VolumeBounds::referenceBorder(aDir);
0252 }
0253
0254 std::vector<double> CylinderVolumeBounds::values() const {
0255   return {m_values.begin(), m_values.end()};
0256 }
0257
0258 void CylinderVolumeBounds::checkConsistency() {
0259   if (get(eMinR) < 0. || get(eMaxR) <= 0.) {
0260     throw std::invalid_argument(
0261         "CylinderVolumeBounds: invalid radial input: minR (" +
0262         std::to_string(get(eMinR)) + ") < 0 or maxR (" +
0263         std::to_string(get(eMaxR)) + ") <= 0");
0264   }
0265   if (get(eMinR) >= get(eMaxR)) {
0266     throw std::invalid_argument(
0267         "CylinderVolumeBounds: invalid radial input: minR (" +
0268         std::to_string(get(eMinR)) + ") >= (" + std::to_string(get(eMaxR)) +
0269         ")");
0270   }
0271   if (get(eHalfLengthZ) <= 0) {
0272     throw std::invalid_argument(
0273         "CylinderVolumeBounds: invalid longitudinal input: hlZ (" +
0274         std::to_string(get(eHalfLengthZ)) + ") <= 0");
0275   }
0276   if (get(eHalfPhiSector) < 0. || get(eHalfPhiSector) > std::numbers::pi) {
0277     throw std::invalid_argument(
0278         "CylinderVolumeBounds: invalid phi sector setup.");
0279   }
0280   if (get(eAveragePhi) != detail::radian_sym(get(eAveragePhi))) {
0281     throw std::invalid_argument(
0282         "CylinderVolumeBounds: invalid phi positioning.");
0283   }
0284   if (get(eBevelMinZ) != detail::radian_sym(get(eBevelMinZ))) {
0285     throw std::invalid_argument("CylinderBounds: invalid bevel at min Z.");
0286   }
0287   if (get(eBevelMaxZ) != detail::radian_sym(get(eBevelMaxZ))) {
0288     throw std::invalid_argument("CylinderBounds: invalid bevel at max Z.");
0289   }
0290 }
0291
0292 void CylinderVolumeBounds::set(BoundValues bValue, double value) {
0293   set({{bValue, value}});
0294 }
0295
0296 void CylinderVolumeBounds::set(
0297     std::initializer_list<std::pair<BoundValues, double>> keyValues) {
0298   std::array<double, eSize> previous = m_values;
0299   for (const auto& [key, value] : keyValues) {
0300     m_values[key] = value;
0301   }
0302   try {
0303     checkConsistency();
0304     buildSurfaceBounds();
0305   } catch (std::invalid_argument& e) {
0306     m_values = previous;
0307     throw e;
0308   }
0309 }
0310
0311 CylinderVolumeBounds::CylinderVolumeBounds(const CylinderVolumeBounds& cylbo) =
0312     default;
0313
0314 }  // namespace Acts