G4EllipticalCone Class Reference

#include <G4EllipticalCone.hh>

Inheritance diagram for G4EllipticalCone:

Collaboration diagram for G4EllipticalCone:

Public Member Functions
	G4EllipticalCone (const G4String &pName, G4double pxSemiAxis, G4double pySemiAxis, G4double zMax, G4double pzTopCut)
virtual	~G4EllipticalCone ()
G4double	GetSemiAxisMax () const
G4double	GetSemiAxisX () const
G4double	GetSemiAxisY () const
G4double	GetZMax () const
G4double	GetZTopCut () const
void	SetSemiAxis (G4double x, G4double y, G4double z)
void	SetZCut (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
G4ThreeVector	GetPointOnSurface () const
std::ostream &	StreamInfo (std::ostream &os) const
G4Polyhedron *	GetPolyhedron () const
void	DescribeYourselfTo (G4VGraphicsScene &scene) const
G4VisExtent	GetExtent () const
G4Polyhedron *	CreatePolyhedron () const
	G4EllipticalCone (__void__ &)
	G4EllipticalCone (const G4EllipticalCone &rhs)
G4EllipticalCone &	operator= (const G4EllipticalCone &rhs)
Public Member Functions inherited from G4VSolid
	G4VSolid (const G4String &name)
virtual	~G4VSolid ()
G4bool	operator== (const G4VSolid &s) const
G4String	GetName () const
void	SetName (const G4String &name)
G4double	GetTolerance () const
virtual void	ComputeDimensions (G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
void	DumpInfo () const
virtual const G4VSolid *	GetConstituentSolid (G4int no) const
virtual G4VSolid *	GetConstituentSolid (G4int no)
virtual const G4DisplacedSolid *	GetDisplacedSolidPtr () const
virtual G4DisplacedSolid *	GetDisplacedSolidPtr ()
	G4VSolid (__void__ &)
	G4VSolid (const G4VSolid &rhs)
G4VSolid &	operator= (const G4VSolid &rhs)
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

Private Attributes
G4double	kRadTolerance
G4double	halfRadTol
G4double	halfCarTol
G4double	fCubicVolume
G4double	fSurfaceArea
G4double	xSemiAxis
G4double	ySemiAxis
G4double	zheight
G4double	semiAxisMax
G4double	zTopCut

Detailed Description

Definition at line 85 of file G4EllipticalCone.hh.

Constructor & Destructor Documentation

G4EllipticalCone() [1/3]

G4EllipticalCone::G4EllipticalCone	(	const G4String &	pName,
		G4double	pxSemiAxis,
		G4double	pySemiAxis,
		G4double	zMax,
		G4double	pzTopCut
	)

Definition at line 70 of file G4EllipticalCone.cc.

  : G4VSolid(pName), fRebuildPolyhedron(false), fpPolyhedron(0),
    fCubicVolume(0.), fSurfaceArea(0.), zTopCut(0.)
{

  kRadTolerance = G4GeometryTolerance::GetInstance()->GetRadialTolerance();

  halfRadTol = 0.5*kRadTolerance;
  halfCarTol = 0.5*kCarTolerance;

  // Check Semi-Axis & Z-cut
  //
  if ( (pxSemiAxis <= 0.) || (pySemiAxis <= 0.) || (pzMax <= 0.) )
  {
     std::ostringstream message;
     message << "Invalid semi-axis or height - " << GetName();
     G4Exception("G4EllipticalCone::G4EllipticalCone()", "GeomSolids0002",
                 FatalErrorInArgument, message);
  }
  if ( pzTopCut <= 0 )
  {
     std::ostringstream message;
     message << "Invalid z-coordinate for cutting plane - " << GetName();
     G4Exception("G4EllipticalCone::G4EllipticalCone()", "InvalidSetup",
                 FatalErrorInArgument, message);
  }

  SetSemiAxis( pxSemiAxis, pySemiAxis, pzMax );
  SetZCut(pzTopCut);
}

Here is the call graph for this function:

Here is the caller graph for this function:

~G4EllipticalCone()

G4EllipticalCone::~G4EllipticalCone ( )

virtual

Definition at line 122 of file G4EllipticalCone.cc.

{
  delete fpPolyhedron; fpPolyhedron = 0;
}

G4EllipticalCone() [2/3]

G4EllipticalCone::G4EllipticalCone

(

__void__ &

a

)

Definition at line 110 of file G4EllipticalCone.cc.

  : G4VSolid(a), fRebuildPolyhedron(false), fpPolyhedron(0),
    kRadTolerance(0.), halfRadTol(0.), halfCarTol(0.), fCubicVolume(0.),
    fSurfaceArea(0.), xSemiAxis(0.), ySemiAxis(0.), zheight(0.),
    semiAxisMax(0.), zTopCut(0.)
{
}

G4EllipticalCone() [3/3]

G4EllipticalCone::G4EllipticalCone

(

const G4EllipticalCone &

rhs

)

Definition at line 131 of file G4EllipticalCone.cc.

  : G4VSolid(rhs),
    fRebuildPolyhedron(false), fpPolyhedron(0),
    kRadTolerance(rhs.kRadTolerance),
    halfRadTol(rhs.halfRadTol), halfCarTol(rhs.halfCarTol), 
    fCubicVolume(rhs.fCubicVolume), fSurfaceArea(rhs.fSurfaceArea),
    xSemiAxis(rhs.xSemiAxis), ySemiAxis(rhs.ySemiAxis), zheight(rhs.zheight),
    semiAxisMax(rhs.semiAxisMax), zTopCut(rhs.zTopCut)
{
}

Member Function Documentation

CalculateExtent()

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

virtual

Implements G4VSolid.

Definition at line 174 of file G4EllipticalCone.cc.

{
  G4SolidExtentList  extentList( axis, voxelLimit );
  
  //
  // We are going to divide up our elliptical face into small pieces
  //
  
  //
  // Choose phi size of our segment(s) based on constants as
  // defined in meshdefs.hh
  //
  G4int numPhi = kMaxMeshSections;
  G4double sigPhi = twopi/numPhi;
  
  //
  // We have to be careful to keep our segments completely outside
  // of the elliptical surface. To do so we imagine we have
  // a simple (unit radius) circular cross section (as in G4Tubs) 
  // and then "stretch" the dimensions as necessary to fit the ellipse.
  //
  G4double rFudge = 1.0/std::cos(0.5*sigPhi);
  G4double dxFudgeBot = xSemiAxis*2.*zheight*rFudge,
           dyFudgeBot = ySemiAxis*2.*zheight*rFudge;
  G4double dxFudgeTop = xSemiAxis*(zheight-zTopCut)*rFudge,
           dyFudgeTop = ySemiAxis*(zheight-zTopCut)*rFudge;
  
  //
  // As we work around the elliptical surface, we build
  // a "phi" segment on the way, and keep track of two
  // additional polygons for the two ends.
  //
  G4ClippablePolygon endPoly1, endPoly2, phiPoly;
  
  G4double phi = 0, 
           cosPhi = std::cos(phi),
           sinPhi = std::sin(phi);
  G4ThreeVector v0( dxFudgeTop*cosPhi, dyFudgeTop*sinPhi, +zTopCut ),
                v1( dxFudgeBot*cosPhi, dyFudgeBot*sinPhi, -zTopCut ),
                w0, w1;
  transform.ApplyPointTransform( v0 );
  transform.ApplyPointTransform( v1 );
  do    // Loop checking, 13.08.2015, G.Cosmo
  {
    phi += sigPhi;
    if (numPhi == 1) phi = 0;  // Try to avoid roundoff
    cosPhi = std::cos(phi), 
    sinPhi = std::sin(phi);
    
    w0 = G4ThreeVector( dxFudgeTop*cosPhi, dyFudgeTop*sinPhi, +zTopCut );
    w1 = G4ThreeVector( dxFudgeBot*cosPhi, dyFudgeBot*sinPhi, -zTopCut );
    transform.ApplyPointTransform( w0 );
    transform.ApplyPointTransform( w1 );
    
    //
    // Add a point to our z ends
    //
    endPoly1.AddVertexInOrder( v0 );
    endPoly2.AddVertexInOrder( v1 );
    
    //
    // Build phi polygon
    //
    phiPoly.ClearAllVertices();
    
    phiPoly.AddVertexInOrder( v0 );
    phiPoly.AddVertexInOrder( v1 );
    phiPoly.AddVertexInOrder( w1 );
    phiPoly.AddVertexInOrder( w0 );
    
    if (phiPoly.PartialClip( voxelLimit, axis ))
    {
      //
      // Get unit normal
      //
      phiPoly.SetNormal( (v1-v0).cross(w0-v0).unit() );
      
      extentList.AddSurface( phiPoly );
    }

    //
    // Next vertex
    //    
    v0 = w0;
    v1 = w1;
  } while( --numPhi > 0 );

  //
  // Process the end pieces
  //
  if (endPoly1.PartialClip( voxelLimit, axis ))
  {
    static const G4ThreeVector normal(0,0,+1);
    endPoly1.SetNormal( transform.TransformAxis(normal) );
    extentList.AddSurface( endPoly1 );
  }
  
  if (endPoly2.PartialClip( voxelLimit, axis ))
  {
    static const G4ThreeVector normal(0,0,-1);
    endPoly2.SetNormal( transform.TransformAxis(normal) );
    extentList.AddSurface( endPoly2 );
  }
  
  //
  // Return min/max value
  //
  return extentList.GetExtent( min, max );
}

Here is the call graph for this function:

Clone()

G4VSolid * G4EllipticalCone::Clone ( ) const

virtual

Reimplemented from G4VSolid.

Definition at line 967 of file G4EllipticalCone.cc.

{
  return new G4EllipticalCone(*this);
}

Here is the call graph for this function:

CreatePolyhedron()

G4Polyhedron * G4EllipticalCone::CreatePolyhedron ( ) const

virtual

Reimplemented from G4VSolid.

Definition at line 1079 of file G4EllipticalCone.cc.

{
  return new G4PolyhedronEllipticalCone(xSemiAxis, ySemiAxis, zheight, zTopCut);
}

Here is the caller graph for this function:

CreateRotatedVertices()

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

protected

DescribeYourselfTo()

void G4EllipticalCone::DescribeYourselfTo ( G4VGraphicsScene &  scene ) const

virtual

Implements G4VSolid.

Definition at line 1061 of file G4EllipticalCone.cc.

{
  scene.AddSolid(*this);
}

Here is the call graph for this function:

DistanceToIn() [1/2]

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

virtual

Implements G4VSolid.

Definition at line 439 of file G4EllipticalCone.cc.

{
  G4double distMin = kInfinity;

  // code from EllipticalTube

  G4double sigz = p.z()+zTopCut;

  //
  // Check z = -dz planer surface
  //

  if (sigz < halfCarTol)
  {
    //
    // We are "behind" the shape in z, and so can
    // potentially hit the rear face. Correct direction?
    //
    if (v.z() <= 0)
    {
      //
      // As long as we are far enough away, we know we
      // can't intersect
      //
      if (sigz < 0) return kInfinity;
      
      //
      // Otherwise, we don't intersect unless we are
      // on the surface of the ellipse
      //

      if ( sqr(p.x()/( xSemiAxis - halfCarTol ))
         + sqr(p.y()/( ySemiAxis - halfCarTol )) <= sqr( zheight+zTopCut ) )
        return kInfinity;

    }
    else
    {
      //
      // How far?
      //
      G4double q = -sigz/v.z();
      
      //
      // Where does that place us?
      //
      G4double xi = p.x() + q*v.x(),
               yi = p.y() + q*v.y();
      
      //
      // Is this on the surface (within ellipse)?
      //
      if ( sqr(xi/xSemiAxis) + sqr(yi/ySemiAxis) <= sqr( zheight + zTopCut ) )
      {
        //
        // Yup. Return q, unless we are on the surface
        //
        return (sigz < -halfCarTol) ? q : 0;
      }
      else if (xi/(xSemiAxis*xSemiAxis)*v.x()
             + yi/(ySemiAxis*ySemiAxis)*v.y() >= 0)
      {
        //
        // Else, if we are traveling outwards, we know
        // we must miss
        //
        //        return kInfinity;
      }
    }
  }

  //
  // Check z = +dz planer surface
  //
  sigz = p.z() - zTopCut;
  
  if (sigz > -halfCarTol)
  {
    if (v.z() >= 0)
    {

      if (sigz > 0) return kInfinity;

      if ( sqr(p.x()/( xSemiAxis - halfCarTol ))
         + sqr(p.y()/( ySemiAxis - halfCarTol )) <= sqr( zheight-zTopCut ) )
        return kInfinity;

    }
    else {
      G4double q = -sigz/v.z();

      G4double xi = p.x() + q*v.x(),
               yi = p.y() + q*v.y();

      if ( sqr(xi/xSemiAxis) + sqr(yi/ySemiAxis) <= sqr( zheight - zTopCut ) )
      {
        return (sigz > -halfCarTol) ? q : 0;
      }
      else if (xi/(xSemiAxis*xSemiAxis)*v.x()
             + yi/(ySemiAxis*ySemiAxis)*v.y() >= 0)
      {
        //        return kInfinity;
      }
    }
  }


#if 0

  // check to see if Z plane is relevant
  //
  if (p.z() < -zTopCut - 0.5*kCarTolerance)
  {
    if (v.z() <= 0.0)
      return distMin; 

    G4double lambda = (-zTopCut - p.z())/v.z();
    
    if ( sqr((lambda*v.x()+p.x())/xSemiAxis) + 
         sqr((lambda*v.y()+p.y())/ySemiAxis) <=
         sqr(zTopCut + zheight + 0.5*kRadTolerance) ) 
    { 
      return distMin = std::fabs(lambda);    
    }
  }

  if (p.z() > zTopCut+0.5*kCarTolerance) 
  {
    if (v.z() >= 0.0)
      { return distMin; }

    G4double lambda  = (zTopCut - p.z()) / v.z();

    if ( sqr((lambda*v.x() + p.x())/xSemiAxis) + 
         sqr((lambda*v.y() + p.y())/ySemiAxis) <=
         sqr(zheight - zTopCut + 0.5*kRadTolerance) )
      {
        return distMin = std::fabs(lambda);
      }
  }
  
  if (p.z() > zTopCut - halfCarTol
   && p.z() < zTopCut + halfCarTol )
  {
    if (v.z() > 0.) 
      { return kInfinity; }

    return distMin = 0.;
  }
  
  if (p.z() < -zTopCut + halfCarTol
   && p.z() > -zTopCut - halfCarTol)
  {
    if (v.z() < 0.)
      { return distMin = kInfinity; }
    
    return distMin = 0.;
  }
  
#endif

  // if we are here then it either intersects or grazes the curved surface 
  // or it does not intersect at all
  //
  G4double A = sqr(v.x()/xSemiAxis) + sqr(v.y()/ySemiAxis) - sqr(v.z());
  G4double B = 2*(v.x()*p.x()/sqr(xSemiAxis) + 
                  v.y()*p.y()/sqr(ySemiAxis) + v.z()*(zheight-p.z()));
  G4double C = sqr(p.x()/xSemiAxis) + sqr(p.y()/ySemiAxis) - 
               sqr(zheight - p.z());
 
  G4double discr = B*B - 4.*A*C;
   
  // if the discriminant is negative it never hits the curved object
  //
  if ( discr < -halfCarTol )
    { return distMin; }
  
  // case below is when it hits or grazes the surface
  //
  if ( (discr >= - halfCarTol ) && (discr < halfCarTol ) )
  {
    return distMin = std::fabs(-B/(2.*A)); 
  }
  
  G4double plus  = (-B+std::sqrt(discr))/(2.*A);
  G4double minus = (-B-std::sqrt(discr))/(2.*A);
 
  // Special case::Point on Surface, Check norm.dot(v)

  if ( ( std::fabs(plus) < halfCarTol )||( std::fabs(minus) < halfCarTol ) )
  {
    G4ThreeVector truenorm(p.x()/(xSemiAxis*xSemiAxis),
                           p.y()/(ySemiAxis*ySemiAxis),
                           -( p.z() - zheight ));
    if ( truenorm*v >= 0)  //  going outside the solid from surface
    {
      return kInfinity;
    }
    else
    {
      return 0;
    }
  }

  // G4double lambda = std::fabs(plus) < std::fabs(minus) ? plus : minus;  
  G4double lambda = 0;

  if ( minus > halfCarTol && minus < distMin ) 
  {
    lambda = minus ;
    // check normal vector   n * v < 0
    G4ThreeVector pin = p + lambda*v;
    if(std::fabs(pin.z())<zTopCut+0.5*kCarTolerance)
    {
      G4ThreeVector truenorm(pin.x()/(xSemiAxis*xSemiAxis),
                             pin.y()/(ySemiAxis*ySemiAxis),
                             - ( pin.z() - zheight ));
      if ( truenorm*v < 0)
      {   // yes, going inside the solid
        distMin = lambda;
      }
    }
  }
  if ( plus > halfCarTol  && plus < distMin )
  {
    lambda = plus ;
    // check normal vector   n * v < 0
    G4ThreeVector pin = p + lambda*v;
    if(std::fabs(pin.z())<zTopCut+0.5*kCarTolerance)
    {
      G4ThreeVector truenorm(pin.x()/(xSemiAxis*xSemiAxis),
                             pin.y()/(ySemiAxis*ySemiAxis),
                             - ( pin.z() - zheight ) );
      if ( truenorm*v < 0)
      {   // yes, going inside the solid
        distMin = lambda;
      }
    }
  }
  if (distMin < halfCarTol) distMin=0.;
  return distMin ;
}

Here is the call graph for this function:

DistanceToIn() [2/2]

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

virtual

Implements G4VSolid.

Definition at line 688 of file G4EllipticalCone.cc.

{
  G4double distR, distR2, distZ, maxDim;
  G4double distRad;  

  // check if the point lies either below z=-zTopCut in bottom elliptical
  // region or on top within cut elliptical region
  //
  if( (p.z() <= -zTopCut) && (sqr(p.x()/xSemiAxis) + sqr(p.y()/ySemiAxis)
                           <= sqr(zTopCut + zheight + 0.5*kCarTolerance )) )
  {  
    //return distZ = std::fabs(zTopCut - p.z());
     return distZ = std::fabs(zTopCut + p.z());
  } 
  
  if( (p.z() >= zTopCut) && (sqr(p.x()/xSemiAxis)+sqr(p.y()/ySemiAxis)
                          <= sqr(zheight - zTopCut + kCarTolerance/2.0 )) )
  {
    return distZ = std::fabs(p.z() - zTopCut);
  } 
  
  // below we use the following approximation: we take the largest of the
  // axes and find the shortest distance to the circular (cut) cone of that
  // radius.  
  //
  maxDim = xSemiAxis >= ySemiAxis ? xSemiAxis:ySemiAxis;
  distRad = std::sqrt(p.x()*p.x()+p.y()*p.y());

  if( p.z() > maxDim*distRad + zTopCut*(1.+maxDim)-sqr(maxDim)*zheight )
  {
    distR2 = sqr(p.z() - zTopCut) + sqr(distRad - maxDim*(zheight - zTopCut));
    return std::sqrt( distR2 );
  } 

  if( distRad > maxDim*( zheight - p.z() ) )
  {
    if( p.z() > maxDim*distRad - (zTopCut*(1.+maxDim)+sqr(maxDim)*zheight) )
    {
      G4double zVal = (p.z()-maxDim*(distRad-maxDim*zheight))/(1.+sqr(maxDim));
      G4double rVal = maxDim*(zheight - zVal);
      return distR  = std::sqrt(sqr(p.z() - zVal) + sqr(distRad - rVal));
    }
  }

  if( distRad <= maxDim*(zheight - p.z()) )
  {
    distR2 = sqr(distRad - maxDim*(zheight + zTopCut)) + sqr(p.z() + zTopCut);
    return std::sqrt( distR2 );    
  }   
  
  return distR = 0;
}

Here is the call graph for this function:

DistanceToOut() [1/2]

G4double G4EllipticalCone::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 746 of file G4EllipticalCone.cc.

{
  G4double distMin, lambda;
  enum surface_e {kPlaneSurf, kCurvedSurf, kNoSurf} surface;
  
  distMin = kInfinity;
  surface = kNoSurf;

  if (v.z() < 0.0)
  {
    lambda = (-p.z() - zTopCut)/v.z();

    if ( (sqr((p.x() + lambda*v.x())/xSemiAxis) + 
          sqr((p.y() + lambda*v.y())/ySemiAxis)) < 
          sqr(zheight + zTopCut + 0.5*kCarTolerance) )
    {
      distMin = std::fabs(lambda);

      if (!calcNorm) { return distMin; }
    } 
    distMin = std::fabs(lambda);
    surface = kPlaneSurf;
  }

  if (v.z() > 0.0)
  {
    lambda = (zTopCut - p.z()) / v.z();

    if ( (sqr((p.x() + lambda*v.x())/xSemiAxis)
        + sqr((p.y() + lambda*v.y())/ySemiAxis) )
       < (sqr(zheight - zTopCut + 0.5*kCarTolerance)) )
    {
      distMin = std::fabs(lambda);
      if (!calcNorm) { return distMin; }
    }
    distMin = std::fabs(lambda);
    surface = kPlaneSurf;
  }
  
  // if we are here then it either intersects or grazes the 
  // curved surface...
  //
  G4double A = sqr(v.x()/xSemiAxis) + sqr(v.y()/ySemiAxis) - sqr(v.z());
  G4double B = 2.*(v.x()*p.x()/sqr(xSemiAxis) +  
                   v.y()*p.y()/sqr(ySemiAxis) + v.z()*(zheight-p.z()));
  G4double C = sqr(p.x()/xSemiAxis) + sqr(p.y()/ySemiAxis)
             - sqr(zheight - p.z());
 
  G4double discr = B*B - 4.*A*C;
  
  if ( discr >= - 0.5*kCarTolerance && discr < 0.5*kCarTolerance )
  { 
    if(!calcNorm) { return distMin = std::fabs(-B/(2.*A)); }
  }

  else if ( discr > 0.5*kCarTolerance )
  {
    G4double plus  = (-B+std::sqrt(discr))/(2.*A);
    G4double minus = (-B-std::sqrt(discr))/(2.*A);

    if ( plus > 0.5*kCarTolerance && minus > 0.5*kCarTolerance )
    {
      // take the shorter distance
      //
      lambda   = std::fabs(plus) < std::fabs(minus) ? plus : minus;
    }
    else
    {
      // at least one solution is close to zero or negative
      // so, take small positive solution or zero 
      //
      lambda   = plus > -0.5*kCarTolerance ? plus : 0;
    }

    if ( std::fabs(lambda) < distMin )
    {
      if( std::fabs(lambda) > 0.5*kCarTolerance)
      {
        distMin  = std::fabs(lambda);
        surface  = kCurvedSurf;
      }
      else  // Point is On the Surface, Check Normal
      {
        G4ThreeVector truenorm(p.x()/(xSemiAxis*xSemiAxis),
                               p.y()/(ySemiAxis*ySemiAxis),
                               -( p.z() - zheight ));
        if( truenorm.dot(v) > 0 )
        {
          distMin  = 0.0;
          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() - zheight ) );
          truenorm /= 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
                  << "   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("G4EllipticalCone::DistanceToOut(p,v,..)",
                      "GeomSolids1002", JustWarning, message);
          break;
      }
    }
  }

  if (distMin<0.5*kCarTolerance) { distMin=0; }

  return distMin;
}

