G4PolyconeSide Class Reference

#include <G4PolyconeSide.hh>

Inheritance diagram for G4PolyconeSide:

G4VCSGface

Public Member Functions

 G4PolyconeSide (const G4PolyconeSideRZ *prevRZ, const G4PolyconeSideRZ *tail, const G4PolyconeSideRZ *head, const G4PolyconeSideRZ *nextRZ, G4double phiStart, G4double deltaPhi, G4bool phiIsOpen, G4bool isAllBehind=false)
virtual ~G4PolyconeSide ()
 G4PolyconeSide (const G4PolyconeSide &source)
G4PolyconeSideoperator= (const G4PolyconeSide &source)
G4bool Intersect (const G4ThreeVector &p, const G4ThreeVector &v, G4bool outgoing, G4double surfTolerance, G4double &distance, G4double &distFromSurface, G4ThreeVector &normal, G4bool &isAllBehind)
G4double Distance (const G4ThreeVector &p, G4bool outgoing)
EInside Inside (const G4ThreeVector &p, G4double tolerance, G4double *bestDistance)
G4ThreeVector Normal (const G4ThreeVector &p, G4double *bestDistance)
G4double Extent (const G4ThreeVector axis)
void CalculateExtent (const EAxis axis, const G4VoxelLimits &voxelLimit, const G4AffineTransform &tranform, G4SolidExtentList &extentList)
G4VCSGfaceClone ()
G4double SurfaceArea ()
G4ThreeVector GetPointOnFace ()
 G4PolyconeSide (__void__ &)

Protected Member Functions

G4double DistanceAway (const G4ThreeVector &p, G4bool opposite, G4double &distOutside2, G4double *rzNorm=0)
G4bool PointOnCone (const G4ThreeVector &hit, G4double normSign, const G4ThreeVector &p, const G4ThreeVector &v, G4ThreeVector &normal)
void CopyStuff (const G4PolyconeSide &source)
G4double GetPhi (const G4ThreeVector &p)

Static Protected Member Functions

static void FindLineIntersect (G4double x1, G4double y1, G4double tx1, G4double ty1, G4double x2, G4double y2, G4double tx2, G4double ty2, G4double &x, G4double &y)

Protected Attributes

G4double r [2]
G4double z [2]
G4double startPhi
G4double deltaPhi
G4bool phiIsOpen
G4bool allBehind
G4IntersectingConecone
G4double rNorm
G4double zNorm
G4double rS
G4double zS
G4double length
G4double prevRS
G4double prevZS
G4double nextRS
G4double nextZS
G4double rNormEdge [2]
G4double zNormEdge [2]
G4int ncorners
G4ThreeVectorcorners

Detailed Description

Definition at line 67 of file G4PolyconeSide.hh.


Constructor & Destructor Documentation

G4PolyconeSide::G4PolyconeSide ( const G4PolyconeSideRZ prevRZ,
const G4PolyconeSideRZ tail,
const G4PolyconeSideRZ head,
const G4PolyconeSideRZ nextRZ,
G4double  phiStart,
G4double  deltaPhi,
G4bool  phiIsOpen,
G4bool  isAllBehind = false 
)

Definition at line 57 of file G4PolyconeSide.cc.

References allBehind, cone, corners, deltaPhi, G4GeometryTolerance::GetInstance(), G4GeometryTolerance::GetSurfaceTolerance(), length, ncorners, nextRS, nextZS, phiIsOpen, prevRS, prevZS, G4PolyconeSideRZ::r, r, rNorm, rNormEdge, rS, startPhi, G4PolyconeSideRZ::z, z, zNorm, zNormEdge, and zS.

Referenced by Clone().

00065   : ncorners(0), corners(0)
00066 {
00067   kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
00068   fSurfaceArea = 0.0;
00069   fPhi.first = G4ThreeVector(0,0,0);
00070   fPhi.second= 0.0;
00071 
00072   //
00073   // Record values
00074   //
00075   r[0] = tail->r; z[0] = tail->z;
00076   r[1] = head->r; z[1] = head->z;
00077   
00078   phiIsOpen = thePhiIsOpen;
00079   if (phiIsOpen)
00080   {
00081     deltaPhi = theDeltaPhi;
00082     startPhi = thePhiStart;
00083 
00084     //
00085     // Set phi values to our conventions
00086     //
00087     while (deltaPhi < 0.0) deltaPhi += twopi;
00088     while (startPhi < 0.0) startPhi += twopi;
00089     
00090     //
00091     // Calculate corner coordinates
00092     //
00093     ncorners = 4;
00094     corners = new G4ThreeVector[ncorners];
00095     
00096     corners[0] = G4ThreeVector( tail->r*std::cos(startPhi),
00097                                 tail->r*std::sin(startPhi), tail->z );
00098     corners[1] = G4ThreeVector( head->r*std::cos(startPhi),
00099                                 head->r*std::sin(startPhi), head->z );
00100     corners[2] = G4ThreeVector( tail->r*std::cos(startPhi+deltaPhi),
00101                                 tail->r*std::sin(startPhi+deltaPhi), tail->z );
00102     corners[3] = G4ThreeVector( head->r*std::cos(startPhi+deltaPhi),
00103                                 head->r*std::sin(startPhi+deltaPhi), head->z );
00104   }
00105   else
00106   {
00107     deltaPhi = twopi;
00108     startPhi = 0.0;
00109   }
00110   
00111   allBehind = isAllBehind;
00112     
00113   //
00114   // Make our intersecting cone
00115   //
00116   cone = new G4IntersectingCone( r, z );
00117   
00118   //
00119   // Calculate vectors in r,z space
00120   //
00121   rS = r[1]-r[0]; zS = z[1]-z[0];
00122   length = std::sqrt( rS*rS + zS*zS);
00123   rS /= length; zS /= length;
00124   
00125   rNorm = +zS;
00126   zNorm = -rS;
00127   
00128   G4double lAdj;
00129   
00130   prevRS = r[0]-prevRZ->r;
00131   prevZS = z[0]-prevRZ->z;
00132   lAdj = std::sqrt( prevRS*prevRS + prevZS*prevZS );
00133   prevRS /= lAdj;
00134   prevZS /= lAdj;
00135 
00136   rNormEdge[0] = rNorm + prevZS;
00137   zNormEdge[0] = zNorm - prevRS;
00138   lAdj = std::sqrt( rNormEdge[0]*rNormEdge[0] + zNormEdge[0]*zNormEdge[0] );
00139   rNormEdge[0] /= lAdj;
00140   zNormEdge[0] /= lAdj;
00141 
00142   nextRS = nextRZ->r-r[1];
00143   nextZS = nextRZ->z-z[1];
00144   lAdj = std::sqrt( nextRS*nextRS + nextZS*nextZS );
00145   nextRS /= lAdj;
00146   nextZS /= lAdj;
00147 
00148   rNormEdge[1] = rNorm + nextZS;
00149   zNormEdge[1] = zNorm - nextRS;
00150   lAdj = std::sqrt( rNormEdge[1]*rNormEdge[1] + zNormEdge[1]*zNormEdge[1] );
00151   rNormEdge[1] /= lAdj;
00152   zNormEdge[1] /= lAdj;
00153 }

