#include <G4Ellipsoid.hh>

Inheritance diagram for G4Ellipsoid:

Public Member Functions
	G4Ellipsoid (const G4String &pName, G4double pxSemiAxis, G4double pySemiAxis, G4double pzSemiAxis, G4double pzBottomCut=0, G4double pzTopCut=0)

virtual	~G4Ellipsoid ()

G4double	GetSemiAxisMax (G4int i) const

G4double	GetZBottomCut () const

G4double	GetZTopCut () const

void	SetSemiAxis (G4double x, G4double y, G4double z)

void	SetZCuts (G4double newzBottomCut, G4double newzTopCut)

G4double	GetCubicVolume ()

G4double	GetSurfaceArea ()

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

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

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=G4bool(false), G4bool validNorm=0, G4ThreeVector n=0) const

G4double	DistanceToOut (const G4ThreeVector &p) const

G4GeometryType	GetEntityType () const

G4VSolid *	Clone () const

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

G4ThreeVector	GetPointOnSurface () const

G4Polyhedron *	GetPolyhedron () const

void	DescribeYourselfTo (G4VGraphicsScene &scene) const

G4VisExtent	GetExtent () const

G4Polyhedron *	CreatePolyhedron () const

	G4Ellipsoid (__void__ &)

	G4Ellipsoid (const G4Ellipsoid &rhs)

