DDG4/src/Geant4ShapeConverter.cpp

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 <DD4hep/Shapes.h>
0016 #include <DD4hep/Printout.h>
0017 #include <DD4hep/detail/ShapesInterna.h>
0018
0019 #include "Geant4ShapeConverter.h"
0020
0021 // ROOT includes
0022 #include <TClass.h>
0023 #include <TGeoMatrix.h>
0024 #include <TGeoBoolNode.h>
0025 #include <TGeoScaledShape.h>
0026
0027 // Geant4 include files
0028 #include <G4Box.hh>
0029 #include <G4Trd.hh>
0030 #include <G4Tubs.hh>
0031 #include <G4Trap.hh>
0032 #include <G4Cons.hh>
0033 #include <G4Hype.hh>
0034 #include <G4Para.hh>
0035 #include <G4Torus.hh>
0036 #include <G4Sphere.hh>
0037 #include <G4CutTubs.hh>
0038 #include <G4Polycone.hh>
0039 #include <G4Polyhedra.hh>
0040 #include <G4Ellipsoid.hh>
0041 #include <G4Paraboloid.hh>
0042 #include <G4TwistedTubs.hh>
0043 #include <G4GenericTrap.hh>
0044 #include <G4ExtrudedSolid.hh>
0045 #include <G4EllipticalTube.hh>
0046 #include <G4TessellatedSolid.hh>
0047 #include <G4TriangularFacet.hh>
0048 #include <G4QuadrangularFacet.hh>
0049
0050 // C/C++ include files
0051
0052 using namespace dd4hep::detail;
0053
0054 /// Namespace for the AIDA detector description toolkit
0055 namespace dd4hep {
0056
0057   /// Namespace for the Geant4 based simulation part of the AIDA detector description toolkit
0058   namespace sim {
0059
0060     static const double CM_2_MM = (CLHEP::centimeter/dd4hep::centimeter);
0061
0062     /// Convert a specific TGeo shape into the geant4 equivalent
0063     template <typename T> G4VSolid* convertShape(const TGeoShape* shape)    {
0064       if ( shape )   {
0065         except("convertShape","Unsupported shape: %s",shape->IsA()->GetName());
0066       }
0067       except("convertShape","Invalid shape conversion requested.");
0068       return 0;
0069     }
0070
0071     template <> G4VSolid* convertShape<TGeoShapeAssembly>(const TGeoShape* /* shape */)  {
0072       return 0;
0073     }
0074
0075     template <> G4VSolid* convertShape<TGeoBBox>(const TGeoShape* shape)  {
0076       const TGeoBBox* sh = (const TGeoBBox*) shape;
0077       return new G4Box(sh->GetName(), sh->GetDX() * CM_2_MM, sh->GetDY() * CM_2_MM, sh->GetDZ() * CM_2_MM);
0078     }
0079
0080     template <> G4VSolid* convertShape<TGeoTube>(const TGeoShape* shape)  {
0081       const TGeoTube* sh = (const TGeoTube*) shape;
0082       return new G4Tubs(sh->GetName(), sh->GetRmin() * CM_2_MM, sh->GetRmax() * CM_2_MM, sh->GetDz() * CM_2_MM, 0, 2. * M_PI);
0083     }
0084
0085     template <> G4VSolid* convertShape<TGeoTubeSeg>(const TGeoShape* shape)  {
0086       const TGeoTubeSeg* sh = (const TGeoTubeSeg*) shape;
0087       return new G4Tubs(sh->GetName(), sh->GetRmin() * CM_2_MM, sh->GetRmax() * CM_2_MM, sh->GetDz() * CM_2_MM,
0088                         sh->GetPhi1() * DEGREE_2_RAD, (sh->GetPhi2()-sh->GetPhi1()) * DEGREE_2_RAD);
0089     }
0090
0091     template <> G4VSolid* convertShape<TGeoCtub>(const TGeoShape* shape)  {
0092       const TGeoCtub* sh = (const TGeoCtub*) shape;
0093       const Double_t* ln = sh->GetNlow();
0094       const Double_t* hn = sh->GetNhigh();
0095       G4ThreeVector   lowNorm (ln[0], ln[1], ln[2]);
0096       G4ThreeVector   highNorm(hn[0], hn[1], hn[2]);
0097       return new G4CutTubs(sh->GetName(),
0098                            sh->GetRmin() * CM_2_MM, sh->GetRmax() * CM_2_MM, sh->GetDz() * CM_2_MM,
0099                            sh->GetPhi1() * DEGREE_2_RAD, (sh->GetPhi2()-sh->GetPhi1()) * DEGREE_2_RAD,
0100                            std::move(lowNorm), std::move(highNorm));
0101     }
0102
0103     template <> G4VSolid* convertShape<TGeoEltu>(const TGeoShape* shape)  {
0104       const TGeoEltu* sh = (const TGeoEltu*) shape;
0105       return new G4EllipticalTube(sh->GetName(),sh->GetA() * CM_2_MM, sh->GetB() * CM_2_MM, sh->GetDz() * CM_2_MM);
0106     }
0107
0108     template <> G4VSolid* convertShape<TwistedTubeObject>(const TGeoShape* shape)  {
0109       const TwistedTubeObject* sh = (const TwistedTubeObject*) shape;
0110       if ( std::fabs(std::fabs(sh->GetNegativeEndZ()) - std::fabs(sh->GetPositiveEndZ())) < 1e-10 )   {
0111         return new G4TwistedTubs(sh->GetName(),sh->GetPhiTwist() * DEGREE_2_RAD,
0112                                  sh->GetRmin() * CM_2_MM, sh->GetRmax() * CM_2_MM,
0113                                  sh->GetPositiveEndZ() * CM_2_MM,
0114                                  sh->GetNsegments(),
0115                                  (sh->GetPhi2()-sh->GetPhi1()) * DEGREE_2_RAD);
0116       }
0117       return new G4TwistedTubs(sh->GetName(),sh->GetPhiTwist() * DEGREE_2_RAD,
0118                                sh->GetRmin() * CM_2_MM, sh->GetRmax() * CM_2_MM,
0119                                sh->GetNegativeEndZ() * CM_2_MM, sh->GetPositiveEndZ() * CM_2_MM,
0120                                sh->GetNsegments(),
0121                                (sh->GetPhi2()-sh->GetPhi1()) * DEGREE_2_RAD);
0122     }
0123
0124     template <> G4VSolid* convertShape<TGeoTrd1>(const TGeoShape* shape)  {
0125       const TGeoTrd1* sh = (const TGeoTrd1*) shape;
0126       return new G4Trd(sh->GetName(),
0127                        sh->GetDx1() * CM_2_MM, sh->GetDx2() * CM_2_MM,
0128                        sh->GetDy() * CM_2_MM, sh->GetDy() * CM_2_MM,
0129                        sh->GetDz() * CM_2_MM);
0130     }
0131
0132     template <> G4VSolid* convertShape<TGeoTrd2>(const TGeoShape* shape)  {
0133       const TGeoTrd2* sh = (const TGeoTrd2*) shape;
0134       return new G4Trd(sh->GetName(),
0135                        sh->GetDx1() * CM_2_MM, sh->GetDx2() * CM_2_MM,
0136                        sh->GetDy1() * CM_2_MM, sh->GetDy2() * CM_2_MM,
0137                        sh->GetDz() * CM_2_MM);
0138     }
0139
0140     template <> G4VSolid* convertShape<TGeoHype>(const TGeoShape* shape)  {
0141       const TGeoHype* sh = (const TGeoHype*) shape;
0142       return new G4Hype(sh->GetName(), sh->GetRmin() * CM_2_MM, sh->GetRmax() * CM_2_MM,
0143                         sh->GetStIn() * DEGREE_2_RAD, sh->GetStOut() * DEGREE_2_RAD,
0144                         sh->GetDz() * CM_2_MM);
0145     }
0146
0147     template <> G4VSolid* convertShape<TGeoArb8>(const TGeoShape* shape)  {
0148       std::vector<G4TwoVector> vertices;
0149       TGeoArb8* sh = (TGeoArb8*) shape;
0150       Double_t* vtx_xy = sh->GetVertices();
0151       for ( std::size_t i=0; i<8; ++i, vtx_xy +=2 )
0152         vertices.emplace_back(vtx_xy[0] * CM_2_MM, vtx_xy[1] * CM_2_MM);
0153       return new G4GenericTrap(sh->GetName(), sh->GetDz() * CM_2_MM, vertices);
0154     }
0155
0156     template <> G4VSolid* convertShape<TGeoPara>(const TGeoShape* shape) {
0157       const auto* sh = static_cast<const TGeoPara*>(shape);
0158       return new G4Para(sh->GetName(),
0159                         sh->GetX() * CM_2_MM, sh->GetY() * CM_2_MM, sh->GetZ() * CM_2_MM,
0160                         sh->GetAlpha() * DEGREE_2_RAD, sh->GetTheta() * DEGREE_2_RAD,
0161                         sh->GetPhi() * DEGREE_2_RAD);
0162     }
0163
0164     template <> G4VSolid* convertShape<TGeoXtru>(const TGeoShape* shape)  {
0165       const TGeoXtru* sh = (const TGeoXtru*) shape;
0166       std::size_t nz = sh->GetNz();
0167       std::vector<G4ExtrudedSolid::ZSection> z;
0168       z.reserve(nz);
0169       for(std::size_t i=0; i<nz; ++i)   {
0170         z.emplace_back(G4ExtrudedSolid::ZSection(sh->GetZ(i) * CM_2_MM, {sh->GetXOffset(i), sh->GetYOffset(i)}, sh->GetScale(i)));
0171       }
0172       std::size_t np = sh->GetNvert();
0173       std::vector<G4TwoVector> polygon;
0174       polygon.reserve(np);
0175       for(std::size_t i=0; i<np; ++i) {
0176         polygon.emplace_back(sh->GetX(i) * CM_2_MM,sh->GetY(i) * CM_2_MM);
0177       }
0178       return new G4ExtrudedSolid(sh->GetName(), polygon, z);
0179     }
0180
0181     template <> G4VSolid* convertShape<TGeoPgon>(const TGeoShape* shape)  {
0182       const TGeoPgon* sh = (const TGeoPgon*) shape;
0183       std::vector<double> rmin, rmax, z;
0184       for (Int_t i = 0; i < sh->GetNz(); ++i) {
0185         rmin.emplace_back(sh->GetRmin(i) * CM_2_MM);
0186         rmax.emplace_back(sh->GetRmax(i) * CM_2_MM);
0187         z.emplace_back(sh->GetZ(i) * CM_2_MM);
0188       }
0189       return new G4Polyhedra(sh->GetName(), sh->GetPhi1() * DEGREE_2_RAD, sh->GetDphi() * DEGREE_2_RAD,
0190                              sh->GetNedges(), sh->GetNz(), &z[0], &rmin[0], &rmax[0]);
0191     }
0192
0193     template <> G4VSolid* convertShape<TGeoPcon>(const TGeoShape* shape)  {
0194       const TGeoPcon* sh = (const TGeoPcon*) shape;
0195       std::vector<double> rmin, rmax, z;
0196       for (Int_t i = 0; i < sh->GetNz(); ++i) {
0197         rmin.emplace_back(sh->GetRmin(i) * CM_2_MM);
0198         rmax.emplace_back(sh->GetRmax(i) * CM_2_MM);
0199         z.emplace_back(sh->GetZ(i) * CM_2_MM);
0200       }
0201       return new G4Polycone(sh->GetName(), sh->GetPhi1() * DEGREE_2_RAD, sh->GetDphi() * DEGREE_2_RAD,
0202                             sh->GetNz(), &z[0], &rmin[0], &rmax[0]);
0203     }
0204
0205     template <> G4VSolid* convertShape<TGeoCone>(const TGeoShape* shape)  {
0206       const TGeoCone* sh = (const TGeoCone*) shape;
0207       return new G4Cons(sh->GetName(), sh->GetRmin1() * CM_2_MM, sh->GetRmax1() * CM_2_MM, sh->GetRmin2() * CM_2_MM,
0208                         sh->GetRmax2() * CM_2_MM, sh->GetDz() * CM_2_MM, 0.0, 2.*M_PI);
0209     }
0210
0211     template <> G4VSolid* convertShape<TGeoConeSeg>(const TGeoShape* shape)  {
0212       const TGeoConeSeg* sh = (const TGeoConeSeg*) shape;
0213       return new G4Cons(sh->GetName(), sh->GetRmin1() * CM_2_MM, sh->GetRmax1() * CM_2_MM,
0214                         sh->GetRmin2() * CM_2_MM, sh->GetRmax2() * CM_2_MM,
0215                         sh->GetDz() * CM_2_MM,
0216                         sh->GetPhi1() * DEGREE_2_RAD, (sh->GetPhi2()-sh->GetPhi1()) * DEGREE_2_RAD);
0217     }
0218
0219     template <> G4VSolid* convertShape<TGeoParaboloid>(const TGeoShape* shape)  {
0220       const TGeoParaboloid* sh = (const TGeoParaboloid*) shape;
0221       return new G4Paraboloid(sh->GetName(), sh->GetDz() * CM_2_MM, sh->GetRlo() * CM_2_MM, sh->GetRhi() * CM_2_MM);
0222     }
0223
0224     template <> G4VSolid* convertShape<TGeoSphere>(const TGeoShape* shape)  {
0225       const TGeoSphere* sh = (const TGeoSphere*) shape;
0226       return new G4Sphere(sh->GetName(), sh->GetRmin() * CM_2_MM, sh->GetRmax() * CM_2_MM,
0227                           sh->GetPhi1() * DEGREE_2_RAD, (sh->GetPhi2()-sh->GetPhi1()) * DEGREE_2_RAD,
0228                           sh->GetTheta1() * DEGREE_2_RAD, (sh->GetTheta2()- sh->GetTheta1()) * DEGREE_2_RAD);
0229     }
0230
0231     template <> G4VSolid* convertShape<TGeoTorus>(const TGeoShape* shape)  {
0232       const TGeoTorus* sh = (const TGeoTorus*) shape;
0233       return new G4Torus(sh->GetName(), sh->GetRmin() * CM_2_MM, sh->GetRmax() * CM_2_MM, sh->GetR() * CM_2_MM,
0234                          sh->GetPhi1() * DEGREE_2_RAD, sh->GetDphi() * DEGREE_2_RAD);
0235     }
0236
0237     template <> G4VSolid* convertShape<TGeoTrap>(const TGeoShape* shape)  {
0238       const TGeoTrap* sh = (const TGeoTrap*) shape;
0239       return new G4Trap(sh->GetName(), sh->GetDz() * CM_2_MM, sh->GetTheta() * DEGREE_2_RAD, sh->GetPhi() * DEGREE_2_RAD,
0240                         sh->GetH1() * CM_2_MM, sh->GetBl1() * CM_2_MM, sh->GetTl1() * CM_2_MM, sh->GetAlpha1() * DEGREE_2_RAD,
0241                         sh->GetH2() * CM_2_MM, sh->GetBl2() * CM_2_MM, sh->GetTl2() * CM_2_MM, sh->GetAlpha2() * DEGREE_2_RAD);
0242     }
0243
0244     template <> G4VSolid* convertShape<G4GenericTrap>(const TGeoShape* shape)  {
0245       std::vector<G4TwoVector> vertices;
0246       TGeoTrap* sh = (TGeoTrap*) shape;
0247       Double_t* vtx_xy = sh->GetVertices();
0248       for ( std::size_t i=0; i<8; ++i, vtx_xy +=2 )
0249         vertices.emplace_back(vtx_xy[0] * CM_2_MM, vtx_xy[1] * CM_2_MM);
0250       return new G4GenericTrap(sh->GetName(), sh->GetDz() * CM_2_MM, vertices);
0251     }
0252
0253     template <> G4VSolid* convertShape<TGeoTessellated>(const TGeoShape* shape)  {
0254       TGeoTessellated*   sh  = (TGeoTessellated*) shape;
0255       G4TessellatedSolid* g4 = new G4TessellatedSolid(sh->GetName());
0256       int num_facet = sh->GetNfacets();
0257
0258       printout(DEBUG,"TessellatedSolid","+++ %s> Converting %d facets", sh->GetName(), num_facet);
0259       for(int i=0; i<num_facet; ++i)  {
0260         const TGeoFacet& facet = sh->GetFacet(i);
0261         int nv = facet.GetNvert();
0262 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,31,1)
0263         const auto& v0 = sh->GetVertex(facet[0]);
0264         const auto& v1 = sh->GetVertex(facet[1]);
0265         const auto& v2 = sh->GetVertex(facet[2]);
0266 #else
0267         const auto& v0 = sh->GetVertex(facet.GetVertexIndex(0));
0268         const auto& v1 = sh->GetVertex(facet.GetVertexIndex(1));
0269         const auto& v2 = sh->GetVertex(facet.GetVertexIndex(2));
0270 #endif
0271         G4VFacet* g4f = 0;
0272         if ( nv == 3 )    {
0273           g4f = new G4TriangularFacet(G4ThreeVector(v0.x() * CM_2_MM, v0.y() * CM_2_MM, v0.z() * CM_2_MM),
0274                                       G4ThreeVector(v1.x() * CM_2_MM, v1.y() * CM_2_MM, v1.z() * CM_2_MM),
0275                                       G4ThreeVector(v2.x() * CM_2_MM, v2.y() * CM_2_MM, v2.z() * CM_2_MM),
0276                                       ABSOLUTE);
0277         }
0278         else if ( nv == 4 )    {
0279 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,31,1)
0280           const auto& v3 = sh->GetVertex(facet[3]);
0281 #else
0282           const auto& v3 = sh->GetVertex(facet.GetVertexIndex(3));
0283 #endif
0284           g4f = new G4QuadrangularFacet(G4ThreeVector(v0.x() * CM_2_MM, v0.y() * CM_2_MM, v0.z() * CM_2_MM),
0285                                         G4ThreeVector(v1.x() * CM_2_MM, v1.y() * CM_2_MM, v1.z() * CM_2_MM),
0286                                         G4ThreeVector(v2.x() * CM_2_MM, v2.y() * CM_2_MM, v2.z() * CM_2_MM),
0287                                         G4ThreeVector(v3.x() * CM_2_MM, v3.y() * CM_2_MM, v3.z() * CM_2_MM),
0288                                         ABSOLUTE);
0289         }
0290         else   {
0291           except("TGeoTessellated", "Tessellated shape [%s] has facet with wrong number of vertices: %d",
0292                  sh->GetName(), nv);
0293         }
0294         g4->AddFacet(g4f);
0295       }
0296       g4->SetSolidClosed(sh->IsClosedBody());
0297       return g4;
0298     }
0299
0300   }    // End namespace sim
0301 }      // End namespace dd4hep