Here is the call graph for this function:

DistanceToOut() [2/2]

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

virtual

Implements G4VSolid.

Definition at line 908 of file G4EllipticalCone.cc.

{
  G4double rds,roo,roo1, distR, distZ, distMin=0.;
  G4double minAxis = xSemiAxis < ySemiAxis ? xSemiAxis : ySemiAxis;

#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
    
  // since we have made the above warning, below we are working assuming p
  // is inside check how close it is to the circular cone with radius equal
  // to the smaller of the axes
  //
  if( sqr(p.x()/minAxis)+sqr(p.y()/minAxis) < sqr(zheight - p.z()) )
  {
    rds     = std::sqrt(sqr(p.x()) + sqr(p.y()));
    roo     = minAxis*(zheight-p.z()); // radius of cone at z= p.z()
    roo1    = minAxis*(zheight-zTopCut); // radius of cone at z=+zTopCut

    distZ=zTopCut - std::fabs(p.z()) ;
    distR=(roo-rds)/(std::sqrt(1+sqr(minAxis)));

    if(rds>roo1)
    {
      distMin=(zTopCut-p.z())*(roo-rds)/(roo-roo1);
      distMin=std::min(distMin,distR);
    }      
    distMin=std::min(distR,distZ);
  }

  return distMin;
}

