src/Geometry/ConeVolumeBounds.cpp

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