EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
TrapezoidVolumeBounds.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file TrapezoidVolumeBounds.cpp
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2016-2018 CERN for the benefit of the Acts project
4 //
5 // This Source Code Form is subject to the terms of the Mozilla Public
6 // License, v. 2.0. If a copy of the MPL was not distributed with this
7 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
8 
9 #include "Acts/Geometry/TrapezoidVolumeBounds.hpp"
10 
11 #include "Acts/Geometry/GeometryStatics.hpp"
12 #include "Acts/Surfaces/PlaneSurface.hpp"
13 #include "Acts/Surfaces/RectangleBounds.hpp"
14 #include "Acts/Surfaces/TrapezoidBounds.hpp"
15 #include "Acts/Utilities/BoundingBox.hpp"
16 
17 #include <cmath>
18 #include <iomanip>
19 #include <iostream>
20 
21 Acts::TrapezoidVolumeBounds::TrapezoidVolumeBounds(double minhalex,
22  double maxhalex,
23  double haley, double halez)
24  : VolumeBounds() {
25  m_values[eHalfLengthXnegY] = minhalex;
26  m_values[eHalfLengthXposY] = maxhalex;
27  m_values[eHalfLengthY] = haley;
28  m_values[eHalfLengthZ] = halez;
29  m_values[eAlpha] = atan2(2 * haley, (maxhalex - minhalex));
30  m_values[eBeta] = M_PI - get(eAlpha);
31  checkConsistency();
32  buildSurfaceBounds();
33 }
34 
35 Acts::TrapezoidVolumeBounds::TrapezoidVolumeBounds(double minhalex,
36  double haley, double halez,
37  double alpha, double beta)
38  : VolumeBounds() {
39  m_values[eHalfLengthXnegY] = minhalex;
40  m_values[eHalfLengthY] = haley;
41  m_values[eHalfLengthZ] = halez;
42  m_values[eAlpha] = alpha;
43  m_values[eBeta] = beta;
44  // now calculate the remaining max half X
45  double gamma = (alpha > beta) ? (alpha - 0.5 * M_PI) : (beta - 0.5 * M_PI);
46  m_values[eHalfLengthXposY] = minhalex + (2. * haley) * tan(gamma);
47 
48  checkConsistency();
49  buildSurfaceBounds();
50 }
51 
52 Acts::OrientedSurfaces Acts::TrapezoidVolumeBounds::orientedSurfaces(
53  const Transform3D& transform) const {
54  OrientedSurfaces oSurfaces;
55  oSurfaces.reserve(6);
56 
57  // Face surfaces xy
58  RotationMatrix3D trapezoidRotation(transform.rotation());
59  Vector3D trapezoidX(trapezoidRotation.col(0));
60  Vector3D trapezoidY(trapezoidRotation.col(1));
61  Vector3D trapezoidZ(trapezoidRotation.col(2));
62  Vector3D trapezoidCenter(transform.translation());
63 
64  // (1) - At negative local z
65  auto nzTransform = transform * Translation3D(0., 0., -get(eHalfLengthZ));
66  auto sf =
67  Surface::makeShared<PlaneSurface>(nzTransform, m_faceXYTrapezoidBounds);
68  oSurfaces.push_back(OrientedSurface(std::move(sf), forward));
69  // (2) - At positive local z
70  auto pzTransform = transform * Translation3D(0., 0., get(eHalfLengthZ));
71  sf = Surface::makeShared<PlaneSurface>(pzTransform, m_faceXYTrapezoidBounds);
72  oSurfaces.push_back(OrientedSurface(std::move(sf), backward));
73 
74  double poshOffset = get(eHalfLengthY) / std::tan(get(eAlpha));
75  double neghOffset = get(eHalfLengthY) / std::tan(get(eBeta));
76  double topShift = poshOffset + neghOffset;
77 
78  // Face surfaces yz
79  // (3) - At point B, attached to beta opening angle
80  Vector3D fbPosition(-get(eHalfLengthXnegY) + neghOffset, 0., 0.);
81  auto fbTransform =
82  transform * Translation3D(fbPosition) *
83  AngleAxis3D(-0.5 * M_PI + get(eBeta), Vector3D(0., 0., 1.)) * s_planeYZ;
84  sf =
85  Surface::makeShared<PlaneSurface>(fbTransform, m_faceBetaRectangleBounds);
86  oSurfaces.push_back(OrientedSurface(std::move(sf), forward));
87 
88  // (4) - At point A, attached to alpha opening angle
89  Vector3D faPosition(get(eHalfLengthXnegY) + poshOffset, 0., 0.);
90  auto faTransform =
91  transform * Translation3D(faPosition) *
92  AngleAxis3D(-0.5 * M_PI + get(eAlpha), Vector3D(0., 0., 1.)) * s_planeYZ;
93  sf = Surface::makeShared<PlaneSurface>(faTransform,
94  m_faceAlphaRectangleBounds);
95  oSurfaces.push_back(OrientedSurface(std::move(sf), backward));
96 
97  // Face surfaces zx
98  // (5) - At negative local y
99  auto nxTransform =
100  transform * Translation3D(0., -get(eHalfLengthY), 0.) * s_planeZX;
101  sf = Surface::makeShared<PlaneSurface>(nxTransform,
102  m_faceZXRectangleBoundsBottom);
103  oSurfaces.push_back(OrientedSurface(std::move(sf), forward));
104  // (6) - At positive local y
105  auto pxTransform =
106  transform * Translation3D(topShift, get(eHalfLengthY), 0.) * s_planeZX;
107  sf = Surface::makeShared<PlaneSurface>(pxTransform,
108  m_faceZXRectangleBoundsTop);
109  oSurfaces.push_back(OrientedSurface(std::move(sf), backward));
110 
111  return oSurfaces;
112 }
113 
114 void Acts::TrapezoidVolumeBounds::buildSurfaceBounds() {
115  m_faceXYTrapezoidBounds = std::make_shared<const TrapezoidBounds>(
116  get(eHalfLengthXnegY), get(eHalfLengthXposY), get(eHalfLengthY));
117 
118  m_faceAlphaRectangleBounds = std::make_shared<const RectangleBounds>(
119  get(eHalfLengthY) / cos(get(eAlpha) - 0.5 * M_PI), get(eHalfLengthZ));
120 
121  m_faceBetaRectangleBounds = std::make_shared<const RectangleBounds>(
122  get(eHalfLengthY) / cos(get(eBeta) - 0.5 * M_PI), get(eHalfLengthZ));
123 
124  m_faceZXRectangleBoundsBottom = std::make_shared<const RectangleBounds>(
125  get(eHalfLengthZ), get(eHalfLengthXnegY));
126 
127  m_faceZXRectangleBoundsTop = std::make_shared<const RectangleBounds>(
128  get(eHalfLengthZ), get(eHalfLengthXposY));
129 }
130 
131 bool Acts::TrapezoidVolumeBounds::inside(const Vector3D& pos,
132  double tol) const {
133  if (std::abs(pos.z()) > get(eHalfLengthZ) + tol) {
134  return false;
135  }
136  if (std::abs(pos.y()) > get(eHalfLengthY) + tol) {
137  return false;
138  }
139  Vector2D locp(pos.x(), pos.y());
140  bool inside(m_faceXYTrapezoidBounds->inside(
141  locp, BoundaryCheck(true, true, tol, tol)));
142  return inside;
143 }
144 
145 std::ostream& Acts::TrapezoidVolumeBounds::toStream(std::ostream& sl) const {
146  return dumpT<std::ostream>(sl);
147 }
148 
149 Acts::Volume::BoundingBox Acts::TrapezoidVolumeBounds::boundingBox(
150  const Acts::Transform3D* trf, const Vector3D& envelope,
151  const Volume* entity) const {
152  double minx = get(eHalfLengthXnegY);
153  double maxx = get(eHalfLengthXposY);
154  double haley = get(eHalfLengthY);
155  double halez = get(eHalfLengthZ);
156 
157  std::array<Vector3D, 8> vertices = {{{-minx, -haley, -halez},
158  {+minx, -haley, -halez},
159  {-maxx, +haley, -halez},
160  {+maxx, +haley, -halez},
161  {-minx, -haley, +halez},
162  {+minx, -haley, +halez},
163  {-maxx, +haley, +halez},
164  {+maxx, +haley, +halez}}};
165 
166  Transform3D transform = Transform3D::Identity();
167  if (trf != nullptr) {
168  transform = *trf;
169  }
170 
171  Vector3D vmin = transform * vertices[0];
172  Vector3D vmax = transform * vertices[0];
173 
174  for (size_t i = 1; i < 8; i++) {
175  const Vector3D vtx = transform * vertices[i];
176  vmin = vmin.cwiseMin(vtx);
177  vmax = vmax.cwiseMax(vtx);
178  }
179 
180  return {entity, vmin - envelope, vmax + envelope};
181 }