src/Geometry/DiamondVolumeBounds.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/DiamondVolumeBounds.hpp"
0010
0011 #include "Acts/Definitions/Direction.hpp"
0012 #include "Acts/Surfaces/BoundaryTolerance.hpp"
0013 #include "Acts/Surfaces/DiamondBounds.hpp"
0014 #include "Acts/Surfaces/PlaneSurface.hpp"
0015 #include "Acts/Surfaces/RectangleBounds.hpp"
0016 #include "Acts/Surfaces/Surface.hpp"
0017
0018 namespace Acts {
0019
0020 DiamondVolumeBounds::DiamondVolumeBounds(double x1, double x2, double x3,
0021                                          double y1, double y2, double halez)
0022
0023     : VolumeBounds() {
0024   m_values[eHalfLengthX1] = x1;
0025   m_values[eHalfLengthX2] = x2;
0026   m_values[eHalfLengthX3] = x3;
0027   m_values[eLengthY1] = y1;
0028   m_values[eLengthY2] = y2;
0029   m_values[eHalfLengthZ] = halez;
0030   m_values[eAlphaAngle] = std::numbers::pi - std::atan2(y1, (x2 - x1));
0031   m_values[eBetaAngle] = std::atan2(y2, (x2 - x3));
0032
0033   checkConsistency();
0034   buildSurfaceBounds();
0035 }
0036
0037 std::vector<double> DiamondVolumeBounds::values() const {
0038   return {m_values.begin(), m_values.end()};
0039 }
0040
0041 std::vector<OrientedSurface> DiamondVolumeBounds::orientedSurfaces(
0042     const Transform3& transform) const {
0043   std::vector<OrientedSurface> surfaces;
0044   surfaces.reserve(8);
0045
0046   // (0) - At negative local z
0047   auto negZTransform = transform * Translation3(0., 0., -get(eHalfLengthZ));
0048   auto sf = Surface::makeShared<PlaneSurface>(negZTransform, m_FaceXYBounds);
0049   surfaces.push_back(OrientedSurface{std::move(sf), Direction::AlongNormal()});
0050
0051   // (1) - At positive local z
0052   auto posZTransform = transform * Translation3(0., 0., get(eHalfLengthZ));
0053   sf = Surface::makeShared<PlaneSurface>(posZTransform, m_FaceXYBounds);
0054   surfaces.push_back(
0055       OrientedSurface{std::move(sf), Direction::OppositeNormal()});
0056
0057   double posXOffset23 = 0.5 * get(eLengthY2) / std::tan(get(eBetaAngle));
0058   double posXOffset12 =
0059       0.5 * get(eLengthY1) / std::tan(std::numbers::pi - get(eAlphaAngle));
0060
0061   // (2) - At negative x face yz12
0062   Vector3 nyz12Position(-get(eHalfLengthX1) - posXOffset12,
0063                         -0.5 * get(eLengthY1), 0.);
0064   auto nyz12Transform =
0065       transform * Translation3(nyz12Position) *
0066       AngleAxis3(-std::numbers::pi / 2. + get(eAlphaAngle), Vector3::UnitZ()) *
0067       s_planeYZ;
0068   sf = Surface::makeShared<PlaneSurface>(nyz12Transform, m_FaceYZ12Bounds);
0069   surfaces.push_back(OrientedSurface{std::move(sf), Direction::AlongNormal()});
0070
0071   // (3) - At positive x face yz12
0072   Vector3 pyz12Position(get(eHalfLengthX1) + posXOffset12,
0073                         -0.5 * get(eLengthY1), 0.);
0074   auto pyz12Transform =
0075       transform * Translation3(pyz12Position) *
0076       AngleAxis3(-std::numbers::pi / 2. + get(eAlphaAngle), -Vector3::UnitZ()) *
0077       s_planeYZ;
0078   sf = Surface::makeShared<PlaneSurface>(pyz12Transform, m_FaceYZ12Bounds);
0079   surfaces.push_back(
0080       OrientedSurface{std::move(sf), Direction::OppositeNormal()});
0081
0082   // (4) - At negative x face yz23
0083   Vector3 nyz23Position(-get(eHalfLengthX3) - posXOffset23,
0084                         0.5 * get(eLengthY2), 0.);
0085   auto nyz23Transform = transform * Translation3(nyz23Position) *
0086                         AngleAxis3(std::numbers::pi / 2. - get(eBetaAngle),
0087                                    Vector3(0., 0., -1.)) *
0088                         s_planeYZ;
0089   sf = Surface::makeShared<PlaneSurface>(nyz23Transform, m_FaceYZ23Bounds);
0090   surfaces.push_back(OrientedSurface{std::move(sf), Direction::AlongNormal()});
0091
0092   // (5) - At positive x face yz23
0093   Vector3 pyz23Position(get(eHalfLengthX3) + posXOffset23, 0.5 * get(eLengthY2),
0094                         0.);
0095   auto pyz23Transform =
0096       transform * Translation3(pyz23Position) *
0097       AngleAxis3(std::numbers::pi / 2. - get(eBetaAngle), Vector3::UnitZ()) *
0098       s_planeYZ;
0099
0100   sf = Surface::makeShared<PlaneSurface>(pyz23Transform, m_FaceYZ23Bounds);
0101   surfaces.push_back(
0102       OrientedSurface{std::move(sf), Direction::OppositeNormal()});
0103
0104   // (6) - At negative y face zx
0105   auto nyTransform =
0106       transform * Translation3(0., -get(eLengthY1), 0.) * s_planeZX;
0107   sf = Surface::makeShared<PlaneSurface>(nyTransform, m_negYFaceZXBounds);
0108   surfaces.push_back(OrientedSurface{std::move(sf), Direction::AlongNormal()});
0109
0110   // (7) - At positive y face zx
0111   auto pyTransform =
0112       transform * Translation3(0., get(eLengthY2), 0.) * s_planeZX;
0113   sf = Surface::makeShared<PlaneSurface>(pyTransform, m_posYFaceZXBounds);
0114   surfaces.push_back(
0115       OrientedSurface{std::move(sf), Direction::OppositeNormal()});
0116
0117   return surfaces;
0118 };
0119
0120 void DiamondVolumeBounds::buildSurfaceBounds() {
0121   m_FaceXYBounds = std::make_shared<Acts::DiamondBounds>(
0122       get(eHalfLengthX1), get(eHalfLengthX2), get(eHalfLengthX3),
0123       get(eLengthY1), get(eLengthY2));
0124
0125   double dx23 =
0126       fastHypot((get(eHalfLengthX2) - get(eHalfLengthX3)), get(eLengthY2));
0127   double dx12 =
0128       fastHypot((get(eHalfLengthX2) - get(eHalfLengthX1)), get(eLengthY1));
0129
0130   m_FaceYZ12Bounds =
0131       std::make_shared<Acts::RectangleBounds>(0.5 * dx12, get(eHalfLengthZ));
0132   m_FaceYZ23Bounds =
0133       std::make_shared<Acts::RectangleBounds>(0.5 * dx23, get(eHalfLengthZ));
0134
0135   m_posYFaceZXBounds = std::make_shared<Acts::RectangleBounds>(
0136       get(eHalfLengthZ), get(eHalfLengthX3));
0137   m_negYFaceZXBounds = std::make_shared<Acts::RectangleBounds>(
0138       get(eHalfLengthZ), get(eHalfLengthX1));
0139 }
0140
0141 Volume::BoundingBox DiamondVolumeBounds::boundingBox(
0142     const Transform3* trf, const Vector3& envelope,
0143     const Volume* entity) const {
0144   double maxX =
0145       std::max({get(eHalfLengthX1), get(eHalfLengthX2), get(eHalfLengthX3)});
0146   double maxY = std::max(get(eLengthY1), get(eLengthY2));
0147   double halez = get(eHalfLengthZ);
0148
0149   std::array<Vector3, 8> vertices = {{{-maxX, -maxY, -halez},
0150                                       {+maxX, -maxY, -halez},
0151                                       {-maxX, +maxY, -halez},
0152                                       {+maxX, +maxY, -halez},
0153                                       {-maxX, -maxY, +halez},
0154                                       {+maxX, -maxY, +halez},
0155                                       {-maxX, +maxY, +halez},
0156                                       {+maxX, +maxY, +halez}}};
0157
0158   Transform3 transform = Transform3::Identity();
0159   if (trf != nullptr) {
0160     transform = *trf;
0161   }
0162
0163   Vector3 vmin = transform * vertices[0];
0164   Vector3 vmax = transform * vertices[0];
0165
0166   for (std::size_t i = 1; i < 8; i++) {
0167     const Vector3 vtx = transform * vertices[i];
0168     vmin = vmin.cwiseMin(vtx);
0169     vmax = vmax.cwiseMax(vtx);
0170   }
0171
0172   return {entity, vmin - envelope, vmax + envelope};
0173 }
0174 bool DiamondVolumeBounds::inside(const Vector3& pos, double tol) const {
0175   if (std::abs(pos.z()) > get(eHalfLengthZ) + tol) {
0176     return false;
0177   }
0178
0179   if (std::abs(pos.y()) > std::max(get(eLengthY1), get(eLengthY2)) + tol) {
0180     return false;
0181   }
0182
0183   Vector2 locPos(pos.x(), pos.y());
0184   return m_FaceXYBounds->inside(locPos,
0185                                 BoundaryTolerance::AbsoluteEuclidean(tol));
0186 }
0187
0188 void DiamondVolumeBounds::checkConsistency() noexcept(false) {
0189   if (get(eHalfLengthX1) <= 0. || get(eHalfLengthX2) <= 0. ||
0190       get(eHalfLengthX3) <= 0.) {
0191     throw std::invalid_argument(
0192         "DiamondVolumeBounds: invalid polygon parameters in x.");
0193   }
0194   if (get(eLengthY1) <= 0. || get(eLengthY2) <= 0.) {
0195     throw std::invalid_argument("DiamondVolumeBounds: invalid y extrusion.");
0196   }
0197   if (get(eHalfLengthZ) <= 0.) {
0198     throw std::invalid_argument("DiamondVolumeBounds: invalid y extrusion.");
0199   }
0200
0201   // make sure this is a convex polygon - do not allow angles > 180 deg
0202   if (get(eHalfLengthX2) < get(eHalfLengthX1) ||
0203       get(eHalfLengthX2) < get(eHalfLengthX3)) {
0204     throw std::invalid_argument(
0205         "DiamondVolumeBounds: invalid polygon shape - not convex.");
0206   }
0207 }
0208
0209 std::ostream& DiamondVolumeBounds::toStream(std::ostream& os) const {
0210   os << std::setiosflags(std::ios::fixed);
0211   os << std::setprecision(5);
0212   os << "DiamondVolumeBounds: (halfX1, halfX2, halfX3, halfY1, halfY2, "
0213         "halfZ) = ";
0214   os << "(" << get(eHalfLengthX1) << ", " << get(eHalfLengthX2) << ", "
0215      << get(eHalfLengthX3) << ", " << get(eLengthY1) << ", " << get(eLengthY2)
0216      << ", " << get(eHalfLengthZ) << ")";
0217   return os;
0218 }
0219
0220 }  // namespace Acts