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