File indexing completed on 2025-03-13 08:19:44
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014 #define _USE_MATH_DEFINES
0015
0016
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
0025 #include <stdexcept>
0026 #include <iomanip>
0027 #include <sstream>
0028
0029
0030 #include <TClass.h>
0031 #include <TGeoMatrix.h>
0032 #include <TGeoBoolNode.h>
0033 #include <TGeoScaledShape.h>
0034 #include <TGeoCompositeShape.h>
0035
0036
0037 namespace dd4hep {
0038
0039
0040 namespace units = dd4hep;
0041
0042
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
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
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
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
0221
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
0603 std::vector<double> get_shape_dimensions(const Handle<TGeoShape>& shape) {
0604 return dimensions<Solid>(shape);
0605 }
0606
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
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
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
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
0957
0958
0959 double boxX = 1.1*rmax + rmax/sin_alpha;
0960 double boxY = rmax;
0961
0962 double boxZ = 1.1 * dz;
0963 double xBox;
0964 if( cutInside )
0965 xBox = r - boxY / sin_alpha;
0966 else
0967 xBox = r + boxY / sin_alpha;
0968
0969
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;
0997 double displacement = 0;
0998 double startPhi = 0;
0999 double halfZ = z;
1000 double halfOpeningAngle = std::asin( x / std::abs( r ))/units::deg;
1001
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
1009
1010 if( r < 0 && std::abs(r) >= x ) {
1011 intersec = true;
1012 h = y1 < y2 ? y2 : y1;
1013 h += h/20.;
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
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
1135 template <> void set_dimensions(const TGeoShape* shape, const std::vector<double>& params) {
1136 set_dimensions<Solid>(Solid(shape), params);
1137 }
1138
1139 void set_shape_dimensions(TGeoShape* shape, const std::vector<double>& params) {
1140 set_dimensions<Solid>(Solid(shape), params);
1141 }
1142
1143 void set_shape_dimensions(Solid solid, const std::vector<double>& params) {
1144 set_dimensions<Solid>(solid, params);
1145 }
1146
1147
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
1295 std::string toStringMesh(const TGeoShape* shape, int prec) {
1296 Solid sol(shape);
1297 std::stringstream os;
1298
1299
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
1333 delete [] points;
1334 return os.str();
1335 }
1336 }
1337