G4PolyconeSide::~G4PolyconeSide (  )  [virtual]

Definition at line 176 of file G4PolyconeSide.cc.

References cone, corners, and phiIsOpen.

00177 {
00178   delete cone;
00179   if (phiIsOpen)  { delete [] corners; }
00180 }

G4PolyconeSide::G4PolyconeSide ( const G4PolyconeSide source  ) 

Definition at line 186 of file G4PolyconeSide.cc.

References CopyStuff().

00187   : G4VCSGface(), ncorners(0), corners(0)
00188 {
00189   CopyStuff( source );
00190 }

G4PolyconeSide::G4PolyconeSide ( __void__ &   ) 

Definition at line 160 of file G4PolyconeSide.cc.

References r, rNormEdge, z, and zNormEdge.

00161   : startPhi(0.), deltaPhi(0.), phiIsOpen(false), allBehind(false),
00162     cone(0), rNorm(0.), zNorm(0.), rS(0.), zS(0.), length(0.),
00163     prevRS(0.), prevZS(0.), nextRS(0.), nextZS(0.), ncorners(0), corners(0),
00164     kCarTolerance(0.), fSurfaceArea(0.)
00165 {
00166   r[0] = r[1] = 0.;
00167   z[0] = z[1] = 0.;
00168   rNormEdge[0]= rNormEdge[1] = 0.;
00169   zNormEdge[0]= zNormEdge[1] = 0.;
00170 }


Member Function Documentation

void G4PolyconeSide::CalculateExtent ( const EAxis  axis,
const G4VoxelLimits voxelLimit,
const G4AffineTransform tranform,
G4SolidExtentList extentList 
) [virtual]

Implements G4VCSGface.

Definition at line 552 of file G4PolyconeSide.cc.

References G4SolidExtentList::AddSurface(), G4ClippablePolygon::AddVertexInOrder(), G4AffineTransform::ApplyPointTransform(), G4ClippablePolygon::ClearAllVertices(), DBL_MIN, deltaPhi, FindLineIntersect(), kMaxMeshSections, kMeshAngleDefault, kMinMeshSections, nextRS, nextZS, G4ClippablePolygon::PartialClip(), phiIsOpen, prevRS, prevZS, r, rNorm, rS, G4ClippablePolygon::SetNormal(), startPhi, G4AffineTransform::TransformAxis(), z, and zS.

