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

928  {
929 
930  double volume = -1.;
931 
932  const CylinderVolumeBounds* cyl = dynamic_cast<const CylinderVolumeBounds*>(&(vol.volumeBounds()));
933  const CuboidVolumeBounds* box = dynamic_cast<const CuboidVolumeBounds*>(&(vol.volumeBounds()));
934  const TrapezoidVolumeBounds* trd = dynamic_cast<const TrapezoidVolumeBounds*>(&(vol.volumeBounds()));
935  const BevelledCylinderVolumeBounds* bcyl =dynamic_cast<const BevelledCylinderVolumeBounds*>(&(vol.volumeBounds()));
936  const PrismVolumeBounds* prism =dynamic_cast<const PrismVolumeBounds*>(&(vol.volumeBounds()));
937  const SimplePolygonBrepVolumeBounds* spb = dynamic_cast<const SimplePolygonBrepVolumeBounds*>(&(vol.volumeBounds()));
938  const CombinedVolumeBounds* comb =dynamic_cast<const CombinedVolumeBounds*>(&(vol.volumeBounds()));
939  const SubtractedVolumeBounds* sub = dynamic_cast<const SubtractedVolumeBounds*>(&(vol.volumeBounds()));
940 
941  if (cyl) {
942  return 2 * cyl->halfPhiSector() * cyl->halflengthZ() *
943  (std::pow(cyl->outerRadius(), 2) -
944  std::pow(cyl->innerRadius(), 2));
945  }
946  if (box) {
947  return 8 * box->halflengthX() * box->halflengthY() * box->halflengthZ();
948  }
949  if (trd) {
950  return 4 * (trd->minHalflengthX() + trd->maxHalflengthX()) *
951  trd->halflengthY() * trd->halflengthZ();
952  }
953  if (bcyl) {
954  int type = bcyl->type();
955  if (type < 1)
956  return 2 * bcyl->halfPhiSector() *
957  (std::pow(bcyl->outerRadius(), 2) -
958  std::pow(bcyl->innerRadius(), 2)) *
959  bcyl->halflengthZ();
960  if (type == 1)
961  return 2 * bcyl->halflengthZ() *
962  (bcyl->halfPhiSector() * std::pow(bcyl->outerRadius(), 2) -
963  std::pow(bcyl->innerRadius(), 2) *
964  std::tan(bcyl->halfPhiSector()));
965  if (type == 2)
966  return 2 * bcyl->halflengthZ() *
967  (-bcyl->halfPhiSector() * std::pow(bcyl->innerRadius(), 2) +
968  std::pow(bcyl->outerRadius(), 2) *
969  std::tan(bcyl->halfPhiSector()));
970  if (type == 3)
971  return 2 * bcyl->halflengthZ() * std::tan(bcyl->halfPhiSector()) *
972  (std::pow(bcyl->outerRadius(), 2) -
973  std::pow(bcyl->innerRadius(), 2));
974  }
975  if (prism) {
976 
977  std::vector<std::pair<double, double>> v = prism->xyVertices();
978  double vv = v[0].first * (v[1].second - v.back().second);
979  for (unsigned int i = 1; i < v.size() - 1; i++) {
980  vv += v[i].first * (v[i + 1].second - v[i - 1].second);
981  }
982  vv += v.back().first * (v[0].second - v[v.size() - 2].second);
983  return vv * prism->halflengthZ();
984  }
985  if (spb) {
986  std::vector<std::pair<double, double>> v = spb->xyVertices();
987  double vv = v[0].first * (v[1].second - v.back().second);
988  for (unsigned int i = 1; i < v.size() - 1; i++) {
989  vv += v[i].first * (v[i + 1].second - v[i - 1].second);
990  }
991  vv += v.back().first * (v[0].second - v[v.size() - 2].second);
992  return vv * spb->halflengthZ();
993  }
994 
995  if (nonBooleanOnly)
996  return volume;
997 
998  if (comb || sub) {
999  return resolveBooleanVolume(vol, precision);
1000  }
1001  return volume;
1002 }

◆ collectMaterial()

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

