src/Geometry/TrapezoidVolumeBounds.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/TrapezoidVolumeBounds.hpp"
0010
0011 #include "Acts/Definitions/Direction.hpp"
0012 #include "Acts/Surfaces/BoundaryTolerance.hpp"
0013 #include "Acts/Surfaces/PlaneSurface.hpp"
0014 #include "Acts/Surfaces/RectangleBounds.hpp"
0015 #include "Acts/Surfaces/Surface.hpp"
0016 #include "Acts/Surfaces/TrapezoidBounds.hpp"
0017 #include "Acts/Utilities/BoundingBox.hpp"
0018 #include "Acts/Utilities/detail/OstreamStateGuard.hpp"
0019
0020 #include <cmath>
0021 #include <cstddef>
0022 #include <iomanip>
0023 #include <numbers>
0024 #include <utility>
0025
0026 namespace Acts {
0027
0028 TrapezoidVolumeBounds::TrapezoidVolumeBounds(double minhalex, double maxhalex,
0029                                              double haley, double halez)
0030     : VolumeBounds() {
0031   m_values[eHalfLengthXnegY] = minhalex;
0032   m_values[eHalfLengthXposY] = maxhalex;
0033   m_values[eHalfLengthY] = haley;
0034   m_values[eHalfLengthZ] = halez;
0035   m_values[eAlpha] = std::atan2(2 * haley, (maxhalex - minhalex));
0036   m_values[eBeta] = std::numbers::pi - get(eAlpha);
0037   checkConsistency();
0038   buildSurfaceBounds();
0039 }
0040
0041 TrapezoidVolumeBounds::TrapezoidVolumeBounds(double minhalex, double haley,
0042                                              double halez, double alpha,
0043                                              double beta)
0044     : VolumeBounds() {
0045   m_values[eHalfLengthXnegY] = minhalex;
0046   m_values[eHalfLengthY] = haley;
0047   m_values[eHalfLengthZ] = halez;
0048   m_values[eAlpha] = alpha;
0049   m_values[eBeta] = beta;
0050   // now calculate the remaining max half X
0051   double gamma = (alpha > beta) ? (alpha - std::numbers::pi / 2.)
0052                                 : (beta - std::numbers::pi / 2.);
0053   m_values[eHalfLengthXposY] = minhalex + (2. * haley) * std::tan(gamma);
0054
0055   checkConsistency();
0056   buildSurfaceBounds();
0057 }
0058
0059 std::vector<double> TrapezoidVolumeBounds::values() const {
0060   return {m_values.begin(), m_values.end()};
0061 }
0062
0063 std::vector<OrientedSurface> TrapezoidVolumeBounds::orientedSurfaces(
0064     const Transform3& transform) const {
0065   std::vector<OrientedSurface> oSurfaces;
0066   oSurfaces.reserve(6);
0067
0068   // Face surfaces xy
0069   RotationMatrix3 trapezoidRotation(transform.rotation());
0070   Vector3 trapezoidX(trapezoidRotation.col(0));
0071   Vector3 trapezoidY(trapezoidRotation.col(1));
0072   Vector3 trapezoidZ(trapezoidRotation.col(2));
0073   Vector3 trapezoidCenter(transform.translation());
0074
0075   //   (1) - At negative local z
0076   auto nzTransform = transform * Translation3(0., 0., -get(eHalfLengthZ));
0077   auto sf =
0078       Surface::makeShared<PlaneSurface>(nzTransform, m_faceXYTrapezoidBounds);
0079   oSurfaces.emplace_back(std::move(sf), Direction::AlongNormal());
0080   //   (2) - At positive local z
0081   auto pzTransform = transform * Translation3(0., 0., get(eHalfLengthZ));
0082   sf = Surface::makeShared<PlaneSurface>(pzTransform, m_faceXYTrapezoidBounds);
0083   oSurfaces.emplace_back(std::move(sf), Direction::OppositeNormal());
0084
0085   double poshOffset = get(eHalfLengthY) / std::tan(get(eAlpha));
0086   double neghOffset = get(eHalfLengthY) / std::tan(get(eBeta));
0087   double topShift = poshOffset + neghOffset;
0088
0089   // Face surfaces yz
0090   // (3) - At point B, attached to beta opening angle
0091   Vector3 fbPosition(-get(eHalfLengthXnegY) + neghOffset, 0., 0.);
0092   auto fbTransform =
0093       transform * Translation3(fbPosition) *
0094       AngleAxis3(-std::numbers::pi / 2. + get(eBeta), Vector3(0., 0., 1.)) *
0095       s_planeYZ;
0096   sf =
0097       Surface::makeShared<PlaneSurface>(fbTransform, m_faceBetaRectangleBounds);
0098   oSurfaces.emplace_back(std::move(sf), Direction::AlongNormal());
0099
0100   // (4) - At point A, attached to alpha opening angle
0101   Vector3 faPosition(get(eHalfLengthXnegY) + poshOffset, 0., 0.);
0102   auto faTransform =
0103       transform * Translation3(faPosition) *
0104       AngleAxis3(-std::numbers::pi / 2. + get(eAlpha), Vector3(0., 0., 1.)) *
0105       s_planeYZ;
0106   sf = Surface::makeShared<PlaneSurface>(faTransform,
0107                                          m_faceAlphaRectangleBounds);
0108   oSurfaces.emplace_back(std::move(sf), Direction::OppositeNormal());
0109
0110   // Face surfaces zx
0111   //   (5) - At negative local y
0112   auto nxTransform =
0113       transform * Translation3(0., -get(eHalfLengthY), 0.) * s_planeZX;
0114   sf = Surface::makeShared<PlaneSurface>(nxTransform,
0115                                          m_faceZXRectangleBoundsBottom);
0116   oSurfaces.emplace_back(std::move(sf), Direction::AlongNormal());
0117   //   (6) - At positive local y
0118   auto pxTransform =
0119       transform * Translation3(topShift, get(eHalfLengthY), 0.) * s_planeZX;
0120   sf = Surface::makeShared<PlaneSurface>(pxTransform,
0121                                          m_faceZXRectangleBoundsTop);
0122   oSurfaces.emplace_back(std::move(sf), Direction::OppositeNormal());
0123
0124   return oSurfaces;
0125 }
0126
0127 void TrapezoidVolumeBounds::buildSurfaceBounds() {
0128   m_faceXYTrapezoidBounds = std::make_shared<const TrapezoidBounds>(
0129       get(eHalfLengthXnegY), get(eHalfLengthXposY), get(eHalfLengthY));
0130
0131   m_faceAlphaRectangleBounds = std::make_shared<const RectangleBounds>(
0132       get(eHalfLengthY) / std::cos(get(eAlpha) - std::numbers::pi / 2.),
0133       get(eHalfLengthZ));
0134
0135   m_faceBetaRectangleBounds = std::make_shared<const RectangleBounds>(
0136       get(eHalfLengthY) / std::cos(get(eBeta) - std::numbers::pi / 2.),
0137       get(eHalfLengthZ));
0138
0139   m_faceZXRectangleBoundsBottom = std::make_shared<const RectangleBounds>(
0140       get(eHalfLengthZ), get(eHalfLengthXnegY));
0141
0142   m_faceZXRectangleBoundsTop = std::make_shared<const RectangleBounds>(
0143       get(eHalfLengthZ), get(eHalfLengthXposY));
0144 }
0145
0146 bool TrapezoidVolumeBounds::inside(const Vector3& pos, double tol) const {
0147   if (std::abs(pos.z()) > get(eHalfLengthZ) + tol) {
0148     return false;
0149   }
0150   if (std::abs(pos.y()) > get(eHalfLengthY) + tol) {
0151     return false;
0152   }
0153   Vector2 locp(pos.x(), pos.y());
0154   return m_faceXYTrapezoidBounds->inside(
0155       locp, BoundaryTolerance::AbsoluteEuclidean(tol));
0156 }
0157
0158 std::ostream& TrapezoidVolumeBounds::toStream(std::ostream& os) const {
0159   detail::OstreamStateGuard guard{os};
0160   os << std::fixed << std::setprecision(5);
0161   os << "TrapezoidVolumeBounds: (halfX @-Y, halfX @+Y, halfY, halfZ, alpha, "
0162         "beta) "
0163         "= ";
0164   os << "(" << get(eHalfLengthXnegY) << ", " << get(eHalfLengthXposY) << ", "
0165      << get(eHalfLengthY) << ", " << get(eHalfLengthZ);
0166   os << ", " << get(eAlpha) << ", " << get(eBeta) << ")";
0167   return os;
0168 }
0169
0170 Volume::BoundingBox TrapezoidVolumeBounds::boundingBox(
0171     const Transform3* trf, const Vector3& envelope,
0172     const Volume* entity) const {
0173   double minx = get(eHalfLengthXnegY);
0174   double maxx = get(eHalfLengthXposY);
0175   double haley = get(eHalfLengthY);
0176   double halez = get(eHalfLengthZ);
0177
0178   std::array<Vector3, 8> vertices = {{{-minx, -haley, -halez},
0179                                       {+minx, -haley, -halez},
0180                                       {-maxx, +haley, -halez},
0181                                       {+maxx, +haley, -halez},
0182                                       {-minx, -haley, +halez},
0183                                       {+minx, -haley, +halez},
0184                                       {-maxx, +haley, +halez},
0185                                       {+maxx, +haley, +halez}}};
0186
0187   Transform3 transform = Transform3::Identity();
0188   if (trf != nullptr) {
0189     transform = *trf;
0190   }
0191
0192   Vector3 vmin = transform * vertices[0];
0193   Vector3 vmax = transform * vertices[0];
0194
0195   for (std::size_t i = 1; i < 8; i++) {
0196     const Vector3 vtx = transform * vertices[i];
0197     vmin = vmin.cwiseMin(vtx);
0198     vmax = vmax.cwiseMax(vtx);
0199   }
0200
0201   return {entity, vmin - envelope, vmax + envelope};
0202 }
0203
0204 void TrapezoidVolumeBounds::checkConsistency() noexcept(false) {
0205   if (get(eHalfLengthXnegY) < 0. || get(eHalfLengthXposY) < 0.) {
0206     throw std::invalid_argument(
0207         "TrapezoidVolumeBounds: invalid trapezoid parameters in x.");
0208   }
0209   if (get(eHalfLengthY) <= 0.) {
0210     throw std::invalid_argument("TrapezoidVolumeBounds: invalid y extrusion.");
0211   }
0212   if (get(eHalfLengthZ) <= 0.) {
0213     throw std::invalid_argument("TrapezoidVolumeBounds: invalid z extrusion.");
0214   }
0215 }
0216
0217 }  // namespace Acts