00556 {
00557   G4ClippablePolygon polygon;
00558   
00559   //
00560   // Here we will approximate (ala G4Cons) and divide our conical section
00561   // into segments, like G4Polyhedra. When doing so, the radius
00562   // is extented far enough such that the segments always lie
00563   // just outside the surface of the conical section we are
00564   // approximating.
00565   //
00566   
00567   //
00568   // Choose phi size of our segment(s) based on constants as
00569   // defined in meshdefs.hh
00570   //
00571   G4int numPhi = (G4int)(deltaPhi/kMeshAngleDefault) + 1;
00572   if (numPhi < kMinMeshSections) 
00573     numPhi = kMinMeshSections;
00574   else if (numPhi > kMaxMeshSections)
00575     numPhi = kMaxMeshSections;
00576     
00577   G4double sigPhi = deltaPhi/numPhi;
00578   
00579   //
00580   // Determine radius factor to keep segments outside
00581   //
00582   G4double rFudge = 1.0/std::cos(0.5*sigPhi);
00583   
00584   //
00585   // Decide which radius to use on each end of the side,
00586   // and whether a transition mesh is required
00587   //
00588   // {r0,z0}  - Beginning of this side
00589   // {r1,z1}  - Ending of this side
00590   // {r2,z0}  - Beginning of transition piece connecting previous
00591   //            side (and ends at beginning of this side)
00592   //
00593   // So, order is 2 --> 0 --> 1.
00594   //                    -------
00595   //
00596   // r2 < 0 indicates that no transition piece is required
00597   //
00598   G4double r0, r1, r2, z0, z1;
00599   
00600   r2 = -1;  // By default: no transition piece
00601   
00602   if (rNorm < -DBL_MIN)
00603   {
00604     //
00605     // This side faces *inward*, and so our mesh has
00606     // the same radius
00607     //
00608     r1 = r[1];
00609     z1 = z[1];
00610     z0 = z[0];
00611     r0 = r[0];
00612     
00613     r2 = -1;
00614     
00615     if (prevZS > DBL_MIN)
00616     {
00617       //
00618       // The previous side is facing outwards
00619       //
00620       if ( prevRS*zS - prevZS*rS > 0 )
00621       {
00622         //
00623         // Transition was convex: build transition piece
00624         //
00625         if (r[0] > DBL_MIN) r2 = r[0]*rFudge;
00626       }
00627       else
00628       {
00629         //
00630         // Transition was concave: short this side
00631         //
00632         FindLineIntersect( z0, r0, zS, rS,
00633                            z0, r0*rFudge, prevZS, prevRS*rFudge, z0, r0 );
00634       }
00635     }
00636     
00637     if ( nextZS > DBL_MIN && (rS*nextZS - zS*nextRS < 0) )
00638     {
00639       //
00640       // The next side is facing outwards, forming a 
00641       // concave transition: short this side
00642       //
00643       FindLineIntersect( z1, r1, zS, rS,
00644                          z1, r1*rFudge, nextZS, nextRS*rFudge, z1, r1 );
00645     }
00646   }
00647   else if (rNorm > DBL_MIN)
00648   {
00649     //
00650     // This side faces *outward* and is given a boost to
00651     // it radius
00652     //
00653     r0 = r[0]*rFudge;
00654     z0 = z[0];
00655     r1 = r[1]*rFudge;
00656     z1 = z[1];
00657     
00658     if (prevZS < -DBL_MIN)
00659     {
00660       //
00661       // The previous side is facing inwards
00662       //
00663       if ( prevRS*zS - prevZS*rS > 0 )
00664       {
00665         //
00666         // Transition was convex: build transition piece
00667         //
00668         if (r[0] > DBL_MIN) r2 = r[0];
00669       }
00670       else
00671       {
00672         //
00673         // Transition was concave: short this side
00674         //
00675         FindLineIntersect( z0, r0, zS, rS*rFudge,
00676                            z0, r[0], prevZS, prevRS, z0, r0 );
00677       }
00678     }
00679     
00680     if ( nextZS < -DBL_MIN && (rS*nextZS - zS*nextRS < 0) )
00681     {
00682       //
00683       // The next side is facing inwards, forming a 
00684       // concave transition: short this side
00685       //
00686       FindLineIntersect( z1, r1, zS, rS*rFudge,
00687                          z1, r[1], nextZS, nextRS, z1, r1 );
00688     }
00689   }
00690   else
00691   {
00692     //
00693     // This side is perpendicular to the z axis (is a disk)
00694     //
00695     // Whether or not r0 needs a rFudge factor depends
00696     // on the normal of the previous edge. Similar with r1
00697     // and the next edge. No transition piece is required.
00698     //
00699     r0 = r[0];
00700     r1 = r[1];
00701     z0 = z[0];
00702     z1 = z[1];
00703     
00704     if (prevZS > DBL_MIN) r0 *= rFudge;
00705     if (nextZS > DBL_MIN) r1 *= rFudge;
00706   }
00707   
00708   //
00709   // Loop
00710   //
00711   G4double phi = startPhi, 
00712            cosPhi = std::cos(phi), 
00713            sinPhi = std::sin(phi);
00714   
00715   G4ThreeVector v0( r0*cosPhi, r0*sinPhi, z0 ),
00716                     v1( r1*cosPhi, r1*sinPhi, z1 ),
00717   v2, w0, w1, w2;
00718   transform.ApplyPointTransform( v0 );
00719   transform.ApplyPointTransform( v1 );
00720   
00721   if (r2 >= 0)
00722   {
00723     v2 = G4ThreeVector( r2*cosPhi, r2*sinPhi, z0 );
00724     transform.ApplyPointTransform( v2 );
00725   }
00726 
00727   do
00728   {
00729     phi += sigPhi;
00730     if (numPhi == 1) phi = startPhi+deltaPhi;  // Try to avoid roundoff
00731     cosPhi = std::cos(phi), 
00732     sinPhi = std::sin(phi);
00733     
00734     w0 = G4ThreeVector( r0*cosPhi, r0*sinPhi, z0 );
00735     w1 = G4ThreeVector( r1*cosPhi, r1*sinPhi, z1 );
00736     transform.ApplyPointTransform( w0 );
00737     transform.ApplyPointTransform( w1 );
00738     
00739     G4ThreeVector deltaV = r0 > r1 ? w0-v0 : w1-v1;
00740     
00741     //
00742     // Build polygon, taking special care to keep the vertices
00743     // in order
00744     //
00745     polygon.ClearAllVertices();
00746 
00747     polygon.AddVertexInOrder( v0 );
00748     polygon.AddVertexInOrder( v1 );
00749     polygon.AddVertexInOrder( w1 );
00750     polygon.AddVertexInOrder( w0 );
00751 
00752     //
00753     // Get extent
00754     //
00755     if (polygon.PartialClip( voxelLimit, axis ))
00756     {
00757       //
00758       // Get dot product of normal with target axis
00759       //
00760       polygon.SetNormal( deltaV.cross(v1-v0).unit() );
00761       
00762       extentList.AddSurface( polygon );
00763     }
00764     
00765     if (r2 >= 0)
00766     {
00767       //
00768       // Repeat, for transition piece
00769       //
00770       w2 = G4ThreeVector( r2*cosPhi, r2*sinPhi, z0 );
00771       transform.ApplyPointTransform( w2 );
00772 
00773       polygon.ClearAllVertices();
00774 
00775       polygon.AddVertexInOrder( v2 );
00776       polygon.AddVertexInOrder( v0 );
00777       polygon.AddVertexInOrder( w0 );
00778       polygon.AddVertexInOrder( w2 );
00779 
00780       if (polygon.PartialClip( voxelLimit, axis ))
00781       {
00782         polygon.SetNormal( deltaV.cross(v0-v2).unit() );
00783         
00784         extentList.AddSurface( polygon );
00785       }
00786       
00787       v2 = w2;
00788     }
00789     
00790     //
00791     // Next vertex
00792     //    
00793     v0 = w0;
00794     v1 = w1;
00795   } while( --numPhi > 0 );
00796   
00797   //
00798   // We are almost done. But, it is important that we leave no
00799   // gaps in the surface of our solid. By using rFudge, however,
00800   // we've done exactly that, if we have a phi segment. 
00801   // Add two additional faces if necessary
00802   //
00803   if (phiIsOpen && rNorm > DBL_MIN)
00804   {    
00805     cosPhi = std::cos(startPhi);
00806     sinPhi = std::sin(startPhi);
00807 
00808     G4ThreeVector a0( r[0]*cosPhi, r[0]*sinPhi, z[0] ),
00809                   a1( r[1]*cosPhi, r[1]*sinPhi, z[1] ),
00810                   b0( r0*cosPhi, r0*sinPhi, z[0] ),
00811                   b1( r1*cosPhi, r1*sinPhi, z[1] );
00812   
00813     transform.ApplyPointTransform( a0 );
00814     transform.ApplyPointTransform( a1 );
00815     transform.ApplyPointTransform( b0 );
00816     transform.ApplyPointTransform( b1 );
00817 
00818     polygon.ClearAllVertices();
00819 
00820     polygon.AddVertexInOrder( a0 );
00821     polygon.AddVertexInOrder( a1 );
00822     polygon.AddVertexInOrder( b0 );
00823     polygon.AddVertexInOrder( b1 );
00824     
00825     if (polygon.PartialClip( voxelLimit , axis))
00826     {
00827       G4ThreeVector normal( sinPhi, -cosPhi, 0 );
00828       polygon.SetNormal( transform.TransformAxis( normal ) );
00829         
00830       extentList.AddSurface( polygon );
00831     }
00832     
00833     cosPhi = std::cos(startPhi+deltaPhi);
00834     sinPhi = std::sin(startPhi+deltaPhi);
00835     
00836     a0 = G4ThreeVector( r[0]*cosPhi, r[0]*sinPhi, z[0] ),
00837     a1 = G4ThreeVector( r[1]*cosPhi, r[1]*sinPhi, z[1] ),
00838     b0 = G4ThreeVector( r0*cosPhi, r0*sinPhi, z[0] ),
00839     b1 = G4ThreeVector( r1*cosPhi, r1*sinPhi, z[1] );
00840     transform.ApplyPointTransform( a0 );
00841     transform.ApplyPointTransform( a1 );
00842     transform.ApplyPointTransform( b0 );
00843     transform.ApplyPointTransform( b1 );
00844 
00845     polygon.ClearAllVertices();
00846 
00847     polygon.AddVertexInOrder( a0 );
00848     polygon.AddVertexInOrder( a1 );
00849     polygon.AddVertexInOrder( b0 );
00850     polygon.AddVertexInOrder( b1 );
00851     
00852     if (polygon.PartialClip( voxelLimit, axis ))
00853     {
00854       G4ThreeVector normal( -sinPhi, cosPhi, 0 );
00855       polygon.SetNormal( transform.TransformAxis( normal ) );
00856         
00857       extentList.AddSurface( polygon );
00858     }
00859   }
00860     
00861   return;
00862 }

G4VCSGface* G4PolyconeSide::Clone (  )  [inline, virtual]

Implements G4VCSGface.

Definition at line 101 of file G4PolyconeSide.hh.

References G4PolyconeSide().

00101 { return new G4PolyconeSide( *this ); }

void G4PolyconeSide::CopyStuff ( const G4PolyconeSide source  )  [protected]

Definition at line 212 of file G4PolyconeSide.cc.

References allBehind, cone, corners, deltaPhi, fSurfaceArea, kCarTolerance, length, ncorners, nextRS, nextZS, phiIsOpen, prevRS, prevZS, r, rNorm, rNormEdge, rS, startPhi, z, zNorm, zNormEdge, and zS.

Referenced by G4PolyconeSide(), and operator=().

00213 {
00214   r[0]    = source.r[0];
00215   r[1]    = source.r[1];
00216   z[0]    = source.z[0];
00217   z[1]    = source.z[1];
00218   
00219   startPhi  = source.startPhi;
00220   deltaPhi  = source.deltaPhi;
00221   phiIsOpen  = source.phiIsOpen;
00222   allBehind  = source.allBehind;
00223 
00224   kCarTolerance = source.kCarTolerance;
00225   fSurfaceArea = source.fSurfaceArea;
00226   
00227   cone    = new G4IntersectingCone( *source.cone );
00228   
00229   rNorm    = source.rNorm;
00230   zNorm    = source.zNorm;
00231   rS    = source.rS;
00232   zS    = source.zS;
00233   length    = source.length;
00234   prevRS    = source.prevRS;
00235   prevZS    = source.prevZS;
00236   nextRS    = source.nextRS;
00237   nextZS    = source.nextZS;
00238   
00239   rNormEdge[0]   = source.rNormEdge[0];
00240   rNormEdge[1]  = source.rNormEdge[1];
00241   zNormEdge[0]  = source.zNormEdge[0];
00242   zNormEdge[1]  = source.zNormEdge[1];
00243   
00244   if (phiIsOpen)
00245   {
00246     ncorners = 4;
00247     corners = new G4ThreeVector[ncorners];
00248     
00249     corners[0] = source.corners[0];
00250     corners[1] = source.corners[1];
00251     corners[2] = source.corners[2];
00252     corners[3] = source.corners[3];
00253   }
00254 }

G4double G4PolyconeSide::Distance ( const G4ThreeVector p,
G4bool  outgoing 
) [virtual]

Implements G4VCSGface.

Definition at line 389 of file G4PolyconeSide.cc.

References DistanceAway().

