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 40 of file VolumeConverter.cxx.

40: 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 937 of file VolumeConverter.cxx.

938 {
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}
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 1045 of file VolumeConverter.cxx.

1047 {
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}
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 1070 of file VolumeConverter.cxx.

1072 {
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{
1101 m_geoShapeConverter.translateGeoShape(sh, transf)};
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}
#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

◆ estimateFraction()

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

the tricky part of volume calculation

Definition at line 1014 of file VolumeConverter.cxx.

1015 {
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}
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 426 of file VolumeConverter.cxx.

428 {
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;
521 VolumeSpan span;
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}
#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
span(T *ptr, std::size_t sz) -> span< T >
A couple needed deduction guides.
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 1145 of file VolumeConverter.cxx.

1145 {
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}
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 163 of file AthMessaging.h.

164{
165 MsgStream* ms = m_msg_tls.get();
166 if (!ms) {
167 if (!m_initialized.test_and_set()) initMessaging();
168 ms = new MsgStream(m_imsg,m_nm);
169 m_msg_tls.reset( ms );
170 }
171
172 ms->setLevel (m_lvl);
173 return *ms;
174}
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 178 of file AthMessaging.h.

179{ 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 (m_lvl <= lvl) {
154 msg() << lvl;
155 return true;
156 } else {
157 return false;
158 }
159}

◆ resolveBooleanVolume()

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

temporary owner of auxiliary volumes

small components

Definition at line 226 of file VolumeConverter.cxx.

227 {
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 =
291 Trk::VolumeIntersection::intersect(**tit, **(tit + 1));
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}
#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 348 of file VolumeConverter.cxx.

349 {
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}
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 42 of file VolumeConverter.cxx.

44 {
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}
static VolumePairVec splitComposedVolume(const Volume &trVol)
Decomposition of volume into set of non-overlapping subtractions from analytically calculable volume.
std::pair< std::shared_ptr< Volume >, std::shared_ptr< Volume > > VolumePair
static constexpr double s_precisionInX0
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)
get a rotation transformation around 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: