ATLAS Offline Software
GeoShapeConverter.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 // Trk
6 #include "TrkDetDescrGeoModelCnv/GeoShapeConverter.h"
7 
8 #include <algorithm>
9 #include <cmath>
10 
11 #include "GeoPrimitives/GeoPrimitivesHelpers.h"
12 #include "GeoPrimitives/GeoPrimitivesToStringConverter.h"
13 #include "TrkVolumes/BoundarySurface.h"
14 #include "TrkVolumes/CombinedVolumeBounds.h"
15 #include "TrkVolumes/CuboidVolumeBounds.h"
16 #include "TrkVolumes/CylinderVolumeBounds.h"
17 #include "TrkVolumes/DoubleTrapezoidVolumeBounds.h"
18 #include "TrkVolumes/SimplePolygonBrepVolumeBounds.h"
19 #include "TrkVolumes/SubtractedVolumeBounds.h"
20 #include "TrkVolumes/TrapezoidVolumeBounds.h"
21 #include "TrkVolumes/Volume.h"
22 #include "TrkVolumes/VolumeExcluder.h"
23 // GeoModelKernel
24 #include "GaudiKernel/SystemOfUnits.h"
25 #include "GeoModelKernel/GeoBox.h"
26 #include "GeoModelKernel/GeoCons.h"
27 #include "GeoModelKernel/GeoDefinitions.h"
28 #include "GeoModelKernel/GeoPara.h"
29 #include "GeoModelKernel/GeoPcon.h"
30 #include "GeoModelKernel/GeoPgon.h"
31 #include "GeoModelKernel/GeoShape.h"
32 #include "GeoModelKernel/GeoShapeIntersection.h"
33 #include "GeoModelKernel/GeoShapeShift.h"
34 #include "GeoModelKernel/GeoShapeSubtraction.h"
35 #include "GeoModelKernel/GeoShapeUnion.h"
36 #include "GeoModelKernel/GeoSimplePolygonBrep.h"
37 #include "GeoModelKernel/GeoTrap.h"
38 #include "GeoModelKernel/GeoTrd.h"
39 #include "GeoModelKernel/GeoTube.h"
40 #include "GeoModelKernel/GeoTubs.h"
41 #include "GeoModelKernel/GeoVolumeCursor.h"
42 
43 namespace {
44 // commonly used angles, ±90°, 180°
45 constexpr double p90deg(90.0 * Gaudi::Units::deg),
46  m90deg(-90.0 * Gaudi::Units::deg),
47  p180deg(180.0 * Gaudi::Units::deg);
48 
49 Amg::Transform3D* makeTransform(const Amg::Transform3D& trf) {
50  return std::make_unique<Amg::Transform3D>(trf).release();
51 }
52 } // namespace
53 
54 namespace Trk {
55 GeoShapeConverter::GeoShapeConverter() : AthMessaging{"GeoShapeConverter"} {}
56 
57 std::unique_ptr<CylinderVolumeBounds> GeoShapeConverter::convert(const GeoTubs* gtubs) {
58  // get the dimensions
59  double rMin = gtubs->getRMin();
60  double rMax = gtubs->getRMax();
61  double halflength = gtubs->getZHalfLength();
62  // create the volumeBounds
63  return std::make_unique<CylinderVolumeBounds>(rMin, rMax, halflength);
64 }
65 
66 std::unique_ptr<CylinderVolumeBounds> GeoShapeConverter::convert(const GeoTube* gtube) {
67  // get the dimensions
68  double rMin = gtube->getRMin();
69  double rMax = gtube->getRMax();
70  double halflength = gtube->getZHalfLength();
71  // create the volumeBounds
72  return std::make_unique<CylinderVolumeBounds>(rMin, rMax, halflength);
73 }
74 
75 std::unique_ptr<CylinderVolumeBounds> GeoShapeConverter::convert(const GeoPcon* gpcon,
76  std::vector<double>& zbounds) {
77 
78  // get the pcon igredients ...
79  unsigned int numberOfPlanes = gpcon->getNPlanes();
80  double rMin{10.e10}, rMax{-10.e10}, zMin{10.e10}, zMax{-10.e10};
81 
82  for (unsigned int iplane = 0; iplane < numberOfPlanes; ++iplane) {
83  zMin = std::min(gpcon->getZPlane(iplane), zMin);
84  zMax = std::max(gpcon->getZPlane(iplane), zMax);
85  rMin = std::min(gpcon->getRMinPlane(iplane), rMin);
86  rMax = std::max(gpcon->getRMaxPlane(iplane), rMax);
87  }
88  zbounds.push_back(zMin);
89  zbounds.push_back(zMax);
90 
91  return std::make_unique<CylinderVolumeBounds>(rMin, rMax,
92  0.5 * (zMax - zMin));
93 }
94 
95 std::unique_ptr<CuboidVolumeBounds> GeoShapeConverter::convert(const GeoBox* gbox) {
96  double halfX = gbox->getXHalfLength();
97  double halfY = gbox->getYHalfLength();
98  double halfZ = gbox->getZHalfLength();
99  return std::make_unique<CuboidVolumeBounds>(halfX, halfY, halfZ);
100 }
101 
102 std::unique_ptr<Volume> GeoShapeConverter::translateGeoShape(const GeoShape* sh,
103  const Amg::Transform3D& transf) const {
104  std::unique_ptr<Volume> vol{};
105  double tol = 0.1;
106 
107  ATH_MSG_DEBUG(" translateGeoShape " << sh->type());
108 
109  if (sh->type() == "Trap") {
110  const GeoTrap* trap = dynamic_cast<const GeoTrap*>(sh);
111  std::unique_ptr<TrapezoidVolumeBounds> volBounds{};
112  if (trap->getDxdyndzp() < trap->getDxdyndzn()) {
113  volBounds = std::make_unique<TrapezoidVolumeBounds>(trap->getDxdyndzp(),
114  trap->getDxdyndzn(),
115  trap->getDydzn(),
116  trap->getZHalfLength());
117  } else {
118  volBounds = std::make_unique<TrapezoidVolumeBounds>(trap->getDxdyndzn(),
119  trap->getDxdyndzp(),
120  trap->getDydzn(),
121  trap->getZHalfLength());
122  }
123  return std::make_unique<Volume>(makeTransform(transf), volBounds.release());
124  } else if (sh->type() == "Pgon") {
125  const GeoPgon* pgon = dynamic_cast<const GeoPgon*>(sh);
126  double hlz = 0.5 * std::abs(pgon->getZPlane(1) - pgon->getZPlane(0));
127  double phiH = pgon->getDPhi() / (2. * pgon->getNSides());
128  double hly = 0.5 * std::cos(phiH) * (pgon->getRMaxPlane(0) - pgon->getRMinPlane(0));
129  double dly = 0.5 * std::cos(phiH) * (pgon->getRMaxPlane(0) + pgon->getRMinPlane(0));
130  double hlxmin = pgon->getRMinPlane(0) * std::sin(phiH);
131  double hlxmax = pgon->getRMaxPlane(0) * std::sin(phiH);
132  if (pgon->getDPhi() == 2 * M_PI) {
133 
134  auto volBounds = std::make_unique<CylinderVolumeBounds>(pgon->getRMaxPlane(0), hlz);
135  auto subBounds = std::make_unique<CuboidVolumeBounds>(hlxmax + tol, hlxmax + tol,
136  hlz + tol);
137  auto volume = std::make_unique<Volume>(makeTransform(transf), volBounds.release());
138  auto bVol = std::make_unique<Volume>(nullptr, subBounds.release());
139  const unsigned int nsides(pgon->getNSides());
140  const double twicePhiH(2.0 * phiH);
141  const double xTranslationDistance = hlxmax + std::cos(phiH) * (pgon->getRMaxPlane(0));
142 
143  const Amg::Translation3D xTranslation{xTranslationDistance, 0., 0.};
144  for (unsigned int i = 0; i < nsides; i++) {
145  const double angle = i * twicePhiH;
146  Amg::Transform3D totalTransform = transf * Amg::getRotateZ3D(angle) * xTranslation;
147  auto volS = std::make_unique<Volume>(*bVol, totalTransform);
148  auto combBounds =std::make_unique<SubtractedVolumeBounds>(volume.release(),
149  volS.release());
150  volume = std::make_unique<Volume>(nullptr, combBounds.release());
151  }
152  return volume;
153  }
154 
155  if (pgon->getNSides() == 1) {
156  auto volBounds = std::make_unique<TrapezoidVolumeBounds>(hlxmin, hlxmax, hly, hlz);
157  Amg::Transform3D totalTransform = transf *
158  Amg::getRotateZ3D(m90deg) *
159  Amg::Translation3D(0., dly, 0.);
160  return std::make_unique<Volume>(makeTransform(totalTransform),
161  volBounds.release());
162  }
163 
164  if (pgon->getNSides() == 2) {
165  auto cylBounds = std::make_unique<CylinderVolumeBounds>(0, dly + hly, hlz);
166  return std::make_unique<Volume>(makeTransform(transf),
167  cylBounds.release());
168  }
169  return vol;
170  } else if (sh->type() == "Trd") {
171  const GeoTrd* trd = dynamic_cast<const GeoTrd*>(sh);
172  //
173  double x1 = trd->getXHalfLength1();
174  double x2 = trd->getXHalfLength2();
175  double y1 = trd->getYHalfLength1();
176  double y2 = trd->getYHalfLength2();
177  double z = trd->getZHalfLength();
178  //
179  ATH_MSG_DEBUG(" Trd x1 " << x1 << " x2 " << x2 << " y1 " << y1 << " y2 "
180  << y2 << " z " << z);
181  // Note this flip comes from the y axis in Tracking -> z axis in
182  // Geomodel
183  if (y1 == y2) {
184  if (x1 <= x2) {
185  auto volBounds = std::make_unique<TrapezoidVolumeBounds>(x1, x2, z, y1);
186  Amg::Transform3D totalTransform = transf * Amg::getRotateX3D(p90deg);
187  ATH_MSG_DEBUG(" Trd new volume case 1 Trapezoid minHalflengthX "
188  << volBounds->minHalflengthX()
189  << " maxHalflengthX() "
190  << volBounds->maxHalflengthX());
191  vol = std::make_unique<Volume>(makeTransform(totalTransform),
192  volBounds.release());
193 
194 
195  } else {
196 
197  auto volBounds = std::make_unique<TrapezoidVolumeBounds>(x2, x1, z, y1);
198  Amg::Transform3D totalTransform = transf *
199  Amg::getRotateY3D(p180deg) *
200  Amg::getRotateZ3D(p180deg);
201 
202  if (msgLvl(MSG::DEBUG)) {
203  const Amg::Vector3D top{-x1, y1, z}, bottom{-x2, y1, -z},
204  top2{x1, y1, z}, bottom2{x2, y1, -z};
205  ATH_MSG_DEBUG(" Trd new volume case 2 Trapezoid minHalflengthX "
206  << volBounds->minHalflengthX() << " maxHalflengthX() "
207  << volBounds->maxHalflengthX());
208  ATH_MSG_DEBUG(" Original topLocal " << Amg::toString(top) << " Radius " << top.perp());
209  ATH_MSG_DEBUG(" Original bottomLocal "<< Amg::toString(bottom) << " Radius " << bottom.perp());
210  Amg::Vector3D topG{transf * top}, bottomG{transf * bottom};
211  ATH_MSG_DEBUG(" top Global " << Amg::toString(topG)<< " Radius " << topG.perp());
212  ATH_MSG_DEBUG(" bottom Global " << Amg::toString(bottomG)<< " Radius "<< bottomG.perp());
213  topG = transf * top2;
214  bottomG = transf * bottom2;
215  ATH_MSG_DEBUG(" top2 Global x " << Amg::toString(topG)<< " Radius " << topG.perp());
216  ATH_MSG_DEBUG(" bottom2 Global " << Amg::toString(bottomG) << " Radius: " << bottomG.perp());
217  const Amg::Vector3D topR{-x2, z, y1}, bottomR{-x1, -z, y1},
218  top2R{x2, z, y1}, bottom2R{x1, -z, y1};
219  topG = totalTransform * topR;
220  bottomG = totalTransform * bottomR;
221  ATH_MSG_DEBUG(" topR Global " << Amg::toString(topG)<< " Radius " << topG.perp());
222  ATH_MSG_DEBUG(" bottomR Global " << Amg::toString(bottomG)<< " Radius "<< bottomG.perp());
223  topG = totalTransform * top2R;
224  bottomG = totalTransform * bottom2R;
225  ATH_MSG_DEBUG(" top2R Global " << Amg::toString(topG) << " Radius "<< topG.perp());
226  ATH_MSG_DEBUG(" bottom2R Global " << Amg::toString(bottomG)<< " Radius "<< bottomG.perp());
227 
228  ATH_MSG_DEBUG(" Original bottomLocal "<< Amg::toString(bottom) << " Radius "<< bottom.perp());
229  ATH_MSG_DEBUG(" topLocal x " << Amg::toString(topR)<< " Radius " << topR.perp());
230  topG = Amg::getRotateY3D(p180deg) * topR;
231  ATH_MSG_DEBUG(" topLocal Y Flip " << Amg::toString(topG)<< " Radius " << topG.perp());
232  topG = Amg::getRotateY3D(p180deg) * topR;
233  ATH_MSG_DEBUG(" topLocal X Flip " << Amg::toString(topG) << " Radius "<< topG.perp());
234  topG = Amg::getRotateZ3D(p90deg) * topR;
235  ATH_MSG_DEBUG(" topLocal XY Flip " << Amg::toString(topG) << " Radius " << topG.perp());
236  topG = Amg::getRotateY3D(p180deg) * topR;
237  ATH_MSG_DEBUG(" topLocal YZ Flip " << Amg::toString(topG) << " Radius " << topG.perp());
238  topG = Amg::getRotateY3D(p90deg) * topR;
239  ATH_MSG_DEBUG(" topLocal XZ Flip " << Amg::toString(topG) << " Radius " << topG.perp());
240  topG = Amg::getRotateY3D(p180deg) * topR;
241  ATH_MSG_DEBUG(" topLocal XY sign Flip x "<< Amg::toString(topG) << " Radius " << topG.perp());
242  }
243  vol = std::make_unique<Volume>(makeTransform(totalTransform),
244  volBounds.release());
245  }
246  return vol;
247  } else if (x1 == x2) {
248  if (y1 < y2) {
249  std::unique_ptr<TrapezoidVolumeBounds> volBounds =
250  std::make_unique<TrapezoidVolumeBounds>(y1, y2, z, x1);
251  ATH_MSG_DEBUG(" Trd new volume case 3 Trapezoid minHalflengthX "
252  << volBounds->minHalflengthX()<< " maxHalflengthX() "
253  << volBounds->maxHalflengthX());
254  Amg::Transform3D totalTransform = transf *
255  Amg::getRotateZ3D(p90deg) *
256  Amg::getRotateX3D(p90deg);
257  vol = std::make_unique<Volume>(makeTransform(totalTransform),
258  volBounds.release());
259 
260  } else {
261  auto volBounds = std::make_unique<TrapezoidVolumeBounds>(y2, y1, z, x1);
262  ATH_MSG_DEBUG(" Trd new volume case 4 Trapezoid minHalflengthX "
263  << volBounds->minHalflengthX()
264  << " maxHalflengthX() "<< volBounds->maxHalflengthX());
265 
266  Amg::Transform3D totalTransform = transf *
267  Amg::getRotateX3D(p180deg) *
268  Amg::getRotateZ3D(p90deg) *
269  Amg::getRotateX3D(p90deg);
270  vol = std::make_unique<Volume>(makeTransform(totalTransform),
271  volBounds.release());
272  }
273  return vol;
274  } else {
275  ATH_MSG_WARNING("PROBLEM: translating trapezoid: not recognized:"
276  << x1 << "," << x2 << "," << y1 << "," << y2 << ","<< z);
277  }
278  } else if (sh->type() == "Box") {
279  const GeoBox* box = dynamic_cast<const GeoBox*>(sh);
280  //
281  double x = box->getXHalfLength();
282  double y = box->getYHalfLength();
283  double z = box->getZHalfLength();
284  std::unique_ptr<CuboidVolumeBounds> volBounds = std::make_unique<CuboidVolumeBounds>(x, y, z);
285  return std::make_unique<Volume>(makeTransform(transf), volBounds.release());
286  } else if (sh->type() == "Para") {
287  const GeoPara* para = dynamic_cast<const GeoPara*>(sh);
288  //
289  double x = para->getXHalfLength();
290  double y = para->getYHalfLength();
291  double z = para->getZHalfLength();
292  auto volBounds = std::make_unique<CuboidVolumeBounds>(x, y, z);
293  return std::make_unique<Volume>(makeTransform(transf), volBounds.release());
294  //
295  } else if (sh->type() == "Tube") {
296  const GeoTube* tube = dynamic_cast<const GeoTube*>(sh);
297  return std::make_unique<Volume>(makeTransform(transf), convert(tube).release());
298  } else if (sh->type() == "Tubs") { // non-trivial case - transform!
299  const GeoTubs* tubs = dynamic_cast<const GeoTubs*>(sh);
300  double rMin = tubs->getRMin();
301  double rMax = tubs->getRMax();
302  double z = tubs->getZHalfLength();
303  double aPhi = tubs->getSPhi();
304  double dPhi = tubs->getDPhi();
305  auto volBounds = std::make_unique<CylinderVolumeBounds>(rMin, rMax, 0.5 * dPhi, z);
306  Amg::Transform3D totalTransform(transf * Amg::getRotateZ3D(aPhi + 0.5 * dPhi));
307  return std::make_unique<Volume>(makeTransform(totalTransform), volBounds.release());
308  } else if (sh->type() == "Cons") {
309  const GeoCons* cons = dynamic_cast<const GeoCons*>(sh);
310  double rMin1 = cons->getRMin1();
311  double rMin2 = cons->getRMin2();
312  double rMax1 = cons->getRMax1();
313  double rMax2 = cons->getRMax2();
314  double z = cons->getDZ();
315  double aPhi = cons->getSPhi();
316  double dPhi = cons->getDPhi();
317  // translate into tube with average radius
318  if (dPhi == 2 * M_PI) {
319  auto volBounds =std::make_unique<CylinderVolumeBounds>(0.5 * (rMin1 + rMin2),
320  0.5 * (rMax1 + rMax2), z);
321  return std::make_unique<Volume>(makeTransform(transf),
322  volBounds.release());
323  } else {
324  auto volBounds =std::make_unique<CylinderVolumeBounds>(0.5 * (rMin1 + rMin2),
325  0.5 * (rMax1 + rMax2),
326  0.5 * dPhi, z);
327  Amg::Transform3D totalTransform = transf * Amg::getRotateZ3D(aPhi + 0.5 * dPhi);
328  return std::make_unique<Volume>(makeTransform(totalTransform),
329  volBounds.release());
330  }
331  } else if (sh->type() == "Pcon") {
332  const GeoPcon* con = dynamic_cast<const GeoPcon*>(sh);
333  std::unique_ptr<CylinderVolumeBounds> volBounds{};
334  double aPhi = con->getSPhi();
335  double dPhi = con->getDPhi();
336  double z1 = con->getZPlane(0);
337  double r1 = con->getRMinPlane(0);
338  double R1 = con->getRMaxPlane(0);
339  std::vector<std::unique_ptr<Volume>> cyls;
340  const unsigned int nPlanes = con->getNPlanes();
341  ATH_MSG_DEBUG(" convert pcon aPhi " << aPhi << " dPhi " << dPhi << " z1 " << z1 << " r1 "
342  << r1 << " R1 " << R1 << " nPlanes " << nPlanes);
343  for (unsigned int iv = 1; iv < nPlanes; iv++) {
344  double z2 = con->getZPlane(iv);
345  double r2 = con->getRMinPlane(iv);
346  double R2 = con->getRMaxPlane(iv);
347  double zshift = 0.5 * (z1 + z2);
348  double hz = 0.5 * std::abs(z1 - z2);
349  double rmin = std::max(r1, r2);
350  double rmax = std::sqrt((R1 * R1 + R1 * R2 + R2 * R2 - r1 * r1 - r1 * r2 - r2 * r2) / 3 + rmin * rmin);
351  ATH_MSG_DEBUG(" iPlane " << iv << " z2 " << z2 << " r2 " << r2 << " R2 " << R2 << " zshift " << zshift
352  << " hz " << hz << " rmin " << rmin << " rmax " << rmax);
353  double dz = con->getZPlane(iv) - con->getZPlane(iv - 1);
354  double drMin = con->getRMinPlane(iv) - con->getRMinPlane(iv - 1);
355  double drMax = con->getRMaxPlane(iv) - con->getRMaxPlane(iv - 1);
356  int nSteps = 1;
357  if (std::abs(dz) > 1 && (std::abs(drMin) > 0.1 * std::abs(dz) ||
358  std::abs(drMax) > 0.1 * std::abs(dz))) {
359  double dMax = std::abs(dz);
360  dMax = std::max(std::abs(drMin), dMax);
361  dMax = std::max(std::abs(drMax), dMax);
362  nSteps = std::clamp(dMax / 50., 2., 20.);
363  ATH_MSG_DEBUG(" Now "
364  << nSteps << " cylinders should be created " << dz
365  << " drMin " << drMin << " drMax " << drMax
366  << " splopeMin " << drMin / dz << " slopeMax "
367  << drMax / dz);
368  }
369  for (int j = 0; j < nSteps; j++) {
371  double zStep = (0.5 + j) * dz / nSteps;
372  if (nSteps > 1) {
373  hz = 0.5 * std::abs(z1 - z2) / nSteps;
374  zshift = z1 + zStep;
375  rmin = r1 + drMin * zStep / dz;
376  rmax = R1 + drMax * zStep / dz;
377  }
378  ATH_MSG_DEBUG(" cylinder " << j << " zshift " << zshift
379  << " rmin " << rmin << " rmax "
380  << rmax << " hz " << hz);
381  // translate into tube sector
382  if (dPhi == 2 * M_PI) {
383  volBounds = std::make_unique<CylinderVolumeBounds>(rmin, rmax, hz);
384  Amg::Transform3D totalTransform = transf * Amg::Translation3D{0., 0., zshift};
385  cyls.emplace_back(std::make_unique<Volume>(makeTransform(totalTransform),
386  volBounds.release()));
387  } else {
388  volBounds = std::make_unique<CylinderVolumeBounds>(rmin, rmax, 0.5 * dPhi, hz);
389  Amg::Transform3D totalTransform = transf * Amg::Translation3D{0., 0., zshift} *
390  Amg::getRotateZ3D(aPhi + 0.5 * dPhi);
391  cyls.emplace_back(std::make_unique<Volume>(makeTransform(totalTransform),
392  volBounds.release()));
393  }
394  } // end loop over steps
395  z1 = z2;
396  r1 = r2;
397  R1 = R2;
398  }
399 
400  if (cyls.size() < 2) {
401  return std::move(cyls[0]);
402  } else {
403  auto comb =std::make_unique<CombinedVolumeBounds>(cyls[0].release(), cyls[1].release(), false);
404  std::unique_ptr<Volume> combVol = std::make_unique<Volume>(nullptr, comb.release());
405  for (unsigned int ic = 2; ic < cyls.size(); ++ic) {
406  comb = std::make_unique<CombinedVolumeBounds>(combVol.release(), cyls[ic].release(), false);
407  combVol = std::make_unique<Volume>(nullptr, comb.release());
408  }
409  return combVol;
410  }
411  } else if (sh->type() == "SimplePolygonBrep") {
412  const GeoSimplePolygonBrep* spb = dynamic_cast<const GeoSimplePolygonBrep*>(sh);
413  unsigned int nv = spb->getNVertices();
414  std::vector<std::pair<double, double>> ivtx(nv);
415  for (unsigned int iv = 0; iv < nv; iv++) {
416  ivtx[iv] = std::make_pair(spb->getXVertex(iv), spb->getYVertex(iv));
417  ATH_MSG_DEBUG(" SimplePolygonBrep x "<< spb->getXVertex(iv) << " y " << spb->getYVertex(iv)
418  << " z " << spb->getDZ());
419  }
420  // translate into trapezoid or double trapezoid if possible
421  if (nv == 4 || nv == 6) {
422  std::vector<double> xstep;
423  std::vector<std::pair<double, double>> ystep;
424  bool trdlike = true;
425  for (unsigned int iv = 0; iv < nv; ++iv) {
426  if (ystep.empty() || spb->getYVertex(iv) > ystep.back().first)
427  ystep.emplace_back(spb->getYVertex(iv),
428  std::abs(spb->getXVertex(iv)));
429  else {
430  std::vector<std::pair<double, double>>::iterator iy = ystep.begin();
431  while (iy + 1 < ystep.end() &&
432  spb->getYVertex(iv) > (*iy).first + 1.e-3) {
433  ++iy;
434  }
435  if (spb->getYVertex(iv) < (*iy).first - 1.e-3) {
436  ystep.insert(iy, std::make_pair(spb->getYVertex(iv), std::abs(spb->getXVertex(iv))));
437  } else if (spb->getYVertex(iv) == (*iy).first &&
438  std::abs(spb->getXVertex(iv)) != (*iy).second) {
439  trdlike = false;
440  }
441  }
442  }
443 
444  if (trdlike) {
445  std::unique_ptr<VolumeBounds> volBounds{};
446  if (nv == 4) {
447  if (ystep[1].second >= ystep[0].second) { // expected ordering
448  volBounds = std::make_unique<TrapezoidVolumeBounds>(ystep[0].second,
449  ystep[1].second,
450  0.5 * (ystep[1].first - ystep[0].first),
451  spb->getDZ());
452  return std::make_unique<Volume>(makeTransform(transf),
453  volBounds.release());
454  }
455  }
456 
457  if (nv == 6) {
458  if (ystep[1].second >= ystep[0].second &&
459  ystep[2].second >= ystep[1].second) { // expected ordering
460  volBounds = std::make_unique<DoubleTrapezoidVolumeBounds>(ystep[0].second,
461  ystep[1].second,
462  ystep[2].second,
463  0.5 * (ystep[1].first - ystep[0].first),
464  0.5 * (ystep[2].first - ystep[1].first),
465  spb->getDZ());
466  return std::make_unique<Volume>(makeTransform(transf * Amg::Translation3D(0., ystep[1].first, 0.)),
467  volBounds.release());
468  }
469  }
470  } // not trd-like
471  }
472  auto newBounds = std::make_unique<SimplePolygonBrepVolumeBounds>(ivtx, spb->getDZ());
473  return std::make_unique<Volume>(makeTransform(transf),
474  newBounds.release());
475  }
476 
477  else if (sh->type() == "Subtraction") {
478  const GeoShapeSubtraction* sub = dynamic_cast<const GeoShapeSubtraction*>(sh);
479 
480  const GeoShape* shA = sub->getOpA();
481  const GeoShape* shB = sub->getOpB();
482  std::unique_ptr<Volume> volA = translateGeoShape(shA, transf);
483  std::unique_ptr<Volume> volB = translateGeoShape(shB, transf);
484  auto volBounds = std::make_unique<SubtractedVolumeBounds>(volA.release(),
485  volB.release());
486  return std::make_unique<Volume>(nullptr, volBounds.release());
487  } else if (sh->type() == "Union") {
488  const GeoShapeUnion* uni = dynamic_cast<const GeoShapeUnion*>(sh);
489  const GeoShape* shA = uni->getOpA();
490  const GeoShape* shB = uni->getOpB();
491  std::unique_ptr<Volume> volA = translateGeoShape(shA, transf);
492  std::unique_ptr<Volume> volB = translateGeoShape(shB, transf);
493  auto volBounds = std::make_unique<CombinedVolumeBounds>(volA.release(),
494  volB.release(), false);
495  return std::make_unique<Volume>(nullptr, volBounds.release());
496  } else if (sh->type() == "Intersection") {
497  const GeoShapeIntersection* intersect = dynamic_cast<const GeoShapeIntersection*>(sh);
498 
499  const GeoShape* shA = intersect->getOpA();
500  const GeoShape* shB = intersect->getOpB();
501  std::unique_ptr<Volume> volA{translateGeoShape(shA, transf)};
502  std::unique_ptr<Volume> volB{translateGeoShape(shB, transf)};
503  auto volBounds = std::make_unique<CombinedVolumeBounds>(volA.release(),
504  volB.release(), true);
505  return std::make_unique<Volume>(nullptr, volBounds.release());
506  }
507 
508  if (sh->type() == "Shift") {
509  const GeoShapeShift* shift = dynamic_cast<const GeoShapeShift*>(sh);
510  return translateGeoShape(shift->getOp(), transf * shift->getX());
511  }
512  ATH_MSG_WARNING("shape " << sh->type() << " not recognized, return 0");
513  return nullptr;
514 }
515 
516 void GeoShapeConverter::decodeShape(const GeoShape* sh) const {
517  ATH_MSG_DEBUG(" ");
518  ATH_MSG_DEBUG("decoding shape:" << sh->type());
519 
520  if (sh->type() == "Pgon") {
521  const GeoPgon* pgon = dynamic_cast<const GeoPgon*>(sh);
522  if (pgon)
523  ATH_MSG_DEBUG("polygon: " << pgon->getNPlanes() << " planes "
524  << pgon->getSPhi() << " "
525  << pgon->getDPhi() << " "
526  << pgon->getNSides());
527  else
528  ATH_MSG_DEBUG("polygon: WARNING: dynamic_cast failed!");
529  }
530 
531  if (sh->type() == "Trd") {
532  const GeoTrd* trd = dynamic_cast<const GeoTrd*>(sh);
533  ATH_MSG_DEBUG("dimensions:" << trd->getXHalfLength1() << ","
534  << trd->getXHalfLength2() << ","
535  << trd->getYHalfLength1() << ","
536  << trd->getYHalfLength2() << ","
537  << trd->getZHalfLength());
538  }
539  if (sh->type() == "Box") {
540  const GeoBox* box = dynamic_cast<const GeoBox*>(sh);
541  ATH_MSG_DEBUG("dimensions:" << box->getXHalfLength() << ","
542  << box->getYHalfLength() << ","
543  << box->getZHalfLength());
544  }
545 
546  if (sh->type() == "Tube") {
547  const GeoTube* tube = dynamic_cast<const GeoTube*>(sh);
548  ATH_MSG_DEBUG("dimensions:" << tube->getRMin() << "," << tube->getRMax()
549  << "," << tube->getZHalfLength());
550  }
551 
552  if (sh->type() == "Tubs") {
553  const GeoTubs* tubs = dynamic_cast<const GeoTubs*>(sh);
554  ATH_MSG_DEBUG("dimensions:" << tubs->getRMin() << "," << tubs->getRMax()
555  << "," << tubs->getZHalfLength() << ","
556  << tubs->getSPhi() << ","
557  << tubs->getDPhi());
558  }
559 
560  if (sh->type() == "Cons") {
561  const GeoCons* cons = dynamic_cast<const GeoCons*>(sh);
562  ATH_MSG_DEBUG("dimensions:"
563  << cons->getRMin1() << "," << cons->getRMin2() << ","
564  << cons->getRMax1() << "," << cons->getRMax2() << ","
565  << cons->getDZ() << "," << cons->getSPhi() << ","
566  << cons->getDPhi());
567  }
568 
569  if (sh->type() == "Subtraction") {
570  const GeoShapeSubtraction* sub =
571  dynamic_cast<const GeoShapeSubtraction*>(sh);
572  const GeoShape* sha = sub->getOpA();
573  const GeoShape* shs = sub->getOpB();
574  ATH_MSG_DEBUG("decoding subtracted shape:");
575  decodeShape(sha);
576  decodeShape(shs);
577  }
578 
579  if (sh->type() == "Union") {
580  const GeoShapeUnion* sub = dynamic_cast<const GeoShapeUnion*>(sh);
581  const GeoShape* shA = sub->getOpA();
582  const GeoShape* shB = sub->getOpB();
583  ATH_MSG_DEBUG("decoding shape A:");
584  decodeShape(shA);
585  ATH_MSG_DEBUG("decoding shape B:");
586  decodeShape(shB);
587  }
588  if (sh->type() == "Shift") {
589  const GeoShapeShift* shift = dynamic_cast<const GeoShapeShift*>(sh);
590  const GeoShape* shA = shift->getOp();
591  const GeoTrf::Transform3D& transf = shift->getX();
592  ATH_MSG_DEBUG("shifted by:transl:"
593  << transf.translation() << ", rot:" << transf(0, 0) << ","
594  << transf(1, 1) << "," << transf(2, 2));
595  decodeShape(shA);
596  }
597 }
598 } // namespace Trk
xAOD::iterator
JetConstituentVector::iterator iterator
Definition: JetConstituentVector.cxx:68
Trk::y
@ y
Definition: ParamDefs.h:56
Muon::makeTransform
Amg::Transform3D * makeTransform(const Amg::Transform3D &trf)
Definition: MuonSpectrometer/MuonDetDescr/MuonTrackingGeometry/MuonTrackingGeometry/Utils.h:14
python.SystemOfUnits.second
int second
Definition: SystemOfUnits.py:120
plotBeamSpotCompare.x1
x1
Definition: plotBeamSpotCompare.py:216
Trk::z
@ z
global position (cartesian)
Definition: ParamDefs.h:57
max
constexpr double max()
Definition: ap_fixedTest.cxx:33
min
constexpr double min()
Definition: ap_fixedTest.cxx:26
plotBeamSpotCompare.x2
x2
Definition: plotBeamSpotCompare.py:218
M_PI
#define M_PI
Definition: ActiveFraction.h:11
deg
#define deg
Definition: SbPolyhedron.cxx:17
drawFromPickle.cos
cos
Definition: drawFromPickle.py:36
Trk::GeoShapeConverter::GeoShapeConverter
GeoShapeConverter()
Definition: GeoShapeConverter.cxx:55
makeTRTBarrelCans.y1
tuple y1
Definition: makeTRTBarrelCans.py:15
TrapezoidVolumeBounds.h
Amg::getRotateZ3D
Amg::Transform3D getRotateZ3D(double angle)
get a rotation transformation around Z-axis
Definition: GeoPrimitivesHelpers.h:270
VolumeExcluder.h
CylinderVolumeBounds.h
Volume.h
Amg::getRotateX3D
Amg::Transform3D getRotateX3D(double angle)
get a rotation transformation around X-axis
Definition: GeoPrimitivesHelpers.h:252
Amg::toString
std::string toString(const Translation3D &translation, int precision=4)
GeoPrimitvesToStringConverter.
Definition: GeoPrimitivesToStringConverter.h:40
AthMessaging::msgLvl
bool msgLvl(const MSG::Level lvl) const
Test the output level.
Definition: AthMessaging.h:151
FullCPAlgorithmsTest_eljob.sh
sh
Definition: FullCPAlgorithmsTest_eljob.py:111
lumiFormat.i
int i
Definition: lumiFormat.py:85
RCU::Shell
Definition: ShellExec.cxx:28
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
angle
double angle(const GeoTrf::Vector2D &a, const GeoTrf::Vector2D &b)
Definition: TRTDetectorFactory_Full.cxx:73
makeTRTBarrelCans.y2
tuple y2
Definition: makeTRTBarrelCans.py:18
TauGNNUtils::Variables::Track::dPhi
bool dPhi(const xAOD::TauJet &tau, const xAOD::TauTrack &track, double &out)
Definition: TauGNNUtils.cxx:538
Amg::Transform3D
Eigen::Affine3d Transform3D
Definition: GeoPrimitives.h:46
AthMessaging
Class to provide easy MsgStream access and capabilities.
Definition: AthMessaging.h:55
checkCorrelInHIST.nSteps
int nSteps
Definition: checkCorrelInHIST.py:458
Trk::GeoShapeConverter::decodeShape
void decodeShape(const GeoShape *) const
Decode and dump arbitrary GeoShape for visual inspection.
Definition: GeoShapeConverter.cxx:516
grepfile.ic
int ic
Definition: grepfile.py:33
Trk
Ensure that the ATLAS eigen extensions are properly loaded.
Definition: FakeTrackBuilder.h:9
SimplePolygonBrepVolumeBounds.h
DoubleTrapezoidVolumeBounds.h
python.EventInfoMgtInit.release
release
Definition: EventInfoMgtInit.py:24
Amg::getRotateY3D
Amg::Transform3D getRotateY3D(double angle)
get a rotation transformation around Y-axis
Definition: GeoPrimitivesHelpers.h:261
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition: GeoPrimitives.h:47
BoundarySurface.h
Trk::GeoShapeConverter::translateGeoShape
std::unique_ptr< Volume > translateGeoShape(const GeoShape *shape, const Amg::Transform3D &trf) const
Convert an arbitrary GeoShape into Trk::Volume.
Definition: GeoShapeConverter.cxx:102
Trk::GeoShapeConverter::convert
static std::unique_ptr< CylinderVolumeBounds > convert(const GeoTubs *gtub)
Convert a tubs.
Definition: GeoShapeConverter.cxx:57
Amg::intersect
std::optional< double > intersect(const AmgVector(N)&posA, const AmgVector(N)&dirA, const AmgVector(N)&posB, const AmgVector(N)&dirB)
Calculates the point B' along the line B that's closest to a second line A.
Definition: GeoPrimitivesHelpers.h:347
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
GeoPrimitivesHelpers.h
DeMoScan.first
bool first
Definition: DeMoScan.py:536
DEBUG
#define DEBUG
Definition: page_access.h:11
Amg::Translation3D
Eigen::Translation< double, 3 > Translation3D
Definition: GeoPrimitives.h:44
GeoPrimitivesToStringConverter.h
GeoShapeConverter.h
CombinedVolumeBounds.h
top
@ top
Definition: TruthClasses.h:64
fitman.hz
def hz
Definition: fitman.py:516
drawFromPickle.sin
sin
Definition: drawFromPickle.py:36
CaloLCW_tf.trf
trf
Definition: CaloLCW_tf.py:20
Trk::x
@ x
Definition: ParamDefs.h:55
calibdata.tube
tube
Definition: calibdata.py:31
CuboidVolumeBounds.h
SubtractedVolumeBounds.h
Generated on Sun Dec 22 2024 21:10:51 for ATLAS Offline Software by   doxygen 1.8.18