00390 {
00391   G4double normSign = outgoing ? -1 : +1;
00392   G4double distFrom, distOut2;
00393   
00394   //
00395   // We have two tries for each hemisphere. Try the closest first.
00396   //
00397   distFrom = normSign*DistanceAway( p, false, distOut2 );
00398   if (distFrom > -0.5*kCarTolerance )
00399   {
00400     //
00401     // Good answer
00402     //
00403     if (distOut2 > 0) 
00404       return std::sqrt( distFrom*distFrom + distOut2 );
00405     else 
00406       return std::fabs(distFrom);
00407   }
00408   
00409   //
00410   // Try second side. 
00411   //
00412   distFrom = normSign*DistanceAway( p,  true, distOut2 );
00413   if (distFrom > -0.5*kCarTolerance)
00414   {
00415 
00416     if (distOut2 > 0) 
00417       return std::sqrt( distFrom*distFrom + distOut2 );
00418     else
00419       return std::fabs(distFrom);
00420   }
00421   
00422   return kInfinity;
00423 }

G4double G4PolyconeSide::DistanceAway ( const G4ThreeVector p,
G4bool  opposite,
G4double distOutside2,
G4double rzNorm = 0 
) [protected]

Definition at line 906 of file G4PolyconeSide.cc.

References deltaPhi, GetPhi(), length, phiIsOpen, r, rNorm, rNormEdge, rS, sqr(), startPhi, z, zNorm, zNormEdge, and zS.

Referenced by Distance(), Inside(), Intersect(), and Normal().

00910 {
00911   //
00912   // Convert our point to r and z
00913   //
00914   G4double rx = p.perp(), zx = p.z();
00915   
00916   //
00917   // Change sign of r if opposite says we should
00918   //
00919   if (opposite) rx = -rx;
00920   
00921   //
00922   // Calculate return value
00923   //
00924   G4double deltaR  = rx - r[0], deltaZ = zx - z[0];
00925   G4double answer = deltaR*rNorm + deltaZ*zNorm;
00926   
00927   //
00928   // Are we off the surface in r,z space?
00929   //
00930   G4double q = deltaR*rS + deltaZ*zS;
00931   if (q < 0)
00932   {
00933     distOutside2 = q*q;
00934     if (edgeRZnorm) *edgeRZnorm = deltaR*rNormEdge[0] + deltaZ*zNormEdge[0];
00935   }
00936   else if (q > length)
00937   {
00938     distOutside2 = sqr( q-length );
00939     if (edgeRZnorm)
00940     {
00941       deltaR = rx - r[1];
00942       deltaZ = zx - z[1];
00943       *edgeRZnorm = deltaR*rNormEdge[1] + deltaZ*zNormEdge[1];
00944     }
00945   }
00946   else
00947   {
00948     distOutside2 = 0;
00949     if (edgeRZnorm) *edgeRZnorm = answer;
00950   }
00951 
00952   if (phiIsOpen)
00953   {
00954     //
00955     // Finally, check phi
00956     //
00957     G4double phi = GetPhi(p);
00958     while( phi < startPhi ) phi += twopi;
00959     
00960     if (phi > startPhi+deltaPhi)
00961     {
00962       //
00963       // Oops. Are we closer to the start phi or end phi?
00964       //
00965       G4double d1 = phi-startPhi-deltaPhi;
00966       while( phi > startPhi ) phi -= twopi;
00967       G4double d2 = startPhi-phi;
00968       
00969       if (d2 < d1) d1 = d2;
00970       
00971       //
00972       // Add result to our distance
00973       //
00974       G4double dist = d1*rx;
00975       
00976       distOutside2 += dist*dist;
00977       if (edgeRZnorm)
00978       {
00979         *edgeRZnorm = std::max(std::fabs(*edgeRZnorm),std::fabs(dist));
00980       }
00981     }
00982   }
00983 
00984   return answer;
00985 }

G4double G4PolyconeSide::Extent ( const G4ThreeVector  axis  )  [virtual]

Implements G4VCSGface.

Definition at line 488 of file G4PolyconeSide.cc.

References cone, DBL_MIN, deltaPhi, GetPhi(), phiIsOpen, r, startPhi, z, G4IntersectingCone::ZHi(), and G4IntersectingCone::ZLo().

00489 {
00490   if (axis.perp2() < DBL_MIN)
00491   {
00492     //
00493     // Special case
00494     //
00495     return axis.z() < 0 ? -cone->ZLo() : cone->ZHi();
00496   }
00497 
00498   //
00499   // Is the axis pointing inside our phi gap?
00500   //
00501   if (phiIsOpen)
00502   {
00503     G4double phi = GetPhi(axis);
00504     while( phi < startPhi ) phi += twopi;
00505     
00506     if (phi > deltaPhi+startPhi)
00507     {
00508       //
00509       // Yeah, looks so. Make four three vectors defining the phi
00510       // opening
00511       //
00512       G4double cosP = std::cos(startPhi), sinP = std::sin(startPhi);
00513       G4ThreeVector a( r[0]*cosP, r[0]*sinP, z[0] );
00514       G4ThreeVector b( r[1]*cosP, r[1]*sinP, z[1] );
00515       cosP = std::cos(startPhi+deltaPhi); sinP = std::sin(startPhi+deltaPhi);
00516       G4ThreeVector c( r[0]*cosP, r[0]*sinP, z[0] );
00517       G4ThreeVector d( r[1]*cosP, r[1]*sinP, z[1] );
00518       
00519       G4double ad = axis.dot(a),
00520          bd = axis.dot(b),
00521          cd = axis.dot(c),
00522          dd = axis.dot(d);
00523       
00524       if (bd > ad) ad = bd;
00525       if (cd > ad) ad = cd;
00526       if (dd > ad) ad = dd;
00527       
00528       return ad;
00529     }
00530   }
00531 
00532   //
00533   // Check either end
00534   //
00535   G4double aPerp = axis.perp();
00536   
00537   G4double a = aPerp*r[0] + axis.z()*z[0];
00538   G4double b = aPerp*r[1] + axis.z()*z[1];
00539   
00540   if (b > a) a = b;
00541   
00542   return a;
00543 }

