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