Here is the call graph for this function:

GetCubicVolume()

G4double G4EllipticalCone::GetCubicVolume ( )

inlinevirtual

Reimplemented from G4VSolid.

GetEntityType()

G4GeometryType G4EllipticalCone::GetEntityType ( ) const

virtual

Implements G4VSolid.

Definition at line 958 of file G4EllipticalCone.cc.

{
  return G4String("G4EllipticalCone");
}

GetExtent()

G4VisExtent G4EllipticalCone::GetExtent ( ) const

virtual

Reimplemented from G4VSolid.

Definition at line 1066 of file G4EllipticalCone.cc.

{
  // Define the sides of the box into which the solid instance would fit.
  //
  G4double maxDim;
  maxDim = xSemiAxis > ySemiAxis ? xSemiAxis : ySemiAxis;
  maxDim = maxDim > zTopCut ? maxDim : zTopCut;
  
  return G4VisExtent (-maxDim, maxDim,
                      -maxDim, maxDim,
                      -maxDim, maxDim);
}

GetPointOnSurface()

G4ThreeVector G4EllipticalCone::GetPointOnSurface ( ) const

virtual

Reimplemented from G4VSolid.

Definition at line 1001 of file G4EllipticalCone.cc.

{

  G4double phi, sinphi, cosphi, aOne, aTwo, aThree,
           chose, zRand, rRand1, rRand2;
  
  G4double rOne = std::sqrt(sqr(xSemiAxis)
                + sqr(ySemiAxis))*(zheight - zTopCut);
  G4double rTwo = std::sqrt(sqr(xSemiAxis)
                + sqr(ySemiAxis))*(zheight + zTopCut);

  G4int it1=0, it2=0;
  
  aOne   = pi*(rOne + rTwo)*std::sqrt(sqr(rOne - rTwo)+sqr(2.*zTopCut));
  aTwo   = pi*xSemiAxis*ySemiAxis*sqr(zheight+zTopCut);
  aThree = pi*xSemiAxis*ySemiAxis*sqr(zheight-zTopCut);  

  phi = G4RandFlat::shoot(0.,twopi);
  cosphi = std::cos(phi);
  sinphi = std::sin(phi);
  
  if(zTopCut >= zheight) aThree = 0.;

  chose = G4RandFlat::shoot(0.,aOne+aTwo+aThree);
  if((chose>=0.) && (chose<aOne))
  {
    zRand = G4RandFlat::shoot(-zTopCut,zTopCut);
    return G4ThreeVector(xSemiAxis*(zheight-zRand)*cosphi,
                         ySemiAxis*(zheight-zRand)*sinphi,zRand);    
  }
  else if((chose>=aOne) && (chose<aOne+aTwo))
  {
    do    // Loop checking, 13.08.2015, G.Cosmo
    {
      rRand1 = G4RandFlat::shoot(0.,1.) ;
      rRand2 = G4RandFlat::shoot(0.,1.) ;
    } while (( rRand2 >= rRand1  ) && (++it1 < 1000)) ;

    //    rRand2 = G4RandFlat::shoot(0.,std::sqrt(1.-sqr(rRand1)));
    return G4ThreeVector(rRand1*xSemiAxis*(zheight+zTopCut)*cosphi,
                         rRand1*ySemiAxis*(zheight+zTopCut)*sinphi, -zTopCut);

  }
  // else
  //

  do    // Loop checking, 13.08.2015, G.Cosmo
  {
    rRand1 = G4RandFlat::shoot(0.,1.) ;
    rRand2 = G4RandFlat::shoot(0.,1.) ;
  } while (( rRand2 >= rRand1  ) && (++it2 < 1000));

  return G4ThreeVector(rRand1*xSemiAxis*(zheight-zTopCut)*cosphi,
                       rRand1*ySemiAxis*(zheight-zTopCut)*sinphi, zTopCut);
}

