ATLAS Offline Software
Public Types | Public Member Functions | Static Public Member Functions | Private Member Functions | Private Attributes | Static Private Attributes | List of all members
Trk::VolumeConverter Class Reference

#include <VolumeConverter.h>

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

Public Types

using VolumePair = std::pair< std::shared_ptr< Volume >, std::shared_ptr< Volume > >
 
using VolumePairVec = std::vector< VolumePair >
 

Public Member Functions

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

Static Public Member Functions

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

Private Member Functions

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

Private Attributes

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

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

using Trk::VolumeConverter::VolumePair = std::pair<std::shared_ptr<Volume>, std::shared_ptr<Volume> >

Definition at line 69 of file VolumeConverter.h.

◆ VolumePairVec

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

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") {}

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 }

◆ 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 }

◆ 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 }

◆ 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 }

◆ 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
773  Amg::Vector3D dir =
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++) {
817  Amg::Vector3D dir =
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++) {
883  Amg::Vector3D dir =
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 }

◆ 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 {
41  m_imsg = Athena::getMessageSvc();
42  m_lvl = m_imsg ?
43  static_cast<MSG::Level>( m_imsg.load()->outputLevel(m_nm) ) :
44  MSG::INFO;
45 }

◆ 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 }

◆ 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 164 of file AthMessaging.h.

165 {
166  MsgStream* ms = m_msg_tls.get();
167  if (!ms) {
168  if (!m_initialized.test_and_set()) initMessaging();
169  ms = new MsgStream(m_imsg,m_nm);
170  m_msg_tls.reset( ms );
171  }
172 
173  ms->setLevel (m_lvl);
174  return *ms;
175 }

◆ 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 179 of file AthMessaging.h.

180 { return msg() << lvl; }

◆ 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_initialized.test_and_set()) initMessaging();
154  if (m_lvl <= lvl) {
155  msg() << lvl;
156  return true;
157  } else {
158  return false;
159  }
160 }

◆ 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 =
313  Trk::VolumeIntersection::intersectApproximative(
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 }

◆ 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 }

◆ 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 }

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.

◆ m_lvl

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

Current logging level.

Definition at line 138 of file AthMessaging.h.

◆ 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

