Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
TGeoSurfaceConverter.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file TGeoSurfaceConverter.cpp
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2020-2023 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/Plugins/TGeo/TGeoSurfaceConverter.hpp"
10 
11 #include "Acts/Definitions/Tolerance.hpp"
12 #include "Acts/Plugins/TGeo/TGeoPrimitivesHelper.hpp"
13 #include "Acts/Surfaces/AnnulusBounds.hpp"
14 #include "Acts/Surfaces/ConvexPolygonBounds.hpp"
15 #include "Acts/Surfaces/CylinderBounds.hpp"
16 #include "Acts/Surfaces/CylinderSurface.hpp"
17 #include "Acts/Surfaces/DiscSurface.hpp"
18 #include "Acts/Surfaces/PlaneSurface.hpp"
19 #include "Acts/Surfaces/RadialBounds.hpp"
20 #include "Acts/Surfaces/RectangleBounds.hpp"
21 #include "Acts/Surfaces/Surface.hpp"
22 #include "Acts/Surfaces/TrapezoidBounds.hpp"
23 #include "Acts/Utilities/Helpers.hpp"
24 
25 #include <algorithm>
26 #include <array>
27 #include <cctype>
28 #include <cstddef>
29 #include <memory>
30 #include <stdexcept>
31 #include <tuple>
32 #include <utility>
33 #include <vector>
34 
35 #include <boost/algorithm/string.hpp>
36 #include <boost/algorithm/string/predicate.hpp>
37 
38 #include "RtypesCore.h"
39 #include "TGeoArb8.h"
40 #include "TGeoBBox.h"
41 #include "TGeoBoolNode.h"
42 #include "TGeoCompositeShape.h"
43 #include "TGeoMatrix.h"
44 #include "TGeoShape.h"
45 #include "TGeoTrd1.h"
46 #include "TGeoTrd2.h"
47 #include "TGeoTube.h"
48 
49 std::tuple<std::shared_ptr<const Acts::CylinderBounds>, const Acts::Transform3,
50  double>
51 Acts::TGeoSurfaceConverter::cylinderComponents(const TGeoShape& tgShape,
52  const Double_t* rotation,
53  const Double_t* translation,
54  const std::string& axes,
55  double scalor) noexcept(false) {
56  std::shared_ptr<const CylinderBounds> bounds = nullptr;
57  Transform3 transform = Transform3::Identity();
58  double thickness = 0.;
59 
60  // Check if it's a tube (segment)
61  auto tube = dynamic_cast<const TGeoTube*>(&tgShape);
62  if (tube != nullptr) {
63  if (not boost::istarts_with(axes, "XY") and
64  not boost::istarts_with(axes, "YX")) {
65  throw std::invalid_argument(
66  "TGeoShape -> CylinderSurface (full): can only be converted with "
67  "'(x/X)(y/Y)(*)' or '(y/Y)(x/X)(*) axes.");
68  }
69 
70  // The sign of the axes
71  int xs = std::islower(axes.at(0)) != 0 ? -1 : 1;
72  int ys = std::islower(axes.at(1)) != 0 ? -1 : 1;
73 
74  // Create translation and rotation
75  Vector3 t(scalor * translation[0], scalor * translation[1],
76  scalor * translation[2]);
77  bool flipxy = not boost::istarts_with(axes, "X");
78  Vector3 ax = flipxy ? xs * Vector3(rotation[1], rotation[4], rotation[7])
79  : xs * Vector3(rotation[0], rotation[3], rotation[6]);
80  Vector3 ay = flipxy ? ys * Vector3(rotation[0], rotation[3], rotation[6])
81  : ys * Vector3(rotation[1], rotation[4], rotation[7]);
82  Vector3 az = ax.cross(ay);
83 
84  double minR = tube->GetRmin() * scalor;
85  double maxR = tube->GetRmax() * scalor;
86  double deltaR = maxR - minR;
87  double medR = 0.5 * (minR + maxR);
88  double halfZ = tube->GetDz() * scalor;
89  if (halfZ > deltaR) {
90  transform = TGeoPrimitivesHelper::makeTransform(ax, ay, az, t);
91  double halfPhi = M_PI;
92  double avgPhi = 0.;
93  // Check if it's a segment
94  auto tubeSeg = dynamic_cast<const TGeoTubeSeg*>(tube);
95  if (tubeSeg != nullptr) {
96  double phi1 = toRadian(tubeSeg->GetPhi1());
97  double phi2 = toRadian(tubeSeg->GetPhi2());
98  if (std::abs(phi2 - phi1) < M_PI * (1. - s_epsilon)) {
99  if (not boost::starts_with(axes, "X")) {
100  throw std::invalid_argument(
101  "TGeoShape -> CylinderSurface (sectorial): can only be "
102  "converted "
103  "with "
104  "'(X)(y/Y)(*)' axes.");
105  }
106  halfPhi = 0.5 * (std::max(phi1, phi2) - std::min(phi1, phi2));
107  avgPhi = 0.5 * (phi1 + phi2);
108  }
109  }
110  bounds = std::make_shared<CylinderBounds>(medR, halfZ, halfPhi, avgPhi);
111  thickness = deltaR;
112  }
113  }
114  return {bounds, transform, thickness};
115 }
116 
117 std::tuple<std::shared_ptr<const Acts::DiscBounds>, const Acts::Transform3,
118  double>
119 Acts::TGeoSurfaceConverter::discComponents(const TGeoShape& tgShape,
120  const Double_t* rotation,
121  const Double_t* translation,
122  const std::string& axes,
123  double scalor) noexcept(false) {
124  using Line2D = Eigen::Hyperplane<double, 2>;
125  std::shared_ptr<const DiscBounds> bounds = nullptr;
126  Transform3 transform = Transform3::Identity();
127 
128  double thickness = 0.;
129  // Special test for composite shape of silicon
130  auto compShape = dynamic_cast<const TGeoCompositeShape*>(&tgShape);
131  if (compShape != nullptr) {
132  if (not boost::istarts_with(axes, "XY")) {
133  throw std::invalid_argument(
134  "TGeoShape -> DiscSurface (Annulus): can only be converted with "
135  "'(x/X)(y/Y)(*)' "
136  "axes");
137  }
138 
139  // Create translation and rotation
140  Vector3 t(scalor * translation[0], scalor * translation[1],
141  scalor * translation[2]);
142  Vector3 ax(rotation[0], rotation[3], rotation[6]);
143  Vector3 ay(rotation[1], rotation[4], rotation[7]);
144  Vector3 az(rotation[2], rotation[5], rotation[8]);
145 
146  transform = TGeoPrimitivesHelper::makeTransform(ax, ay, az, t);
147 
148  auto interNode = dynamic_cast<TGeoIntersection*>(compShape->GetBoolNode());
149  if (interNode != nullptr) {
150  auto baseTube = dynamic_cast<TGeoTubeSeg*>(interNode->GetLeftShape());
151  if (baseTube != nullptr) {
152  double rMin = baseTube->GetRmin() * scalor;
153  double rMax = baseTube->GetRmax() * scalor;
154  auto maskShape = dynamic_cast<TGeoArb8*>(interNode->GetRightShape());
155  if (maskShape != nullptr) {
156  auto maskTransform = interNode->GetRightMatrix();
157  // Get the only vertices
158  const Double_t* polyVrt = maskShape->GetVertices();
159  // Apply the whole transformation stored for the
160  // polyhedron, since there is a translation and
161  // also a side flip that needs to be applied.
162  // @TODO check that 3rd coordinate is not altered by
163  // the transformation ?
164  std::vector<Vector2> vertices;
165  for (unsigned int v = 0; v < 8; v += 2) {
166  std::array<double, 3> local{polyVrt[v + 0], polyVrt[v + 1], 0.};
167  std::array<double, 3> global{};
168  maskTransform->LocalToMaster(local.data(), global.data());
169  Vector2 vtx = Vector2(global[0] * scalor, global[1] * scalor);
170  vertices.push_back(vtx);
171  }
172 
173  std::vector<std::pair<Vector2, Vector2>> boundLines;
174  for (size_t i = 0; i < vertices.size(); ++i) {
175  Vector2 a = vertices.at(i);
176  Vector2 b = vertices.at((i + 1) % vertices.size());
177  Vector2 ab = b - a;
178  double phi = VectorHelpers::phi(ab);
179 
180  if (std::abs(phi) > 3 * M_PI / 4. || std::abs(phi) < M_PI / 4.) {
181  if (a.norm() < b.norm()) {
182  boundLines.push_back(std::make_pair(a, b));
183  } else {
184  boundLines.push_back(std::make_pair(b, a));
185  }
186  }
187  }
188 
189  if (boundLines.size() != 2) {
190  throw std::logic_error(
191  "Input DiscPoly bounds type does not have sensible edges.");
192  }
193 
194  Line2D lA =
195  Line2D::Through(boundLines[0].first, boundLines[0].second);
196  Line2D lB =
197  Line2D::Through(boundLines[1].first, boundLines[1].second);
198  Vector2 ix = lA.intersection(lB);
199 
200  const Eigen::Translation3d originTranslation(ix.x(), ix.y(), 0.);
201  const Vector2 originShift = -ix;
202 
203  // Update transform by prepending the origin shift translation
204  transform = transform * originTranslation;
205  // Transform phi line point to new origin and get phi
206  double phi1 =
207  VectorHelpers::phi(boundLines[0].second - boundLines[0].first);
208  double phi2 =
209  VectorHelpers::phi(boundLines[1].second - boundLines[1].first);
210  double phiMax = std::max(phi1, phi2);
211  double phiMin = std::min(phi1, phi2);
212  double phiShift = 0.;
213 
214  // Create the bounds
215  auto annulusBounds = std::make_shared<const AnnulusBounds>(
216  rMin, rMax, phiMin, phiMax, originShift, phiShift);
217 
218  thickness = maskShape->GetDZ() * scalor;
219  bounds = annulusBounds;
220  }
221  }
222  }
223  } else {
224  // Check if it's a tube
225  auto tube = dynamic_cast<const TGeoTube*>(&tgShape);
226  if (tube != nullptr) {
227  if (not boost::istarts_with(axes, "XY") and
228  not boost::istarts_with(axes, "YX")) {
229  throw std::invalid_argument(
230  "TGeoShape -> DiscSurface: can only be converted with "
231  "'(x/X)(y/Y)(*)' or '(y/Y)(x/X)(*) axes.");
232  }
233 
234  // The sign of the axes
235  int xs = std::islower(axes.at(0)) != 0 ? -1 : 1;
236  int ys = std::islower(axes.at(1)) != 0 ? -1 : 1;
237 
238  // Create translation and rotation
239  Vector3 t(scalor * translation[0], scalor * translation[1],
240  scalor * translation[2]);
241  Vector3 ax = xs * Vector3(rotation[0], rotation[3], rotation[6]);
242  Vector3 ay = ys * Vector3(rotation[1], rotation[4], rotation[7]);
243  Vector3 az = ax.cross(ay);
244  transform = TGeoPrimitivesHelper::makeTransform(ax, ay, az, t);
245 
246  double minR = tube->GetRmin() * scalor;
247  double maxR = tube->GetRmax() * scalor;
248  double halfZ = tube->GetDz() * scalor;
249  double halfPhi = M_PI;
250  double avgPhi = 0.;
251  // Check if it's a segment
252  auto tubeSeg = dynamic_cast<const TGeoTubeSeg*>(tube);
253  if (tubeSeg != nullptr) {
254  double phi1 = toRadian(tubeSeg->GetPhi1());
255  double phi2 = toRadian(tubeSeg->GetPhi2());
256  if (std::abs(phi2 - phi1) < 2 * M_PI * (1. - s_epsilon)) {
257  if (not boost::starts_with(axes, "X")) {
258  throw std::invalid_argument(
259  "TGeoShape -> CylinderSurface (sectorial): can only be "
260  "converted "
261  "with "
262  "'(X)(y/Y)(*)' axes.");
263  }
264  halfPhi = 0.5 * (std::max(phi1, phi2) - std::min(phi1, phi2));
265  avgPhi = 0.5 * (phi1 + phi2);
266  }
267  }
268  bounds = std::make_shared<RadialBounds>(minR, maxR, halfPhi, avgPhi);
269  thickness = 2 * halfZ;
270  }
271  }
272  return {bounds, transform, thickness};
273 }
274 
275 std::tuple<std::shared_ptr<const Acts::PlanarBounds>, const Acts::Transform3,
276  double>
277 Acts::TGeoSurfaceConverter::planeComponents(const TGeoShape& tgShape,
278  const Double_t* rotation,
279  const Double_t* translation,
280  const std::string& axes,
281  double scalor) noexcept(false) {
282  // Create translation and rotation
283  Vector3 t(scalor * translation[0], scalor * translation[1],
284  scalor * translation[2]);
285  Vector3 ax(rotation[0], rotation[3], rotation[6]);
286  Vector3 ay(rotation[1], rotation[4], rotation[7]);
287  Vector3 az(rotation[2], rotation[5], rotation[8]);
288 
289  std::shared_ptr<const PlanarBounds> bounds = nullptr;
290 
291  // Check if it's a box - always true, hence last ressort
292  const TGeoBBox* box = dynamic_cast<const TGeoBBox*>(&tgShape);
293 
294  // Check if it's a trapezoid2
295  const TGeoTrd1* trapezoid1 = dynamic_cast<const TGeoTrd1*>(&tgShape);
296  if ((trapezoid1 != nullptr) and not boost::istarts_with(axes, "XZ")) {
297  throw std::invalid_argument(
298  "TGeoTrd1 -> PlaneSurface: can only be converted with '(x/X)(z/Z)(*)' "
299  "axes");
300  }
301 
302  // Check if it's a trapezoid2
303  const TGeoTrd2* trapezoid2 = dynamic_cast<const TGeoTrd2*>(&tgShape);
304  if (trapezoid2 != nullptr) {
305  if (not boost::istarts_with(axes, "X") and
306  std::abs(trapezoid2->GetDx1() - trapezoid2->GetDx2()) > s_epsilon) {
307  throw std::invalid_argument(
308  "TGeoTrd2 -> PlaneSurface: dx1 must be be equal to dx2 if not taken "
309  "as trapezoidal side.");
310  } else if (not boost::istarts_with(axes, "Y") and
311  std::abs(trapezoid2->GetDy1() - trapezoid2->GetDy2()) >
312  s_epsilon) {
313  throw std::invalid_argument(
314  "TGeoTrd2 -> PlaneSurface: dy1 must be be equal to dy2 if not taken "
315  "as trapezoidal side.");
316  }
317  // Not allowed
318  if (boost::istarts_with(axes, "XY") or boost::istarts_with(axes, "YX")) {
319  throw std::invalid_argument(
320  "TGeoTrd2 -> PlaneSurface: only works with (x/X)(z/Z) and "
321  "(y/Y)(z/Z).");
322  }
323  }
324 
325  // Check if it's a Arb8
326  const TGeoArb8* polygon8c = dynamic_cast<const TGeoArb8*>(&tgShape);
327  TGeoArb8* polygon8 = nullptr;
328  if (polygon8c != nullptr) {
329  // Needed otherwise you can access GetVertices
330  polygon8 = const_cast<TGeoArb8*>(polygon8c);
331  }
332 
333  if ((polygon8c != nullptr) and
334  not(boost::istarts_with(axes, "XY") or boost::istarts_with(axes, "YX"))) {
335  throw std::invalid_argument(
336  "TGeoArb8 -> PlaneSurface: dz must be normal component of Surface.");
337  }
338 
339  // The thickness will be filled
340  double thickness = 0.;
341 
342  // The sign of the axes
343  int xs = std::islower(axes.at(0)) != 0 ? -1 : 1;
344  int ys = std::islower(axes.at(1)) != 0 ? -1 : 1;
345 
346  // Set up the columns : only cyclic iterations are allowed
347  Vector3 cx = xs * ax;
348  Vector3 cy = ys * ay;
349  if (boost::istarts_with(axes, "XY")) {
350  if (trapezoid2 != nullptr) {
351  double dx1 = (ys < 0) ? trapezoid1->GetDx2() : trapezoid1->GetDx1();
352  double dx2 = (ys < 0) ? trapezoid1->GetDx1() : trapezoid1->GetDx2();
353  bounds = std::make_shared<const TrapezoidBounds>(
354  scalor * dx1, scalor * dx2, scalor * trapezoid2->GetDy1());
355  thickness = 2 * scalor * trapezoid2->GetDz();
356  } else if (polygon8 != nullptr) {
357  Double_t* tgverts = polygon8->GetVertices();
358  std::vector<Vector2> pVertices;
359  for (unsigned int ivtx = 0; ivtx < 4; ++ivtx) {
360  pVertices.push_back(Vector2(scalor * xs * tgverts[ivtx * 2],
361  scalor * ys * tgverts[ivtx * 2 + 1]));
362  }
363  bounds = std::make_shared<ConvexPolygonBounds<4>>(pVertices);
364  thickness = 2 * scalor * polygon8->GetDz();
365  } else if (box != nullptr) {
366  bounds = std::make_shared<const RectangleBounds>(scalor * box->GetDX(),
367  scalor * box->GetDY());
368  thickness = 2 * scalor * box->GetDZ();
369  }
370  } else if (boost::istarts_with(axes, "YZ")) {
371  cx = xs * ay;
372  cy = ys * az;
373  if (trapezoid1 != nullptr) {
374  throw std::invalid_argument(
375  "TGeoTrd1 can only be converted with '(x/X)(z/Z)(y/Y)' axes");
376  } else if (trapezoid2 != nullptr) {
377  double dx1 = (ys < 0) ? trapezoid2->GetDy2() : trapezoid2->GetDy1();
378  double dx2 = (ys < 0) ? trapezoid2->GetDy1() : trapezoid2->GetDy2();
379  bounds = std::make_shared<const TrapezoidBounds>(
380  scalor * dx1, scalor * dx2, scalor * trapezoid2->GetDz());
381  thickness = 2 * scalor * trapezoid2->GetDx1();
382  } else if (box != nullptr) {
383  bounds = std::make_shared<const RectangleBounds>(scalor * box->GetDY(),
384  scalor * box->GetDZ());
385  thickness = 2 * scalor * box->GetDX();
386  }
387  } else if (boost::istarts_with(axes, "ZX")) {
388  cx = xs * az;
389  cy = ys * ax;
390  if (box != nullptr) {
391  bounds = std::make_shared<const RectangleBounds>(scalor * box->GetDZ(),
392  scalor * box->GetDX());
393  thickness = 2 * scalor * box->GetDY();
394  }
395  } else if (boost::istarts_with(axes, "XZ")) {
396  cx = xs * ax;
397  cy = ys * az;
398  if (trapezoid1 != nullptr) {
399  double dx1 = (ys < 0) ? trapezoid1->GetDx2() : trapezoid1->GetDx1();
400  double dx2 = (ys < 0) ? trapezoid1->GetDx1() : trapezoid1->GetDx2();
401  bounds = std::make_shared<const TrapezoidBounds>(
402  scalor * dx1, scalor * dx2, scalor * trapezoid1->GetDz());
403  thickness = 2 * scalor * trapezoid1->GetDy();
404  } else if (trapezoid2 != nullptr) {
405  double dx1 = (ys < 0) ? trapezoid2->GetDx2() : trapezoid2->GetDx1();
406  double dx2 = (ys < 0) ? trapezoid2->GetDx1() : trapezoid2->GetDx2();
407  bounds = std::make_shared<const TrapezoidBounds>(
408  scalor * dx1, scalor * dx2, scalor * trapezoid2->GetDz());
409  thickness = 2 * scalor * trapezoid2->GetDy1();
410  } else if (box != nullptr) {
411  bounds = std::make_shared<const RectangleBounds>(scalor * box->GetDX(),
412  scalor * box->GetDZ());
413  thickness = 2 * scalor * box->GetDY();
414  }
415  } else if (boost::istarts_with(axes, "YX")) {
416  cx = xs * ay;
417  cy = ys * ax;
418  if (trapezoid2 != nullptr) {
419  double dx1 = (ys < 0) ? trapezoid2->GetDy2() : trapezoid2->GetDy1();
420  double dx2 = (ys < 0) ? trapezoid2->GetDy1() : trapezoid2->GetDy2();
421  bounds = std::make_shared<const TrapezoidBounds>(
422  scalor * dx1, scalor * dx2, scalor * trapezoid2->GetDx1());
423  thickness = 2 * scalor * trapezoid2->GetDz();
424  } else if (polygon8 != nullptr) {
425  const Double_t* tgverts = polygon8->GetVertices();
426  std::vector<Vector2> pVertices;
427  for (unsigned int ivtx = 0; ivtx < 4; ++ivtx) {
428  pVertices.push_back(Vector2(scalor * xs * tgverts[ivtx * 2 + 1],
429  scalor * ys * tgverts[ivtx * 2]));
430  }
431  bounds = std::make_shared<ConvexPolygonBounds<4>>(pVertices);
432  thickness = 2 * scalor * polygon8->GetDz();
433  } else if (box != nullptr) {
434  bounds = std::make_shared<const RectangleBounds>(scalor * box->GetDY(),
435  scalor * box->GetDX());
436  thickness = 2 * scalor * box->GetDZ();
437  }
438  } else if (boost::istarts_with(axes, "ZY")) {
439  cx = xs * az;
440  cy = ys * ay;
441  if (box != nullptr) {
442  bounds = std::make_shared<const RectangleBounds>(scalor * box->GetDZ(),
443  scalor * box->GetDY());
444  thickness = 2 * scalor * box->GetDX();
445  }
446  } else {
447  throw std::invalid_argument(
448  "TGeoConverter: axes definition must be permutation of "
449  "'(x/X)(y/Y)(z/Z)'");
450  }
451 
452  // Create the normal vector & the transform
453  auto cz = cx.cross(cy);
454  auto transform = TGeoPrimitivesHelper::makeTransform(cx, cy, cz, t);
455 
456  return {bounds, transform, thickness};
457 }
458 
459 std::tuple<std::shared_ptr<Acts::Surface>, Acts::ActsScalar>
460 Acts::TGeoSurfaceConverter::toSurface(const TGeoShape& tgShape,
461  const TGeoMatrix& tgMatrix,
462  const std::string& axes,
463  double scalor) noexcept(false) {
464  // Get the placement and orientation in respect to its mother
465  const Double_t* rotation = tgMatrix.GetRotationMatrix();
466  const Double_t* translation = tgMatrix.GetTranslation();
467 
468  auto [cBounds, cTransform, cThickness] =
469  cylinderComponents(tgShape, rotation, translation, axes, scalor);
470  if (cBounds != nullptr) {
471  return {Surface::makeShared<CylinderSurface>(cTransform, cBounds),
472  cThickness};
473  }
474 
475  auto [dBounds, dTransform, dThickness] =
476  discComponents(tgShape, rotation, translation, axes, scalor);
477  if (dBounds != nullptr) {
478  return {Surface::makeShared<DiscSurface>(dTransform, dBounds), dThickness};
479  }
480 
481  auto [pBounds, pTransform, pThickness] =
482  planeComponents(tgShape, rotation, translation, axes, scalor);
483  if (pBounds != nullptr) {
484  return {Surface::makeShared<PlaneSurface>(pTransform, pBounds), pThickness};
485  }
486 
487  return {nullptr, 0.};
488 }