#include <G4PolyhedraSide.hh>

Inheritance diagram for G4PolyhedraSide:

Data Structures
struct	sG4PolyhedraSideEdge

struct	sG4PolyhedraSideVec

Public Member Functions
	G4PolyhedraSide (const G4PolyhedraSideRZ prevRZ, const G4PolyhedraSideRZ tail, const G4PolyhedraSideRZ head, const G4PolyhedraSideRZ nextRZ, G4int numSide, G4double phiStart, G4double phiTotal, G4bool phiIsOpen, G4bool isAllBehind=false)

virtual	~G4PolyhedraSide ()

	G4PolyhedraSide (const G4PolyhedraSide &source)

G4PolyhedraSide &	operator= (const G4PolyhedraSide &source)

G4bool	Intersect (const G4ThreeVector &p, const G4ThreeVector &v, G4bool outgoing, G4double surfTolerance, G4double &distance, G4double &distFromSurface, G4ThreeVector &normal, G4bool &allBehind)

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)

G4VCSGface *	Clone ()

G4double	SurfaceTriangle (G4ThreeVector p1, G4ThreeVector p2, G4ThreeVector p3, G4ThreeVector *p4)

G4ThreeVector	GetPointOnPlane (G4ThreeVector p0, G4ThreeVector p1, G4ThreeVector p2, G4ThreeVector p3, G4double *Area)

G4double	SurfaceArea ()

G4ThreeVector	GetPointOnFace ()

	G4PolyhedraSide (__void__ &)

G4int	GetInstanceID () const

Public Member Functions inherited from G4VCSGface
	G4VCSGface ()

virtual	~G4VCSGface ()

Static Public Member Functions
static const G4PhSideManager &	GetSubInstanceManager ()

Protected Types
typedef struct G4PolyhedraSide::sG4PolyhedraSideEdge	G4PolyhedraSideEdge

typedef struct G4PolyhedraSide::sG4PolyhedraSideVec	G4PolyhedraSideVec

Protected Member Functions
G4bool	IntersectSidePlane (const G4ThreeVector &p, const G4ThreeVector &v, const G4PolyhedraSideVec &vec, G4double normSign, G4double surfTolerance, G4double &distance, G4double &distFromSurface)

G4int	LineHitsSegments (const G4ThreeVector &p, const G4ThreeVector &v, G4int i1, G4int i2)

G4int	ClosestPhiSegment (G4double phi)

G4int	PhiSegment (G4double phi)

G4double	GetPhi (const G4ThreeVector &p)

G4double	DistanceToOneSide (const G4ThreeVector &p, const G4PolyhedraSideVec &vec, G4double *normDist)

G4double	DistanceAway (const G4ThreeVector &p, const G4PolyhedraSideVec &vec, G4double *normDist)

void	CopyStuff (const G4PolyhedraSide &source)

Protected Attributes
G4int	numSide

G4double	r [2]

G4double	z [2]

G4double	startPhi

G4double	deltaPhi

G4double	endPhi

G4bool	phiIsOpen

G4bool	allBehind

G4IntersectingCone *	cone

G4PolyhedraSideVec *	vecs

G4PolyhedraSideEdge *	edges

G4double	lenRZ

G4double	lenPhi [2]

G4double	edgeNorm

Friends
struct	sG4PolyhedraSideVec

Detailed Description

Definition at line 102 of file G4PolyhedraSide.hh.

Member Typedef Documentation

typedef struct G4PolyhedraSide::sG4PolyhedraSideEdge G4PolyhedraSide::G4PolyhedraSideEdge

protected

typedef struct G4PolyhedraSide::sG4PolyhedraSideVec G4PolyhedraSide::G4PolyhedraSideVec

protected

Constructor & Destructor Documentation

G4PolyhedraSide::G4PolyhedraSide	(	const G4PolyhedraSideRZ *	prevRZ,
		const G4PolyhedraSideRZ *	tail,
		const G4PolyhedraSideRZ *	head,
		const G4PolyhedraSideRZ *	nextRZ,
		G4int	numSide,
		G4double	phiStart,
		G4double	phiTotal,
		G4bool	phiIsOpen,
		G4bool	isAllBehind = `false`
	)

Definition at line 74 of file G4PolyhedraSide.cc.

