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