#include <G4GenericTrap.hh>

Inheritance diagram for G4GenericTrap:

Public Member Functions
	G4GenericTrap (const G4String &name, G4double halfZ, const std::vector< G4TwoVector > &vertices)

	~G4GenericTrap ()

G4double	GetZHalfLength () const

G4int	GetNofVertices () const

G4TwoVector	GetVertex (G4int index) const

const std::vector< G4TwoVector > &	GetVertices () const

G4double	GetTwistAngle (G4int index) const

G4bool	IsTwisted () const

G4int	GetVisSubdivisions () const

void	SetVisSubdivisions (G4int subdiv)

EInside	Inside (const G4ThreeVector &p) const

G4ThreeVector	SurfaceNormal (const G4ThreeVector &p) const

G4double	DistanceToIn (const G4ThreeVector &p, const G4ThreeVector &v) const

G4double	DistanceToIn (const G4ThreeVector &p) const

G4double	DistanceToOut (const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=false, G4bool validNorm=0, G4ThreeVector n=0) const

G4double	DistanceToOut (const G4ThreeVector &p) const

G4bool	CalculateExtent (const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pmin, G4double &pmax) const

G4GeometryType	GetEntityType () const

G4VSolid *	Clone () const

std::ostream &	StreamInfo (std::ostream &os) const

G4ThreeVector	GetPointOnSurface () const

G4double	GetCubicVolume ()

G4double	GetSurfaceArea ()

G4Polyhedron *	GetPolyhedron () const

void	DescribeYourselfTo (G4VGraphicsScene &scene) const

G4VisExtent	GetExtent () const

G4Polyhedron *	CreatePolyhedron () const

	G4GenericTrap (__void__ &)

	G4GenericTrap (const G4GenericTrap &rhs)

G4GenericTrap &	operator= (const G4GenericTrap &rhs)

Public Member Functions inherited from G4VSolid
	G4VSolid (const G4String &name)

virtual	~G4VSolid ()

G4bool	operator== (const G4VSolid &s) const

G4String	GetName () const

void	SetName (const G4String &name)

G4double	GetTolerance () const

virtual void	ComputeDimensions (G4VPVParameterisation p, const G4int n, const G4VPhysicalVolume pRep)

void	DumpInfo () const

virtual const G4VSolid *	GetConstituentSolid (G4int no) const

virtual G4VSolid *	GetConstituentSolid (G4int no)

virtual const G4DisplacedSolid *	GetDisplacedSolidPtr () const

virtual G4DisplacedSolid *	GetDisplacedSolidPtr ()

	G4VSolid (__void__ &)

	G4VSolid (const G4VSolid &rhs)

G4VSolid &	operator= (const G4VSolid &rhs)

Protected Attributes
G4Polyhedron *	fpPolyhedron

Protected Attributes inherited from G4VSolid
G4double	kCarTolerance

Additional Inherited Members
Protected Member Functions inherited from G4VSolid
void	CalculateClippedPolygonExtent (G4ThreeVectorList &pPolygon, const G4VoxelLimits &pVoxelLimit, const EAxis pAxis, G4double &pMin, G4double &pMax) const

void	ClipCrossSection (G4ThreeVectorList *pVertices, const G4int pSectionIndex, const G4VoxelLimits &pVoxelLimit, const EAxis pAxis, G4double &pMin, G4double &pMax) const

void	ClipBetweenSections (G4ThreeVectorList *pVertices, const G4int pSectionIndex, const G4VoxelLimits &pVoxelLimit, const EAxis pAxis, G4double &pMin, G4double &pMax) const

void	ClipPolygon (G4ThreeVectorList &pPolygon, const G4VoxelLimits &pVoxelLimit, const EAxis pAxis) const

G4double	EstimateCubicVolume (G4int nStat, G4double epsilon) const

G4double	EstimateSurfaceArea (G4int nStat, G4double ell) const

Detailed Description

Definition at line 79 of file G4GenericTrap.hh.

Constructor & Destructor Documentation

G4GenericTrap::G4GenericTrap	(	const G4String &	name,
		G4double	halfZ,
		const std::vector< G4TwoVector > &	vertices
	)

Definition at line 66 of file G4GenericTrap.cc.

References FatalErrorInArgument, G4endl, G4Exception(), JustWarning, and G4VSolid::kCarTolerance.

Referenced by Clone().

   : G4VSolid(name),
     fpPolyhedron(0),
     fDz(halfZ),
     fVertices(),
     fIsTwisted(false),
     fTessellatedSolid(0),
     fMinBBoxVector(G4ThreeVector(0,0,0)),
     fMaxBBoxVector(G4ThreeVector(0,0,0)),
     fVisSubdivisions(0),
     fSurfaceArea(0.),
     fCubicVolume(0.)
    
 {
   // General constructor
   const G4double min_length=5*1.e-6;
   G4double length = 0.;
   G4int k=0;
   G4String errorDescription = "InvalidSetup in \" ";
   errorDescription += name;
   errorDescription += "\""; 
 
   halfCarTolerance = kCarTolerance*0.5;
 
   // Check vertices size
 
   if ( G4int(vertices.size()) != fgkNofVertices )
   {
     G4Exception("G4GenericTrap::G4GenericTrap()", "GeomSolids0002",
                 FatalErrorInArgument, "Number of vertices != 8");
   }            
   
   // Check dZ
   // 
   if (halfZ < kCarTolerance)
   {
      G4Exception("G4GenericTrap::G4GenericTrap()", "GeomSolids0002",
                 FatalErrorInArgument, "dZ is too small or negative");
   }           
  
   // Check Ordering and Copy vertices 
   //
   if(CheckOrder(vertices))
   {
     for (G4int i=0; i<fgkNofVertices; ++i) {fVertices.push_back(vertices[i]);}
   }
   else
   { 
     for (G4int i=0; i <4; ++i) {fVertices.push_back(vertices[3-i]);}
     for (G4int i=0; i <4; ++i) {fVertices.push_back(vertices[7-i]);}
   }
 
    // Check length of segments and Adjust
   // 
   for (G4int j=0; j < 2; j++)
   {
     for (G4int i=1; i<4; ++i)
     {
       k = j*4+i;
       length = (fVertices[k]-fVertices[k-1]).mag();
       if ( ( length < min_length) && ( length > kCarTolerance ) )
       {
         std::ostringstream message;
         message << "Length segment is too small." << G4endl
                 << "Distance between " << fVertices[k-1] << " and "
                 << fVertices[k] << " is only " << length << " mm !"; 
         G4Exception("G4GenericTrap::G4GenericTrap()", "GeomSolids1001",
                     JustWarning, message, "Vertices will be collapsed.");
         fVertices[k]=fVertices[k-1];
       }
     }
   }
 
   // Compute Twist
   //
   for( G4int i=0; i<4; i++) { fTwist[i]=0.; }
   fIsTwisted = ComputeIsTwisted();
 
   // Compute Bounding Box 
   //
   ComputeBBox();
   
   // If not twisted - create tessellated solid 
   // (an alternative implementation for testing)
   //
 #ifdef G4TESS_TEST
    if ( !fIsTwisted )  { fTessellatedSolid = CreateTessellatedSolid(); }
 #endif
 }