material collection for layers

Definition at line 1035 of file VolumeConverter.cxx.

1037  {
1038  // sf is the area of the layer collecting the material
1039 
1040  // solution relying on GeoModel
1041  // currently involves hit&miss on-fly calculation of boolean volumes
1042  // GeoModelTools::MaterialComponent mat =
1043  // gm_materialHelper.collectMaterial(pv); Material newMP =
1044  // convert(mat.first); double d = mat.second / sf; layMat.addMaterial(newMP,
1045  // d / newMP.x0()); return;
1046 
1047  std::vector<MaterialComponent> materialContent;
1048  collectMaterialContent(pv, materialContent);
1049 
1050  for (const auto& mat : materialContent) {
1051  if (mat.second < 0)
1052  continue; // protection unsolved booleans
1053  double d = sf > 0 ? mat.second / sf : 0.;
1054  if (d > 0)
1055  layMat.addMaterial(mat.first,
1056  (mat.first.X0 > 0 ? d / mat.first.X0 : 0.));
1057  }
1058 }

◆ collectMaterialContent()

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

material collection for volumes

Definition at line 1060 of file VolumeConverter.cxx.

1062  {
1063 
1064  // solution relying on GeoModel
1065  // currently involves hit&miss on-fly calculation of boolean volumes
1066  // GeoModelTools::MaterialComponent mat =
1067  // gm_materialHelper.collectMaterial(pv); Material newMP =
1068  // convert(mat.first); materialContent.push_back( MaterialComponent( newMP,
1069  // mat.second) ); return;
1070 
1071  const GeoLogVol* lv = gv->getLogVol();
1072  Material mat = Trk::GeoMaterialConverter::convert(lv->getMaterial());
1073 
1074  double motherVolume = 0.;
1075 
1076  // skip volume calculation for dummy material configuration
1077  if (!Trk::GeoMaterialConverter::dummy_material(lv->getMaterial())) {
1078  const GeoShape* sh = lv->getShape();
1079  while (sh && sh->type() == "Shift") {
1080  const GeoShapeShift* shift = dynamic_cast<const GeoShapeShift*>(sh);
1081  sh = shift ? shift->getOp() : nullptr;
1082  }
1083 
1084  bool isBoolean =
1085  sh && (sh->type() == "Subtraction" || sh->type() == "Union" ||
1086  sh->type() == "Intersection");
1087 
1088  if (isBoolean) {
1089  Amg::Transform3D transf{Amg::Transform3D::Identity()};
1090  std::unique_ptr<Volume> vol{
1091  m_geoShapeConverter.translateGeoShape(sh, transf)};
1092  motherVolume =
1093  calculateVolume(*vol, false, std::pow(1.e-3 * mat.X0, 3));
1094  if (motherVolume < 0) {
1095  // m_geoShapeConverter.decodeShape(sh);
1096  }
1097  } else
1098  motherVolume = lv->getShape()->volume();
1099  }
1100 
1101  double childVol = 0;
1102  std::string cPrevious = " ";
1103  size_t nIdentical = 0;
1104  std::vector<Trk::MaterialComponent> childMat;
1105  std::vector<const GeoVPhysVol*> children = geoGetVolumesNoXform (gv);
1106  for (const GeoVPhysVol* cv : children) {
1107  std::string cname = cv->getLogVol()->getName();
1108  if (cname == cPrevious)
1109  nIdentical++; // assuming identity for identical name and branching
1110  // history
1111  else { // scale and collect material from previous item
1112  for (const auto& cmat : childMat) {
1113  materialContent.emplace_back(cmat.first, nIdentical * cmat.second);
1114  childVol += materialContent.back().second;
1115  }
1116  childMat.clear(); // reset
1117  nIdentical = 1; // current
1118  collectMaterialContent(cv, childMat);
1119  }
1120  }
1121  for (const auto& cmat : childMat) {
1122  materialContent.emplace_back(cmat.first, nIdentical * cmat.second);
1123  childVol += materialContent.back().second;
1124  }
1125  if (motherVolume > 0 && childVol > 0)
1126  motherVolume += -1. * childVol;
1127 
1128  ATH_MSG_DEBUG("collected material:" << lv->getName() << ":made of:"
1129  << lv->getMaterial()->getName()
1130  << ":density(g/mm3)" << mat.rho
1131  << ":mass:" << mat.rho * motherVolume);
1132  materialContent.emplace_back(mat, motherVolume);
1133 }

