ATLAS Offline Software
Loading...
Searching...
No Matches
Trk::VolumeConverter Class Reference

A Simple Helper Class that collects methods for material simplification. More...

#include <VolumeConverter.h>

Inheritance diagram for Trk::VolumeConverter:
Collaboration diagram for Trk::VolumeConverter:

Public Types

using VolumePair
using VolumePairVec = std::vector<VolumePair>

Public Member Functions

 VolumeConverter ()
std::unique_ptr< TrackingVolumetranslate (const GeoVPhysVol *gv, bool simplify, bool blend, double blendMassLimit) const
 translation of GeoVPhysVol to Trk::TrackingVolume
double resolveBooleanVolume (const Volume &trVol, double tolerance) const
std::unique_ptr< VolumeSpanfindVolumeSpan (const VolumeBounds &volBounds, const Amg::Transform3D &transform, double zTol, double phiTol) const
 Estimation of the geometrical volume span.
double calculateVolume (const Volume &vol, bool nonBooleanOnly=false, double precision=1.e-3) const
 Volume calculation : by default return analytical solution only.
double estimateFraction (const VolumePair &sub, double precision) const
 the tricky part of volume calculation
void collectMaterial (const GeoVPhysVol *pv, Trk::MaterialProperties &layMat, double sf) const
 material collection for layers
void collectMaterialContent (const GeoVPhysVol *gv, std::vector< Trk::MaterialComponent > &materialContent) const
 material collection for volumes
bool msgLvl (const MSG::Level lvl) const
 Test the output level.
MsgStream & msg () const
 The standard message stream.
MsgStream & msg (const MSG::Level lvl) const
 The standard message stream.
void setLevel (MSG::Level lvl)
 Change the current logging level.

Static Public Member Functions

static VolumePairVec splitComposedVolume (const Volume &trVol)
 Decomposition of volume into set of non-overlapping subtractions from analytically calculable volume.

Private Member Functions

double leadingVolume (const GeoShape *sh) const
void initMessaging () const
 Initialize our message level and MessageSvc.

Private Attributes

Trk::GeoShapeConverter m_geoShapeConverter
 shape converter
std::string m_nm
 Message source name.
boost::thread_specific_ptr< MsgStream > m_msg_tls
 MsgStream instance (a std::cout like with print-out levels).
std::atomic< IMessageSvc * > m_imsg { nullptr }
 MessageSvc pointer.
std::atomic< MSG::Level > m_lvl { MSG::NIL }
 Current logging level.
std::atomic_flag m_initialized ATLAS_THREAD_SAFE = ATOMIC_FLAG_INIT
 Messaging initialized (initMessaging).

Static Private Attributes

static constexpr double s_precisionInX0

Detailed Description

A Simple Helper Class that collects methods for material simplification.

Author
sarka.nosp@m..tod.nosp@m.orova.nosp@m.@cer.nosp@m.n.ch

Definition at line 60 of file VolumeConverter.h.

Member Typedef Documentation

◆ VolumePair

Initial value:
std::pair<std::shared_ptr<Volume>, std::shared_ptr<Volume>>

Definition at line 69 of file VolumeConverter.h.

◆ VolumePairVec

Definition at line 71 of file VolumeConverter.h.

Constructor & Destructor Documentation

◆ VolumeConverter()

Trk::VolumeConverter::VolumeConverter ( )

Definition at line 44 of file VolumeConverter.cxx.

44: AthMessaging("VolumeConverter") {}
AthMessaging()
Default constructor:

Member Function Documentation

◆ calculateVolume()

double Trk::VolumeConverter::calculateVolume ( const Volume & vol,
bool nonBooleanOnly = false,
double precision = 1.e-3 ) const

Volume calculation : by default return analytical solution only.

Definition at line 941 of file VolumeConverter.cxx.

942 {
943
944 double volume = -1.;
945
946 const CylinderVolumeBounds* cyl = dynamic_cast<const CylinderVolumeBounds*>(&(vol.volumeBounds()));
947 const CuboidVolumeBounds* box = dynamic_cast<const CuboidVolumeBounds*>(&(vol.volumeBounds()));
948 const TrapezoidVolumeBounds* trd = dynamic_cast<const TrapezoidVolumeBounds*>(&(vol.volumeBounds()));
949 const BevelledCylinderVolumeBounds* bcyl =dynamic_cast<const BevelledCylinderVolumeBounds*>(&(vol.volumeBounds()));
950 const PrismVolumeBounds* prism =dynamic_cast<const PrismVolumeBounds*>(&(vol.volumeBounds()));
951 const SimplePolygonBrepVolumeBounds* spb = dynamic_cast<const SimplePolygonBrepVolumeBounds*>(&(vol.volumeBounds()));
952 const CombinedVolumeBounds* comb =dynamic_cast<const CombinedVolumeBounds*>(&(vol.volumeBounds()));
953 const SubtractedVolumeBounds* sub = dynamic_cast<const SubtractedVolumeBounds*>(&(vol.volumeBounds()));
954
955 if (cyl) {
956 return 2 * cyl->halfPhiSector() * cyl->halflengthZ() *
957 (std::pow(cyl->outerRadius(), 2) -
958 std::pow(cyl->innerRadius(), 2));
959 }
960 if (box) {
961 return 8 * box->halflengthX() * box->halflengthY() * box->halflengthZ();
962 }
963 if (trd) {
964 return 4 * (trd->minHalflengthX() + trd->maxHalflengthX()) *
965 trd->halflengthY() * trd->halflengthZ();
966 }
967 if (bcyl) {
968 int type = bcyl->type();
969 if (type < 1)
970 return 2 * bcyl->halfPhiSector() *
971 (std::pow(bcyl->outerRadius(), 2) -
972 std::pow(bcyl->innerRadius(), 2)) *
973 bcyl->halflengthZ();
974 if (type == 1)
975 return 2 * bcyl->halflengthZ() *
976 (bcyl->halfPhiSector() * std::pow(bcyl->outerRadius(), 2) -
977 std::pow(bcyl->innerRadius(), 2) *
978 std::tan(bcyl->halfPhiSector()));
979 if (type == 2)
980 return 2 * bcyl->halflengthZ() *
981 (-bcyl->halfPhiSector() * std::pow(bcyl->innerRadius(), 2) +
982 std::pow(bcyl->outerRadius(), 2) *
983 std::tan(bcyl->halfPhiSector()));
984 if (type == 3)
985 return 2 * bcyl->halflengthZ() * std::tan(bcyl->halfPhiSector()) *
986 (std::pow(bcyl->outerRadius(), 2) -
987 std::pow(bcyl->innerRadius(), 2));
988 }
989 if (prism) {
990
991 std::vector<std::pair<double, double>> v = prism->xyVertices();
992 double vv = v[0].first * (v[1].second - v.back().second);
993 for (unsigned int i = 1; i < v.size() - 1; i++) {
994 vv += v[i].first * (v[i + 1].second - v[i - 1].second);
995 }
996 vv += v.back().first * (v[0].second - v[v.size() - 2].second);
997 return vv * prism->halflengthZ();
998 }
999 if (spb) {
1000 std::vector<std::pair<double, double>> v = spb->xyVertices();
1001 double vv = v[0].first * (v[1].second - v.back().second);
1002 for (unsigned int i = 1; i < v.size() - 1; i++) {
1003 vv += v[i].first * (v[i + 1].second - v[i - 1].second);
1004 }
1005 vv += v.back().first * (v[0].second - v[v.size() - 2].second);
1006 return vv * spb->halflengthZ();
1007 }
1008
1009 if (nonBooleanOnly)
1010 return volume;
1011
1012 if (comb || sub) {
1013 return resolveBooleanVolume(vol, precision);
1014 }
1015 return volume;
1016}
double resolveBooleanVolume(const Volume &trVol, double tolerance) const
@ v
Definition ParamDefs.h:78

◆ collectMaterial()

void Trk::VolumeConverter::collectMaterial ( const GeoVPhysVol * pv,
Trk::MaterialProperties & layMat,
double sf ) const

material collection for layers

Definition at line 1049 of file VolumeConverter.cxx.

