00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042 #include "G4Box.hh"
00043
00044 #include "G4SystemOfUnits.hh"
00045 #include "G4VoxelLimits.hh"
00046 #include "G4AffineTransform.hh"
00047 #include "Randomize.hh"
00048
00049 #include "G4VPVParameterisation.hh"
00050
00051 #include "G4VGraphicsScene.hh"
00052 #include "G4Polyhedron.hh"
00053 #include "G4NURBS.hh"
00054 #include "G4NURBSbox.hh"
00055 #include "G4VisExtent.hh"
00056
00058
00059
00060
00061 G4Box::G4Box(const G4String& pName,
00062 G4double pX,
00063 G4double pY,
00064 G4double pZ)
00065 : G4CSGSolid(pName), fDx(pX), fDy(pY), fDz(pZ)
00066 {
00067 if ( (pX < 2*kCarTolerance)
00068 && (pY < 2*kCarTolerance)
00069 && (pZ < 2*kCarTolerance) )
00070 {
00071 std::ostringstream message;
00072 message << "Dimensions too small for Solid: " << GetName() << "!" << G4endl
00073 << " hX, hY, hZ = " << pX << ", " << pY << ", " << pZ;
00074 G4Exception("G4Box::G4Box()", "GeomSolids0002", FatalException, message);
00075 }
00076 }
00077
00079
00080
00081
00082
00083 G4Box::G4Box( __void__& a )
00084 : G4CSGSolid(a), fDx(0.), fDy(0.), fDz(0.)
00085 {
00086 }
00087
00089
00090
00091
00092 G4Box::~G4Box()
00093 {
00094 }
00095
00097
00098
00099
00100 G4Box::G4Box(const G4Box& rhs)
00101 : G4CSGSolid(rhs), fDx(rhs.fDx), fDy(rhs.fDy), fDz(rhs.fDz)
00102 {
00103 }
00104
00106
00107
00108
00109 G4Box& G4Box::operator = (const G4Box& rhs)
00110 {
00111
00112
00113 if (this == &rhs) { return *this; }
00114
00115
00116
00117 G4CSGSolid::operator=(rhs);
00118
00119
00120
00121 fDx = rhs.fDx;
00122 fDy = rhs.fDy;
00123 fDz = rhs.fDz;
00124
00125 return *this;
00126 }
00127
00129
00130 void G4Box::SetXHalfLength(G4double dx)
00131 {
00132 if(dx > 2*kCarTolerance)
00133 {
00134 fDx = dx;
00135 }
00136 else
00137 {
00138 std::ostringstream message;
00139 message << "Dimension X too small for solid: " << GetName() << "!"
00140 << G4endl
00141 << " hX = " << dx;
00142 G4Exception("G4Box::SetXHalfLength()", "GeomSolids0002",
00143 FatalException, message);
00144 }
00145 fCubicVolume= 0.;
00146 fSurfaceArea= 0.;
00147 fpPolyhedron = 0;
00148 }
00149
00150 void G4Box::SetYHalfLength(G4double dy)
00151 {
00152 if(dy > 2*kCarTolerance)
00153 {
00154 fDy = dy;
00155 }
00156 else
00157 {
00158 std::ostringstream message;
00159 message << "Dimension Y too small for solid: " << GetName() << "!"
00160 << G4endl
00161 << " hY = " << dy;
00162 G4Exception("G4Box::SetYHalfLength()", "GeomSolids0002",
00163 FatalException, message);
00164 }
00165 fCubicVolume= 0.;
00166 fSurfaceArea= 0.;
00167 fpPolyhedron = 0;
00168 }
00169
00170 void G4Box::SetZHalfLength(G4double dz)
00171 {
00172 if(dz > 2*kCarTolerance)
00173 {
00174 fDz = dz;
00175 }
00176 else
00177 {
00178 std::ostringstream message;
00179 message << "Dimension Z too small for solid: " << GetName() << "!"
00180 << G4endl
00181 << " hZ = " << dz;
00182 G4Exception("G4Box::SetZHalfLength()", "GeomSolids0002",
00183 FatalException, message);
00184 }
00185 fCubicVolume= 0.;
00186 fSurfaceArea= 0.;
00187 fpPolyhedron = 0;
00188 }
00189
00191
00192
00193
00194
00195 void G4Box::ComputeDimensions(G4VPVParameterisation* p,
00196 const G4int n,
00197 const G4VPhysicalVolume* pRep)
00198 {
00199 p->ComputeDimensions(*this,n,pRep);
00200 }
00201
00203
00204
00205
00206 G4bool G4Box::CalculateExtent(const EAxis pAxis,
00207 const G4VoxelLimits& pVoxelLimit,
00208 const G4AffineTransform& pTransform,
00209 G4double& pMin, G4double& pMax) const
00210 {
00211 if (!pTransform.IsRotated())
00212 {
00213
00214
00215
00216
00217 G4double xoffset,xMin,xMax;
00218 G4double yoffset,yMin,yMax;
00219 G4double zoffset,zMin,zMax;
00220
00221 xoffset = pTransform.NetTranslation().x() ;
00222 xMin = xoffset - fDx ;
00223 xMax = xoffset + fDx ;
00224
00225 if (pVoxelLimit.IsXLimited())
00226 {
00227 if ((xMin > pVoxelLimit.GetMaxXExtent()+kCarTolerance) ||
00228 (xMax < pVoxelLimit.GetMinXExtent()-kCarTolerance)) { return false ; }
00229 else
00230 {
00231 xMin = std::max(xMin, pVoxelLimit.GetMinXExtent());
00232 xMax = std::min(xMax, pVoxelLimit.GetMaxXExtent());
00233 }
00234 }
00235 yoffset = pTransform.NetTranslation().y() ;
00236 yMin = yoffset - fDy ;
00237 yMax = yoffset + fDy ;
00238
00239 if (pVoxelLimit.IsYLimited())
00240 {
00241 if ((yMin > pVoxelLimit.GetMaxYExtent()+kCarTolerance) ||
00242 (yMax < pVoxelLimit.GetMinYExtent()-kCarTolerance)) { return false ; }
00243 else
00244 {
00245 yMin = std::max(yMin, pVoxelLimit.GetMinYExtent());
00246 yMax = std::min(yMax, pVoxelLimit.GetMaxYExtent());
00247 }
00248 }
00249 zoffset = pTransform.NetTranslation().z() ;
00250 zMin = zoffset - fDz ;
00251 zMax = zoffset + fDz ;
00252
00253 if (pVoxelLimit.IsZLimited())
00254 {
00255 if ((zMin > pVoxelLimit.GetMaxZExtent()+kCarTolerance) ||
00256 (zMax < pVoxelLimit.GetMinZExtent()-kCarTolerance)) { return false ; }
00257 else
00258 {
00259 zMin = std::max(zMin, pVoxelLimit.GetMinZExtent());
00260 zMax = std::min(zMax, pVoxelLimit.GetMaxZExtent());
00261 }
00262 }
00263 switch (pAxis)
00264 {
00265 case kXAxis:
00266 pMin = xMin ;
00267 pMax = xMax ;
00268 break ;
00269 case kYAxis:
00270 pMin=yMin;
00271 pMax=yMax;
00272 break;
00273 case kZAxis:
00274 pMin=zMin;
00275 pMax=zMax;
00276 break;
00277 default:
00278 break;
00279 }
00280 pMin -= kCarTolerance ;
00281 pMax += kCarTolerance ;
00282
00283 return true;
00284 }
00285 else
00286 {
00287 G4bool existsAfterClip = false ;
00288 G4ThreeVectorList* vertices ;
00289
00290 pMin = +kInfinity ;
00291 pMax = -kInfinity ;
00292
00293
00294
00295 vertices = CreateRotatedVertices(pTransform) ;
00296 ClipCrossSection(vertices,0,pVoxelLimit,pAxis,pMin,pMax) ;
00297 ClipCrossSection(vertices,4,pVoxelLimit,pAxis,pMin,pMax) ;
00298 ClipBetweenSections(vertices,0,pVoxelLimit,pAxis,pMin,pMax) ;
00299
00300 if (pVoxelLimit.IsLimited(pAxis) == false)
00301 {
00302 if ( (pMin != kInfinity) || (pMax != -kInfinity) )
00303 {
00304 existsAfterClip = true ;
00305
00306
00307
00308 pMin -= kCarTolerance;
00309 pMax += kCarTolerance;
00310 }
00311 }
00312 else
00313 {
00314 G4ThreeVector clipCentre(
00315 ( pVoxelLimit.GetMinXExtent()+pVoxelLimit.GetMaxXExtent())*0.5,
00316 ( pVoxelLimit.GetMinYExtent()+pVoxelLimit.GetMaxYExtent())*0.5,
00317 ( pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5);
00318
00319 if ( (pMin != kInfinity) || (pMax != -kInfinity) )
00320 {
00321 existsAfterClip = true ;
00322
00323
00324
00325
00326 clipCentre(pAxis) = pVoxelLimit.GetMinExtent(pAxis);
00327
00328 if (Inside(pTransform.Inverse().TransformPoint(clipCentre)) != kOutside)
00329 {
00330 pMin = pVoxelLimit.GetMinExtent(pAxis);
00331 }
00332 else
00333 {
00334 pMin -= kCarTolerance;
00335 }
00336 clipCentre(pAxis) = pVoxelLimit.GetMaxExtent(pAxis);
00337
00338 if (Inside(pTransform.Inverse().TransformPoint(clipCentre)) != kOutside)
00339 {
00340 pMax = pVoxelLimit.GetMaxExtent(pAxis);
00341 }
00342 else
00343 {
00344 pMax += kCarTolerance;
00345 }
00346 }
00347
00348
00349
00350
00351
00352
00353 else if (Inside(pTransform.Inverse().TransformPoint(clipCentre))
00354 != kOutside)
00355 {
00356 existsAfterClip = true ;
00357 pMin = pVoxelLimit.GetMinExtent(pAxis) ;
00358 pMax = pVoxelLimit.GetMaxExtent(pAxis) ;
00359 }
00360 }
00361 delete vertices;
00362 return existsAfterClip;
00363 }
00364 }
00365
00367
00368
00369
00370 EInside G4Box::Inside(const G4ThreeVector& p) const
00371 {
00372 static const G4double delta=0.5*kCarTolerance;
00373 EInside in = kOutside ;
00374 G4ThreeVector q(std::fabs(p.x()), std::fabs(p.y()), std::fabs(p.z()));
00375
00376 if ( q.x() <= (fDx - delta) )
00377 {
00378 if (q.y() <= (fDy - delta) )
00379 {
00380 if ( q.z() <= (fDz - delta) ) { in = kInside ; }
00381 else if ( q.z() <= (fDz + delta) ) { in = kSurface ; }
00382 }
00383 else if ( q.y() <= (fDy + delta) )
00384 {
00385 if ( q.z() <= (fDz + delta) ) { in = kSurface ; }
00386 }
00387 }
00388 else if ( q.x() <= (fDx + delta) )
00389 {
00390 if ( q.y() <= (fDy + delta) )
00391 {
00392 if ( q.z() <= (fDz + delta) ) { in = kSurface ; }
00393 }
00394 }
00395 return in ;
00396 }
00397
00399
00400
00401
00402
00403
00404 G4ThreeVector G4Box::SurfaceNormal( const G4ThreeVector& p) const
00405 {
00406 G4double distx, disty, distz ;
00407 G4ThreeVector norm(0.,0.,0.);
00408
00409
00410
00411 distx = std::fabs(std::fabs(p.x()) - fDx) ;
00412 disty = std::fabs(std::fabs(p.y()) - fDy) ;
00413 distz = std::fabs(std::fabs(p.z()) - fDz) ;
00414
00415
00416
00417
00418 static const G4double delta = 0.5*kCarTolerance;
00419 const G4ThreeVector nX = G4ThreeVector( 1.0, 0,0 );
00420 const G4ThreeVector nmX = G4ThreeVector(-1.0, 0,0 );
00421 const G4ThreeVector nY = G4ThreeVector( 0, 1.0,0 );
00422 const G4ThreeVector nmY = G4ThreeVector( 0,-1.0,0 );
00423 const G4ThreeVector nZ = G4ThreeVector( 0, 0, 1.0);
00424 const G4ThreeVector nmZ = G4ThreeVector( 0, 0,- 1.0);
00425
00426 G4ThreeVector normX(0.,0.,0.), normY(0.,0.,0.), normZ(0.,0.,0.);
00427 G4ThreeVector sumnorm(0., 0., 0.);
00428 G4int noSurfaces=0;
00429
00430 if (distx <= delta)
00431 {
00432 noSurfaces ++;
00433 if ( p.x() >= 0. ) { normX= nX ; }
00434 else { normX= nmX; }
00435 sumnorm= normX;
00436 }
00437
00438 if (disty <= delta)
00439 {
00440 noSurfaces ++;
00441 if ( p.y() >= 0. ) { normY= nY; }
00442 else { normY= nmY; }
00443 sumnorm += normY;
00444 }
00445
00446 if (distz <= delta)
00447 {
00448 noSurfaces ++;
00449 if ( p.z() >= 0. ) { normZ= nZ; }
00450 else { normZ= nmZ; }
00451 sumnorm += normZ;
00452 }
00453
00454 static const G4double invSqrt2 = 1.0 / std::sqrt(2.0);
00455 static const G4double invSqrt3 = 1.0 / std::sqrt(3.0);
00456
00457 if( noSurfaces > 0 )
00458 {
00459 if( noSurfaces == 1 )
00460 {
00461 norm= sumnorm;
00462 }
00463 else
00464 {
00465
00466 if( noSurfaces == 2 )
00467 {
00468
00469 norm = invSqrt2 * sumnorm;
00470 }
00471 else
00472 {
00473
00474 norm = invSqrt3 * sumnorm;
00475 }
00476 }
00477 }
00478 else
00479 {
00480 #ifdef G4CSGDEBUG
00481 G4Exception("G4Box::SurfaceNormal(p)", "Notification", JustWarning,
00482 "Point p is not on surface !?" );
00483 #endif
00484 norm = ApproxSurfaceNormal(p);
00485 }
00486
00487 return norm;
00488 }
00489
00491
00492
00493
00494
00495 G4ThreeVector G4Box::ApproxSurfaceNormal( const G4ThreeVector& p ) const
00496 {
00497 G4double distx, disty, distz ;
00498 G4ThreeVector norm(0.,0.,0.);
00499
00500
00501
00502 distx = std::fabs(std::fabs(p.x()) - fDx) ;
00503 disty = std::fabs(std::fabs(p.y()) - fDy) ;
00504 distz = std::fabs(std::fabs(p.z()) - fDz) ;
00505
00506 if ( distx <= disty )
00507 {
00508 if ( distx <= distz )
00509 {
00510 if ( p.x() < 0 ) { norm = G4ThreeVector(-1.0,0,0) ; }
00511 else { norm = G4ThreeVector( 1.0,0,0) ; }
00512 }
00513 else
00514 {
00515 if ( p.z() < 0 ) { norm = G4ThreeVector(0,0,-1.0) ; }
00516 else { norm = G4ThreeVector(0,0, 1.0) ; }
00517 }
00518 }
00519 else
00520 {
00521 if ( disty <= distz )
00522 {
00523 if ( p.y() < 0 ) { norm = G4ThreeVector(0,-1.0,0) ; }
00524 else { norm = G4ThreeVector(0, 1.0,0) ; }
00525 }
00526 else
00527 {
00528 if ( p.z() < 0 ) { norm = G4ThreeVector(0,0,-1.0) ; }
00529 else { norm = G4ThreeVector(0,0, 1.0) ; }
00530 }
00531 }
00532 return norm;
00533 }
00534
00536
00537
00538
00539
00540
00541
00542
00543
00544
00545
00546
00547
00548
00549
00550
00551
00552
00553
00554
00555
00556 G4double G4Box::DistanceToIn(const G4ThreeVector& p,
00557 const G4ThreeVector& v) const
00558 {
00559 G4double safx, safy, safz ;
00560 G4double smin=0.0, sminy, sminz ;
00561 G4double smax=kInfinity, smaxy, smaxz ;
00562 G4double stmp ;
00563 G4double sOut=kInfinity, sOuty=kInfinity, sOutz=kInfinity ;
00564
00565 static const G4double delta = 0.5*kCarTolerance;
00566
00567 safx = std::fabs(p.x()) - fDx ;
00568 safy = std::fabs(p.y()) - fDy ;
00569 safz = std::fabs(p.z()) - fDz ;
00570
00571
00572
00573
00574
00575
00576 if ( ((p.x()*v.x() >= 0.0) && (safx > -delta))
00577 || ((p.y()*v.y() >= 0.0) && (safy > -delta))
00578 || ((p.z()*v.z() >= 0.0) && (safz > -delta)) )
00579 {
00580 return kInfinity ;
00581 }
00582
00583
00584
00585
00586 if ( v.x() )
00587 {
00588 stmp = 1.0/std::fabs(v.x()) ;
00589
00590 if (safx >= 0.0)
00591 {
00592 smin = safx*stmp ;
00593 smax = (fDx+std::fabs(p.x()))*stmp ;
00594 }
00595 else
00596 {
00597 if (v.x() < 0) { sOut = (fDx + p.x())*stmp ; }
00598 else { sOut = (fDx - p.x())*stmp ; }
00599 }
00600 }
00601
00602
00603
00604 if ( v.y() )
00605 {
00606 stmp = 1.0/std::fabs(v.y()) ;
00607
00608 if (safy >= 0.0)
00609 {
00610 sminy = safy*stmp ;
00611 smaxy = (fDy+std::fabs(p.y()))*stmp ;
00612
00613 if (sminy > smin) { smin=sminy ; }
00614 if (smaxy < smax) { smax=smaxy ; }
00615
00616 if (smin >= (smax-delta))
00617 {
00618 return kInfinity ;
00619 }
00620 }
00621 else
00622 {
00623 if (v.y() < 0) { sOuty = (fDy + p.y())*stmp ; }
00624 else { sOuty = (fDy - p.y())*stmp ; }
00625 if( sOuty < sOut ) { sOut = sOuty ; }
00626 }
00627 }
00628
00629
00630
00631 if ( v.z() )
00632 {
00633 stmp = 1.0/std::fabs(v.z()) ;
00634
00635 if ( safz >= 0.0 )
00636 {
00637 sminz = safz*stmp ;
00638 smaxz = (fDz+std::fabs(p.z()))*stmp ;
00639
00640 if (sminz > smin) { smin = sminz ; }
00641 if (smaxz < smax) { smax = smaxz ; }
00642
00643 if (smin >= (smax-delta))
00644 {
00645 return kInfinity ;
00646 }
00647 }
00648 else
00649 {
00650 if (v.z() < 0) { sOutz = (fDz + p.z())*stmp ; }
00651 else { sOutz = (fDz - p.z())*stmp ; }
00652 if( sOutz < sOut ) { sOut = sOutz ; }
00653 }
00654 }
00655
00656 if (sOut <= (smin + delta))
00657 {
00658 return kInfinity ;
00659 }
00660 if (smin < delta) { smin = 0.0 ; }
00661
00662 return smin ;
00663 }
00664
00666
00667
00668
00669
00670
00671
00672 G4double G4Box::DistanceToIn(const G4ThreeVector& p) const
00673 {
00674 G4double safex, safey, safez, safe = 0.0 ;
00675
00676 safex = std::fabs(p.x()) - fDx ;
00677 safey = std::fabs(p.y()) - fDy ;
00678 safez = std::fabs(p.z()) - fDz ;
00679
00680 if (safex > safe) { safe = safex ; }
00681 if (safey > safe) { safe = safey ; }
00682 if (safez > safe) { safe = safez ; }
00683
00684 return safe ;
00685 }
00686
00688
00689
00690
00691
00692
00693
00694
00695 G4double G4Box::DistanceToOut( const G4ThreeVector& p,const G4ThreeVector& v,
00696 const G4bool calcNorm,
00697 G4bool *validNorm,G4ThreeVector *n) const
00698 {
00699 ESide side = kUndefined ;
00700 G4double pdist,stmp,snxt=kInfinity;
00701
00702 static const G4double delta = 0.5*kCarTolerance;
00703
00704 if (calcNorm) { *validNorm = true ; }
00705
00706 if (v.x() > 0)
00707 {
00708 pdist = fDx - p.x() ;
00709
00710 if (pdist > delta)
00711 {
00712 snxt = pdist/v.x() ;
00713 side = kPX ;
00714 }
00715 else
00716 {
00717 if (calcNorm) { *n = G4ThreeVector(1,0,0) ; }
00718 return snxt = 0 ;
00719 }
00720 }
00721 else if (v.x() < 0)
00722 {
00723 pdist = fDx + p.x() ;
00724
00725 if (pdist > delta)
00726 {
00727 snxt = -pdist/v.x() ;
00728 side = kMX ;
00729 }
00730 else
00731 {
00732 if (calcNorm) { *n = G4ThreeVector(-1,0,0) ; }
00733 return snxt = 0 ;
00734 }
00735 }
00736
00737 if (v.y() > 0)
00738 {
00739 pdist = fDy-p.y();
00740
00741 if (pdist > delta)
00742 {
00743 stmp = pdist/v.y();
00744
00745 if (stmp < snxt)
00746 {
00747 snxt = stmp;
00748 side = kPY;
00749 }
00750 }
00751 else
00752 {
00753 if (calcNorm) { *n = G4ThreeVector(0,1,0) ; }
00754 return snxt = 0 ;
00755 }
00756 }
00757 else if (v.y() < 0)
00758 {
00759 pdist = fDy + p.y() ;
00760
00761 if (pdist > delta)
00762 {
00763 stmp = -pdist/v.y();
00764
00765 if ( stmp < snxt )
00766 {
00767 snxt = stmp;
00768 side = kMY;
00769 }
00770 }
00771 else
00772 {
00773 if (calcNorm) { *n = G4ThreeVector(0,-1,0) ; }
00774 return snxt = 0 ;
00775 }
00776 }
00777
00778 if (v.z() > 0)
00779 {
00780 pdist = fDz-p.z();
00781
00782 if ( pdist > delta )
00783 {
00784 stmp = pdist/v.z();
00785
00786 if ( stmp < snxt )
00787 {
00788 snxt = stmp;
00789 side = kPZ;
00790 }
00791 }
00792 else
00793 {
00794 if (calcNorm) { *n = G4ThreeVector(0,0,1) ; }
00795 return snxt = 0 ;
00796 }
00797 }
00798 else if (v.z() < 0)
00799 {
00800 pdist = fDz + p.z();
00801
00802 if ( pdist > delta )
00803 {
00804 stmp = -pdist/v.z();
00805
00806 if ( stmp < snxt )
00807 {
00808 snxt = stmp;
00809 side = kMZ;
00810 }
00811 }
00812 else
00813 {
00814 if (calcNorm) { *n = G4ThreeVector(0,0,-1) ; }
00815 return snxt = 0 ;
00816 }
00817 }
00818
00819 if (calcNorm)
00820 {
00821 switch (side)
00822 {
00823 case kPX:
00824 *n=G4ThreeVector(1,0,0);
00825 break;
00826 case kMX:
00827 *n=G4ThreeVector(-1,0,0);
00828 break;
00829 case kPY:
00830 *n=G4ThreeVector(0,1,0);
00831 break;
00832 case kMY:
00833 *n=G4ThreeVector(0,-1,0);
00834 break;
00835 case kPZ:
00836 *n=G4ThreeVector(0,0,1);
00837 break;
00838 case kMZ:
00839 *n=G4ThreeVector(0,0,-1);
00840 break;
00841 default:
00842 G4cout << G4endl;
00843 DumpInfo();
00844 std::ostringstream message;
00845 G4int oldprc = message.precision(16);
00846 message << "Undefined side for valid surface normal to solid."
00847 << G4endl
00848 << "Position:" << G4endl << G4endl
00849 << "p.x() = " << p.x()/mm << " mm" << G4endl
00850 << "p.y() = " << p.y()/mm << " mm" << G4endl
00851 << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl
00852 << "Direction:" << G4endl << G4endl
00853 << "v.x() = " << v.x() << G4endl
00854 << "v.y() = " << v.y() << G4endl
00855 << "v.z() = " << v.z() << G4endl << G4endl
00856 << "Proposed distance :" << G4endl << G4endl
00857 << "snxt = " << snxt/mm << " mm" << G4endl;
00858 message.precision(oldprc);
00859 G4Exception("G4Box::DistanceToOut(p,v,..)", "GeomSolids1002",
00860 JustWarning, message);
00861 break;
00862 }
00863 }
00864 return snxt;
00865 }
00866
00868
00869
00870
00871
00872 G4double G4Box::DistanceToOut(const G4ThreeVector& p) const
00873 {
00874 G4double safx1,safx2,safy1,safy2,safz1,safz2,safe=0.0;
00875
00876 #ifdef G4CSGDEBUG
00877 if( Inside(p) == kOutside )
00878 {
00879 G4int oldprc = G4cout.precision(16) ;
00880 G4cout << G4endl ;
00881 DumpInfo();
00882 G4cout << "Position:" << G4endl << G4endl ;
00883 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ;
00884 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ;
00885 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ;
00886 G4cout.precision(oldprc) ;
00887 G4Exception("G4Box::DistanceToOut(p)", "GeomSolids1002",
00888 JustWarning, "Point p is outside !?" );
00889 }
00890 #endif
00891
00892 safx1 = fDx - p.x() ;
00893 safx2 = fDx + p.x() ;
00894 safy1 = fDy - p.y() ;
00895 safy2 = fDy + p.y() ;
00896 safz1 = fDz - p.z() ;
00897 safz2 = fDz + p.z() ;
00898
00899
00900
00901 if (safx2 < safx1) { safe = safx2; }
00902 else { safe = safx1; }
00903 if (safy1 < safe) { safe = safy1; }
00904 if (safy2 < safe) { safe = safy2; }
00905 if (safz1 < safe) { safe = safz1; }
00906 if (safz2 < safe) { safe = safz2; }
00907
00908 if (safe < 0) { safe = 0 ; }
00909 return safe ;
00910 }
00911
00913
00914
00915
00916
00917
00918
00919
00920
00921 G4ThreeVectorList*
00922 G4Box::CreateRotatedVertices(const G4AffineTransform& pTransform) const
00923 {
00924 G4ThreeVectorList* vertices = new G4ThreeVectorList();
00925
00926 if (vertices)
00927 {
00928 vertices->reserve(8);
00929 G4ThreeVector vertex0(-fDx,-fDy,-fDz) ;
00930 G4ThreeVector vertex1(fDx,-fDy,-fDz) ;
00931 G4ThreeVector vertex2(fDx,fDy,-fDz) ;
00932 G4ThreeVector vertex3(-fDx,fDy,-fDz) ;
00933 G4ThreeVector vertex4(-fDx,-fDy,fDz) ;
00934 G4ThreeVector vertex5(fDx,-fDy,fDz) ;
00935 G4ThreeVector vertex6(fDx,fDy,fDz) ;
00936 G4ThreeVector vertex7(-fDx,fDy,fDz) ;
00937
00938 vertices->push_back(pTransform.TransformPoint(vertex0));
00939 vertices->push_back(pTransform.TransformPoint(vertex1));
00940 vertices->push_back(pTransform.TransformPoint(vertex2));
00941 vertices->push_back(pTransform.TransformPoint(vertex3));
00942 vertices->push_back(pTransform.TransformPoint(vertex4));
00943 vertices->push_back(pTransform.TransformPoint(vertex5));
00944 vertices->push_back(pTransform.TransformPoint(vertex6));
00945 vertices->push_back(pTransform.TransformPoint(vertex7));
00946 }
00947 else
00948 {
00949 DumpInfo();
00950 G4Exception("G4Box::CreateRotatedVertices()",
00951 "GeomSolids0003", FatalException,
00952 "Error in allocation of vertices. Out of memory !");
00953 }
00954 return vertices;
00955 }
00956
00958
00959
00960
00961 G4GeometryType G4Box::GetEntityType() const
00962 {
00963 return G4String("G4Box");
00964 }
00965
00967
00968
00969
00970 std::ostream& G4Box::StreamInfo(std::ostream& os) const
00971 {
00972 G4int oldprc = os.precision(16);
00973 os << "-----------------------------------------------------------\n"
00974 << " *** Dump for solid - " << GetName() << " ***\n"
00975 << " ===================================================\n"
00976 << " Solid type: G4Box\n"
00977 << " Parameters: \n"
00978 << " half length X: " << fDx/mm << " mm \n"
00979 << " half length Y: " << fDy/mm << " mm \n"
00980 << " half length Z: " << fDz/mm << " mm \n"
00981 << "-----------------------------------------------------------\n";
00982 os.precision(oldprc);
00983
00984 return os;
00985 }
00986
00988
00989
00990
00991
00992
00993
00994 G4ThreeVector G4Box::GetPointOnSurface() const
00995 {
00996 G4double px, py, pz, select, sumS;
00997 G4double Sxy = fDx*fDy, Sxz = fDx*fDz, Syz = fDy*fDz;
00998
00999 sumS = Sxy + Sxz + Syz;
01000 select = sumS*G4UniformRand();
01001
01002 if( select < Sxy )
01003 {
01004 px = -fDx +2*fDx*G4UniformRand();
01005 py = -fDy +2*fDy*G4UniformRand();
01006
01007 if(G4UniformRand() > 0.5) { pz = fDz; }
01008 else { pz = -fDz; }
01009 }
01010 else if ( ( select - Sxy ) < Sxz )
01011 {
01012 px = -fDx +2*fDx*G4UniformRand();
01013 pz = -fDz +2*fDz*G4UniformRand();
01014
01015 if(G4UniformRand() > 0.5) { py = fDy; }
01016 else { py = -fDy; }
01017 }
01018 else
01019 {
01020 py = -fDy +2*fDy*G4UniformRand();
01021 pz = -fDz +2*fDz*G4UniformRand();
01022
01023 if(G4UniformRand() > 0.5) { px = fDx; }
01024 else { px = -fDx; }
01025 }
01026 return G4ThreeVector(px,py,pz);
01027 }
01028
01030
01031
01032
01033 G4VSolid* G4Box::Clone() const
01034 {
01035 return new G4Box(*this);
01036 }
01037
01039
01040
01041
01042 void G4Box::DescribeYourselfTo (G4VGraphicsScene& scene) const
01043 {
01044 scene.AddSolid (*this);
01045 }
01046
01047 G4VisExtent G4Box::GetExtent() const
01048 {
01049 return G4VisExtent (-fDx, fDx, -fDy, fDy, -fDz, fDz);
01050 }
01051
01052 G4Polyhedron* G4Box::CreatePolyhedron () const
01053 {
01054 return new G4PolyhedronBox (fDx, fDy, fDz);
01055 }
01056
01057 G4NURBS* G4Box::CreateNURBS () const
01058 {
01059 return new G4NURBSbox (fDx, fDy, fDz);
01060 }