EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
GenericCuboidVolumeBounds.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file GenericCuboidVolumeBounds.cpp
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2019 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 
11 #include "Acts/Geometry/Volume.hpp"
18 
19 #include <array>
20 #include <ostream>
21 
23  const std::array<Acts::Vector3D, 8>& vertices) noexcept(false)
24  : m_vertices(vertices) {
25  construct();
26 }
27 
29  const std::array<double, GenericCuboidVolumeBounds::eSize>&
30  values) noexcept(false)
31  : m_vertices() {
32  for (size_t iv = 0; iv < 8; ++iv) {
33  m_vertices[iv] =
34  Vector3D(values[iv * 3], values[iv * 3 + 1], values[iv * 3 + 2]);
35  }
36  construct();
37 }
38 
40  double tol) const {
41  constexpr std::array<size_t, 6> vtxs = {0, 4, 0, 1, 2, 1};
42  // needs to be on same side, get ref
43  bool ref = std::signbit((gpos - m_vertices[vtxs[0]]).dot(m_normals[0]));
44  for (size_t i = 1; i < 6; i++) {
45  double dot = (gpos - m_vertices[vtxs[i]]).dot(m_normals[i]);
46  if (std::signbit(dot) != ref) {
47  // technically outside, but how far?
48  if (std::abs(dot) > tol) {
49  // distance greater than tol
50  return false;
51  }
52  // distance smaller than tol, ignore
53  }
54  }
55  return true;
56 }
57 
59  const Transform3D& transform) const {
60  OrientedSurfaces oSurfaces;
61 
62  // approximate cog of the volume
63  Vector3D cog(0, 0, 0);
64 
65  for (size_t i = 0; i < 8; i++) {
66  cog += m_vertices[i];
67  }
68 
69  cog *= 0.125; // 1/8.
70 
71  auto make_surface = [&](const auto& a, const auto& b, const auto& c,
72  const auto& d) -> void {
73  // calculate centroid of these points
74  Vector3D ctrd = (a + b + c + d) / 4.;
75  // create normal
76  const Vector3D ab = b - a, ac = c - a;
77  Vector3D normal = ab.cross(ac).normalized();
78 
79  NavigationDirection nDir = ((cog - d).dot(normal) < 0) ? backward : forward;
80 
81  // build transform from z unit to normal
82  // z is normal in local coordinates
83  // Volume local to surface local
84  Transform3D vol2srf;
85  vol2srf = (Eigen::Quaternion<double>().setFromTwoVectors(
86  normal, Vector3D::UnitZ()));
87 
88  vol2srf = vol2srf * Translation3D(-ctrd);
89 
90  // now calculate position of vertices in surface local frame
91  Vector3D a_l, b_l, c_l, d_l;
92  a_l = vol2srf * a;
93  b_l = vol2srf * b;
94  c_l = vol2srf * c;
95  d_l = vol2srf * d;
96 
97  std::vector<Vector2D> vertices({{a_l.x(), a_l.y()},
98  {b_l.x(), b_l.y()},
99  {c_l.x(), c_l.y()},
100  {d_l.x(), d_l.y()}});
101 
102  auto polyBounds = std::make_shared<const ConvexPolygonBounds<4>>(vertices);
103  auto srfTrf = transform * vol2srf.inverse();
104  auto srf = Surface::makeShared<PlaneSurface>(srfTrf, polyBounds);
105 
106  oSurfaces.push_back(OrientedSurface(std::move(srf), nDir));
107  };
108 
109  make_surface(m_vertices[0], m_vertices[1], m_vertices[2], m_vertices[3]);
110  make_surface(m_vertices[4], m_vertices[5], m_vertices[6], m_vertices[7]);
111  make_surface(m_vertices[0], m_vertices[3], m_vertices[7], m_vertices[4]);
112  make_surface(m_vertices[1], m_vertices[2], m_vertices[6], m_vertices[5]);
113  make_surface(m_vertices[2], m_vertices[3], m_vertices[7], m_vertices[6]);
114  make_surface(m_vertices[1], m_vertices[0], m_vertices[4], m_vertices[5]);
115 
116  return oSurfaces;
117 }
118 
120  std::ostream& sl) const {
121  sl << "Acts::GenericCuboidVolumeBounds: vertices (x, y, z) =\n";
122  for (size_t i = 0; i < 8; i++) {
123  if (i > 0) {
124  sl << ",\n";
125  }
126  sl << "[" << m_vertices[i].transpose() << "]";
127  }
128  return sl;
129 }
130 
132  // calculate approximate center of gravity first, so we can make sure
133  // the normals point inwards
134  Vector3D cog(0, 0, 0);
135 
136  for (size_t i = 0; i < 8; i++) {
137  cog += m_vertices[i];
138  }
139 
140  cog *= 0.125; // 1/8.
141 
142  size_t idx = 0;
143 
144  auto handle_face = [&](const auto& a, const auto& b, const auto& c,
145  const auto& d) {
146  // we assume a b c d to be counter clockwise
147  const Vector3D ab = b - a, ac = c - a;
148  Vector3D normal = ab.cross(ac).normalized();
149 
150  if ((cog - a).dot(normal) < 0) {
151  // normal points outwards, flip normal
152  normal *= -1.;
153  }
154 
155  // get rid of -0 values if present
156  normal += Vector3D::Zero();
157 
158  // check if d is on the surface
159  if (std::abs((a - d).dot(normal)) > 1e-6) {
160  throw(std::invalid_argument(
161  "GenericCuboidBounds: Four points do not lie on the same plane!"));
162  }
163 
164  m_normals[idx] = normal;
165  idx++;
166  };
167 
168  // handle faces
169  handle_face(m_vertices[0], m_vertices[1], m_vertices[2], m_vertices[3]);
170  handle_face(m_vertices[4], m_vertices[5], m_vertices[6], m_vertices[7]);
171  handle_face(m_vertices[0], m_vertices[3], m_vertices[7], m_vertices[4]);
172  handle_face(m_vertices[1], m_vertices[2], m_vertices[6], m_vertices[5]);
173  handle_face(m_vertices[2], m_vertices[3], m_vertices[7], m_vertices[6]);
174  handle_face(m_vertices[1], m_vertices[0], m_vertices[4], m_vertices[5]);
175 }
176 
177 std::vector<double> Acts::GenericCuboidVolumeBounds::values() const {
178  std::vector<double> rvalues;
179  rvalues.reserve(eSize);
180  for (size_t iv = 0; iv < 8; ++iv) {
181  for (size_t ic = 0; ic < 3; ++ic) {
182  rvalues.push_back(m_vertices[iv][ic]);
183  }
184  }
185  return rvalues;
186 }
187 
189  const Acts::Transform3D* trf, const Vector3D& envelope,
190  const Volume* entity) const {
191  Vector3D vmin, vmax;
192 
193  Transform3D transform = Transform3D::Identity();
194  if (trf != nullptr) {
195  transform = *trf;
196  }
197 
198  vmin = transform * m_vertices[0];
199  vmax = transform * m_vertices[0];
200 
201  for (size_t i = 1; i < 8; i++) {
202  Vector3D vtx = transform * m_vertices[i];
203  vmin = vmin.cwiseMin(vtx);
204  vmax = vmax.cwiseMax(vtx);
205  }
206 
207  return {entity, vmin - envelope, vmax + envelope};
208 }
209 
211  const Transform3D& transform) const {
212  auto draw_face = [&](const auto& a, const auto& b, const auto& c,
213  const auto& d) {
214  helper.face(std::vector<Vector3D>(
215  {transform * a, transform * b, transform * c, transform * d}));
216  };
217 
218  draw_face(m_vertices[0], m_vertices[1], m_vertices[2], m_vertices[3]);
219  draw_face(m_vertices[4], m_vertices[5], m_vertices[6], m_vertices[7]);
220  draw_face(m_vertices[0], m_vertices[3], m_vertices[7], m_vertices[4]);
221  draw_face(m_vertices[1], m_vertices[2], m_vertices[6], m_vertices[5]);
222  draw_face(m_vertices[2], m_vertices[3], m_vertices[7], m_vertices[6]);
223  draw_face(m_vertices[1], m_vertices[0], m_vertices[4], m_vertices[5]);
224 }