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