src/Geometry/CylinderVolumeBounds.cpp

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