constexpr 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:
xAOD::iterator
JetConstituentVector::iterator iterator
Definition: JetConstituentVector.cxx:68
AthMessaging::m_lvl
std::atomic< MSG::Level > m_lvl
Current logging level.
Definition: AthMessaging.h:138
AllowedVariables::e
e
Definition: AsgElectronSelectorTool.cxx:37
CxxUtils::span
span(T *ptr, std::size_t sz) -> span< T >
A couple needed deduction guides.
ReadCellNoiseFromCoolCompare.s1
s1
Definition: ReadCellNoiseFromCoolCompare.py:378
PlotCalibFromCool.norm
norm
Definition: PlotCalibFromCool.py:100
Trk::VolumePartVec
std::vector< VolumePart > VolumePartVec
Definition: VolumeConverter.h:50
Trk::VolumeConverter::m_geoShapeConverter
Trk::GeoShapeConverter m_geoShapeConverter
shape converter
Definition: VolumeConverter.h:102
Muon::makeTransform
std::unique_ptr< Amg::Transform3D > makeTransform(const Amg::Transform3D &trf)
Definition: MuonSpectrometer/MuonDetDescr/MuonTrackingGeometry/MuonTrackingGeometry/Utils.h:14
Trk::VolumeConverter::calculateVolume
double calculateVolume(const Volume &vol, bool nonBooleanOnly=false, double precision=1.e-3) const
Volume calculation : by default return analytical solution only.
Definition: VolumeConverter.cxx:937
perp
Scalar perp() const
perp method - perpenticular length
Definition: AmgMatrixBasePlugin.h:44
Trk::VolumeConverter::VolumePairVec
std::vector< VolumePair > VolumePairVec
Definition: VolumeConverter.h:71
hist_file_dump.d
d
Definition: hist_file_dump.py:142
max
constexpr double max()
Definition: ap_fixedTest.cxx:33
Trk::VolumeConverter::s_precisionInX0
static constexpr double s_precisionInX0
Definition: VolumeConverter.h:104
mat
GeoMaterial * mat
Definition: LArDetectorConstructionTBEC.cxx:55
min
constexpr double min()
Definition: ap_fixedTest.cxx:26
geoGetVolumesNoXform
std::vector< const GeoVPhysVol * > geoGetVolumesNoXform(const GeoGraphNode *node, int depthLimit=1, int sizeHint=20)
Return the child volumes.
Definition: GeoVisitVolumes.cxx:240
Trk::VolumeConverter::findVolumeSpan
std::unique_ptr< VolumeSpan > findVolumeSpan(const VolumeBounds &volBounds, const Amg::Transform3D &transform, double zTol, double phiTol) const
Estimation of the geometrical volume span.
Definition: VolumeConverter.cxx:426
python.SystemOfUnits.second
float second
Definition: SystemOfUnits.py:135
Trk::GeoMaterialConverter::convert
static Material convert(const GeoMaterial *gm)
Single conversion , input type GeoMaterial - output type Trk::MaterialProperties.
Definition: GeoMaterialConverter.cxx:18
M_PI
#define M_PI
Definition: ActiveFraction.h:11
drawFromPickle.cos
cos
Definition: drawFromPickle.py:36
ATH_MSG_VERBOSE
#define ATH_MSG_VERBOSE(x)
Definition: AthMsgStreamMacros.h:28
AthMessaging::m_imsg
std::atomic< IMessageSvc * > m_imsg
MessageSvc pointer.
Definition: AthMessaging.h:135
Trk::VolumeIntersection::intersect
static std::pair< bool, std::unique_ptr< Trk::Volume > > intersect(const Volume &volA, const Volume &volB)
Definition: VolumeIntersection.cxx:38
python.CaloAddPedShiftConfig.type
type
Definition: CaloAddPedShiftConfig.py:42
Athena::getMessageSvc
IMessageSvc * getMessageSvc(bool quiet=false)
Definition: getMessageSvc.cxx:20
Amg::getRotateZ3D
Amg::Transform3D getRotateZ3D(double angle)
get a rotation transformation around Z-axis
Definition: GeoPrimitivesHelpers.h:270
PlotCalibFromCool.ie
ie
Definition: PlotCalibFromCool.py:420
AthMessaging::AthMessaging
AthMessaging()
Default constructor:
TileDCSDataPlotter.tit
tit
Definition: TileDCSDataPlotter.py:892
Trk::MaterialProperties::addMaterial
void addMaterial(const Material &mp, float dInX0)
Material averaging.
Definition: MaterialProperties.cxx:46
Trk::VolumeConverter::estimateFraction
double estimateFraction(const VolumePair &sub, double precision) const
the tricky part of volume calculation
Definition: VolumeConverter.cxx:1014
TrigConf::MSGTC::Level
Level
Definition: Trigger/TrigConfiguration/TrigConfBase/TrigConfBase/MsgStream.h:21
FullCPAlgorithmsTest_eljob.sh
sh
Definition: FullCPAlgorithmsTest_eljob.py:114
lumiFormat.i
int i
Definition: lumiFormat.py:85
Trk::VolumeConverter::collectMaterialContent
void collectMaterialContent(const GeoVPhysVol *gv, std::vector< Trk::MaterialComponent > &materialContent) const
material collection for volumes
Definition: VolumeConverter.cxx:1070
RCU::Shell
Definition: ShellExec.cxx:28
Trk::VolumeConverter::leadingVolume
double leadingVolume(const GeoShape *sh) const
Definition: VolumeConverter.cxx:1145
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
Amg::Transform3D
Eigen::Affine3d Transform3D
Definition: GeoPrimitives.h:46
Amg::transform
Amg::Vector3D transform(Amg::Vector3D &v, Amg::Transform3D &tr)
Transform a point from a Trasformation3D.
Definition: GeoPrimitivesHelpers.h:156
Trk::VolumeIntersection::intersectApproximative
static std::pair< bool, std::unique_ptr< Trk::Volume > > intersectApproximative(const Volume &volA, const Volume &volB)
Definition: VolumeIntersection.cxx:416
Trk::Volume::clone
virtual Volume * clone() const
polymorpic deep copy
Definition: Volume.cxx:60
drawFromPickle.tan
tan
Definition: drawFromPickle.py:36
python.SystemOfUnits.rad
float rad
Definition: SystemOfUnits.py:126
Trk::VolumeConverter::resolveBooleanVolume
double resolveBooleanVolume(const Volume &trVol, double tolerance) const
Definition: VolumeConverter.cxx:226
AthMessaging::msg
MsgStream & msg() const
The standard message stream.
Definition: AthMessaging.h:164
Trk::SubtractedVolumeBounds::outer
const Volume * outer() const
This method returns the outer Volume.
Definition: SubtractedVolumeBounds.h:111
beamspotman.dir
string dir
Definition: beamspotman.py:619
tolerance
Definition: suep_shower.h:17
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:240
TRT::Hit::ident
@ ident
Definition: HitInfo.h:77
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition: GeoPrimitives.h:47
Trk::SubtractedVolumeBounds::clone
SubtractedVolumeBounds * clone() const override final
Virtual constructor.
Definition: SubtractedVolumeBounds.h:102
Trk::GeoShapeConverter::translateGeoShape
std::unique_ptr< Volume > translateGeoShape(const GeoShape *shape, const Amg::Transform3D &trf) const
Convert an arbitrary GeoShape into Trk::Volume.
Definition: GeoShapeConverter.cxx:105
mapkey::sf
@ sf
Definition: TElectronEfficiencyCorrectionTool.cxx:38
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
python.Constants.INFO
int INFO
Definition: Control/AthenaCommon/python/Constants.py:15
unit
const PlainObject unit() const
This is a plugin that makes Eigen look like CLHEP & defines some convenience methods.
Definition: AmgMatrixBasePlugin.h:21
AthMessaging::m_nm
std::string m_nm
Message source name.
Definition: AthMessaging.h:129
Trk::SubtractedVolumeBounds
Definition: SubtractedVolumeBounds.h:40
ReadCellNoiseFromCoolCompare.s2
s2
Definition: ReadCellNoiseFromCoolCompare.py:379
Amg::Translation3D
Eigen::Translation< double, 3 > Translation3D
Definition: GeoPrimitives.h:44
Trk::VolumeConverter::VolumePair
std::pair< std::shared_ptr< Volume >, std::shared_ptr< Volume > > VolumePair
Definition: VolumeConverter.h:70
TauGNNUtils::Variables::Track::dPhi
bool dPhi(const xAOD::TauJet &tau, const xAOD::TauTrack &track, float &out)
Definition: TauGNNUtils.cxx:412
python.changerun.pv
pv
Definition: changerun.py:79
python.DecayParser.children
children
Definition: DecayParser.py:31
Trk::phi
@ phi
Definition: ParamDefs.h:75
Trk::VolumeConverter::splitComposedVolume
static VolumePairVec splitComposedVolume(const Volume &trVol)
Decomposition of volume into set of non-overlapping subtractions from analytically calculable volume.
Definition: VolumeConverter.cxx:348
AthMessaging::initMessaging
void initMessaging() const
Initialize our message level and MessageSvc.
Definition: AthMessaging.cxx:39
drawFromPickle.sin
sin
Definition: drawFromPickle.py:36
Trk::GeoMaterialConverter::dummy_material
static bool dummy_material(const GeoMaterial *)
hardcoded dummy materials : TODO : find generic criterium ( density ? radiation length ?...
Definition: GeoMaterialConverter.cxx:39
Trk::CombinedVolumeBounds
Definition: CombinedVolumeBounds.h:44
AthMessaging::m_msg_tls
boost::thread_specific_ptr< MsgStream > m_msg_tls
MsgStream instance (a std::cout like with print-out levels)
Definition: AthMessaging.h:132
pow
constexpr int pow(int base, int exp) noexcept
Definition: ap_fixedTest.cxx:15
PlotCalibFromCool.vv
vv
Definition: PlotCalibFromCool.py:716
Material
@ Material
Definition: MaterialTypes.h:8
WriteBchToCool.update
update
Definition: WriteBchToCool.py:67
Trk::SubtractedVolumeBounds::inner
const Volume * inner() const
This method returns the inner Volume.
Definition: SubtractedVolumeBounds.h:114
Trk::v
@ v
Definition: ParamDefs.h:78
python.SystemOfUnits.ms
float ms
Definition: SystemOfUnits.py:148
Generated on Mon Sep 29 2025 22:15:55 for ATLAS Offline Software by   doxygen 1.8.18