ATLAS Offline Software
VolumeConverter.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
3 */
4 
6 
8 // Trk
20 // GeoModel
21 #include "GeoModelKernel/GeoShapeIntersection.h"
22 #include "GeoModelKernel/GeoShapeShift.h"
23 #include "GeoModelKernel/GeoShapeSubtraction.h"
24 #include "GeoModelKernel/GeoShapeUnion.h"
25 #include "GeoModelKernel/GeoTrd.h"
27 
28 // STL
29 #include <algorithm>
30 #include <iostream>
31 namespace {
32  const Trk::Material dummyMaterial{1.e10, 1.e10, 0., 0., 0.};
33 
34 std::unique_ptr<Amg::Transform3D> makeTransform(const Amg::Transform3D& trf) {
35  return std::make_unique<Amg::Transform3D>(trf);
36  }
37 }
38 
39 namespace Trk {
41 
42 std::unique_ptr<TrackingVolume> VolumeConverter::translate(const GeoVPhysVol* gv,
43  bool simplify, bool blend,
44  double blendMassLimit) const {
45 
46  const std::string name = gv->getLogVol()->getName();
47 
48  Amg::Transform3D ident{Amg::Transform3D::Identity()};
49  std::unique_ptr<Volume> volGeo{m_geoShapeConverter.translateGeoShape(
50  gv->getLogVol()->getShape(), ident)};
51 
52  // resolve volume structure into a set of non-overlapping subtractions from
53  // analytically calculable shapes
54  VolumePairVec constituents = splitComposedVolume(*volGeo);
55 
56  // material properties
57  Material mat = Trk::GeoMaterialConverter::convert(gv->getLogVol()->getMaterial());
58 
59  // calculate precision of volume estimate taking into account material
60  // properties
61  double precision = s_precisionInX0 * mat.X0; // required precision in mm
62 
63  // volume estimate from GeoShape
64  double volumeFromGeoShape =
65  -1; // replace with database info when available
66 
67  // volume estimate from resolveBoolean
68  // double volumeBoolean = calculateVolume(volGeo,false,pow(precision,3)); //
69  // TODO : test on inert material
70 
71  // volume estimate from Volume
72  double volume = volumeFromGeoShape >= 0 ? volumeFromGeoShape : -1.;
73  if ((simplify || blend) && volumeFromGeoShape < 0) {
74  double fraction = 0.;
75  volume = 0.;
76  for (const VolumePair& cs : constituents) {
77  fraction = estimateFraction(cs, precision);
78  if (fraction < 0) {
79  volume = -1;
80  break;
81  } else
82  volume += fraction * calculateVolume(*cs.first);
83  }
84  }
85 
86  // evaluate complexity of shape
87  // simple case
88  if (constituents.size() == 1 && !constituents[0].second) {
89  return std::make_unique<TrackingVolume>(*volGeo, mat, nullptr, nullptr,
90  name);
91  }
92  // build envelope
93  std::unique_ptr<Volume> envelope{};
94  std::string envName = name;
95 
96  std::unique_ptr<TrackingVolume> trEnv{};
97 
98  bool blended = false;
99 
100  if (constituents.size() == 1) {
101 
102  envelope = std::make_unique<Volume>(*(constituents.front().first),
103  volGeo->transform());
104  double volEnv = calculateVolume(*constituents.front().first);
105 
106  if (blend && volume > 0 && volEnv > 0 &&
107  volume * mat.rho < blendMassLimit)
108  blended = true;
109 
110  if ((simplify || blended) && volume > 0 && volEnv > 0) {
111  // simplified material rescales X0, l0 and density
112  double fraction = volume / volEnv;
113  Material matScaled(mat.X0 / fraction, mat.L0 / fraction, mat.A,
114  mat.Z, fraction * mat.rho);
115  if (blend && !blended)
116  envName = envName + "_PERM";
117  trEnv = std::make_unique<TrackingVolume>(*envelope, matScaled,
118  nullptr, nullptr, envName);
119  } else {
120  auto confinedVols =std::make_unique<std::vector<TrackingVolume*>>();
121  confinedVols->push_back(std::make_unique<TrackingVolume>(*volGeo, mat, nullptr, nullptr, name).release());
122  envName = name + "_envelope";
123  trEnv = std::make_unique<TrackingVolume>(*envelope, dummyMaterial, std::move(confinedVols), envName);
124  }
125 
126  return trEnv;
127  }
128 
129  // composed shapes : derive envelope from span
130  Amg::Transform3D transf = volGeo->transform();
131  std::unique_ptr<VolumeSpan> span =
132  findVolumeSpan(volGeo->volumeBounds(), transf, 0., 0.);
133 
134  bool isCyl = false;
135  for (const auto& fv : constituents) {
136  const CylinderVolumeBounds* cyl =
137  dynamic_cast<const CylinderVolumeBounds*>(
138  &(fv.first->volumeBounds()));
139  if (cyl) {
140  isCyl = true;
141  break;
142  }
143  }
144 
145  ATH_MSG_DEBUG(__FILE__ << ":" << __LINE__
146  << "envelope estimate: object contains cylinder:"
147  << name << ":" << isCyl);
148  ATH_MSG_DEBUG(__FILE__ << ":" << __LINE__
149  << "complex volume span for envelope:" << name
150  << ":x range:" << (*span).xMin << ","
151  << (*span).xMax);
152  ATH_MSG_DEBUG(__FILE__ << ":" << __LINE__
153  << "complex volume span for envelope:" << name
154  << ":y range:" << (*span).yMin << ","
155  << (*span).yMax);
156  ATH_MSG_DEBUG(__FILE__ << ":" << __LINE__
157  << "complex volume span for envelope:" << name
158  << ":z range:" << (*span).zMin << ","
159  << (*span).zMax);
160  ATH_MSG_DEBUG(__FILE__ << ":" << __LINE__
161  << "complex volume span for envelope:" << name
162  << ":R range:" << (*span).rMin << ","
163  << (*span).rMax);
164  ATH_MSG_DEBUG(__FILE__ << ":" << __LINE__
165  << "complex volume span for envelope:" << name
166  << ":phi range:" << (*span).phiMin << ","
167  << (*span).phiMax);
168 
169  if (!isCyl) { // cuboid envelope
170  Amg::Transform3D cylTrf{
171  transf * Amg::Translation3D{0.5 * ((*span).xMin + (*span).xMax),
172  0.5 * ((*span).yMin + (*span).yMax),
173  0.5 * ((*span).zMin + (*span).zMax)}};
174 
175  std::shared_ptr<VolumeBounds> bounds =
176  std::make_shared<CuboidVolumeBounds>(
177  0.5 * ((*span).xMax - (*span).xMin),
178  0.5 * ((*span).yMax - (*span).yMin),
179  0.5 * ((*span).zMax - (*span).zMin));
180  envelope = std::make_unique<Volume>(
181  makeTransform(cylTrf), std::move(bounds));
182  } else {
183  double dPhi = (*span).phiMin > (*span).phiMax
184  ? (*span).phiMax - (*span).phiMin + 2 * M_PI
185  : (*span).phiMax - (*span).phiMin;
186  std::shared_ptr<VolumeBounds> cylBounds{};
187  Amg::Transform3D cylTrf{transf};
188  if (dPhi < 2 * M_PI) {
189  double aPhi = 0.5 * ((*span).phiMax + (*span).phiMin);
190  cylBounds = std::make_shared<CylinderVolumeBounds>(
191  (*span).rMin, (*span).rMax, 0.5 * dPhi,
192  0.5 * ((*span).zMax - (*span).zMin));
193  cylTrf = cylTrf * Amg::getRotateZ3D(aPhi);
194  } else {
195  cylBounds = std::make_shared<CylinderVolumeBounds>(
196  (*span).rMin, (*span).rMax,
197  0.5 * ((*span).zMax - (*span).zMin));
198  }
199  envelope = std::make_unique<Volume>(
200  makeTransform(cylTrf), std::move(cylBounds));
201  }
202 
203  double volEnv = calculateVolume(*envelope);
204 
205  if (blend && volume > 0 && volEnv > 0 && volume * mat.rho < blendMassLimit)
206  blended = true;
207 
208  if ((simplify || blended) && volume > 0 && volEnv > 0) {
209  double fraction = volume / volEnv;
210  Material matScaled(mat.X0 / fraction, mat.L0 / fraction, mat.A, mat.Z,
211  fraction * mat.rho);
212  if (blend && !blended)
213  envName = envName + "_PERM";
214  trEnv = std::make_unique<TrackingVolume>(*envelope, mat, nullptr,
215  nullptr, envName);
216  } else {
217  auto confinedVols = std::make_unique<std::vector<TrackingVolume*>>();
218  confinedVols->push_back( std::make_unique<TrackingVolume>(*volGeo, mat, nullptr, nullptr, name).release());
219  envName = envName + "_envelope";
220  trEnv = std::make_unique<TrackingVolume>(*envelope, dummyMaterial, std::move(confinedVols), envName);
221  }
222 
223  return trEnv;
224 }
225 
227  double tolerance) const {
228 
229  VolumePartVec constituents{};
230  VolumePart inputVol{};
231  inputVol.parts.push_back(std::make_unique<Volume>(trVol));
232  constituents.push_back(std::move(inputVol));
233  VolumePartVec::iterator sIter = constituents.begin();
234 
236  double volume = 0;
237  while (sIter != constituents.end()) {
238  bool update = false;
239  for (unsigned int ii = 0; ii < (*sIter).parts.size(); ++ii) {
240  const VolumeBounds& bounds{((*sIter).parts[ii]->volumeBounds())};
241  const CombinedVolumeBounds* comb =
242  dynamic_cast<const CombinedVolumeBounds*>(&bounds);
243  const SubtractedVolumeBounds* sub =
244  dynamic_cast<const SubtractedVolumeBounds*>(&bounds);
245  if (comb) {
246  (*sIter).parts[ii].reset(comb->first()->clone());
247  VolumePart vp(*sIter); //copy here
248  constituents.push_back(vp); //inser copy the iter can be invalidated
249  constituents.back().parts[ii].reset(comb->second()->clone()); //modify
250  constituents.push_back(vp); //push copy
251  constituents.back().parts.emplace_back(comb->second()->clone());//modify
252  constituents.back().sign = -1. * constituents.back().sign;
253  update = true;
254  break;
255  } else if (sub) {
256  (*sIter).parts[ii].reset(sub->outer()->clone());
258  double volSub = calculateVolume(*sub->inner(), true, tolerance);
259  if (volSub < tolerance) {
260  volume += -1. * (*sIter).sign * volSub;
261  } else {
262  constituents.emplace_back(*sIter);
263  constituents.back().parts.emplace_back(
264  sub->inner()->clone());
265  constituents.back().sign = -1. * constituents.back().sign;
266  }
267  update = true;
268  break;
269  } else {
270  // component small, below tolerance
271  double volSmall = calculateVolume(*(*sIter).parts[ii]);
272  if (volSmall < tolerance) {
273  sIter=constituents.erase(sIter);
274  update = true;
275  break;
276  }
277  }
278  } //
279  if (update)
280  sIter = constituents.begin();
281  else if ((*sIter).parts.size() == 1) {
282  double volSingle = calculateVolume(*(*sIter).parts[0]);
283  volume += (*sIter).sign * volSingle;
284  sIter=constituents.erase(sIter);
285  } else {
286  std::vector<std::shared_ptr<Volume>>::iterator tit =
287  (*sIter).parts.begin();
288  bool noovrlp = false;
289  while (tit + 1 != (*sIter).parts.end()) {
290  std::pair<bool, std::unique_ptr<Volume>> overlap =
292  if (overlap.first && !overlap.second) {
293  sIter=constituents.erase(sIter);
294  noovrlp = true;
295  break;
296  } // no intersection
297  else if (overlap.first && overlap.second) {
298  (*sIter).parts.erase(tit, tit + 2);
299  (*sIter).parts.push_back(std::move(overlap.second));
300  tit = (*sIter).parts.begin();
301  } else {
302  if (calculateVolume(**tit) < tolerance) {
303  sIter=constituents.erase(sIter);
304  noovrlp = true;
305  break;
306  }
307  if (calculateVolume(**(tit + 1)) < tolerance) {
308  sIter=constituents.erase(sIter);
309  noovrlp = true;
310  break;
311  }
312  std::pair<bool, std::unique_ptr<Volume>> overlap =
314  **tit, **(tit + 1));
315  if (overlap.first) {
316  if (overlap.second) {
317  (*sIter).parts.erase(tit, tit + 2);
318  (*sIter).parts.push_back(std::move(overlap.second));
319  tit = (*sIter).parts.begin();
320  } else {
321  sIter=constituents.erase(sIter);
322  noovrlp = true;
323  break; // no intersection
324  }
325  } else
326  ++tit;
327  }
328  }
329  if (noovrlp) {
330  } else if ((*sIter).parts.size() == 1) {
331  double volSingle = calculateVolume(*(*sIter).parts[0]);
332  volume += (*sIter).sign * volSingle;
333  sIter=constituents.erase(sIter);
334  } else {
335  ++sIter;
336  }
337  }
338  }
339 
340  if (!constituents.empty()) {
341  ATH_MSG_VERBOSE("boolean volume resolved to "
342  << constituents.size() << " items "
343  << ":volume estimate:" << volume);
344  }
345  return volume;
346 }
347 
349  const Trk::Volume& trVol) {
350 
351  VolumePairVec constituents;
352  constituents.emplace_back(std::make_unique<Volume>(trVol), nullptr);
353  VolumePairVec::iterator sIter = constituents.begin();
354  std::shared_ptr<VolumeBounds> newBounds{};
355  while (sIter != constituents.end()) {
357  const CombinedVolumeBounds* comb =
358  dynamic_cast<const Trk::CombinedVolumeBounds*>(
359  &((*sIter).first->volumeBounds()));
360  const Trk::SubtractedVolumeBounds* sub =
361  dynamic_cast<const Trk::SubtractedVolumeBounds*>(
362  &((*sIter).first->volumeBounds()));
364  if (comb) {
365  std::shared_ptr<Volume> subVol = (*sIter).second;
366  sIter = constituents.erase(sIter);
367  std::shared_ptr<Volume> combFirst{comb->first()->clone()};
368  std::shared_ptr<Volume> combSecond{comb->second()->clone()};
369  if (comb->intersection()) {
370  newBounds = std::make_shared<Trk::SubtractedVolumeBounds>(
371  std::unique_ptr<Trk::Volume>(combFirst->clone()), std::unique_ptr<Trk::Volume>(combSecond->clone()));
372  std::unique_ptr<Trk::Volume> newSubVol =
373  std::make_unique<Volume>(nullptr, std::move(newBounds));
374  if (subVol) {
375  newBounds = std::make_shared<CombinedVolumeBounds>(
376  std::unique_ptr<Trk::Volume>(subVol->clone()), std::move(newSubVol), false);
377  std::shared_ptr<Volume> newCSubVol =
378  std::make_unique<Volume>(nullptr, std::move(newBounds));
379  constituents.insert(sIter,
380  std::make_pair(combFirst, newCSubVol));
381  } else {
382  constituents.insert(
383  sIter, std::make_pair(combFirst, std::move(newSubVol)));
384  }
385  } else {
386  constituents.insert(sIter, std::make_pair(combFirst, subVol));
387  if (subVol) {
388  newBounds = std::make_shared<CombinedVolumeBounds>(
389  std::unique_ptr<Trk::Volume>(subVol->clone()),
390  std::unique_ptr<Trk::Volume>(combFirst->clone()), false);
391  std::unique_ptr<Trk::Volume> newSubVol =
392  std::make_unique<Volume>(nullptr, std::move(newBounds));
393  constituents.insert(
394  sIter,
395  std::make_pair(combSecond, std::move(newSubVol)));
396  } else {
397  constituents.insert(sIter,
398  std::make_pair(combSecond, combFirst));
399  }
400  }
401  sIter = constituents.begin();
402  } else if (sub) {
403  std::shared_ptr<Volume> subVol = (*sIter).second;
404  sIter = constituents.erase(sIter);
405  std::shared_ptr<Volume> innerVol{sub->inner()->clone()};
406  std::shared_ptr<Volume> outerVol{sub->outer()->clone()};
407  if (subVol) {
408  newBounds = std::make_shared<CombinedVolumeBounds>(
409  std::unique_ptr<Trk::Volume>(subVol->clone()),
410  std::unique_ptr<Trk::Volume>(innerVol->clone()), false);
411  std::unique_ptr<Volume> newSubVol =
412  std::make_unique<Trk::Volume>(nullptr, newBounds);
413  constituents.insert(
414  sIter, std::make_pair(outerVol, std::move(newSubVol)));
415  } else {
416  constituents.insert(sIter, std::make_pair(outerVol, innerVol));
417  }
418  sIter = constituents.begin();
419  } else {
420  ++sIter;
421  }
422  }
423  return constituents;
424 }
425 
426 std::unique_ptr<VolumeSpan> VolumeConverter::findVolumeSpan(
427  const VolumeBounds& volBounds, const Amg::Transform3D& transform,
428  double zTol, double phiTol) const {
429  // volume shape
430  const CuboidVolumeBounds* box =
431  dynamic_cast<const CuboidVolumeBounds*>(&volBounds);
432  const TrapezoidVolumeBounds* trd =
433  dynamic_cast<const TrapezoidVolumeBounds*>(&volBounds);
434  const DoubleTrapezoidVolumeBounds* dtrd =
435  dynamic_cast<const DoubleTrapezoidVolumeBounds*>(&volBounds);
436  const BevelledCylinderVolumeBounds* bcyl =
437  dynamic_cast<const BevelledCylinderVolumeBounds*>(&volBounds);
438  const CylinderVolumeBounds* cyl =
439  dynamic_cast<const CylinderVolumeBounds*>(&volBounds);
440  const SubtractedVolumeBounds* sub =
441  dynamic_cast<const SubtractedVolumeBounds*>(&volBounds);
442  const CombinedVolumeBounds* comb =
443  dynamic_cast<const CombinedVolumeBounds*>(&volBounds);
444  const SimplePolygonBrepVolumeBounds* spb =
445  dynamic_cast<const SimplePolygonBrepVolumeBounds*>(&volBounds);
446  const PrismVolumeBounds* prism =
447  dynamic_cast<const PrismVolumeBounds*>(&volBounds);
448 
449  double dPhi = 0.;
450 
451  if (sub) {
452  return findVolumeSpan(sub->outer()->volumeBounds(),
453  transform * sub->outer()->transform(), zTol,
454  phiTol);
455  }
456 
457  if (comb) {
458  std::unique_ptr<VolumeSpan> s1 = findVolumeSpan(
459  comb->first()->volumeBounds(),
460  transform * comb->first()->transform(), zTol, phiTol);
461  std::unique_ptr<VolumeSpan> s2 = findVolumeSpan(
462  comb->second()->volumeBounds(),
463  transform * comb->second()->transform(), zTol, phiTol);
464 
465  VolumeSpan scomb;
466  scomb.rMin = std::min((*s1).rMin, (*s2).rMin);
467  scomb.rMax = std::max((*s1).rMax, (*s2).rMax);
468  scomb.xMin = std::min((*s1).xMin, (*s2).xMin);
469  scomb.xMax = std::max((*s1).xMax, (*s2).xMax);
470  scomb.yMin = std::min((*s1).yMin, (*s2).yMin);
471  scomb.yMax = std::max((*s1).yMax, (*s2).yMax);
472  scomb.zMin = std::min((*s1).zMin, (*s2).zMin);
473  scomb.zMax = std::max((*s1).zMax, (*s2).zMax);
474  if ((*s1).phiMin < (*s1).phiMax && (*s2).phiMin < (*s2).phiMax) {
475  scomb.phiMin = std::min((*s1).phiMin, (*s2).phiMin);
476  scomb.phiMax = std::max((*s1).phiMax, (*s2).phiMax);
477  } else if ((*s1).phiMin < (*s1).phiMax && (*s2).phiMin > (*s2).phiMax) {
478  if ((*s1).phiMin > (*s2).phiMax) {
479  scomb.phiMin = std::min((*s1).phiMin, (*s2).phiMin);
480  scomb.phiMax = (*s2).phiMax;
481  } else if ((*s1).phiMax < (*s2).phiMin) {
482  scomb.phiMin = (*s2).phiMin;
483  scomb.phiMax = std::max((*s1).phiMax, (*s2).phiMax);
484  } else {
485  scomb.phiMin = 0.;
486  scomb.phiMax = 2 * M_PI;
487  }
488  } else if ((*s1).phiMin > (*s1).phiMax && (*s2).phiMin < (*s2).phiMax) {
489  if ((*s2).phiMin > (*s1).phiMax) {
490  scomb.phiMin = std::min((*s1).phiMin, (*s2).phiMin);
491  scomb.phiMax = (*s1).phiMax;
492  } else if ((*s2).phiMax < (*s1).phiMin) {
493  scomb.phiMin = (*s1).phiMin;
494  scomb.phiMax = std::max((*s1).phiMax, (*s2).phiMax);
495  } else {
496  scomb.phiMin = 0.;
497  scomb.phiMax = 2 * M_PI;
498  }
499  } else {
500  scomb.phiMin = std::min((*s1).phiMin, (*s2).phiMin);
501  scomb.phiMax = std::max((*s1).phiMax, (*s2).phiMax);
502  }
503  return std::make_unique<VolumeSpan>(scomb);
504  }
505 
506  //
507  double minZ{1.e6};
508  double maxZ{-1.e6};
509  double minPhi{2 * M_PI};
510  double maxPhi{0.};
511  double minR{1.e6};
512  double maxR{0.};
513  double minX{1.e6};
514  double maxX{-1.e6};
515  double minY{1.e6};
516  double maxY{-1.e6};
517 
518  // defined vertices and edges
519  std::vector<Amg::Vector3D> vtx;
520  std::vector<std::pair<int, int>> edges;
522 
523  if (box) {
524  vtx.emplace_back(box->halflengthX(), box->halflengthY(),
525  box->halflengthZ());
526  vtx.emplace_back(-box->halflengthX(), box->halflengthY(),
527  box->halflengthZ());
528  vtx.emplace_back(box->halflengthX(), -box->halflengthY(),
529  box->halflengthZ());
530  vtx.emplace_back(-box->halflengthX(), -box->halflengthY(),
531  box->halflengthZ());
532  vtx.emplace_back(box->halflengthX(), box->halflengthY(),
533  -box->halflengthZ());
534  vtx.emplace_back(-box->halflengthX(), box->halflengthY(),
535  -box->halflengthZ());
536  vtx.emplace_back(box->halflengthX(), -box->halflengthY(),
537  -box->halflengthZ());
538  vtx.emplace_back(-box->halflengthX(), -box->halflengthY(),
539  -box->halflengthZ());
540  edges.emplace_back(0, 1);
541  edges.emplace_back(0, 2);
542  edges.emplace_back(1, 3);
543  edges.emplace_back(2, 3);
544  edges.emplace_back(4, 5);
545  edges.emplace_back(4, 6);
546  edges.emplace_back(5, 7);
547  edges.emplace_back(6, 7);
548  edges.emplace_back(0, 4);
549  edges.emplace_back(1, 5);
550  edges.emplace_back(2, 6);
551  edges.emplace_back(3, 7);
552  }
553  if (trd) {
554  vtx.emplace_back(trd->maxHalflengthX(), trd->halflengthY(),
555  trd->halflengthZ());
556  vtx.emplace_back(-trd->maxHalflengthX(), trd->halflengthY(),
557  trd->halflengthZ());
558  vtx.emplace_back(trd->minHalflengthX(), -trd->halflengthY(),
559  trd->halflengthZ());
560  vtx.emplace_back(-trd->minHalflengthX(), -trd->halflengthY(),
561  trd->halflengthZ());
562  vtx.emplace_back(trd->maxHalflengthX(), trd->halflengthY(),
563  -trd->halflengthZ());
564  vtx.emplace_back(-trd->maxHalflengthX(), trd->halflengthY(),
565  -trd->halflengthZ());
566  vtx.emplace_back(trd->minHalflengthX(), -trd->halflengthY(),
567  -trd->halflengthZ());
568  vtx.emplace_back(-trd->minHalflengthX(), -trd->halflengthY(),
569  -trd->halflengthZ());
570  edges.emplace_back(0, 1);
571  edges.emplace_back(0, 2);
572  edges.emplace_back(1, 3);
573  edges.emplace_back(2, 3);
574  edges.emplace_back(4, 5);
575  edges.emplace_back(4, 6);
576  edges.emplace_back(5, 7);
577  edges.emplace_back(6, 7);
578  edges.emplace_back(0, 4);
579  edges.emplace_back(1, 5);
580  edges.emplace_back(2, 6);
581  edges.emplace_back(3, 7);
582  }
583  if (dtrd) {
584  vtx.emplace_back(dtrd->maxHalflengthX(), 2 * dtrd->halflengthY2(),
585  dtrd->halflengthZ());
586  vtx.emplace_back(-dtrd->maxHalflengthX(), 2 * dtrd->halflengthY2(),
587  dtrd->halflengthZ());
588  vtx.emplace_back(dtrd->medHalflengthX(), 0., dtrd->halflengthZ());
589  vtx.emplace_back(-dtrd->medHalflengthX(), 0., dtrd->halflengthZ());
590  vtx.emplace_back(dtrd->minHalflengthX(), -2 * dtrd->halflengthY1(),
591  dtrd->halflengthZ());
592  vtx.emplace_back(-dtrd->minHalflengthX(), -2 * dtrd->halflengthY1(),
593  dtrd->halflengthZ());
594  vtx.emplace_back(dtrd->maxHalflengthX(), 2 * dtrd->halflengthY2(),
595  -dtrd->halflengthZ());
596  vtx.emplace_back(-dtrd->maxHalflengthX(), 2 * dtrd->halflengthY2(),
597  -dtrd->halflengthZ());
598  vtx.emplace_back(dtrd->medHalflengthX(), 0., -dtrd->halflengthZ());
599  vtx.emplace_back(-dtrd->medHalflengthX(), 0., -dtrd->halflengthZ());
600  vtx.emplace_back(dtrd->minHalflengthX(), -2 * dtrd->halflengthY1(),
601  -dtrd->halflengthZ());
602  vtx.emplace_back(-dtrd->minHalflengthX(), -2 * dtrd->halflengthY1(),
603  -dtrd->halflengthZ());
604  edges.emplace_back(0, 1);
605  edges.emplace_back(0, 2);
606  edges.emplace_back(1, 3);
607  edges.emplace_back(2, 4);
608  edges.emplace_back(3, 5);
609  edges.emplace_back(4, 5);
610  edges.emplace_back(6, 7);
611  edges.emplace_back(6, 8);
612  edges.emplace_back(7, 9);
613  edges.emplace_back(8, 10);
614  edges.emplace_back(9, 11);
615  edges.emplace_back(10, 11);
616  edges.emplace_back(0, 6);
617  edges.emplace_back(1, 7);
618  edges.emplace_back(2, 8);
619  edges.emplace_back(3, 9);
620  edges.emplace_back(4, 10);
621  edges.emplace_back(5, 11);
622  }
623  if (bcyl) {
624  dPhi = bcyl->halfPhiSector();
625  vtx.emplace_back(0., 0., bcyl->halflengthZ());
626  vtx.emplace_back(0., 0., -bcyl->halflengthZ());
627  edges.emplace_back(0, 1);
628  if (dPhi < M_PI) {
629  const double cosDphi = std::cos(dPhi);
630  const double sinDphi = std::sin(dPhi);
631  vtx.emplace_back(bcyl->outerRadius() * cosDphi,
632  bcyl->outerRadius() * sinDphi,
633  bcyl->halflengthZ());
634  vtx.emplace_back(bcyl->innerRadius() * cosDphi,
635  bcyl->innerRadius() * sinDphi,
636  bcyl->halflengthZ());
637  vtx.emplace_back(bcyl->outerRadius() * cosDphi,
638  -bcyl->outerRadius() * sinDphi,
639  bcyl->halflengthZ());
640  vtx.emplace_back(bcyl->innerRadius() * cosDphi,
641  -bcyl->innerRadius() * sinDphi,
642  bcyl->halflengthZ());
643  vtx.emplace_back(bcyl->outerRadius() * cosDphi,
644  bcyl->outerRadius() * sinDphi,
645  -bcyl->halflengthZ());
646  vtx.emplace_back(bcyl->innerRadius() * cosDphi,
647  bcyl->innerRadius() * sinDphi,
648  -bcyl->halflengthZ());
649  vtx.emplace_back(bcyl->outerRadius() * cosDphi,
650  -bcyl->outerRadius() * sinDphi,
651  -bcyl->halflengthZ());
652  vtx.emplace_back(bcyl->innerRadius() * cosDphi,
653  -bcyl->innerRadius() * sinDphi,
654  -bcyl->halflengthZ());
655  vtx.emplace_back(bcyl->outerRadius(), 0.,
656  0.); // to distinguish phi intervals for cylinders
657  // aligned with z axis
658  edges.emplace_back(2, 3);
659  edges.emplace_back(4, 5);
660  edges.emplace_back(6, 7);
661  edges.emplace_back(8, 9);
662  if (bcyl->type() == 1 || bcyl->type() == 3) {
663  edges.emplace_back(3, 5);
664  edges.emplace_back(7, 9);
665  }
666  if (bcyl->type() == 2 || bcyl->type() == 3) {
667  edges.emplace_back(2, 4);
668  edges.emplace_back(6, 8);
669  }
670  }
671  }
672  if (cyl) {
673  dPhi = cyl->halfPhiSector();
674  vtx.emplace_back(0., 0., cyl->halflengthZ());
675  vtx.emplace_back(0., 0., -cyl->halflengthZ());
676  edges.emplace_back(0, 1);
677  if (dPhi < M_PI) {
678  const double cosDphi = std::cos(dPhi);
679  const double sinDphi = std::sin(dPhi);
680  vtx.emplace_back(cyl->outerRadius() * cosDphi,
681  cyl->outerRadius() * sinDphi, cyl->halflengthZ());
682  vtx.emplace_back(cyl->innerRadius() * cosDphi,
683  cyl->innerRadius() * sinDphi, cyl->halflengthZ());
684  vtx.emplace_back(cyl->outerRadius() * cosDphi,
685  -cyl->outerRadius() * sinDphi, cyl->halflengthZ());
686  vtx.emplace_back(cyl->outerRadius() * cosDphi,
687  -cyl->outerRadius() * sinDphi, cyl->halflengthZ());
688  vtx.emplace_back(cyl->outerRadius() * cosDphi,
689  cyl->outerRadius() * sinDphi, -cyl->halflengthZ());
690  vtx.emplace_back(cyl->innerRadius() * cosDphi,
691  cyl->innerRadius() * sinDphi, -cyl->halflengthZ());
692  vtx.emplace_back(cyl->outerRadius() * cosDphi,
693  -cyl->outerRadius() * sinDphi,
694  -cyl->halflengthZ());
695  vtx.emplace_back(cyl->outerRadius() * cosDphi,
696  -cyl->outerRadius() * sinDphi,
697  -cyl->halflengthZ());
698  vtx.emplace_back(cyl->outerRadius(), 0.,
699  0.); // to distinguish phi intervals for cylinders
700  // aligned with z axis
701  edges.emplace_back(2, 3);
702  edges.emplace_back(4, 5);
703  edges.emplace_back(6, 7);
704  edges.emplace_back(8, 9);
705  }
706  }
707 
708  if (spb) {
709  const std::vector<std::pair<double, double>> vtcs = spb->xyVertices();
710  for (const auto& vtc : vtcs) {
711  vtx.emplace_back(vtc.first, vtc.second, spb->halflengthZ());
712  vtx.emplace_back(vtc.first, vtc.second, -spb->halflengthZ());
713  edges.emplace_back(vtx.size() - 2, vtx.size() - 1);
714  if (vtx.size() > 2) {
715  edges.emplace_back(
716  vtx.size() - 4, vtx.size() - 2);
717  edges.emplace_back(
718  vtx.size() - 3, vtx.size() - 1);
719  }
720  if (vtx.size() > 4) { // some diagonals
721  edges.emplace_back(vtx.size() - 2, 1);
722  edges.emplace_back(vtx.size() - 1, 0);
723  }
724  }
725  edges.emplace_back(0, vtx.size() - 2);
726  edges.emplace_back(1, vtx.size() - 1);
727  }
728 
729  if (prism) {
730  const std::vector<std::pair<double, double>> vtcs = prism->xyVertices();
731  for (const auto& vtc : vtcs) {
732  vtx.emplace_back(vtc.first, vtc.second, prism->halflengthZ());
733  vtx.emplace_back(vtc.first, vtc.second, -prism->halflengthZ());
734  edges.emplace_back(vtx.size() - 2, vtx.size() - 1);
735  if (vtx.size() > 2) {
736  edges.emplace_back(
737  vtx.size() - 4, vtx.size() - 2);
738  edges.emplace_back(
739  vtx.size() - 3, vtx.size() - 1);
740  }
741  }
742  edges.emplace_back(0, vtx.size() - 2);
743  edges.emplace_back(1, vtx.size() - 1);
744  }
745 
746  std::vector<Amg::Vector3D> vtxt;
747 
748  for (unsigned int ie = 0; ie < vtx.size(); ie++) {
749  Amg::Vector3D gp = transform * vtx[ie];
750  vtxt.push_back(gp);
751 
752  double phi = gp.phi() + M_PI;
753  double rad = gp.perp();
754 
755  // collect limits from vertices
756  minX = std::min(minX, gp[0]);
757  maxX = std::max(maxX, gp[0]);
758  minY = std::min(minY, gp[1]);
759  maxY = std::max(maxY, gp[1]);
760  minZ = std::min(minZ, gp[2]);
761  maxZ = std::max(maxZ, gp[2]);
762  minR = std::min(minR, rad);
763  maxR = std::max(maxR, rad);
764  maxPhi = std::max(maxPhi, phi);
765  minPhi = std::min(minPhi, phi);
766  }
767 
768  if (cyl || bcyl) {
769 
770  double ro = cyl ? cyl->outerRadius() : bcyl->outerRadius();
771  double ri = cyl ? cyl->innerRadius() : bcyl->innerRadius();
772  // z span corrected for theta inclination
774  (vtxt[edges[0].first] - vtxt[edges[0].second]).unit();
775  maxZ += ro * sin(dir.theta());
776  minZ += -ro * sin(dir.theta());
777  // azimuthal & radial extent
778  if (ro < minR) { // excentric object, phi span driven by z-R extent
779  // calculate point of closest approach
780  PerigeeSurface peri;
781  Intersection closest = peri.straightLineIntersection(vtxt[1], dir);
782  double le = (vtxt[0] - vtxt[1]).norm();
783  if ((closest.position - vtxt[0]).norm() < le &&
784  (closest.position - vtxt[1]).norm() < le) {
785  if (minR > closest.position.perp() - ro)
786  minR = std::max(0., closest.position.perp() - ro);
787  // use for phi check
788  double phiClosest = closest.position.phi() + M_PI;
789  if (phiClosest < minPhi || phiClosest > maxPhi) {
790  double phiTmp = minPhi;
791  minPhi = maxPhi;
792  maxPhi = phiTmp;
793  }
794  } else
795  minR = std::max(0., minR - ro * std::abs(dir.z()));
796 
797  const double aTan = std::atan2(ro, minR);
798  minPhi += -aTan;
799  maxPhi += aTan;
800  if (minPhi < 0)
801  minPhi += 2 * M_PI;
802  if (maxPhi > 2 * M_PI)
803  maxPhi += -2 * M_PI;
804 
805  maxR += ro * std::abs(cos(dir.theta()));
806  } else {
807 
808  double rAx = std::max(vtxt[0].perp(), vtxt[1].perp());
809  if (rAx < ri)
810  minR = ri - rAx;
811  else
812  minR = std::max(0., minR - ro * std::abs(cos(dir.theta())));
813 
814  // loop over edges to check inner radial extent
815  PerigeeSurface peri;
816  for (unsigned int ie = 0; ie < edges.size(); ie++) {
818  (vtxt[edges[ie].first] - vtxt[edges[ie].second]).unit();
819  Intersection closest =
820  peri.straightLineIntersection(vtxt[edges[ie].second], dir);
821  double le =
822  (vtxt[edges[ie].first] - vtxt[edges[ie].second]).norm();
823  if ((closest.position - vtxt[edges[ie].first]).norm() < le &&
824  (closest.position - vtxt[edges[ie].second]).norm() < le)
825  if (minR > closest.position.perp())
826  minR = closest.position.perp();
827  }
828 
829  if (vtxt.size() > 10) { // cylindrical section
830  // find spread of phi extent at section (-) boundary : vertices
831  // 4,5,8,9
832  double phiSecLmin = std::min(
833  std::min(vtxt[4].phi() + M_PI, vtxt[5].phi() + M_PI),
834  std::min(vtxt[8].phi() + M_PI, vtxt[9].phi() + M_PI));
835  double phiSecLmax = std::max(
836  std::max(vtxt[4].phi() + M_PI, vtxt[5].phi() + M_PI),
837  std::max(vtxt[8].phi() + M_PI, vtxt[9].phi() + M_PI));
838  // find spread of phi extent at section (+) boundary : vertices
839  // 2,3,6,7
840  double phiSecUmin = std::min(
841  std::min(vtxt[2].phi() + M_PI, vtxt[3].phi() + M_PI),
842  std::min(vtxt[6].phi() + M_PI, vtxt[7].phi() + M_PI));
843  double phiSecUmax = std::max(
844  std::max(vtxt[2].phi() + M_PI, vtxt[3].phi() + M_PI),
845  std::max(vtxt[6].phi() + M_PI, vtxt[7].phi() + M_PI));
846  minPhi = std::min(std::min(phiSecLmin, phiSecLmax),
847  std::min(phiSecUmin, phiSecUmax));
848  maxPhi = std::max(std::max(phiSecLmin, phiSecLmax),
849  std::max(phiSecUmin, phiSecUmax));
850  if (vtxt[10].phi() + M_PI < minPhi ||
851  vtxt[10].phi() + M_PI > maxPhi) {
852  minPhi = 3 * M_PI;
853  maxPhi = 0.;
854  double phiTmp;
855  for (unsigned int iv = 2; iv < vtxt.size(); iv++) {
856  phiTmp = vtxt[iv].phi() + M_PI;
857  if (phiTmp < M_PI)
858  phiTmp += 2 * M_PI;
859  minPhi = phiTmp < minPhi ? phiTmp : minPhi;
860  maxPhi = phiTmp > maxPhi ? phiTmp : maxPhi;
861  }
862  if (minPhi > 2 * M_PI)
863  minPhi += -2 * M_PI;
864  if (maxPhi > 2 * M_PI)
865  maxPhi += -2 * M_PI;
866  }
867  } else {
868  minPhi = 0.;
869  maxPhi = 2 * M_PI;
870  maxR += ro * std::abs(std::cos(dir.theta()));
871  }
872  if (minPhi >= maxPhi && (minPhi - maxPhi) < M_PI) {
873  minPhi = 0.;
874  maxPhi = 2 * M_PI;
875  }
876  }
877  } // end cyl & bcyl
878 
879  if (!cyl && !bcyl) {
880  // loop over edges to check inner radial extent
881  PerigeeSurface peri;
882  for (unsigned int ie = 0; ie < edges.size(); ie++) {
884  (vtxt[edges[ie].first] - vtxt[edges[ie].second]).unit();
885  Intersection closest =
886  peri.straightLineIntersection(vtxt[edges[ie].second], dir);
887  double le = (vtxt[edges[ie].first] - vtxt[edges[ie].second]).norm();
888  if ((closest.position - vtxt[edges[ie].first]).norm() < le &&
889  (closest.position - vtxt[edges[ie].second]).norm() < le)
890  if (minR > closest.position.perp())
891  minR = closest.position.perp();
892  } // end loop over edges
893  // verify phi span - may run across step
894  if (std::abs(maxPhi - minPhi) > M_PI) {
895  double phiTmp = minPhi;
896  minPhi = 3 * M_PI;
897  maxPhi = 0.; // redo the search
898  for (unsigned int iv = 0; iv < vtxt.size(); iv++) {
899  phiTmp = vtxt[iv].phi() + M_PI;
900  if (phiTmp < M_PI)
901  phiTmp += 2 * M_PI;
902  minPhi = phiTmp < minPhi ? phiTmp : minPhi;
903  maxPhi = phiTmp > maxPhi ? phiTmp : maxPhi;
904  }
905  if (minPhi > 2 * M_PI)
906  minPhi += -2 * M_PI;
907  if (maxPhi > 2 * M_PI)
908  maxPhi += -2 * M_PI;
909  if (minPhi >= maxPhi && (minPhi - maxPhi) < M_PI) {
910  minPhi = 0.;
911  maxPhi = 2 * M_PI;
912  }
913  }
914  }
915 
916  if (cyl || bcyl || box || trd || dtrd || spb || prism) {
917  span.zMin = minZ - zTol;
918  span.zMax = maxZ - +zTol;
919  minPhi = (minPhi - phiTol) < 0 ? minPhi - phiTol + 2 * M_PI
920  : minPhi - phiTol;
921  span.phiMin = minPhi;
922  maxPhi = (maxPhi + phiTol) > 2 * M_PI ? maxPhi + phiTol - 2 * M_PI
923  : maxPhi + phiTol;
924  span.phiMax = maxPhi;
925  span.rMin = std::max(70.001, minR - zTol);
926  span.rMax = maxR + zTol;
927  span.xMin = minX - zTol;
928  span.xMax = maxX - +zTol;
929  span.yMin = minY - zTol;
930  span.yMax = maxY - +zTol;
931  } else {
932  ATH_MSG_WARNING("VolumeConverter::volume shape not recognized ");
933  }
934  return std::make_unique<VolumeSpan>(span);
935 }
936 
937 double VolumeConverter::calculateVolume(const Volume& vol, bool nonBooleanOnly,
938  double precision) const {
939 
940  double volume = -1.;
941 
942  const CylinderVolumeBounds* cyl = dynamic_cast<const CylinderVolumeBounds*>(&(vol.volumeBounds()));
943  const CuboidVolumeBounds* box = dynamic_cast<const CuboidVolumeBounds*>(&(vol.volumeBounds()));
944  const TrapezoidVolumeBounds* trd = dynamic_cast<const TrapezoidVolumeBounds*>(&(vol.volumeBounds()));
945  const BevelledCylinderVolumeBounds* bcyl =dynamic_cast<const BevelledCylinderVolumeBounds*>(&(vol.volumeBounds()));
946  const PrismVolumeBounds* prism =dynamic_cast<const PrismVolumeBounds*>(&(vol.volumeBounds()));
947  const SimplePolygonBrepVolumeBounds* spb = dynamic_cast<const SimplePolygonBrepVolumeBounds*>(&(vol.volumeBounds()));
948  const CombinedVolumeBounds* comb =dynamic_cast<const CombinedVolumeBounds*>(&(vol.volumeBounds()));
949  const SubtractedVolumeBounds* sub = dynamic_cast<const SubtractedVolumeBounds*>(&(vol.volumeBounds()));
950 
951  if (cyl) {
952  return 2 * cyl->halfPhiSector() * cyl->halflengthZ() *
953  (std::pow(cyl->outerRadius(), 2) -
954  std::pow(cyl->innerRadius(), 2));
955  }
956  if (box) {
957  return 8 * box->halflengthX() * box->halflengthY() * box->halflengthZ();
958  }
959  if (trd) {
960  return 4 * (trd->minHalflengthX() + trd->maxHalflengthX()) *
961  trd->halflengthY() * trd->halflengthZ();
962  }
963  if (bcyl) {
964  int type = bcyl->type();
965  if (type < 1)
966  return 2 * bcyl->halfPhiSector() *
967  (std::pow(bcyl->outerRadius(), 2) -
968  std::pow(bcyl->innerRadius(), 2)) *
969  bcyl->halflengthZ();
970  if (type == 1)
971  return 2 * bcyl->halflengthZ() *
972  (bcyl->halfPhiSector() * std::pow(bcyl->outerRadius(), 2) -
973  std::pow(bcyl->innerRadius(), 2) *
974  std::tan(bcyl->halfPhiSector()));
975  if (type == 2)
976  return 2 * bcyl->halflengthZ() *
977  (-bcyl->halfPhiSector() * std::pow(bcyl->innerRadius(), 2) +
978  std::pow(bcyl->outerRadius(), 2) *
979  std::tan(bcyl->halfPhiSector()));
980  if (type == 3)
981  return 2 * bcyl->halflengthZ() * std::tan(bcyl->halfPhiSector()) *
982  (std::pow(bcyl->outerRadius(), 2) -
983  std::pow(bcyl->innerRadius(), 2));
984  }
985  if (prism) {
986 
987  std::vector<std::pair<double, double>> v = prism->xyVertices();
988  double vv = v[0].first * (v[1].second - v.back().second);
989  for (unsigned int i = 1; i < v.size() - 1; i++) {
990  vv += v[i].first * (v[i + 1].second - v[i - 1].second);
991  }
992  vv += v.back().first * (v[0].second - v[v.size() - 2].second);
993  return vv * prism->halflengthZ();
994  }
995  if (spb) {
996  std::vector<std::pair<double, double>> v = spb->xyVertices();
997  double vv = v[0].first * (v[1].second - v.back().second);
998  for (unsigned int i = 1; i < v.size() - 1; i++) {
999  vv += v[i].first * (v[i + 1].second - v[i - 1].second);
1000  }
1001  vv += v.back().first * (v[0].second - v[v.size() - 2].second);
1002  return vv * spb->halflengthZ();
1003  }
1004 
1005  if (nonBooleanOnly)
1006  return volume;
1007 
1008  if (comb || sub) {
1009  return resolveBooleanVolume(vol, precision);
1010  }
1011  return volume;
1012 }
1013 
1015  double precision) const {
1016 
1017  if (!sub.first)
1018  return 0.;
1019 
1020  if (sub.first && !sub.second)
1021  return 1.;
1022  double fraction = -1.;
1023 
1024  std::pair<bool, std::unique_ptr<Volume>> overlap =
1025  Trk::VolumeIntersection::intersect(*sub.first, *sub.second);
1026 
1027  if (overlap.first && !overlap.second)
1028  return fraction = 1.;
1029  else if (overlap.first && overlap.second) {
1030  fraction = 1. - calculateVolume(*overlap.second, true, precision) /
1031  calculateVolume(*sub.first, true, precision);
1032  return fraction;
1033  }
1034  // resolve embedded volumes
1035 
1036  // trivial within required precision
1037  double volA = calculateVolume(*sub.first, true, precision);
1038  double volB = calculateVolume(*sub.second, true, precision);
1039  if ((volA > 0 && volA < precision) || (volB > 0 && volB < precision))
1040  return 1.;
1041 
1042  return fraction;
1043 }
1044 
1045 void VolumeConverter::collectMaterial(const GeoVPhysVol* pv,
1046  MaterialProperties& layMat,
1047  double sf) const {
1048  // sf is the area of the layer collecting the material
1049 
1050  // solution relying on GeoModel
1051  // currently involves hit&miss on-fly calculation of boolean volumes
1052  // GeoModelTools::MaterialComponent mat =
1053  // gm_materialHelper.collectMaterial(pv); Material newMP =
1054  // convert(mat.first); double d = mat.second / sf; layMat.addMaterial(newMP,
1055  // d / newMP.x0()); return;
1056 
1057  std::vector<MaterialComponent> materialContent;
1058  collectMaterialContent(pv, materialContent);
1059 
1060  for (const auto& mat : materialContent) {
1061  if (mat.second < 0)
1062  continue; // protection unsolved booleans
1063  double d = sf > 0 ? mat.second / sf : 0.;
1064  if (d > 0)
1065  layMat.addMaterial(mat.first,
1066  (mat.first.X0 > 0 ? d / mat.first.X0 : 0.));
1067  }
1068 }
1069 
1071  const GeoVPhysVol* gv,
1072  std::vector<MaterialComponent>& materialContent) const {
1073 
1074  // solution relying on GeoModel
1075  // currently involves hit&miss on-fly calculation of boolean volumes
1076  // GeoModelTools::MaterialComponent mat =
1077  // gm_materialHelper.collectMaterial(pv); Material newMP =
1078  // convert(mat.first); materialContent.push_back( MaterialComponent( newMP,
1079  // mat.second) ); return;
1080 
1081  const GeoLogVol* lv = gv->getLogVol();
1082  Material mat = Trk::GeoMaterialConverter::convert(lv->getMaterial());
1083 
1084  double motherVolume = 0.;
1085 
1086  // skip volume calculation for dummy material configuration
1087  if (!Trk::GeoMaterialConverter::dummy_material(lv->getMaterial())) {
1088  const GeoShape* sh = lv->getShape();
1089  while (sh && sh->type() == "Shift") {
1090  const GeoShapeShift* shift = dynamic_cast<const GeoShapeShift*>(sh);
1091  sh = shift ? shift->getOp() : nullptr;
1092  }
1093 
1094  bool isBoolean =
1095  sh && (sh->type() == "Subtraction" || sh->type() == "Union" ||
1096  sh->type() == "Intersection");
1097 
1098  if (isBoolean) {
1099  Amg::Transform3D transf{Amg::Transform3D::Identity()};
1100  std::unique_ptr<Volume> vol{
1102  motherVolume =
1103  calculateVolume(*vol, false, std::pow(1.e-3 * mat.X0, 3));
1104  if (motherVolume < 0) {
1105  // m_geoShapeConverter.decodeShape(sh);
1106  }
1107  } else
1108  motherVolume = lv->getShape()->volume();
1109  }
1110 
1111  double childVol = 0;
1112  std::string cPrevious = " ";
1113  size_t nIdentical = 0;
1114  std::vector<Trk::MaterialComponent> childMat;
1115  std::vector<const GeoVPhysVol*> children = geoGetVolumesNoXform (gv);
1116  for (const GeoVPhysVol* cv : children) {
1117  std::string cname = cv->getLogVol()->getName();
1118  if (cname == cPrevious)
1119  nIdentical++; // assuming identity for identical name and branching
1120  // history
1121  else { // scale and collect material from previous item
1122  for (const auto& cmat : childMat) {
1123  materialContent.emplace_back(cmat.first, nIdentical * cmat.second);
1124  childVol += materialContent.back().second;
1125  }
1126  childMat.clear(); // reset
1127  nIdentical = 1; // current
1128  collectMaterialContent(cv, childMat);
1129  }
1130  }
1131  for (const auto& cmat : childMat) {
1132  materialContent.emplace_back(cmat.first, nIdentical * cmat.second);
1133  childVol += materialContent.back().second;
1134  }
1135  if (motherVolume > 0 && childVol > 0)
1136  motherVolume += -1. * childVol;
1137 
1138  ATH_MSG_DEBUG("collected material:" << lv->getName() << ":made of:"
1139  << lv->getMaterial()->getName()
1140  << ":density(g/mm3)" << mat.rho
1141  << ":mass:" << mat.rho * motherVolume);
1142  materialContent.emplace_back(mat, motherVolume);
1143 }
1144 
1145 double VolumeConverter::leadingVolume(const GeoShape* sh) const {
1146 
1147  if (sh->type() == "Subtraction") {
1148  const GeoShapeSubtraction* sub =
1149  dynamic_cast<const GeoShapeSubtraction*>(sh);
1150  if (sub)
1151  return leadingVolume(sub->getOpA());
1152  }
1153  if (sh->type() == "Union") {
1154  const GeoShapeUnion* uni = dynamic_cast<const GeoShapeUnion*>(sh);
1155  if (uni)
1156  return leadingVolume(uni->getOpA()) + leadingVolume(uni->getOpB());
1157  }
1158  if (sh->type() == "Intersection") {
1159  const GeoShapeIntersection* intr =
1160  dynamic_cast<const GeoShapeIntersection*>(sh);
1161  if (intr)
1162  return std::min(leadingVolume(intr->getOpA()),
1163  leadingVolume(intr->getOpB()));
1164  }
1165  if (sh->type() == "Shift") {
1166  const GeoShapeShift* shift = dynamic_cast<const GeoShapeShift*>(sh);
1167  if (shift)
1168  return leadingVolume(shift->getOp());
1169  }
1170 
1171  return sh->volume();
1172 }
1173 } // namespace Trk
xAOD::iterator
JetConstituentVector::iterator iterator
Definition: JetConstituentVector.cxx:68
Trk::CuboidVolumeBounds::halflengthZ
double halflengthZ() const
This method returns the halflength in local z.
Definition: CuboidVolumeBounds.h:136
AllowedVariables::e
e
Definition: AsgElectronSelectorTool.cxx:37
Trk::TrapezoidVolumeBounds::maxHalflengthX
double maxHalflengthX() const
This method returns the maximal halflength in local x.
Definition: TrapezoidVolumeBounds.h:160
Trk::BevelledCylinderVolumeBounds::outerRadius
double outerRadius() const
This method returns the outer radius.
Definition: BevelledCylinderVolumeBounds.h:238
CxxUtils::span
span(T *ptr, std::size_t sz) -> span< T >
A couple needed deduction guides.
Trk::CombinedVolumeBounds::intersection
bool intersection() const
This method distinguishes between Union(0) and Intersection(1)
Definition: CombinedVolumeBounds.h:127
PlotCalibFromCool.norm
norm
Definition: PlotCalibFromCool.py:100
Trk::SimplePolygonBrepVolumeBounds::halflengthZ
double halflengthZ() const
This method returns the halflength in local z.
Definition: SimplePolygonBrepVolumeBounds.h:151
Trk::TrapezoidVolumeBounds::minHalflengthX
double minHalflengthX() const
This method returns the minimal halflength in local x.
Definition: TrapezoidVolumeBounds.h:157
Trk::VolumePartVec
std::vector< VolumePart > VolumePartVec
Definition: VolumeConverter.h:50
Trk::VolumeConverter::m_geoShapeConverter
Trk::GeoShapeConverter m_geoShapeConverter
shape converter
Definition: VolumeConverter.h:102
Trk::Intersection
Definition: Intersection.h:24
Muon::makeTransform
std::unique_ptr< Amg::Transform3D > makeTransform(const Amg::Transform3D &trf)
Definition: MuonSpectrometer/MuonDetDescr/MuonTrackingGeometry/MuonTrackingGeometry/Utils.h:14
Trk::VolumeSpan::xMax
double xMax
Definition: VolumeConverter.h:39
Trk::BevelledCylinderVolumeBounds
Definition: BevelledCylinderVolumeBounds.h:100
Trk::BevelledCylinderVolumeBounds::type
int type() const
This method returns the type.
Definition: BevelledCylinderVolumeBounds.h:245
Trk::SimplePolygonBrepVolumeBounds
Definition: SimplePolygonBrepVolumeBounds.h:44
Trk::VolumeConverter::calculateVolume
double calculateVolume(const Volume &vol, bool nonBooleanOnly=false, double precision=1.e-3) const
Volume calculation : by default return analytical solution only.
Definition: VolumeConverter.cxx:937
perp
Scalar perp() const
perp method - perpenticular length
Definition: AmgMatrixBasePlugin.h:44
Trk::VolumeConverter::VolumePairVec
std::vector< VolumePair > VolumePairVec
Definition: VolumeConverter.h:71
Trk::PerigeeSurface
Definition: PerigeeSurface.h:43
Trk::CuboidVolumeBounds
Definition: CuboidVolumeBounds.h:52
hist_file_dump.d
d
Definition: hist_file_dump.py:142
max
constexpr double max()
Definition: ap_fixedTest.cxx:33
Trk::VolumeConverter::s_precisionInX0
static constexpr double s_precisionInX0
Definition: VolumeConverter.h:104
mat
GeoMaterial * mat
Definition: LArDetectorConstructionTBEC.cxx:55
min
constexpr double min()
Definition: ap_fixedTest.cxx:26
geoGetVolumesNoXform
std::vector< const GeoVPhysVol * > geoGetVolumesNoXform(const GeoGraphNode *node, int depthLimit=1, int sizeHint=20)
Return the child volumes.
Definition: GeoVisitVolumes.cxx:240
Trk::VolumeConverter::findVolumeSpan
std::unique_ptr< VolumeSpan > findVolumeSpan(const VolumeBounds &volBounds, const Amg::Transform3D &transform, double zTol, double phiTol) const
Estimation of the geometrical volume span.
Definition: VolumeConverter.cxx:426
python.SystemOfUnits.second
float second
Definition: SystemOfUnits.py:135
Trk::GeoMaterialConverter::convert
static Material convert(const GeoMaterial *gm)
Single conversion , input type GeoMaterial - output type Trk::MaterialProperties.
Definition: GeoMaterialConverter.cxx:18
M_PI
#define M_PI
Definition: ActiveFraction.h:11
Trk::VolumeSpan::phiMax
double phiMax
Definition: VolumeConverter.h:35
Trk::CuboidVolumeBounds::halflengthX
double halflengthX() const
This method returns the halflength in local x.
Definition: CuboidVolumeBounds.h:132
drawFromPickle.cos
cos
Definition: drawFromPickle.py:36
Trk::VolumeConverter::collectMaterial
void collectMaterial(const GeoVPhysVol *pv, Trk::MaterialProperties &layMat, double sf) const
material collection for layers
Definition: VolumeConverter.cxx:1045
ATH_MSG_VERBOSE
#define ATH_MSG_VERBOSE(x)
Definition: AthMsgStreamMacros.h:28
Trk::VolumeSpan::rMax
double rMax
Definition: VolumeConverter.h:37
Trk::CombinedVolumeBounds::first
const Volume * first() const
This method returns the first VolumeBounds.
Definition: CombinedVolumeBounds.h:119
Trk::VolumeIntersection::intersect
static std::pair< bool, std::unique_ptr< Trk::Volume > > intersect(const Volume &volA, const Volume &volB)
Definition: VolumeIntersection.cxx:38
python.CaloAddPedShiftConfig.type
type
Definition: CaloAddPedShiftConfig.py:42
Trk::VolumeSpan::yMin
double yMin
Definition: VolumeConverter.h:40
TrapezoidVolumeBounds.h
Trk::DoubleTrapezoidVolumeBounds::halflengthZ
double halflengthZ() const
This method returns the halflength in local z.
Definition: DoubleTrapezoidVolumeBounds.h:200
Trk::TrapezoidVolumeBounds::halflengthY
double halflengthY() const
This method returns the halflength in local y.
Definition: TrapezoidVolumeBounds.h:163
Amg::getRotateZ3D
Amg::Transform3D getRotateZ3D(double angle)
get a rotation transformation around Z-axis
Definition: GeoPrimitivesHelpers.h:270
PlotCalibFromCool.ie
ie
Definition: PlotCalibFromCool.py:420
Trk::DoubleTrapezoidVolumeBounds::maxHalflengthX
double maxHalflengthX() const
This method returns the X halflength at maximal Y (local coordinates)
Definition: DoubleTrapezoidVolumeBounds.h:188
PrismVolumeBounds.h
Trk::CuboidVolumeBounds::halflengthY
double halflengthY() const
This method returns the halflength in local y.
Definition: CuboidVolumeBounds.h:134
GeoVisitVolumes.h
Visitor to process all volumes under a GeoModel node.
Trk::VolumeBounds
Definition: VolumeBounds.h:45
CylinderVolumeBounds.h
TileDCSDataPlotter.tit
tit
Definition: TileDCSDataPlotter.py:892
Trk::MaterialProperties::addMaterial
void addMaterial(const Material &mp, float dInX0)
Material averaging.
Definition: MaterialProperties.cxx:46
Trk::PerigeeSurface::straightLineIntersection
virtual Intersection straightLineIntersection(const Amg::Vector3D &pos, const Amg::Vector3D &dir, bool forceDir=false, Trk::BoundaryCheck bchk=false) const override final
fast straight line intersection schema - standard: provides closest intersection and (signed) path le...
Definition: PerigeeSurface.cxx:226
Trk::VolumeConverter::estimateFraction
double estimateFraction(const VolumePair &sub, double precision) const
the tricky part of volume calculation
Definition: VolumeConverter.cxx:1014
FullCPAlgorithmsTest_eljob.sh
sh
Definition: FullCPAlgorithmsTest_eljob.py:114
lumiFormat.i
int i
Definition: lumiFormat.py:85
Trk::VolumePart
Definition: VolumeConverter.h:46
Trk::DoubleTrapezoidVolumeBounds::halflengthY1
double halflengthY1() const
This method returns the halflength1 in local y.
Definition: DoubleTrapezoidVolumeBounds.h:192
Trk::VolumeConverter::collectMaterialContent
void collectMaterialContent(const GeoVPhysVol *gv, std::vector< Trk::MaterialComponent > &materialContent) const
material collection for volumes
Definition: VolumeConverter.cxx:1070
RCU::Shell
Definition: ShellExec.cxx:28
Trk::VolumeConverter::leadingVolume
double leadingVolume(const GeoShape *sh) const
Definition: VolumeConverter.cxx:1145
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
Trk::PrismVolumeBounds::halflengthZ
double halflengthZ() const
This method returns the halflength in local z.
Definition: PrismVolumeBounds.h:125
TauGNNUtils::Variables::Track::dPhi
bool dPhi(const xAOD::TauJet &tau, const xAOD::TauTrack &track, double &out)
Definition: TauGNNUtils.cxx:549
BevelledCylinderVolumeBounds.h
Amg::Transform3D
Eigen::Affine3d Transform3D
Definition: GeoPrimitives.h:46
Amg::transform
Amg::Vector3D transform(Amg::Vector3D &v, Amg::Transform3D &tr)
Transform a point from a Trasformation3D.
Definition: GeoPrimitivesHelpers.h:156
Trk::CylinderVolumeBounds::halflengthZ
double halflengthZ() const
This method returns the halflengthZ.
Definition: CylinderVolumeBounds.h:207
AthMessaging
Class to provide easy MsgStream access and capabilities.
Definition: AthMessaging.h:55
Trk::VolumeIntersection::intersectApproximative
static std::pair< bool, std::unique_ptr< Trk::Volume > > intersectApproximative(const Volume &volA, const Volume &volB)
Definition: VolumeIntersection.cxx:416
Trk::Volume::clone
virtual Volume * clone() const
polymorpic deep copy
Definition: Volume.cxx:60
drawFromPickle.tan
tan
Definition: drawFromPickle.py:36
python.SystemOfUnits.rad
float rad
Definition: SystemOfUnits.py:126
Trk::PrismVolumeBounds
Definition: PrismVolumeBounds.h:44
Trk::VolumeConverter::resolveBooleanVolume
double resolveBooleanVolume(const Volume &trVol, double tolerance) const
Definition: VolumeConverter.cxx:226
Trk::Intersection::position
Amg::Vector3D position
Definition: Intersection.h:25
Trk::DoubleTrapezoidVolumeBounds::halflengthY2
double halflengthY2() const
This method returns the halflength2 in local y.
Definition: DoubleTrapezoidVolumeBounds.h:196
Trk::SubtractedVolumeBounds::outer
const Volume * outer() const
This method returns the outer Volume.
Definition: SubtractedVolumeBounds.h:111
beamspotman.dir
string dir
Definition: beamspotman.py:621
tolerance
Definition: suep_shower.h:17
Trk::CylinderVolumeBounds
Definition: CylinderVolumeBounds.h:70
Trk
Ensure that the ATLAS eigen extensions are properly loaded.
Definition: FakeTrackBuilder.h:9
SimplePolygonBrepVolumeBounds.h
DoubleTrapezoidVolumeBounds.h
Trk::PrismVolumeBounds::xyVertices
const std::vector< std::pair< double, double > > & xyVertices() const
This method returns the set of xy generating vertices.
Definition: PrismVolumeBounds.h:120
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:240
Trk::CylinderVolumeBounds::outerRadius
double outerRadius() const
This method returns the outer radius.
Definition: CylinderVolumeBounds.h:191
Trk::VolumeSpan::phiMin
double phiMin
Definition: VolumeConverter.h:34
Trk::Volume::transform
const Amg::Transform3D & transform() const
Return methods for geometry transform.
Definition: Volume.h:83
Trk::TrapezoidVolumeBounds
Definition: TrapezoidVolumeBounds.h:57
TRT::Hit::ident
@ ident
Definition: HitInfo.h:77
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition: GeoPrimitives.h:47
Trk::CylinderVolumeBounds::halfPhiSector
double halfPhiSector() const
This method returns the halfPhiSector angle.
Definition: CylinderVolumeBounds.h:203
Trk::BevelledCylinderVolumeBounds::innerRadius
double innerRadius() const
This method returns the inner radius.
Definition: BevelledCylinderVolumeBounds.h:237
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:105
Trk::VolumeSpan::zMax
double zMax
Definition: VolumeConverter.h:43
Trk::VolumeSpan::yMax
double yMax
Definition: VolumeConverter.h:41
Trk::MaterialProperties
Definition: MaterialProperties.h:40
TrackingVolume.h
Trk::DoubleTrapezoidVolumeBounds::minHalflengthX
double minHalflengthX() const
This method returns the X halflength at minimal Y.
Definition: DoubleTrapezoidVolumeBounds.h:180
mapkey::sf
@ sf
Definition: TElectronEfficiencyCorrectionTool.cxx:38
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
Trk::CylinderVolumeBounds::innerRadius
double innerRadius() const
This method returns the inner radius.
Definition: CylinderVolumeBounds.h:187
Trk::DoubleTrapezoidVolumeBounds::medHalflengthX
double medHalflengthX() const
This method returns the (maximal) halflength in local x.
Definition: DoubleTrapezoidVolumeBounds.h:184
GeoPrimitivesHelpers.h
unit
const PlainObject unit() const
This is a plugin that makes Eigen look like CLHEP & defines some convenience methods.
Definition: AmgMatrixBasePlugin.h:21
Trk::SubtractedVolumeBounds
Definition: SubtractedVolumeBounds.h:40
ReadCellNoiseFromCoolCompare.s2
s2
Definition: ReadCellNoiseFromCoolCompare.py:379
Amg::Translation3D
Eigen::Translation< double, 3 > Translation3D
Definition: GeoPrimitives.h:44
Trk::VolumeConverter::VolumePair
std::pair< std::shared_ptr< Volume >, std::shared_ptr< Volume > > VolumePair
Definition: VolumeConverter.h:70
Trk::Volume::volumeBounds
const VolumeBounds & volumeBounds() const
returns the volumeBounds()
Definition: Volume.h:96
Trk::DoubleTrapezoidVolumeBounds
Definition: DoubleTrapezoidVolumeBounds.h:66
python.changerun.pv
pv
Definition: changerun.py:79
CombinedVolumeBounds.h
VolumeConverter.h
python.DecayParser.children
children
Definition: DecayParser.py:32
Trk::VolumePart::parts
std::vector< std::shared_ptr< Volume > > parts
Definition: VolumeConverter.h:47
Trk::VolumeSpan::zMin
double zMin
Definition: VolumeConverter.h:42
Trk::VolumeConverter::translate
std::unique_ptr< TrackingVolume > translate(const GeoVPhysVol *gv, bool simplify, bool blend, double blendMassLimit) const
translation of GeoVPhysVol to Trk::TrackingVolume
Definition: VolumeConverter.cxx:42
Trk::phi
@ phi
Definition: ParamDefs.h:75
Trk::TrapezoidVolumeBounds::halflengthZ
double halflengthZ() const
This method returns the halflength in local z.
Definition: TrapezoidVolumeBounds.h:164
Trk::VolumeConverter::splitComposedVolume
static VolumePairVec splitComposedVolume(const Volume &trVol)
Decomposition of volume into set of non-overlapping subtractions from analytically calculable volume.
Definition: VolumeConverter.cxx:348
Trk::VolumeSpan
Definition: VolumeConverter.h:33
Trk::Material
Definition: Material.h:117
drawFromPickle.sin
sin
Definition: drawFromPickle.py:36
Trk::GeoMaterialConverter::dummy_material
static bool dummy_material(const GeoMaterial *)
hardcoded dummy materials : TODO : find generic criterium ( density ? radiation length ?...
Definition: GeoMaterialConverter.cxx:39
CaloLCW_tf.trf
trf
Definition: CaloLCW_tf.py:20
Trk::VolumeSpan::rMin
double rMin
Definition: VolumeConverter.h:36
Trk::CombinedVolumeBounds
Definition: CombinedVolumeBounds.h:44
Trk::BevelledCylinderVolumeBounds::halflengthZ
double halflengthZ() const
This method returns the halflengthZ.
Definition: BevelledCylinderVolumeBounds.h:242
Trk::CombinedVolumeBounds::second
const Volume * second() const
This method returns the second VolumeBounds.
Definition: CombinedVolumeBounds.h:123
Trk::BevelledCylinderVolumeBounds::halfPhiSector
double halfPhiSector() const
This method returns the halfPhiSector angle.
Definition: BevelledCylinderVolumeBounds.h:241
Trk::Volume
Definition: Volume.h:36
pow
constexpr int pow(int base, int exp) noexcept
Definition: ap_fixedTest.cxx:15
PlotCalibFromCool.vv
vv
Definition: PlotCalibFromCool.py:716
Trk::SimplePolygonBrepVolumeBounds::xyVertices
const std::vector< std::pair< double, double > > & xyVertices() const
This method returns the set of xy generating vertices.
Definition: SimplePolygonBrepVolumeBounds.h:147
Trk::VolumeConverter::VolumeConverter
VolumeConverter()
Definition: VolumeConverter.cxx:40
CuboidVolumeBounds.h
Trk::SubtractedVolumeBounds::inner
const Volume * inner() const
This method returns the inner Volume.
Definition: SubtractedVolumeBounds.h:114
Trk::v
@ v
Definition: ParamDefs.h:78
Trk::VolumeSpan::xMin
double xMin
Definition: VolumeConverter.h:38
SubtractedVolumeBounds.h