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