u4/tests/G4VSolid_Test.cc

0001
0002 // ./G4VSolid_Test.sh
0003
0004 #include <iostream>
0005 #include <iomanip>
0006 #include <sstream>
0007 #include <vector>
0008 #include <string>
0009 #include <cstring>
0010 #include <cassert>
0011
0012 #include "G4ThreeVector.hh"
0013 #include "G4Polycone.hh"
0014 #include "G4Ellipsoid.hh"
0015 #include "G4UnionSolid.hh"
0016
0017 #include "sgeomdefs.h"
0018 #include "sfmt.h"
0019 #include "ssolid.h"
0020
0021 #include <CLHEP/Units/SystemOfUnits.h>
0022 using CLHEP::mm ;
0023 using CLHEP::deg ;
0024
0025
0026 G4VSolid* MakePolycone( double r, double z0 , double z1 )
0027 {
0028     G4double phiStart = 0.00*deg ;
0029     G4double phiTotal = 360.00*deg ;
0030     G4int numZPlanes = 2 ;
0031     G4double zPlane[] = { z0    , z1 } ;
0032     G4double rInner[] = {  0.0  , 0.0   } ;
0033     G4double rOuter[] = {  r    , r  } ;
0034
0035     G4VSolid* solid_1_2 = new G4Polycone(
0036                                "_1_2",
0037                                phiStart,
0038                                phiTotal,
0039                                numZPlanes,
0040                                zPlane,
0041                                rInner,
0042                                rOuter
0043                                );
0044
0045     return solid_1_2 ;
0046 }
0047
0048 G4VSolid* MakeUpperHemiEllipsoid( double P_I_R, double P_I_H )
0049 {
0050     G4VSolid* solid_I = new G4Ellipsoid(
0051                     "_I",
0052                     P_I_R,
0053                     P_I_R,
0054                     P_I_H,
0055                     0, // pzBottomCut -> equator
0056                     P_I_H // pzTopCut -> top
0057                     );
0058     return solid_I ;
0059 }
0060
0061 G4VSolid* MakeLowerHemiEllipsoid( double P_I_R, double P_I_H )
0062 {
0063     G4VSolid* solid_III = new G4Ellipsoid(
0064                       "_III",
0065                       P_I_R,
0066                       P_I_R,
0067                       P_I_H,
0068                       -P_I_H,
0069                       0);
0070
0071     return solid_III ;
0072 }
0073
0074
0075 G4VSolid* MakeInner1()
0076 {
0077     G4double P_I_R = 254.*mm ;
0078     G4double P_I_H = 190.*mm ;
0079     return MakeUpperHemiEllipsoid( P_I_R, P_I_H );
0080 }
0081
0082
0083 G4VSolid* MakeInner2(const char* style)
0084 {
0085     bool hemi_only = strcmp(style, "NNVT") == 0 ;
0086
0087     G4double P_I_R = 254.*mm ;
0088     G4double P_I_H = 190.*mm ;
0089     G4double m2_h  = 5.*mm ;
0090
0091     G4VSolid* solid_1_2 = MakePolycone( P_I_R, -m2_h, 0. ) ;
0092     G4VSolid* solid_III = MakeLowerHemiEllipsoid( P_I_R, P_I_H );
0093
0094     G4VSolid* solid_1_3 = new G4UnionSolid(
0095                  "_1_3",
0096                  solid_1_2,
0097                  solid_III,
0098                  0,
0099                  G4ThreeVector(0,0,-m2_h )
0100                  );
0101
0102     return hemi_only ? solid_III : solid_1_3 ;
0103 }
0104
0105
0106
0107 void DumpDist(const char* style, const char* desc)
0108 {
0109     G4ThreeVector dir(0,0,1) ;  // +Z
0110     G4VSolid* sol[2] ;
0111     sol[0] = MakeInner1() ;
0112     sol[1] = MakeInner2(style) ;
0113
0114     EInside  inn[2] ;
0115     G4double d2i[2] ;
0116     G4double d2o[2] ;
0117     G4double dis[2] ;
0118
0119     double zstep = 5. ;
0120     const char* spacer = "    " ;
0121
0122     std::cout
0123         << std::setw(55) << ""
0124         << spacer
0125         << std::setw(10) << style
0126         << std::setw(20) << desc
0127         << std::endl
0128         ;
0129
0130     std::cout
0131         << std::setw(15) << "pos"
0132         << std::setw(10) << "INNER1"
0133         << std::setw(10) << "DistToIn"
0134         << std::setw(10) << "DistToOut"
0135         << std::setw(10) << "Dist"
0136         << spacer
0137         << std::setw(10) << "INNER2"
0138         << std::setw(10) << "DistToIn"
0139         << std::setw(10) << "DistToOut"
0140         << std::setw(10) << "Dist"
0141         << spacer
0142         << std::setw(5) << "trig"
0143         << std::endl
0144         ;
0145
0146     for(double z=20 ; z >= -300 ; z -= zstep )
0147     {
0148         zstep = z < 10 && z > -10 ? 1. : 5. ;
0149         G4ThreeVector pos(0,0,z) ;
0150
0151         for(int i=0 ; i < 2 ; i++)
0152         {
0153             d2i[i] = sol[i]->DistanceToIn(pos, dir);
0154             d2o[i] = sol[i]->DistanceToOut(pos, dir);
0155             dis[i] = ssolid::Distance_( sol[i], pos, dir, inn[i] );
0156         }
0157
0158         double dist1 = d2i[0] ;
0159         double dist2 = d2i[1] ;
0160
0161         bool trig = false ;
0162         if(dist1 == kInfinity)
0163         {
0164             trig = false;
0165         }
0166         else if(dist1>dist2)
0167         {
0168             trig = false;
0169         }
0170         else
0171         {
0172             trig = true;
0173         }
0174
0175
0176         std::cout
0177             << std::setw(15) << sfmt::Format(&pos)
0178             << std::setw(10) << sgeomdefs::EInside_(inn[0])
0179             << sfmt::Format(d2i[0])
0180             << sfmt::Format(d2o[0])
0181             << sfmt::Format(dis[0])
0182             << spacer
0183             << std::setw(10) << sgeomdefs::EInside_(inn[1])
0184             << sfmt::Format(d2i[1])
0185             << sfmt::Format(d2o[1])
0186             << sfmt::Format(dis[1])
0187             << spacer
0188             << std::setw(5) << ( trig ? "YES" : "NO " )
0189             << std::endl
0190             ;
0191     }
0192
0193
0194 }
0195
0196
0197
0198 int main()
0199 {
0200     DumpDist("NNVT", " INNER2: Just Lower-Hemi-Ellipsoid");
0201     DumpDist("HAMA", " INNER2: Union of Polycone and Lower-Hemi-Ellipsoid" );
0202     return 0 ;
0203 }
0204
0205
0206 /**
0207
0208 320 // Approximate nearest distance from the point p to the union of
0209 321 // two solids
0210 322
0211 323 G4double
0212 324 G4UnionSolid::DistanceToIn( const G4ThreeVector& p) const
0213 325 {
0214 326 #ifdef G4BOOLDEBUG
0215 327   if( Inside(p) == kInside )
0216 328   {
0217 329     G4cout << "WARNING - Invalid call in "
0218 330            << "G4UnionSolid::DistanceToIn(p)" << G4endl
0219 331            << "  Point p is inside !" << G4endl;
0220 332     G4cout << "          p = " << p << G4endl;
0221 333     G4cerr << "WARNING - Invalid call in "
0222 334            << "G4UnionSolid::DistanceToIn(p)" << G4endl
0223 335            << "  Point p is inside !" << G4endl;
0224 336     G4cerr << "          p = " << p << G4endl;
0225 337   }
0226 338 #endif
0227 339   G4double distA = fPtrSolidA->DistanceToIn(p) ;
0228 340   G4double distB = fPtrSolidB->DistanceToIn(p) ;
0229 341   G4double safety = std::min(distA,distB) ;
0230 342   if(safety < 0.0) safety = 0.0 ;
0231 343   return safety ;
0232 344 }
0233 345
0234
0235
0236 // note that fPtrSolidB will be a G4DisplacedSolid when there is a RHS transform
0237
0238 **/
0239
0240
0241