#include <G4Ellipsoid.hh>

Inheritance diagram for G4Ellipsoid:

Collaboration 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)

G4double	EstimateCubicVolume (G4int nStat, G4double epsilon) const

G4double	EstimateSurfaceArea (G4int nStat, G4double ell) const

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

Protected Attributes
G4bool	fRebuildPolyhedron

G4Polyhedron *	fpPolyhedron

Protected Attributes inherited from G4VSolid
G4double	kCarTolerance

Detailed Description

Definition at line 60 of file G4Ellipsoid.hh.

Constructor & Destructor Documentation

◆ G4Ellipsoid() [1/3]

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

Definition at line 69 of file G4Ellipsoid.cc.

   : G4VSolid(pName), fRebuildPolyhedron(false), 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);
   }
 }

Here is the call graph for this function:

Here is the caller graph for this function:

◆ ~G4Ellipsoid()

G4Ellipsoid::~G4Ellipsoid ( )

virtual

Definition at line 131 of file G4Ellipsoid.cc.

 {
   delete fpPolyhedron; fpPolyhedron = 0;
 }

◆ G4Ellipsoid() [2/3]

G4Ellipsoid::G4Ellipsoid ( __void__ & a )

Definition at line 119 of file G4Ellipsoid.cc.

   : G4VSolid(a), fRebuildPolyhedron(false), 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() [3/3]

G4Ellipsoid::G4Ellipsoid ( const G4Ellipsoid & rhs )

Definition at line 140 of file G4Ellipsoid.cc.

   : G4VSolid(rhs),
     fRebuildPolyhedron(false), 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

◆ CalculateExtent()

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

virtual

Implements G4VSolid.

Definition at line 199 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;
   }
 }

Here is the call graph for this function:

◆ Clone()

G4VSolid * G4Ellipsoid::Clone ( ) const

virtual

Reimplemented from G4VSolid.

Definition at line 971 of file G4Ellipsoid.cc.

 {
   return new G4Ellipsoid(*this);
 }

Here is the call graph for this function:

◆ ComputeDimensions()

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

virtual

Reimplemented from G4VSolid.

Definition at line 187 of file G4Ellipsoid.cc.

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

Here is the call graph for this function:

◆ CreatePolyhedron()

G4Polyhedron * G4Ellipsoid::CreatePolyhedron ( ) const

virtual

Reimplemented from G4VSolid.

Definition at line 1087 of file G4Ellipsoid.cc.

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

Here is the caller graph for this function:

◆ CreateRotatedVertices()

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

protected

Definition at line 846 of file G4Ellipsoid.cc.

 {
   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;
 }

Here is the call graph for this function:

Here is the caller graph for this function:

◆ DescribeYourselfTo()

void G4Ellipsoid::DescribeYourselfTo ( G4VGraphicsScene & scene ) const

virtual

Implements G4VSolid.

Definition at line 1073 of file G4Ellipsoid.cc.

 {
   scene.AddSolid(*this);
 }

Here is the call graph for this function:

◆ DistanceToIn() [1/2]

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

virtual

Implements G4VSolid.

Definition at line 521 of file G4Ellipsoid.cc.

 {
   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;
 } 

Here is the call graph for this function:

◆ DistanceToIn() [2/2]

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

virtual

Implements G4VSolid.

Definition at line 614 of file G4Ellipsoid.cc.

 {
   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;
   }
 }

Here is the call graph for this function:

◆ DistanceToOut() [1/2]

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 657 of file G4Ellipsoid.cc.

 {
   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;
 }

Here is the call graph for this function:

◆ DistanceToOut() [2/2]

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

virtual

Implements G4VSolid.

Definition at line 780 of file G4Ellipsoid.cc.

 {
   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;
   }
 }

Here is the call graph for this function:

◆ GetCubicVolume()

G4double G4Ellipsoid::GetCubicVolume ( )

inlinevirtual

Reimplemented from G4VSolid.

◆ GetEntityType()

G4GeometryType G4Ellipsoid::GetEntityType ( ) const

virtual

Implements G4VSolid.