Here is the call graph for this function:

GetPolyhedron()

G4Polyhedron * G4EllipticalCone::GetPolyhedron ( ) const

virtual

Reimplemented from G4VSolid.

Definition at line 1084 of file G4EllipticalCone.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 G4EllipticalCone::GetSemiAxisMax ( ) const

inline

Here is the caller graph for this function:

GetSemiAxisX()

G4double G4EllipticalCone::GetSemiAxisX ( ) const

inline

Here is the caller graph for this function:

GetSemiAxisY()

G4double G4EllipticalCone::GetSemiAxisY ( ) const

inline

Here is the caller graph for this function:

GetSurfaceArea()

G4double G4EllipticalCone::GetSurfaceArea ( )

inlinevirtual

Reimplemented from G4VSolid.

GetZMax()

G4double G4EllipticalCone::GetZMax ( ) const

inline

Here is the caller graph for this function:

GetZTopCut()

G4double G4EllipticalCone::GetZTopCut ( ) const

inline

Here is the caller graph for this function:

Inside()

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

virtual

Implements G4VSolid.

Definition at line 293 of file G4EllipticalCone.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() < -zTopCut - halfCarTol)
    || (p.z() > zTopCut + halfCarTol ) )
  {
    return in = kOutside; 
  }

  rad2oo= sqr(p.x()/( xSemiAxis + halfRadTol ))
        + sqr(p.y()/( ySemiAxis + halfRadTol ));

  if ( rad2oo > sqr( zheight-p.z() ) )
  {
    return in = kOutside; 
  }

  //  rad2oi= sqr( p.x()*(1.0 + 0.5*kRadTolerance/(xSemiAxis*xSemiAxis)) )
  //      + sqr( p.y()*(1.0 + 0.5*kRadTolerance/(ySemiAxis*ySemiAxis)) );
  rad2oi = sqr(p.x()/( xSemiAxis - halfRadTol ))
        + sqr(p.y()/( ySemiAxis - halfRadTol ));
     
  if (rad2oi < sqr( zheight-p.z() ) )
  {
    in = ( ( p.z() < -zTopCut + halfRadTol )
        || ( p.z() >  zTopCut - halfRadTol ) ) ? kSurface : kInside;
  }
  else 
  {
    in = kSurface;
  }

  return in;
}