Referenced by Clone().

 {
 
   instanceID = subInstanceManager.CreateSubInstance();
 
   kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
   fSurfaceArea=0.;
   G4MT_phphi.first = G4ThreeVector(0,0,0);
   G4MT_phphi.second= 0.0;
 
   //
   // Record values
   //
   r[0] = tail->r; z[0] = tail->z;
   r[1] = head->r; z[1] = head->z;
   
   G4double phiTotal;
   
   //
   // Set phi to our convention
   //
   startPhi = thePhiStart;
   while (startPhi < 0.0) startPhi += twopi;
   
   phiIsOpen = thePhiIsOpen;
   phiTotal = (phiIsOpen) ? thePhiTotal : twopi;
   
   allBehind = isAllBehind;
     
   //
   // Make our intersecting cone
   //
   cone = new G4IntersectingCone( r, z );
   
   //
   // Construct side plane vector set
   //
   numSide = theNumSide;
   deltaPhi = phiTotal/theNumSide;
   endPhi = startPhi+phiTotal;
   
   vecs = new G4PolyhedraSideVec[numSide];
   
   edges = new G4PolyhedraSideEdge[phiIsOpen ? numSide+1 : numSide];
   
   //
   // ...this is where we start
   //
   G4double phi = startPhi;
   G4ThreeVector a1( r[0]*std::cos(phi), r[0]*std::sin(phi), z[0] ),
           b1( r[1]*std::cos(phi), r[1]*std::sin(phi), z[1] ),
           c1( prevRZ->r*std::cos(phi), prevRZ->r*std::sin(phi), prevRZ->z ),
           d1( nextRZ->r*std::cos(phi), nextRZ->r*std::sin(phi), nextRZ->z ),
           a2, b2, c2, d2;
   G4PolyhedraSideEdge *edge = edges;
           
   G4PolyhedraSideVec *vec = vecs;
   do
   {
     //
     // ...this is where we are going
     //
     phi += deltaPhi;
     a2 = G4ThreeVector( r[0]*std::cos(phi), r[0]*std::sin(phi), z[0] );
     b2 = G4ThreeVector( r[1]*std::cos(phi), r[1]*std::sin(phi), z[1] );
     c2 = G4ThreeVector( prevRZ->r*std::cos(phi), prevRZ->r*std::sin(phi), prevRZ->z );
     d2 = G4ThreeVector( nextRZ->r*std::cos(phi), nextRZ->r*std::sin(phi), nextRZ->z );
     
     G4ThreeVector tt;  
     
     //
     // ...build some relevant vectors.
     //    the point is to sacrifice a little memory with precalcs 
     //    to gain speed
     //
     vec->center = 0.25*( a1 + a2 + b1 + b2 );
     
     tt = b2 + b1 - a2 - a1;
     vec->surfRZ = tt.unit();
     if (vec==vecs) lenRZ = 0.25*tt.mag();
     
     tt = b2 - b1 + a2 - a1;
     vec->surfPhi = tt.unit();
     if (vec==vecs)
     {
       lenPhi[0] = 0.25*tt.mag();
       tt = b2 - b1;
       lenPhi[1] = (0.5*tt.mag()-lenPhi[0])/lenRZ;
     }
     
     tt = vec->surfPhi.cross(vec->surfRZ);
     vec->normal = tt.unit();
     
     //
     // ...edge normals are the average of the normals of
     //    the two faces they connect.
     //
     // ...edge normals are necessary if we are to accurately
     //    decide if a point is "inside" a face. For non-convex
     //    shapes, it is absolutely necessary to know information
     //    on adjacent faces to accurate determine this.
     //
     // ...we don't need them for the phi edges, since that
     //    information is taken care of internally. The r/z edges,
     //    however, depend on the adjacent G4PolyhedraSide.
     //
     G4ThreeVector a12, adj;
     
     a12 = a2-a1;
 
     adj = 0.5*(c1+c2-a1-a2);
     adj = adj.cross(a12);  
     adj = adj.unit() + vec->normal;       
     vec->edgeNorm[0] = adj.unit();
     
     a12 = b1-b2;
     adj = 0.5*(d1+d2-b1-b2);
     adj = adj.cross(a12);  
     adj = adj.unit() + vec->normal;       
     vec->edgeNorm[1] = adj.unit();
     
     //
     // ...the corners are crucial. It is important that
     //    they are calculated consistently for adjacent
     //    G4PolyhedraSides, to avoid gaps caused by roundoff.
     //
     vec->edges[0] = edge;
     edge->corner[0] = a1;
     edge->corner[1] = b1;
     edge++;
     vec->edges[1] = edge;
 
     a1 = a2;
     b1 = b2;
     c1 = c2;
     d1 = d2;
   } while( ++vec < vecs+numSide );
   
   //
   // Clean up hanging edge
   //
   if (phiIsOpen)
   {
     edge->corner[0] = a2;
     edge->corner[1] = b2;
   }
   else
   {
     vecs[numSide-1].edges[1] = edges;
   }
   
   //
   // Go back and fill in remaining fields in edges
   //
   vec = vecs;
   G4PolyhedraSideVec *prev = vecs+numSide-1;
   do
   {
     edge = vec->edges[0];    // The edge between prev and vec
     
     //
     // Okay: edge normal is average of normals of adjacent faces
     //
     G4ThreeVector eNorm = vec->normal + prev->normal;
     edge->normal = eNorm.unit();  
     
     //
     // Vertex normal is average of norms of adjacent surfaces (all four)
     // However, vec->edgeNorm is unit vector in some direction
     // as the sum of normals of adjacent PolyhedraSide with vec.
     // The normalization used for this vector should be the same
     // for vec and prev.
     //
     eNorm = vec->edgeNorm[0] + prev->edgeNorm[0];
     edge->cornNorm[0] = eNorm.unit();
   
     eNorm = vec->edgeNorm[1] + prev->edgeNorm[1];
     edge->cornNorm[1] = eNorm.unit();
   } while( prev=vec, ++vec < vecs + numSide );
   
   if (phiIsOpen)
   {
     // G4double rFact = std::cos(0.5*deltaPhi);
     //
     // If phi is open, we need to patch up normals of the
     // first and last edges and their corresponding
     // vertices.
     //
     // We use vectors that are in the plane of the
     // face. This should be safe.
     //
     vec = vecs;
     
     G4ThreeVector normvec = vec->edges[0]->corner[0]
                           - vec->edges[0]->corner[1];
     normvec = normvec.cross(vec->normal);
     if (normvec.dot(vec->surfPhi) > 0) normvec = -normvec;
 
     vec->edges[0]->normal = normvec.unit();
     
     vec->edges[0]->cornNorm[0] = (vec->edges[0]->corner[0]
                                 - vec->center).unit();
     vec->edges[0]->cornNorm[1] = (vec->edges[0]->corner[1]
                                 - vec->center).unit();
     
     //
     // Repeat for ending phi
     //
     vec = vecs + numSide - 1;
     
     normvec = vec->edges[1]->corner[0] - vec->edges[1]->corner[1];
     normvec = normvec.cross(vec->normal);
     if (normvec.dot(vec->surfPhi) < 0) normvec = -normvec;
 
     vec->edges[1]->normal = normvec.unit();
     
     vec->edges[1]->cornNorm[0] = (vec->edges[1]->corner[0]
                                 - vec->center).unit();
     vec->edges[1]->cornNorm[1] = (vec->edges[1]->corner[1]
                                 - vec->center).unit();
   }
   
   //
   // edgeNorm is the factor one multiplies the distance along vector phi
   // on the surface of one of our sides in order to calculate the distance
   // from the edge. (see routine DistanceAway)
   //
   edgeNorm = 1.0/std::sqrt( 1.0 + lenPhi[1]*lenPhi[1] );
 }

G4PolyhedraSide::~G4PolyhedraSide ( )

virtual

Definition at line 331 of file G4PolyhedraSide.cc.

References cone, edges, and vecs.

 {
   delete cone;
   delete [] vecs;
   delete [] edges;
 }

G4PolyhedraSide::G4PolyhedraSide ( const G4PolyhedraSide & source )

Definition at line 342 of file G4PolyhedraSide.cc.

References CopyStuff(), and G4GeomSplitter< T >::CreateSubInstance().

   : G4VCSGface()
 {
   instanceID = subInstanceManager.CreateSubInstance();
 
   CopyStuff( source );
 }

G4PolyhedraSide::G4PolyhedraSide ( __void__ & )

Definition at line 317 of file G4PolyhedraSide.cc.

References lenPhi, r, and z.

   : numSide(0), startPhi(0.), deltaPhi(0.), endPhi(0.),
     phiIsOpen(false), allBehind(false), cone(0), vecs(0), edges(0),
     lenRZ(0.), edgeNorm(0.), kCarTolerance(0.), fSurfaceArea(0.), instanceID(0)
 {
   r[0] = r[1] = 0.;
   z[0] = z[1] = 0.;
   lenPhi[0] = lenPhi[1] = 0.;
 }

Member Function Documentation

void G4PolyhedraSide::CalculateExtent	(	const EAxis	axis,
		const G4VoxelLimits &	voxelLimit,
		const G4AffineTransform &	tranform,
		G4SolidExtentList &	extentList
	)

virtual

Implements G4VCSGface.

Definition at line 731 of file G4PolyhedraSide.cc.

References G4SolidExtentList::AddSurface(), G4ClippablePolygon::AddVertexInOrder(), G4PolyhedraSide::sG4PolyhedraSideEdge::corner, G4PolyhedraSide::sG4PolyhedraSideVec::edges, G4PolyhedraSide::sG4PolyhedraSideVec::normal, numSide, G4ClippablePolygon::PartialClip(), G4ClippablePolygon::SetNormal(), G4AffineTransform::TransformAxis(), and vecs.

 {
   //
   // Loop over all sides
   //
   G4PolyhedraSideVec *vec = vecs;
   do
   {
     //
     // Fill our polygon with the four corners of
     // this side, after the specified transformation
     //
     G4ClippablePolygon polygon;
     
     polygon.AddVertexInOrder(transform.
                              TransformPoint(vec->edges[0]->corner[0]));
     polygon.AddVertexInOrder(transform.
                              TransformPoint(vec->edges[0]->corner[1]));
     polygon.AddVertexInOrder(transform.
                              TransformPoint(vec->edges[1]->corner[1]));
     polygon.AddVertexInOrder(transform.
                              TransformPoint(vec->edges[1]->corner[0]));
     
     //
     // Get extent
     //  
     if (polygon.PartialClip( voxelLimit, axis ))
     {
       //
       // Get dot product of normal along target axis
       //
       polygon.SetNormal( transform.TransformAxis(vec->normal) );
 
       extentList.AddSurface( polygon );
     }
   } while( ++vec < vecs+numSide );
   
   return;
 }

G4VCSGface* G4PolyhedraSide::Clone ( )

inlinevirtual

Implements G4VCSGface.

Definition at line 138 of file G4PolyhedraSide.hh.

References G4PolyhedraSide().

138 { return new G4PolyhedraSide( *this ); }

G4PolyhedraSide::G4PolyhedraSide

G4PolyhedraSide(const G4PolyhedraSideRZ *prevRZ, const G4PolyhedraSideRZ *tail, const G4PolyhedraSideRZ *head, const G4PolyhedraSideRZ *nextRZ, G4int numSide, G4double phiStart, G4double phiTotal, G4bool phiIsOpen, G4bool isAllBehind=false)

Definition: G4PolyhedraSide.cc:74

G4int G4PolyhedraSide::ClosestPhiSegment ( G4double phi )

protected

Definition at line 943 of file G4PolyhedraSide.cc.

References endPhi, numSide, PhiSegment(), startPhi, and python.hepunit::twopi.

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

 {
   G4int iPhi = PhiSegment( phi0 );
   if (iPhi >= 0) return iPhi;
   
   //
   // Boogers! The points falls inside the phi segment.
   // Look for the closest point: the start, or  end
   //
   G4double phi = phi0;
   
   while( phi < startPhi ) phi += twopi;
   G4double d1 = phi-endPhi;
 
   while( phi > startPhi ) phi -= twopi;
   G4double d2 = startPhi-phi;
   
   return (d2 < d1) ? 0 : numSide-1;
 }

void G4PolyhedraSide::CopyStuff ( const G4PolyhedraSide & source )

protected

Definition at line 371 of file G4PolyhedraSide.cc.

