Geometry/test/testVolumes.cpp

0001 #include <DetectorDescription/Core/interface/DDSolid.h>
0002 #include <DetectorDescription/Core/interface/DDSolidShapes.h>
0003
0004 #include <DetectorDescription/Core/interface/Box.h>
0005 #include <DetectorDescription/Core/interface/Cons.h>
0006 #include <DetectorDescription/Core/interface/CutTubs.h>
0007 #include <DetectorDescription/Core/interface/EllipticalTube.h>
0008 #include <DetectorDescription/Core/interface/ExtrudedPolygon.h>
0009 #include <DetectorDescription/Core/interface/Polycone.h>
0010 #include <DetectorDescription/Core/interface/Polyhedra.h>
0011 #include <DetectorDescription/Core/interface/Sphere.h>
0012 #include <DetectorDescription/Core/interface/Torus.h>
0013 #include <DetectorDescription/Core/interface/Trap.h>
0014 #include <DetectorDescription/Core/interface/Tubs.h>
0015
0016 #include <CLHEP/Units/SystemOfUnits.h>
0017 #include <DataFormats/GeometryVector/interface/Pi.h>
0018 #include <G4Box.hh>
0019 #include <G4Cons.hh>
0020 #include <G4CutTubs.hh>
0021 #include <G4Ellipsoid.hh>
0022 #include <G4EllipticalTube.hh>
0023 #include <G4ExtrudedSolid.hh>
0024 #include <G4Orb.hh>
0025 #include <G4Para.hh>
0026 #include <G4Polycone.hh>
0027 #include <G4Polyhedra.hh>
0028 #include <G4Sphere.hh>
0029 #include <G4Torus.hh>
0030 #include <G4Trap.hh>
0031 #include <G4Trd.hh>
0032 #include <G4Tubs.hh>
0033 #include <string>
0034
0035 using CLHEP::cm;
0036 using CLHEP::cm3;
0037 using CLHEP::deg;
0038 //
0039 // See Geant4 documentation for more details:
0040 //
0041 // http://geant4.web.cern.ch/geant4/UserDocumentation/UsersGuides/ForApplicationDeveloper/html/ch04.html#sect.Geom.Solids
0042 //
0043 //
0044 // This test verifies convertion of the DDD solids to Geant4 Constructed Solid
0045 // Geometry (CSG) Solids
0046 //
0047
0048 //
0049 // 1. Box:
0050 //
0051 // G4Box(const G4String& pName,
0052 //       G4double  pX,
0053 //       G4double  pY,
0054 //       G4double  pZ)
0055 void doBox(const std::string &name, double xHalfLength, double yHalfLength, double zHalfLength) {
0056   G4Box g4(name, xHalfLength, yHalfLength, zHalfLength);
0057   DDI::Box dd(xHalfLength, yHalfLength, zHalfLength);
0058   DDBox dds = DDSolidFactory::box(name, xHalfLength, yHalfLength, zHalfLength);
0059   dd.stream(std::cout);
0060   std::cout << std::endl;
0061   std::cout << "\tg4 volume = " << g4.GetCubicVolume() / cm3 << " cm3" << std::endl;
0062   std::cout << "\tdd volume = " << dd.volume() / cm3 << " cm3" << std::endl;
0063   std::cout << "\tDD Information: " << dds << " vol= " << dds.volume() << std::endl;
0064 }
0065
0066 //
0067 // 2. Cylindrical Section or Tube:
0068 //
0069 // G4Tubs(const G4String& pName,
0070 //        G4double  pRMin,
0071 //        G4double  pRMax,
0072 //        G4double  pDz,
0073 //        G4double  pSPhi,
0074 //        G4double  pDPhi)
0075 void doTubs(const std::string &name, double rIn, double rOut, double zhalf, double startPhi, double deltaPhi) {
0076   G4Tubs g4(name, rIn, rOut, zhalf, startPhi, deltaPhi);
0077   DDI::Tubs dd(zhalf, rIn, rOut, startPhi, deltaPhi);
0078   DDTubs dds = DDSolidFactory::tubs(name, zhalf, rIn, rOut, startPhi, deltaPhi);
0079   dd.stream(std::cout);
0080   std::cout << std::endl;
0081   std::cout << "\tg4 volume = " << g4.GetCubicVolume() / cm3 << " cm3" << std::endl;
0082   std::cout << "\tdd volume = " << dd.volume() / cm3 << " cm3" << std::endl;
0083   std::cout << "\tDD Information: " << dds << " vol= " << dds.volume() << std::endl;
0084 }
0085
0086 //
0087 // 3. Cone or Conical section:
0088 //
0089 // G4Cons(const G4String& pName,
0090 //        G4double  pRmin1,
0091 //        G4double  pRmax1,
0092 //        G4double  pRmin2,
0093 //        G4double  pRmax2,
0094 //        G4double  pDz,
0095 //        G4double  pSPhi,
0096 //        G4double  pDPhi)
0097 void doCons(const std::string &name,
0098             double rIn1,
0099             double rOut1,
0100             double rIn2,
0101             double rOut2,
0102             double zhalf,
0103             double startPhi,
0104             double deltaPhi) {
0105   G4Cons g4(name, rIn1, rOut1, rIn2, rOut2, zhalf, startPhi, deltaPhi);
0106   DDI::Cons dd(zhalf, rIn1, rOut1, rIn2, rOut2, startPhi, deltaPhi);
0107   DDCons dds = DDSolidFactory::cons(name, zhalf, rIn1, rOut1, rIn2, rOut2, startPhi, deltaPhi);
0108   dd.stream(std::cout);
0109   std::cout << std::endl;
0110   std::cout << "\tg4 volume = " << g4.GetCubicVolume() / cm3 << " cm3" << std::endl;
0111   std::cout << "\tdd volume = " << dd.volume() / cm3 << " cm3" << std::endl;
0112   std::cout << "\tDD Information: " << dds << " vol= " << dds.volume() << std::endl;
0113 }
0114
0115 //
0116 // 5. Trapezoid:
0117 //
0118 // G4Trd(const G4String& pName,
0119 //       G4double  dx1,
0120 //       G4double  dx2,
0121 //       G4double  dy1,
0122 //       G4double  dy2,
0123 //       G4double  dz)
0124 void doTrd(const std::string &name, double dx1, double dx2, double dy1, double dy2, double dz) {
0125   G4Trd g4(name, dx1, dx2, dy1, dy2, dz);
0126   /////////////////////////////////////////
0127   // DDD does not have direct implementation of Trd.
0128   // Use generic trapezoid instead.
0129   DDI::Trap dd(dz, 0.0 /* pTheta */, 0.0 /* pPhi */, dy1, dx1, dx1, 0.0 /* pAlp1 */, dy2, dx2, dx2, 0.0 /* pAlp2 */);
0130   DDTrap dds = DDSolidFactory::trap(
0131       name, dz, 0.0 /* pTheta */, 0.0 /* pPhi */, dy1, dx1, dx1, 0.0 /* pAlp1 */, dy2, dx2, dx2, 0.0 /* pAlp2 */);
0132   dd.stream(std::cout);
0133   std::cout << std::endl;
0134   std::cout << "\tg4 volume = " << g4.GetCubicVolume() / cm3 << " cm3" << std::endl;
0135   std::cout << "\tdd volume = " << dd.volume() / cm3 << " cm3" << std::endl;
0136   std::cout << "\tDD Information: " << dds << " vol= " << dds.volume() << std::endl;
0137 }
0138
0139 //
0140 // 6. Generic Trapezoid:
0141 //
0142 // G4Trap(const G4String& pName,
0143 //        G4double   pZ,
0144 //        G4double   pY,
0145 //        G4double   pX,
0146 //        G4double   pLTX)
0147 //
0148 // G4Trap(const G4String& pName,
0149 //        G4double   pDz,   G4double   pTheta,
0150 //        G4double   pPhi,  G4double   pDy1,
0151 //        G4double   pDx1,  G4double   pDx2,
0152 //        G4double   pAlp1, G4double   pDy2,
0153 //        G4double   pDx3,  G4double   pDx4,
0154 //        G4double   pAlp2)
0155 void doTrap(const std::string &name,
0156             double dz,
0157             double pTheta,
0158             double pPhi,
0159             double pDy1,
0160             double pDx1,
0161             double pDx2,
0162             double pAlp1,
0163             double pDy2,
0164             double pDx3,
0165             double pDx4,
0166             double pAlp2) {
0167   G4Trap g4(name, dz, pTheta, pPhi, pDy1, pDx1, pDx2, pAlp1, pDy2, pDx3, pDx4, pAlp2);
0168
0169   // Note, the order of parameters is different:
0170   DDI::Trap dd(dz, pTheta, pPhi, pDy2, pDx3, pDx4, pAlp1, pDy1, pDx1, pDx2, pAlp2);
0171   DDTrap dds = DDSolidFactory::trap(name, dz, pTheta, pPhi, pDy2, pDx3, pDx4, pAlp1, pDy1, pDx1, pDx2, pAlp2);
0172   dd.stream(std::cout);
0173   std::cout << std::endl;
0174   std::cout << "\tg4 volume = " << g4.GetCubicVolume() / cm3 << " cm3" << std::endl;
0175   std::cout << "\tdd volume = " << dd.volume() / cm3 << " cm3" << std::endl;
0176   std::cout << "\tDD Information: " << dds << " vol= " << dds.volume() << std::endl;
0177 }
0178
0179 //
0180 // 7. Sphere or Spherical Shell Section:
0181 //
0182 // G4Sphere(const G4String& pName,
0183 //      G4double   pRmin,
0184 //      G4double   pRmax,
0185 //      G4double   pSPhi,
0186 //      G4double   pDPhi,
0187 //      G4double   pSTheta,
0188 //      G4double   pDTheta )
0189 void doSphere(const std::string &name,
0190               double innerRadius,
0191               double outerRadius,
0192               double startPhi,
0193               double deltaPhi,
0194               double startTheta,
0195               double deltaTheta) {
0196   G4Sphere g4(name, innerRadius, outerRadius, startPhi, deltaPhi, startTheta, deltaTheta);
0197   DDI::Sphere dd(innerRadius, outerRadius, startPhi, deltaPhi, startTheta, deltaTheta);
0198   DDSphere dds = DDSolidFactory::sphere(name, innerRadius, outerRadius, startPhi, deltaPhi, startTheta, deltaTheta);
0199   dd.stream(std::cout);
0200   std::cout << std::endl;
0201   std::cout << "\tg4 volume = " << g4.GetCubicVolume() / cm3 << " cm3" << std::endl;
0202   std::cout << "\tdd volume = " << dd.volume() / cm3 << " cm3" << std::endl;
0203   std::cout << "\tDD Information: " << dds << " vol= " << dds.volume() << std::endl;
0204 }
0205
0206 //
0207 // 9. Torus:
0208 //
0209 // G4Torus(const G4String& pName,
0210 //     G4double   pRmin,
0211 //     G4double   pRmax,
0212 //     G4double   pRtor,
0213 //     G4double   pSPhi,
0214 //     G4double   pDPhi)
0215 void doTorus(const std::string &name, double rMin, double rMax, double radius, double sPhi, double dPhi) {
0216   G4Torus g4(name, rMin, rMax, radius, sPhi, dPhi);
0217   DDI::Torus dd(rMin, rMax, radius, sPhi, dPhi);
0218   DDTorus dds = DDSolidFactory::torus(name, rMin, rMax, radius, sPhi, dPhi);
0219   dd.stream(std::cout);
0220   std::cout << std::endl;
0221   std::cout << "\tg4 volume = " << g4.GetCubicVolume() / cm3 << " cm3" << std::endl;
0222   std::cout << "\tdd volume = " << dd.volume() / cm3 << " cm3" << std::endl;
0223   std::cout << "\tDD Information: " << dds << " vol= " << dds.volume() << std::endl;
0224 }
0225
0226 // Specific CSG Solids
0227 //
0228
0229 //
0230 // 10. Polycons:
0231 //
0232 // G4Polycone(const G4String& pName,
0233 //     G4double   phiStart,
0234 //     G4double   phiTotal,
0235 //     G4int      numZPlanes,
0236 //     const G4double   zPlane[],
0237 //     const G4double   rInner[],
0238 //     const G4double   rOuter[])
0239 //
0240 // G4Polycone(const G4String& pName,
0241 //     G4double   phiStart,
0242 //     G4double   phiTotal,
0243 //     G4int      numRZ,
0244 //     const G4double  r[],
0245 //     const G4double  z[])
0246 void doPolycone1(const std::string &name,
0247                  double phiStart,
0248                  double phiTotal,
0249                  const std::vector<double> &z,
0250                  const std::vector<double> &rInner,
0251                  const std::vector<double> &rOuter) {
0252   G4Polycone g4(name, phiStart, phiTotal, z.size(), &z[0], &rInner[0], &rOuter[0]);
0253   DDI::Polycone dd(phiStart, phiTotal, z, rInner, rOuter);
0254   DDPolycone dds = DDSolidFactory::polycone(name, phiStart, phiTotal, z, rInner, rOuter);
0255   dd.stream(std::cout);
0256   std::cout << std::endl;
0257   std::cout << "\tg4 volume = " << g4.GetCubicVolume() / cm3 << " cm3" << std::endl;
0258   std::cout << "\tdd volume = " << dd.volume() / cm3 << " cm3" << std::endl;
0259   std::cout << "\tDD Information: " << dds << " vol= " << dds.volume() << std::endl;
0260 }
0261
0262 void doPolycone2(const std::string &name,
0263                  double phiStart,
0264                  double phiTotal,
0265                  const std::vector<double> &z,
0266                  const std::vector<double> &r) {
0267   std::cout << "### doPolycone_RZ: "
0268             << "phi1=" << phiStart / deg << " phi2=" << phiTotal / deg << " N= " << z.size() << std::endl;
0269   for (size_t i = 0; i < z.size(); ++i) {
0270     std::cout << " R= " << r[i] << " Z= " << z[i] << std::endl;
0271   }
0272   G4Polycone g4(name, phiStart, phiTotal, z.size(), &r[0], &z[0]);
0273   DDI::Polycone dd(phiStart, phiTotal, z, r);
0274   DDPolycone dds = DDSolidFactory::polycone(name, phiStart, phiTotal, z, r);
0275   dd.stream(std::cout);
0276   std::cout << std::endl;
0277   std::cout << "\tg4 volume = " << g4.GetCubicVolume() / cm3 << " cm3" << std::endl;
0278   std::cout << "\tdd volume = " << dd.volume() / cm3 << " cm3" << std::endl;
0279   std::cout << "\tDD Information: " << dds << " vol= " << dds.volume() << std::endl;
0280 }
0281
0282 //
0283 // 11. Polyhedra (PGON):
0284 //
0285 // G4Polyhedra(const G4String& pName,
0286 //      G4double  phiStart,
0287 //      G4double  phiTotal,
0288 //      G4int     numSide,
0289 //      G4int     numZPlanes,
0290 //      const G4double  zPlane[],
0291 //      const G4double  rInner[],
0292 //      const G4double  rOuter[] )
0293 //
0294 // G4Polyhedra(const G4String& pName,
0295 //      G4double  phiStart,
0296 //      G4double  phiTotal,
0297 //      G4int     numSide,
0298 //      G4int     numRZ,
0299 //      const G4double  r[],
0300 //      const G4double  z[] )
0301 void doPolyhedra1(const std::string &name,
0302                   int sides,
0303                   double phiStart,
0304                   double phiTotal,
0305                   const std::vector<double> &z,
0306                   const std::vector<double> &rInner,
0307                   const std::vector<double> &rOuter) {
0308   G4Polyhedra g4(name, phiStart, phiTotal, sides, z.size(), &z[0], &rInner[0], &rOuter[0]);
0309   DDI::Polyhedra dd(sides, phiStart, phiTotal, z, rInner, rOuter);
0310   DDPolyhedra dds = DDSolidFactory::polyhedra(name, sides, phiStart, phiTotal, z, rInner, rOuter);
0311   dd.stream(std::cout);
0312   std::cout << std::endl;
0313   std::cout << "\tg4 volume = " << g4.GetCubicVolume() / cm3 << " cm3" << std::endl;
0314   std::cout << "\tdd volume = " << dd.volume() / cm3 << " cm3" << std::endl;
0315   std::cout << "\tDD Information: " << dds << " vol= " << dds.volume() << std::endl;
0316 }
0317
0318 void doPolyhedra2(const std::string &name,
0319                   int sides,
0320                   double phiStart,
0321                   double phiTotal,
0322                   const std::vector<double> &z,
0323                   const std::vector<double> &r) {
0324   G4Polyhedra g4(name, phiStart, phiTotal, sides, z.size(), &r[0], &z[0]);
0325   DDI::Polyhedra dd(sides, phiStart, phiTotal, z, r);
0326   DDPolyhedra dds = DDSolidFactory::polyhedra(name, sides, phiStart, phiTotal, z, r);
0327   dd.stream(std::cout);
0328   std::cout << std::endl;
0329   std::cout << "\tg4 volume = " << g4.GetCubicVolume() / cm3 << " cm3" << std::endl;
0330   std::cout << "\tdd volume = " << dd.volume() / cm3 << " cm3" << std::endl;
0331   std::cout << "\tDD Information: " << dds << " vol= " << dds.volume() << std::endl;
0332 }
0333
0334 //
0335 // 12. Tube with an elliptical cross section:
0336 //
0337 // G4EllipticalTube(const G4String& pName,
0338 //       G4double  Dx,
0339 //       G4double  Dy,
0340 //       G4double  Dz)
0341 void doEllipticalTube(const std::string &name, double xSemiaxis, double ySemiAxis, double zHeight) {
0342   G4EllipticalTube g4t(name, xSemiaxis, ySemiAxis, zHeight);
0343   DDI::EllipticalTube ddt(xSemiaxis, ySemiAxis, zHeight);
0344   DDEllipticalTube ddet = DDSolidFactory::ellipticalTube(name, xSemiaxis, ySemiAxis, zHeight);
0345   ddt.stream(std::cout);
0346   std::cout << std::endl;
0347   std::cout << "\tg4 volume = " << g4t.GetCubicVolume() / cm3 << " cm3" << std::endl;
0348   std::cout << "\tdd volume = " << ddt.volume() / cm3 << " cm3" << std::endl;
0349   std::cout << "\tcalc volume = " << 2 * zHeight * Geom::pi() * ySemiAxis * xSemiaxis / cm3 << " cm3 " << std::endl;
0350   std::cout << "\tDD Information: ";
0351   std::cout << ddet << " vol= " << ddet.volume() << std::endl;
0352 }
0353
0354 //
0355 // 14. Cone with Elliptical Cross Section:
0356 //
0357 // G4EllipticalCone(const G4String& pName,
0358 //       G4double  pxSemiAxis,
0359 //       G4double  pySemiAxis,
0360 //       G4double  zMax,
0361 //       G4double  pzTopCut)
0362
0363 //
0364 // 15. Paraboloid, a solid with parabolic profile:
0365 //
0366 // G4Paraboloid(const G4String& pName,
0367 //       G4double  Dz,
0368 //       G4double  R1,
0369 //       G4double  R2)
0370
0371 //
0372 // 16. Tube with Hyperbolic Profile:
0373 //
0374 // G4Hype(const G4String& pName,
0375 //        G4double  innerRadius,
0376 //        G4double  outerRadius,
0377 //        G4double  innerStereo,
0378 //        G4double  outerStereo,
0379 //        G4double  halfLenZ)
0380
0381 //
0382 // 17. Tetrahedra:
0383 //
0384 // G4Tet(const G4String& pName,
0385 //       G4ThreeVector  anchor,
0386 //       G4ThreeVector  p2,
0387 //       G4ThreeVector  p3,
0388 //       G4ThreeVector  p4,
0389 //       G4bool         *degeneracyFlag=0)
0390
0391 //
0392 // 18. Extruded Polygon:
0393 //
0394 // G4ExtrudedSolid(const G4String&                pName,
0395 //      std::vector<G4TwoVector> polygon,
0396 //      std::vector<ZSection>    zsections)
0397 //
0398 // G4ExtrudedSolid(const G4String&                pName,
0399 //      std::vector<G4TwoVector> polygon,
0400 //      G4double                 hz,
0401 //      G4TwoVector off1, G4double scale1,
0402 //      G4TwoVector off2, G4double scale2)
0403
0404 //
0405 // 19. Box Twisted:
0406 //
0407 // G4TwistedBox(const G4String& pName,
0408 //       G4double  twistedangle,
0409 //       G4double  pDx,
0410 //       G4double  pDy,
0411 //       G4double  pDz)
0412
0413 //
0414 // 20. Trapezoid Twisted along One Axis:
0415 //
0416 // G4TwistedTrap(const G4String& pName,
0417 //        G4double  twistedangle,
0418 //        G4double  pDxx1,
0419 //        G4double  pDxx2,
0420 //        G4double  pDy,
0421 //        G4double   pDz)
0422 //
0423 // G4TwistedTrap(const G4String& pName,
0424 //        G4double  twistedangle,
0425 //        G4double  pDz,
0426 //        G4double  pTheta,
0427 //        G4double  pPhi,
0428 //        G4double  pDy1,
0429 //        G4double  pDx1,
0430 //        G4double  pDx2,
0431 //        G4double  pDy2,
0432 //        G4double  pDx3,
0433 //        G4double  pDx4,
0434 //        G4double  pAlph)
0435
0436 //
0437 // 21. Twisted Trapezoid with x and y dimensions varying along z:
0438 //
0439 // G4TwistedTrd(const G4String& pName,
0440 //       G4double  pDx1,
0441 //       G4double  pDx2,
0442 //       G4double  pDy1,
0443 //       G4double  pDy2,
0444 //       G4double  pDz,
0445 //       G4double  twistedangle)
0446
0447 //
0448 // 22. Generic trapezoid with optionally collapsing vertices:
0449 //
0450 // G4GenericTrap(const G4String& pName,
0451 //        G4double  pDz,
0452 //        const std::vector<G4TwoVector>& vertices)
0453
0454 //
0455 // 23. Tube Section Twisted along Its Axis:
0456 //
0457 // G4TwistedTubs(const G4String& pName,
0458 //        G4double  twistedangle,
0459 //        G4double  endinnerrad,
0460 //        G4double  endouterrad,
0461 //        G4double  halfzlen,
0462 //        G4double  dphi)
0463
0464 //
0465 // 24. Cylindrical Cut Section or Cut Tube:
0466 //
0467 // G4CutTubs(const G4String& pName,
0468 //           G4double  pRMin,
0469 //           G4double  pRMax,
0470 //           G4double  pDz,
0471 //           G4double  pSPhi,
0472 //           G4double  pDPhi,
0473 //           G4ThreeVector pLowNorm,
0474 //           G4ThreeVector pHighNorm)
0475 void doCutTubs(const std::string &name,
0476                double rIn,
0477                double rOut,
0478                double zhalf,
0479                double startPhi,
0480                double deltaPhi,
0481                std::array<double, 3> lowNorm,
0482                std::array<double, 3> highNorm) {
0483   G4CutTubs g4(name,
0484                rIn,
0485                rOut,
0486                zhalf,
0487                startPhi,
0488                deltaPhi,
0489                G4ThreeVector(lowNorm[0], lowNorm[1], lowNorm[2]),
0490                G4ThreeVector(highNorm[0], highNorm[1], highNorm[2]));
0491   DDI::CutTubs dd(
0492       zhalf, rIn, rOut, startPhi, deltaPhi, lowNorm[0], lowNorm[1], lowNorm[2], highNorm[0], highNorm[1], highNorm[2]);
0493   DDCutTubs dds = DDSolidFactory::cuttubs(name,
0494                                           zhalf,
0495                                           rIn,
0496                                           rOut,
0497                                           startPhi,
0498                                           deltaPhi,
0499                                           lowNorm[0],
0500                                           lowNorm[1],
0501                                           lowNorm[2],
0502                                           highNorm[0],
0503                                           highNorm[1],
0504                                           highNorm[2]);
0505   dd.stream(std::cout);
0506   std::cout << std::endl;
0507   std::cout << "\tg4 volume = " << g4.GetCubicVolume() / cm3 << " cm3" << std::endl;
0508   std::cout << "\tdd volume = " << dd.volume() / cm3 << " cm3" << std::endl;
0509   std::cout << "\tDD Information: " << dds << " vol= " << dds.volume() << std::endl;
0510 }
0511
0512 //
0513 // 25. Extruded Polygon:
0514 //
0515 // The extrusion of an arbitrary polygon (extruded solid)
0516 // with fixed outline in the defined Z sections can be defined as follows
0517 // (in a general way, or as special construct with two Z sections):
0518 //
0519 //    G4ExtrudedSolid(const G4String& pName,
0520 //                    std::vector<G4TwoVector> polygon,
0521 //                    std::vector<ZSection> zsections)
0522 //
0523 void doExtrudedPgon(const std::string &name,
0524                     const std::vector<double> x,
0525                     const std::vector<double> y,
0526                     const std::vector<double> z,
0527                     const std::vector<double> zx,
0528                     const std::vector<double> zy,
0529                     const std::vector<double> zscale) {
0530   std::vector<G4TwoVector> polygon;
0531   std::vector<G4ExtrudedSolid::ZSection> zsections;
0532   for (unsigned int it = 0; it < x.size(); ++it)
0533     polygon.emplace_back(x[it], y[it]);
0534   for (unsigned int it = 0; it < z.size(); ++it)
0535     zsections.emplace_back(z[it], G4TwoVector(zx[it], zy[it]), zscale[it]);
0536   G4ExtrudedSolid g4(name, polygon, zsections);
0537   DDI::ExtrudedPolygon dd(x, y, z, zx, zy, zscale);
0538   DDExtrudedPolygon dds = DDSolidFactory::extrudedpolygon(name, x, y, z, zx, zy, zscale);
0539
0540   dd.stream(std::cout);
0541   std::cout << std::endl;
0542   std::cout << "\tg4 volume = " << g4.GetCubicVolume() / cm3 << " cm3" << std::endl;
0543   std::cout << "\tdd volume = " << dd.volume() / cm3 << " cm3" << std::endl;
0544   std::cout << "\tDD Information: " << dds << " vol= " << dds.volume() << std::endl;
0545 }
0546
0547 int main(int argc, char *argv[]) {
0548   double xSemiaxis(2. * cm);
0549   double ySemiAxis(2. * cm);
0550   double zHeight(2. * cm);
0551   std::string name("fred1");
0552
0553   //
0554   // 1. Box:
0555   //
0556   std::cout << "\n\nBox tests\n" << std::endl;
0557   doBox(name, xSemiaxis, ySemiAxis, zHeight);
0558   std::cout << std::endl;
0559
0560   //
0561   // 2. Cylindrical Section or Tube:
0562   //
0563   std::cout << "\n\nTub tests\n" << std::endl;
0564   double rIn = 10. * cm;
0565   double rOut = 15. * cm;
0566   double zhalf = 20. * cm;
0567   double startPhi = 0. * deg;
0568   double deltaPhi = 90. * deg;
0569   doTubs(name, rIn, rOut, zhalf, startPhi, deltaPhi);
0570   std::cout << std::endl;
0571
0572   //
0573   // 3. Cone or Conical section:
0574   //
0575   std::cout << "\n\nCons tests\n" << std::endl;
0576   double rIn2 = 20. * cm;
0577   double rOut2 = 25. * cm;
0578   doCons(name, rIn, rOut, rIn2, rOut2, zhalf, startPhi, deltaPhi);
0579   std::cout << std::endl;
0580
0581   //
0582   // 5. Trapezoid:
0583   //
0584   std::cout << "\n\nTrapezoid tests\n" << std::endl;
0585   double dx1 = 10. * cm;
0586   double dx2 = 30. * cm;
0587   double dy1 = 15. * cm;
0588   double dy2 = 30. * cm;
0589   double dz = 60. * cm;
0590   doTrd(name, dx1, dx2, dy1, dy2, dz);
0591   std::cout << std::endl;
0592
0593   //
0594   // 6. Generic Trapezoid:
0595   //
0596   std::cout << "\n\nGeneric Trapezoid tests\n" << std::endl;
0597   double pTheta = 0. * deg;
0598   double pPhi = 0. * deg;
0599   double pDy1 = 30. * cm;
0600   double pDx1 = 30. * cm;
0601   double pDx2 = 30. * cm;
0602   double pAlp1 = 0. * deg;
0603   double pDy2 = 15. * cm;
0604   double pDx3 = 10. * cm;
0605   double pDx4 = 10. * cm;
0606   double pAlp2 = 0. * deg;
0607   doTrap(name, dz, pTheta, pPhi, pDy1, pDx1, pDx2, pAlp1, pDy2, pDx3, pDx4, pAlp2);
0608   std::cout << std::endl;
0609
0610   //
0611   // 7. Sphere or Spherical Shell Section:
0612   //
0613   std::cout << "\n\nSphere tests\n" << std::endl;
0614   std::cout << "This next should be the same as a 2cm ball: " << std::endl;
0615   doSphere("fred1", 0.0 * cm, 2.0 * cm, 0. * deg, 360. * deg, 0., 180. * deg);
0616   std::cout << "Manual computation gives: " << 4. / 3. * Geom::pi() * 2.0 * cm * 2.0 * cm * 2.0 * cm / cm3 << std::endl;
0617   std::cout << "If you mess up phi and theta you get: " << std::endl;
0618   doSphere("fred1", 0.0 * cm, 2.0 * cm, 0. * deg, 180. * deg, 0., 360. * deg);
0619   std::cout << "\n1 cm thick shell: " << std::endl;
0620   doSphere("fred1", 2.0 * cm, 3.0 * cm, 0. * deg, 360. * deg, 0., 180. * deg);
0621   std::cout << "Manual computation gives: "
0622             << 4. / 3. * Geom::pi() * 3.0 * cm * 3.0 * cm * 3.0 * cm / cm3 -
0623                    4. / 3. * Geom::pi() * 2.0 * cm * 2.0 * cm * 2.0 * cm / cm3
0624             << std::endl;
0625   std::cout << "\nHALF of the above 1 cm thick shell: " << std::endl;
0626   doSphere("fred1", 2.0 * cm, 3.0 * cm, 0. * deg, 180. * deg, 0., 180. * deg);
0627   std::cout << "Manual computation gives: "
0628             << (4. / 3. * Geom::pi() * 3.0 * cm * 3.0 * cm * 3.0 * cm / cm3 -
0629                 4. / 3. * Geom::pi() * 2.0 * cm * 2.0 * cm * 2.0 * cm / cm3) /
0630                    2.
0631             << std::endl;
0632   std::cout << "\n30 degree span in theta; full phi \"top\" hemisphere" << std::endl;
0633   doSphere("fred1", 2.0 * cm, 3.0 * cm, 0. * deg, 360. * deg, 10. * deg, 30. * deg);
0634   std::cout << "\n30 degree span in theta; full phi \"bottom\" hemisphere; "
0635                "mirror of above, so should be same."
0636             << std::endl;
0637   doSphere("fred1", 2.0 * cm, 3.0 * cm, 0. * deg, 360. * deg, 140. * deg, 30. * deg);
0638   std::cout << "\n30 degree span in theta; full phi around equator (should be "
0639                "bigger than above)"
0640             << std::endl;
0641   doSphere("fred1", 2.0 * cm, 3.0 * cm, 0. * deg, 360. * deg, 75. * deg, 30. * deg);
0642
0643   //
0644   // 9. Torus:
0645   //
0646   std::cout << "\n\nTorus tests\n" << std::endl;
0647   double radius = 200. * cm;
0648   doTorus(name, rIn, rOut, radius, startPhi, deltaPhi);
0649   std::cout << std::endl;
0650
0651   //
0652   // 10. Polycons:
0653   //
0654   std::cout << "\n\nPolycons tests\n" << std::endl;
0655   double phiStart = 45. * deg;
0656   double phiTotal = 325. * deg;
0657   double inner[] = {0, 0, 0, 0, 0, 0, 0, 0, 0};
0658   std::vector<double> rInner(inner, inner + sizeof(inner) / sizeof(double));
0659   double outer[] = {0, 10, 10, 5, 5, 10, 10, 2, 2};
0660   std::vector<double> rOuter(outer, outer + sizeof(outer) / sizeof(double));
0661   double pl[] = {5, 7, 9, 11, 25, 27, 29, 31, 35};
0662   std::vector<double> z(pl, pl + sizeof(pl) / sizeof(double));
0663   doPolycone1(name, phiStart, phiTotal, z, rInner, rOuter);
0664   std::cout << std::endl;
0665
0666   doPolycone2(name, phiStart, phiTotal, z, rOuter);
0667   std::cout << std::endl;
0668
0669   //
0670   // 11. Polyhedra (PGON):
0671   //
0672   std::cout << "\n\nPolyhedra tests\n" << std::endl;
0673   int sides = 3;
0674   doPolyhedra1(name, sides, phiStart, phiTotal, z, rInner, rOuter);
0675   std::cout << std::endl;
0676
0677   doPolyhedra2(name, sides, phiStart, phiTotal, z, rOuter);
0678   std::cout << std::endl;
0679
0680   //
0681   // 12. Tube with an elliptical cross section:
0682   //
0683
0684   std::cout << "\n\nElliptical Tube tests\n" << std::endl;
0685   doEllipticalTube(name, xSemiaxis, ySemiAxis, zHeight);
0686   std::cout << std::endl;
0687   ySemiAxis = 3. * cm;
0688   doEllipticalTube(name, xSemiaxis, ySemiAxis, zHeight);
0689   std::cout << std::endl;
0690   xSemiaxis = 3. * cm;
0691   ySemiAxis = 2. * cm;
0692   zHeight = 10. * cm;
0693   doEllipticalTube(name, xSemiaxis, ySemiAxis, zHeight);
0694   std::cout << std::endl;
0695   xSemiaxis = 300. * cm;
0696   ySemiAxis = 400. * cm;
0697   zHeight = 3000. * cm;
0698   doEllipticalTube(name, xSemiaxis, ySemiAxis, zHeight);
0699
0700   //
0701   // 14. Cone with Elliptical Cross Section:
0702   //
0703
0704   //
0705   // 15. Paraboloid, a solid with parabolic profile:
0706   //
0707
0708   //
0709   // 16. Tube with Hyperbolic Profile:
0710   //
0711
0712   //
0713   // 17. Tetrahedra:
0714   //
0715
0716   //
0717   // 18. Extruded Polygon:
0718   //
0719
0720   //
0721   // 19. Box Twisted:
0722   //
0723
0724   //
0725   // 20. Trapezoid Twisted along One Axis:
0726   //
0727
0728   //
0729   // 21. Twisted Trapezoid with x and y dimensions varying along z:
0730   //
0731
0732   //
0733   // 22. Generic trapezoid with optionally collapsing vertices:
0734   //
0735
0736   //
0737   // 23. Tube Section Twisted along Its Axis:
0738   //
0739
0740   //
0741   // 24. Cylindrical Cut Section or Cut Tube:
0742   //
0743   std::cout << "\n\nCutTub tests\n" << std::endl;
0744   std::array<double, 3> lowNorm = {{0, -0.7, -0.71}};
0745   std::array<double, 3> highNorm = {{0.7, 0, 0.71}};
0746
0747   doCutTubs(name, rIn, rOut, zhalf, startPhi, deltaPhi, lowNorm, highNorm);
0748   std::cout << std::endl;
0749
0750   //
0751   // 25. Extruded Polygon:
0752   //
0753   // The extrusion of an arbitrary polygon (extruded solid)
0754   // with fixed outline in the defined Z sections can be defined as follows
0755   // (in a general way, or as special construct with two Z sections):
0756   //
0757   //    G4ExtrudedSolid(const G4String& pName,
0758   //                    std::vector<G4TwoVector> polygon,
0759   //                    std::vector<ZSection> zsections)
0760   //
0761   std::cout << "\n\nExtruded Polygon tests\n" << std::endl;
0762   std::vector<double> x = {-300, -300, 300, 300, 150, 150, -150, -150};
0763   std::vector<double> y = {-300, 300, 300, -300, -300, 150, 150, -300};
0764   std::vector<double> epz = {-600, -150, 100, 600};
0765   std::vector<double> zx = {0, 0, 0, 0};
0766   std::vector<double> zy = {300, -300, 0, 300};
0767   std::vector<double> zscale = {8, 10, 6, 12};
0768
0769   doExtrudedPgon(name, x, y, epz, zx, zy, zscale);
0770   std::cout << std::endl;
0771
0772   return EXIT_SUCCESS;
0773 }