Here is the call graph for this function:

Here is the caller graph for this function:

operator=()

G4EllipticalCone & G4EllipticalCone::operator=

(

const G4EllipticalCone &

rhs

)

Definition at line 146 of file G4EllipticalCone.cc.

{
   // Check assignment to self
   //
   if (this == &rhs)  { return *this; }

   // Copy base class data
   //
   G4VSolid::operator=(rhs);

   // Copy data
   //
   kRadTolerance = rhs.kRadTolerance;
   halfRadTol = rhs.halfRadTol; halfCarTol = rhs.halfCarTol;
   fCubicVolume = rhs.fCubicVolume; fSurfaceArea = rhs.fSurfaceArea;
   xSemiAxis = rhs.xSemiAxis; ySemiAxis = rhs.ySemiAxis;
   zheight = rhs.zheight; semiAxisMax = rhs.semiAxisMax; zTopCut = rhs.zTopCut;
   fRebuildPolyhedron = false;
   delete fpPolyhedron; fpPolyhedron = 0;

   return *this;
}

Here is the call graph for this function:

SetSemiAxis()

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

inline

Here is the caller graph for this function:

SetZCut()

void G4EllipticalCone::SetZCut ( G4double  newzTopCut )

inline

Here is the caller graph for this function:

StreamInfo()

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

