G4Box.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 //
00027 // $Id: G4Box.cc 69788 2013-05-15 12:06:57Z gcosmo $
00028 //
00029 // 
00030 //
00031 // Implementation for G4Box class
00032 //
00033 //  24.06.98 - V.Grichine: insideEdge in DistanceToIn(p,v)
00034 //  20.09.98 - V.Grichine: new algorithm of DistanceToIn(p,v)
00035 //  07.05.00 - V.Grichine: d= DistanceToIn(p,v), if d<e/2, d=0
00036 //  09.06.00 - V.Grichine: safety in DistanceToIn(p) against Inside(p)=kOutside
00037 //             and information before exception in DistanceToOut(p,v,...)
00038 //  15.11.00 - D.Williams, V.Grichine: bug fixed in CalculateExtent - change
00039 //                                     algorithm for rotated vertices
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 // Constructor - check & set half widths
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) )  // limit to thickness of surfaces
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 // Fake default constructor - sets only member data and allocates memory
00081 //                            for usage restricted to object persistency.
00082 
00083 G4Box::G4Box( __void__& a )
00084   : G4CSGSolid(a), fDx(0.), fDy(0.), fDz(0.)
00085 {
00086 }
00087 
00089 //
00090 // Destructor
00091 
00092 G4Box::~G4Box()
00093 {
00094 }
00095 
00097 //
00098 // Copy constructor
00099 
00100 G4Box::G4Box(const G4Box& rhs)
00101   : G4CSGSolid(rhs), fDx(rhs.fDx), fDy(rhs.fDy), fDz(rhs.fDz)
00102 {
00103 }
00104 
00106 //
00107 // Assignment operator
00108 
00109 G4Box& G4Box::operator = (const G4Box& rhs) 
00110 {
00111    // Check assignment to self
00112    //
00113    if (this == &rhs)  { return *this; }
00114 
00115    // Copy base class data
00116    //
00117    G4CSGSolid::operator=(rhs);
00118 
00119    // Copy data
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)  // limit to thickness of surfaces
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)  // limit to thickness of surfaces
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)  // limit to thickness of surfaces
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 // Dispatch to parameterisation for replication mechanism dimension
00193 // computation & modification.
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 // Calculate extent under transform and specified limit
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     // Special case handling for unrotated boxes
00214     // Compute x/y/z mins and maxs respecting limits, with early returns
00215     // if outside limits. Then switch() on pAxis
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  // General rotated case - create and clip mesh to boundaries
00286   {
00287     G4bool existsAfterClip = false ;
00288     G4ThreeVectorList* vertices ;
00289 
00290     pMin = +kInfinity ;
00291     pMax = -kInfinity ;
00292 
00293     // Calculate rotated vertex coordinates
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         // Add 2*tolerance to avoid precision troubles
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         // Check to see if endpoints are in the solid
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       // Check for case where completely enveloping clipping volume
00349       // If point inside then we are confident that the solid completely
00350       // envelopes the clipping volume. Hence set min/max extents according
00351       // to clipping volume extents along the specified axis.
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 // Return whether point inside/outside/on surface, using tolerance
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 // Calculate side nearest to p, and return normal
00401 // If two sides are equidistant, normal of first side (x/y/z) 
00402 // encountered returned
00403 
00404 G4ThreeVector G4Box::SurfaceNormal( const G4ThreeVector& p) const
00405 {
00406   G4double distx, disty, distz ;
00407   G4ThreeVector norm(0.,0.,0.);
00408 
00409   // Calculate distances as if in 1st octant
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   // New code for particle on surface including edges and corners with specific
00416   // normals
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)         // on X/mX surface and around
00431   {
00432     noSurfaces ++; 
00433     if ( p.x() >= 0. )  { normX= nX ; }       // on +X surface : (1,0,0)
00434     else                { normX= nmX; }       //                 (-1,0,0)
00435     sumnorm= normX; 
00436   }
00437 
00438   if (disty <= delta)    // on one of the +Y or -Y surfaces
00439   {
00440     noSurfaces ++; 
00441     if ( p.y() >= 0. )  { normY= nY;  }       // on +Y surface
00442     else                { normY= nmY; }
00443     sumnorm += normY; 
00444   }
00445 
00446   if (distz <= delta)    // on one of the +Z or -Z surfaces
00447   {
00448     noSurfaces ++; 
00449     if ( p.z() >= 0. )  { normZ= nZ;  }       // on +Z surface
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       // norm = sumnorm . unit(); 
00466       if( noSurfaces == 2 )
00467       { 
00468         // 2 surfaces -> on edge 
00469         norm = invSqrt2 * sumnorm; 
00470       }
00471       else
00472       { 
00473         // 3 surfaces (on corner)
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 // Algorithm for SurfaceNormal() following the original specification
00493 // for points not on the surface
00494 
00495 G4ThreeVector G4Box::ApproxSurfaceNormal( const G4ThreeVector& p ) const
00496 {
00497   G4double distx, disty, distz ;
00498   G4ThreeVector norm(0.,0.,0.);
00499 
00500   // Calculate distances as if in 1st octant
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 )     // Closest to X
00509     {
00510       if ( p.x() < 0 ) { norm = G4ThreeVector(-1.0,0,0) ; }
00511       else             { norm = G4ThreeVector( 1.0,0,0) ; }
00512     }
00513     else                      // Closest to Z
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 )      // Closest to Y
00522     {
00523       if ( p.y() < 0 ) { norm = G4ThreeVector(0,-1.0,0) ; }
00524       else             { norm = G4ThreeVector(0, 1.0,0) ; }
00525     }
00526     else                       // Closest to Z
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 // Calculate distance to box from an outside point
00538 // - return kInfinity if no intersection.
00539 //
00540 // ALGORITHM:
00541 //
00542 // Check that if point lies outside x/y/z extent of box, travel is towards
00543 // the box (ie. there is a possibility of an intersection)
00544 //
00545 // Calculate pairs of minimum and maximum distances for x/y/z travel for
00546 // intersection with the box's x/y/z extent.
00547 // If there is a valid intersection, it is given by the maximum min distance
00548 // (ie. distance to satisfy x/y/z intersections) *if* <= minimum max distance
00549 // (ie. distance after which 1+ of x/y/z intersections not satisfied)
00550 //
00551 // NOTE:
00552 //
00553 // `Inside' safe - meaningful answers given if point is inside the exact
00554 // shape.
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 ; // , sminx ;
00561   G4double smax=kInfinity, smaxy, smaxz ; // , smaxx ;  // they always > 0
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 ;     // minimum distance to x surface of shape
00568   safy = std::fabs(p.y()) - fDy ;
00569   safz = std::fabs(p.z()) - fDz ;
00570 
00571   // Will we intersect?
00572   // If safx/y/z is >-tol/2 the point is outside/on the box's x/y/z extent.
00573   // If both p.x/y/z and v.x/y/z repectively are both positive/negative,
00574   // travel is in a direction away from the shape.
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 ;  // travel away or parallel within tolerance
00581   }
00582 
00583   // Compute min / max distances for x/y/z travel:
00584   // X Planes
00585 
00586   if ( v.x() )  // != 0
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   // Y Planes
00603 
00604   if ( v.y() )  // != 0
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 ;  // touch XY corner
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   // Z planes
00630 
00631   if ( v.z() )  // != 0
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 ;    // touch ZX or ZY corners
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)) // travel over edge
00657   {
00658     return kInfinity ;
00659   }
00660   if (smin < delta)  { smin = 0.0 ; }
00661 
00662   return smin ;
00663 }
00664 
00666 // 
00667 // Appoximate distance to box.
00668 // Returns largest perpendicular distance to the closest x/y/z sides of
00669 // the box, which is the most fast estimation of the shortest distance to box
00670 // - If inside return 0
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 // Calculate distance to surface of box from inside
00690 // by calculating distances to box's x/y/z planes.
00691 // Smallest distance is exact distance to exiting.
00692 // - Eliminate one side of each pair by considering direction of v
00693 // - when leaving a surface & v.close, return 0
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 ; }  // All normals are valid
00705 
00706   if (v.x() > 0)   // X planes
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)   // Y planes
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)        // Z planes
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 // Calculate exact shortest distance to any boundary from inside
00870 // - If outside return 0
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   // shortest Dist to any boundary now MIN(safx1,safx2,safy1..)
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 // Create a List containing the transformed vertices
00915 // Ordering [0-3] -fDz cross section
00916 //          [4-7] +fDz cross section such that [0] is below [4],
00917 //                                             [1] below [5] etc.
00918 // Note:
00919 //  Caller has deletion resposibility
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 // GetEntityType
00960 
00961 G4GeometryType G4Box::GetEntityType() const
00962 {
00963   return G4String("G4Box");
00964 }
00965 
00967 //
00968 // Stream object contents to an output stream
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 // GetPointOnSurface
00990 //
00991 // Return a point (G4ThreeVector) randomly and uniformly selected
00992 // on the solid surface
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 // Make a clone of the object
01032 //
01033 G4VSolid* G4Box::Clone() const
01034 {
01035   return new G4Box(*this);
01036 }
01037 
01039 //
01040 // Methods for visualisation
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 }

Generated on Mon May 27 17:47:46 2013 for Geant4 by  doxygen 1.4.7