ECCE @ 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 
10 // TrapezoidVolumeBounds.cpp, Acts project
12 
13 #include "Acts/Geometry/TrapezoidVolumeBounds.hpp"
14 
15 #include "Acts/Geometry/GeometryStatics.hpp"
16 #include "Acts/Surfaces/PlaneSurface.hpp"
17 #include "Acts/Surfaces/RectangleBounds.hpp"
18 #include "Acts/Surfaces/TrapezoidBounds.hpp"
19 #include "Acts/Utilities/BoundingBox.hpp"
20 
21 #include <cmath>
22 #include <iomanip>
23 #include <iostream>
24 
25 Acts::TrapezoidVolumeBounds::TrapezoidVolumeBounds(double minhalex,
26  double maxhalex,
27  double haley, double halez)
28  : VolumeBounds() {
29  m_values[eHalfLengthXnegY] = minhalex;
30  m_values[eHalfLengthXposY] = maxhalex;
31  m_values[eHalfLengthY] = haley;
32  m_values[eHalfLengthZ] = halez;
33  m_values[eAlpha] =
34  atan((m_values[eHalfLengthXposY] - m_values[eHalfLengthXnegY]) / 2 /
35  m_values[eHalfLengthY]) +
36  0.5 * M_PI;
37  m_values[eBeta] = m_values[eAlpha];
38  checkConsistency();
39  buildSurfaceBounds();
40 }
41 
42 Acts::TrapezoidVolumeBounds::TrapezoidVolumeBounds(double minhalex,
43  double haley, double halez,
44  double alpha, double beta)
45  : VolumeBounds() {
46  m_values[eHalfLengthXnegY] = minhalex;
47  m_values[eHalfLengthY] = haley;
48  m_values[eHalfLengthZ] = halez;
49  m_values[eAlpha] = alpha;
50  m_values[eBeta] = beta;
51  // now calculate the remaining max half X
52  double gamma = (alpha > beta) ? (alpha - 0.5 * M_PI) : (beta - 0.5 * M_PI);
53  m_values[eHalfLengthXposY] = minhalex + (2. * haley) * tan(gamma);
54 
55  checkConsistency();
56  buildSurfaceBounds();
57 }
58 
59 std::vector<std::shared_ptr<const Acts::Surface>>
60 Acts::TrapezoidVolumeBounds::decomposeToSurfaces(
61  const Transform3D* transformPtr) const {
62  std::vector<std::shared_ptr<const Surface>> rSurfaces;
63 
64  Transform3D transform =
65  (transformPtr == nullptr) ? Transform3D::Identity() : (*transformPtr);
66  const Transform3D* tTransform = nullptr;
67  // face surfaces xy
68  RotationMatrix3D trapezoidRotation(transform.rotation());
69  Vector3D trapezoidX(trapezoidRotation.col(0));
70  Vector3D trapezoidY(trapezoidRotation.col(1));
71  Vector3D trapezoidZ(trapezoidRotation.col(2));
72  Vector3D trapezoidCenter(transform.translation());
73 
74  // (1) - at negative local z
75  tTransform =
76  new Transform3D(transform * AngleAxis3D(M_PI, Vector3D(0., 1., 0.)) *
77  Translation3D(Vector3D(0., 0., get(eHalfLengthZ))));
78 
79  rSurfaces.push_back(Surface::makeShared<PlaneSurface>(
80  std::shared_ptr<const Transform3D>(tTransform), m_faceXYTrapezoidBounds));
81  // (2) - at positive local z
82  tTransform = new Transform3D(
83  transform * Translation3D(Vector3D(0., 0., get(eHalfLengthZ))));
84  rSurfaces.push_back(Surface::makeShared<PlaneSurface>(
85  std::shared_ptr<const Transform3D>(tTransform), m_faceXYTrapezoidBounds));
86 
87  // face surfaces yz
88  // transmute cyclical
89  // (3) - at point A, attached to alpha opening angle
90  Vector3D A(get(eHalfLengthXnegY), get(eHalfLengthY), trapezoidCenter.z());
91  RotationMatrix3D alphaZRotation =
92  (s_idRotation *
93  AngleAxis3D(get(eAlpha) - 0.5 * M_PI, Vector3D(0., 0., 1.)))
94  .toRotationMatrix();
95  RotationMatrix3D faceAlphaRotation;
96  faceAlphaRotation.col(0) = alphaZRotation.col(1);
97  faceAlphaRotation.col(1) = -alphaZRotation.col(2);
98  faceAlphaRotation.col(2) = -alphaZRotation.col(0);
99 
100  // Vector3D
101  // faceAlphaPosition(A+faceAlphaRotation.colX()*m_faceAlphaRectangleBounds->halflengthX());
102  Vector3D faceAlphaPosition0(
103  -0.5 * (get(eHalfLengthXnegY) + get(eHalfLengthXposY)), 0., 0.);
104  Vector3D faceAlphaPosition = transform * faceAlphaPosition0;
105  tTransform = new Transform3D(Translation3D(faceAlphaPosition) *
106  (trapezoidRotation * faceAlphaRotation));
107  rSurfaces.push_back(Surface::makeShared<PlaneSurface>(
108  std::shared_ptr<const Transform3D>(tTransform),
109  m_faceAlphaRectangleBounds));
110 
111  // (4) - at point B, attached to beta opening angle
112  Vector3D B(get(eHalfLengthXnegY), -get(eHalfLengthY), trapezoidCenter.z());
113  RotationMatrix3D betaZRotation =
114  (s_idRotation *
115  AngleAxis3D(-(get(eBeta) - 0.5 * M_PI), Vector3D(0., 0., 1.)))
116  .toRotationMatrix();
117  RotationMatrix3D faceBetaRotation;
118  faceBetaRotation.col(0) = betaZRotation.col(1);
119  faceBetaRotation.col(1) = betaZRotation.col(2);
120  faceBetaRotation.col(2) = betaZRotation.col(0);
121  // Vector3D
122  // faceBetaPosition(B+faceBetaRotation.colX()*m_faceBetaRectangleBounds->halflengthX());
123  Vector3D faceBetaPosition0(
124  0.5 * (get(eHalfLengthXnegY) + get(eHalfLengthXposY)), 0., 0.);
125  Vector3D faceBetaPosition = transform * faceBetaPosition0;
126  tTransform = new Transform3D(Translation3D(faceBetaPosition) *
127  (trapezoidRotation * faceBetaRotation));
128  rSurfaces.push_back(Surface::makeShared<PlaneSurface>(
129  std::shared_ptr<const Transform3D>(tTransform),
130  m_faceBetaRectangleBounds));
131 
132  // face surfaces zx
133  // (5) - at negative local y
134  tTransform =
135  new Transform3D(transform * AngleAxis3D(M_PI, Vector3D(1., 0., 0.)) *
136  Translation3D(Vector3D(0., get(eHalfLengthY), 0.)) *
137  AngleAxis3D(-0.5 * M_PI, Vector3D(0., 1., 0.)) *
138  AngleAxis3D(-0.5 * M_PI, Vector3D(1., 0., 0.)));
139  rSurfaces.push_back(Surface::makeShared<PlaneSurface>(
140  std::shared_ptr<const Transform3D>(tTransform),
141  std::shared_ptr<const PlanarBounds>(m_faceZXRectangleBoundsBottom)));
142  // (6) - at positive local y
143  tTransform = new Transform3D(
144  transform * Translation3D(Vector3D(0., get(eHalfLengthY), 0.)) *
145  AngleAxis3D(-0.5 * M_PI, Vector3D(0., 1., 0.)) *
146  AngleAxis3D(-0.5 * M_PI, Vector3D(1., 0., 0.)));
147  rSurfaces.push_back(Surface::makeShared<PlaneSurface>(
148  std::shared_ptr<const Transform3D>(tTransform),
149  std::shared_ptr<const PlanarBounds>(m_faceZXRectangleBoundsTop)));
150 
151  return rSurfaces;
152 }
153 
154 void Acts::TrapezoidVolumeBounds::buildSurfaceBounds() {
155  m_faceXYTrapezoidBounds = std::make_shared<const TrapezoidBounds>(
156  get(eHalfLengthXnegY), get(eHalfLengthXposY), get(eHalfLengthY));
157 
158  m_faceAlphaRectangleBounds = std::make_shared<const RectangleBounds>(
159  get(eHalfLengthY) / cos(get(eAlpha) - 0.5 * M_PI), get(eHalfLengthZ));
160 
161  m_faceBetaRectangleBounds = std::make_shared<const RectangleBounds>(
162  get(eHalfLengthY) / cos(get(eBeta) - 0.5 * M_PI), get(eHalfLengthZ));
163 
164  m_faceZXRectangleBoundsBottom = std::make_shared<const RectangleBounds>(
165  get(eHalfLengthZ), get(eHalfLengthXnegY));
166 
167  m_faceZXRectangleBoundsTop = std::make_shared<const RectangleBounds>(
168  get(eHalfLengthZ), get(eHalfLengthXposY));
169 }
170 
171 bool Acts::TrapezoidVolumeBounds::inside(const Vector3D& pos,
172  double tol) const {
173  if (std::abs(pos.z()) > get(eHalfLengthZ) + tol) {
174  return false;
175  }
176  if (std::abs(pos.y()) > get(eHalfLengthY) + tol) {
177  return false;
178  }
179  Vector2D locp(pos.x(), pos.y());
180  bool inside(m_faceXYTrapezoidBounds->inside(
181  locp, BoundaryCheck(true, true, tol, tol)));
182  return inside;
183 }
184 
185 std::ostream& Acts::TrapezoidVolumeBounds::toStream(std::ostream& sl) const {
186  return dumpT<std::ostream>(sl);
187 }
188 
189 Acts::Volume::BoundingBox Acts::TrapezoidVolumeBounds::boundingBox(
190  const Acts::Transform3D* trf, const Vector3D& envelope,
191  const Volume* entity) const {
192  double minx = get(eHalfLengthXnegY);
193  double maxx = get(eHalfLengthXposY);
194  double haley = get(eHalfLengthY);
195  double halez = get(eHalfLengthZ);
196 
197  std::array<Vector3D, 8> vertices = {{{-minx, -haley, -halez},
198  {+minx, -haley, -halez},
199  {-maxx, +haley, -halez},
200  {+maxx, +haley, -halez},
201  {-minx, -haley, +halez},
202  {+minx, -haley, +halez},
203  {-maxx, +haley, +halez},
204  {+maxx, +haley, +halez}}};
205 
206  Transform3D transform = Transform3D::Identity();
207  if (trf != nullptr) {
208  transform = *trf;
209  }
210 
211  Vector3D vmin = transform * vertices[0];
212  Vector3D vmax = transform * vertices[0];
213 
214  for (size_t i = 1; i < 8; i++) {
215  const Vector3D vtx = transform * vertices[i];
216  vmin = vmin.cwiseMin(vtx);
217  vmax = vmax.cwiseMax(vtx);
218  }
219 
220  return {entity, vmin - envelope, vmax + envelope};
221 }