virtual

Implements G4VSolid.

Definition at line 976 of file G4EllipticalCone.cc.

{
  G4int oldprc = os.precision(16);
  os << "-----------------------------------------------------------\n"
     << "    *** Dump for solid - " << GetName() << " ***\n"
     << "    ===================================================\n"
     << " Solid type: G4EllipticalCone\n"
     << " Parameters: \n"

     << "    semi-axis x: " << xSemiAxis/mm << " mm \n"
     << "    semi-axis y: " << ySemiAxis/mm << " mm \n"
     << "    height    z: " << zheight/mm << " mm \n"
     << "    half length in  z: " << zTopCut/mm << " mm \n"
     << "-----------------------------------------------------------\n";
  os.precision(oldprc);

  return os;
}

Here is the call graph for this function:

SurfaceNormal()

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

virtual

Implements G4VSolid.

Definition at line 339 of file G4EllipticalCone.cc.

{

  G4double rx = sqr(p.x()/xSemiAxis), 
           ry = sqr(p.y()/ySemiAxis);

  G4double rds = std::sqrt(rx + ry); 

  G4ThreeVector norm;

  if( (p.z() < -zTopCut) && ((rx+ry) < sqr(zTopCut + zheight)) )
  {
    return G4ThreeVector( 0., 0., -1. ); 
  }

  if( (p.z() > (zheight > zTopCut ? zheight : zTopCut)) &&
      ((rx+ry) < sqr(zheight-zTopCut)) )
  {
    return G4ThreeVector( 0., 0., 1. );
  }

  if( p.z() > rds + 2.*zTopCut - zheight ) 
  {
    if ( p.z() > zTopCut )
    {
      if( p.x() == 0. ) 
      {
        norm = G4ThreeVector( 0., p.y() < 0. ? -1. : 1., 1. ); 
        return norm /= norm.mag();
      } 
      if( p.y() == 0. )
      {
        norm = G4ThreeVector( p.x() < 0. ? -1. : 1., 0., 1. ); 
        return norm /= norm.mag();
      } 
      
      G4double k =  std::fabs(p.x()/p.y());
      G4double c2 = sqr(zheight-zTopCut)/(1./sqr(xSemiAxis)+sqr(k/ySemiAxis));
      G4double x  = std::sqrt(c2);
      G4double y  = k*x;
        
      x /= sqr(xSemiAxis);
      y /= sqr(ySemiAxis);
      
      norm = G4ThreeVector( p.x() < 0. ? -x : x, 
                            p.y() < 0. ? -y : y,
                            - ( zheight - zTopCut ) );
      norm /= norm.mag();
      norm += G4ThreeVector( 0., 0., 1. );
      return norm /= norm.mag();      
    }
    
    return G4ThreeVector( 0., 0., 1. );    
  }
  
  if( p.z() < rds - 2.*zTopCut - zheight )
  {
    if( p.x() == 0. ) 
    {
      norm = G4ThreeVector( 0., p.y() < 0. ? -1. : 1., -1. ); 
      return norm /= norm.mag();
    } 
    if( p.y() == 0. )
    {
      norm = G4ThreeVector( p.x() < 0. ? -1. : 1., 0., -1. ); 
      return norm /= norm.mag();
    } 
    
    G4double k =  std::fabs(p.x()/p.y());
    G4double c2 = sqr(zheight+zTopCut)/(1./sqr(xSemiAxis)+sqr(k/ySemiAxis));
    G4double x  = std::sqrt(c2);
    G4double y  = k*x;
    
    x /= sqr(xSemiAxis);
    y /= sqr(ySemiAxis);
    
    norm = G4ThreeVector( p.x() < 0. ? -x : x, 
                          p.y() < 0. ? -y : y,
                          - ( zheight - zTopCut ) );
    norm /= norm.mag();
    norm += G4ThreeVector( 0., 0., -1. );
    return norm /= norm.mag();      
  }
    
  norm  = G4ThreeVector(p.x()/sqr(xSemiAxis), p.y()/sqr(ySemiAxis), rds);
   
  G4double k = std::tan(pi/8.);
  G4double c = -zTopCut - k*(zTopCut + zheight);

  if( p.z() < -k*rds + c )
    return G4ThreeVector (0.,0.,-1.);

  return norm /= norm.mag();
}