◆ estimateFraction()

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

the tricky part of volume calculation

Definition at line 1004 of file VolumeConverter.cxx.

1005  {
1006 
1007  if (!sub.first)
1008  return 0.;
1009 
1010  if (sub.first && !sub.second)
1011  return 1.;
1012  double fraction = -1.;
1013 
1014  std::pair<bool, std::unique_ptr<Volume>> overlap =
1015  Trk::VolumeIntersection::intersect(*sub.first, *sub.second);
1016 
1017  if (overlap.first && !overlap.second)
1018  return fraction = 1.;
1019  else if (overlap.first && overlap.second) {
1020  fraction = 1. - calculateVolume(*overlap.second, true, precision) /
1021  calculateVolume(*sub.first, true, precision);
1022  return fraction;
1023  }
1024  // resolve embedded volumes
1025 
1026  // trivial within required precision
1027  double volA = calculateVolume(*sub.first, true, precision);
1028  double volB = calculateVolume(*sub.second, true, precision);
1029  if ((volA > 0 && volA < precision) || (volB > 0 && volB < precision))
1030  return 1.;
1031 
1032  return fraction;
1033 }

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

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

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

1135  {
1136 
1137  if (sh->type() == "Subtraction") {
1138  const GeoShapeSubtraction* sub =
1139  dynamic_cast<const GeoShapeSubtraction*>(sh);
1140  if (sub)
1141  return leadingVolume(sub->getOpA());
1142  }
1143  if (sh->type() == "Union") {
1144  const GeoShapeUnion* uni = dynamic_cast<const GeoShapeUnion*>(sh);
1145  if (uni)
1146  return leadingVolume(uni->getOpA()) + leadingVolume(uni->getOpB());
1147  }
1148  if (sh->type() == "Intersection") {
1149  const GeoShapeIntersection* intr =
1150  dynamic_cast<const GeoShapeIntersection*>(sh);
1151  if (intr)
1152  return std::min(leadingVolume(intr->getOpA()),
1153  leadingVolume(intr->getOpB()));
1154  }
1155  if (sh->type() == "Shift") {
1156  const GeoShapeShift* shift = dynamic_cast<const GeoShapeShift*>(sh);
1157  if (shift)
1158  return leadingVolume(shift->getOp());
1159  }
1160 
1161  return sh->volume();
1162 }

◆ 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);
248  constituents.push_back(vp);
249  constituents.back().parts[ii].reset(comb->second()->clone());
250  constituents.push_back(vp);
251  constituents.back().parts.emplace_back(comb->second()->clone());
252  constituents.back().sign = -1. * constituents.back().sign;
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::unique_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_unique<Trk::SubtractedVolumeBounds>(
371  combFirst->clone(), combSecond->clone());
372  std::unique_ptr<Trk::Volume> newSubVol =
373  std::make_unique<Volume>(nullptr, newBounds.release());
374  if (subVol) {
375  newBounds = std::make_unique<CombinedVolumeBounds>(
376  subVol->clone(), newSubVol.release(), false);
377  std::shared_ptr<Volume> newCSubVol =
378  std::make_unique<Volume>(nullptr, newBounds.release());
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_unique<CombinedVolumeBounds>(
389  subVol->clone(), combFirst->clone(), false);
390  std::unique_ptr<Trk::Volume> newSubVol =
391  std::make_unique<Volume>(nullptr, newBounds.release());
392  constituents.insert(
393  sIter,
394  std::make_pair(combSecond, std::move(newSubVol)));
395  } else {
396  constituents.insert(sIter,
397  std::make_pair(combSecond, combFirst));
398  }
399  }
400  sIter = constituents.begin();
401  } else if (sub) {
402  std::shared_ptr<Volume> subVol = (*sIter).second;
403  sIter = constituents.erase(sIter);
404  std::shared_ptr<Volume> innerVol{sub->inner()->clone()};
405  std::shared_ptr<Volume> outerVol{sub->outer()->clone()};
406  if (subVol) {
407  newBounds = std::make_unique<CombinedVolumeBounds>(
408  subVol->clone(), innerVol->clone(), false);
409  std::unique_ptr<Volume> newSubVol =
410  std::make_unique<Trk::Volume>(nullptr, newBounds.release());
411  constituents.insert(
412  sIter, std::make_pair(outerVol, std::move(newSubVol)));
413  } else {
414  constituents.insert(sIter, std::make_pair(outerVol, innerVol));
415  }
416  sIter = constituents.begin();
417  } else {
418  ++sIter;
419  }
420  }
421  return constituents;
422 }