G4GenericTrap::~G4GenericTrap ( )

Definition at line 181 of file G4GenericTrap.cc.

 {
   // Destructor
   delete fTessellatedSolid;
 }

G4GenericTrap::G4GenericTrap ( __void__ & a )

Definition at line 159 of file G4GenericTrap.cc.

   : G4VSolid(a),
     fpPolyhedron(0),
     halfCarTolerance(0.),
     fDz(0.),
     fVertices(),
     fIsTwisted(false),
     fTessellatedSolid(0),
     fMinBBoxVector(G4ThreeVector(0,0,0)),
     fMaxBBoxVector(G4ThreeVector(0,0,0)),
     fVisSubdivisions(0),
     fSurfaceArea(0.),
     fCubicVolume(0.)
 {
   // Fake default constructor - sets only member data and allocates memory
   //                            for usage restricted to object persistency.
 
   for (size_t i=0; i<4; ++i)  { fTwist[i]=0.; }
 }

G4GenericTrap::G4GenericTrap ( const G4GenericTrap & rhs )

Definition at line 189 of file G4GenericTrap.cc.

   : G4VSolid(rhs),
     fpPolyhedron(0), halfCarTolerance(rhs.halfCarTolerance),
     fDz(rhs.fDz), fVertices(rhs.fVertices),
     fIsTwisted(rhs.fIsTwisted), fTessellatedSolid(0),
     fMinBBoxVector(rhs.fMinBBoxVector), fMaxBBoxVector(rhs.fMaxBBoxVector),
     fVisSubdivisions(rhs.fVisSubdivisions),
     fSurfaceArea(rhs.fSurfaceArea), fCubicVolume(rhs.fCubicVolume) 
 {
    for (size_t i=0; i<4; ++i)  { fTwist[i] = rhs.fTwist[i]; }
 #ifdef G4TESS_TEST
    if (rhs.fTessellatedSolid && !fIsTwisted )
    { fTessellatedSolid = CreateTessellatedSolid(); } 
 #endif
 }

Member Function Documentation

G4bool G4GenericTrap::CalculateExtent	(	const EAxis	pAxis,
		const G4VoxelLimits &	pVoxelLimit,
		const G4AffineTransform &	pTransform,
		G4double &	pmin,
		G4double &	pmax
	)		const

virtual

Implements G4VSolid.

Definition at line 1193 of file G4GenericTrap.cc.

 {
 #ifdef G4TESS_TEST
   if ( fTessellatedSolid )
   {
     return fTessellatedSolid->CalculateExtent(pAxis, pVoxelLimit,
                                               pTransform, pMin, pMax);
   }     
 #endif 
 
   // Computes bounding vectors for a shape
   //
   G4double Dx,Dy;
   G4ThreeVector minVec = GetMinimumBBox();
   G4ThreeVector maxVec = GetMaximumBBox();
   Dx = 0.5*(maxVec.x()- minVec.x());
   Dy = 0.5*(maxVec.y()- minVec.y());
 
   if (!pTransform.IsRotated())
   {
     // Special case handling for unrotated shapes
     // Compute x/y/z mins and maxs respecting limits, with early returns
     // if outside limits. Then switch() on pAxis
     //
     G4double xoffset,xMin,xMax;
     G4double yoffset,yMin,yMax;
     G4double zoffset,zMin,zMax;
 
     xoffset=pTransform.NetTranslation().x();
     xMin=xoffset-Dx;
     xMax=xoffset+Dx;
     if (pVoxelLimit.IsXLimited())
     {
       if ( (xMin>pVoxelLimit.GetMaxXExtent()+kCarTolerance)
         || (xMax<pVoxelLimit.GetMinXExtent()-kCarTolerance) )
       {
         return false;
       }
       else
       {
         if (xMin<pVoxelLimit.GetMinXExtent())
         {
           xMin=pVoxelLimit.GetMinXExtent();
         }
         if (xMax>pVoxelLimit.GetMaxXExtent())
         {
           xMax=pVoxelLimit.GetMaxXExtent();
         }
       }
     }
 
     yoffset=pTransform.NetTranslation().y();
     yMin=yoffset-Dy;
     yMax=yoffset+Dy;
     if (pVoxelLimit.IsYLimited())
     {
       if ( (yMin>pVoxelLimit.GetMaxYExtent()+kCarTolerance)
         || (yMax<pVoxelLimit.GetMinYExtent()-kCarTolerance) )
       {
         return false;
       }
       else
       {
         if (yMin<pVoxelLimit.GetMinYExtent())
         {
           yMin=pVoxelLimit.GetMinYExtent();
         }
         if (yMax>pVoxelLimit.GetMaxYExtent())
         {
           yMax=pVoxelLimit.GetMaxYExtent();
         }
       }
     }
 
     zoffset=pTransform.NetTranslation().z();
     zMin=zoffset-fDz;
     zMax=zoffset+fDz;
     if (pVoxelLimit.IsZLimited())
     {
       if ( (zMin>pVoxelLimit.GetMaxZExtent()+kCarTolerance)
         || (zMax<pVoxelLimit.GetMinZExtent()-kCarTolerance) )
       {
         return false;
       }
       else
       {
         if (zMin<pVoxelLimit.GetMinZExtent())
         {
           zMin=pVoxelLimit.GetMinZExtent();
         }
         if (zMax>pVoxelLimit.GetMaxZExtent())
         {
           zMax=pVoxelLimit.GetMaxZExtent();
         }
       }
     }
 
     switch (pAxis)
     {
       case kXAxis:
         pMin = xMin;
         pMax = xMax;
         break;
       case kYAxis:
         pMin = yMin;
         pMax = yMax;
         break;
       case kZAxis:
         pMin = zMin;
         pMax = zMax;
         break;
       default:
         break;
     }
     pMin-=kCarTolerance;
     pMax+=kCarTolerance;
 
     return true;
   }
   else
   {
     // General rotated case - create and clip mesh to boundaries
 
     G4bool existsAfterClip=false;
     G4ThreeVectorList *vertices;
 
     pMin=+kInfinity;
     pMax=-kInfinity;
 
     // Calculate rotated vertex coordinates
     //
     vertices=CreateRotatedVertices(pTransform);
     ClipCrossSection(vertices,0,pVoxelLimit,pAxis,pMin,pMax);
     ClipCrossSection(vertices,4,pVoxelLimit,pAxis,pMin,pMax);
     ClipBetweenSections(vertices,0,pVoxelLimit,pAxis,pMin,pMax);
 
     if ( (pMin!=kInfinity) || (pMax!=-kInfinity) )
     {
       existsAfterClip=true;
     
       // Add 2*tolerance to avoid precision troubles
       //
       pMin-=kCarTolerance;
       pMax+=kCarTolerance;
     }
     else
     {
       // Check for case where completely enveloping clipping volume.
       // If point inside then we are confident that the solid completely
       // envelopes the clipping volume. Hence set min/max extents according
       // to clipping volume extents along the specified axis.
       //
       G4ThreeVector clipCentre(
               (pVoxelLimit.GetMinXExtent()+pVoxelLimit.GetMaxXExtent())*0.5,
               (pVoxelLimit.GetMinYExtent()+pVoxelLimit.GetMaxYExtent())*0.5,
               (pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5);
 
       if (Inside(pTransform.Inverse().TransformPoint(clipCentre))!=kOutside)
       {
         existsAfterClip=true;
         pMin=pVoxelLimit.GetMinExtent(pAxis);
         pMax=pVoxelLimit.GetMaxExtent(pAxis);
       }
     }
     delete vertices;
     return existsAfterClip;
   }
 }

G4VSolid * G4GenericTrap::Clone ( ) const

virtual

Reimplemented from G4VSolid.

Definition at line 1420 of file G4GenericTrap.cc.

References G4GenericTrap().

 {
   return new G4GenericTrap(*this);
 }

G4Polyhedron * G4GenericTrap::CreatePolyhedron ( ) const

virtual

Reimplemented from G4VSolid.

Definition at line 2147 of file G4GenericTrap.cc.

References G4PolyhedronArbitrary::AddFacet(), G4PolyhedronArbitrary::AddVertex(), G4TessellatedSolid::CreatePolyhedron(), GetTwistAngle(), GetVisSubdivisions(), G4PolyhedronArbitrary::InvertFacets(), G4PolyhedronArbitrary::SetReferences(), test::x, CLHEP::Hep2Vector::x(), CLHEP::Hep3Vector::x(), CLHEP::Hep2Vector::y(), and CLHEP::Hep3Vector::y().

Referenced by GetPolyhedron().

 {
 
 #ifdef G4TESS_TEST
   if ( fTessellatedSolid )
   { 
     return fTessellatedSolid->CreatePolyhedron();
   }  
 #endif 
   
   // Approximation of Twisted Side
   // Construct extra Points, if Twisted Side
   //
   G4PolyhedronArbitrary* polyhedron;
   size_t nVertices, nFacets;
 
   G4int subdivisions=0;
   G4int i;
   if(fIsTwisted)
   {
     if ( GetVisSubdivisions()!= 0 )
     {
       subdivisions=GetVisSubdivisions();
     }
     else
     {
       // Estimation of Number of Subdivisions for smooth visualisation
       //
       G4double maxTwist=0.;
       for(i=0; i<4; i++)
       {
         if(GetTwistAngle(i)>maxTwist) { maxTwist=GetTwistAngle(i); }
       }
 
       // Computes bounding vectors for the shape
       //
       G4double Dx,Dy;
       G4ThreeVector minVec = GetMinimumBBox();
       G4ThreeVector maxVec = GetMaximumBBox();
       Dx = 0.5*(maxVec.x()- minVec.y());
       Dy = 0.5*(maxVec.y()- minVec.y());
       if (Dy > Dx)  { Dx=Dy; }
     
       subdivisions=8*G4int(maxTwist/(Dx*Dx*Dx)*fDz);
       if (subdivisions<4)  { subdivisions=4; }
       if (subdivisions>30) { subdivisions=30; }
     }
   }
   G4int sub4=4*subdivisions;
   nVertices = 8+subdivisions*4;
   nFacets = 6+subdivisions*4;
   G4double cf=1./(subdivisions+1);
   polyhedron = new G4PolyhedronArbitrary (nVertices, nFacets);
 
   // Add Vertex
   //
   for (i=0;i<4;i++)
   {
     polyhedron->AddVertex(G4ThreeVector(fVertices[i].x(),
                                         fVertices[i].y(),-fDz));
   }
   for( i=0;i<subdivisions;i++)
   {
     for(G4int j=0;j<4;j++)
     {
       G4TwoVector u=fVertices[j]+cf*(i+1)*( fVertices[j+4]-fVertices[j]);
       polyhedron->AddVertex(G4ThreeVector(u.x(),u.y(),-fDz+cf*2*fDz*(i+1)));
     }    
   }
   for (i=4;i<8;i++)
   {
     polyhedron->AddVertex(G4ThreeVector(fVertices[i].x(),
                                         fVertices[i].y(),fDz));
   }
 
   // Add Facets
   //
   polyhedron->AddFacet(1,4,3,2);  //Z-plane
   for (i=0;i<subdivisions+1;i++)
   {
     G4int is=i*4;
     polyhedron->AddFacet(5+is,8+is,4+is,1+is);
     polyhedron->AddFacet(8+is,7+is,3+is,4+is);
     polyhedron->AddFacet(7+is,6+is,2+is,3+is);
     polyhedron->AddFacet(6+is,5+is,1+is,2+is); 
   }
   polyhedron->AddFacet(5+sub4,6+sub4,7+sub4,8+sub4);  //Z-plane
 
   polyhedron->SetReferences();
   polyhedron->InvertFacets();
 
   return (G4Polyhedron*) polyhedron;
 }

void G4GenericTrap::DescribeYourselfTo ( G4VGraphicsScene & scene ) const

virtual

Implements G4VSolid.

Definition at line 2110 of file G4GenericTrap.cc.

References G4VGraphicsScene::AddSolid(), and G4TessellatedSolid::DescribeYourselfTo().

 {
 
 #ifdef G4TESS_TEST
   if ( fTessellatedSolid )
   { 
     return fTessellatedSolid->DescribeYourselfTo(scene);
   }
 #endif  
   
   scene.AddSolid(*this);
 }

G4double G4GenericTrap::DistanceToIn	(	const G4ThreeVector &	p,
		const G4ThreeVector &	v
	)		const

virtual

Implements G4VSolid.

Definition at line 749 of file G4GenericTrap.cc.

References G4TessellatedSolid::DistanceToIn(), CLHEP::Hep3Vector::dot(), Inside(), kOutside, n, test::v, and CLHEP::Hep3Vector::z().

 {
 #ifdef G4TESS_TEST
   if ( fTessellatedSolid )
   { 
     return fTessellatedSolid->DistanceToIn(p, v);
   }  
 #endif
     
   G4double dist[5];
   G4ThreeVector n;
 
   // Check lateral faces
   //
   G4int i;
   for (i=0; i<4; i++)
   {
     dist[i]=DistToPlane(p, v, i);  
   }
 
   // Check Z planes
   //
   dist[4]=kInfinity;
   if (std::fabs(p.z())>fDz-halfCarTolerance)
   {
     if (v.z())
     {
       G4ThreeVector pt;
       if (p.z()>0)
       {
         dist[4] = (fDz-p.z())/v.z();
       }
       else
       {   
         dist[4] = (-fDz-p.z())/v.z();
       }   
       if (dist[4]<-halfCarTolerance)
       {
         dist[4]=kInfinity;
       }
       else
       {
         if(dist[4]<halfCarTolerance)
         {
           if(p.z()>0)  { n=G4ThreeVector(0,0,1); }
           else         { n=G4ThreeVector(0,0,-1); }
           if (n.dot(v)<0) { dist[4]=0.; }
           else            { dist[4]=kInfinity; }
         }
         pt=p+dist[4]*v;
         if (Inside(pt)==kOutside)  { dist[4]=kInfinity; }
       }
     }
   }   
   G4double distmin = dist[0];
   for (i=1;i<5;i++)
   {
     if (dist[i] < distmin)  { distmin = dist[i]; }
   }
 
   if (distmin<halfCarTolerance)  { distmin=0.; }
 
   return distmin;
 }    

G4double G4GenericTrap::DistanceToIn ( const G4ThreeVector & p ) const

virtual

Implements G4VSolid.

Definition at line 817 of file G4GenericTrap.cc.

References G4TessellatedSolid::DistanceToIn(), and CLHEP::Hep3Vector::z().

 {
   // Computes the closest distance from given point to this shape
 
 #ifdef G4TESS_TEST
   if ( fTessellatedSolid )
   {
     return fTessellatedSolid->DistanceToIn(p);
   }  
 #endif
  
   G4double safz = std::fabs(p.z())-fDz;
   if(safz<0) { safz=0; }
 
   G4int iseg;
   G4double safe  = safz;
   G4double safxy = safz;
  
   for (iseg=0; iseg<4; iseg++)
   { 
     safxy = SafetyToFace(p,iseg);
     if (safxy>safe)  { safe=safxy; }
   }
 
   return safe;
 }

G4double G4GenericTrap::DistanceToOut	(	const G4ThreeVector &	p,
		const G4ThreeVector &	v,
		const G4bool	calcNorm = `false`,
		G4bool *	validNorm = `0`,
		G4ThreeVector *	n = `0`
	)		const

virtual

Implements G4VSolid.

Definition at line 921 of file G4GenericTrap.cc.

References test::b, test::c, G4TessellatedSolid::DistanceToOut(), CLHEP::Hep3Vector::dot(), G4VSolid::DumpInfo(), G4endl, G4Exception(), JustWarning, G4VSolid::kCarTolerance, kOutside, python.hepunit::mm, test::v, CLHEP::Hep3Vector::x(), CLHEP::Hep3Vector::y(), and CLHEP::Hep3Vector::z().

 {
 #ifdef G4TESS_TEST
   if ( fTessellatedSolid )
   { 
     return fTessellatedSolid->DistanceToOut(p, v, calcNorm, validNorm, n);
   }
 #endif
 
   G4double distmin;
   G4bool lateral_cross = false;
   ESide side = kUndefined;
  
   if (calcNorm)  { *validNorm=true; } // All normals are valid
 
   if (v.z() < 0)
   {
     distmin=(-fDz-p.z())/v.z();
     if (calcNorm) { side=kMZ; *n=G4ThreeVector(0,0,-1); }
   }
   else
   {
     if (v.z() > 0)
     { 
       distmin = (fDz-p.z())/v.z(); 
       if (calcNorm) { side=kPZ; *n=G4ThreeVector(0,0,1); } 
     }
     else            { distmin = kInfinity; }
   }      
 
   G4double dz2 =0.5/fDz;
   G4double xa,xb,xc,xd;
   G4double ya,yb,yc,yd;
 
   for (G4int ipl=0; ipl<4; ipl++)
   {
     G4int j = (ipl+1)%4;
     xa=fVertices[ipl].x();
     ya=fVertices[ipl].y();
     xb=fVertices[ipl+4].x();
     yb=fVertices[ipl+4].y();
     xc=fVertices[j].x();
     yc=fVertices[j].y();
     xd=fVertices[4+j].x();
     yd=fVertices[4+j].y();
 
     if ( ((std::fabs(xb-xd)+std::fabs(yb-yd))<halfCarTolerance)
       || ((std::fabs(xa-xc)+std::fabs(ya-yc))<halfCarTolerance) )
     {
       G4double q=DistToTriangle(p,v,ipl) ;
       if ( (q>=0) && (q<distmin) )
       {
         distmin=q;
         lateral_cross=true;
         side=ESide(ipl+1);
       }
       continue;
     }
     G4double tx1 =dz2*(xb-xa);
     G4double ty1 =dz2*(yb-ya);
     G4double tx2 =dz2*(xd-xc);
     G4double ty2 =dz2*(yd-yc);
     G4double dzp =fDz+p.z();
     G4double xs1 =xa+tx1*dzp;
     G4double ys1 =ya+ty1*dzp;
     G4double xs2 =xc+tx2*dzp;
     G4double ys2 =yc+ty2*dzp;
     G4double dxs =xs2-xs1;
     G4double dys =ys2-ys1;
     G4double dtx =tx2-tx1;
     G4double dty =ty2-ty1;
     G4double a = (dtx*v.y()-dty*v.x()+(tx1*ty2-tx2*ty1)*v.z())*v.z();
     G4double b = dxs*v.y()-dys*v.x()+(dtx*p.y()-dty*p.x()+ty2*xs1-ty1*xs2
                + tx1*ys2-tx2*ys1)*v.z();
     G4double c=dxs*p.y()-dys*p.x()+xs1*ys2-xs2*ys1;
     G4double q=kInfinity;
 
     if (std::fabs(a) < kCarTolerance)
     {           
       if (std::fabs(b) < kCarTolerance) { continue; }
       q=-c/b;
          
       // Check for Point on the Surface
       //
       if ((q > -halfCarTolerance) && (q < distmin))
       {
         if (q < halfCarTolerance)
         {
           if (NormalToPlane(p,ipl).dot(v)<0.)  { continue; }
         }
         distmin =q;
         lateral_cross=true;
         side=ESide(ipl+1);
       }   
       continue;
     }
     G4double d=b*b-4*a*c;
     if (d >= 0.)
     {
       if (a > 0)  { q=0.5*(-b-std::sqrt(d))/a; }
       else        { q=0.5*(-b+std::sqrt(d))/a; }
 
       // Check for Point on the Surface
       //
       if (q > -halfCarTolerance )
       {
         if (q < distmin)
         {
           if(q < halfCarTolerance)
           {
             if (NormalToPlane(p,ipl).dot(v)<0.)  // Check second root
             {
               if (a > 0)  { q=0.5*(-b+std::sqrt(d))/a; }
               else        { q=0.5*(-b-std::sqrt(d))/a; }
               if (( q > halfCarTolerance) && (q < distmin))
               {
                 distmin=q;
                 lateral_cross = true;
                 side=ESide(ipl+1);
               }
               continue;
             }
           }
           distmin = q;
           lateral_cross = true;
           side=ESide(ipl+1);
         }
       }
       else
       {
         if (a > 0)  { q=0.5*(-b+std::sqrt(d))/a; }
         else        { q=0.5*(-b-std::sqrt(d))/a; }
 
         // Check for Point on the Surface
         //
         if ((q > -halfCarTolerance) && (q < distmin))
         {
           if (q < halfCarTolerance)
           {
             if (NormalToPlane(p,ipl).dot(v)<0.)  // Check second root
             {
               if (a > 0)  { q=0.5*(-b-std::sqrt(d))/a; }
               else        { q=0.5*(-b+std::sqrt(d))/a; }
               if ( ( q > halfCarTolerance) && (q < distmin) )
               {
                 distmin=q;
                 lateral_cross = true;
                 side=ESide(ipl+1);
               }
               continue;
             }  
           }
           distmin =q;
           lateral_cross = true;
           side=ESide(ipl+1);
         }   
       }
     }
   }
   if (!lateral_cross)  // Make sure that track crosses the top or bottom
   {
     if (distmin >= kInfinity)  { distmin=kCarTolerance; }
     G4ThreeVector pt=p+distmin*v;
     
     // Check if propagated point is in the polygon
     //
     G4int i=0;
     if (v.z()>0.) { i=4; }
     std::vector<G4TwoVector> xy;
     for ( G4int j=0; j<4; j++)  { xy.push_back(fVertices[i+j]); }
 
     // Check Inside
     //
     if (InsidePolygone(pt,xy)==kOutside)
     { 
       if(calcNorm)
       {
         if (v.z()>0) {side= kPZ; *n = G4ThreeVector(0,0,1);}
         else         { side=kMZ; *n = G4ThreeVector(0,0,-1);}
       } 
       return 0.;
     }
     else
     {
       if(v.z()>0) {side=kPZ;}
       else        {side=kMZ;}
     }
   }
 
   if (calcNorm)
   {
     G4ThreeVector pt=p+v*distmin;     
     switch (side)
     {
       case kXY0:
         *n=NormalToPlane(pt,0);
         break;
       case kXY1:
         *n=NormalToPlane(pt,1);
         break;
       case kXY2:
         *n=NormalToPlane(pt,2);
         break;
       case kXY3:
         *n=NormalToPlane(pt,3);
         break;
       case kPZ:
         *n=G4ThreeVector(0,0,1);
         break;
       case kMZ:
         *n=G4ThreeVector(0,0,-1);
         break;
       default:
         DumpInfo();
         std::ostringstream message;
         G4int oldprc = message.precision(16);
         message << "Undefined side for valid surface normal to solid." << G4endl
                 << "Position:" << G4endl
                 << "  p.x() = "   << p.x()/mm << " mm" << G4endl
                 << "  p.y() = "   << p.y()/mm << " mm" << G4endl
                 << "  p.z() = "   << p.z()/mm << " mm" << G4endl
                 << "Direction:" << G4endl
                 << "  v.x() = "   << v.x() << G4endl
                 << "  v.y() = "   << v.y() << G4endl
                 << "  v.z() = "   << v.z() << G4endl
                 << "Proposed distance :" << G4endl
                 << "  distmin = " << distmin/mm << " mm";
         message.precision(oldprc);
         G4Exception("G4GenericTrap::DistanceToOut(p,v,..)",
                     "GeomSolids1002", JustWarning, message);
         break;
      }
   }
  
   if (distmin<halfCarTolerance)  { distmin=0.; }
  
   return distmin;
 }    

G4double G4GenericTrap::DistanceToOut ( const G4ThreeVector & p ) const

virtual

Implements G4VSolid.

Definition at line 1166 of file G4GenericTrap.cc.

References G4TessellatedSolid::DistanceToOut(), and CLHEP::Hep3Vector::z().

 {
 
 #ifdef G4TESS_TEST
   if ( fTessellatedSolid )
   { 
     return fTessellatedSolid->DistanceToOut(p);
   }  
 #endif
 
   G4double safz = fDz-std::fabs(p.z());
   if (safz<0) { safz = 0; }
 
   G4double safe  = safz;
   G4double safxy = safz;
  
   for (G4int iseg=0; iseg<4; iseg++)
   { 
     safxy = std::fabs(SafetyToFace(p,iseg));
     if (safxy < safe)  { safe = safxy; }
   }
 
   return safe;
 }    

G4double G4GenericTrap::GetCubicVolume ( )

virtual

Reimplemented from G4VSolid.

Definition at line 1538 of file G4GenericTrap.cc.

References G4VSolid::GetCubicVolume().

 {
   if(fCubicVolume != 0.) {;}
   else   { fCubicVolume = G4VSolid::GetCubicVolume(); }
   return fCubicVolume;
 }

G4GeometryType G4GenericTrap::GetEntityType ( ) const

virtual

Implements G4VSolid.

Definition at line 1413 of file G4GenericTrap.cc.

Referenced by StreamInfo().

 {
   return G4String("G4GenericTrap");
 }

G4VisExtent G4GenericTrap::GetExtent ( ) const

virtual

Reimplemented from G4VSolid.

Definition at line 2125 of file G4GenericTrap.cc.

References G4TessellatedSolid::GetExtent(), CLHEP::Hep3Vector::x(), and CLHEP::Hep3Vector::y().

 {
   // Computes bounding vectors for the shape
 
 #ifdef G4TESS_TEST
   if ( fTessellatedSolid )
   { 
     return fTessellatedSolid->GetExtent();
   } 
 #endif
    
   G4double Dx,Dy;
   G4ThreeVector minVec = GetMinimumBBox();
   G4ThreeVector maxVec = GetMaximumBBox();
   Dx = 0.5*(maxVec.x()- minVec.x());
   Dy = 0.5*(maxVec.y()- minVec.y());
 
   return G4VisExtent (-Dx, Dx, -Dy, Dy, -fDz, fDz); 
 }    

G4int G4GenericTrap::GetNofVertices ( ) const

inline

G4ThreeVector G4GenericTrap::GetPointOnSurface ( ) const

virtual

Reimplemented from G4VSolid.

Definition at line 1450 of file G4GenericTrap.cc.

References G4UniformRand, G4TessellatedSolid::GetPointOnSurface(), test::v, test::x, CLHEP::Hep2Vector::x(), and CLHEP::Hep2Vector::y().

 {
 
 #ifdef G4TESS_TEST
   if ( fTessellatedSolid )
   { 
     return fTessellatedSolid->GetPointOnSurface();
   }  
 #endif
 
   G4ThreeVector point;
   G4TwoVector u,v,w;
   G4double rand,area,chose,cf,lambda0,lambda1,alfa,beta,zp;
   G4int ipl,j;
  
   std::vector<G4ThreeVector> vertices;
   for (G4int i=0; i<4;i++)
   {
     vertices.push_back(G4ThreeVector(fVertices[i].x(),fVertices[i].y(),-fDz));
   }
   for (G4int i=4; i<8;i++)
   {
     vertices.push_back(G4ThreeVector(fVertices[i].x(),fVertices[i].y(),fDz));
   }
 
   // Surface Area of Planes(only estimation for twisted)
   //
   G4double Surface0=GetFaceSurfaceArea(vertices[0],vertices[1],
                                        vertices[2],vertices[3]);//-fDz plane 
   G4double Surface1=GetFaceSurfaceArea(vertices[0],vertices[1],
                                        vertices[5],vertices[4]);// Lat plane
   G4double Surface2=GetFaceSurfaceArea(vertices[3],vertices[0],
                                        vertices[4],vertices[7]);// Lat plane
   G4double Surface3=GetFaceSurfaceArea(vertices[2],vertices[3],
                                        vertices[7],vertices[6]);// Lat plane
   G4double Surface4=GetFaceSurfaceArea(vertices[2],vertices[1],
                                        vertices[5],vertices[6]);// Lat plane
   G4double Surface5=GetFaceSurfaceArea(vertices[4],vertices[5],
                                        vertices[6],vertices[7]);// fDz plane
   rand = G4UniformRand();
   area = Surface0+Surface1+Surface2+Surface3+Surface4+Surface5;
   chose = rand*area;
 
   if ( ( chose < Surface0)
     || ( chose > (Surface0+Surface1+Surface2+Surface3+Surface4)) )
   {                                        // fDz or -fDz Plane
     ipl = G4int(G4UniformRand()*4);
     j = (ipl+1)%4;
     if(chose < Surface0)
     { 
       zp = -fDz;
       u = fVertices[ipl]; v = fVertices[j];
       w = fVertices[(ipl+3)%4]; 
     }
     else
     {
       zp = fDz;
       u = fVertices[ipl+4]; v = fVertices[j+4];
       w = fVertices[(ipl+3)%4+4]; 
     }
     alfa = G4UniformRand();
     beta = G4UniformRand();
     lambda1=alfa*beta;
     lambda0=alfa-lambda1;
     v = v-u;
     w = w-u;
     v = u+lambda0*v+lambda1*w;
   }
   else                                     // Lateral Plane Twisted or Not
   {
     if (chose < Surface0+Surface1) { ipl=0; }
     else if (chose < Surface0+Surface1+Surface2) { ipl=1; }
     else if (chose < Surface0+Surface1+Surface2+Surface3) { ipl=2; }
     else { ipl=3; }
     j = (ipl+1)%4;
     zp = -fDz+G4UniformRand()*2*fDz;
     cf = 0.5*(fDz-zp)/fDz;
     u = fVertices[ipl+4]+cf*( fVertices[ipl]-fVertices[ipl+4]);
     v = fVertices[j+4]+cf*(fVertices[j]-fVertices[j+4]); 
     v = u+(v-u)*G4UniformRand();
   }
   point=G4ThreeVector(v.x(),v.y(),zp);
 
   return point;
 }

G4Polyhedron * G4GenericTrap::GetPolyhedron ( ) const

virtual

Reimplemented from G4VSolid.

Definition at line 2088 of file G4GenericTrap.cc.

References CreatePolyhedron(), fpPolyhedron, HepPolyhedron::GetNumberOfRotationSteps(), G4Polyhedron::GetNumberOfRotationStepsAtTimeOfCreation(), and G4TessellatedSolid::GetPolyhedron().

 {
 
 #ifdef G4TESS_TEST
   if ( fTessellatedSolid )
   { 
     return fTessellatedSolid->GetPolyhedron();
   }
 #endif  
   
   if ( (!fpPolyhedron)
     || (fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() !=
         fpPolyhedron->GetNumberOfRotationSteps()) )
   {
     delete fpPolyhedron;
     fpPolyhedron = CreatePolyhedron();
   }
   return fpPolyhedron;
 }    

G4double G4GenericTrap::GetSurfaceArea ( )

virtual

Reimplemented from G4VSolid.

Definition at line 1547 of file G4GenericTrap.cc.

References G4VSolid::GetSurfaceArea(), and test::x.

 {
   if(fSurfaceArea != 0.) {;}
   else
   {
     std::vector<G4ThreeVector> vertices;
     for (G4int i=0; i<4;i++)
     {
       vertices.push_back(G4ThreeVector(fVertices[i].x(),fVertices[i].y(),-fDz));
     }
     for (G4int i=4; i<8;i++)
     {
       vertices.push_back(G4ThreeVector(fVertices[i].x(),fVertices[i].y(),fDz));
     }
 
     // Surface Area of Planes(only estimation for twisted)
     //
     G4double fSurface0=GetFaceSurfaceArea(vertices[0],vertices[1],
                                           vertices[2],vertices[3]);//-fDz plane
     G4double fSurface1=GetFaceSurfaceArea(vertices[0],vertices[1],
                                           vertices[5],vertices[4]);// Lat plane
     G4double fSurface2=GetFaceSurfaceArea(vertices[3],vertices[0],
                                           vertices[4],vertices[7]);// Lat plane 
     G4double fSurface3=GetFaceSurfaceArea(vertices[2],vertices[3],
                                           vertices[7],vertices[6]);// Lat plane
     G4double fSurface4=GetFaceSurfaceArea(vertices[2],vertices[1],
                                           vertices[5],vertices[6]);// Lat plane
     G4double fSurface5=GetFaceSurfaceArea(vertices[4],vertices[5],
                                           vertices[6],vertices[7]);// fDz plane 
 
     // Total Surface Area
     //
     if(!fIsTwisted)
     {
       fSurfaceArea = fSurface0+fSurface1+fSurface2
                    + fSurface3+fSurface4+fSurface5;
     }
     else
     {
       fSurfaceArea = G4VSolid::GetSurfaceArea();
     }
   }
   return fSurfaceArea;
 }

G4double G4GenericTrap::GetTwistAngle ( G4int index ) const

inline

Referenced by CreatePolyhedron(), and SurfaceNormal().

G4TwoVector G4GenericTrap::GetVertex ( G4int index ) const

inline

const std::vector<G4TwoVector>& G4GenericTrap::GetVertices ( ) const

inline

Referenced by G4GDMLWriteSolids::GenTrapWrite().

G4int G4GenericTrap::GetVisSubdivisions ( ) const

inline

Referenced by CreatePolyhedron().

G4double G4GenericTrap::GetZHalfLength ( ) const

inline

Referenced by G4GDMLWriteSolids::GenTrapWrite().

EInside G4GenericTrap::Inside ( const G4ThreeVector & p ) const

virtual

Implements G4VSolid.

Definition at line 339 of file G4GenericTrap.cc.

References G4TessellatedSolid::Inside(), kInside, kOutside, kSurface, and CLHEP::Hep3Vector::z().

Referenced by CalculateExtent(), and DistanceToIn().

 {
   // Test if point is inside this shape
 
 #ifdef G4TESS_TEST
    if ( fTessellatedSolid )
    { 
      return fTessellatedSolid->Inside(p);
    }
 #endif  
 
   EInside innew=kOutside;
   std::vector<G4TwoVector> xy;
  
   if (std::fabs(p.z()) <= fDz+halfCarTolerance)  // First check Z range
   {
     // Compute intersection between Z plane containing point and the shape
     //
     G4double cf = 0.5*(fDz-p.z())/fDz;
     for (G4int i=0; i<4; i++)
     {
       xy.push_back(fVertices[i+4]+cf*( fVertices[i]-fVertices[i+4]));
     }
 
     innew=InsidePolygone(p,xy);
 
     if( (innew==kInside) || (innew==kSurface) )
     { 
       if(std::fabs(p.z()) > fDz-halfCarTolerance)  { innew=kSurface; }
     }
   }
   return innew;    
 } 

G4bool G4GenericTrap::IsTwisted ( ) const

inline

G4GenericTrap & G4GenericTrap::operator= ( const G4GenericTrap & rhs )

Definition at line 207 of file G4GenericTrap.cc.

References fpPolyhedron, and G4VSolid::operator=().

 {
    // Check assignment to self
    //
    if (this == &rhs)  { return *this; }
 
    // Copy base class data
    //
    G4VSolid::operator=(rhs);
 
    // Copy data
    //
    fpPolyhedron = 0; halfCarTolerance = rhs.halfCarTolerance;
    fDz = rhs.fDz; fVertices = rhs.fVertices;
    fIsTwisted = rhs.fIsTwisted; fTessellatedSolid = 0;
    fMinBBoxVector = rhs.fMinBBoxVector; fMaxBBoxVector = rhs.fMaxBBoxVector;
    fVisSubdivisions = rhs.fVisSubdivisions;
    fSurfaceArea = rhs.fSurfaceArea; fCubicVolume = rhs.fCubicVolume;
 
    for (size_t i=0; i<4; ++i)  { fTwist[i] = rhs.fTwist[i]; }
 #ifdef G4TESS_TEST
    if (rhs.fTessellatedSolid && !fIsTwisted )
    { delete fTessellatedSolid; fTessellatedSolid = CreateTessellatedSolid(); } 
 #endif
 
    return *this;
 }

void G4GenericTrap::SetVisSubdivisions ( G4int subdiv )

inline

std::ostream & G4GenericTrap::StreamInfo ( std::ostream & os ) const

virtual

Implements G4VSolid.

Definition at line 1427 of file G4GenericTrap.cc.

References G4endl, GetEntityType(), G4VSolid::GetName(), and python.hepunit::mm.

 {
   G4int oldprc = os.precision(16);
   os << "-----------------------------------------------------------\n"
      << "    *** Dump for solid - " << GetName() << " *** \n"
      << "    =================================================== \n"
      << " Solid geometry type: " << GetEntityType()  << G4endl
      << "   half length Z: " << fDz/mm << " mm \n"
      << "   list of vertices:\n";
      
   for ( G4int i=0; i<fgkNofVertices; ++i )
   {
     os << std::setw(5) << "#" << i 
        << "   vx = " << fVertices[i].x()/mm << " mm" 
        << "   vy = " << fVertices[i].y()/mm << " mm" << G4endl;
   }
   os.precision(oldprc);
 
   return os;
 } 

G4ThreeVector G4GenericTrap::SurfaceNormal ( const G4ThreeVector & p ) const

virtual

Implements G4VSolid.

Definition at line 375 of file G4GenericTrap.cc.

References CLHEP::Hep3Vector::cross(), G4Exception(), GetTwistAngle(), JustWarning, G4VSolid::kCarTolerance, G4TessellatedSolid::SurfaceNormal(), CLHEP::Hep3Vector::unit(), test::x, and CLHEP::Hep3Vector::z().

 {
   // Calculate side nearest to p, and return normal
   // If two sides are equidistant, sum of the Normal is returned
 
 #ifdef G4TESS_TEST
   if ( fTessellatedSolid )
   {
     return fTessellatedSolid->SurfaceNormal(p);
   }  
 #endif   
  
   G4ThreeVector lnorm, sumnorm(0.,0.,0.), apprnorm(0.,0.,1.),
                 p0, p1, p2, r1, r2, r3, r4;
   G4int noSurfaces = 0; 
   G4double distxy,distz;
   G4bool zPlusSide=false;
 
   distz = fDz-std::fabs(p.z());
   if (distz < halfCarTolerance)
   {
     if(p.z()>0)
     {
       zPlusSide=true;
       sumnorm=G4ThreeVector(0,0,1);
     }
     else
     {
       sumnorm=G4ThreeVector(0,0,-1);
     }
     noSurfaces ++;
   } 
 
   // Check lateral planes
   //
   std:: vector<G4TwoVector> vertices;  
   G4double cf = 0.5*(fDz-p.z())/fDz;
   for (G4int i=0; i<4; i++)
   {
     vertices.push_back(fVertices[i+4]+cf*(fVertices[i]-fVertices[i+4]));
   }
 
   // Compute distance for lateral planes
   //
   for (G4int q=0; q<4; q++)
   {
     p0=G4ThreeVector(vertices[q].x(),vertices[q].y(),p.z());
     if(zPlusSide)
     {
       p1=G4ThreeVector(fVertices[q].x(),fVertices[q].y(),-fDz);
     }
     else
     {
       p1=G4ThreeVector(fVertices[q+4].x(),fVertices[q+4].y(),fDz); 
     }
     p2=G4ThreeVector(vertices[(q+1)%4].x(),vertices[(q+1)%4].y(),p.z());
 
     // Collapsed vertices
     //
     if ( (p2-p0).mag2() < kCarTolerance )
     {
       if ( std::fabs(p.z()+fDz) > kCarTolerance )
       {
         p2=G4ThreeVector(fVertices[(q+1)%4].x(),fVertices[(q+1)%4].y(),-fDz);
       }
       else
       {
         p2=G4ThreeVector(fVertices[(q+1)%4+4].x(),fVertices[(q+1)%4+4].y(),fDz);
       }
     }
     lnorm = (p1-p0).cross(p2-p0);
     lnorm = lnorm.unit();
     if(zPlusSide)  { lnorm=-lnorm; }
 
     // Adjust Normal for Twisted Surface
     //
     if ( (fIsTwisted) && (GetTwistAngle(q)!=0) )
     {
       G4double normP=(p2-p0).mag();
       if(normP)
       {
         G4double proj=(p-p0).dot(p2-p0)/normP;
         if(proj<0)     { proj=0; }
         if(proj>normP) { proj=normP; }
         G4int j=(q+1)%4;
         r1=G4ThreeVector(fVertices[q+4].x(),fVertices[q+4].y(),fDz);
         r2=G4ThreeVector(fVertices[j+4].x(),fVertices[j+4].y(),fDz);
         r3=G4ThreeVector(fVertices[q].x(),fVertices[q].y(),-fDz);
         r4=G4ThreeVector(fVertices[j].x(),fVertices[j].y(),-fDz);
         r1=r1+proj*(r2-r1)/normP;
         r3=r3+proj*(r4-r3)/normP;
         r2=r1-r3;
         r4=r2.cross(p2-p0); r4=r4.unit();
         lnorm=r4;
       }
     }   // End if fIsTwisted
 
     distxy=std::fabs((p0-p).dot(lnorm));
     if ( distxy<halfCarTolerance )
     {
       noSurfaces ++;
 
       // Negative sign for Normal is taken for Outside Normal
       //
       sumnorm=sumnorm+lnorm;
     }
 
     // For ApproxSurfaceNormal
     //
     if (distxy<distz)
     {
       distz=distxy;
       apprnorm=lnorm;
     }      
   }  // End for loop
 
   // Calculate final Normal, add Normal in the Corners and Touching Sides
   //
   if ( noSurfaces == 0 )
   {
     G4Exception("G4GenericTrap::SurfaceNormal(p)", "GeomSolids1002",
                 JustWarning, "Point p is not on surface !?" );
     sumnorm=apprnorm;
     // Add Approximative Surface Normal Calculation?
   }
   else if ( noSurfaces == 1 )  { sumnorm = sumnorm; }
   else                         { sumnorm = sumnorm.unit(); }
 
   return sumnorm ; 
 }    

Field Documentation

G4Polyhedron* G4GenericTrap::fpPolyhedron

mutableprotected

Definition at line 191 of file G4GenericTrap.hh.

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

The documentation for this class was generated from the following files:

Public Member Functions

Protected Attributes

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Field Documentation