1051 {
1052 // sf is the area of the layer collecting the material
1053
1054 // solution relying on GeoModel
1055 // currently involves hit&miss on-fly calculation of boolean volumes
1056 // GeoModelTools::MaterialComponent mat =
1057 // gm_materialHelper.collectMaterial(pv); Material newMP =
1058 // convert(mat.first); double d = mat.second / sf; layMat.addMaterial(newMP,
1059 // d / newMP.x0()); return;
1060
1061 std::vector<MaterialComponent> materialContent;
1062 collectMaterialContent(pv, materialContent);
1063
1064 for (const auto& mat : materialContent) {
1065 if (mat.second < 0)
1066 continue; // protection unsolved booleans
1067 double d = sf > 0 ? mat.second / sf : 0.;
1068 if (d > 0)
1069 layMat.addMaterial(mat.first,
1070 (mat.first.X0 > 0 ? d / mat.first.X0 : 0.));
1071 }
1072}
void addMaterial(const Material &mp, float dInX0)
Material averaging.
void collectMaterialContent(const GeoVPhysVol *gv, std::vector< Trk::MaterialComponent > &materialContent) const
material collection for volumes

◆ collectMaterialContent()

void Trk::VolumeConverter::collectMaterialContent ( const GeoVPhysVol * gv,
std::vector< Trk::MaterialComponent > & materialContent ) const

material collection for volumes

Definition at line 1074 of file VolumeConverter.cxx.