References allBehind, cone, deltaPhi, edgeNorm, G4PolyhedraSide::sG4PolyhedraSideVec::edges, edges, endPhi, lenPhi, lenRZ, numSide, phiIsOpen, r, startPhi, vecs, and z.

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

 {
   //
   // The simple stuff
   //
   numSide    = source.numSide;
   r[0]    = source.r[0];
   r[1]    = source.r[1];
   z[0]    = source.z[0];
   z[1]    = source.z[1];
   startPhi  = source.startPhi;
   deltaPhi  = source.deltaPhi;
   endPhi    = source.endPhi;
   phiIsOpen = source.phiIsOpen;
   allBehind = source.allBehind;
   
   lenRZ     = source.lenRZ;
   lenPhi[0] = source.lenPhi[0];
   lenPhi[1] = source.lenPhi[1];
   edgeNorm  = source.edgeNorm;
 
   kCarTolerance = source.kCarTolerance;
   fSurfaceArea = source.fSurfaceArea;
 
   cone = new G4IntersectingCone( *source.cone );
 
   //
   // Duplicate edges
   //
   G4int  numEdges = phiIsOpen ? numSide+1 : numSide;
   edges = new G4PolyhedraSideEdge[numEdges];
   
   G4PolyhedraSideEdge *edge = edges,
           *sourceEdge = source.edges;
   do
   {
     *edge = *sourceEdge;
   } while( ++sourceEdge, ++edge < edges + numEdges);
 
   //
   // Duplicate vecs
   //
   vecs = new G4PolyhedraSideVec[numSide];
   
   G4PolyhedraSideVec *vec = vecs,
          *sourceVec = source.vecs;
   do
   {
     *vec = *sourceVec;
     vec->edges[0] = edges + (sourceVec->edges[0] - source.edges);
     vec->edges[1] = edges + (sourceVec->edges[1] - source.edges);
   } while( ++sourceVec, ++vec < vecs + numSide );
 }

G4double G4PolyhedraSide::Distance	(	const G4ThreeVector &	p,
		G4bool	outgoing
	)

virtual

Implements G4VCSGface.

Definition at line 586 of file G4PolyhedraSide.cc.

References G4PolyhedraSide::sG4PolyhedraSideVec::center, ClosestPhiSegment(), DistanceAway(), CLHEP::Hep3Vector::dot(), GetPhi(), and vecs.

 {
   G4double normSign = outgoing ? -1 : +1;
   
   //
   // Try the closest phi segment first
   //
   G4int iPhi = ClosestPhiSegment( GetPhi(p) );
   
   G4ThreeVector pdotc = p - vecs[iPhi].center;
   G4double normDist = pdotc.dot(vecs[iPhi].normal);
   
   if (normSign*normDist > -0.5*kCarTolerance)
   {
     return DistanceAway( p, vecs[iPhi], &normDist );
   }
 
   //
   // Now we have an interesting problem... do we try to find the
   // closest facing side??
   //
   // Considered carefully, the answer is no. We know that if we
   // are asking for the distance out, we are supposed to be inside,
   // and vice versa.
   //
   
   return kInfinity;
 }

G4double G4PolyhedraSide::DistanceAway	(	const G4ThreeVector &	p,
		const G4PolyhedraSideVec &	vec,
		G4double *	normDist
	)

protected

Definition at line 1060 of file G4PolyhedraSide.cc.

References G4PolyhedraSide::sG4PolyhedraSideVec::center, G4PolyhedraSide::sG4PolyhedraSideEdge::corner, G4PolyhedraSide::sG4PolyhedraSideEdge::cornNorm, CLHEP::Hep3Vector::dot(), G4PolyhedraSide::sG4PolyhedraSideVec::edgeNorm, edgeNorm, G4PolyhedraSide::sG4PolyhedraSideVec::edges, lenPhi, lenRZ, G4PolyhedraSide::sG4PolyhedraSideEdge::normal, G4PolyhedraSide::sG4PolyhedraSideVec::surfPhi, and G4PolyhedraSide::sG4PolyhedraSideVec::surfRZ.