Here is the call graph for this function:

Member Data Documentation

fCubicVolume

G4double G4EllipticalCone::fCubicVolume

private

Definition at line 173 of file G4EllipticalCone.hh.

fpPolyhedron

G4Polyhedron* G4EllipticalCone::fpPolyhedron

mutableprotected

Definition at line 166 of file G4EllipticalCone.hh.

fRebuildPolyhedron

G4bool G4EllipticalCone::fRebuildPolyhedron

mutableprotected

Definition at line 165 of file G4EllipticalCone.hh.

fSurfaceArea

G4double G4EllipticalCone::fSurfaceArea

private

Definition at line 174 of file G4EllipticalCone.hh.

halfCarTol

G4double G4EllipticalCone::halfCarTol

private

Definition at line 171 of file G4EllipticalCone.hh.

halfRadTol

G4double G4EllipticalCone::halfRadTol

private

Definition at line 171 of file G4EllipticalCone.hh.

kRadTolerance

G4double G4EllipticalCone::kRadTolerance

private

Definition at line 170 of file G4EllipticalCone.hh.

semiAxisMax

G4double G4EllipticalCone::semiAxisMax

private

Definition at line 175 of file G4EllipticalCone.hh.

xSemiAxis

G4double G4EllipticalCone::xSemiAxis

private

Definition at line 175 of file G4EllipticalCone.hh.

ySemiAxis

G4double G4EllipticalCone::ySemiAxis

private

Definition at line 175 of file G4EllipticalCone.hh.

zheight

G4double G4EllipticalCone::zheight

private

Definition at line 175 of file G4EllipticalCone.hh.

zTopCut

G4double G4EllipticalCone::zTopCut

private

Definition at line 175 of file G4EllipticalCone.hh.

---

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

• G4EllipticalCone.hh
• G4EllipticalCone.cc