Definition at line 962 of file G4Ellipsoid.cc.

 {
   return G4String("G4Ellipsoid");
 }

◆ GetExtent()

G4VisExtent G4Ellipsoid::GetExtent ( ) const

virtual

Reimplemented from G4VSolid.

Definition at line 1078 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);
 }

◆ GetPointOnSurface()

G4ThreeVector G4Ellipsoid::GetPointOnSurface ( ) const

virtual

Reimplemented from G4VSolid.

Definition at line 1005 of file G4Ellipsoid.cc.

 {
   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   = G4RandFlat::shoot(0.,twopi);
   
   cosphi = std::cos(phi);   sinphi = std::sin(phi);
   costheta = G4RandFlat::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 = G4RandFlat::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 = G4RandFlat::shoot(-1.,1.)*xSemiAxis
           * std::sqrt(1-sqr(zTopCut/zSemiAxis));
     yRand = G4RandFlat::shoot(-1.,1.)*ySemiAxis
           * std::sqrt(1.-sqr(zTopCut/zSemiAxis)-sqr(xRand/xSemiAxis));
     zRand = zTopCut;
     return G4ThreeVector (xRand,yRand,zRand);
   }
   else
   {
     xRand = G4RandFlat::shoot(-1.,1.)*xSemiAxis
           * std::sqrt(1-sqr(zBottomCut/zSemiAxis));
     yRand = G4RandFlat::shoot(-1.,1.)*ySemiAxis
           * std::sqrt(1.-sqr(zBottomCut/zSemiAxis)-sqr(xRand/xSemiAxis)); 
     zRand = zBottomCut;
     return G4ThreeVector (xRand,yRand,zRand);
   }
 }

Here is the call graph for this function:

◆ GetPolyhedron()

G4Polyhedron * G4Ellipsoid::GetPolyhedron ( ) const

virtual

Reimplemented from G4VSolid.

Definition at line 1093 of file G4Ellipsoid.cc.

 {
   if (!fpPolyhedron ||
       fRebuildPolyhedron ||
       fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() !=
       fpPolyhedron->GetNumberOfRotationSteps())
     {
       G4AutoLock l(&polyhedronMutex);
       delete fpPolyhedron;
       fpPolyhedron = CreatePolyhedron();
       fRebuildPolyhedron = false;
       l.unlock();
     }
   return fpPolyhedron;
 }

Here is the call graph for this function:

◆ GetSemiAxisMax()

G4double G4Ellipsoid::GetSemiAxisMax ( G4int i ) const

inline

Here is the caller graph for this function:

◆ GetSurfaceArea()

G4double G4Ellipsoid::GetSurfaceArea ( )

inlinevirtual

Reimplemented from G4VSolid.

◆ GetZBottomCut()

G4double G4Ellipsoid::GetZBottomCut ( ) const

inline

Here is the caller graph for this function:

◆ GetZTopCut()

G4double G4Ellipsoid::GetZTopCut ( ) const

inline

Here is the caller graph for this function:

◆ Inside()

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

virtual

Implements G4VSolid.

Definition at line 445 of file G4Ellipsoid.cc.

 {
   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;
 
 }

Here is the call graph for this function:

Here is the caller graph for this function:

◆ operator=()

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

Definition at line 157 of file G4Ellipsoid.cc.

 {
    // Check assignment to self
    //
    if (this == &rhs)  { return *this; }
 
    // Copy base class data
    //
    G4VSolid::operator=(rhs);
 
    // Copy data
    //
    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;
    fRebuildPolyhedron = false;
    delete fpPolyhedron; fpPolyhedron = 0;
 
    return *this;
 }

Here is the call graph for this function:

◆ SetSemiAxis()

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

inline

Here is the caller graph for this function:

◆ SetZCuts()

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

inline

Here is the caller graph for this function:

◆ StreamInfo()

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

virtual

Implements G4VSolid.

Definition at line 980 of file G4Ellipsoid.cc.

 {
   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;
 }

Here is the call graph for this function:

◆ SurfaceNormal()

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

virtual

Implements G4VSolid.