1076 {
1077
1078 // solution relying on GeoModel
1079 // currently involves hit&miss on-fly calculation of boolean volumes
1080 // GeoModelTools::MaterialComponent mat =
1081 // gm_materialHelper.collectMaterial(pv); Material newMP =
1082 // convert(mat.first); materialContent.push_back( MaterialComponent( newMP,
1083 // mat.second) ); return;
1084
1085 const GeoLogVol* lv = gv->getLogVol();
1086 Material mat = Trk::GeoMaterialConverter::convert(lv->getMaterial());
1087
1088 double motherVolume = 0.;
1089
1090 // skip volume calculation for dummy material configuration
1091 if (!Trk::GeoMaterialConverter::dummy_material(lv->getMaterial())) {
1092 const GeoShape* sh = lv->getShape();
1093 while (sh && sh->type() == "Shift") {
1094 const GeoShapeShift* shift = dynamic_cast<const GeoShapeShift*>(sh);
1095 sh = shift ? shift->getOp() : nullptr;
1096 }
1097
1098 bool isBoolean =
1099 sh && (sh->type() == "Subtraction" || sh->type() == "Union" ||
1100 sh->type() == "Intersection");
1101
1102 if (isBoolean) {
1103 Amg::Transform3D transf{Amg::Transform3D::Identity()};
1104 std::unique_ptr<Volume> vol{
1105 m_geoShapeConverter.translateGeoShape(sh, transf)};
1106 motherVolume =
1107 calculateVolume(*vol, false, std::pow(1.e-3 * mat.X0, 3));
1108 if (motherVolume < 0) {
1109 // m_geoShapeConverter.decodeShape(sh);
1110 }
1111 } else
1112 motherVolume = lv->getShape()->volume();
1113 }
1114
1115 double childVol = 0;
1116 std::string cPrevious = " ";
1117 size_t nIdentical = 0;
1118 std::vector<Trk::MaterialComponent> childMat;
1119 std::vector<const GeoVPhysVol*> children = geoGetVolumesNoXform (gv);
1120 for (const GeoVPhysVol* cv : children) {
1121 std::string cname = cv->getLogVol()->getName();
1122 if (cname == cPrevious)
1123 nIdentical++; // assuming identity for identical name and branching
1124 // history
1125 else { // scale and collect material from previous item
1126 for (const auto& cmat : childMat) {
1127 materialContent.emplace_back(cmat.first, nIdentical * cmat.second);
1128 childVol += materialContent.back().second;
1129 }
1130 childMat.clear(); // reset
1131 nIdentical = 1; // current
1132 collectMaterialContent(cv, childMat);
1133 }
1134 }
1135 for (const auto& cmat : childMat) {
1136 materialContent.emplace_back(cmat.first, nIdentical * cmat.second);
1137 childVol += materialContent.back().second;
1138 }
1139 if (motherVolume > 0 && childVol > 0)
1140 motherVolume += -1. * childVol;
1141
1142 ATH_MSG_DEBUG("collected material:" << lv->getName() << ":made of:"
1143 << lv->getMaterial()->getName()
1144 << ":density(g/mm3)" << mat.rho
1145 << ":mass:" << mat.rho * motherVolume);
1146 materialContent.emplace_back(mat, motherVolume);
1147}
#define ATH_MSG_DEBUG(x)
std::vector< const GeoVPhysVol * > geoGetVolumesNoXform(const GeoGraphNode *node, int depthLimit=1, int sizeHint=20)
Return the child volumes.
@ Material
static bool dummy_material(const GeoMaterial *)
hardcoded dummy materials : TODO : find generic criterium ( density ?
static Material convert(const GeoMaterial *gm)
Single conversion , input type GeoMaterial - output type Trk::MaterialProperties.
Trk::GeoShapeConverter m_geoShapeConverter
shape converter
double calculateVolume(const Volume &vol, bool nonBooleanOnly=false, double precision=1.e-3) const
Volume calculation : by default return analytical solution only.
Eigen::Affine3d Transform3D
virtual void shift(size_t pos, ptrdiff_t offs) override
Shift the elements of the container.

◆ estimateFraction()

double Trk::VolumeConverter::estimateFraction ( const VolumePair & sub,
double precision ) const

the tricky part of volume calculation

Definition at line 1018 of file VolumeConverter.cxx.

1019 {
1020
1021 if (!sub.first)
1022 return 0.;
1023
1024 if (sub.first && !sub.second)
1025 return 1.;
1026 double fraction = -1.;
1027
1028 std::pair<bool, std::unique_ptr<Volume>> overlap =
1029 Trk::VolumeIntersection::intersect(*sub.first, *sub.second);
1030
1031 if (overlap.first && !overlap.second)
1032 return fraction = 1.;
1033 else if (overlap.first && overlap.second) {
1034 fraction = 1. - calculateVolume(*overlap.second, true, precision) /
1035 calculateVolume(*sub.first, true, precision);
1036 return fraction;
1037 }
1038 // resolve embedded volumes
1039
1040 // trivial within required precision
1041 double volA = calculateVolume(*sub.first, true, precision);
1042 double volB = calculateVolume(*sub.second, true, precision);
1043 if ((volA > 0 && volA < precision) || (volB > 0 && volB < precision))
1044 return 1.;
1045
1046 return fraction;
1047}
static std::pair< bool, std::unique_ptr< Trk::Volume > > intersect(const Volume &volA, const Volume &volB)

◆ findVolumeSpan()

std::unique_ptr< VolumeSpan > Trk::VolumeConverter::findVolumeSpan ( const VolumeBounds & volBounds,
const Amg::Transform3D & transform,
double zTol,
double phiTol ) const

Estimation of the geometrical volume span.

Definition at line 430 of file VolumeConverter.cxx.

432 {
433 // volume shape
434 const CuboidVolumeBounds* box =
435 dynamic_cast<const CuboidVolumeBounds*>(&volBounds);
436 const TrapezoidVolumeBounds* trd =
437 dynamic_cast<const TrapezoidVolumeBounds*>(&volBounds);
438 const DoubleTrapezoidVolumeBounds* dtrd =
439 dynamic_cast<const DoubleTrapezoidVolumeBounds*>(&volBounds);
440 const BevelledCylinderVolumeBounds* bcyl =
441 dynamic_cast<const BevelledCylinderVolumeBounds*>(&volBounds);
442 const CylinderVolumeBounds* cyl =
443 dynamic_cast<const CylinderVolumeBounds*>(&volBounds);
444 const SubtractedVolumeBounds* sub =
445 dynamic_cast<const SubtractedVolumeBounds*>(&volBounds);
446 const CombinedVolumeBounds* comb =
447 dynamic_cast<const CombinedVolumeBounds*>(&volBounds);
448 const SimplePolygonBrepVolumeBounds* spb =
449 dynamic_cast<const SimplePolygonBrepVolumeBounds*>(&volBounds);
450 const PrismVolumeBounds* prism =
451 dynamic_cast<const PrismVolumeBounds*>(&volBounds);
452
453 double dPhi = 0.;
454
455 if (sub) {
456 return findVolumeSpan(sub->outer()->volumeBounds(),
457 transform * sub->outer()->transform(), zTol,
458 phiTol);
459 }
460
461 if (comb) {
462 std::unique_ptr<VolumeSpan> s1 = findVolumeSpan(
463 comb->first()->volumeBounds(),
464 transform * comb->first()->transform(), zTol, phiTol);
465 std::unique_ptr<VolumeSpan> s2 = findVolumeSpan(
466 comb->second()->volumeBounds(),
467 transform * comb->second()->transform(), zTol, phiTol);
468
469 VolumeSpan scomb;
470 scomb.rMin = std::min((*s1).rMin, (*s2).rMin);
471 scomb.rMax = std::max((*s1).rMax, (*s2).rMax);
472 scomb.xMin = std::min((*s1).xMin, (*s2).xMin);
473 scomb.xMax = std::max((*s1).xMax, (*s2).xMax);
474 scomb.yMin = std::min((*s1).yMin, (*s2).yMin);
475 scomb.yMax = std::max((*s1).yMax, (*s2).yMax);
476 scomb.zMin = std::min((*s1).zMin, (*s2).zMin);
477 scomb.zMax = std::max((*s1).zMax, (*s2).zMax);
478 if ((*s1).phiMin < (*s1).phiMax && (*s2).phiMin < (*s2).phiMax) {
479 scomb.phiMin = std::min((*s1).phiMin, (*s2).phiMin);
480 scomb.phiMax = std::max((*s1).phiMax, (*s2).phiMax);
481 } else if ((*s1).phiMin < (*s1).phiMax && (*s2).phiMin > (*s2).phiMax) {
482 if ((*s1).phiMin > (*s2).phiMax) {
483 scomb.phiMin = std::min((*s1).phiMin, (*s2).phiMin);
484 scomb.phiMax = (*s2).phiMax;
485 } else if ((*s1).phiMax < (*s2).phiMin) {
486 scomb.phiMin = (*s2).phiMin;
487 scomb.phiMax = std::max((*s1).phiMax, (*s2).phiMax);
488 } else {
489 scomb.phiMin = 0.;
490 scomb.phiMax = 2 * M_PI;
491 }
492 } else if ((*s1).phiMin > (*s1).phiMax && (*s2).phiMin < (*s2).phiMax) {
493 if ((*s2).phiMin > (*s1).phiMax) {
494 scomb.phiMin = std::min((*s1).phiMin, (*s2).phiMin);
495 scomb.phiMax = (*s1).phiMax;
496 } else if ((*s2).phiMax < (*s1).phiMin) {
497 scomb.phiMin = (*s1).phiMin;
498 scomb.phiMax = std::max((*s1).phiMax, (*s2).phiMax);
499 } else {
500 scomb.phiMin = 0.;
501 scomb.phiMax = 2 * M_PI;
502 }
503 } else {
504 scomb.phiMin = std::min((*s1).phiMin, (*s2).phiMin);
505 scomb.phiMax = std::max((*s1).phiMax, (*s2).phiMax);
506 }
507 return std::make_unique<VolumeSpan>(scomb);
508 }
509
510 //
511 double minZ{1.e6};
512 double maxZ{-1.e6};
513 double minPhi{2 * M_PI};
514 double maxPhi{0.};
515 double minR{1.e6};
516 double maxR{0.};
517 double minX{1.e6};
518 double maxX{-1.e6};
519 double minY{1.e6};
520 double maxY{-1.e6};
521
522 // defined vertices and edges
523 std::vector<Amg::Vector3D> vtx;
524 std::vector<std::pair<int, int>> edges;
525 VolumeSpan span;
526
527 if (box) {
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 vtx.emplace_back(box->halflengthX(), -box->halflengthY(),
541 -box->halflengthZ());
542 vtx.emplace_back(-box->halflengthX(), -box->halflengthY(),
543 -box->halflengthZ());
544 edges.emplace_back(0, 1);
545 edges.emplace_back(0, 2);
546 edges.emplace_back(1, 3);
547 edges.emplace_back(2, 3);
548 edges.emplace_back(4, 5);
549 edges.emplace_back(4, 6);
550 edges.emplace_back(5, 7);
551 edges.emplace_back(6, 7);
552 edges.emplace_back(0, 4);
553 edges.emplace_back(1, 5);
554 edges.emplace_back(2, 6);
555 edges.emplace_back(3, 7);
556 }
557 if (trd) {
558 vtx.emplace_back(trd->maxHalflengthX(), trd->halflengthY(),
559 trd->halflengthZ());
560 vtx.emplace_back(-trd->maxHalflengthX(), trd->halflengthY(),
561 trd->halflengthZ());
562 vtx.emplace_back(trd->minHalflengthX(), -trd->halflengthY(),
563 trd->halflengthZ());
564 vtx.emplace_back(-trd->minHalflengthX(), -trd->halflengthY(),
565 trd->halflengthZ());
566 vtx.emplace_back(trd->maxHalflengthX(), trd->halflengthY(),
567 -trd->halflengthZ());
568 vtx.emplace_back(-trd->maxHalflengthX(), trd->halflengthY(),
569 -trd->halflengthZ());
570 vtx.emplace_back(trd->minHalflengthX(), -trd->halflengthY(),
571 -trd->halflengthZ());
572 vtx.emplace_back(-trd->minHalflengthX(), -trd->halflengthY(),
573 -trd->halflengthZ());
574 edges.emplace_back(0, 1);
575 edges.emplace_back(0, 2);
576 edges.emplace_back(1, 3);
577 edges.emplace_back(2, 3);
578 edges.emplace_back(4, 5);
579 edges.emplace_back(4, 6);
580 edges.emplace_back(5, 7);
581 edges.emplace_back(6, 7);
582 edges.emplace_back(0, 4);
583 edges.emplace_back(1, 5);
584 edges.emplace_back(2, 6);
585 edges.emplace_back(3, 7);
586 }
587 if (dtrd) {
588 vtx.emplace_back(dtrd->maxHalflengthX(), 2 * dtrd->halflengthY2(),
589 dtrd->halflengthZ());
590 vtx.emplace_back(-dtrd->maxHalflengthX(), 2 * dtrd->halflengthY2(),
591 dtrd->halflengthZ());
592 vtx.emplace_back(dtrd->medHalflengthX(), 0., dtrd->halflengthZ());
593 vtx.emplace_back(-dtrd->medHalflengthX(), 0., dtrd->halflengthZ());
594 vtx.emplace_back(dtrd->minHalflengthX(), -2 * dtrd->halflengthY1(),
595 dtrd->halflengthZ());
596 vtx.emplace_back(-dtrd->minHalflengthX(), -2 * dtrd->halflengthY1(),
597 dtrd->halflengthZ());
598 vtx.emplace_back(dtrd->maxHalflengthX(), 2 * dtrd->halflengthY2(),
599 -dtrd->halflengthZ());
600 vtx.emplace_back(-dtrd->maxHalflengthX(), 2 * dtrd->halflengthY2(),
601 -dtrd->halflengthZ());
602 vtx.emplace_back(dtrd->medHalflengthX(), 0., -dtrd->halflengthZ());
603 vtx.emplace_back(-dtrd->medHalflengthX(), 0., -dtrd->halflengthZ());
604 vtx.emplace_back(dtrd->minHalflengthX(), -2 * dtrd->halflengthY1(),
605 -dtrd->halflengthZ());
606 vtx.emplace_back(-dtrd->minHalflengthX(), -2 * dtrd->halflengthY1(),
607 -dtrd->halflengthZ());
608 edges.emplace_back(0, 1);
609 edges.emplace_back(0, 2);
610 edges.emplace_back(1, 3);
611 edges.emplace_back(2, 4);
612 edges.emplace_back(3, 5);
613 edges.emplace_back(4, 5);
614 edges.emplace_back(6, 7);
615 edges.emplace_back(6, 8);
616 edges.emplace_back(7, 9);
617 edges.emplace_back(8, 10);
618 edges.emplace_back(9, 11);
619 edges.emplace_back(10, 11);
620 edges.emplace_back(0, 6);
621 edges.emplace_back(1, 7);
622 edges.emplace_back(2, 8);
623 edges.emplace_back(3, 9);
624 edges.emplace_back(4, 10);
625 edges.emplace_back(5, 11);
626 }
627 if (bcyl) {
628 dPhi = bcyl->halfPhiSector();
629 vtx.emplace_back(0., 0., bcyl->halflengthZ());
630 vtx.emplace_back(0., 0., -bcyl->halflengthZ());
631 edges.emplace_back(0, 1);
632 if (dPhi < M_PI) {
633 const double cosDphi = std::cos(dPhi);
634 const double sinDphi = std::sin(dPhi);
635 vtx.emplace_back(bcyl->outerRadius() * cosDphi,
636 bcyl->outerRadius() * sinDphi,
637 bcyl->halflengthZ());
638 vtx.emplace_back(bcyl->innerRadius() * cosDphi,
639 bcyl->innerRadius() * sinDphi,
640 bcyl->halflengthZ());
641 vtx.emplace_back(bcyl->outerRadius() * cosDphi,
642 -bcyl->outerRadius() * sinDphi,
643 bcyl->halflengthZ());
644 vtx.emplace_back(bcyl->innerRadius() * cosDphi,
645 -bcyl->innerRadius() * sinDphi,
646 bcyl->halflengthZ());
647 vtx.emplace_back(bcyl->outerRadius() * cosDphi,
648 bcyl->outerRadius() * sinDphi,
649 -bcyl->halflengthZ());
650 vtx.emplace_back(bcyl->innerRadius() * cosDphi,
651 bcyl->innerRadius() * sinDphi,
652 -bcyl->halflengthZ());
653 vtx.emplace_back(bcyl->outerRadius() * cosDphi,
654 -bcyl->outerRadius() * sinDphi,
655 -bcyl->halflengthZ());
656 vtx.emplace_back(bcyl->innerRadius() * cosDphi,
657 -bcyl->innerRadius() * sinDphi,
658 -bcyl->halflengthZ());
659 vtx.emplace_back(bcyl->outerRadius(), 0.,
660 0.); // to distinguish phi intervals for cylinders
661 // aligned with z axis
662 edges.emplace_back(2, 3);
663 edges.emplace_back(4, 5);
664 edges.emplace_back(6, 7);
665 edges.emplace_back(8, 9);
666 if (bcyl->type() == 1 || bcyl->type() == 3) {
667 edges.emplace_back(3, 5);
668 edges.emplace_back(7, 9);
669 }
670 if (bcyl->type() == 2 || bcyl->type() == 3) {
671 edges.emplace_back(2, 4);
672 edges.emplace_back(6, 8);
673 }
674 }
675 }
676 if (cyl) {
677 dPhi = cyl->halfPhiSector();
678 vtx.emplace_back(0., 0., cyl->halflengthZ());
679 vtx.emplace_back(0., 0., -cyl->halflengthZ());
680 edges.emplace_back(0, 1);
681 if (dPhi < M_PI) {
682 const double cosDphi = std::cos(dPhi);
683 const double sinDphi = std::sin(dPhi);
684 vtx.emplace_back(cyl->outerRadius() * cosDphi,
685 cyl->outerRadius() * sinDphi, cyl->halflengthZ());
686 vtx.emplace_back(cyl->innerRadius() * cosDphi,
687 cyl->innerRadius() * sinDphi, cyl->halflengthZ());
688 vtx.emplace_back(cyl->outerRadius() * cosDphi,
689 -cyl->outerRadius() * sinDphi, cyl->halflengthZ());
690 vtx.emplace_back(cyl->outerRadius() * cosDphi,
691 -cyl->outerRadius() * sinDphi, cyl->halflengthZ());
692 vtx.emplace_back(cyl->outerRadius() * cosDphi,
693 cyl->outerRadius() * sinDphi, -cyl->halflengthZ());
694 vtx.emplace_back(cyl->innerRadius() * cosDphi,
695 cyl->innerRadius() * sinDphi, -cyl->halflengthZ());
696 vtx.emplace_back(cyl->outerRadius() * cosDphi,
697 -cyl->outerRadius() * sinDphi,
698 -cyl->halflengthZ());
699 vtx.emplace_back(cyl->outerRadius() * cosDphi,
700 -cyl->outerRadius() * sinDphi,
701 -cyl->halflengthZ());
702 vtx.emplace_back(cyl->outerRadius(), 0.,
703 0.); // to distinguish phi intervals for cylinders
704 // aligned with z axis
705 edges.emplace_back(2, 3);
706 edges.emplace_back(4, 5);
707 edges.emplace_back(6, 7);
708 edges.emplace_back(8, 9);
709 }
710 }
711
712 if (spb) {
713 const std::vector<std::pair<double, double>> vtcs = spb->xyVertices();
714 for (const auto& vtc : vtcs) {
715 vtx.emplace_back(vtc.first, vtc.second, spb->halflengthZ());
716 vtx.emplace_back(vtc.first, vtc.second, -spb->halflengthZ());
717 edges.emplace_back(vtx.size() - 2, vtx.size() - 1);
718 if (vtx.size() > 2) {
719 edges.emplace_back(
720 vtx.size() - 4, vtx.size() - 2);
721 edges.emplace_back(
722 vtx.size() - 3, vtx.size() - 1);
723 }
724 if (vtx.size() > 4) { // some diagonals
725 edges.emplace_back(vtx.size() - 2, 1);
726 edges.emplace_back(vtx.size() - 1, 0);
727 }
728 }
729 edges.emplace_back(0, vtx.size() - 2);
730 edges.emplace_back(1, vtx.size() - 1);
731 }
732
733 if (prism) {
734 const std::vector<std::pair<double, double>> vtcs = prism->xyVertices();
735 for (const auto& vtc : vtcs) {
736 vtx.emplace_back(vtc.first, vtc.second, prism->halflengthZ());
737 vtx.emplace_back(vtc.first, vtc.second, -prism->halflengthZ());
738 edges.emplace_back(vtx.size() - 2, vtx.size() - 1);
739 if (vtx.size() > 2) {
740 edges.emplace_back(
741 vtx.size() - 4, vtx.size() - 2);
742 edges.emplace_back(
743 vtx.size() - 3, vtx.size() - 1);
744 }
745 }
746 edges.emplace_back(0, vtx.size() - 2);
747 edges.emplace_back(1, vtx.size() - 1);
748 }
749
750 std::vector<Amg::Vector3D> vtxt;
751
752 for (unsigned int ie = 0; ie < vtx.size(); ie++) {
753 Amg::Vector3D gp = transform * vtx[ie];
754 vtxt.push_back(gp);
755
756 double phi = gp.phi() + M_PI;
757 double rad = gp.perp();
758
759 // collect limits from vertices
760 minX = std::min(minX, gp[0]);
761 maxX = std::max(maxX, gp[0]);
762 minY = std::min(minY, gp[1]);
763 maxY = std::max(maxY, gp[1]);
764 minZ = std::min(minZ, gp[2]);
765 maxZ = std::max(maxZ, gp[2]);
766 minR = std::min(minR, rad);
767 maxR = std::max(maxR, rad);
768 maxPhi = std::max(maxPhi, phi);
769 minPhi = std::min(minPhi, phi);
770 }
771
772 if (cyl || bcyl) {
773
774 double ro = cyl ? cyl->outerRadius() : bcyl->outerRadius();
775 double ri = cyl ? cyl->innerRadius() : bcyl->innerRadius();
776 // z span corrected for theta inclination
778 (vtxt[edges[0].first] - vtxt[edges[0].second]).unit();
779 maxZ += ro * sin(dir.theta());
780 minZ += -ro * sin(dir.theta());
781 // azimuthal & radial extent
782 if (ro < minR) { // excentric object, phi span driven by z-R extent
783 // calculate point of closest approach
784 PerigeeSurface peri;
785 Intersection closest = peri.straightLineIntersection(vtxt[1], dir);
786 double le = (vtxt[0] - vtxt[1]).norm();
787 if ((closest.position - vtxt[0]).norm() < le &&
788 (closest.position - vtxt[1]).norm() < le) {
789 if (minR > closest.position.perp() - ro)
790 minR = std::max(0., closest.position.perp() - ro);
791 // use for phi check
792 double phiClosest = closest.position.phi() + M_PI;
793 if (phiClosest < minPhi || phiClosest > maxPhi) {
794 double phiTmp = minPhi;
795 minPhi = maxPhi;
796 maxPhi = phiTmp;
797 }
798 } else
799 minR = std::max(0., minR - ro * std::abs(dir.z()));
800
801 const double aTan = std::atan2(ro, minR);
802 minPhi += -aTan;
803 maxPhi += aTan;
804 if (minPhi < 0)
805 minPhi += 2 * M_PI;
806 if (maxPhi > 2 * M_PI)
807 maxPhi += -2 * M_PI;
808
809 maxR += ro * std::abs(cos(dir.theta()));
810 } else {
811
812 double rAx = std::max(vtxt[0].perp(), vtxt[1].perp());
813 if (rAx < ri)
814 minR = ri - rAx;
815 else
816 minR = std::max(0., minR - ro * std::abs(cos(dir.theta())));
817
818 // loop over edges to check inner radial extent
819 PerigeeSurface peri;
820 for (unsigned int ie = 0; ie < edges.size(); ie++) {
822 (vtxt[edges[ie].first] - vtxt[edges[ie].second]).unit();
823 Intersection closest =
824 peri.straightLineIntersection(vtxt[edges[ie].second], dir);
825 double le =
826 (vtxt[edges[ie].first] - vtxt[edges[ie].second]).norm();
827 if ((closest.position - vtxt[edges[ie].first]).norm() < le &&
828 (closest.position - vtxt[edges[ie].second]).norm() < le)
829 if (minR > closest.position.perp())
830 minR = closest.position.perp();
831 }
832
833 if (vtxt.size() > 10) { // cylindrical section
834 // find spread of phi extent at section (-) boundary : vertices
835 // 4,5,8,9
836 double phiSecLmin = std::min(
837 std::min(vtxt[4].phi() + M_PI, vtxt[5].phi() + M_PI),
838 std::min(vtxt[8].phi() + M_PI, vtxt[9].phi() + M_PI));
839 double phiSecLmax = std::max(
840 std::max(vtxt[4].phi() + M_PI, vtxt[5].phi() + M_PI),
841 std::max(vtxt[8].phi() + M_PI, vtxt[9].phi() + M_PI));
842 // find spread of phi extent at section (+) boundary : vertices
843 // 2,3,6,7
844 double phiSecUmin = std::min(
845 std::min(vtxt[2].phi() + M_PI, vtxt[3].phi() + M_PI),
846 std::min(vtxt[6].phi() + M_PI, vtxt[7].phi() + M_PI));
847 double phiSecUmax = std::max(
848 std::max(vtxt[2].phi() + M_PI, vtxt[3].phi() + M_PI),
849 std::max(vtxt[6].phi() + M_PI, vtxt[7].phi() + M_PI));
850 minPhi = std::min(std::min(phiSecLmin, phiSecLmax),
851 std::min(phiSecUmin, phiSecUmax));
852 maxPhi = std::max(std::max(phiSecLmin, phiSecLmax),
853 std::max(phiSecUmin, phiSecUmax));
854 if (vtxt[10].phi() + M_PI < minPhi ||
855 vtxt[10].phi() + M_PI > maxPhi) {
856 minPhi = 3 * M_PI;
857 maxPhi = 0.;
858 double phiTmp;
859 for (unsigned int iv = 2; iv < vtxt.size(); iv++) {
860 phiTmp = vtxt[iv].phi() + M_PI;
861 if (phiTmp < M_PI)
862 phiTmp += 2 * M_PI;
863 minPhi = phiTmp < minPhi ? phiTmp : minPhi;
864 maxPhi = phiTmp > maxPhi ? phiTmp : maxPhi;
865 }
866 if (minPhi > 2 * M_PI)
867 minPhi += -2 * M_PI;
868 if (maxPhi > 2 * M_PI)
869 maxPhi += -2 * M_PI;
870 }
871 } else {
872 minPhi = 0.;
873 maxPhi = 2 * M_PI;
874 maxR += ro * std::abs(std::cos(dir.theta()));
875 }
876 if (minPhi >= maxPhi && (minPhi - maxPhi) < M_PI) {
877 minPhi = 0.;
878 maxPhi = 2 * M_PI;
879 }
880 }
881 } // end cyl & bcyl
882
883 if (!cyl && !bcyl) {
884 // loop over edges to check inner radial extent
885 PerigeeSurface peri;
886 for (unsigned int ie = 0; ie < edges.size(); ie++) {
888 (vtxt[edges[ie].first] - vtxt[edges[ie].second]).unit();
889 Intersection closest =
890 peri.straightLineIntersection(vtxt[edges[ie].second], dir);
891 double le = (vtxt[edges[ie].first] - vtxt[edges[ie].second]).norm();
892 if ((closest.position - vtxt[edges[ie].first]).norm() < le &&
893 (closest.position - vtxt[edges[ie].second]).norm() < le)
894 if (minR > closest.position.perp())
895 minR = closest.position.perp();
896 } // end loop over edges
897 // verify phi span - may run across step
898 if (std::abs(maxPhi - minPhi) > M_PI) {
899 double phiTmp = minPhi;
900 minPhi = 3 * M_PI;
901 maxPhi = 0.; // redo the search
902 for (unsigned int iv = 0; iv < vtxt.size(); iv++) {
903 phiTmp = vtxt[iv].phi() + M_PI;
904 if (phiTmp < M_PI)
905 phiTmp += 2 * M_PI;
906 minPhi = phiTmp < minPhi ? phiTmp : minPhi;
907 maxPhi = phiTmp > maxPhi ? phiTmp : maxPhi;
908 }
909 if (minPhi > 2 * M_PI)
910 minPhi += -2 * M_PI;
911 if (maxPhi > 2 * M_PI)
912 maxPhi += -2 * M_PI;
913 if (minPhi >= maxPhi && (minPhi - maxPhi) < M_PI) {
914 minPhi = 0.;
915 maxPhi = 2 * M_PI;
916 }
917 }
918 }
919
920 if (cyl || bcyl || box || trd || dtrd || spb || prism) {
921 span.zMin = minZ - zTol;
922 span.zMax = maxZ - +zTol;
923 minPhi = (minPhi - phiTol) < 0 ? minPhi - phiTol + 2 * M_PI
924 : minPhi - phiTol;
925 span.phiMin = minPhi;
926 maxPhi = (maxPhi + phiTol) > 2 * M_PI ? maxPhi + phiTol - 2 * M_PI
927 : maxPhi + phiTol;
928 span.phiMax = maxPhi;
929 span.rMin = std::max(70.001, minR - zTol);
930 span.rMax = maxR + zTol;
931 span.xMin = minX - zTol;
932 span.xMax = maxX - +zTol;
933 span.yMin = minY - zTol;
934 span.yMax = maxY - +zTol;
935 } else {
936 ATH_MSG_WARNING("VolumeConverter::volume shape not recognized ");
937 }
938 return std::make_unique<VolumeSpan>(span);
939}
#define M_PI
Scalar perp() const
perp method - perpendicular length
const PlainObject unit() const
This is a plugin that makes Eigen look like CLHEP & defines some convenience methods.
#define ATH_MSG_WARNING(x)
std::unique_ptr< VolumeSpan > findVolumeSpan(const VolumeBounds &volBounds, const Amg::Transform3D &transform, double zTol, double phiTol) const
Estimation of the geometrical volume span.
Amg::Vector3D transform(Amg::Vector3D &v, Amg::Transform3D &tr)
Transform a point from a Trasformation3D.
Eigen::Matrix< double, 3, 1 > Vector3D
CxxUtils::transform_view_with_at< PackedLinkVector_span, ELSpanConverter > span
bool dPhi(const xAOD::TauJet &tau, const xAOD::CaloVertexedTopoCluster &cluster, float &out)
@ phi
Definition ParamDefs.h:75

◆ initMessaging()

void AthMessaging::initMessaging ( ) const
privateinherited

Initialize our message level and MessageSvc.

This method should only be called once.

Definition at line 39 of file AthMessaging.cxx.

40{
42 // If user did not set an explicit level, set a default
43 if (m_lvl == MSG::NIL) {
44 m_lvl = m_imsg ?
45 static_cast<MSG::Level>( m_imsg.load()->outputLevel(m_nm) ) :
46 MSG::INFO;
47 }
48}
std::string m_nm
Message source name.
std::atomic< IMessageSvc * > m_imsg
MessageSvc pointer.
std::atomic< MSG::Level > m_lvl
Current logging level.
IMessageSvc * getMessageSvc(bool quiet=false)

◆ leadingVolume()

double Trk::VolumeConverter::leadingVolume ( const GeoShape * sh) const
private

Definition at line 1149 of file VolumeConverter.cxx.

1149 {
1150
1151 if (sh->type() == "Subtraction") {
1152 const GeoShapeSubtraction* sub =
1153 dynamic_cast<const GeoShapeSubtraction*>(sh);
1154 if (sub)
1155 return leadingVolume(sub->getOpA());
1156 }
1157 if (sh->type() == "Union") {
1158 const GeoShapeUnion* uni = dynamic_cast<const GeoShapeUnion*>(sh);
1159 if (uni)
1160 return leadingVolume(uni->getOpA()) + leadingVolume(uni->getOpB());
1161 }
1162 if (sh->type() == "Intersection") {
1163 const GeoShapeIntersection* intr =
1164 dynamic_cast<const GeoShapeIntersection*>(sh);
1165 if (intr)
1166 return std::min(leadingVolume(intr->getOpA()),
1167 leadingVolume(intr->getOpB()));
1168 }
1169 if (sh->type() == "Shift") {
1170 const GeoShapeShift* shift = dynamic_cast<const GeoShapeShift*>(sh);
1171 if (shift)
1172 return leadingVolume(shift->getOp());
1173 }
1174
1175 return sh->volume();
1176}
double leadingVolume(const GeoShape *sh) const

◆ msg() [1/2]

MsgStream & AthMessaging::msg ( ) const
inlineinherited

The standard message stream.

Returns a reference to the default message stream May not be invoked before sysInitialize() has been invoked.

Definition at line 167 of file AthMessaging.h.

168{
169 MsgStream* ms = m_msg_tls.get();
170 if (!ms) {
171 if (!m_initialized.test_and_set()) initMessaging();
172 ms = new MsgStream(m_imsg,m_nm);
173 m_msg_tls.reset( ms );
174 }
175
176 ms->setLevel (m_lvl);
177 return *ms;
178}
boost::thread_specific_ptr< MsgStream > m_msg_tls
MsgStream instance (a std::cout like with print-out levels).
void initMessaging() const
Initialize our message level and MessageSvc.

◆ msg() [2/2]

MsgStream & AthMessaging::msg ( const MSG::Level lvl) const
inlineinherited

The standard message stream.

Returns a reference to the default message stream May not be invoked before sysInitialize() has been invoked.

Definition at line 182 of file AthMessaging.h.

183{ return msg() << lvl; }
MsgStream & msg() const
The standard message stream.

◆ msgLvl()

bool AthMessaging::msgLvl ( const MSG::Level lvl) const
inlineinherited

Test the output level.

Parameters
lvlThe message level to test against
Returns
boolean Indicating if messages at given level will be printed
Return values
trueMessages at level "lvl" will be printed

Definition at line 151 of file AthMessaging.h.

152{
153 // If user did not set explicit message level we have to initialize
154 // the messaging and retrieve the default via the MessageSvc.
155 if (m_lvl==MSG::NIL && !m_initialized.test_and_set()) initMessaging();
156
157 if (m_lvl <= lvl) {
158 msg() << lvl;
159 return true;
160 } else {
161 return false;
162 }
163}

◆ resolveBooleanVolume()

double Trk::VolumeConverter::resolveBooleanVolume ( const Volume & trVol,
double tolerance ) const

temporary owner of auxiliary volumes

small components

Definition at line 230 of file VolumeConverter.cxx.

231 {
232
233 VolumePartVec constituents{};
234 VolumePart inputVol{};
235 inputVol.parts.push_back(std::make_unique<Volume>(trVol));
236 constituents.push_back(std::move(inputVol));
237 VolumePartVec::iterator sIter = constituents.begin();
238
240 double volume = 0;
241 while (sIter != constituents.end()) {
242 bool update = false;
243 for (unsigned int ii = 0; ii < (*sIter).parts.size(); ++ii) {
244 const VolumeBounds& bounds{((*sIter).parts[ii]->volumeBounds())};
245 const CombinedVolumeBounds* comb =
246 dynamic_cast<const CombinedVolumeBounds*>(&bounds);
247 const SubtractedVolumeBounds* sub =
248 dynamic_cast<const SubtractedVolumeBounds*>(&bounds);
249 if (comb) {
250 (*sIter).parts[ii].reset(comb->first()->clone());
251 VolumePart vp(*sIter); //copy here
252 constituents.push_back(vp); //inser copy the iter can be invalidated
253 constituents.back().parts[ii].reset(comb->second()->clone()); //modify
254 constituents.push_back(vp); //push copy
255 constituents.back().parts.emplace_back(comb->second()->clone());//modify
256 constituents.back().sign = -1. * constituents.back().sign;
257 update = true;
258 break;
259 } else if (sub) {
260 (*sIter).parts[ii].reset(sub->outer()->clone());
262 double volSub = calculateVolume(*sub->inner(), true, tolerance);
263 if (volSub < tolerance) {
264 volume += -1. * (*sIter).sign * volSub;
265 } else {
266 constituents.emplace_back(*sIter);
267 constituents.back().parts.emplace_back(
268 sub->inner()->clone());
269 constituents.back().sign = -1. * constituents.back().sign;
270 }
271 update = true;
272 break;
273 } else {
274 // component small, below tolerance
275 double volSmall = calculateVolume(*(*sIter).parts[ii]);
276 if (volSmall < tolerance) {
277 sIter=constituents.erase(sIter);
278 update = true;
279 break;
280 }
281 }
282 } //
283 if (update)
284 sIter = constituents.begin();
285 else if ((*sIter).parts.size() == 1) {
286 double volSingle = calculateVolume(*(*sIter).parts[0]);
287 volume += (*sIter).sign * volSingle;
288 sIter=constituents.erase(sIter);
289 } else {
290 std::vector<std::shared_ptr<Volume>>::iterator tit =
291 (*sIter).parts.begin();
292 bool noovrlp = false;
293 while (tit + 1 != (*sIter).parts.end()) {
294 std::pair<bool, std::unique_ptr<Volume>> overlap =
295 Trk::VolumeIntersection::intersect(**tit, **(tit + 1));
296 if (overlap.first && !overlap.second) {
297 sIter=constituents.erase(sIter);
298 noovrlp = true;
299 break;
300 } // no intersection
301 else if (overlap.first && overlap.second) {
302 (*sIter).parts.erase(tit, tit + 2);
303 (*sIter).parts.push_back(std::move(overlap.second));
304 tit = (*sIter).parts.begin();
305 } else {
306 if (calculateVolume(**tit) < tolerance) {
307 sIter=constituents.erase(sIter);
308 noovrlp = true;
309 break;
310 }
311 if (calculateVolume(**(tit + 1)) < tolerance) {
312 sIter=constituents.erase(sIter);
313 noovrlp = true;
314 break;
315 }
316 std::pair<bool, std::unique_ptr<Volume>> overlap =
318 **tit, **(tit + 1));
319 if (overlap.first) {
320 if (overlap.second) {
321 (*sIter).parts.erase(tit, tit + 2);
322 (*sIter).parts.push_back(std::move(overlap.second));
323 tit = (*sIter).parts.begin();
324 } else {
325 sIter=constituents.erase(sIter);
326 noovrlp = true;
327 break; // no intersection
328 }
329 } else
330 ++tit;
331 }
332 }
333 if (noovrlp) {
334 } else if ((*sIter).parts.size() == 1) {
335 double volSingle = calculateVolume(*(*sIter).parts[0]);
336 volume += (*sIter).sign * volSingle;
337 sIter=constituents.erase(sIter);
338 } else {
339 ++sIter;
340 }
341 }
342 }
343
344 if (!constituents.empty()) {
345 ATH_MSG_VERBOSE("boolean volume resolved to "
346 << constituents.size() << " items "
347 << ":volume estimate:" << volume);
348 }
349 return volume;
350}
#define ATH_MSG_VERBOSE(x)
static std::pair< bool, std::unique_ptr< Trk::Volume > > intersectApproximative(const Volume &volA, const Volume &volB)
std::vector< VolumePart > VolumePartVec
constexpr double tolerance

