src/Geometry/ConeVolumeBounds.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/ConeVolumeBounds.hpp"
0010
0011 #include "Acts/Definitions/Direction.hpp"
0012 #include "Acts/Definitions/Tolerance.hpp"
0013 #include "Acts/Surfaces/ConeBounds.hpp"
0014 #include "Acts/Surfaces/ConeSurface.hpp"
0015 #include "Acts/Surfaces/ConvexPolygonBounds.hpp"
0016 #include "Acts/Surfaces/CylinderBounds.hpp"
0017 #include "Acts/Surfaces/CylinderSurface.hpp"
0018 #include "Acts/Surfaces/DiscSurface.hpp"
0019 #include "Acts/Surfaces/PlaneSurface.hpp"
0020 #include "Acts/Surfaces/RadialBounds.hpp"
0021 #include "Acts/Surfaces/Surface.hpp"
0022 #include "Acts/Utilities/BoundingBox.hpp"
0023 #include "Acts/Utilities/detail/OstreamStateGuard.hpp"
0024 #include "Acts/Utilities/detail/periodic.hpp"
0025
0026 #include <algorithm>
0027 #include <cmath>
0028 #include <iomanip>
0029 #include <numbers>
0030 #include <stdexcept>
0031 #include <utility>
0032
0033 namespace Acts {
0034
0035 ConeVolumeBounds::ConeVolumeBounds(double innerAlpha, double innerOffsetZ,
0036                                    double outerAlpha, double outerOffsetZ,
0037                                    double halflengthZ, double averagePhi,
0038                                    double halfPhiSector) noexcept(false)
0039     : VolumeBounds(), m_values() {
0040   m_values[eInnerAlpha] = innerAlpha;
0041   m_values[eInnerOffsetZ] = innerOffsetZ;
0042   m_values[eOuterAlpha] = outerAlpha;
0043   m_values[eOuterOffsetZ] = outerOffsetZ;
0044   m_values[eHalfLengthZ] = halflengthZ;
0045   m_values[eAveragePhi] = averagePhi;
0046   m_values[eHalfPhiSector] = halfPhiSector;
0047   buildSurfaceBounds();
0048   checkConsistency();
0049 }
0050
0051 ConeVolumeBounds::ConeVolumeBounds(double cylinderR, double alpha,
0052                                    double offsetZ, double halflengthZ,
0053                                    double averagePhi,
0054                                    double halfPhiSector) noexcept(false)
0055     : VolumeBounds(), m_values() {
0056   m_values[eInnerAlpha] = 0.;
0057   m_values[eInnerOffsetZ] = 0.;
0058   m_values[eOuterAlpha] = 0.;
0059   m_values[eOuterOffsetZ] = 0.;
0060   m_values[eHalfLengthZ] = halflengthZ;
0061   m_values[eAveragePhi] = averagePhi;
0062   m_values[eHalfPhiSector] = halfPhiSector;
0063
0064   // Cone parameters
0065   double tanAlpha = std::tan(alpha);
0066   double zmin = offsetZ - halflengthZ;
0067   double zmax = offsetZ + halflengthZ;
0068   double rmin = std::abs(zmin) * tanAlpha;
0069   double rmax = std::abs(zmax) * tanAlpha;
0070
0071   if (rmin >= cylinderR) {
0072     // Cylindrical cut-out of a cone
0073     m_innerRmin = cylinderR;
0074     m_innerRmax = cylinderR;
0075     m_outerTanAlpha = tanAlpha;
0076     m_outerRmin = rmin;
0077     m_outerRmax = rmax;
0078     m_values[eOuterAlpha] = alpha;
0079     m_values[eOuterOffsetZ] = offsetZ;
0080   } else if (rmax <= cylinderR) {
0081     // Conical cut-out of a cylinder
0082     m_outerRmin = cylinderR;
0083     m_outerRmax = cylinderR;
0084     m_innerTanAlpha = tanAlpha;
0085     m_innerRmin = rmin;
0086     m_innerRmax = rmax;
0087     m_values[eInnerAlpha] = alpha;
0088     m_values[eInnerOffsetZ] = offsetZ;
0089   } else {
0090     throw std::domain_error(
0091         "Cylinder and Cone are intersecting, not possible.");
0092   }
0093   buildSurfaceBounds();
0094   checkConsistency();
0095 }
0096
0097 std::vector<double> ConeVolumeBounds::values() const {
0098   return {m_values.begin(), m_values.end()};
0099 }
0100
0101 std::vector<Acts::OrientedSurface> Acts::ConeVolumeBounds::orientedSurfaces(
0102     const Transform3& transform) const {
0103   std::vector<OrientedSurface> oSurfaces;
0104   oSurfaces.reserve(6);
0105
0106   // Create an inner Cone
0107   if (m_innerConeBounds != nullptr) {
0108     auto innerConeTrans = transform * Translation3(0., 0., -get(eInnerOffsetZ));
0109     auto innerCone =
0110         Surface::makeShared<ConeSurface>(innerConeTrans, m_innerConeBounds);
0111     oSurfaces.emplace_back(std::move(innerCone), Direction::AlongNormal());
0112   } else if (m_innerCylinderBounds != nullptr) {
0113     // Or alternatively the inner Cylinder
0114     auto innerCylinder =
0115         Surface::makeShared<CylinderSurface>(transform, m_innerCylinderBounds);
0116     oSurfaces.emplace_back(std::move(innerCylinder), Direction::AlongNormal());
0117   }
0118
0119   // Create an outer Cone
0120   if (m_outerConeBounds != nullptr) {
0121     auto outerConeTrans = transform * Translation3(0., 0., -get(eOuterOffsetZ));
0122     auto outerCone =
0123         Surface::makeShared<ConeSurface>(outerConeTrans, m_outerConeBounds);
0124     oSurfaces.emplace_back(std::move(outerCone), Direction::OppositeNormal());
0125   } else if (m_outerCylinderBounds != nullptr) {
0126     // or alternatively an outer Cylinder
0127     auto outerCylinder =
0128         Surface::makeShared<CylinderSurface>(transform, m_outerCylinderBounds);
0129     oSurfaces.emplace_back(std::move(outerCylinder),
0130                            Direction::OppositeNormal());
0131   }
0132
0133   // Set a disc at Zmin
0134   if (m_negativeDiscBounds != nullptr) {
0135     auto negativeDiscTrans =
0136         transform * Translation3(0., 0., -get(eHalfLengthZ));
0137     auto negativeDisc = Surface::makeShared<DiscSurface>(negativeDiscTrans,
0138                                                          m_negativeDiscBounds);
0139     oSurfaces.emplace_back(std::move(negativeDisc), Direction::AlongNormal());
0140   }
0141
0142   // Set a disc at Zmax
0143   auto positiveDiscTrans = transform * Translation3(0., 0., get(eHalfLengthZ));
0144   auto positiveDisc =
0145       Surface::makeShared<DiscSurface>(positiveDiscTrans, m_positiveDiscBounds);
0146   oSurfaces.emplace_back(std::move(positiveDisc), Direction::OppositeNormal());
0147
0148   if (m_sectorBounds) {
0149     RotationMatrix3 sectorRotation;
0150     sectorRotation.col(0) = Vector3::UnitZ();
0151     sectorRotation.col(1) = Vector3::UnitX();
0152     sectorRotation.col(2) = Vector3::UnitY();
0153
0154     Transform3 negSectorRelTrans{sectorRotation};
0155     negSectorRelTrans.prerotate(
0156         AngleAxis3(get(eAveragePhi) - get(eHalfPhiSector), Vector3::UnitZ()));
0157     auto negSectorAbsTrans = transform * negSectorRelTrans;
0158     auto negSectorPlane =
0159         Surface::makeShared<PlaneSurface>(negSectorAbsTrans, m_sectorBounds);
0160     oSurfaces.emplace_back(std::move(negSectorPlane), Direction::AlongNormal());
0161
0162     Transform3 posSectorRelTrans{sectorRotation};
0163     posSectorRelTrans.prerotate(
0164         AngleAxis3(get(eAveragePhi) + get(eHalfPhiSector), Vector3::UnitZ()));
0165     auto posSectorAbsTrans = transform * posSectorRelTrans;
0166     auto posSectorPlane =
0167         Surface::makeShared<PlaneSurface>(posSectorAbsTrans, m_sectorBounds);
0168
0169     oSurfaces.emplace_back(std::move(posSectorPlane),
0170                            Direction::OppositeNormal());
0171   }
0172   return oSurfaces;
0173 }
0174
0175 void ConeVolumeBounds::checkConsistency() noexcept(false) {
0176   if (innerRmin() > outerRmin() || innerRmax() > outerRmax()) {
0177     throw std::invalid_argument("ConeVolumeBounds: invalid radial input.");
0178   }
0179   if (get(eHalfLengthZ) <= 0) {
0180     throw std::invalid_argument(
0181         "ConeVolumeBounds: invalid longitudinal input.");
0182   }
0183   if (get(eHalfPhiSector) < 0. || get(eHalfPhiSector) > std::numbers::pi) {
0184     throw std::invalid_argument("ConeVolumeBounds: invalid phi sector setup.");
0185   }
0186   if (get(eAveragePhi) != detail::radian_sym(get(eAveragePhi))) {
0187     throw std::invalid_argument("ConeVolumeBounds: invalid phi positioning.");
0188   }
0189   if (get(eInnerAlpha) == 0. && get(eOuterAlpha) == 0.) {
0190     throw std::invalid_argument(
0191         "ConeVolumeBounds: neither inner nor outer cone.");
0192   }
0193 }
0194
0195 bool ConeVolumeBounds::inside(const Vector3& pos, double tol) const {
0196   double z = pos.z();
0197   double zmin = z + tol;
0198   double zmax = z - tol;
0199   // Quick check outside z
0200   if (zmin < -get(eHalfLengthZ) || zmax > get(eHalfLengthZ)) {
0201     return false;
0202   }
0203   double r = VectorHelpers::perp(pos);
0204   if (std::abs(get(eHalfPhiSector) - std::numbers::pi) > s_onSurfaceTolerance) {
0205     // need to check the phi sector - approximate phi tolerance
0206     double phitol = tol / r;
0207     double phi = VectorHelpers::phi(pos);
0208     double phimin = phi - phitol;
0209     double phimax = phi + phitol;
0210     if (phimin < get(eAveragePhi) - get(eHalfPhiSector) ||
0211         phimax > get(eAveragePhi) + get(eHalfPhiSector)) {
0212       return false;
0213     }
0214   }
0215   // We are within phi sector check box r quickly
0216   double rmin = r + tol;
0217   double rmax = r - tol;
0218   if (rmin > innerRmax() && rmax < outerRmin()) {
0219     return true;
0220   }
0221   // Finally we need to check the cone
0222   if (m_innerConeBounds != nullptr) {
0223     double innerConeR = m_innerConeBounds->r(std::abs(z + get(eInnerOffsetZ)));
0224     if (innerConeR > rmin) {
0225       return false;
0226     }
0227   } else if (innerRmax() > rmin) {
0228     return false;
0229   }
0230   // And the outer cone
0231   if (m_outerConeBounds != nullptr) {
0232     double outerConeR = m_outerConeBounds->r(std::abs(z + get(eOuterOffsetZ)));
0233     if (outerConeR < rmax) {
0234       return false;
0235     }
0236   } else if (outerRmax() < rmax) {
0237     return false;
0238   }
0239   return true;
0240 }
0241
0242 void ConeVolumeBounds::buildSurfaceBounds() {
0243   // Build inner cone or inner cylinder
0244   if (get(eInnerAlpha) > s_epsilon) {
0245     m_innerTanAlpha = std::tan(get(eInnerAlpha));
0246     double innerZmin = get(eInnerOffsetZ) - get(eHalfLengthZ);
0247     double innerZmax = get(eInnerOffsetZ) + get(eHalfLengthZ);
0248     m_innerRmin = std::abs(innerZmin) * m_innerTanAlpha;
0249     m_innerRmax = std::abs(innerZmax) * m_innerTanAlpha;
0250     m_innerConeBounds =
0251         std::make_shared<ConeBounds>(get(eInnerAlpha), innerZmin, innerZmax,
0252                                      get(eHalfPhiSector), get(eAveragePhi));
0253   } else if (m_innerRmin == m_innerRmax && m_innerRmin > s_epsilon) {
0254     m_innerCylinderBounds = std::make_shared<CylinderBounds>(
0255         m_innerRmin, get(eHalfLengthZ), get(eHalfPhiSector), get(eAveragePhi));
0256   }
0257
0258   if (get(eOuterAlpha) > s_epsilon) {
0259     m_outerTanAlpha = std::tan(get(eOuterAlpha));
0260     double outerZmin = get(eOuterOffsetZ) - get(eHalfLengthZ);
0261     double outerZmax = get(eOuterOffsetZ) + get(eHalfLengthZ);
0262     m_outerRmin = std::abs(outerZmin) * m_outerTanAlpha;
0263     m_outerRmax = std::abs(outerZmax) * m_outerTanAlpha;
0264     m_outerConeBounds =
0265         std::make_shared<ConeBounds>(get(eOuterAlpha), outerZmin, outerZmax,
0266                                      get(eHalfPhiSector), get(eAveragePhi));
0267
0268   } else if (m_outerRmin == m_outerRmax) {
0269     m_outerCylinderBounds = std::make_shared<CylinderBounds>(
0270         m_outerRmax, get(eHalfLengthZ), get(eHalfPhiSector), get(eAveragePhi));
0271   }
0272
0273   if (get(eHalfLengthZ) < std::max(get(eInnerOffsetZ), get(eOuterOffsetZ))) {
0274     m_negativeDiscBounds = std::make_shared<RadialBounds>(
0275         m_innerRmin, m_outerRmin, get(eHalfPhiSector), get(eAveragePhi));
0276   }
0277
0278   m_positiveDiscBounds = std::make_shared<RadialBounds>(
0279       m_innerRmax, m_outerRmax, get(eHalfPhiSector), get(eAveragePhi));
0280
0281   // Create the sector bounds
0282   if (std::abs(get(eHalfPhiSector) - std::numbers::pi) > s_epsilon) {
0283     // The 4 points building the sector
0284     std::vector<Vector2> polyVertices = {{-get(eHalfLengthZ), m_innerRmin},
0285                                          {get(eHalfLengthZ), m_innerRmax},
0286                                          {get(eHalfLengthZ), m_outerRmax},
0287                                          {-get(eHalfLengthZ), m_outerRmin}};
0288     m_sectorBounds =
0289         std::make_shared<ConvexPolygonBounds<4>>(std::move(polyVertices));
0290   }
0291 }
0292
0293 std::ostream& ConeVolumeBounds::toStream(std::ostream& os) const {
0294   detail::OstreamStateGuard guard{os};
0295   os << std::fixed << std::setprecision(5);
0296   os << "Acts::ConeVolumeBounds : (innerAlpha, innerOffsetZ, outerAlpha,";
0297   os << "  outerOffsetZ, halflenghZ, averagePhi, halfPhiSector) = ";
0298   os << get(eInnerAlpha) << ", " << get(eInnerOffsetZ) << ", ";
0299   os << get(eOuterAlpha) << ", " << get(eOuterOffsetZ) << ", ";
0300   os << get(eHalfLengthZ) << ", " << get(eAveragePhi) << std::endl;
0301   return os;
0302 }
0303
0304 Volume::BoundingBox ConeVolumeBounds::boundingBox(const Transform3* trf,
0305                                                   const Vector3& envelope,
0306                                                   const Volume* entity) const {
0307   Vector3 vmin(-outerRmax(), -outerRmax(), -0.5 * get(eHalfLengthZ));
0308   Vector3 vmax(outerRmax(), outerRmax(), 0.5 * get(eHalfLengthZ));
0309   Volume::BoundingBox box(entity, vmin - envelope, vmax + envelope);
0310   return trf == nullptr ? box : box.transformed(*trf);
0311 }
0312
0313 double ConeVolumeBounds::innerRmin() const {
0314   return m_innerRmin;
0315 }
0316
0317 double ConeVolumeBounds::innerRmax() const {
0318   return m_innerRmax;
0319 }
0320
0321 double ConeVolumeBounds::innerTanAlpha() const {
0322   return m_innerTanAlpha;
0323 }
0324
0325 double ConeVolumeBounds::outerRmin() const {
0326   return m_outerRmin;
0327 }
0328
0329 double ConeVolumeBounds::outerRmax() const {
0330   return m_outerRmax;
0331 }
0332
0333 double ConeVolumeBounds::outerTanAlpha() const {
0334   return m_outerTanAlpha;
0335 }
0336
0337 }  // namespace Acts