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