◆ setLevel()

void AthMessaging::setLevel ( MSG::Level lvl)
inherited

Change the current logging level.

Use this rather than msg().setLevel() for proper operation with MT.

Definition at line 28 of file AthMessaging.cxx.

29{
30 m_lvl = lvl;
31}

◆ splitComposedVolume()

VolumeConverter::VolumePairVec Trk::VolumeConverter::splitComposedVolume ( const Volume & trVol)
static

Decomposition of volume into set of non-overlapping subtractions from analytically calculable volume.

Check whether the first operand in the iterator is a composite one

Combined one --> Union or Intersecion

Definition at line 352 of file VolumeConverter.cxx.

353 {
354
355 VolumePairVec constituents;
356 constituents.emplace_back(std::make_unique<Volume>(trVol), nullptr);
357 VolumePairVec::iterator sIter = constituents.begin();
358 std::shared_ptr<VolumeBounds> newBounds{};
359 while (sIter != constituents.end()) {
361 const CombinedVolumeBounds* comb =
362 dynamic_cast<const Trk::CombinedVolumeBounds*>(
363 &((*sIter).first->volumeBounds()));
364 const Trk::SubtractedVolumeBounds* sub =
365 dynamic_cast<const Trk::SubtractedVolumeBounds*>(
366 &((*sIter).first->volumeBounds()));
368 if (comb) {
369 std::shared_ptr<Volume> subVol = (*sIter).second;
370 sIter = constituents.erase(sIter);
371 std::shared_ptr<Volume> combFirst{comb->first()->clone()};
372 std::shared_ptr<Volume> combSecond{comb->second()->clone()};
373 if (comb->intersection()) {
374 newBounds = std::make_shared<Trk::SubtractedVolumeBounds>(
375 std::unique_ptr<Trk::Volume>(combFirst->clone()), std::unique_ptr<Trk::Volume>(combSecond->clone()));
376 std::unique_ptr<Trk::Volume> newSubVol =
377 std::make_unique<Volume>(nullptr, std::move(newBounds));
378 if (subVol) {
379 newBounds = std::make_shared<CombinedVolumeBounds>(
380 std::unique_ptr<Trk::Volume>(subVol->clone()), std::move(newSubVol), false);
381 std::shared_ptr<Volume> newCSubVol =
382 std::make_unique<Volume>(nullptr, std::move(newBounds));
383 constituents.insert(sIter,
384 std::make_pair(combFirst, newCSubVol));
385 } else {
386 constituents.insert(
387 sIter, std::make_pair(combFirst, std::move(newSubVol)));
388 }
389 } else {
390 constituents.insert(sIter, std::make_pair(combFirst, subVol));
391 if (subVol) {
392 newBounds = std::make_shared<CombinedVolumeBounds>(
393 std::unique_ptr<Trk::Volume>(subVol->clone()),
394 std::unique_ptr<Trk::Volume>(combFirst->clone()), false);
395 std::unique_ptr<Trk::Volume> newSubVol =
396 std::make_unique<Volume>(nullptr, std::move(newBounds));
397 constituents.insert(
398 sIter,
399 std::make_pair(combSecond, std::move(newSubVol)));
400 } else {
401 constituents.insert(sIter,
402 std::make_pair(combSecond, combFirst));
403 }
404 }
405 sIter = constituents.begin();
406 } else if (sub) {
407 std::shared_ptr<Volume> subVol = (*sIter).second;
408 sIter = constituents.erase(sIter);
409 std::shared_ptr<Volume> innerVol{sub->inner()->clone()};
410 std::shared_ptr<Volume> outerVol{sub->outer()->clone()};
411 if (subVol) {
412 newBounds = std::make_shared<CombinedVolumeBounds>(
413 std::unique_ptr<Trk::Volume>(subVol->clone()),
414 std::unique_ptr<Trk::Volume>(innerVol->clone()), false);
415 std::unique_ptr<Volume> newSubVol =
416 std::make_unique<Trk::Volume>(nullptr, newBounds);
417 constituents.insert(
418 sIter, std::make_pair(outerVol, std::move(newSubVol)));
419 } else {
420 constituents.insert(sIter, std::make_pair(outerVol, innerVol));
421 }
422 sIter = constituents.begin();
423 } else {
424 ++sIter;
425 }
426 }
427 return constituents;
428}
const Volume * inner() const
This method returns the inner Volume.
const Volume * outer() const
This method returns the outer Volume.
std::vector< VolumePair > VolumePairVec
virtual Volume * clone() const
polymorpic deep copy
Definition Volume.cxx:60