Referenced by Distance(), and DistanceToOneSide().

 {
   G4double distOut2;
   G4ThreeVector pct = p - vec.center;
   G4double distFaceNorm = *normDist;
   
   //
   // Okay, are we inside bounds?
   //
   G4double pcDotRZ  = pct.dot(vec.surfRZ);
   G4double pcDotPhi = pct.dot(vec.surfPhi);
   
   //
   // Go through all permutations.
   //                                                   Phi
   //               |              |                     ^
   //           B   |      H       |   E                 |
   //        ------[1]------------[3]-----               |
   //               |XXXXXXXXXXXXXX|                     +----> RZ
   //           C   |XXXXXXXXXXXXXX|   F
   //               |XXXXXXXXXXXXXX|
   //        ------[0]------------[2]----
   //           A   |      G       |   D
   //               |              |
   //
   // It's real messy, but at least it's quick
   //
   
   if (pcDotRZ < -lenRZ)
   {
     G4double lenPhiZ = lenPhi[0] - lenRZ*lenPhi[1];
     G4double distOutZ = pcDotRZ+lenRZ;
     //
     // Below in RZ
     //
     if (pcDotPhi < -lenPhiZ)
     {
       //
       // ...and below in phi. Find distance to point (A)
       //
       G4double distOutPhi = pcDotPhi+lenPhiZ;
       distOut2 = distOutPhi*distOutPhi + distOutZ*distOutZ;
       G4ThreeVector pa = p - vec.edges[0]->corner[0];
       *normDist = pa.dot(vec.edges[0]->cornNorm[0]);
     }
     else if (pcDotPhi > lenPhiZ)
     {
       //
       // ...and above in phi. Find distance to point (B)
       //
       G4double distOutPhi = pcDotPhi-lenPhiZ;
       distOut2 = distOutPhi*distOutPhi + distOutZ*distOutZ;
       G4ThreeVector pb = p - vec.edges[1]->corner[0];
       *normDist = pb.dot(vec.edges[1]->cornNorm[0]);
     }
     else
     {
       //
       // ...and inside in phi. Find distance to line (C)
       //
       G4ThreeVector pa = p - vec.edges[0]->corner[0];
       distOut2 = distOutZ*distOutZ;
       *normDist = pa.dot(vec.edgeNorm[0]);
     }
   }
   else if (pcDotRZ > lenRZ)
   {
     G4double lenPhiZ = lenPhi[0] + lenRZ*lenPhi[1];
     G4double distOutZ = pcDotRZ-lenRZ;
     //
     // Above in RZ
     //
     if (pcDotPhi < -lenPhiZ)
     {
       //
       // ...and below in phi. Find distance to point (D)
       //
       G4double distOutPhi = pcDotPhi+lenPhiZ;
       distOut2 = distOutPhi*distOutPhi + distOutZ*distOutZ;
       G4ThreeVector pd = p - vec.edges[0]->corner[1];
       *normDist = pd.dot(vec.edges[0]->cornNorm[1]);
     }
     else if (pcDotPhi > lenPhiZ)
     {
       //
       // ...and above in phi. Find distance to point (E)
       //
       G4double distOutPhi = pcDotPhi-lenPhiZ;
       distOut2 = distOutPhi*distOutPhi + distOutZ*distOutZ;
       G4ThreeVector pe = p - vec.edges[1]->corner[1];
       *normDist = pe.dot(vec.edges[1]->cornNorm[1]);
     }
     else
     {
       //
       // ...and inside in phi. Find distance to line (F)
       //
       distOut2 = distOutZ*distOutZ;
       G4ThreeVector pd = p - vec.edges[0]->corner[1];
       *normDist = pd.dot(vec.edgeNorm[1]);
     }
   }
   else
   {
     G4double lenPhiZ = lenPhi[0] + pcDotRZ*lenPhi[1];
     //
     // We are inside RZ bounds
     // 
     if (pcDotPhi < -lenPhiZ)
     {
       //
       // ...and below in phi. Find distance to line (G)
       //
       G4double distOut = edgeNorm*(pcDotPhi+lenPhiZ);
       distOut2 = distOut*distOut;
       G4ThreeVector pd = p - vec.edges[0]->corner[1];
       *normDist = pd.dot(vec.edges[0]->normal);
     }
     else if (pcDotPhi > lenPhiZ)
     {
       //
       // ...and above in phi. Find distance to line (H)
       //
       G4double distOut = edgeNorm*(pcDotPhi-lenPhiZ);
       distOut2 = distOut*distOut;
       G4ThreeVector pe = p - vec.edges[1]->corner[1];
       *normDist = pe.dot(vec.edges[1]->normal);
     }
     else
     {
       //
       // Inside bounds! No penalty.
       //
       return std::fabs(distFaceNorm);
     }
   }
   return std::sqrt( distFaceNorm*distFaceNorm + distOut2 );
 }

G4double G4PolyhedraSide::DistanceToOneSide	(	const G4ThreeVector &	p,
		const G4PolyhedraSideVec &	vec,
		G4double *	normDist
	)

protected

Definition at line 1036 of file G4PolyhedraSide.cc.

References G4PolyhedraSide::sG4PolyhedraSideVec::center, DistanceAway(), CLHEP::Hep3Vector::dot(), and G4PolyhedraSide::sG4PolyhedraSideVec::normal.

Referenced by Inside(), and Normal().

 {
   G4ThreeVector pct = p - vec.center;
   
   //
   // Get normal distance
   //
   *normDist = vec.normal.dot(pct);
 
   //
   // Add edge penalty
   //
   return DistanceAway( p, vec, normDist );
 }

G4double G4PolyhedraSide::Extent ( const G4ThreeVector axis )

virtual

Implements G4VCSGface.

Definition at line 671 of file G4PolyhedraSide.cc.

References cone, G4PolyhedraSide::sG4PolyhedraSideEdge::corner, DBL_MIN, CLHEP::Hep3Vector::dot(), G4PolyhedraSide::sG4PolyhedraSideVec::edges, GetPhi(), numSide, CLHEP::Hep3Vector::perp2(), PhiSegment(), vecs, CLHEP::Hep3Vector::z(), G4IntersectingCone::ZHi(), and G4IntersectingCone::ZLo().

 {
   if (axis.perp2() < DBL_MIN)
   {
     //
     // Special case
     //
     return axis.z() < 0 ? -cone->ZLo() : cone->ZHi();
   }
 
   G4int iPhi, i1, i2;
   G4double best;
   G4ThreeVector *list[4];
   
   //
   // Which phi segment, if any, does the axis belong to
   //
   iPhi = PhiSegment( GetPhi(axis) );
   
   if (iPhi < 0)
   {
     //
     // No phi segment? Check front edge of first side and
     // last edge of second side
     //
     i1 = 0; i2 = numSide-1;
   }
   else
   {
     //
     // Check all corners of matching phi side
     //
     i1 = iPhi; i2 = iPhi;
   }
   
   list[0] = vecs[i1].edges[0]->corner;
   list[1] = vecs[i1].edges[0]->corner+1;
   list[2] = vecs[i2].edges[1]->corner;
   list[3] = vecs[i2].edges[1]->corner+1;
         
   //
   // Who's biggest?
   //
   best = -kInfinity;
   G4ThreeVector **vec = list;
   do
   {
     G4double answer = (*vec)->dot(axis);
     if (answer > best) best = answer;
   } while( ++vec < list+4 );
   
   return best;
 }

G4int G4PolyhedraSide::GetInstanceID ( ) const

inline

Definition at line 161 of file G4PolyhedraSide.hh.

161 { return instanceID; }

G4double G4PolyhedraSide::GetPhi ( const G4ThreeVector & p )

protected

Definition at line 1009 of file G4PolyhedraSide.cc.

