Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-03-13 08:19:44

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 #define _USE_MATH_DEFINES
0015 
0016 // Framework include files
0017 #include <DD4hep/Shapes.h>
0018 #include <DD4hep/MatrixHelpers.h>
0019 #include <DD4hep/DD4hepUnits.h>
0020 #include <DD4hep/ShapeTags.h>
0021 #include <DD4hep/Printout.h>
0022 #include <DD4hep/detail/ShapesInterna.h>
0023 
0024 // C/C++ include files
0025 #include <stdexcept>
0026 #include <iomanip>
0027 #include <sstream>
0028 
0029 // ROOT includes
0030 #include <TClass.h>
0031 #include <TGeoMatrix.h>
0032 #include <TGeoBoolNode.h>
0033 #include <TGeoScaledShape.h>
0034 #include <TGeoCompositeShape.h>
0035 
0036 /// Namespace for the AIDA detector description toolkit
0037 namespace dd4hep {
0038 
0039   // name units differently
0040   namespace units = dd4hep;
0041 
0042   // Local, anonymous stuff
0043   namespace {
0044     template <typename T> inline T* get_ptr(const TGeoShape* shape)    {
0045       if ( shape && shape->IsA() == T::Class() ) return (T*)shape;
0046       except("Dimension","Invalid shape pointer!");
0047       return 0;
0048     }
0049     inline void invalidSetDimensionCall(const TGeoShape* sh, const std::vector<double>& params)   {
0050       except("Solid","+++ Shape:%s setDimension: Invalid number of parameters: %ld",
0051              (sh ? typeName(typeid(*sh)) : typeName(typeid(sh))).c_str(), params.size());
0052     }
0053     template <typename T> inline bool check_shape_type(const Handle<TGeoShape>& solid)   {
0054       return solid.isValid() && solid->IsA() == T::Class();
0055     }
0056   }
0057 
0058   /// Type check of various shapes. Result like dynamic_cast. Compare with python's isinstance(obj,type)
0059   template <typename SOLID> bool isInstance(const Handle<TGeoShape>& solid)   {
0060     return check_shape_type<typename SOLID::Object>(solid);
0061   }
0062   template bool isInstance<Box>               (const Handle<TGeoShape>& solid);
0063   template bool isInstance<Scale>             (const Handle<TGeoShape>& solid);
0064   template bool isInstance<ShapelessSolid>    (const Handle<TGeoShape>& solid);
0065   template bool isInstance<HalfSpace>         (const Handle<TGeoShape>& solid);
0066   template bool isInstance<ConeSegment>       (const Handle<TGeoShape>& solid);
0067   template bool isInstance<CutTube>           (const Handle<TGeoShape>& solid);
0068   template bool isInstance<EllipticalTube>    (const Handle<TGeoShape>& solid);
0069   template bool isInstance<Trap>              (const Handle<TGeoShape>& solid);
0070   template bool isInstance<Trd1>              (const Handle<TGeoShape>& solid);
0071   template bool isInstance<Trd2>              (const Handle<TGeoShape>& solid);
0072   template bool isInstance<Torus>             (const Handle<TGeoShape>& solid);
0073   template bool isInstance<Sphere>            (const Handle<TGeoShape>& solid);
0074   template bool isInstance<Paraboloid>        (const Handle<TGeoShape>& solid);
0075   template bool isInstance<Hyperboloid>       (const Handle<TGeoShape>& solid);
0076   template bool isInstance<PolyhedraRegular>  (const Handle<TGeoShape>& solid);
0077   template bool isInstance<Polyhedra>         (const Handle<TGeoShape>& solid);
0078   template bool isInstance<ExtrudedPolygon>   (const Handle<TGeoShape>& solid);
0079   template bool isInstance<TessellatedSolid>   (const Handle<TGeoShape>& solid);
0080   template bool isInstance<BooleanSolid>      (const Handle<TGeoShape>& solid);
0081 
0082   template <> bool isInstance<Cone>(const Handle<TGeoShape>& solid)  {
0083     return check_shape_type<TGeoConeSeg>(solid) || check_shape_type<TGeoCone>(solid);
0084   }
0085   template <> bool isInstance<Tube>(const Handle<TGeoShape>& solid)  {
0086     return check_shape_type<TGeoTubeSeg>(solid)
0087       || check_shape_type<TGeoCtub>(solid)
0088       || check_shape_type<TwistedTubeObject>(solid);
0089   }
0090   template <> bool isInstance<Polycone>(const Handle<TGeoShape>& solid)   {
0091     return check_shape_type<TGeoPcon>(solid)    || check_shape_type<TGeoPgon>(solid);
0092   }
0093   template <> bool isInstance<EightPointSolid>(const Handle<TGeoShape>& solid)   {
0094     if ( solid.isValid() )   {
0095       TClass* c = solid->IsA();
0096       return c==TGeoArb8::Class() || c==TGeoTrap::Class() || c==TGeoGtra::Class();
0097     }
0098     return false;
0099   }
0100   template <> bool isInstance<TwistedTube>(const Handle<TGeoShape>& solid)  {
0101     return check_shape_type<TwistedTubeObject>(solid);
0102   }
0103   template <> bool isInstance<TruncatedTube>(const Handle<TGeoShape>& solid)   {
0104     return check_shape_type<TGeoCompositeShape>(solid)
0105       &&   ::strcmp(solid->GetTitle(), TRUNCATEDTUBE_TAG) == 0;
0106   }
0107   template <> bool isInstance<PseudoTrap>(const Handle<TGeoShape>& solid)   {
0108     return check_shape_type<TGeoCompositeShape>(solid)
0109       &&   ::strcmp(solid->GetTitle(), PSEUDOTRAP_TAG) == 0;
0110   }
0111   template <> bool isInstance<SubtractionSolid>(const Handle<TGeoShape>& solid)   {
0112     TGeoCompositeShape* sh = (TGeoCompositeShape*)solid.ptr();
0113     return sh && sh->IsA() == TGeoCompositeShape::Class()
0114       &&   sh->GetBoolNode()->GetBooleanOperator() == TGeoBoolNode::kGeoSubtraction;
0115   }
0116   template <> bool isInstance<UnionSolid>(const Handle<TGeoShape>& solid)   {
0117     TGeoCompositeShape* sh = (TGeoCompositeShape*)solid.ptr();
0118     return sh && sh->IsA() == TGeoCompositeShape::Class()
0119       &&   sh->GetBoolNode()->GetBooleanOperator() == TGeoBoolNode::kGeoUnion;
0120   }
0121   template <> bool isInstance<IntersectionSolid>(const Handle<TGeoShape>& solid)   {
0122     TGeoCompositeShape* sh = (TGeoCompositeShape*)solid.ptr();
0123     return sh && sh->IsA() == TGeoCompositeShape::Class()
0124       &&   sh->GetBoolNode()->GetBooleanOperator() == TGeoBoolNode::kGeoIntersection;
0125   }
0126 
0127   /// Type check of various shapes. Do not allow for polymorphism. Types must match exactly
0128   template <typename SOLID> bool isA(const Handle<TGeoShape>& solid)   {
0129     return check_shape_type<typename SOLID::Object>(solid);
0130   }
0131   template bool isA<Box>(const Handle<TGeoShape>& solid);
0132   template bool isA<Scale>(const Handle<TGeoShape>& solid);
0133   template bool isA<ShapelessSolid>(const Handle<TGeoShape>& solid);
0134   template bool isA<HalfSpace>(const Handle<TGeoShape>& solid);
0135   template bool isA<Cone>(const Handle<TGeoShape>& solid);
0136   template bool isA<ConeSegment>(const Handle<TGeoShape>& solid);
0137   template bool isA<Tube>(const Handle<TGeoShape>& solid);
0138   template bool isA<CutTube>(const Handle<TGeoShape>& solid);
0139   template bool isA<EllipticalTube>(const Handle<TGeoShape>& solid);
0140   template bool isA<Trap>(const Handle<TGeoShape>& solid);
0141   template bool isA<Trd1>(const Handle<TGeoShape>& solid);
0142   template bool isA<Trd2>(const Handle<TGeoShape>& solid);
0143   template bool isA<Torus>(const Handle<TGeoShape>& solid);
0144   template bool isA<Sphere>(const Handle<TGeoShape>& solid);
0145   template bool isA<Paraboloid>(const Handle<TGeoShape>& solid);
0146   template bool isA<Hyperboloid>(const Handle<TGeoShape>& solid);
0147   template bool isA<PolyhedraRegular>(const Handle<TGeoShape>& solid);
0148   template bool isA<Polyhedra>(const Handle<TGeoShape>& solid);
0149   template bool isA<ExtrudedPolygon>(const Handle<TGeoShape>& solid);
0150   template bool isA<Polycone>(const Handle<TGeoShape>& solid);
0151   template bool isA<EightPointSolid>(const Handle<TGeoShape>& solid);
0152   template bool isA<TessellatedSolid>(const Handle<TGeoShape>& solid);
0153   template <> bool isA<TwistedTube>(const Handle<TGeoShape>& solid)   {
0154     return check_shape_type<TwistedTubeObject>(solid)
0155       &&   ::strcmp(solid->GetTitle(), TWISTEDTUBE_TAG) == 0;
0156   }
0157   template <> bool isA<TruncatedTube>(const Handle<TGeoShape>& solid)   {
0158     return check_shape_type<TGeoCompositeShape>(solid)
0159       &&   ::strcmp(solid->GetTitle(), TRUNCATEDTUBE_TAG) == 0;
0160   }
0161   template <> bool isA<PseudoTrap>(const Handle<TGeoShape>& solid)   {
0162     return check_shape_type<TGeoCompositeShape>(solid)
0163       &&   ::strcmp(solid->GetTitle(), PSEUDOTRAP_TAG) == 0;
0164   }
0165   template <> bool isA<SubtractionSolid>(const Handle<TGeoShape>& solid)   {
0166     TGeoCompositeShape* sh = (TGeoCompositeShape*)solid.ptr();
0167     return sh && sh->IsA() == TGeoCompositeShape::Class()
0168       &&   sh->GetBoolNode()->GetBooleanOperator() == TGeoBoolNode::kGeoSubtraction
0169       &&   ::strcmp(sh->GetTitle(), SUBTRACTION_TAG) == 0;
0170   }
0171   template <> bool isA<UnionSolid>(const Handle<TGeoShape>& solid)   {
0172     TGeoCompositeShape* sh = (TGeoCompositeShape*)solid.ptr();
0173     return sh && sh->IsA() == TGeoCompositeShape::Class()
0174       &&   sh->GetBoolNode()->GetBooleanOperator() == TGeoBoolNode::kGeoUnion
0175       &&   ::strcmp(sh->GetTitle(), UNION_TAG) == 0;
0176   }
0177   template <> bool isA<IntersectionSolid>(const Handle<TGeoShape>& solid)   {
0178     TGeoCompositeShape* sh = (TGeoCompositeShape*)solid.ptr();
0179     return sh && sh->IsA() == TGeoCompositeShape::Class()
0180       &&   sh->GetBoolNode()->GetBooleanOperator() == TGeoBoolNode::kGeoIntersection
0181       &&   ::strcmp(sh->GetTitle(), INTERSECTION_TAG) == 0;
0182   }
0183 
0184   /// Retrieve tag name from shape type
0185   std::string get_shape_tag(const TGeoShape* sh)    {
0186     if ( sh )    {
0187       TClass* cl = sh->IsA();
0188       if ( cl == TGeoShapeAssembly::Class() )
0189         return SHAPELESS_TAG;
0190       else if ( cl == TGeoBBox::Class() )
0191         return BOX_TAG;
0192       else if ( cl == TGeoHalfSpace::Class() )
0193         return HALFSPACE_TAG;
0194       else if ( cl == TGeoPgon::Class() )
0195         return POLYHEDRA_TAG;
0196       else if ( cl == TGeoPcon::Class() )
0197         return POLYCONE_TAG;
0198       else if ( cl == TGeoConeSeg::Class() )
0199         return CONESEGMENT_TAG;
0200       else if ( cl == TGeoCone::Class() )
0201         return CONE_TAG;
0202       else if ( cl == TGeoTube::Class() )
0203         return TUBE_TAG;
0204       else if ( cl == TGeoTubeSeg::Class() )
0205         return TUBE_TAG;
0206       else if ( cl == TGeoCtub::Class() )
0207         return CUTTUBE_TAG;
0208       else if ( cl == TGeoEltu::Class() )
0209         return ELLIPTICALTUBE_TAG;
0210       else if ( cl == TwistedTubeObject::Class() )
0211         return TWISTEDTUBE_TAG;
0212       else if ( cl == TGeoTrd1::Class() )
0213         return TRD1_TAG;
0214       else if ( cl == TGeoTrd2::Class() )
0215         return TRD2_TAG;
0216       else if ( cl == TGeoParaboloid::Class() )
0217         return PARABOLOID_TAG;
0218       else if ( cl == TGeoHype::Class() )
0219         return HYPERBOLOID_TAG;
0220       //else if ( cl == TGeoGtra::Class() )
0221       //  return 
0222       else if ( cl == TGeoTrap::Class() )
0223         return TRAP_TAG;
0224       else if ( cl == TGeoArb8::Class() )
0225         return EIGHTPOINTSOLID_TAG;
0226       else if ( cl == TGeoSphere::Class() )
0227         return SPHERE_TAG;
0228       else if ( cl == TGeoTorus::Class() )
0229         return TORUS_TAG;
0230       else if ( cl == TGeoXtru::Class() )
0231         return EXTRUDEDPOLYGON_TAG;
0232       else if ( cl == TGeoScaledShape::Class() )
0233         return SCALE_TAG;
0234       else if ( cl == TGeoTessellated::Class() )
0235         return TESSELLATEDSOLID_TAG;
0236       else if (isA<TruncatedTube>(sh) )
0237         return TRUNCATEDTUBE_TAG;
0238       else if (isA<PseudoTrap>(sh) )
0239         return PSEUDOTRAP_TAG;
0240       else if ( isInstance<SubtractionSolid>(sh) )
0241     return SUBTRACTION_TAG;
0242       else if ( isInstance<UnionSolid>(sh) )
0243     return UNION_TAG;
0244       else if ( isInstance<IntersectionSolid>(sh) )
0245     return INTERSECTION_TAG;
0246       else
0247         except("Solid","Failed to type of shape: %s [%s]",
0248            sh->GetName(), cl->GetName() );
0249       return "";
0250     }
0251     except("Solid","Failed to access dimensions [Invalid handle].");
0252     return "";
0253   }
0254 
0255   template <typename T> std::vector<double> dimensions(const TGeoShape* shape)    {
0256     std::stringstream str;
0257     if ( shape )
0258       str << "Shape: dimension(" << typeName(typeid(*shape)) << "): Invalid call!";
0259     else
0260       str << "Shape: dimension<TGeoShape): Invalid call && pointer!";
0261     throw std::runtime_error(str.str());
0262   }
0263   template <> std::vector<double> dimensions<TGeoShapeAssembly>(const TGeoShape* shape)    {
0264     const auto* sh = get_ptr<TGeoShapeAssembly>(shape);
0265     return { sh->GetDX(), sh->GetDY(), sh->GetDZ() };
0266   }
0267   template <> std::vector<double> dimensions<TGeoBBox>(const TGeoShape* shape)    {
0268     const auto* sh = get_ptr<TGeoBBox>(shape);
0269     return { sh->GetDX(), sh->GetDY(), sh->GetDZ() };
0270   }
0271   template <> std::vector<double> dimensions<TGeoHalfSpace>(const TGeoShape* shape)    {
0272     auto* sh = get_ptr<TGeoHalfSpace>(shape);
0273     return { sh->GetPoint()[0], sh->GetPoint()[1], sh->GetPoint()[2],
0274         sh->GetNorm()[0], sh->GetNorm()[1], sh->GetNorm()[2] };
0275   }
0276   template <> std::vector<double> dimensions<TGeoPcon>(const TGeoShape* shape)    {
0277     const TGeoPcon* sh = get_ptr<TGeoPcon>(shape);
0278     std::vector<double> pars { sh->GetPhi1()*units::deg, sh->GetDphi()*units::deg, double(sh->GetNz()) };
0279     pars.reserve(3+3*sh->GetNz());
0280     for (Int_t i = 0; i < sh->GetNz(); ++i) {
0281       pars.emplace_back(sh->GetZ(i));
0282       pars.emplace_back(sh->GetRmin(i));
0283       pars.emplace_back(sh->GetRmax(i));
0284     }
0285     return pars;
0286   }
0287   template <> std::vector<double> dimensions<TGeoConeSeg>(const TGeoShape* shape)    {
0288     const TGeoConeSeg* sh = get_ptr<TGeoConeSeg>(shape);
0289     return { sh->GetDz(), sh->GetRmin1(), sh->GetRmax1(), sh->GetRmin2(), sh->GetRmax2(),
0290         sh->GetPhi1()*units::deg, sh->GetPhi2()*units::deg };
0291   }
0292   template <> std::vector<double> dimensions<TGeoCone>(const TGeoShape* shape)    {
0293     const TGeoCone* sh = get_ptr<TGeoCone>(shape);
0294     return { sh->GetDz(), sh->GetRmin1(), sh->GetRmax1(), sh->GetRmin2(), sh->GetRmax2() };
0295   }  
0296   template <> std::vector<double> dimensions<TGeoTube>(const TGeoShape* shape)    {
0297     const TGeoTube* sh = get_ptr<TGeoTube>(shape);
0298     return { sh->GetRmin(), sh->GetRmax(), sh->GetDz() };
0299   }
0300   template <> std::vector<double> dimensions<TGeoTubeSeg>(const TGeoShape* shape)    {
0301     const TGeoTubeSeg* sh = get_ptr<TGeoTubeSeg>(shape);
0302     return { sh->GetRmin(), sh->GetRmax(), sh->GetDz(), sh->GetPhi1()*units::deg, sh->GetPhi2()*units::deg };
0303   }
0304   template <> std::vector<double> dimensions<TwistedTubeObject>(const TGeoShape* shape)    {
0305     const TwistedTubeObject* sh = get_ptr<TwistedTubeObject>(shape);
0306     return { sh->GetPhiTwist(), sh->GetRmin(), sh->GetRmax(),
0307         sh->GetNegativeEndZ(), sh->GetPositiveEndZ(),
0308         double(sh->GetNsegments()), sh->GetPhi2()*units::deg };
0309   }
0310   template <> std::vector<double> dimensions<TGeoCtub>(const TGeoShape* shape)    {
0311     const TGeoCtub* sh = get_ptr<TGeoCtub>(shape);
0312     const Double_t* lo = sh->GetNlow();
0313     const Double_t* hi = sh->GetNhigh();
0314     return { sh->GetRmin(), sh->GetRmax(), sh->GetDz(),
0315         sh->GetPhi1()*units::deg, sh->GetPhi2()*units::deg,
0316         lo[0], lo[1], lo[2], hi[0], hi[1], hi[2] };
0317   }
0318   template <> std::vector<double> dimensions<TGeoEltu>(const TGeoShape* shape)    {
0319     const TGeoEltu* sh = get_ptr<TGeoEltu>(shape);
0320     return { sh->GetA(), sh->GetB(), sh->GetDz() };
0321   }
0322   template <> std::vector<double> dimensions<TGeoTrd1>(const TGeoShape* shape)    {
0323     const TGeoTrd1* sh = get_ptr<TGeoTrd1>(shape);
0324     return { sh->GetDx1(), sh->GetDx2(), sh->GetDy(), sh->GetDz() };
0325   }
0326   template <> std::vector<double> dimensions<TGeoTrd2>(const TGeoShape* shape)    {
0327     const TGeoTrd2* sh = get_ptr<TGeoTrd2>(shape);
0328     return { sh->GetDx1(), sh->GetDx2(), sh->GetDy1(), sh->GetDy2(), sh->GetDz() };
0329   }
0330   template <> std::vector<double> dimensions<TGeoParaboloid>(const TGeoShape* shape)    {
0331     const TGeoParaboloid* sh = get_ptr<TGeoParaboloid>(shape);
0332     return { sh->GetRlo(), sh->GetRhi(), sh->GetDz() };
0333   }
0334   template <> std::vector<double> dimensions<TGeoHype>(const TGeoShape* shape)    {
0335     const TGeoHype* sh = get_ptr<TGeoHype>(shape);
0336     return { sh->GetDz(),
0337         sh->GetRmin(), sh->GetStIn()*units::deg,
0338         sh->GetRmax(), sh->GetStOut()*units::deg };
0339   }
0340   template <> std::vector<double> dimensions<TGeoSphere>(const TGeoShape* shape)    {
0341     const TGeoSphere* sh = get_ptr<TGeoSphere>(shape);
0342     return { sh->GetRmin(), sh->GetRmax(),
0343         sh->GetTheta1()*units::deg, sh->GetTheta2()*units::deg,
0344         sh->GetPhi1()*units::deg,   sh->GetPhi2()*units::deg };
0345   }
0346   template <> std::vector<double> dimensions<TGeoTorus>(const TGeoShape* shape)    {
0347     const TGeoTorus* sh = get_ptr<TGeoTorus>(shape);
0348     return { sh->GetR(), sh->GetRmin(), sh->GetRmax(),
0349         sh->GetPhi1()*units::deg, sh->GetDphi()*units::deg };
0350   }
0351   template <> std::vector<double> dimensions<TGeoPgon>(const TGeoShape* shape)    {
0352     const TGeoPgon* sh = get_ptr<TGeoPgon>(shape);
0353     std::vector<double> pars { sh->GetPhi1()*units::deg, sh->GetDphi()*units::deg, double(sh->GetNedges()), double(sh->GetNz()) };
0354     pars.reserve(4+3*sh->GetNz());
0355     for (Int_t i = 0; i < sh->GetNz(); ++i) {
0356       pars.emplace_back(sh->GetZ(i));
0357       pars.emplace_back(sh->GetRmin(i));
0358       pars.emplace_back(sh->GetRmax(i));
0359     }
0360     return pars;
0361   }
0362   template <> std::vector<double> dimensions<TGeoXtru>(const TGeoShape* shape)    {
0363     const TGeoXtru* sh = get_ptr<TGeoXtru>(shape);
0364     Int_t nz = sh->GetNz();
0365     std::vector<double> pars { double(nz) };
0366     pars.reserve(1+4*nz);
0367     for(Int_t i=0; i<nz; ++i)   {
0368       pars.emplace_back(sh->GetZ(i));
0369       pars.emplace_back(sh->GetXOffset(i));
0370       pars.emplace_back(sh->GetYOffset(i));
0371       pars.emplace_back(sh->GetScale(i));
0372     }
0373     return pars;
0374   }
0375   template <> std::vector<double> dimensions<TGeoArb8>(const TGeoShape* shape)    {
0376     TGeoArb8* sh = get_ptr<TGeoArb8>(shape);
0377     struct _V { double xy[8][2]; } *vtx = (_V*)sh->GetVertices();
0378     std::vector<double> pars { sh->GetDz() };
0379     for ( std::size_t i=0; i<8; ++i )   {
0380       pars.emplace_back(vtx->xy[i][0]);
0381       pars.emplace_back(vtx->xy[i][1]);
0382     }
0383     return pars;
0384   }
0385   template <> std::vector<double> dimensions<TGeoTrap>(const TGeoShape* shape)    {
0386     const TGeoTrap* sh = get_ptr<TGeoTrap>(shape);
0387     return { sh->GetDz(), sh->GetTheta()*units::deg, sh->GetPhi()*units::deg,
0388         sh->GetH1(), sh->GetBl1(), sh->GetTl1(), sh->GetAlpha1()*units::deg,
0389         sh->GetH2(), sh->GetBl2(), sh->GetTl2(), sh->GetAlpha2()*units::deg };
0390   }
0391   template <> std::vector<double> dimensions<TGeoGtra>(const TGeoShape* shape)    {
0392     const TGeoGtra* sh = get_ptr<TGeoGtra>(shape);
0393     return { sh->GetDz(), sh->GetTheta()*units::deg, sh->GetPhi()*units::deg,
0394         sh->GetH1(), sh->GetBl1(), sh->GetTl1(), sh->GetAlpha1()*units::deg,
0395         sh->GetH2(), sh->GetBl2(), sh->GetTl2(), sh->GetAlpha2()*units::deg,
0396         sh->GetTwistAngle()*units::deg };
0397   }
0398   template <> std::vector<double> dimensions<TGeoScaledShape>(const TGeoShape* shape)    {
0399     TGeoScaledShape* sh = get_ptr<TGeoScaledShape>(shape);
0400     TGeoShape*       s_sh = sh->GetShape();
0401     const Double_t*  scale = sh->GetScale()->GetScale();
0402     std::vector<double>   pars {scale[0],scale[1],scale[2]};
0403     std::vector<double>   s_pars = get_shape_dimensions(s_sh);
0404     for(auto p : s_pars) pars.push_back(p);
0405     return pars;
0406   }
0407   
0408   template <> std::vector<double> dimensions<TGeoTessellated>(const TGeoShape* shape)    {
0409     TGeoTessellated* sh = get_ptr<TGeoTessellated>(shape);
0410     int num_facet = sh->GetNfacets();
0411     int num_vtx   = sh->GetNvertices();
0412     std::vector<double> pars;
0413 
0414     printout(DEBUG,"TessellatedSolid","+++ Saving %d vertices, %d facets",num_vtx, num_facet);
0415     pars.reserve(num_facet*5+num_vtx*3+2);
0416     pars.emplace_back(num_vtx);
0417     pars.emplace_back(num_facet);
0418     for(int i=0; i<num_vtx; ++i)  {
0419       const auto& v = sh->GetVertex(i);
0420       pars.emplace_back(v.x());
0421       pars.emplace_back(v.y());
0422       pars.emplace_back(v.z());
0423     }
0424     for(int i=0; i<num_facet; ++i)  {
0425       const TGeoFacet& f = sh->GetFacet(i);
0426       pars.emplace_back(double(f.GetNvert()));
0427       for(int j=0, n=f.GetNvert(); j<n; ++j)   {
0428 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,31,1)
0429     int idx = f[j];
0430         pars.emplace_back(double(idx));
0431 #else
0432         pars.emplace_back(double(f.GetVertexIndex(j)));
0433 #endif
0434       }
0435     }
0436     return pars;
0437   }
0438   
0439   template <> std::vector<double> dimensions<TGeoCompositeShape>(const TGeoShape* shape)    {
0440     const TGeoCompositeShape* sh = get_ptr<TGeoCompositeShape>(shape);
0441     const TGeoBoolNode*  boolean = sh->GetBoolNode();
0442     TGeoMatrix*    right_matrix  = boolean->GetRightMatrix();
0443     TGeoShape*     left_solid    = boolean->GetLeftShape();
0444     TGeoShape*     right_solid   = boolean->GetRightShape();
0445 
0446     TGeoBoolNode::EGeoBoolType oper = boolean->GetBooleanOperator();
0447     TGeoMatrix*   left_matrix   = boolean->GetRightMatrix();
0448     const Double_t*   left_tr   = left_matrix->GetTranslation();
0449     const Double_t*   left_rot  = left_matrix->GetRotationMatrix();
0450     const Double_t*   right_tr  = right_matrix->GetTranslation();
0451     const Double_t*   right_rot = right_matrix->GetRotationMatrix();
0452 
0453     std::vector<double>    pars { double(oper) };
0454     std::vector<double>    left_par  = Solid(left_solid).dimensions();
0455     std::vector<double>    right_par = Solid(right_solid).dimensions();
0456 
0457     pars.insert(pars.end(), left_par.begin(), left_par.end());
0458     pars.insert(pars.end(), left_rot, left_rot+9);
0459     pars.insert(pars.end(), left_tr, left_tr+3);
0460 
0461     pars.insert(pars.end(), right_par.begin(), right_par.end());
0462     pars.insert(pars.end(), right_rot, right_rot+9);
0463     pars.insert(pars.end(), right_tr, right_tr+3);
0464     return pars;
0465   }
0466 
0467   template <typename T> std::vector<double> dimensions(const Handle<TGeoShape>& shape)
0468   {  return dimensions<typename T::Object>(get_ptr<typename T::Object>(shape.ptr()));  }
0469   template std::vector<double> dimensions<ShapelessSolid>   (const Handle<TGeoShape>& shape);
0470   template std::vector<double> dimensions<Box>              (const Handle<TGeoShape>& shape);
0471   template std::vector<double> dimensions<Scale>            (const Handle<TGeoShape>& shape);
0472   template std::vector<double> dimensions<HalfSpace>        (const Handle<TGeoShape>& shape);
0473   template std::vector<double> dimensions<Polycone>         (const Handle<TGeoShape>& shape);
0474   template std::vector<double> dimensions<ConeSegment>      (const Handle<TGeoShape>& shape);
0475   template std::vector<double> dimensions<Tube>             (const Handle<TGeoShape>& shape);
0476   template std::vector<double> dimensions<CutTube>          (const Handle<TGeoShape>& shape);
0477   template std::vector<double> dimensions<TwistedTube>      (const Handle<TGeoShape>& shape);
0478   template std::vector<double> dimensions<EllipticalTube>   (const Handle<TGeoShape>& shape);
0479   template std::vector<double> dimensions<Cone>             (const Handle<TGeoShape>& shape);
0480   template std::vector<double> dimensions<Trap>             (const Handle<TGeoShape>& shape);
0481   template std::vector<double> dimensions<Trd1>             (const Handle<TGeoShape>& shape);
0482   template std::vector<double> dimensions<Trd2>             (const Handle<TGeoShape>& shape);
0483   template std::vector<double> dimensions<Torus>            (const Handle<TGeoShape>& shape);
0484   template std::vector<double> dimensions<Sphere>           (const Handle<TGeoShape>& shape);
0485   template std::vector<double> dimensions<Paraboloid>       (const Handle<TGeoShape>& shape);
0486   template std::vector<double> dimensions<Hyperboloid>      (const Handle<TGeoShape>& shape);
0487   template std::vector<double> dimensions<PolyhedraRegular> (const Handle<TGeoShape>& shape);
0488   template std::vector<double> dimensions<Polyhedra>        (const Handle<TGeoShape>& shape);
0489   template std::vector<double> dimensions<ExtrudedPolygon>  (const Handle<TGeoShape>& shape);
0490   template std::vector<double> dimensions<EightPointSolid>  (const Handle<TGeoShape>& shape); 
0491   template std::vector<double> dimensions<TessellatedSolid>  (const Handle<TGeoShape>& shape); 
0492   template std::vector<double> dimensions<BooleanSolid>     (const Handle<TGeoShape>& shape);
0493   template std::vector<double> dimensions<SubtractionSolid> (const Handle<TGeoShape>& shape);
0494   template std::vector<double> dimensions<UnionSolid>       (const Handle<TGeoShape>& shape);
0495   template std::vector<double> dimensions<IntersectionSolid>(const Handle<TGeoShape>& shape);
0496 
0497   template <> std::vector<double> dimensions<PseudoTrap>(const Handle<TGeoShape>& shape)   {
0498     const TGeoCompositeShape* sh = get_ptr<TGeoCompositeShape>(shape.ptr());
0499     TGeoMatrix*    right_matrix  = sh->GetBoolNode()->GetRightMatrix();
0500     std::stringstream params(right_matrix->GetTitle());
0501     std::vector<double> pars;
0502     pars.reserve(7);
0503 #ifdef DIMENSION_DEBUG
0504     cout << "dimensions: [" << PSEUDOTRAP_TAG << "]" << endl
0505          << right_matrix->GetTitle() << endl;
0506 #endif
0507     for(std::size_t i=0; i<7; ++i)   {
0508       double val;
0509       params >> val;
0510       pars.emplace_back(val);
0511       if ( !params.good() )   {
0512         except("PseudoTrap","+++ dimensions: Invalid parameter stream.");
0513       }
0514     }
0515     return pars;
0516   }
0517 
0518   template <> std::vector<double> dimensions<TruncatedTube>(const Handle<TGeoShape>& shape)   {
0519     const TGeoCompositeShape* sh = get_ptr<TGeoCompositeShape>(shape.ptr());
0520     TGeoMatrix* right_matrix  = sh->GetBoolNode()->GetRightMatrix();
0521     std::stringstream params(right_matrix->GetTitle());
0522     std::vector<double> pars;
0523     pars.reserve(8);
0524     for(std::size_t i=0; i<8; ++i)   {
0525       double val;
0526       params >> val;
0527       pars.emplace_back(val);
0528       if ( !params.good() )   {
0529         except("TruncatedTube","+++ dimensions: Invalid parameter stream.");
0530       }
0531     }
0532     return pars;
0533   }
0534   
0535   template <> std::vector<double> dimensions<Solid>(const Handle<TGeoShape>& shape)   {
0536     if ( shape.ptr() )   {
0537       TClass* cl = shape->IsA();
0538       if ( cl == TGeoShapeAssembly::Class() )
0539         return dimensions<TGeoShapeAssembly>(shape.ptr() );
0540       else if ( cl == TGeoBBox::Class() )
0541         return dimensions<TGeoBBox>(shape.ptr() );
0542       else if ( cl == TGeoHalfSpace::Class() )
0543         return dimensions<TGeoHalfSpace>(shape.ptr() );
0544       else if ( cl == TGeoPgon::Class() )
0545         return dimensions<TGeoPgon>(shape.ptr() );
0546       else if ( cl == TGeoPcon::Class() )
0547         return dimensions<TGeoPcon>(shape.ptr() );
0548       else if ( cl == TGeoConeSeg::Class() )
0549         return dimensions<TGeoConeSeg>(shape.ptr() );
0550       else if ( cl == TGeoCone::Class() )
0551         return dimensions<TGeoCone>(shape.ptr() );
0552       else if ( cl == TGeoTube::Class() )
0553         return dimensions<TGeoTube>(shape.ptr() );
0554       else if ( cl == TGeoTubeSeg::Class() )
0555         return dimensions<TGeoTubeSeg>(shape.ptr() );
0556       else if ( cl == TGeoCtub::Class() )
0557         return dimensions<TGeoCtub>(shape.ptr() );
0558       else if ( cl == TGeoEltu::Class() )
0559         return dimensions<TGeoEltu>(shape.ptr() );
0560       else if ( cl == TwistedTubeObject::Class() )
0561         return dimensions<TwistedTubeObject>(shape.ptr() );
0562       else if ( cl == TGeoTrd1::Class() )
0563         return dimensions<TGeoTrd1>(shape.ptr() );
0564       else if ( cl == TGeoTrd2::Class() )
0565         return dimensions<TGeoTrd2>(shape.ptr() );
0566       else if ( cl == TGeoParaboloid::Class() )
0567         return dimensions<TGeoParaboloid>(shape.ptr() );
0568       else if ( cl == TGeoHype::Class() )
0569         return dimensions<TGeoHype>(shape.ptr() );
0570       else if ( cl == TGeoGtra::Class() )
0571         return dimensions<TGeoGtra>(shape.ptr() );
0572       else if ( cl == TGeoTrap::Class() )
0573         return dimensions<TGeoTrap>(shape.ptr() );
0574       else if ( cl == TGeoArb8::Class() )
0575         return dimensions<TGeoArb8>(shape.ptr() );
0576       else if ( cl == TGeoSphere::Class() )
0577         return dimensions<TGeoSphere>(shape.ptr() );
0578       else if ( cl == TGeoTorus::Class() )
0579         return dimensions<TGeoTorus>(shape.ptr() );
0580       else if ( cl == TGeoXtru::Class() )
0581         return dimensions<TGeoXtru>(shape.ptr() );
0582       else if ( cl == TGeoScaledShape::Class() )
0583         return dimensions<TGeoScaledShape>(shape.ptr() );
0584       else if ( cl == TGeoTessellated::Class() )
0585         return dimensions<TGeoTessellated>(shape.ptr() );
0586       else if (isA<TruncatedTube>(shape.ptr() ))
0587         return dimensions<TruncatedTube>(shape);
0588       else if (isA<PseudoTrap>(shape.ptr() ))
0589         return dimensions<PseudoTrap>(shape);
0590       else if ( cl == TGeoCompositeShape::Class() )
0591         return dimensions<TGeoCompositeShape>(shape.ptr() );
0592       else  {
0593         printout(ERROR,"Solid","Failed to access dimensions for shape of type:%s.",
0594                  cl->GetName());
0595       }
0596       return {};
0597     }
0598     except("Solid","Failed to access dimensions [Invalid handle].");
0599     return {};
0600   }
0601 
0602   /// Access Shape dimension parameters (As in TGeo, but angles in radians rather than degrees)
0603   std::vector<double> get_shape_dimensions(const Handle<TGeoShape>& shape)    {
0604     return dimensions<Solid>(shape);
0605   }
0606   /// Access Shape dimension parameters (As in TGeo, but angles in radians rather than degrees)
0607   std::vector<double> get_shape_dimensions(TGeoShape* shape)   {
0608     return dimensions<Solid>(Solid(shape));
0609   }
0610 
0611   template <typename T> void set_dimensions(T shape, const std::vector<double>& )    {
0612     std::stringstream str;
0613     if ( shape )
0614       str << "Shape: set_dimension(" << typeName(typeid(*shape)) << "): Invalid call!";
0615     else
0616       str << "Shape: set_dimension<TGeoShape): Invalid call && pointer!";
0617     throw std::runtime_error(str.str());
0618   }
0619 
0620   template <> void set_dimensions(TGeoShapeAssembly* , const std::vector<double>& )    {
0621     printout(WARNING,"ShapelessSolid","+++ setDimensions is a compatibility call. Nothing implemented.");
0622   }
0623   template <> void set_dimensions(TGeoBBox* sh, const std::vector<double>& params)    {
0624     auto pars = params;
0625     if ( pars.size() != 3 )   {
0626       invalidSetDimensionCall(sh,params);
0627     }
0628     Solid(sh)._setDimensions(&pars[0]);
0629   }
0630   template <> void set_dimensions(TGeoHalfSpace* sh, const std::vector<double>& params)   {
0631     auto pars = params;
0632     if ( params.size() != 6 )   {
0633       invalidSetDimensionCall(sh,params);
0634     }
0635     Solid(sh)._setDimensions(&pars[0]);
0636   }
0637   template <> void set_dimensions(TGeoPcon* sh, const std::vector<double>& params)   {
0638     if ( params.size() < 3 )   {
0639       invalidSetDimensionCall(sh,params);
0640     }
0641     std::size_t nz = std::size_t(params[2]);
0642     if ( params.size() != 3 + 3*nz )   {
0643       invalidSetDimensionCall(sh,params);
0644     }
0645     std::vector<double> pars = params;
0646     pars[0] /= units::deg;
0647     pars[1] /= units::deg;
0648     Solid(sh)._setDimensions(&pars[0]);
0649   }
0650   template <> void set_dimensions(TGeoConeSeg* sh, const std::vector<double>& params)   {
0651     if ( params.size() != 7 )   {
0652       invalidSetDimensionCall(sh,params);
0653     }
0654     auto pars = params;
0655     pars[5] /= units::deg;
0656     pars[6] /= units::deg;
0657     Solid(sh)._setDimensions(&pars[0]);
0658   }
0659   template <> void set_dimensions(TGeoCone* sh, const std::vector<double>& params)   {
0660     auto pars = params;
0661     if ( params.size() != 5 )   {
0662       invalidSetDimensionCall(sh,params);
0663     }
0664     Solid(sh)._setDimensions(&pars[0]);
0665   }
0666   template <> void set_dimensions(TGeoTube* sh, const std::vector<double>& params)   {
0667     if ( params.size() != 3 )   {
0668       invalidSetDimensionCall(sh,params);
0669     }
0670     auto pars = params;
0671     Solid(sh)._setDimensions(&pars[0]);
0672   }
0673   template <> void set_dimensions(TGeoTubeSeg* sh, const std::vector<double>& params)   {
0674     if ( params.size() != 5 )   {
0675       invalidSetDimensionCall(sh,params);
0676     }
0677     auto pars = params;
0678     pars[3] /= units::deg;
0679     pars[4] /= units::deg;
0680     Solid(sh)._setDimensions(&pars[0]);
0681   }
0682   template <> void set_dimensions(TwistedTubeObject* sh, const std::vector<double>& params)   {
0683     if ( params.size() != 7 )   {
0684       invalidSetDimensionCall(sh,params);
0685     }
0686     auto pars = params;
0687     sh->fPhiTwist = pars[0]/units::deg;
0688     sh->fNegativeEndz = pars[3];
0689     sh->fPositiveEndz = pars[4];
0690     sh->fNsegments    = (int)pars[5];
0691 
0692     pars[0] = pars[1];
0693     pars[1] = pars[2];
0694     pars[2] = (-pars[3]+pars[4])/2.0;
0695     pars[3] = 0.0;
0696     pars[4] = pars[6];
0697     pars.resize(5);
0698     set_dimensions((TGeoTubeSeg*)sh, pars);
0699   }
0700   template <> void set_dimensions(TGeoCtub* sh, const std::vector<double>& params)   {
0701     if ( params.size() != 11 )   {
0702       invalidSetDimensionCall(sh,params);
0703     }
0704     auto pars = params;
0705     pars[3] /= units::deg;
0706     pars[4] /= units::deg;
0707     Solid(sh)._setDimensions(&pars[0]);
0708   }
0709   template <> void set_dimensions(TGeoEltu* sh, const std::vector<double>& params)   {
0710     auto pars = params;
0711     if ( params.size() != 3 )   {
0712       invalidSetDimensionCall(sh,params);
0713     }
0714     Solid(sh)._setDimensions(&pars[0]);
0715   }
0716   template <> void set_dimensions(TGeoTrd1* sh, const std::vector<double>& params)   {
0717     auto pars = params;
0718     if ( params.size() != 4 )   {
0719       invalidSetDimensionCall(sh,params);
0720     }
0721     Solid(sh)._setDimensions(&pars[0]);
0722   }
0723   template <> void set_dimensions(TGeoTrd2* sh, const std::vector<double>& params)   {
0724     auto pars = params;
0725     if ( params.size() != 5 )   {
0726       invalidSetDimensionCall(sh,params);
0727     }
0728     Solid(sh)._setDimensions(&pars[0]);
0729   }
0730   template <> void set_dimensions(TGeoParaboloid* sh, const std::vector<double>& params)   {
0731     auto pars = params;
0732     if ( params.size() != 3 )   {
0733       invalidSetDimensionCall(sh,params);
0734     }
0735     Solid(sh)._setDimensions(&pars[0]);
0736   }
0737   template <> void set_dimensions(TGeoHype* sh, const std::vector<double>& params)   {
0738     if ( params.size() != 5 )   {
0739       invalidSetDimensionCall(sh,params);
0740     }
0741     double pars[] = { params[0], params[1], params[2]/units::deg, params[3], params[4]/units::deg };
0742     Solid(sh)._setDimensions(pars);
0743   }
0744   template <> void set_dimensions(TGeoSphere* sh, const std::vector<double>& params)   {
0745     if ( params.size() < 2 )   {
0746       invalidSetDimensionCall(sh,params);
0747     }
0748     double pars[] = { params[0], params[1], 0.0, M_PI/units::deg, 0.0, 2*M_PI/units::deg };
0749     if (params.size() > 2) pars[2] = params[2]/units::deg;
0750     if (params.size() > 3) pars[3] = params[3]/units::deg;
0751     if (params.size() > 4) pars[4] = params[4]/units::deg;
0752     if (params.size() > 5) pars[5] = params[5]/units::deg;
0753     Sphere(sh).access()->SetDimensions(pars, params.size());
0754     sh->ComputeBBox();
0755   }
0756   template <> void set_dimensions(TGeoTorus* sh, const std::vector<double>& params)   {
0757     if ( params.size() != 5 )   {
0758       invalidSetDimensionCall(sh,params);
0759     }
0760     double pars[] = { params[0], params[1], params[2], params[3]/units::deg, params[4]/units::deg };
0761     Solid(sh)._setDimensions(&pars[0]);
0762   }
0763   template <> void set_dimensions(TGeoPgon* sh, const std::vector<double>& params)   {
0764     auto pars = params;
0765     if ( params.size() < 4 || params.size() != 4 + 3*std::size_t(params[3]) )   {
0766       invalidSetDimensionCall(sh,params);
0767     }
0768     pars[0]  /= units::deg;
0769     pars[1]  /= units::deg;
0770     Solid(sh)._setDimensions(&pars[0]);
0771   }
0772   template <> void set_dimensions(TGeoXtru* sh, const std::vector<double>& params)   {
0773     auto pars = params;
0774     if ( params.size() < 1 || params.size() != 1 + 4*std::size_t(params[0]) )   {
0775       invalidSetDimensionCall(sh,params);
0776     }
0777     Solid(sh)._setDimensions(&pars[0]);
0778   }
0779   template <> void set_dimensions(TGeoArb8* sh, const std::vector<double>& params)   {
0780     if ( params.size() != 17 )   {
0781       invalidSetDimensionCall(sh,params);
0782     }
0783     auto pars = params;
0784     Solid(sh)._setDimensions(&pars[0]);
0785   }
0786   template <> void set_dimensions(TGeoTrap* sh, const std::vector<double>& params)   {
0787     auto pars = params;
0788     if ( params.size() != 11 )   {
0789       invalidSetDimensionCall(sh,params);
0790     }
0791     pars[1]  /= units::deg;
0792     pars[2]  /= units::deg;
0793     pars[6]  /= units::deg;
0794     pars[10] /= units::deg;
0795     Solid(sh)._setDimensions(&pars[0]);
0796   }
0797   template <> void set_dimensions(TGeoGtra* sh, const std::vector<double>& params)   {
0798     auto pars = params;
0799     if ( params.size() != 12 )   {
0800       invalidSetDimensionCall(sh,params);
0801     }
0802     pars[1]  /= units::deg;
0803     pars[2]  /= units::deg;
0804     pars[6]  /= units::deg;
0805     pars[10] /= units::deg;
0806     pars[11] /= units::deg;
0807     Solid(sh)._setDimensions(&pars[0]);
0808   }
0809   template <> void set_dimensions(TGeoScaledShape* sh, const std::vector<double>& params)   {
0810     Solid            s_sh(sh->GetShape());
0811     sh->GetScale()->SetScale(params[0], params[1], params[2]);
0812     auto pars = params;
0813     s_sh.access()->SetDimensions(&pars[3]);
0814   }
0815   
0816   template <> void set_dimensions(TGeoTessellated* sh, const std::vector<double>& params)    {
0817     int num_vtx   = params[0];
0818     int num_facet = params[1];
0819     std::vector<TessellatedSolid::Vertex> vertices;
0820     std::size_t i_par = 1;
0821     printout(DEBUG,"TessellatedSolid","+++ Loading %d vertices, %d facets",num_vtx, num_facet);
0822     for (int i=0; i<num_vtx; ++i)   {
0823       double x = params[++i_par];
0824       double y = params[++i_par];
0825       double z = params[++i_par];
0826       vertices.emplace_back(x, y, z);
0827     }
0828     std::string nam = sh->GetName();
0829     std::string tit = sh->GetTitle();
0830     sh->~TGeoTessellated();
0831     new(sh) TGeoTessellated(nam.c_str(), vertices);
0832     sh->SetTitle(tit.c_str());
0833     for (int i=0; i<num_facet; ++i)   {
0834       int i0, i1, i2, i3;
0835       int n_vtx = params[++i_par];
0836       i0 = n_vtx>0 ? params[++i_par] : -1;
0837       i1 = n_vtx>1 ? params[++i_par] : -1;
0838       i2 = n_vtx>2 ? params[++i_par] : -1;
0839       i3 = n_vtx>3 ? params[++i_par] : -1;
0840       if ( n_vtx == 3 )
0841         sh->AddFacet(i0,i1,i2);
0842       else if ( n_vtx == 4 )
0843         sh->AddFacet(i0,i1,i2,i3);
0844     }
0845     sh->CloseShape(true, true, false);
0846   }
0847   
0848   template <> void set_dimensions(TGeoCompositeShape*, const std::vector<double>&)   {
0849     // In general TGeoCompositeShape instances have an empty SetDimension call
0850 #if 0
0851     TGeoCompositeShape* sh      = (TGeoCompositeShape*) shape;
0852     TGeoBoolNode* boolean       = sh->GetBoolNode();
0853     TGeoMatrix*   right_matrix  = boolean->GetRightMatrix();
0854     TGeoShape*    left_solid    = boolean->GetLeftShape();
0855     TGeoShape*    right_solid   = boolean->GetRightShape();
0856     auto pars = params;
0857     Solid(sh)._setDimensions(&pars[0]);
0858 #endif
0859 #ifdef DIMENSION_DEBUG
0860     throw std::runtime_error("Composite shape. setDimensions is not implemented!");
0861 #endif
0862   }
0863 
0864   template <> void set_dimensions(ShapelessSolid shape, const std::vector<double>& params)
0865   {  set_dimensions(shape.ptr(), params);   }
0866   template <> void set_dimensions(Box shape, const std::vector<double>& params)
0867   {  set_dimensions(shape.ptr(), params);   }
0868   template <> void set_dimensions(Scale shape, const std::vector<double>& params)
0869   {  set_dimensions(shape.ptr(), params);   }
0870   template <> void set_dimensions(HalfSpace shape, const std::vector<double>& params)
0871   {  set_dimensions(shape.ptr(), params);   }
0872   template <> void set_dimensions(Polycone shape, const std::vector<double>& params)
0873   {  set_dimensions(shape.ptr(), params);   }
0874   template <> void set_dimensions(Cone shape, const std::vector<double>& params)
0875   {  set_dimensions(shape.ptr(), params);   }
0876   template <> void set_dimensions(ConeSegment shape, const std::vector<double>& params)
0877   {  set_dimensions(shape.ptr(), params);   }
0878   template <> void set_dimensions(Tube shape, const std::vector<double>& params)
0879   {  set_dimensions(shape.ptr(), params);   }
0880   template <> void set_dimensions(CutTube shape, const std::vector<double>& params)
0881   {  set_dimensions(shape.ptr(), params);   }
0882   template <> void set_dimensions(TwistedTube shape, const std::vector<double>& params)
0883   {  set_dimensions((TwistedTubeObject*)shape.ptr(), params);   }
0884   template <> void set_dimensions(EllipticalTube shape, const std::vector<double>& params)
0885   {  set_dimensions(shape.ptr(), params);   }
0886   template <> void set_dimensions(Trap shape, const std::vector<double>& params)
0887   {  set_dimensions(shape.ptr(), params);   }
0888   template <> void set_dimensions(Trd1 shape, const std::vector<double>& params)
0889   {  set_dimensions(shape.ptr(), params);   }
0890   template <> void set_dimensions(Trd2 shape, const std::vector<double>& params)
0891   {  set_dimensions(shape.ptr(), params);   }
0892   template <> void set_dimensions(Paraboloid shape, const std::vector<double>& params)
0893   {  set_dimensions(shape.ptr(), params);   }
0894   template <> void set_dimensions(Hyperboloid shape, const std::vector<double>& params)
0895   {  set_dimensions(shape.ptr(), params);   }
0896   template <> void set_dimensions(Sphere shape, const std::vector<double>& params)
0897   {  set_dimensions(shape.ptr(), params);   }
0898   template <> void set_dimensions(Torus shape, const std::vector<double>& params)
0899   {  set_dimensions(shape.ptr(), params);   }
0900   template <> void set_dimensions(PolyhedraRegular shape, const std::vector<double>& params)
0901   {  set_dimensions(shape.ptr(), params);   }
0902   template <> void set_dimensions(Polyhedra shape, const std::vector<double>& params)
0903   {  set_dimensions(shape.ptr(), params);   }
0904   template <> void set_dimensions(ExtrudedPolygon shape, const std::vector<double>& params)
0905   {  set_dimensions(shape.ptr(), params);   }
0906   template <> void set_dimensions(EightPointSolid shape, const std::vector<double>& params)
0907   {  set_dimensions(shape.ptr(), params);   }
0908   template <> void set_dimensions(TessellatedSolid shape, const std::vector<double>& params)
0909   {  set_dimensions(shape.ptr(), params);   }
0910   template <> void set_dimensions(BooleanSolid shape, const std::vector<double>& params)
0911   {  set_dimensions(shape.ptr(), params);   }
0912   template <> void set_dimensions(SubtractionSolid shape, const std::vector<double>& params)
0913   {  set_dimensions(shape.ptr(), params);   }
0914   template <> void set_dimensions(UnionSolid shape, const std::vector<double>& params)
0915   {  set_dimensions(shape.ptr(), params);   }
0916   template <> void set_dimensions(IntersectionSolid shape, const std::vector<double>& params)
0917   {  set_dimensions(shape.ptr(), params);   }
0918 
0919   template <> void set_dimensions(TruncatedTube shape, const std::vector<double>& params)  {
0920     TGeoCompositeShape* sh = (TGeoCompositeShape*) shape.ptr();
0921     TGeoBoolNode* boolean = sh->GetBoolNode();
0922     TGeoMatrix* right_matrix  = boolean->GetRightMatrix();
0923     TGeoShape*  left_solid    = boolean->GetLeftShape();
0924     TGeoShape*  right_solid   = boolean->GetRightShape();
0925     TGeoTubeSeg* tubs = (TGeoTubeSeg*)left_solid;
0926     TGeoBBox*    box  = (TGeoBBox*)right_solid;
0927     double dz    = params[0];
0928     double rmin  = params[1];
0929     double rmax  = params[2];
0930     double startPhi = params[3]/units::deg;
0931     double deltaPhi = params[4]/units::deg;
0932     double cutAtStart = params[5];
0933     double cutAtDelta = params[6];
0934     bool   cutInside  = params[7] > 0.5;
0935 #ifdef DIMENSION_DEBUG
0936     cout << "setDimensions: [" << TRUNCATEDTUBE_TAG << "]" << endl
0937          << right_matrix->GetTitle() << endl;
0938 #endif
0939     // check the parameters
0940     if( rmin <= 0 || rmax <= 0 || cutAtStart <= 0 || cutAtDelta <= 0 )
0941       except(TRUNCATEDTUBE_TAG,"++ 0 <= rIn,cutAtStart,rOut,cutAtDelta,rOut violated!");
0942     else if( rmin >= rmax )
0943       except(TRUNCATEDTUBE_TAG,"++ rIn<rOut violated!");
0944     else if( startPhi != 0. )
0945       except(TRUNCATEDTUBE_TAG,"++ startPhi != 0 not supported!");
0946 
0947     double r         = cutAtStart;
0948     double R         = cutAtDelta;
0949     // angle of the box w.r.t. tubs
0950     double cath      = r - R * std::cos( deltaPhi*units::deg );
0951     double hypo      = std::sqrt( r*r + R*R - 2.*r*R * std::cos( deltaPhi*units::deg ));
0952     double cos_alpha = cath / hypo;
0953     double alpha     = std::acos( cos_alpha );
0954     double sin_alpha = std::sin( std::fabs(alpha) );
0955   
0956     // exaggerate dimensions - does not matter, it's subtracted!
0957     // If we don't, the **edge** of the box would cut into the tube segment
0958     // for larger delta-phi values
0959     double boxX      = 1.1*rmax + rmax/sin_alpha; // Need to adjust for move!
0960     double boxY      = rmax;
0961     // width of the box > width of the tubs
0962     double boxZ      = 1.1 * dz;
0963     double xBox;      // center point of the box
0964     if( cutInside )
0965       xBox = r - boxY / sin_alpha;
0966     else
0967       xBox = r + boxY / sin_alpha;
0968 
0969     // rotationmatrix of box w.r.t. tubs
0970     TGeoRotation rot;
0971     rot.RotateZ( -alpha/units::deg );
0972     double box_dim[] = {boxX, boxY, boxZ};
0973     double tub_dim[] = {rmin, rmax, dz, startPhi, deltaPhi};
0974     box->SetDimensions(box_dim);
0975     tubs->SetDimensions(tub_dim);
0976     TGeoCombiTrans* combi = (TGeoCombiTrans*)right_matrix;
0977     combi->SetRotation(rot);
0978     combi->SetTranslation(xBox, 0, 0);
0979     return;
0980   }
0981   template <> void set_dimensions(PseudoTrap shape, const std::vector<double>& params)  {
0982     TGeoCompositeShape* sh = (TGeoCompositeShape*) shape.ptr();
0983     TGeoBoolNode* boolean = sh->GetBoolNode();
0984     TGeoMatrix* right_matrix  = boolean->GetRightMatrix();
0985     TGeoShape*  left_solid    = boolean->GetLeftShape();
0986     TGeoShape*  right_solid   = boolean->GetRightShape();
0987     double x1 = params[0];
0988     double x2 = params[1];
0989     double y1 = params[2];
0990     double y2 = params[3];
0991     double z  = params[4];
0992     double r  = params[5];
0993     bool   atMinusZ = params[6] > 0.5;
0994     double x            = atMinusZ ? x1 : x2;
0995     double h            = 0;
0996     bool   intersec     = false; // union or intersection solid
0997     double displacement = 0;
0998     double startPhi     = 0;
0999     double halfZ        = z;
1000     double halfOpeningAngle = std::asin( x / std::abs( r ))/units::deg;
1001     /// calculate the displacement of the tubs w.r.t. to the trap, determine the opening angle of the tubs
1002     double delta        = std::sqrt( r * r - x * x );
1003 
1004 #ifdef DIMENSION_DEBUG
1005     cout << "setDimensions: [" << PSEUDOTRAP_TAG << "]" << endl
1006          << right_matrix->GetTitle() << endl;
1007 #endif
1008     // Implementation from :
1009     // https://cmssdt.cern.ch/lxr/source/Fireworks/Geometry/src/TGeoMgrFromDdd.cc#0538
1010     if( r < 0 && std::abs(r) >= x )  {
1011       intersec = true;       // intersection solid
1012       h = y1 < y2 ? y2 : y1; // tubs half height
1013       h += h/20.;            // enlarge a bit - for subtraction solid
1014       if( atMinusZ )    {
1015         displacement = - halfZ - delta; 
1016         startPhi     =  90.0 - halfOpeningAngle;
1017       }
1018       else    {
1019         displacement =   halfZ + delta;
1020         startPhi     = -90.0 - halfOpeningAngle;
1021       }
1022     }
1023     else if( r > 0 && std::abs(r) >= x )  {
1024       if( atMinusZ )    {
1025         displacement = - halfZ + delta;
1026         startPhi     = 270.0 - halfOpeningAngle;
1027         h = y1;
1028       }
1029       else
1030       {
1031         displacement =   halfZ - delta; 
1032         startPhi     =  90.0 - halfOpeningAngle;
1033         h = y2;
1034       }    
1035     }
1036     else  {
1037       except(PSEUDOTRAP_TAG,"Check parameters of the PseudoTrap!");   
1038     }
1039     double trd2_dim[] = { x1, x2, y1, y2, halfZ };
1040     double tube_dim[] = { 0.0, std::abs(r), h, startPhi, startPhi + halfOpeningAngle*2. };
1041     if( intersec && left_solid->IsA() == TGeoTrd2::Class() )   {
1042       left_solid->SetDimensions(trd2_dim);
1043       right_solid->SetDimensions(tube_dim);
1044     }
1045     else if ( right_solid->IsA() == TGeoCompositeShape::Class() )   {
1046       double box_dim[] = { 1.1*x, 1.1*h, std::sqrt(r*r-x*x) };
1047       TGeoCompositeShape* comp = (TGeoCompositeShape*)right_solid;
1048       TGeoBoolNode* b_s = comp->GetBoolNode();
1049       b_s->GetLeftShape()->SetDimensions(tube_dim);
1050       b_s->GetRightShape()->SetDimensions(box_dim);
1051       left_solid->SetDimensions(trd2_dim);
1052     }
1053     else  {
1054       except("PseudoTrap","+++ Incompatible change of parameters.");
1055     }
1056     ((TGeoTranslation*)right_matrix)->SetTranslation(0,0,displacement);
1057     std::stringstream str;
1058     str << x1 << " " << x2 << " " << y1 << " " << y2 << " " << z << " "
1059         << r << " " << char(atMinusZ ? '1' : '0') << " ";
1060     right_matrix->SetTitle(str.str().c_str());
1061   }
1062 
1063   /// Set the shape dimensions. As for the TGeo shape, but angles in rad rather than degrees.
1064   template <> void set_dimensions(Solid shape, const std::vector<double>& params)  {
1065     if ( shape.ptr() ) {
1066       TClass* cl = shape->IsA();
1067       Solid solid(shape);
1068       if ( cl == TGeoShapeAssembly::Class() )
1069         set_dimensions(ShapelessSolid(shape), params);
1070       else if ( cl == TGeoBBox::Class() )
1071         set_dimensions(Box(shape), params);
1072       else if ( cl == TGeoHalfSpace::Class() )
1073         set_dimensions(HalfSpace(shape), params);
1074       else if ( cl == TGeoPcon::Class() )
1075         set_dimensions(Polycone(shape), params);
1076       else if ( cl == TGeoConeSeg::Class() )
1077         set_dimensions(ConeSegment(shape), params);
1078       else if ( cl == TGeoCone::Class() )
1079         set_dimensions(Cone(shape), params);
1080       else if ( cl == TGeoTube::Class() )
1081         set_dimensions(Tube(shape), params);
1082       else if ( cl == TGeoTubeSeg::Class() )
1083         set_dimensions(Tube(shape), params);
1084       else if ( cl == TGeoCtub::Class() )
1085         set_dimensions(CutTube(shape), params);
1086       else if ( cl == TGeoEltu::Class() )
1087         set_dimensions(EllipticalTube(shape), params);
1088       else if ( cl == TwistedTubeObject::Class() )
1089         set_dimensions(TwistedTube(shape), params);
1090       else if ( cl == TGeoTrd1::Class() )
1091         set_dimensions(Trd1(shape), params);
1092       else if ( cl == TGeoTrd2::Class() )
1093         set_dimensions(Trd2(shape), params);
1094       else if ( cl == TGeoParaboloid::Class() )
1095         set_dimensions(Paraboloid(shape), params);
1096       else if ( cl == TGeoHype::Class() )
1097         set_dimensions(Hyperboloid(shape), params);
1098       else if ( cl == TGeoSphere::Class() )
1099         set_dimensions(Sphere(shape), params);
1100       else if ( cl == TGeoTorus::Class() )
1101         set_dimensions(Torus(shape), params);
1102       else if ( cl == TGeoTrap::Class() )
1103         set_dimensions(Trap(shape), params);
1104       else if ( cl == TGeoPgon::Class() )
1105         set_dimensions(PolyhedraRegular(shape), params);
1106       else if ( cl == TGeoXtru::Class() )
1107         set_dimensions(ExtrudedPolygon(shape), params);
1108       else if ( cl == TGeoArb8::Class() )
1109         set_dimensions(EightPointSolid(shape), params);
1110       else if ( cl == TGeoTessellated::Class() )
1111         set_dimensions(TessellatedSolid(shape), params);
1112       else if ( cl == TGeoScaledShape::Class() )  {
1113         TGeoScaledShape* sh = (TGeoScaledShape*) shape.ptr();
1114         set_dimensions(sh, params);
1115       }
1116       else if ( isA<TruncatedTube>(shape) )
1117         set_dimensions(TruncatedTube(shape), params);
1118       else if ( isA<PseudoTrap>(shape) )
1119         set_dimensions(PseudoTrap(shape), params);
1120       else if ( isA<SubtractionSolid>(solid) )
1121         set_dimensions(SubtractionSolid(shape), params);
1122       else if ( isA<UnionSolid>(solid) )
1123         set_dimensions(UnionSolid(shape), params);
1124       else if ( isA<IntersectionSolid>(solid) )
1125         set_dimensions(IntersectionSolid(shape), params);
1126       else if ( isA<BooleanSolid>(solid) )
1127         set_dimensions(BooleanSolid(shape), params);
1128       else
1129         printout(ERROR,"Solid","Failed to set dimensions for shape of type:%s.",cl->GetName() );
1130       return;
1131     }
1132     except("Solid","set_shape_dimensions: Failed to set dimensions [Invalid handle].");
1133   }
1134   /// Set the shape dimensions. As for the TGeo shape, but angles in rad rather than degrees.
1135   template <> void set_dimensions(const TGeoShape* shape, const std::vector<double>& params)  {
1136     set_dimensions<Solid>(Solid(shape), params);
1137   }
1138   /// Set the shape dimensions. As for the TGeo shape, but angles in rad rather than degrees.
1139   void set_shape_dimensions(TGeoShape* shape, const std::vector<double>& params)   {
1140     set_dimensions<Solid>(Solid(shape), params);
1141   }
1142   /// Set the shape dimensions. As for the TGeo shape, but angles in rad rather than degrees.
1143   void set_shape_dimensions(Solid solid, const std::vector<double>& params)   {
1144     set_dimensions<Solid>(solid, params);
1145   }
1146 
1147   /// Pretty print of solid attributes
1148   std::string toStringSolid(const TGeoShape* shape, int precision)   {
1149     if ( !shape )  {
1150       return "[Invalid shape]";
1151     }
1152 
1153     std::stringstream log;
1154     TClass* cl = shape->IsA();
1155     int prec = log.precision();
1156     log.setf(std::ios::fixed,std::ios::floatfield);
1157     log << std::setprecision(precision);
1158     log << cl->GetName();
1159     if ( cl == TGeoBBox::Class() )   {
1160       TGeoBBox* sh = (TGeoBBox*) shape;
1161       log << " x:" << sh->GetDX()
1162           << " y:" << sh->GetDY()
1163           << " z:" << sh->GetDZ();
1164     }
1165     else if ( cl == TGeoHalfSpace::Class() ) {
1166       TGeoHalfSpace* sh = (TGeoHalfSpace*)(const_cast<TGeoShape*>(shape));
1167       log << " Point:  (" << sh->GetPoint()[0] << ", " << sh->GetPoint()[1] << ", " << sh->GetPoint()[2] << ") " 
1168           << " Normal: (" << sh->GetNorm()[0]  << ", " << sh->GetNorm()[1]  << ", " << sh->GetNorm()[2]  << ") ";
1169     }
1170     else if ( cl == TGeoTube::Class() ) {
1171       const TGeoTube* sh = (const TGeoTube*) shape;
1172       log << " rmin:" << sh->GetRmin() << " rmax:" << sh->GetRmax() << " dz:" << sh->GetDz();
1173     }
1174     else if ( cl == TGeoTubeSeg::Class() ) {
1175       const TGeoTubeSeg* sh = (const TGeoTubeSeg*) shape;
1176       log << " rmin:" << sh->GetRmin() << " rmax:" << sh->GetRmax() << " dz:" << sh->GetDz()
1177           << " Phi1:" << sh->GetPhi1() << " Phi2:" << sh->GetPhi2();
1178     }
1179     else if ( cl == TGeoCtub::Class() ) {
1180       const TGeoCtub* sh = (const TGeoCtub*) shape;
1181       const Double_t*   hi = sh->GetNhigh();
1182       const Double_t*   lo = sh->GetNlow();
1183       log << " rmin:" << sh->GetRmin() << " rmax:" << sh->GetRmax() << " dz:" << sh->GetDz()
1184           << " Phi1:" << sh->GetPhi1() << " Phi2:" << sh->GetPhi2()
1185           << " lx:" << lo[0] << " ly:" << lo[1] << " lz:" << lo[2]
1186           << " hx:" << hi[0] << " hy:" << hi[1] << " hz:" << hi[2];
1187     }
1188     else if ( cl == TGeoEltu::Class() ) {
1189       const TGeoEltu* sh = (const TGeoEltu*) shape;
1190       log << " A:" << sh->GetA() << " B:" << sh->GetB() << " dz:" << sh->GetDz();
1191     }
1192     else if ( cl == TGeoTrd1::Class() ) {
1193       const TGeoTrd1* sh = (const TGeoTrd1*) shape;
1194       log << " x1:" << sh->GetDx1() << " x2:" << sh->GetDx2() << " y:" << sh->GetDy() << " z:" << sh->GetDz();
1195     }
1196     else if ( cl == TGeoTrd2::Class() ) {
1197       const TGeoTrd2* sh = (const TGeoTrd2*) shape;
1198       log << " x1:" << sh->GetDx1() << " x2:" << sh->GetDx2()
1199           << " y1:" << sh->GetDy1() << " y2:" << sh->GetDy2() << " z:" << sh->GetDz();
1200     }
1201     else if ( cl == TGeoTrap::Class() )   {
1202       const TGeoTrap* sh = (const TGeoTrap*) shape;
1203       log << " dz:" << sh->GetDz() << " Theta:" << sh->GetTheta() << " Phi:" << sh->GetPhi()
1204           << " H1:" << sh->GetH1() << " Bl1:"   << sh->GetBl1()   << " Tl1:" << sh->GetTl1() << " Alpha1:" << sh->GetAlpha1()
1205           << " H2:" << sh->GetH2() << " Bl2:"   << sh->GetBl2()   << " Tl2:" << sh->GetTl2() << " Alpha2:" << sh->GetAlpha2();
1206     }
1207     else if ( cl == TGeoHype::Class() )   {
1208       const TGeoHype* sh = (const TGeoHype*) shape;
1209       log << " rmin:" << sh->GetRmin() << " rmax:"  << sh->GetRmax() << " dz:" << sh->GetDz()
1210           << " StIn:" << sh->GetStIn() << " StOut:" << sh->GetStOut();
1211     }
1212     else if ( cl == TGeoPgon::Class() )   {
1213       const TGeoPgon* sh = (const TGeoPgon*) shape;
1214       log << " Phi1:"   << sh->GetPhi1()   << " dPhi:" << sh->GetDphi()
1215           << " NEdges:" << sh->GetNedges() << " Nz:" << sh->GetNz();
1216       for(int i=0, n=sh->GetNz(); i<n; ++i)  {
1217         log << " i=" << i << " z:" << sh->GetZ(i)
1218             << " r:[" << sh->GetRmin(i) << "," << sh->GetRmax(i) << "]";
1219       }
1220     }
1221     else if ( cl == TGeoPcon::Class() )  {
1222       const TGeoPcon* sh = (const TGeoPcon*) shape;
1223       log << " Phi1:" << sh->GetPhi1() << " dPhi:" << sh->GetDphi() << " Nz:" << sh->GetNz();
1224       for(int i=0, n=sh->GetNz(); i<n; ++i)  {
1225         log << " i=" << i << " z:" << sh->GetZ(i)
1226             << " r:[" << sh->GetRmin(i) << "," << sh->GetRmax(i) << "]";
1227       }
1228     }
1229     else if ( cl == TGeoCone::Class() )  {
1230       const TGeoCone* sh = (const TGeoCone*) shape;
1231       log << " rmin1:" << sh->GetRmin1() << " rmax1:" << sh->GetRmax1()
1232           << " rmin2:" << sh->GetRmin2() << " rmax2:" << sh->GetRmax2()
1233           << " dz:"    << sh->GetDz();
1234     }
1235     else if ( cl == TGeoConeSeg::Class() )  {
1236       const TGeoConeSeg* sh = (const TGeoConeSeg*) shape;
1237       log << " rmin1:" << sh->GetRmin1() << " rmax1:" << sh->GetRmax1()
1238           << " rmin2:" << sh->GetRmin2() << " rmax2:" << sh->GetRmax2()
1239           << " dz:"    << sh->GetDz()
1240           << " Phi1:"  << sh->GetPhi1() << " Phi2:" << sh->GetPhi2();
1241     }
1242     else if ( cl == TGeoParaboloid::Class() )  {
1243       const TGeoParaboloid* sh = (const TGeoParaboloid*) shape;
1244       log << " dz:" << sh->GetDz() << " RLo:" << sh->GetRlo() << " Rhi:" << sh->GetRhi();
1245     }
1246     else if ( cl == TGeoSphere::Class() )   {
1247       const TGeoSphere* sh = (const TGeoSphere*) shape;
1248       log << " rmin:" << sh->GetRmin() << " rmax:" << sh->GetRmax()
1249           << " Phi1:" << sh->GetPhi1() << " Phi2:" << sh->GetPhi2()
1250           << " Theta1:" << sh->GetTheta1() << " Theta2:" << sh->GetTheta2();
1251     }
1252     else if ( cl == TGeoTorus::Class() )   {
1253       const TGeoTorus* sh = (const TGeoTorus*) shape;
1254       log << " rmin:" << sh->GetRmin() << " rmax:" << sh->GetRmax() << " r:" << sh->GetR()
1255           << " Phi1:" << sh->GetPhi1() << " dPhi:" << sh->GetDphi();
1256     }
1257     else if ( cl == TGeoArb8::Class() )   {
1258       TGeoArb8* sh = (TGeoArb8*) shape;
1259       const Double_t* v = sh->GetVertices();
1260       log << " dz:" << sh->GetDz();
1261       for(int i=0; i<8; ++i) {
1262         log << " x[" << i << "]:" << *v; ++v;
1263         log << " y[" << i << "]:" << *v; ++v;
1264       }
1265     }
1266     else if ( cl == TGeoXtru::Class() )   {
1267       const TGeoXtru* sh = (const TGeoXtru*) shape;
1268       log << " X:   " << sh->GetX(0)
1269           << " Y:   " << sh->GetY(0)
1270           << " Z:   " << sh->GetZ(0)
1271           << " Xoff:" << sh->GetXOffset(0)
1272           << " Yoff:" << sh->GetYOffset(0)
1273           << " Nvtx:" << sh->GetNvert()
1274           << " Nz:  " << sh->GetNz();
1275     }
1276     else if (shape->IsA() == TGeoCompositeShape::Class() )   {
1277       const TGeoCompositeShape* sh = (const TGeoCompositeShape*) shape;
1278       const TGeoBoolNode* boolean = sh->GetBoolNode();
1279       const TGeoShape* left  = boolean->GetLeftShape();
1280       const TGeoShape* right = boolean->GetRightShape();
1281       TGeoBoolNode::EGeoBoolType oper = boolean->GetBooleanOperator();
1282       if (oper == TGeoBoolNode::kGeoSubtraction)
1283         log << "Subtraction: ";
1284       else if (oper == TGeoBoolNode::kGeoUnion)
1285         log << "Union: ";
1286       else if (oper == TGeoBoolNode::kGeoIntersection)
1287         log << "Intersection: ";
1288       log << " Left:" << toStringSolid(left) << " Right:" << toStringSolid(right);
1289     }
1290     log << std::setprecision(prec);
1291     return log.str();
1292   }
1293 
1294   /// Output mesh vertices to string
1295   std::string toStringMesh(const TGeoShape* shape, int prec)  {
1296     Solid sol(shape);
1297     std::stringstream os;
1298 
1299     // Prints shape parameters
1300     Int_t nvert = 0, nsegs = 0, npols = 0;
1301 
1302     sol->GetMeshNumbers(nvert, nsegs, npols);
1303     Double_t* points = new Double_t[3*nvert];
1304     sol->SetPoints(points);
1305 
1306     os << std::setw(16) << std::left << sol->IsA()->GetName()
1307        << " " << nvert << " Mesh-points:" << std::endl;
1308     os << std::setw(16) << std::left << sol->IsA()->GetName() << " " << sol->GetName()
1309        << " N(mesh)=" << sol->GetNmeshVertices()
1310        << "  N(vert)=" << nvert << "  N(seg)=" << nsegs << "  N(pols)=" << npols << std::endl;
1311     
1312     for(Int_t i=0; i<nvert; ++i)   {
1313       Double_t* p = points + 3*i;
1314       os << std::setw(16) << std::left << sol->IsA()->GetName() << " " << std::setw(3) << std::left << i
1315          << " Local  ("  << std::setw(7) << std::setprecision(prec) << std::fixed << std::right << p[0]/dd4hep::cm
1316          << ", "         << std::setw(7) << std::setprecision(prec) << std::fixed << std::right << p[1]/dd4hep::cm
1317          << ", "         << std::setw(7) << std::setprecision(prec) << std::fixed << std::right << p[2]/dd4hep::cm
1318          << ")" << std::endl;
1319     }
1320     Box box = sol;
1321     const Double_t* org = box->GetOrigin();
1322     os << std::setw(16) << std::left << sol->IsA()->GetName()
1323        << " Bounding box: "
1324        << " dx="        << std::setw(7) << std::setprecision(prec) << std::fixed << std::right << box->GetDX()/dd4hep::cm
1325        << " dy="        << std::setw(7) << std::setprecision(prec) << std::fixed << std::right << box->GetDY()/dd4hep::cm
1326        << " dz="        << std::setw(7) << std::setprecision(prec) << std::fixed << std::right << box->GetDZ()/dd4hep::cm
1327        << " Origin: x=" << std::setw(7) << std::setprecision(prec) << std::fixed << std::right << org[0]/dd4hep::cm
1328        << " y="         << std::setw(7) << std::setprecision(prec) << std::fixed << std::right << org[1]/dd4hep::cm
1329        << " z="         << std::setw(7) << std::setprecision(prec) << std::fixed << std::right << org[2]/dd4hep::cm
1330        << std::endl;
1331   
1332     /// -------------------- DONE --------------------
1333     delete [] points;
1334     return os.str();
1335   }
1336 }
1337