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"
35 return std::make_unique<Amg::Transform3D>(
trf);
43 bool simplify,
bool blend,
44 double blendMassLimit)
const {
46 const std::string
name = gv->getLogVol()->getName();
50 gv->getLogVol()->getShape(),
ident)};
64 double volumeFromGeoShape =
72 double volume = volumeFromGeoShape >= 0 ? volumeFromGeoShape : -1.;
73 if ((simplify || blend) && volumeFromGeoShape < 0) {
88 if (constituents.size() == 1 && !constituents[0].second) {
89 return std::make_unique<TrackingVolume>(*volGeo,
mat,
nullptr,
nullptr,
93 std::unique_ptr<Volume> envelope{};
94 std::string envName =
name;
96 std::unique_ptr<TrackingVolume> trEnv{};
100 if (constituents.size() == 1) {
102 envelope = std::make_unique<Volume>(*(constituents.front().first),
103 volGeo->transform());
106 if (blend && volume > 0 && volEnv > 0 &&
107 volume *
mat.rho < blendMassLimit)
110 if ((simplify || blended) && volume > 0 && volEnv > 0) {
112 double fraction = volume / volEnv;
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);
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);
131 std::unique_ptr<VolumeSpan>
span =
135 for (
const auto& fv : constituents) {
138 &(fv.first->volumeBounds()));
146 <<
"envelope estimate: object contains cylinder:"
147 <<
name <<
":" << isCyl);
149 <<
"complex volume span for envelope:" <<
name
150 <<
":x range:" << (*span).xMin <<
","
153 <<
"complex volume span for envelope:" <<
name
154 <<
":y range:" << (*span).yMin <<
","
157 <<
"complex volume span for envelope:" <<
name
158 <<
":z range:" << (*span).zMin <<
","
161 <<
"complex volume span for envelope:" <<
name
162 <<
":R range:" << (*span).rMin <<
","
165 <<
"complex volume span for envelope:" <<
name
166 <<
":phi range:" << (*span).phiMin <<
","
172 0.5 * ((*span).yMin + (*span).yMax),
173 0.5 * ((*span).zMin + (*span).zMax)}};
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>(
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{};
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));
195 cylBounds = std::make_shared<CylinderVolumeBounds>(
196 (*span).rMin, (*span).rMax,
197 0.5 * ((*span).zMax - (*span).zMin));
199 envelope = std::make_unique<Volume>(
205 if (blend && volume > 0 && volEnv > 0 && volume *
mat.rho < blendMassLimit)
208 if ((simplify || blended) && volume > 0 && volEnv > 0) {
209 double fraction = volume / volEnv;
212 if (blend && !blended)
213 envName = envName +
"_PERM";
214 trEnv = std::make_unique<TrackingVolume>(*envelope,
mat,
nullptr,
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);
231 inputVol.
parts.push_back(std::make_unique<Volume>(trVol));
232 constituents.push_back(std::move(inputVol));
237 while (sIter != constituents.end()) {
239 for (
unsigned int ii = 0; ii < (*sIter).parts.size(); ++ii) {
240 const VolumeBounds& bounds{((*sIter).parts[ii]->volumeBounds())};
246 (*sIter).parts[ii].reset(comb->
first()->
clone());
248 constituents.push_back(vp);
249 constituents.back().parts[ii].reset(comb->
second()->
clone());
250 constituents.push_back(vp);
251 constituents.back().parts.emplace_back(comb->
second()->
clone());
252 constituents.back().sign = -1. * constituents.back().sign;
256 (*sIter).parts[ii].reset(sub->
outer()->
clone());
260 volume += -1. * (*sIter).sign * volSub;
262 constituents.emplace_back(*sIter);
263 constituents.back().parts.emplace_back(
265 constituents.back().sign = -1. * constituents.back().sign;
273 sIter=constituents.erase(sIter);
280 sIter = constituents.begin();
281 else if ((*sIter).parts.size() == 1) {
283 volume += (*sIter).sign * volSingle;
284 sIter=constituents.erase(sIter);
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);
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();
303 sIter=constituents.erase(sIter);
308 sIter=constituents.erase(sIter);
312 std::pair<bool, std::unique_ptr<Volume>> overlap =
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();
321 sIter=constituents.erase(sIter);
330 }
else if ((*sIter).parts.size() == 1) {
332 volume += (*sIter).sign * volSingle;
333 sIter=constituents.erase(sIter);
340 if (!constituents.empty()) {
342 << constituents.size() <<
" items "
343 <<
":volume estimate:" << volume);
352 constituents.emplace_back(std::make_unique<Volume>(trVol),
nullptr);
354 std::shared_ptr<VolumeBounds> newBounds{};
355 while (sIter != constituents.end()) {
359 &((*sIter).first->volumeBounds()));
362 &((*sIter).first->volumeBounds()));
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()};
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));
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));
383 sIter, std::make_pair(combFirst, std::move(newSubVol)));
386 constituents.insert(sIter, std::make_pair(combFirst, 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));
395 std::make_pair(combSecond, std::move(newSubVol)));
397 constituents.insert(sIter,
398 std::make_pair(combSecond, combFirst));
401 sIter = constituents.begin();
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()};
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);
414 sIter, std::make_pair(outerVol, std::move(newSubVol)));
416 constituents.insert(sIter, std::make_pair(outerVol, innerVol));
418 sIter = constituents.begin();
428 double zTol,
double phiTol)
const {
474 if ((*s1).phiMin < (*s1).phiMax && (*s2).phiMin < (*s2).phiMax) {
477 }
else if ((*s1).phiMin < (*s1).phiMax && (*s2).phiMin > (*s2).phiMax) {
478 if ((*s1).phiMin > (*s2).phiMax) {
480 scomb.
phiMax = (*s2).phiMax;
481 }
else if ((*s1).phiMax < (*s2).phiMin) {
482 scomb.
phiMin = (*s2).phiMin;
488 }
else if ((*s1).phiMin > (*s1).phiMax && (*s2).phiMin < (*s2).phiMax) {
489 if ((*s2).phiMin > (*s1).phiMax) {
491 scomb.
phiMax = (*s1).phiMax;
492 }
else if ((*s2).phiMax < (*s1).phiMin) {
493 scomb.
phiMin = (*s1).phiMin;
503 return std::make_unique<VolumeSpan>(scomb);
509 double minPhi{2 *
M_PI};
519 std::vector<Amg::Vector3D> vtx;
520 std::vector<std::pair<int, int>> edges;
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);
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);
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);
627 edges.emplace_back(0, 1);
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);
666 if (bcyl->
type() == 2 || bcyl->
type() == 3) {
667 edges.emplace_back(2, 4);
668 edges.emplace_back(6, 8);
676 edges.emplace_back(0, 1);
701 edges.emplace_back(2, 3);
702 edges.emplace_back(4, 5);
703 edges.emplace_back(6, 7);
704 edges.emplace_back(8, 9);
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) {
716 vtx.size() - 4, vtx.size() - 2);
718 vtx.size() - 3, vtx.size() - 1);
720 if (vtx.size() > 4) {
721 edges.emplace_back(vtx.size() - 2, 1);
722 edges.emplace_back(vtx.size() - 1, 0);
725 edges.emplace_back(0, vtx.size() - 2);
726 edges.emplace_back(1, vtx.size() - 1);
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) {
737 vtx.size() - 4, vtx.size() - 2);
739 vtx.size() - 3, vtx.size() - 1);
742 edges.emplace_back(0, vtx.size() - 2);
743 edges.emplace_back(1, vtx.size() - 1);
746 std::vector<Amg::Vector3D> vtxt;
748 for (
unsigned int ie = 0;
ie < vtx.size();
ie++) {
753 double rad = gp.perp();
774 (vtxt[edges[0].first] - vtxt[edges[0].second]).
unit();
775 maxZ += ro *
sin(
dir.theta());
776 minZ += -ro *
sin(
dir.theta());
782 double le = (vtxt[0] - vtxt[1]).
norm();
783 if ((closest.
position - vtxt[0]).norm() < le &&
785 if (minR > closest.
position.perp() - ro)
789 if (phiClosest < minPhi || phiClosest > maxPhi) {
790 double phiTmp = minPhi;
795 minR =
std::max(0., minR - ro * std::abs(
dir.z()));
797 const double aTan = std::atan2(ro, minR);
802 if (maxPhi > 2 *
M_PI)
805 maxR += ro * std::abs(
cos(
dir.theta()));
816 for (
unsigned int ie = 0;
ie < edges.size();
ie++) {
818 (vtxt[edges[
ie].first] - vtxt[edges[
ie].second]).
unit();
822 (vtxt[edges[
ie].first] - vtxt[edges[
ie].second]).
norm();
823 if ((closest.
position - vtxt[edges[
ie].first]).norm() < le &&
829 if (vtxt.size() > 10) {
850 if (vtxt[10].
phi() +
M_PI < minPhi ||
851 vtxt[10].
phi() +
M_PI > maxPhi) {
855 for (
unsigned int iv = 2; iv < vtxt.size(); iv++) {
856 phiTmp = vtxt[iv].phi() +
M_PI;
859 minPhi = phiTmp < minPhi ? phiTmp : minPhi;
860 maxPhi = phiTmp > maxPhi ? phiTmp : maxPhi;
862 if (minPhi > 2 *
M_PI)
864 if (maxPhi > 2 *
M_PI)
872 if (minPhi >= maxPhi && (minPhi - maxPhi) <
M_PI) {
882 for (
unsigned int ie = 0;
ie < edges.size();
ie++) {
884 (vtxt[edges[
ie].first] - vtxt[edges[
ie].second]).
unit();
887 double le = (vtxt[edges[
ie].first] - vtxt[edges[
ie].second]).
norm();
888 if ((closest.
position - vtxt[edges[
ie].first]).norm() < le &&
894 if (std::abs(maxPhi - minPhi) >
M_PI) {
895 double phiTmp = minPhi;
898 for (
unsigned int iv = 0; iv < vtxt.size(); iv++) {
899 phiTmp = vtxt[iv].phi() +
M_PI;
902 minPhi = phiTmp < minPhi ? phiTmp : minPhi;
903 maxPhi = phiTmp > maxPhi ? phiTmp : maxPhi;
905 if (minPhi > 2 *
M_PI)
907 if (maxPhi > 2 *
M_PI)
909 if (minPhi >= maxPhi && (minPhi - maxPhi) <
M_PI) {
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
921 span.phiMin = minPhi;
922 maxPhi = (maxPhi + phiTol) > 2 *
M_PI ? maxPhi + phiTol - 2 *
M_PI
924 span.phiMax = maxPhi;
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;
934 return std::make_unique<VolumeSpan>(
span);
938 double precision)
const {
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);
992 vv +=
v.back().first * (
v[0].second -
v[
v.size() - 2].second);
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);
1001 vv +=
v.back().first * (
v[0].second -
v[
v.size() - 2].second);
1015 double precision)
const {
1020 if (sub.first && !sub.second)
1022 double fraction = -1.;
1024 std::pair<bool, std::unique_ptr<Volume>> overlap =
1027 if (overlap.first && !overlap.second)
1028 return fraction = 1.;
1029 else if (overlap.first && overlap.second) {
1039 if ((volA > 0 && volA < precision) || (volB > 0 && volB < precision))
1057 std::vector<MaterialComponent> materialContent;
1060 for (
const auto&
mat : materialContent) {
1063 double d =
sf > 0 ?
mat.second /
sf : 0.;
1066 (
mat.first.X0 > 0 ?
d /
mat.first.X0 : 0.));
1071 const GeoVPhysVol* gv,
1072 std::vector<MaterialComponent>& materialContent)
const {
1081 const GeoLogVol* lv = gv->getLogVol();
1084 double motherVolume = 0.;
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;
1095 sh && (
sh->type() ==
"Subtraction" ||
sh->type() ==
"Union" ||
1096 sh->type() ==
"Intersection");
1100 std::unique_ptr<Volume> vol{
1104 if (motherVolume < 0) {
1108 motherVolume = lv->getShape()->volume();
1111 double childVol = 0;
1112 std::string cPrevious =
" ";
1113 size_t nIdentical = 0;
1114 std::vector<Trk::MaterialComponent> childMat;
1116 for (
const GeoVPhysVol* cv :
children) {
1117 std::string cname = cv->getLogVol()->getName();
1118 if (cname == cPrevious)
1122 for (
const auto& cmat : childMat) {
1123 materialContent.emplace_back(cmat.first, nIdentical * cmat.second);
1124 childVol += materialContent.back().second;
1131 for (
const auto& cmat : childMat) {
1132 materialContent.emplace_back(cmat.first, nIdentical * cmat.second);
1133 childVol += materialContent.back().second;
1135 if (motherVolume > 0 && childVol > 0)
1136 motherVolume += -1. * childVol;
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);
1147 if (
sh->type() ==
"Subtraction") {
1148 const GeoShapeSubtraction* sub =
1149 dynamic_cast<const GeoShapeSubtraction*
>(
sh);
1153 if (
sh->type() ==
"Union") {
1154 const GeoShapeUnion* uni =
dynamic_cast<const GeoShapeUnion*
>(
sh);
1158 if (
sh->type() ==
"Intersection") {
1159 const GeoShapeIntersection* intr =
1160 dynamic_cast<const GeoShapeIntersection*
>(
sh);
1165 if (
sh->type() ==
"Shift") {
1166 const GeoShapeShift* shift =
dynamic_cast<const GeoShapeShift*
>(
sh);
1171 return sh->volume();