◆ translate()

std::unique_ptr< TrackingVolume > Trk::VolumeConverter::translate ( const GeoVPhysVol * gv,
bool simplify,
bool blend,
double blendMassLimit ) const

translation of GeoVPhysVol to Trk::TrackingVolume

Definition at line 46 of file VolumeConverter.cxx.

48 {
49
50 const std::string name = gv->getLogVol()->getName();
51
52 Amg::Transform3D ident{Amg::Transform3D::Identity()};
53 std::unique_ptr<Volume> volGeo{m_geoShapeConverter.translateGeoShape(
54 gv->getLogVol()->getShape(), ident)};
55
56 // resolve volume structure into a set of non-overlapping subtractions from
57 // analytically calculable shapes
58 VolumePairVec constituents = splitComposedVolume(*volGeo);
59
60 // material properties
61 Material mat = Trk::GeoMaterialConverter::convert(gv->getLogVol()->getMaterial());
62
63 // calculate precision of volume estimate taking into account material
64 // properties
65 double precision = s_precisionInX0 * mat.X0; // required precision in mm
66
67 // volume estimate from GeoShape
68 double volumeFromGeoShape =
69 -1; // replace with database info when available
70
71 // volume estimate from resolveBoolean
72 // double volumeBoolean = calculateVolume(volGeo,false,pow(precision,3)); //
73 // TODO : test on inert material
74
75 // volume estimate from Volume
76 double volume = volumeFromGeoShape >= 0 ? volumeFromGeoShape : -1.;
77 if ((simplify || blend) && volumeFromGeoShape < 0) {
78 double fraction = 0.;
79 volume = 0.;
80 for (const VolumePair& cs : constituents) {
81 fraction = estimateFraction(cs, precision);
82 if (fraction < 0) {
83 volume = -1;
84 break;
85 } else
86 volume += fraction * calculateVolume(*cs.first);
87 }
88 }
89
90 // evaluate complexity of shape
91 // simple case
92 if (constituents.size() == 1 && !constituents[0].second) {
93 return std::make_unique<TrackingVolume>(*volGeo, mat, nullptr, nullptr,
94 name);
95 }
96 // build envelope
97 std::unique_ptr<Volume> envelope{};
98 std::string envName = name;
99
100 std::unique_ptr<TrackingVolume> trEnv{};
101
102 bool blended = false;
103
104 if (constituents.size() == 1) {
105
106 envelope = std::make_unique<Volume>(*(constituents.front().first),
107 volGeo->transform());
108 double volEnv = calculateVolume(*constituents.front().first);
109
110 if (blend && volume > 0 && volEnv > 0 &&
111 volume * mat.rho < blendMassLimit)
112 blended = true;
113
114 if ((simplify || blended) && volume > 0 && volEnv > 0) {
115 // simplified material rescales X0, l0 and density
116 double fraction = volume / volEnv;
117 Material matScaled(mat.X0 / fraction, mat.L0 / fraction, mat.A,
118 mat.Z, fraction * mat.rho);
119 if (blend && !blended)
120 envName = envName + "_PERM";
121 trEnv = std::make_unique<TrackingVolume>(*envelope, matScaled,
122 nullptr, nullptr, envName);
123 } else {
124 auto confinedVols =std::make_unique<std::vector<TrackingVolume*>>();
125 confinedVols->push_back(std::make_unique<TrackingVolume>(*volGeo, mat, nullptr, nullptr, name).release());
126 envName = name + "_envelope";
127 trEnv = std::make_unique<TrackingVolume>(*envelope, dummyMaterial, std::move(confinedVols), envName);
128 }
129
130 return trEnv;
131 }
132
133 // composed shapes : derive envelope from span
134 Amg::Transform3D transf = volGeo->transform();
135 std::unique_ptr<VolumeSpan> span =
136 findVolumeSpan(volGeo->volumeBounds(), transf, 0., 0.);
137
138 bool isCyl = false;
139 for (const auto& fv : constituents) {
140 const CylinderVolumeBounds* cyl =
141 dynamic_cast<const CylinderVolumeBounds*>(
142 &(fv.first->volumeBounds()));
143 if (cyl) {
144 isCyl = true;
145 break;
146 }
147 }
148
149 ATH_MSG_DEBUG(__FILE__ << ":" << __LINE__
150 << "envelope estimate: object contains cylinder:"
151 << name << ":" << isCyl);
152 ATH_MSG_DEBUG(__FILE__ << ":" << __LINE__
153 << "complex volume span for envelope:" << name
154 << ":x range:" << (*span).xMin << ","
155 << (*span).xMax);
156 ATH_MSG_DEBUG(__FILE__ << ":" << __LINE__
157 << "complex volume span for envelope:" << name
158 << ":y range:" << (*span).yMin << ","
159 << (*span).yMax);
160 ATH_MSG_DEBUG(__FILE__ << ":" << __LINE__
161 << "complex volume span for envelope:" << name
162 << ":z range:" << (*span).zMin << ","
163 << (*span).zMax);
164 ATH_MSG_DEBUG(__FILE__ << ":" << __LINE__
165 << "complex volume span for envelope:" << name
166 << ":R range:" << (*span).rMin << ","
167 << (*span).rMax);
168 ATH_MSG_DEBUG(__FILE__ << ":" << __LINE__
169 << "complex volume span for envelope:" << name
170 << ":phi range:" << (*span).phiMin << ","
171 << (*span).phiMax);
172
173 if (!isCyl) { // cuboid envelope
174 Amg::Transform3D cylTrf{
175 transf * Amg::Translation3D{0.5 * ((*span).xMin + (*span).xMax),
176 0.5 * ((*span).yMin + (*span).yMax),
177 0.5 * ((*span).zMin + (*span).zMax)}};
178
179 std::shared_ptr<VolumeBounds> bounds =
180 std::make_shared<CuboidVolumeBounds>(
181 0.5 * ((*span).xMax - (*span).xMin),
182 0.5 * ((*span).yMax - (*span).yMin),
183 0.5 * ((*span).zMax - (*span).zMin));
184 envelope = std::make_unique<Volume>(
185 makeTransform(cylTrf), std::move(bounds));
186 } else {
187 double dPhi = (*span).phiMin > (*span).phiMax
188 ? (*span).phiMax - (*span).phiMin + 2 * M_PI
189 : (*span).phiMax - (*span).phiMin;
190 std::shared_ptr<VolumeBounds> cylBounds{};
191 Amg::Transform3D cylTrf{transf};
192 if (dPhi < 2 * M_PI) {
193 double aPhi = 0.5 * ((*span).phiMax + (*span).phiMin);
194 cylBounds = std::make_shared<CylinderVolumeBounds>(
195 (*span).rMin, (*span).rMax, 0.5 * dPhi,
196 0.5 * ((*span).zMax - (*span).zMin));
197 cylTrf = cylTrf * Amg::getRotateZ3D(aPhi);
198 } else {
199 cylBounds = std::make_shared<CylinderVolumeBounds>(
200 (*span).rMin, (*span).rMax,
201 0.5 * ((*span).zMax - (*span).zMin));
202 }
203 envelope = std::make_unique<Volume>(
204 makeTransform(cylTrf), std::move(cylBounds));
205 }
206
207 double volEnv = calculateVolume(*envelope);
208
209 if (blend && volume > 0 && volEnv > 0 && volume * mat.rho < blendMassLimit)
210 blended = true;
211
212 if ((simplify || blended) && volume > 0 && volEnv > 0) {
213 double fraction = volume / volEnv;
214 Material matScaled(mat.X0 / fraction, mat.L0 / fraction, mat.A, mat.Z,
215 fraction * mat.rho);
216 if (blend && !blended)
217 envName = envName + "_PERM";
218 trEnv = std::make_unique<TrackingVolume>(*envelope, mat, nullptr,
219 nullptr, envName);
220 } else {
221 auto confinedVols = std::make_unique<std::vector<TrackingVolume*>>();
222 confinedVols->push_back( std::make_unique<TrackingVolume>(*volGeo, mat, nullptr, nullptr, name).release());
223 envName = envName + "_envelope";
224 trEnv = std::make_unique<TrackingVolume>(*envelope, dummyMaterial, std::move(confinedVols), envName);
225 }
226
227 return trEnv;
228}
static VolumePairVec splitComposedVolume(const Volume &trVol)
Decomposition of volume into set of non-overlapping subtractions from analytically calculable volume.
static constexpr double s_precisionInX0
std::pair< std::shared_ptr< Volume >, std::shared_ptr< Volume > > VolumePair
double estimateFraction(const VolumePair &sub, double precision) const
the tricky part of volume calculation
static std::string release
Definition computils.h:50
Amg::Transform3D getRotateZ3D(double angle)
Rotate the coordinate system by an angle around the z-axis.
Eigen::Translation< double, 3 > Translation3D
std::unique_ptr< Amg::Transform3D > makeTransform(const Amg::Transform3D &trf)
@ ident
Definition HitInfo.h:77