G4Ellipsoid &	operator= (const G4Ellipsoid &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

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 Member Functions
G4ThreeVectorList *	CreateRotatedVertices (const G4AffineTransform &pT, G4int &noPV) const

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

Protected Attributes
G4Polyhedron *	fpPolyhedron

Protected Attributes inherited from G4VSolid
G4double	kCarTolerance

Detailed Description

Definition at line 59 of file G4Ellipsoid.hh.

Constructor & Destructor Documentation

G4Ellipsoid::G4Ellipsoid	(	const G4String &	pName,
		G4double	pxSemiAxis,
		G4double	pySemiAxis,
		G4double	pzSemiAxis,
		G4double	pzBottomCut = `0`,
		G4double	pzTopCut = `0`
	)

Definition at line 63 of file G4Ellipsoid.cc.

References FatalErrorInArgument, G4Exception(), G4GeometryTolerance::GetInstance(), G4VSolid::GetName(), G4GeometryTolerance::GetRadialTolerance(), G4VSolid::kCarTolerance, SetSemiAxis(), and SetZCuts().

Referenced by Clone().

   : G4VSolid(pName), fpPolyhedron(0), fCubicVolume(0.), fSurfaceArea(0.),
     zBottomCut(0.), zTopCut(0.)
 {
   // note: for users that want to use the full ellipsoid it is useful
   // to include a default for the cuts 
 
   kRadTolerance = G4GeometryTolerance::GetInstance()->GetRadialTolerance();
 
   halfCarTolerance = kCarTolerance*0.5;
   halfRadTolerance = kRadTolerance*0.5;
 
   // Check Semi-Axis
   if ( (pxSemiAxis<=0.) || (pySemiAxis<=0.) || (pzSemiAxis<=0.) )
   {
      std::ostringstream message;
      message << "Invalid semi-axis - " << GetName();
      G4Exception("G4Ellipsoid::G4Ellipsoid()", "GeomSolids0002",
                  FatalErrorInArgument, message);
   }
   SetSemiAxis(pxSemiAxis, pySemiAxis, pzSemiAxis);
 
   if ( pzBottomCut == 0 && pzTopCut == 0 )
   {
      SetZCuts(-pzSemiAxis, pzSemiAxis);
   }
   else if ( (pzBottomCut < pzSemiAxis) && (pzTopCut > -pzSemiAxis)
          && (pzBottomCut < pzTopCut) )
   {
      SetZCuts(pzBottomCut, pzTopCut);
   }
   else
   {
      std::ostringstream message;
      message << "Invalid z-coordinate for cutting plane - " << GetName();
      G4Exception("G4Ellipsoid::G4Ellipsoid()", "GeomSolids0002",
                  FatalErrorInArgument, message);
   }
 }

G4Ellipsoid::~G4Ellipsoid ( )

virtual

Definition at line 125 of file G4Ellipsoid.cc.

126 {

127 }

G4Ellipsoid::G4Ellipsoid ( __void__ & a )

Definition at line 113 of file G4Ellipsoid.cc.

   : G4VSolid(a), fpPolyhedron(0), kRadTolerance(0.),
     halfCarTolerance(0.), halfRadTolerance(0.), fCubicVolume(0.),
     fSurfaceArea(0.), xSemiAxis(0.), ySemiAxis(0.), zSemiAxis(0.),
     semiAxisMax(0.), zBottomCut(0.), zTopCut(0.)
 {
 }

G4Ellipsoid::G4Ellipsoid ( const G4Ellipsoid & rhs )

Definition at line 133 of file G4Ellipsoid.cc.

   : G4VSolid(rhs),
     fpPolyhedron(0), kRadTolerance(rhs.kRadTolerance),
     halfCarTolerance(rhs.halfCarTolerance),
     halfRadTolerance(rhs.halfRadTolerance),
     fCubicVolume(rhs.fCubicVolume), fSurfaceArea(rhs.fSurfaceArea),
     xSemiAxis(rhs.xSemiAxis), ySemiAxis(rhs.ySemiAxis),
     zSemiAxis(rhs.zSemiAxis), semiAxisMax(rhs.semiAxisMax),
     zBottomCut(rhs.zBottomCut), zTopCut(rhs.zTopCut)
 {
 }

Member Function Documentation

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

virtual

Implements G4VSolid.

Definition at line 189 of file G4Ellipsoid.cc.

 {
   if (!pTransform.IsRotated())
   {
     // Special case handling for unrotated solid ellipsoid
     // Compute x/y/z mins and maxs for bounding box respecting limits,
     // with early returns if outside limits. Then switch() on pAxis,
     // and compute exact x and y limit for x/y case
 
     G4double xoffset,xMin,xMax;
     G4double yoffset,yMin,yMax;
     G4double zoffset,zMin,zMax;
 
     G4double maxDiff,newMin,newMax;
     G4double xoff,yoff;
 
     xoffset=pTransform.NetTranslation().x();
     xMin=xoffset - xSemiAxis;
     xMax=xoffset + xSemiAxis;
     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 - ySemiAxis;
     yMax=yoffset + ySemiAxis;
     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 + (-zSemiAxis > zBottomCut ? -zSemiAxis : zBottomCut);
     zMax=zoffset + ( zSemiAxis < zTopCut ? zSemiAxis : zTopCut);
     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();
         }
       }
     }
 
     // if here, then known to cut bounding box around ellipsoid
     //
     xoff = (xoffset < xMin) ? (xMin-xoffset)
          : (xoffset > xMax) ? (xoffset-xMax) : 0.0;
     yoff = (yoffset < yMin) ? (yMin-yoffset)
          : (yoffset > yMax) ? (yoffset-yMax) : 0.0;
 
     // detailed calculations
     // NOTE: does not use X or Y offsets to adjust Z range,
     // and does not use Z offset to adjust X or Y range,
     // which is consistent with G4Sphere::CalculateExtent behavior
     //
     switch (pAxis)
     {
       case kXAxis:
         if (yoff==0.)
         {
           // YZ limits cross max/min x => no change
           //
           pMin=xMin;
           pMax=xMax;
         }
         else
         {
           // YZ limits don't cross max/min x => compute max delta x,
           // hence new mins/maxs
           //
           maxDiff= 1.0-sqr(yoff/ySemiAxis);
           if (maxDiff < 0.0) { return false; }
           maxDiff= xSemiAxis * std::sqrt(maxDiff);
           newMin=xoffset-maxDiff;
           newMax=xoffset+maxDiff;
           pMin=(newMin<xMin) ? xMin : newMin;
           pMax=(newMax>xMax) ? xMax : newMax;
         }
         break;
       case kYAxis:
         if (xoff==0.)
         {
           // XZ limits cross max/min y => no change
           //
           pMin=yMin;
           pMax=yMax;
         }
         else
         {
           // XZ limits don't cross max/min y => compute max delta y,
           // hence new mins/maxs
           //
           maxDiff= 1.0-sqr(xoff/xSemiAxis);
           if (maxDiff < 0.0) { return false; }
           maxDiff= ySemiAxis * std::sqrt(maxDiff);
           newMin=yoffset-maxDiff;
           newMax=yoffset+maxDiff;
           pMin=(newMin<yMin) ? yMin : newMin;
           pMax=(newMax>yMax) ? yMax : newMax;
         }
         break;
       case kZAxis:
         pMin=zMin;
         pMax=zMax;
         break;
       default:
         break;
     }
   
     pMin-=kCarTolerance;
     pMax+=kCarTolerance;
     return true;
   }
   else  // not rotated
   {
     G4int i,j,noEntries,noBetweenSections;
     G4bool existsAfterClip=false;
 
     // Calculate rotated vertex coordinates
 
     G4int noPolygonVertices=0;
     G4ThreeVectorList* vertices =
       CreateRotatedVertices(pTransform,noPolygonVertices);
 
     pMin=+kInfinity;
     pMax=-kInfinity;
 
     noEntries=vertices->size(); // noPolygonVertices*noPhiCrossSections
     noBetweenSections=noEntries-noPolygonVertices;
     
     G4ThreeVectorList ThetaPolygon;
     for (i=0;i<noEntries;i+=noPolygonVertices)
     {
       for(j=0;j<(noPolygonVertices/2)-1;j++)
       {
         ThetaPolygon.push_back((*vertices)[i+j]);  
         ThetaPolygon.push_back((*vertices)[i+j+1]);  
         ThetaPolygon.push_back((*vertices)[i+noPolygonVertices-2-j]);
         ThetaPolygon.push_back((*vertices)[i+noPolygonVertices-1-j]);
         CalculateClippedPolygonExtent(ThetaPolygon,pVoxelLimit,pAxis,pMin,pMax);
         ThetaPolygon.clear();
       }
     }
     for (i=0;i<noBetweenSections;i+=noPolygonVertices)
     {
       for(j=0;j<noPolygonVertices-1;j++)
       {
         ThetaPolygon.push_back((*vertices)[i+j]);  
         ThetaPolygon.push_back((*vertices)[i+j+1]);  
         ThetaPolygon.push_back((*vertices)[i+noPolygonVertices+j+1]);
         ThetaPolygon.push_back((*vertices)[i+noPolygonVertices+j]);
         CalculateClippedPolygonExtent(ThetaPolygon,pVoxelLimit,pAxis,pMin,pMax);
         ThetaPolygon.clear();
       }
       ThetaPolygon.push_back((*vertices)[i+noPolygonVertices-1]);
       ThetaPolygon.push_back((*vertices)[i]);
       ThetaPolygon.push_back((*vertices)[i+noPolygonVertices]);
       ThetaPolygon.push_back((*vertices)[i+2*noPolygonVertices-1]);
       CalculateClippedPolygonExtent(ThetaPolygon,pVoxelLimit,pAxis,pMin,pMax);
       ThetaPolygon.clear();
     }
     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 * G4Ellipsoid::Clone ( ) const

virtual

Reimplemented from G4VSolid.

Definition at line 961 of file G4Ellipsoid.cc.

References G4Ellipsoid().

 {
   return new G4Ellipsoid(*this);
 }

void G4Ellipsoid::ComputeDimensions	(	G4VPVParameterisation *	p,
		const G4int	n,
		const G4VPhysicalVolume *	pRep
	)

virtual

Reimplemented from G4VSolid.

Definition at line 177 of file G4Ellipsoid.cc.

References G4VPVParameterisation::ComputeDimensions().

 {
   p->ComputeDimensions(*this,n,pRep);
 }

G4Polyhedron * G4Ellipsoid::CreatePolyhedron ( ) const

virtual

Reimplemented from G4VSolid.

Definition at line 1077 of file G4Ellipsoid.cc.

Referenced by GetPolyhedron().

 {
   return new G4PolyhedronEllipsoid(xSemiAxis, ySemiAxis, zSemiAxis,
                                    zBottomCut, zTopCut);
 }

G4ThreeVectorList * G4Ellipsoid::CreateRotatedVertices	(	const G4AffineTransform &	pT,
		G4int &	noPV
	)		const

protected

Definition at line 836 of file G4Ellipsoid.cc.

References G4VSolid::DumpInfo(), FatalException, G4Exception(), kMeshAngleDefault, python.hepunit::pi, G4AffineTransform::TransformPoint(), and python.hepunit::twopi.

Referenced by CalculateExtent().

 {
   G4ThreeVectorList *vertices;
   G4ThreeVector vertex;
   G4double meshAnglePhi, meshRMaxFactor,
            crossAnglePhi, coscrossAnglePhi, sincrossAnglePhi, sAnglePhi;
   G4double meshTheta, crossTheta, startTheta;
   G4double rMaxX, rMaxY, rMaxZ, rMaxMax, rx, ry, rz;
   G4int crossSectionPhi, noPhiCrossSections, crossSectionTheta, noThetaSections;
 
   // Phi cross sections
   //
   noPhiCrossSections=G4int (twopi/kMeshAngleDefault)+1;  // = 9!
     
 /*
   if (noPhiCrossSections<kMinMeshSections)        // <3
   {
     noPhiCrossSections=kMinMeshSections;
   }
   else if (noPhiCrossSections>kMaxMeshSections)   // >37
   {
     noPhiCrossSections=kMaxMeshSections;
   }
 */
   meshAnglePhi=twopi/(noPhiCrossSections-1);
     
   // Set start angle such that mesh will be at fRMax
   // on the x axis. Will give better extent calculations when not rotated.
     
   sAnglePhi = -meshAnglePhi*0.5;
 
   // Theta cross sections
     
   noThetaSections = G4int(pi/kMeshAngleDefault)+3;  //  = 7!
 
 /*
   if (noThetaSections<kMinMeshSections)       // <3
   {
     noThetaSections=kMinMeshSections;
   }
   else if (noThetaSections>kMaxMeshSections)  // >37
   {
     noThetaSections=kMaxMeshSections;
   }
 */
   meshTheta= pi/(noThetaSections-2);
     
   // Set start angle such that mesh will be at fRMax
   // on the z axis. Will give better extent calculations when not rotated.
     
   startTheta = -meshTheta*0.5;
 
   meshRMaxFactor =  1.0/std::cos(0.5*
                     std::sqrt(meshAnglePhi*meshAnglePhi+meshTheta*meshTheta));
   rMaxMax= (xSemiAxis > ySemiAxis ? xSemiAxis : ySemiAxis);
   if (zSemiAxis > rMaxMax) rMaxMax= zSemiAxis;
   rMaxX= xSemiAxis + rMaxMax*(meshRMaxFactor-1.0);
   rMaxY= ySemiAxis + rMaxMax*(meshRMaxFactor-1.0);
   rMaxZ= zSemiAxis + rMaxMax*(meshRMaxFactor-1.0);
   G4double* cosCrossTheta = new G4double[noThetaSections];
   G4double* sinCrossTheta = new G4double[noThetaSections];    
   vertices=new G4ThreeVectorList(noPhiCrossSections*noThetaSections);
   if (vertices && cosCrossTheta && sinCrossTheta)
   {
     for (crossSectionTheta=0; crossSectionTheta<noThetaSections;
          crossSectionTheta++)
     {
       // Compute sine and cosine table (for historical reasons)
       //
       crossTheta=startTheta+crossSectionTheta*meshTheta;
       cosCrossTheta[crossSectionTheta]=std::cos(crossTheta);
       sinCrossTheta[crossSectionTheta]=std::sin(crossTheta);
     }
     for (crossSectionPhi=0; crossSectionPhi<noPhiCrossSections;
          crossSectionPhi++)
     {
       crossAnglePhi=sAnglePhi+crossSectionPhi*meshAnglePhi;
       coscrossAnglePhi=std::cos(crossAnglePhi);
       sincrossAnglePhi=std::sin(crossAnglePhi);
       for (crossSectionTheta=0; crossSectionTheta<noThetaSections;
            crossSectionTheta++)
       {
         // Compute coordinates of cross section at section crossSectionPhi
         //
         rx= sinCrossTheta[crossSectionTheta]*coscrossAnglePhi*rMaxX;
         ry= sinCrossTheta[crossSectionTheta]*sincrossAnglePhi*rMaxY;
         rz= cosCrossTheta[crossSectionTheta]*rMaxZ;
         if (rz < zBottomCut)
           { rz= zBottomCut; }
         if (rz > zTopCut)
           { rz= zTopCut; }
         vertex= G4ThreeVector(rx,ry,rz);
         vertices->push_back(pTransform.TransformPoint(vertex));
       }    // Theta forward     
     }    // Phi
     noPolygonVertices = noThetaSections ;
   }
   else
   {
     DumpInfo();
     G4Exception("G4Ellipsoid::CreateRotatedVertices()",
                 "GeomSolids0003", FatalException,
                 "Error in allocation of vertices. Out of memory !");
   }
 
   delete[] cosCrossTheta;
   delete[] sinCrossTheta;
 
   return vertices;
 }

void G4Ellipsoid::DescribeYourselfTo ( G4VGraphicsScene & scene ) const

virtual

Implements G4VSolid.

Definition at line 1063 of file G4Ellipsoid.cc.

References G4VGraphicsScene::AddSolid().

 {
   scene.AddSolid(*this);
 }

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

virtual

Implements G4VSolid.

Definition at line 511 of file G4Ellipsoid.cc.

References CLHEP::Hep3Vector::dot(), Inside(), kOutside, G4INCL::Math::min(), sqr(), SurfaceNormal(), CLHEP::Hep3Vector::x(), CLHEP::Hep3Vector::y(), and CLHEP::Hep3Vector::z().

 {
   G4double distMin = std::min(xSemiAxis,ySemiAxis);
   const G4double dRmax = 100.*std::min(distMin,zSemiAxis);
   distMin= kInfinity;
 
   // check to see if Z plane is relevant
   if (p.z() <= zBottomCut+halfCarTolerance)
   {
     if (v.z() <= 0.0) { return distMin; }
     G4double distZ = (zBottomCut - p.z()) / v.z();
 
     if ( (distZ > -halfRadTolerance) && (Inside(p+distZ*v) != kOutside) )
     {
       // early exit since can't intercept curved surface if we reach here
       if ( std::fabs(distZ) < halfRadTolerance ) { distZ=0.; }
       return distMin= distZ;
     }
   }
   if (p.z() >= zTopCut-halfCarTolerance)
   {
     if (v.z() >= 0.0) { return distMin;}
     G4double distZ = (zTopCut - p.z()) / v.z();
     if ( (distZ > -halfRadTolerance) && (Inside(p+distZ*v) != kOutside) )
     {
       // early exit since can't intercept curved surface if we reach here
       if ( std::fabs(distZ) < halfRadTolerance ) { distZ=0.; }
       return distMin= distZ;
     }
   }
   // if fZCut1 <= p.z() <= fZCut2, then must hit curved surface
 
   // now check curved surface intercept
   G4double A,B,C;
 
   A= sqr(v.x()/xSemiAxis) + sqr(v.y()/ySemiAxis) + sqr(v.z()/zSemiAxis);
   C= sqr(p.x()/xSemiAxis) + sqr(p.y()/ySemiAxis) + sqr(p.z()/zSemiAxis) - 1.0;
   B= 2.0 * ( p.x()*v.x()/(xSemiAxis*xSemiAxis)
            + p.y()*v.y()/(ySemiAxis*ySemiAxis)
            + p.z()*v.z()/(zSemiAxis*zSemiAxis) );
 
   C= B*B - 4.0*A*C;
   if (C > 0.0)
   {    
     G4double distR= (-B - std::sqrt(C)) / (2.0*A);
     G4double intZ = p.z()+distR*v.z();
     if ( (distR > halfRadTolerance)
       && (intZ >= zBottomCut-halfRadTolerance)
       && (intZ <= zTopCut+halfRadTolerance) )
     { 
       distMin = distR;
     }
     else if( (distR >- halfRadTolerance)
         && (intZ >= zBottomCut-halfRadTolerance)
         && (intZ <= zTopCut+halfRadTolerance) )
     {
       // p is on the curved surface, DistanceToIn returns 0 or kInfinity:
       // DistanceToIn returns 0, if second root is positive (means going inside)
       // If second root is negative, DistanceToIn returns kInfinity (outside)
       //
       distR = (-B + std::sqrt(C) ) / (2.0*A);
       if(distR>0.) { distMin=0.; }
     }
     else
     {
       distR= (-B + std::sqrt(C)) / (2.0*A);
       intZ = p.z()+distR*v.z();
       if ( (distR > halfRadTolerance)
         && (intZ >= zBottomCut-halfRadTolerance)
         && (intZ <= zTopCut+halfRadTolerance) )
       {
         G4ThreeVector norm=SurfaceNormal(p);
         if (norm.dot(v)<0.) { distMin = distR; }
       }
     }
     if ( (distMin!=kInfinity) && (distMin>dRmax) ) 
     {                    // Avoid rounding errors due to precision issues on
                          // 64 bits systems. Split long distances and recompute
       G4double fTerm = distMin-std::fmod(distMin,dRmax);
       distMin = fTerm + DistanceToIn(p+fTerm*v,v);
     }
   }
   
   if (std::fabs(distMin)<halfRadTolerance) { distMin=0.; }
   return distMin;
 } 

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

virtual

Implements G4VSolid.

Definition at line 604 of file G4Ellipsoid.cc.

References CLHEP::Hep3Vector::mag(), CLHEP::Hep3Vector::x(), CLHEP::Hep3Vector::y(), and CLHEP::Hep3Vector::z().

 {
   G4double distR, distZ;
 
   // normal vector:  parallel to normal, magnitude 1/(characteristic radius)
   //
   G4ThreeVector norm(p.x()/(xSemiAxis*xSemiAxis),
                      p.y()/(ySemiAxis*ySemiAxis),
                      p.z()/(zSemiAxis*zSemiAxis));
   G4double radius= 1.0/norm.mag();
 
   // approximate distance to curved surface ( <= actual distance )
   //
   distR= (p*norm - 1.0) * radius / 2.0;
 
   // Distance to z-cut plane
   //
   distZ= zBottomCut - p.z();
   if (distZ < 0.0)
   {
     distZ = p.z() - zTopCut;
   }
 
   // Distance to closest surface from outside
   //
   if (distZ < 0.0)
   {
     return (distR < 0.0) ? 0.0 : distR;
   }
   else if (distR < 0.0)
   {
     return distZ;
   }
   else
   {
     return (distZ < distR) ? distZ : distR;
   }
 }

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

virtual

Implements G4VSolid.

Definition at line 647 of file G4Ellipsoid.cc.

References G4VSolid::DumpInfo(), G4endl, G4Exception(), JustWarning, CLHEP::Hep3Vector::mag(), python.hepunit::mm, sqr(), test::v, CLHEP::Hep3Vector::x(), CLHEP::Hep3Vector::y(), and CLHEP::Hep3Vector::z().

 {
   G4double distMin;
   enum surface_e {kPlaneSurf, kCurvedSurf, kNoSurf} surface;
   
   distMin= kInfinity;
   surface= kNoSurf;
 
   // check to see if Z plane is relevant
   //
   if (v.z() < 0.0)
   {
     G4double distZ = (zBottomCut - p.z()) / v.z();
     if (distZ < 0.0)
     {
       distZ= 0.0;
       if (!calcNorm) {return 0.0;}
     }
     distMin= distZ;
     surface= kPlaneSurf;
   }
   if (v.z() > 0.0)
   {
     G4double distZ = (zTopCut - p.z()) / v.z();
     if (distZ < 0.0)
     {
       distZ= 0.0;
       if (!calcNorm) {return 0.0;}
     }
     distMin= distZ;
     surface= kPlaneSurf;
   }
 
   // normal vector:  parallel to normal, magnitude 1/(characteristic radius)
   //
   G4ThreeVector nearnorm(p.x()/(xSemiAxis*xSemiAxis),
                          p.y()/(ySemiAxis*ySemiAxis),
                          p.z()/(zSemiAxis*zSemiAxis));
   
   // now check curved surface intercept
   //
   G4double A,B,C;
   
   A= sqr(v.x()/xSemiAxis) + sqr(v.y()/ySemiAxis) + sqr(v.z()/zSemiAxis);
   C= (p * nearnorm) - 1.0;
   B= 2.0 * (v * nearnorm);
 
   C= B*B - 4.0*A*C;
   if (C > 0.0)
   {
     G4double distR= (-B + std::sqrt(C) ) / (2.0*A);
     if (distR < 0.0)
     {
       distR= 0.0;
       if (!calcNorm) {return 0.0;}
     }
     if (distR < distMin)
     {
       distMin= distR;
       surface= kCurvedSurf;
     }
   }
 
   // set normal if requested
   //
   if (calcNorm)
   {
     if (surface == kNoSurf)
     {
       *validNorm = false;
     }
     else
     {
       *validNorm = true;
       switch (surface)
       {
         case kPlaneSurf:
           *n= G4ThreeVector(0.,0.,(v.z() > 0.0 ? 1. : -1.));
           break;
         case kCurvedSurf:
         {
           G4ThreeVector pexit= p + distMin*v;
           G4ThreeVector truenorm(pexit.x()/(xSemiAxis*xSemiAxis),
                                  pexit.y()/(ySemiAxis*ySemiAxis),
                                  pexit.z()/(zSemiAxis*zSemiAxis));
           truenorm *= 1.0/truenorm.mag();
           *n= truenorm;
         } break;
         default:           // Should never reach this case ...
           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 << 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("G4Ellipsoid::DistanceToOut(p,v,..)",
                       "GeomSolids1002", JustWarning, message);
           break;
       }
     }
   }
    
   return distMin;
 }

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

virtual

Implements G4VSolid.

Definition at line 770 of file G4Ellipsoid.cc.

References G4VSolid::DumpInfo(), G4endl, G4Exception(), Inside(), JustWarning, kOutside, CLHEP::Hep3Vector::mag(), python.hepunit::mm, CLHEP::Hep3Vector::x(), CLHEP::Hep3Vector::y(), and CLHEP::Hep3Vector::z().

 {
   G4double distR, distZ;
 
 #ifdef G4SPECSDEBUG
   if( Inside(p) == kOutside )
   {
      DumpInfo();
      std::ostringstream message;
      G4int oldprc = message.precision(16);
      message << "Point p is outside !?" << G4endl
              << "Position:"  << G4endl
              << "   p.x() = "   << p.x()/mm << " mm" << G4endl
              << "   p.y() = "   << p.y()/mm << " mm" << G4endl
              << "   p.z() = "   << p.z()/mm << " mm";
      message.precision(oldprc) ;
      G4Exception("G4Ellipsoid::DistanceToOut(p)", "GeomSolids1002",
                  JustWarning, message);
   }
 #endif
 
   // Normal vector:  parallel to normal, magnitude 1/(characteristic radius)
   //
   G4ThreeVector norm(p.x()/(xSemiAxis*xSemiAxis),
                      p.y()/(ySemiAxis*ySemiAxis),
                      p.z()/(zSemiAxis*zSemiAxis));
 
   // the following is a safe inlined "radius= min(1.0/norm.mag(),p.mag())
   //
   G4double radius= p.mag();
   G4double tmp= norm.mag();
   if ( (tmp > 0.0) && (1.0 < radius*tmp) ) {radius = 1.0/tmp;}
 
   // Approximate distance to curved surface ( <= actual distance )
   //
   distR = (1.0 - p*norm) * radius / 2.0;
     
   // Distance to z-cut plane
   //
   distZ = p.z() - zBottomCut;
   if (distZ < 0.0) {distZ= zTopCut - p.z();}
 
   // Distance to closest surface from inside
   //
   if ( (distZ < 0.0) || (distR < 0.0) )
   {
     return 0.0;
   }
   else
   {
     return (distZ < distR) ? distZ : distR;
   }
 }

G4double G4Ellipsoid::GetCubicVolume ( )

inlinevirtual

Reimplemented from G4VSolid.

G4GeometryType G4Ellipsoid::GetEntityType ( ) const

virtual

Implements G4VSolid.

Definition at line 952 of file G4Ellipsoid.cc.

 {
   return G4String("G4Ellipsoid");
 }

G4VisExtent G4Ellipsoid::GetExtent ( ) const

virtual

Reimplemented from G4VSolid.

Definition at line 1068 of file G4Ellipsoid.cc.

 {
   // Define the sides of the box into which the G4Ellipsoid instance would fit.
   //
   return G4VisExtent (-semiAxisMax, semiAxisMax,
                       -semiAxisMax, semiAxisMax,
                       -semiAxisMax, semiAxisMax);
 }

G4ThreeVector G4Ellipsoid::GetPointOnSurface ( ) const

virtual

Reimplemented from G4VSolid.

Definition at line 995 of file G4Ellipsoid.cc.

References python.hepunit::pi, G4INCL::DeJongSpin::shoot(), sqr(), and python.hepunit::twopi.

 {
   G4double aTop, aBottom, aCurved, chose, xRand, yRand, zRand, phi;
   G4double cosphi, sinphi, costheta, sintheta, alpha, beta, max1, max2, max3;
 
   max1  = xSemiAxis > ySemiAxis ? xSemiAxis : ySemiAxis;
   max1  = max1 > zSemiAxis ? max1 : zSemiAxis;
   if (max1 == xSemiAxis)      { max2 = ySemiAxis; max3 = zSemiAxis; }
   else if (max1 == ySemiAxis) { max2 = xSemiAxis; max3 = zSemiAxis; }
   else                        { max2 = xSemiAxis; max3 = ySemiAxis; }
 
   phi   = RandFlat::shoot(0.,twopi);
   
   cosphi = std::cos(phi);   sinphi = std::sin(phi);
   costheta = RandFlat::shoot(zBottomCut,zTopCut)/zSemiAxis;
   sintheta = std::sqrt(1.-sqr(costheta));
   
   alpha = 1.-sqr(max2/max1); beta  = 1.-sqr(max3/max1);
   
   aTop    = pi*xSemiAxis*ySemiAxis*(1 - sqr(zTopCut/zSemiAxis));
   aBottom = pi*xSemiAxis*ySemiAxis*(1 - sqr(zBottomCut/zSemiAxis));
   
   // approximation
   // from:" http://www.citr.auckland.ac.nz/techreports/2004/CITR-TR-139.pdf"
   aCurved = 4.*pi*max1*max2*(1.-1./6.*(alpha+beta)-
                             1./120.*(3.*sqr(alpha)+2.*alpha*beta+3.*sqr(beta)));
 
   aCurved *= 0.5*(1.2*zTopCut/zSemiAxis - 1.2*zBottomCut/zSemiAxis);
   
   if( ( zTopCut >= zSemiAxis && zBottomCut <= -1.*zSemiAxis )
    || ( zTopCut == 0 && zBottomCut ==0 ) )
   {
     aTop = 0; aBottom = 0;
   }
   
   chose = RandFlat::shoot(0.,aTop + aBottom + aCurved); 
   
   if(chose < aCurved)
   { 
     xRand = xSemiAxis*sintheta*cosphi;
     yRand = ySemiAxis*sintheta*sinphi;
     zRand = zSemiAxis*costheta;
     return G4ThreeVector (xRand,yRand,zRand); 
   }
   else if(chose >= aCurved && chose < aCurved + aTop)
   {
     xRand = RandFlat::shoot(-1.,1.)*xSemiAxis
           * std::sqrt(1-sqr(zTopCut/zSemiAxis));
     yRand = RandFlat::shoot(-1.,1.)*ySemiAxis
           * std::sqrt(1.-sqr(zTopCut/zSemiAxis)-sqr(xRand/xSemiAxis));
     zRand = zTopCut;
     return G4ThreeVector (xRand,yRand,zRand);
   }
   else
   {
     xRand = RandFlat::shoot(-1.,1.)*xSemiAxis
           * std::sqrt(1-sqr(zBottomCut/zSemiAxis));
     yRand = RandFlat::shoot(-1.,1.)*ySemiAxis
           * std::sqrt(1.-sqr(zBottomCut/zSemiAxis)-sqr(xRand/xSemiAxis)); 
     zRand = zBottomCut;
     return G4ThreeVector (xRand,yRand,zRand);
   }
 }

G4Polyhedron * G4Ellipsoid::GetPolyhedron ( ) const

virtual

Reimplemented from G4VSolid.

Definition at line 1083 of file G4Ellipsoid.cc.

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

 {
   if (!fpPolyhedron ||
       fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() !=
       fpPolyhedron->GetNumberOfRotationSteps())
     {
       delete fpPolyhedron;
       fpPolyhedron = CreatePolyhedron();
     }
   return fpPolyhedron;
 }

G4double G4Ellipsoid::GetSemiAxisMax ( G4int i ) const

inline

Referenced by G4GDMLWriteParamvol::Ellipsoid_dimensionsWrite(), G4GDMLWriteSolids::EllipsoidWrite(), export_G4Ellipsoid(), and G4tgbGeometryDumper::GetSolidParams().

G4double G4Ellipsoid::GetSurfaceArea ( )

inlinevirtual

Reimplemented from G4VSolid.

G4double G4Ellipsoid::GetZBottomCut ( ) const

inline

Referenced by G4GDMLWriteParamvol::Ellipsoid_dimensionsWrite(), G4GDMLWriteSolids::EllipsoidWrite(), export_G4Ellipsoid(), and G4tgbGeometryDumper::GetSolidParams().

G4double G4Ellipsoid::GetZTopCut ( ) const

inline

Referenced by G4GDMLWriteParamvol::Ellipsoid_dimensionsWrite(), G4GDMLWriteSolids::EllipsoidWrite(), export_G4Ellipsoid(), and G4tgbGeometryDumper::GetSolidParams().

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

virtual

Implements G4VSolid.

Definition at line 435 of file G4Ellipsoid.cc.

References kInside, kOutside, kSurface, sqr(), CLHEP::Hep3Vector::x(), CLHEP::Hep3Vector::y(), and CLHEP::Hep3Vector::z().

Referenced by CalculateExtent(), DistanceToIn(), and DistanceToOut().

 {
   G4double rad2oo,  // outside surface outer tolerance
            rad2oi;  // outside surface inner tolerance
   EInside in;
 
   // check this side of z cut first, because that's fast
   //
   if (p.z() < zBottomCut-halfRadTolerance) { return in=kOutside; }
   if (p.z() > zTopCut+halfRadTolerance)    { return in=kOutside; }
 
   rad2oo= sqr(p.x()/(xSemiAxis+halfRadTolerance))
         + sqr(p.y()/(ySemiAxis+halfRadTolerance))
         + sqr(p.z()/(zSemiAxis+halfRadTolerance));
 
   if (rad2oo > 1.0)  { return in=kOutside; }
     
   rad2oi= sqr(p.x()*(1.0+halfRadTolerance/xSemiAxis)/xSemiAxis)
       + sqr(p.y()*(1.0+halfRadTolerance/ySemiAxis)/ySemiAxis)
       + sqr(p.z()*(1.0+halfRadTolerance/zSemiAxis)/zSemiAxis);
 
   // Check radial surfaces
   //  sets `in' (already checked for rad2oo > 1.0)
   //
   if (rad2oi < 1.0)
   {
     in = ( (p.z() < zBottomCut+halfRadTolerance)
         || (p.z() > zTopCut-halfRadTolerance) ) ? kSurface : kInside;
     if ( rad2oi > 1.0-halfRadTolerance )  { in=kSurface; }
   }
   else 
   {
     in = kSurface;
   }
   return in;
 
 }

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

Definition at line 149 of file G4Ellipsoid.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; kRadTolerance = rhs.kRadTolerance;
    halfCarTolerance = rhs.halfCarTolerance;
    halfRadTolerance = rhs.halfRadTolerance;
    fCubicVolume = rhs.fCubicVolume; fSurfaceArea = rhs.fSurfaceArea;
    xSemiAxis = rhs.xSemiAxis; ySemiAxis = rhs.ySemiAxis;
    zSemiAxis = rhs.zSemiAxis; semiAxisMax = rhs.semiAxisMax;
    zBottomCut = rhs.zBottomCut; zTopCut = rhs.zTopCut;
 
    return *this;
 }

void G4Ellipsoid::SetSemiAxis	(	G4double	x,
		G4double	y,
		G4double	z
	)

inline

Referenced by export_G4Ellipsoid(), and G4Ellipsoid().

void G4Ellipsoid::SetZCuts	(	G4double	newzBottomCut,
		G4double	newzTopCut
	)

inline

Referenced by export_G4Ellipsoid(), and G4Ellipsoid().

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

virtual

Implements G4VSolid.

Definition at line 970 of file G4Ellipsoid.cc.

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

 {
   G4int oldprc = os.precision(16);
   os << "-----------------------------------------------------------\n"
      << "    *** Dump for solid - " << GetName() << " ***\n"
      << "    ===================================================\n"
      << " Solid type: G4Ellipsoid\n"
      << " Parameters: \n"
 
      << "    semi-axis x: " << xSemiAxis/mm << " mm \n"
      << "    semi-axis y: " << ySemiAxis/mm << " mm \n"
      << "    semi-axis z: " << zSemiAxis/mm << " mm \n"
      << "    max semi-axis: " << semiAxisMax/mm << " mm \n"
      << "    lower cut plane level z: " << zBottomCut/mm << " mm \n"
      << "    upper cut plane level z: " << zTopCut/mm << " mm \n"
      << "-----------------------------------------------------------\n";
   os.precision(oldprc);
 
   return os;
 }

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

virtual

Implements G4VSolid.

Definition at line 477 of file G4Ellipsoid.cc.

References CLHEP::Hep3Vector::mag(), CLHEP::Hep3Vector::x(), CLHEP::Hep3Vector::y(), and CLHEP::Hep3Vector::z().

Referenced by DistanceToIn().

 {
   G4double distR, distZBottom, distZTop;
 
   // normal vector with special magnitude:  parallel to normal, units 1/length
   // norm*p == 1.0 if on surface, >1.0 if outside, <1.0 if inside
   //
   G4ThreeVector norm(p.x()/(xSemiAxis*xSemiAxis),
                      p.y()/(ySemiAxis*ySemiAxis),
                      p.z()/(zSemiAxis*zSemiAxis));
   G4double radius = 1.0/norm.mag();
 
   // approximate distance to curved surface
   //
   distR = std::fabs( (p*norm - 1.0) * radius ) / 2.0;
 
   // Distance to z-cut plane
   //
   distZBottom = std::fabs( p.z() - zBottomCut );
   distZTop = std::fabs( p.z() - zTopCut );
 
   if ( (distZBottom < distR) || (distZTop < distR) )
   {
     return G4ThreeVector(0.,0.,(distZBottom < distZTop) ? -1.0 : 1.0);
   }
   return ( norm *= radius );
 }

Field Documentation

G4Polyhedron* G4Ellipsoid::fpPolyhedron

mutableprotected

Definition at line 135 of file G4Ellipsoid.hh.

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

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

Public Member Functions

Protected Member Functions

Protected Attributes

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation

Field Documentation