References G4MT_phphi, and CLHEP::Hep3Vector::phi().

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

 {
   G4double val=0.;
 
   if (G4MT_phphi.first != p)
   {
     val = p.phi();
     G4MT_phphi.first = p;
     G4MT_phphi.second = val;
   }
   else
   {
     val = G4MT_phphi.second;
   }
   return val;
 }

G4ThreeVector G4PolyhedraSide::GetPointOnFace ( )

virtual

Implements G4VCSGface.

Definition at line 1286 of file G4PolyhedraSide.cc.

References G4PolyhedraSide::sG4PolyhedraSideEdge::corner, G4PolyhedraSide::sG4PolyhedraSideVec::edges, G4UniformRand, GetPointOnPlane(), numSide, and vecs.

 {
   // Define the variables
   //
   std::vector<G4double>areas;
   std::vector<G4ThreeVector>points;
   G4double area=0;
   G4double result1;
   G4ThreeVector point1;
   G4ThreeVector v1,v2,v3,v4; 
   G4PolyhedraSideVec *vec = vecs;
 
   // Do a loop on all SideEdge
   //
   do
   {
     // Define 4points for a Plane or Triangle
     //
     v1=vec->edges[0]->corner[0];
     v2=vec->edges[0]->corner[1];
     v3=vec->edges[1]->corner[1];
     v4=vec->edges[1]->corner[0];
     point1=GetPointOnPlane(v1,v2,v3,v4,&result1);
     points.push_back(point1);
     areas.push_back(result1);
     area+=result1;
   } while( ++vec < vecs+numSide );
 
   // Choose randomly one of the surfaces and point on it
   //
   G4double chose = area*G4UniformRand();
   G4double Achose1,Achose2;
   Achose1=0;Achose2=0.; 
   G4int i=0;
   do 
   {
     Achose2+=areas[i];
     if(chose>=Achose1 && chose<Achose2)
     {
       point1=points[i] ; break;     
     }
     i++; Achose1=Achose2;
   } while( i<numSide );
  
   return point1;
 }

G4ThreeVector G4PolyhedraSide::GetPointOnPlane	(	G4ThreeVector	p0,
		G4ThreeVector	p1,
		G4ThreeVector	p2,
		G4ThreeVector	p3,
		G4double *	Area
	)

Definition at line 1229 of file G4PolyhedraSide.cc.

References G4UniformRand, and SurfaceTriangle().

Referenced by GetPointOnFace(), and SurfaceArea().

 {
   G4double chose,aOne,aTwo;
   G4ThreeVector point1,point2;
   aOne = SurfaceTriangle(p0,p1,p2,&point1);
   aTwo = SurfaceTriangle(p2,p3,p0,&point2);
   *Area= aOne+aTwo;
 
   chose = G4UniformRand()*(aOne+aTwo);
   if( (chose>=0.) && (chose < aOne) )
   {
    return (point1);    
   }
   return (point2);
 }

const G4PhSideManager & G4PolyhedraSide::GetSubInstanceManager ( )

static

Definition at line 63 of file G4PolyhedraSide.cc.

Referenced by G4SolidsWorkspace::G4SolidsWorkspace().

 {
   return subInstanceManager;
 }

EInside G4PolyhedraSide::Inside	(	const G4ThreeVector &	p,
		G4double	tolerance,
		G4double *	bestDistance
	)

virtual

Implements G4VCSGface.

Definition at line 619 of file G4PolyhedraSide.cc.

References ClosestPhiSegment(), DistanceToOneSide(), GetPhi(), kInside, kOutside, kSurface, and vecs.

 {
   //
   // Which phi segment is closest to this point?
   //
   G4int iPhi = ClosestPhiSegment( GetPhi(p) );
   
   G4double norm;
   
   //
   // Get distance to this segment
   //
   *bestDistance = DistanceToOneSide( p, vecs[iPhi], &norm );
   
   //
   // Use distance along normal to decide return value
   //
   if ( (std::fabs(norm) < tolerance) && (*bestDistance < 2.0*tolerance) )
     return kSurface;
   else if (norm < 0)
     return kInside;
   else  
     return kOutside;
 }

G4bool G4PolyhedraSide::Intersect	(	const G4ThreeVector &	p,
		const G4ThreeVector &	v,
		G4bool	outgoing,
		G4double	surfTolerance,
		G4double &	distance,
		G4double &	distFromSurface,
		G4ThreeVector &	normal,
		G4bool &	allBehind
	)

virtual

Implements G4VCSGface.

Definition at line 467 of file G4PolyhedraSide.cc.