Member Data Documentation

◆ ATLAS_THREAD_SAFE

std::atomic_flag m_initialized AthMessaging::ATLAS_THREAD_SAFE = ATOMIC_FLAG_INIT
mutableprivateinherited

Messaging initialized (initMessaging).

Definition at line 141 of file AthMessaging.h.

◆ m_geoShapeConverter

Trk::GeoShapeConverter Trk::VolumeConverter::m_geoShapeConverter
private

shape converter

Definition at line 102 of file VolumeConverter.h.

◆ m_imsg

std::atomic<IMessageSvc*> AthMessaging::m_imsg { nullptr }
mutableprivateinherited

MessageSvc pointer.

Definition at line 135 of file AthMessaging.h.

135{ nullptr };

◆ m_lvl

std::atomic<MSG::Level> AthMessaging::m_lvl { MSG::NIL }
mutableprivateinherited

Current logging level.

Definition at line 138 of file AthMessaging.h.

138{ MSG::NIL };

◆ m_msg_tls

boost::thread_specific_ptr<MsgStream> AthMessaging::m_msg_tls
mutableprivateinherited

MsgStream instance (a std::cout like with print-out levels).

Definition at line 132 of file AthMessaging.h.

◆ m_nm

std::string AthMessaging::m_nm
privateinherited

Message source name.

Definition at line 129 of file AthMessaging.h.

◆ s_precisionInX0

double Trk::VolumeConverter::s_precisionInX0
staticconstexprprivate
Initial value:
=
1.e-3

Definition at line 104 of file VolumeConverter.h.


The documentation for this class was generated from the following files: