104 std::unique_ptr<Volume> vol{};
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(),
116 trap->getZHalfLength());
118 volBounds = std::make_unique<TrapezoidVolumeBounds>(trap->getDxdyndzn(),
121 trap->getZHalfLength());
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) {
134 auto volBounds = std::make_unique<CylinderVolumeBounds>(pgon->getRMaxPlane(0), hlz);
135 auto subBounds = std::make_unique<CuboidVolumeBounds>(hlxmax + tol, hlxmax + 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));
144 for (
unsigned int i = 0;
i < nsides;
i++) {
145 const double angle =
i * twicePhiH;
147 auto volS = std::make_unique<Volume>(*bVol, totalTransform);
148 auto combBounds =std::make_unique<SubtractedVolumeBounds>(volume.release(),
150 volume = std::make_unique<Volume>(
nullptr, combBounds.release());
155 if (pgon->getNSides() == 1) {
156 auto volBounds = std::make_unique<TrapezoidVolumeBounds>(hlxmin, hlxmax, hly, hlz);
160 return std::make_unique<Volume>(
makeTransform(totalTransform),
161 volBounds.release());
164 if (pgon->getNSides() == 2) {
165 auto cylBounds = std::make_unique<CylinderVolumeBounds>(0, dly + hly, hlz);
167 cylBounds.release());
170 }
else if (
sh->type() ==
"Trd") {
171 const GeoTrd* trd =
dynamic_cast<const GeoTrd*
>(
sh);
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();
180 <<
y2 <<
" z " <<
z);
185 auto volBounds = std::make_unique<TrapezoidVolumeBounds>(
x1,
x2,
z,
y1);
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());
197 auto volBounds = std::make_unique<TrapezoidVolumeBounds>(
x2,
x1,
z,
y1);
205 ATH_MSG_DEBUG(
" Trd new volume case 2 Trapezoid minHalflengthX "
206 << volBounds->minHalflengthX() <<
" maxHalflengthX() "
207 << volBounds->maxHalflengthX());
213 topG = transf * top2;
214 bottomG = transf * bottom2;
219 topG = totalTransform * topR;
220 bottomG = totalTransform * bottomR;
223 topG = totalTransform * top2R;
224 bottomG = totalTransform * bottom2R;
243 vol = std::make_unique<Volume>(
makeTransform(totalTransform),
244 volBounds.release());
247 }
else if (
x1 ==
x2) {
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());
257 vol = std::make_unique<Volume>(
makeTransform(totalTransform),
258 volBounds.release());
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());
270 vol = std::make_unique<Volume>(
makeTransform(totalTransform),
271 volBounds.release());
276 <<
x1 <<
"," <<
x2 <<
"," <<
y1 <<
"," <<
y2 <<
","<<
z);
278 }
else if (
sh->type() ==
"Box") {
279 const GeoBox* box =
dynamic_cast<const GeoBox*
>(
sh);
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);
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());
295 }
else if (
sh->type() ==
"Tube") {
296 const GeoTube*
tube =
dynamic_cast<const GeoTube*
>(
sh);
298 }
else if (
sh->type() ==
"Tubs") {
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);
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();
319 auto volBounds =std::make_unique<CylinderVolumeBounds>(0.5 * (rMin1 + rMin2),
320 0.5 * (rMax1 + rMax2),
z);
322 volBounds.release());
324 auto volBounds =std::make_unique<CylinderVolumeBounds>(0.5 * (rMin1 + rMin2),
325 0.5 * (rMax1 + rMax2),
328 return std::make_unique<Volume>(
makeTransform(totalTransform),
329 volBounds.release());
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);
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);
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.);
364 <<
nSteps <<
" cylinders should be created " << dz
365 <<
" drMin " << drMin <<
" drMax " << drMax
366 <<
" splopeMin " << drMin / dz <<
" slopeMax "
369 for (
int j = 0; j <
nSteps; j++) {
371 double zStep = (0.5 + j) * dz /
nSteps;
373 hz = 0.5 * std::abs(z1 - z2) /
nSteps;
375 rmin =
r1 + drMin * zStep / dz;
376 rmax = R1 + drMax * zStep / dz;
379 <<
" rmin " << rmin <<
" rmax "
380 << rmax <<
" hz " <<
hz);
383 volBounds = std::make_unique<CylinderVolumeBounds>(rmin, rmax,
hz);
385 cyls.emplace_back(std::make_unique<Volume>(
makeTransform(totalTransform),
386 volBounds.release()));
388 volBounds = std::make_unique<CylinderVolumeBounds>(rmin, rmax, 0.5 *
dPhi,
hz);
391 cyls.emplace_back(std::make_unique<Volume>(
makeTransform(totalTransform),
392 volBounds.release()));
400 if (cyls.size() < 2) {
401 return std::move(cyls[0]);
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());
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());
421 if (nv == 4 || nv == 6) {
422 std::vector<double> xstep;
423 std::vector<std::pair<double, double>> ystep;
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)));
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) {
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) {
445 std::unique_ptr<VolumeBounds> volBounds{};
448 volBounds = std::make_unique<TrapezoidVolumeBounds>(ystep[0].
second,
453 volBounds.release());
460 volBounds = std::make_unique<DoubleTrapezoidVolumeBounds>(ystep[0].
second,
467 volBounds.release());
472 auto newBounds = std::make_unique<SimplePolygonBrepVolumeBounds>(ivtx, spb->getDZ());
474 newBounds.release());
477 else if (
sh->type() ==
"Subtraction") {
478 const GeoShapeSubtraction* sub =
dynamic_cast<const GeoShapeSubtraction*
>(
sh);
480 const GeoShape* shA = sub->getOpA();
481 const GeoShape* shB = sub->getOpB();
484 auto volBounds = std::make_unique<SubtractedVolumeBounds>(volA.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();
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);
499 const GeoShape* shA =
intersect->getOpA();
500 const GeoShape* shB =
intersect->getOpB();
503 auto volBounds = std::make_unique<CombinedVolumeBounds>(volA.release(),
504 volB.release(),
true);
505 return std::make_unique<Volume>(
nullptr, volBounds.release());
508 if (
sh->type() ==
"Shift") {
509 const GeoShapeShift* shift =
dynamic_cast<const GeoShapeShift*
>(
sh);