Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-01-30 09:16:59 

0001 //==========================================================================
0002 //  AIDA Detector description implementation 
0003 //--------------------------------------------------------------------------
0004 // Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
0005 // All rights reserved.
0006 //
0007 // For the licensing terms see $DD4hepINSTALL/LICENSE.
0008 // For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
0009 //
0010 // Author     : M.Frank
0011 //
0012 //==========================================================================
0013 
0014 // Framework include files
0015 #include "GeometryTreeDump.h"
0016 #include <DD4hep/Detector.h>
0017 
0018 // ROOT includes
0019 #include <TROOT.h>
0020 #include <TColor.h>
0021 #include <TGeoShape.h>
0022 #include <TGeoCone.h>
0023 #include <TGeoParaboloid.h>
0024 #include <TGeoPcon.h>
0025 #include <TGeoPgon.h>
0026 #include <TGeoSphere.h>
0027 #include <TGeoTorus.h>
0028 #include <TGeoTube.h>
0029 #include <TGeoTrd1.h>
0030 #include <TGeoTrd2.h>
0031 #include <TGeoArb8.h>
0032 #include <TGeoMatrix.h>
0033 #include <TGeoBoolNode.h>
0034 #include <TGeoCompositeShape.h>
0035 #include <TClass.h>
0036 #include <TGeoManager.h>
0037 #include <TMath.h>
0038 
0039 // C/C++ include files
0040 #include <iostream>
0041 
0042 using namespace dd4hep;
0043 
0044 namespace {
0045   std::string indent = "";
0046 
0047   /// Reference to output stream
0048   std::ostream& m_output = std::cout;
0049 
0050   void getAngles(const Double_t* m, Double_t &phi, Double_t &theta, Double_t &psi) {
0051     // Retreive Euler angles.
0052     // Check if theta is 0 or 180.
0053     if (TMath::Abs(1. - TMath::Abs(m[8])) < 1.e-9) {
0054       theta = TMath::ACos(m[8]) * RAD_2_DEGREE;
0055       phi = TMath::ATan2(-m[8] * m[1], m[0]) * RAD_2_DEGREE;
0056       psi = 0.;   // convention, phi+psi matters
0057       return;
0058     }
0059     // sin(theta) != 0
0060     phi = TMath::ATan2(m[2], -m[5]);
0061     Double_t sphi = TMath::Sin(phi);
0062     if (TMath::Abs(sphi) < 1.e-9)
0063       theta = -TMath::ASin(m[5] / TMath::Cos(phi)) * RAD_2_DEGREE;
0064     else
0065       theta = TMath::ASin(m[2] / sphi) * RAD_2_DEGREE;
0066     phi *= RAD_2_DEGREE;
0067     psi = TMath::ATan2(m[6], m[7]) * RAD_2_DEGREE;
0068   }
0069 }
0070 
0071 /// Dump logical volume in GDML format to output stream
0072 void* detail::GeometryTreeDump::handleVolume(const std::string& name, Volume vol) const {
0073   VisAttr vis = vol.visAttributes();
0074   TGeoShape* shape = vol->GetShape();
0075   TGeoMedium* medium = vol->GetMedium();
0076   int num = vol->GetNdaughters();
0077 
0078   m_output << "\t\t<volume name=\"" << name << "\">" << std::endl;
0079   m_output << "\t\t\t<solidref ref=\"" << shape->GetName() << "\"/>" << std::endl;
0080   m_output << "\t\t\t<materialref ref=\"" << medium->GetName() << "\"/>" << std::endl;
0081   if (vis.isValid()) {
0082     m_output << "\t\t\t<visref ref=\"" << vis.name() << "\"/>" << std::endl;
0083   }
0084   if (num > 0) {
0085     for (int i = 0; i < num; ++i) {
0086       TGeoNode*   geo_nod = vol.ptr()->GetNode(vol->GetNode(i)->GetName());
0087       TGeoVolume* geo_vol = geo_nod->GetVolume();
0088       TGeoMatrix* geo_mat = geo_nod->GetMatrix();
0089       m_output << "\t\t\t<physvol>" << std::endl;
0090       m_output << "\t\t\t\t<volumeref ref=\"" << geo_vol->GetName() << "\"/>" << std::endl;
0091       if (geo_mat) {
0092         if (geo_mat->IsTranslation()) {
0093           m_output << "\t\t\t\t<positionref ref=\"" << geo_nod->GetName() << "_pos\"/>" << std::endl;
0094         }
0095         if (geo_mat->IsRotation()) {
0096           m_output << "\t\t\t\t<rotationref ref=\"" << geo_nod->GetName() << "_rot\"/>" << std::endl;
0097         }
0098       }
0099       m_output << "\t\t\t</physvol>" << std::endl;
0100     }
0101   }
0102   m_output << "\t\t</volume>" << std::endl;
0103   return 0;
0104 }
0105 
0106 /// Dump solid in GDML format to output stream
0107 void* detail::GeometryTreeDump::handleSolid(const std::string& name, const TGeoShape* shape) const {
0108   if (shape) {
0109     if (shape->IsA() == TGeoBBox::Class()) {
0110       const TGeoBBox* sh = (const TGeoBBox*) shape;
0111       m_output << "\t\t<box name=\"" << name << "_shape\" x=\"" << sh->GetDX() << "\" y=\"" << sh->GetDY() << "\" z=\""
0112                << sh->GetDZ() << "\" lunit=\"cm\"/>" << std::endl;
0113     }
0114     else if (shape->IsA() == TGeoTube::Class()) {
0115       const TGeoTube* sh = (const TGeoTube*) shape;
0116       m_output << "\t\t<tube name=\"" << name << "_shape\" rmin=\"" << sh->GetRmin() << "\" rmax=\"" << sh->GetRmax() << "\" z=\""
0117                << sh->GetDz() << "\" startphi=\"0.0\" deltaphi=\"360.0\" aunit=\"deg\" lunit=\"cm\"/>" << std::endl;
0118     }
0119     else if (shape->IsA() == TGeoTubeSeg::Class()) {
0120       const TGeoTubeSeg* sh = (const TGeoTubeSeg*) shape;
0121       m_output << "\t\t<tube name=\"" << name << "_shape\" rmin=\"" << sh->GetRmin() << "\" rmax=\"" << sh->GetRmax() << "\" z=\""
0122                << sh->GetDz() << "\" startphi=\"" << sh->GetPhi1() << "\" deltaphi=\"" << sh->GetPhi2()
0123                << "\" aunit=\"deg\" lunit=\"cm\"/>" << std::endl;
0124     }
0125     else if (shape->IsA() == TGeoTrd1::Class()) {
0126       const TGeoTrd1* sh = (const TGeoTrd1*) shape;
0127       m_output << "\t\t<tube name=\"" << name << "_shape\" x1=\"" << sh->GetDx1() << "\" x2=\"" << sh->GetDx2() << "\" y1=\""
0128                << sh->GetDy() << "\" y2=\"" << sh->GetDy() << "\" z=\"" << sh->GetDz() << "\" lunit=\"cm\"/>" << std::endl;
0129     }
0130     else if (shape->IsA() == TGeoTrd2::Class()) {
0131       const TGeoTrd2* sh = (const TGeoTrd2*) shape;
0132       m_output << "\t\t<tube name=\"" << name << "_shape\" x1=\"" << sh->GetDx1() << "\" x2=\"" << sh->GetDx2() << "\" y1=\""
0133                << sh->GetDy1() << "\" y2=\"" << sh->GetDy2() << "\" z=\"" << sh->GetDz() << "\" lunit=\"cm\"/>" << std::endl;
0134     }
0135     else if (shape->IsA() == TGeoPgon::Class()) {
0136       const TGeoPgon* sh = (const TGeoPgon*) shape;
0137       m_output << "\t\t<polyhedra name=\"" << name << "_shape\" startphi=\"" << sh->GetPhi1() << "\" deltaphi=\"" << sh->GetDphi()
0138                << "\" numsides=\"" << sh->GetNedges() << "\" aunit=\"deg\" lunit=\"cm\">" << std::endl;
0139       for (int i = 0; i < sh->GetNz(); ++i) {
0140         m_output << "\t\t\t<zplane z=\"" << sh->GetZ(i) << "\" rmin=\"" << sh->GetRmin(i) << "\" rmax=\"" << sh->GetRmax(i)
0141                  << "\" lunit=\"cm\"/>" << std::endl;
0142       }
0143       m_output << "\t\t</polyhedra>" << std::endl;
0144     }
0145     else if (shape->IsA() == TGeoPcon::Class()) {
0146       const TGeoPcon* sh = (const TGeoPcon*) shape;
0147       m_output << "\t\t<polycone name=\"" << name << "_shape\" startphi=\"" << sh->GetPhi1() << "\" deltaphi=\"" << sh->GetDphi()
0148                << "\" aunit=\"deg\" lunit=\"cm\">" << std::endl;
0149       for (int i = 0; i < sh->GetNz(); ++i) {
0150         m_output << "\t\t\t<zplane z=\"" << sh->GetZ(i) << "\" rmin=\"" << sh->GetRmin(i) << "\" rmax=\"" << sh->GetRmax(i)
0151                  << "\" lunit=\"cm\"/>" << std::endl;
0152       }
0153       m_output << "\t\t</polycone>" << std::endl;
0154     }
0155     else if (shape->IsA() == TGeoCompositeShape::Class()) {
0156       std::string nn = name;
0157       const TGeoCompositeShape* sh = (const TGeoCompositeShape*) shape;
0158       const TGeoBoolNode* boolean = sh->GetBoolNode();
0159       TGeoBoolNode::EGeoBoolType oper = boolean->GetBooleanOperator();
0160 
0161       handleSolid(name + "_left", boolean->GetLeftShape());
0162       handleSolid(name + "_right", boolean->GetRightShape());
0163 
0164       if (oper == TGeoBoolNode::kGeoSubtraction)
0165         m_output << "\t\t<subtraction name=\"" << nn << "\">" << std::endl;
0166       else if (oper == TGeoBoolNode::kGeoUnion)
0167         m_output << "\t\t<union name=\"" << nn << "\">" << std::endl;
0168       else if (oper == TGeoBoolNode::kGeoIntersection)
0169         m_output << "\t\t<intersection name=\"" << nn << "\">" << std::endl;
0170 
0171       m_output << "\t\t\t<first ref=\"" << nn << "_left\"/>" << std::endl;
0172       m_output << "\t\t\t<second ref=\"" << nn << "_right\"/>" << std::endl;
0173       indent = "\t";
0174       handleTransformation("", boolean->GetRightMatrix());
0175       indent = "";
0176 
0177       if (oper == TGeoBoolNode::kGeoSubtraction)
0178         m_output << "\t\t</subtraction>" << std::endl;
0179       else if (oper == TGeoBoolNode::kGeoUnion)
0180         m_output << "\t\t</union>" << std::endl;
0181       else if (oper == TGeoBoolNode::kGeoIntersection)
0182         m_output << "\t\t</intersection>" << std::endl;
0183     }
0184     else {
0185       std::cerr << "Failed to handle unknwon solid shape:" << shape->IsA()->GetName() << std::endl;
0186     }
0187   }
0188   return 0;
0189 }
0190 
0191 /// Dump structure information in GDML format to output stream
0192 void detail::GeometryTreeDump::handleStructure(const std::set<Volume>& volset) const {
0193   m_output << "\t<structure>" << std::endl;
0194   for (const auto& vol : volset)
0195     handleVolume(vol->GetName(), vol);
0196   m_output << "\t</structure>" << std::endl;
0197 }
0198 
0199 /// Dump single volume transformation in GDML format to output stream
0200 void* detail::GeometryTreeDump::handleTransformation(const std::string& name, const TGeoMatrix* mat) const {
0201   if (mat) {
0202     if (mat->IsTranslation()) {
0203       const Double_t* f = mat->GetTranslation();
0204       m_output << indent << "\t\t<position ";
0205       if (!name.empty())
0206         m_output << "name=\"" << name << "_pos\" ";
0207       m_output << "x=\"" << f[0] << "\" " << "y=\"" << f[1] << "\" " << "z=\"" << f[2] << "\" unit=\"cm\"/>" << std::endl;
0208     }
0209     if (mat->IsRotation()) {
0210       const Double_t* matrix = mat->GetRotationMatrix();
0211       Double_t theta = 0.0, phi = 0.0, psi = 0.0;
0212       getAngles(matrix, theta, phi, psi);
0213       m_output << indent << "\t\t<rotation ";
0214       if (!name.empty())
0215         m_output << "name=\"" << name << "_rot\" ";
0216       m_output << "x=\"" << theta << "\" " << "y=\"" << psi << "\" " << "z=\"" << phi << "\" unit=\"deg\"/>" << std::endl;
0217     }
0218   }
0219   return 0;
0220 }
0221 
0222 /// Dump Transformations in GDML format to output stream
0223 void detail::GeometryTreeDump::handleTransformations(const std::vector<std::pair<std::string, TGeoMatrix*> >& trafos) const {
0224   m_output << "\t<define>" << std::endl;
0225   for ( const auto& t : trafos )
0226     handleTransformation(t.first, t.second);
0227   m_output << "\t</define>" << std::endl;
0228 }
0229 
0230 /// Dump all solids in GDML format to output stream
0231 void detail::GeometryTreeDump::handleSolids(const std::set<TGeoShape*>& solids) const {
0232   m_output << "\t<solids>" << std::endl;
0233   for ( const auto sh : solids )
0234     handleSolid(sh->GetName(), sh);
0235   m_output << "\t</solids>" << std::endl;
0236 }
0237 
0238 /// Dump all constants in GDML format to output stream
0239 void detail::GeometryTreeDump::handleDefines(const Detector::HandleMap& defs) const {
0240   m_output << "\t<define>" << std::endl;
0241   for ( const auto& def : defs )
0242     m_output << "\t\t<constant name=\"" << def.second->name << "\" value=\"" << def.second->type << "\" />"
0243              << std::endl;
0244   m_output << "\t</define>" << std::endl;
0245 }
0246 
0247 /// Dump all visualisation specs in Detector format to output stream
0248 void detail::GeometryTreeDump::handleVisualisation(const Detector::HandleMap&) const {
0249 }
0250 
0251 static std::string _path = "";
0252 static void dumpStructure(PlacedVolume pv, int level) {
0253   TGeoNode* current = pv.ptr();
0254   TGeoVolume* volume = current->GetVolume();
0255   TObjArray* nodes = volume->GetNodes();
0256   int num_children = nodes ? nodes->GetEntriesFast() : 0;
0257   char fmt[128];
0258 
0259   _path += "/";
0260   _path += current->GetName();
0261   ::snprintf(fmt, sizeof(fmt), "%%4d %%%ds %%7s %%s\n", level * 2 + 5);
0262   ::printf(fmt, level, "", "  ->LV:  ", volume->GetName());
0263   ::printf(fmt, level, "", "  ->PV:  ", current->GetName());
0264   ::printf(fmt, level, "", "  ->path:", _path.c_str());
0265   if (num_children > 0) {
0266     for (int i = 0; i < num_children; ++i) {
0267       TGeoNode* node = static_cast<TGeoNode*>(nodes->At(i));
0268       dumpStructure(PlacedVolume(node), level + 1);
0269     }
0270   }
0271   _path = _path.substr(0, _path.length() - 1 - strlen(current->GetName()));
0272 }
0273 
0274 static void dumpDetectors(DetElement parent, int level) {
0275   const DetElement::Children& children = parent.children();
0276   PlacedVolume pl = parent.placement();
0277   char fmt[128];
0278 
0279   _path += "/";
0280   _path += parent.name();
0281   ::snprintf(fmt, sizeof(fmt), "%%4d %%%ds %%7s %%s\n", level * 2 + 5);
0282   ::printf(fmt, level, "", "->path:", _path.c_str());
0283   if (pl.isValid()) {
0284     ::printf(fmt, level, "", "   ->placement:", parent.placementPath().c_str());
0285     ::printf(fmt, level, "", "   ->logvol:   ", pl->GetVolume()->GetName());
0286     ::printf(fmt, level, "", "   ->shape:    ", pl->GetVolume()->GetShape()->GetName());
0287   }
0288   for (const auto& c : children)
0289     dumpDetectors(c.second, level + 1);
0290 
0291   _path = _path.substr(0, _path.length() - 1 - strlen(parent.name()));
0292 }
0293 
0294 void detail::GeometryTreeDump::create(DetElement top) {
0295   dumpDetectors(top, 0);
0296   dumpStructure(top.placement(), 0);
0297 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!