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