References allBehind, G4PolyhedraSide::sG4PolyhedraSideVec::center, G4PolyhedraSide::sG4PolyhedraSideEdge::corner, CLHEP::Hep3Vector::cross(), CLHEP::Hep3Vector::dot(), G4PolyhedraSide::sG4PolyhedraSideVec::edges, lenPhi, lenRZ, G4PolyhedraSide::sG4PolyhedraSideVec::normal, numSide, G4InuclParticleNames::pp, r, G4PolyhedraSide::sG4PolyhedraSideVec::surfPhi, G4PolyhedraSide::sG4PolyhedraSideVec::surfRZ, test::v, and vecs.

 {
   G4double normSign = outgoing ? +1 : -1;
   
   //
   // ------------------TO BE IMPLEMENTED---------------------
   // Testing the intersection of individual phi faces is
   // pretty straight forward. The simple thing therefore is to
   // form a loop and check them all in sequence.
   //
   // But, I worry about one day someone making
   // a polygon with a thousands sides. A linear search
   // would not be ideal in such a case.
   //
   // So, it would be nice to be able to quickly decide
   // which face would be intersected. One can make a very
   // good guess by using the intersection with a cone.
   // However, this is only reliable in 99% of the cases.
   //
   // My solution: make a decent guess as to the one or
   // two potential faces might get intersected, and then
   // test them. If we have the wrong face, use the test
   // to make a better guess.
   //
   // Since we might have two guesses, form a queue of
   // potential intersecting faces. Keep an array of 
   // already tested faces to avoid doing one more than
   // once.
   //
   // Result: at worst, an iterative search. On average,
   // a little more than two tests would be required.
   //
   G4ThreeVector q = p + v;
   
   G4int face = 0;
   G4PolyhedraSideVec *vec = vecs;
   do
   {
     //
     // Correct normal?
     //
     G4double dotProd = normSign*v.dot(vec->normal);
     if (dotProd <= 0) continue;
   
     //
     // Is this face in front of the point along the trajectory?
     //
     G4ThreeVector delta = p - vec->center;
     distFromSurface = -normSign*delta.dot(vec->normal);
     
     if (distFromSurface < -surfTolerance) continue;
     
     //
     //                            phi
     //      c -------- d           ^
     //      |          |           |
     //      a -------- b           +---> r/z
     //
     //
     // Do we remain on this particular segment?
     //
     G4ThreeVector qc = q - vec->edges[1]->corner[0];
     G4ThreeVector qd = q - vec->edges[1]->corner[1];
     
     if (normSign*qc.cross(qd).dot(v) < 0) continue;
     
     G4ThreeVector qa = q - vec->edges[0]->corner[0];
     G4ThreeVector qb = q - vec->edges[0]->corner[1];
     
     if (normSign*qa.cross(qb).dot(v) > 0) continue;
     
     //
     // We found the one and only segment we might be intersecting.
     // Do we remain within r/z bounds?
     //
     
     if (r[0] > 1/kInfinity && normSign*qa.cross(qc).dot(v) < 0) return false;
     if (r[1] > 1/kInfinity && normSign*qb.cross(qd).dot(v) > 0) return false;
     
     //
     // We allow the face to be slightly behind the trajectory
     // (surface tolerance) only if the point p is within
     // the vicinity of the face
     //
     if (distFromSurface < 0)
     {
       G4ThreeVector ps = p - vec->center; 
       
       G4double rz = ps.dot(vec->surfRZ);
       if (std::fabs(rz) > lenRZ+surfTolerance) return false; 
 
       G4double pp = ps.dot(vec->surfPhi);
       if (std::fabs(pp) > lenPhi[0] + lenPhi[1]*rz + surfTolerance) return false;
     }
       
 
     //
     // Intersection found. Return answer.
     //
     distance = distFromSurface/dotProd;
     normal = vec->normal;
     isAllBehind = allBehind;
     return true;
   } while( ++vec, ++face < numSide );
 
   //
   // Oh well. Better luck next time.
   //
   return false;
 }

G4bool G4PolyhedraSide::IntersectSidePlane	(	const G4ThreeVector &	p,
		const G4ThreeVector &	v,
		const G4PolyhedraSideVec &	vec,
		G4double	normSign,
		G4double	surfTolerance,
		G4double &	distance,
		G4double &	distFromSurface
	)

protected

Definition at line 802 of file G4PolyhedraSide.cc.

References G4PolyhedraSide::sG4PolyhedraSideVec::center, G4PolyhedraSide::sG4PolyhedraSideEdge::corner, CLHEP::Hep3Vector::cross(), CLHEP::Hep3Vector::dot(), G4PolyhedraSide::sG4PolyhedraSideVec::edges, lenRZ, G4PolyhedraSide::sG4PolyhedraSideVec::normal, r, G4PolyhedraSide::sG4PolyhedraSideVec::surfRZ, and test::v.

 {
   //
   // Correct normal? Here we have straight sides, and can safely ignore
   // intersections where the dot product with the normal is zero.
   //
   G4double dotProd = normSign*v.dot(vec.normal);
   
   if (dotProd <= 0) return false;
   
   //
   // Calculate distance to surface. If the side is too far
   // behind the point, we must reject it.
   //
   G4ThreeVector delta = p - vec.center;
   distFromSurface = -normSign*delta.dot(vec.normal);
     
   if (distFromSurface < -surfTolerance) return false;
 
   //
   // Calculate precise distance to intersection with the side
   // (along the trajectory, not normal to the surface)
   //
   distance = distFromSurface/dotProd;
   
   //
   // Do we fall off the r/z extent of the segment?
   //
   // Calculate this very, very carefully! Why?
   //         1. If a RZ end is at R=0, you can't miss!
   //         2. If you just fall off in RZ, the answer must
   //            be consistent with adjacent G4PolyhedraSide faces.
   // (2) implies that only variables used by other G4PolyhedraSide
   // faces may be used, which includes only: p, v, and the edge corners.
   // It also means that one side is a ">" or "<", which the other
   // must be ">=" or "<=". Fortunately, this isn't a new problem.
   // The solution below I borrowed from Joseph O'Rourke,
   // "Computational Geometry in C (Second Edition)"
   // See: http://cs.smith.edu/~orourke/
   //
   G4ThreeVector ic = p + distance*v - vec.center;
   G4double atRZ = vec.surfRZ.dot(ic);
   
   if (atRZ < 0)
   {
     if (r[0]==0) return true;    // Can't miss!
     
     if (atRZ < -lenRZ*1.2) return false;  // Forget it! Missed by a mile.
     
     G4ThreeVector q = p + v;    
     G4ThreeVector qa = q - vec.edges[0]->corner[0],
                   qb = q - vec.edges[1]->corner[0];
     G4ThreeVector qacb = qa.cross(qb);
     if (normSign*qacb.dot(v) < 0) return false;
     
     if (distFromSurface < 0)
     {
       if (atRZ < -lenRZ-surfTolerance) return false;
     }
   }
   else if (atRZ > 0)
   {
     if (r[1]==0) return true;    // Can't miss!
     
     if (atRZ > lenRZ*1.2) return false;  // Missed by a mile
     
     G4ThreeVector q = p + v;    
     G4ThreeVector qa = q - vec.edges[0]->corner[1],
                   qb = q - vec.edges[1]->corner[1];
     G4ThreeVector qacb = qa.cross(qb);
     if (normSign*qacb.dot(v) >= 0) return false;
     
     if (distFromSurface < 0)
     {
       if (atRZ > lenRZ+surfTolerance) return false;
     }
   }
 
   return true;
 }

G4int G4PolyhedraSide::LineHitsSegments	(	const G4ThreeVector &	p,
		const G4ThreeVector &	v,
		G4int *	i1,
		G4int *	i2
	)

protected

Definition at line 897 of file G4PolyhedraSide.cc.

References cone, G4IntersectingCone::LineHitsCone(), n, PhiSegment(), CLHEP::Hep3Vector::x(), and CLHEP::Hep3Vector::y().

 {
   G4double s1, s2;
   //
   // First, decide if and where the line intersects the cone
   //
   G4int n = cone->LineHitsCone( p, v, &s1, &s2 );
   
   if (n==0) return 0;
   
   //
   // Try first intersection.
   //
   *i1 = PhiSegment( std::atan2( p.y() + s1*v.y(), p.x() + s1*v.x() ) );
   if (n==1)
   {
     return (*i1 < 0) ? 0 : 1;
   }
   
   //
   // Try second intersection
   //
   *i2 = PhiSegment( std::atan2( p.y() + s2*v.y(), p.x() + s2*v.x() ) );
   if (*i1 == *i2) return 0;
   
   if (*i1 < 0)
   {
     if (*i2 < 0) return 0;
     *i1 = *i2;
     return 1;
   }
 
   if (*i2 < 0) return 1;
   
   return 2;
 }

G4ThreeVector G4PolyhedraSide::Normal	(	const G4ThreeVector &	p,
		G4double *	bestDistance
	)

virtual

Implements G4VCSGface.

Definition at line 650 of file G4PolyhedraSide.cc.

References ClosestPhiSegment(), DistanceToOneSide(), GetPhi(), G4PolyhedraSide::sG4PolyhedraSideVec::normal, and vecs.

 {
   //
   // Which phi segment is closest to this point?
   //
   G4int iPhi = ClosestPhiSegment( GetPhi(p) );
 
   //
   // Get distance to this segment
   //
   G4double norm;
   *bestDistance = DistanceToOneSide( p, vecs[iPhi], &norm );
 
   return vecs[iPhi].normal;
 }

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

Definition at line 354 of file G4PolyhedraSide.cc.

References cone, CopyStuff(), edges, and vecs.

 {
   if (this == &source) return *this;
   
   delete cone;
   delete [] vecs;
   delete [] edges;
   
   CopyStuff( source );
 
   return *this;
 }

G4int G4PolyhedraSide::PhiSegment ( G4double phi )

protected

Definition at line 971 of file G4PolyhedraSide.cc.

References deltaPhi, numSide, phiIsOpen, startPhi, and python.hepunit::twopi.

Referenced by ClosestPhiSegment(), Extent(), and LineHitsSegments().

 {
   //
   // How far are we from phiStart? Come up with a positive answer
   // that is less than 2*PI
   //
   G4double phi = phi0 - startPhi;
   while( phi < 0      ) phi += twopi;
   while( phi > twopi ) phi -= twopi;
 
   //
   // Divide
   //
   G4int answer = (G4int)(phi/deltaPhi);
   
   if (answer >= numSide)
   {
     if (phiIsOpen)
     {
       return -1;  // Looks like we missed
     }
     else
     {
       answer = numSide-1;  // Probably just roundoff
     }
   }
   
   return answer;
 }

G4double G4PolyhedraSide::SurfaceArea ( )

virtual

Implements G4VCSGface.

Definition at line 1251 of file G4PolyhedraSide.cc.

References G4PolyhedraSide::sG4PolyhedraSideEdge::corner, G4PolyhedraSide::sG4PolyhedraSideVec::edges, GetPointOnPlane(), numSide, and vecs.

 {
   if( fSurfaceArea==0. )
   { 
     // Define the variables
     //
     G4double area,areas;
     G4ThreeVector point1;
     G4ThreeVector v1,v2,v3,v4; 
     G4PolyhedraSideVec *vec = vecs;
     areas=0.;
 
     // Do a loop on all SideEdge
     //
     do
     {
       // Define 4points for a Plane or Triangle
       //
       v1=vec->edges[0]->corner[0];
       v2=vec->edges[0]->corner[1];
       v3=vec->edges[1]->corner[1];
       v4=vec->edges[1]->corner[0];
       point1=GetPointOnPlane(v1,v2,v3,v4,&area);
       areas+=area;
     } while( ++vec < vecs + numSide);
 
     fSurfaceArea=areas;
   }
   return fSurfaceArea;
 }

G4double G4PolyhedraSide::SurfaceTriangle	(	G4ThreeVector	p1,
		G4ThreeVector	p2,
		G4ThreeVector	p3,
		G4ThreeVector *	p4
	)

Definition at line 1206 of file G4PolyhedraSide.cc.

References G4UniformRand, and test::v.

Referenced by GetPointOnPlane().

 {
   G4ThreeVector v, w;
   
   v = p3 - p1;
   w = p1 - p2;
   G4double lambda1 = G4UniformRand();
   G4double lambda2 = lambda1*G4UniformRand();
  
   *p4=p2 + lambda1*w + lambda2*v;
   return 0.5*(v.cross(w)).mag();
 }