#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 ()

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

G4NURBS *	CreateNURBS () 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

virtual G4bool	CalculateExtent (const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pMin, G4double &pMax) const =0

virtual EInside	Inside (const G4ThreeVector &p) const =0

virtual G4ThreeVector	SurfaceNormal (const G4ThreeVector &p) const =0

virtual G4double	DistanceToIn (const G4ThreeVector &p, const G4ThreeVector &v) const =0

virtual G4double	DistanceToIn (const G4ThreeVector &p) const =0

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

virtual G4double	DistanceToOut (const G4ThreeVector &p) const =0

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

virtual G4double	GetCubicVolume ()

virtual G4double	GetSurfaceArea ()

virtual G4GeometryType	GetEntityType () const =0

virtual G4ThreeVector	GetPointOnSurface () const

virtual G4VSolid *	Clone () const

virtual std::ostream &	StreamInfo (std::ostream &os) const =0

void	DumpInfo () const

virtual void	DescribeYourselfTo (G4VGraphicsScene &scene) const =0

virtual G4VisExtent	GetExtent () const

virtual G4Polyhedron *	CreatePolyhedron () const

virtual G4NURBS *	CreateNURBS () const

virtual G4Polyhedron *	GetPolyhedron () 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() [1/3]

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

Definition at line 64 of file G4Ellipsoid.cc.

  : 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();
 
  // 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::~G4Ellipsoid ( )

virtual

Definition at line 122 of file G4Ellipsoid.cc.

123{

124}

◆ G4Ellipsoid() [2/3]

G4Ellipsoid::G4Ellipsoid ( __void__ & a )

Definition at line 111 of file G4Ellipsoid.cc.

  : G4VSolid(a), fpPolyhedron(0), kRadTolerance(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 130 of file G4Ellipsoid.cc.

  : G4VSolid(rhs),
    fpPolyhedron(0), kRadTolerance(rhs.kRadTolerance),
    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 170 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;
  }
}

◆ Clone()

G4VSolid * G4Ellipsoid::Clone ( ) const

virtual

Reimplemented from G4VSolid.

Definition at line 947 of file G4Ellipsoid.cc.

{
  return new G4Ellipsoid(*this);
}

◆ CreateNURBS()

G4NURBS * G4Ellipsoid::CreateNURBS ( ) const

virtual

Reimplemented from G4VSolid.

Definition at line 1063 of file G4Ellipsoid.cc.

{
  // Box for now!!!
  //
  return new G4NURBSbox(semiAxisMax, semiAxisMax, semiAxisMax);
}

◆ CreatePolyhedron()

G4Polyhedron * G4Ellipsoid::CreatePolyhedron ( ) const

virtual

Reimplemented from G4VSolid.

Definition at line 1070 of file G4Ellipsoid.cc.

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

Referenced by GetPolyhedron().

◆ CreateRotatedVertices()

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

protected

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

Referenced by CalculateExtent().

◆ DescribeYourselfTo()

void G4Ellipsoid::DescribeYourselfTo ( G4VGraphicsScene & scene ) const

virtual

Implements G4VSolid.

Definition at line 1049 of file G4Ellipsoid.cc.

{
  scene.AddSolid(*this);
}

◆ DistanceToIn() [1/2]

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

virtual

Implements G4VSolid.

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

◆ DistanceToIn() [2/2]

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

virtual

Implements G4VSolid.

Definition at line 494 of file G4Ellipsoid.cc.

{
  static const G4double halfCarTolerance=kCarTolerance*0.5;
  static const G4double halfRadTolerance=kRadTolerance*0.5;
 
  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;
} 

Referenced by DistanceToIn().

◆ DistanceToOut() [1/2]

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

virtual

Implements G4VSolid.

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

◆ DistanceToOut() [2/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 633 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;
}

◆ GetCubicVolume()

G4double G4Ellipsoid::GetCubicVolume ( )

inlinevirtual

Reimplemented from G4VSolid.

◆ GetEntityType()

G4GeometryType G4Ellipsoid::GetEntityType ( ) const

virtual

Implements G4VSolid.

Definition at line 938 of file G4Ellipsoid.cc.

{
  return G4String("G4Ellipsoid");
}

◆ GetExtent()

G4VisExtent G4Ellipsoid::GetExtent ( ) const

virtual

Reimplemented from G4VSolid.

Definition at line 1054 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 981 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   = 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);
  }
}

◆ GetPolyhedron()

G4Polyhedron * G4Ellipsoid::GetPolyhedron ( ) const

virtual

Reimplemented from G4VSolid.

Definition at line 1076 of file G4Ellipsoid.cc.

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

◆ GetSemiAxisMax()

G4double G4Ellipsoid::GetSemiAxisMax ( G4int i ) const

inline

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

◆ GetSurfaceArea()

G4double G4Ellipsoid::GetSurfaceArea ( )

inlinevirtual

Reimplemented from G4VSolid.

◆ GetZBottomCut()

G4double G4Ellipsoid::GetZBottomCut ( ) const

inline

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

◆ GetZTopCut()

G4double G4Ellipsoid::GetZTopCut ( ) const

inline

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

◆ Inside()

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

virtual

Implements G4VSolid.

Definition at line 416 of file G4Ellipsoid.cc.

{
  G4double rad2oo,  // outside surface outer tolerance
           rad2oi;  // outside surface inner tolerance
  EInside in;
 
  static const G4double halfRadTolerance=kRadTolerance*0.5;
 
  // 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;
 
}

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

◆ operator=()

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

Definition at line 144 of file G4Ellipsoid.cc.

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

◆ SetSemiAxis()

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

inline

Referenced by G4Ellipsoid().

◆ SetZCuts()

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

inline

Referenced by G4Ellipsoid().

◆ StreamInfo()

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

virtual

Implements G4VSolid.

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

◆ SurfaceNormal()

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

virtual

Implements G4VSolid.

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

Referenced by DistanceToIn().

Member Data Documentation

◆ fpPolyhedron

G4Polyhedron* G4Ellipsoid::fpPolyhedron

mutableprotected

Definition at line 132 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

◆ G4Ellipsoid() [1/3]

◆ ~G4Ellipsoid()

◆ G4Ellipsoid() [2/3]

◆ G4Ellipsoid() [3/3]

Member Function Documentation

◆ CalculateExtent()

◆ Clone()

◆ CreateNURBS()

◆ 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

◆ fpPolyhedron