void G4PolyconeSide::FindLineIntersect ( G4double  x1,
G4double  y1,
G4double  tx1,
G4double  ty1,
G4double  x2,
G4double  y2,
G4double  tx2,
G4double  ty2,
G4double x,
G4double y 
) [static, protected]

Definition at line 1062 of file G4PolyconeSide.cc.

Referenced by CalculateExtent().

01067 {
01068   //
01069   // The solution is a simple linear equation
01070   //
01071   G4double deter = tx1*ty2 - tx2*ty1;
01072   
01073   G4double s1 = ((x2-x1)*ty2 - tx2*(y2-y1))/deter;
01074   G4double s2 = ((x2-x1)*ty1 - tx1*(y2-y1))/deter;
01075 
01076   //
01077   // We want the answer to not depend on which order the
01078   // lines were specified. Take average.
01079   //
01080   x = 0.5*( x1+s1*tx1 + x2+s2*tx2 );
01081   y = 0.5*( y1+s1*ty1 + y2+s2*ty2 );
01082 }

G4double G4PolyconeSide::GetPhi ( const G4ThreeVector p  )  [protected]

Definition at line 871 of file G4PolyconeSide.cc.

Referenced by DistanceAway(), Extent(), and PointOnCone().

00872 {
00873   G4double val=0.;
00874 
00875   if (fPhi.first != p)
00876   {
00877     val = p.phi();
00878     fPhi.first = p;
00879     fPhi.second = val;
00880   }
00881   else
00882   {
00883     val = fPhi.second;
00884   }
00885   return val;
00886 }

G4ThreeVector G4PolyconeSide::GetPointOnFace (  )  [virtual]

Implements G4VCSGface.

Definition at line 1100 of file G4PolyconeSide.cc.

References deltaPhi, G4UniformRand, r, startPhi, and z.

01101 {
01102   G4double x,y,zz;
01103   G4double rr,phi,dz,dr;
01104   dr=r[1]-r[0];dz=z[1]-z[0];
01105   phi=startPhi+deltaPhi*G4UniformRand();
01106   rr=r[0]+dr*G4UniformRand();
01107  
01108   x=rr*std::cos(phi);
01109   y=rr*std::sin(phi);
01110 
01111   // PolyconeSide has a Ring Form
01112   //
01113   if (dz==0.)
01114   {
01115     zz=z[0];
01116   }
01117   else
01118   {
01119     if(dr==0.)  // PolyconeSide has a Tube Form
01120     {
01121       zz = z[0]+dz*G4UniformRand();
01122     }
01123     else
01124     {
01125       zz = z[0]+(rr-r[0])*dz/dr;
01126     }
01127   }
01128 
01129   return G4ThreeVector(x,y,zz);
01130 }

EInside G4PolyconeSide::Inside ( const G4ThreeVector p,
G4double  tolerance,
G4double bestDistance 
) [virtual]

Implements G4VCSGface.

Definition at line 429 of file G4PolyconeSide.cc.

References DistanceAway(), kInside, kOutside, and kSurface.

00432 {
00433   //
00434   // Check both sides
00435   //
00436   G4double distFrom[2], distOut2[2], dist2[2];
00437   G4double edgeRZnorm[2];
00438      
00439   distFrom[0] =  DistanceAway( p, false, distOut2[0], edgeRZnorm   );
00440   distFrom[1] =  DistanceAway( p,  true, distOut2[1], edgeRZnorm+1 );
00441   
00442   dist2[0] = distFrom[0]*distFrom[0] + distOut2[0];
00443   dist2[1] = distFrom[1]*distFrom[1] + distOut2[1];
00444   
00445   //
00446   // Who's closest?
00447   //
00448   G4int i = std::fabs(dist2[0]) < std::fabs(dist2[1]) ? 0 : 1;
00449   
00450   *bestDistance = std::sqrt( dist2[i] );
00451   
00452   //
00453   // Okay then, inside or out?
00454   //
00455   if ( (std::fabs(edgeRZnorm[i]) < tolerance)
00456     && (distOut2[i] < tolerance*tolerance) )
00457     return kSurface;
00458   else if (edgeRZnorm[i] < 0)
00459     return kInside;
00460   else
00461     return kOutside;
00462 }

G4bool G4PolyconeSide::Intersect ( const G4ThreeVector p,
const G4ThreeVector v,
G4bool  outgoing,
G4double  surfTolerance,
G4double distance,
G4double distFromSurface,
G4ThreeVector normal,
G4bool isAllBehind 
) [virtual]

Implements G4VCSGface.

Definition at line 260 of file G4PolyconeSide.cc.

References allBehind, cone, DBL_MIN, DistanceAway(), G4IntersectingCone::LineHitsCone(), PointOnCone(), rNorm, and zNorm.