Definition at line 487 of file G4Ellipsoid.cc.

 {
   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 );
 }

Here is the call graph for this function:

Here is the caller graph for this function:

Member Data Documentation

◆ fCubicVolume

G4double G4Ellipsoid::fCubicVolume

private

Definition at line 144 of file G4Ellipsoid.hh.

◆ fpPolyhedron

G4Polyhedron* G4Ellipsoid::fpPolyhedron

mutableprotected

Definition at line 137 of file G4Ellipsoid.hh.

◆ fRebuildPolyhedron

G4bool G4Ellipsoid::fRebuildPolyhedron

mutableprotected

Definition at line 136 of file G4Ellipsoid.hh.

◆ fSurfaceArea

G4double G4Ellipsoid::fSurfaceArea

private

Definition at line 145 of file G4Ellipsoid.hh.

◆ halfCarTolerance

G4double G4Ellipsoid::halfCarTolerance

private

Definition at line 142 of file G4Ellipsoid.hh.

◆ halfRadTolerance

G4double G4Ellipsoid::halfRadTolerance

private

Definition at line 142 of file G4Ellipsoid.hh.

◆ kRadTolerance

G4double G4Ellipsoid::kRadTolerance

private

Definition at line 141 of file G4Ellipsoid.hh.

◆ semiAxisMax

G4double G4Ellipsoid::semiAxisMax

private

Definition at line 146 of file G4Ellipsoid.hh.

◆ xSemiAxis

G4double G4Ellipsoid::xSemiAxis

private

Definition at line 146 of file G4Ellipsoid.hh.

◆ ySemiAxis

G4double G4Ellipsoid::ySemiAxis

private

Definition at line 146 of file G4Ellipsoid.hh.

◆ zBottomCut

G4double G4Ellipsoid::zBottomCut

private

Definition at line 146 of file G4Ellipsoid.hh.

◆ zSemiAxis

G4double G4Ellipsoid::zSemiAxis

private

Definition at line 146 of file G4Ellipsoid.hh.

◆ zTopCut

G4double G4Ellipsoid::zTopCut

private

Definition at line 146 of file G4Ellipsoid.hh.

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

Private Attributes
G4double	kRadTolerance

G4double	halfCarTolerance

G4double	halfRadTolerance

G4double	fCubicVolume

G4double	fSurfaceArea

G4double	xSemiAxis

G4double	ySemiAxis

G4double	zSemiAxis

G4double	semiAxisMax

G4double	zBottomCut

G4double	zTopCut

Public Member Functions

Protected Member Functions

Protected Attributes

Private Attributes

Detailed Description

Constructor & Destructor Documentation

◆ G4Ellipsoid() [1/3]

◆ ~G4Ellipsoid()

◆ G4Ellipsoid() [2/3]

◆ G4Ellipsoid() [3/3]

Member Function Documentation

◆ CalculateExtent()

◆ Clone()

◆ ComputeDimensions()

◆ CreatePolyhedron()

◆ CreateRotatedVertices()

◆ DescribeYourselfTo()

◆ DistanceToIn() [1/2]

◆ DistanceToIn() [2/2]

◆ DistanceToOut() [1/2]

◆ DistanceToOut() [2/2]

◆ GetCubicVolume()

◆ GetEntityType()

◆ GetExtent()

◆ GetPointOnSurface()

◆ GetPolyhedron()

◆ GetSemiAxisMax()

◆ GetSurfaceArea()

◆ GetZBottomCut()

◆ GetZTopCut()

◆ Inside()

◆ operator=()

◆ SetSemiAxis()

◆ SetZCuts()

◆ StreamInfo()

◆ SurfaceNormal()

Member Data Documentation

◆ fCubicVolume

◆ fpPolyhedron

◆ fRebuildPolyhedron

◆ fSurfaceArea

◆ halfCarTolerance

◆ halfRadTolerance

◆ kRadTolerance

◆ semiAxisMax

◆ xSemiAxis

◆ ySemiAxis

◆ zBottomCut

◆ zSemiAxis

◆ zTopCut