◆ 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, confinedVols.release(), 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::unique_ptr<VolumeBounds> bounds =
176  std::make_unique<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), bounds.release());
182  } else {
183  double dPhi = (*span).phiMin > (*span).phiMax
184  ? (*span).phiMax - (*span).phiMin + 2 * M_PI
185  : (*span).phiMax - (*span).phiMin;
186  std::unique_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_unique<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_unique<CylinderVolumeBounds>(
196  (*span).rMin, (*span).rMax,
197  0.5 * ((*span).zMax - (*span).zMin));
198  }
199  envelope = std::make_unique<Volume>(
200  makeTransform(cylTrf), cylBounds.release());
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, confinedVols.release(), 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
Muon::makeTransform
Amg::Transform3D * makeTransform(const Amg::Transform3D &trf)
Definition: MuonSpectrometer/MuonDetDescr/MuonTrackingGeometry/MuonTrackingGeometry/Utils.h:14
python.SystemOfUnits.second
int second
Definition: SystemOfUnits.py:120
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
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:927
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:137
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:424
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
python.SystemOfUnits.ms
int ms
Definition: SystemOfUnits.py:132
Trk::VolumeIntersection::intersect
static std::pair< bool, std::unique_ptr< Trk::Volume > > intersect(const Volume &volA, const Volume &volB)
Definition: VolumeIntersection.cxx:38
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:890
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:1004
TrigConf::MSGTC::Level
Level
Definition: Trigger/TrigConfiguration/TrigConfBase/TrigConfBase/MsgStream.h:21
FullCPAlgorithmsTest_eljob.sh
sh
Definition: FullCPAlgorithmsTest_eljob.py:111
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:1060
RCU::Shell
Definition: ShellExec.cxx:28
Trk::VolumeConverter::leadingVolume
double leadingVolume(const GeoShape *sh) const
Definition: VolumeConverter.cxx:1135
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
TauGNNUtils::Variables::Track::dPhi
bool dPhi(const xAOD::TauJet &tau, const xAOD::TauTrack &track, double &out)
Definition: TauGNNUtils.cxx:538
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
Pseudo-constructor.
Definition: Volume.cxx:78
drawFromPickle.tan
tan
Definition: drawFromPickle.py:36
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:113
beamspotman.dir
string dir
Definition: beamspotman.py:623
tolerance
Definition: suep_shower.h:17
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:228
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:104
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:102
mapkey::sf
@ sf
Definition: TElectronEfficiencyCorrectionTool.cxx:38
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
unit
const PlainObject unit() const
This is a plugin that makes Eigen look like CLHEP & defines some convenience methods.
Definition: AmgMatrixBasePlugin.h:21
python.CaloScaleNoiseConfig.type
type
Definition: CaloScaleNoiseConfig.py:78
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
python.changerun.pv
pv
Definition: changerun.py:81
python.DecayParser.children
children
Definition: DecayParser.py:32
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:42
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:116
python.SystemOfUnits.rad
int rad
Definition: SystemOfUnits.py:111
Trk::v
@ v
Definition: ParamDefs.h:78
Generated on Sun Dec 22 2024 22:10:37 for ATLAS Offline Software by   doxygen 1.8.18