00268 {
00269   G4double s1, s2;
00270   G4double normSign = outgoing ? +1 : -1;
00271   
00272   isAllBehind = allBehind;
00273 
00274   //
00275   // Check for two possible intersections
00276   //
00277   G4int nside = cone->LineHitsCone( p, v, &s1, &s2 );
00278   if (nside == 0) return false;
00279     
00280   //
00281   // Check the first side first, since it is (supposed to be) closest
00282   //
00283   G4ThreeVector hit = p + s1*v;
00284   
00285   if (PointOnCone( hit, normSign, p, v, normal ))
00286   {
00287     //
00288     // Good intersection! What about the normal? 
00289     //
00290     if (normSign*v.dot(normal) > 0)
00291     {
00292       //
00293       // We have a valid intersection, but it could very easily
00294       // be behind the point. To decide if we tolerate this,
00295       // we have to see if the point p is on the surface near
00296       // the intersecting point.
00297       //
00298       // What does it mean exactly for the point p to be "near"
00299       // the intersection? It means that if we draw a line from
00300       // p to the hit, the line remains entirely within the
00301       // tolerance bounds of the cone. To test this, we can
00302       // ask if the normal is correct near p.
00303       //
00304       G4double pr = p.perp();
00305       if (pr < DBL_MIN) pr = DBL_MIN;
00306       G4ThreeVector pNormal( rNorm*p.x()/pr, rNorm*p.y()/pr, zNorm );
00307       if (normSign*v.dot(pNormal) > 0)
00308       {
00309         //
00310         // p and intersection in same hemisphere
00311         //
00312         G4double distOutside2;
00313         distFromSurface = -normSign*DistanceAway( p, false, distOutside2 );
00314         if (distOutside2 < surfTolerance*surfTolerance)
00315         {
00316           if (distFromSurface > -surfTolerance)
00317           {
00318             //
00319             // We are just inside or away from the
00320             // surface. Accept *any* value of distance.
00321             //
00322             distance = s1;
00323             return true;
00324           }
00325         }
00326       }
00327       else 
00328         distFromSurface = s1;
00329       
00330       //
00331       // Accept positive distances
00332       //
00333       if (s1 > 0)
00334       {
00335         distance = s1;
00336         return true;
00337       }
00338     }
00339   }  
00340   
00341   if (nside==1) return false;
00342   
00343   //
00344   // Well, try the second hit
00345   //  
00346   hit = p + s2*v;
00347   
00348   if (PointOnCone( hit, normSign, p, v, normal ))
00349   {
00350     //
00351     // Good intersection! What about the normal? 
00352     //
00353     if (normSign*v.dot(normal) > 0)
00354     {
00355       G4double pr = p.perp();
00356       if (pr < DBL_MIN) pr = DBL_MIN;
00357       G4ThreeVector pNormal( rNorm*p.x()/pr, rNorm*p.y()/pr, zNorm );
00358       if (normSign*v.dot(pNormal) > 0)
00359       {
00360         G4double distOutside2;
00361         distFromSurface = -normSign*DistanceAway( p, false, distOutside2 );
00362         if (distOutside2 < surfTolerance*surfTolerance)
00363         {
00364           if (distFromSurface > -surfTolerance)
00365           {
00366             distance = s2;
00367             return true;
00368           }
00369         }
00370       }
00371       else 
00372         distFromSurface = s2;
00373       
00374       if (s2 > 0)
00375       {
00376         distance = s2;
00377         return true;
00378       }
00379     }
00380   }  
00381 
00382   //
00383   // Better luck next time
00384   //
00385   return false;
00386 }

G4ThreeVector G4PolyconeSide::Normal ( const G4ThreeVector p,
G4double bestDistance 
) [virtual]

Implements G4VCSGface.

Definition at line 468 of file G4PolyconeSide.cc.

References DistanceAway(), rNorm, and zNorm.

00470 {
00471   if (p == G4ThreeVector(0.,0.,0.))  { return p; }
00472 
00473   G4double dFrom, dOut2;
00474   
00475   dFrom = DistanceAway( p, false, dOut2 );
00476   
00477   *bestDistance = std::sqrt( dFrom*dFrom + dOut2 );
00478   
00479   G4double rds = p.perp();
00480   if (rds!=0.) { return G4ThreeVector(rNorm*p.x()/rds,rNorm*p.y()/rds,zNorm); }
00481   return G4ThreeVector( 0.,0., zNorm ).unit();
00482 }

G4PolyconeSide & G4PolyconeSide::operator= ( const G4PolyconeSide source  ) 

Definition at line 196 of file G4PolyconeSide.cc.

References cone, CopyStuff(), corners, and phiIsOpen.

00197 {
00198   if (this == &source)  { return *this; }
00199 
00200   delete cone;
00201   if (phiIsOpen)  { delete [] corners; }
00202   
00203   CopyStuff( source );
00204   
00205   return *this;
00206 }

G4bool G4PolyconeSide::PointOnCone ( const G4ThreeVector hit,
G4double  normSign,
const G4ThreeVector p,
const G4ThreeVector v,
G4ThreeVector normal 
) [protected]

Definition at line 993 of file G4PolyconeSide.cc.

References cone, corners, DBL_MIN, deltaPhi, GetPhi(), G4IntersectingCone::HitOn(), phiIsOpen, rNorm, startPhi, and zNorm.

Referenced by Intersect().

00998 {
00999   G4double rx = hit.perp();
01000   //
01001   // Check radial/z extent, as appropriate
01002   //
01003   if (!cone->HitOn( rx, hit.z() )) return false;
01004   
01005   if (phiIsOpen)
01006   {
01007     G4double phiTolerant = 2.0*kCarTolerance/(rx+kCarTolerance);
01008     //
01009     // Check phi segment. Here we have to be careful
01010     // to use the standard method consistent with
01011     // PolyPhiFace. See PolyPhiFace::InsideEdgesExact
01012     //
01013     G4double phi = GetPhi(hit);
01014     while( phi < startPhi-phiTolerant ) phi += twopi;
01015     
01016     if (phi > startPhi+deltaPhi+phiTolerant) return false;
01017     
01018     if (phi > startPhi+deltaPhi-phiTolerant)
01019     {
01020       //
01021       // Exact treatment
01022       //
01023       G4ThreeVector qx = p + v;
01024       G4ThreeVector qa = qx - corners[2],
01025               qb = qx - corners[3];
01026       G4ThreeVector qacb = qa.cross(qb);
01027       
01028       if (normSign*qacb.dot(v) < 0) return false;
01029     }
01030     else if (phi < phiTolerant)
01031     {
01032       G4ThreeVector qx = p + v;
01033       G4ThreeVector qa = qx - corners[1],
01034               qb = qx - corners[0];
01035       G4ThreeVector qacb = qa.cross(qb);
01036       
01037       if (normSign*qacb.dot(v) < 0) return false;
01038     }
01039   }
01040   
01041   //
01042   // We have a good hit! Calculate normal
01043   //
01044   if (rx < DBL_MIN) 
01045     normal = G4ThreeVector( 0, 0, zNorm < 0 ? -1 : 1 );
01046   else
01047     normal = G4ThreeVector( rNorm*hit.x()/rx, rNorm*hit.y()/rx, zNorm );
01048   return true;
01049 }

G4double G4PolyconeSide::SurfaceArea (  )  [virtual]

Implements G4VCSGface.

Definition at line 1087 of file G4PolyconeSide.cc.

References deltaPhi, r, sqr(), and z.

01088 { 
01089   if(fSurfaceArea==0)
01090   {
01091     fSurfaceArea = (r[0]+r[1])* std::sqrt(sqr(r[0]-r[1])+sqr(z[0]-z[1]));
01092     fSurfaceArea *= 0.5*(deltaPhi);
01093   }  
01094   return fSurfaceArea;
01095 }


Field Documentation

G4bool G4PolyconeSide::allBehind [protected]

Definition at line 138 of file G4PolyconeSide.hh.

Referenced by CopyStuff(), G4PolyconeSide(), and Intersect().

G4IntersectingCone* G4PolyconeSide::cone [protected]

Definition at line 140 of file G4PolyconeSide.hh.

Referenced by CopyStuff(), Extent(), G4PolyconeSide(), Intersect(), operator=(), PointOnCone(), and ~G4PolyconeSide().

G4ThreeVector* G4PolyconeSide::corners [protected]

Definition at line 154 of file G4PolyconeSide.hh.

Referenced by CopyStuff(), G4PolyconeSide(), operator=(), PointOnCone(), and ~G4PolyconeSide().

G4double G4PolyconeSide::deltaPhi [protected]

Definition at line 135 of file G4PolyconeSide.hh.

Referenced by CalculateExtent(), CopyStuff(), DistanceAway(), Extent(), G4PolyconeSide(), GetPointOnFace(), PointOnCone(), and SurfaceArea().

G4double G4PolyconeSide::length [protected]

Definition at line 144 of file G4PolyconeSide.hh.

Referenced by CopyStuff(), DistanceAway(), and G4PolyconeSide().

G4int G4PolyconeSide::ncorners [protected]

Definition at line 153 of file G4PolyconeSide.hh.

Referenced by CopyStuff(), and G4PolyconeSide().

G4double G4PolyconeSide::nextRS [protected]

Definition at line 147 of file G4PolyconeSide.hh.

Referenced by CalculateExtent(), CopyStuff(), and G4PolyconeSide().

G4double G4PolyconeSide::nextZS [protected]

Definition at line 147 of file G4PolyconeSide.hh.

Referenced by CalculateExtent(), CopyStuff(), and G4PolyconeSide().

G4bool G4PolyconeSide::phiIsOpen [protected]

Definition at line 137 of file G4PolyconeSide.hh.

Referenced by CalculateExtent(), CopyStuff(), DistanceAway(), Extent(), G4PolyconeSide(), operator=(), PointOnCone(), and ~G4PolyconeSide().

G4double G4PolyconeSide::prevRS [protected]

Definition at line 145 of file G4PolyconeSide.hh.

Referenced by CalculateExtent(), CopyStuff(), and G4PolyconeSide().

G4double G4PolyconeSide::prevZS [protected]

Definition at line 145 of file G4PolyconeSide.hh.

Referenced by CalculateExtent(), CopyStuff(), and G4PolyconeSide().

G4double G4PolyconeSide::r[2] [protected]

Definition at line 134 of file G4PolyconeSide.hh.

Referenced by CalculateExtent(), CopyStuff(), DistanceAway(), Extent(), G4PolyconeSide(), GetPointOnFace(), and SurfaceArea().

G4double G4PolyconeSide::rNorm [protected]

Definition at line 142 of file G4PolyconeSide.hh.

Referenced by CalculateExtent(), CopyStuff(), DistanceAway(), G4PolyconeSide(), Intersect(), Normal(), and PointOnCone().

G4double G4PolyconeSide::rNormEdge[2] [protected]

Definition at line 150 of file G4PolyconeSide.hh.

Referenced by CopyStuff(), DistanceAway(), and G4PolyconeSide().

G4double G4PolyconeSide::rS [protected]

Definition at line 143 of file G4PolyconeSide.hh.

Referenced by CalculateExtent(), CopyStuff(), DistanceAway(), and G4PolyconeSide().

G4double G4PolyconeSide::startPhi [protected]

Definition at line 135 of file G4PolyconeSide.hh.

Referenced by CalculateExtent(), CopyStuff(), DistanceAway(), Extent(), G4PolyconeSide(), GetPointOnFace(), and PointOnCone().

G4double G4PolyconeSide::z[2] [protected]

Definition at line 134 of file G4PolyconeSide.hh.

Referenced by CalculateExtent(), CopyStuff(), DistanceAway(), Extent(), G4PolyconeSide(), GetPointOnFace(), and SurfaceArea().

G4double G4PolyconeSide::zNorm [protected]

Definition at line 142 of file G4PolyconeSide.hh.

Referenced by CopyStuff(), DistanceAway(), G4PolyconeSide(), Intersect(), Normal(), and PointOnCone().

G4double G4PolyconeSide::zNormEdge[2] [protected]

Definition at line 150 of file G4PolyconeSide.hh.

Referenced by CopyStuff(), DistanceAway(), and G4PolyconeSide().

G4double G4PolyconeSide::zS [protected]

Definition at line 143 of file G4PolyconeSide.hh.

Referenced by CalculateExtent(), CopyStuff(), DistanceAway(), and G4PolyconeSide().


The documentation for this class was generated from the following files:
Generated on Mon May 27 17:52:57 2013 for Geant4 by  doxygen 1.4.7