G4Cons Class Reference

#include <G4Cons.hh>

Inheritance diagram for G4Cons:

Collaboration diagram for G4Cons:

Public Member Functions
	G4Cons (const G4String &pName, G4double pRmin1, G4double pRmax1, G4double pRmin2, G4double pRmax2, G4double pDz, G4double pSPhi, G4double pDPhi)
	~G4Cons ()
G4double	GetInnerRadiusMinusZ () const
G4double	GetOuterRadiusMinusZ () const
G4double	GetInnerRadiusPlusZ () const
G4double	GetOuterRadiusPlusZ () const
G4double	GetZHalfLength () const
G4double	GetStartPhiAngle () const
G4double	GetDeltaPhiAngle () const
void	SetInnerRadiusMinusZ (G4double Rmin1)
void	SetOuterRadiusMinusZ (G4double Rmax1)
void	SetInnerRadiusPlusZ (G4double Rmin2)
void	SetOuterRadiusPlusZ (G4double Rmax2)
void	SetZHalfLength (G4double newDz)
void	SetStartPhiAngle (G4double newSPhi, G4bool trig=true)
void	SetDeltaPhiAngle (G4double newDPhi)
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
G4ThreeVector	GetPointOnSurface () const
G4VSolid *	Clone () const
std::ostream &	StreamInfo (std::ostream &os) const
void	DescribeYourselfTo (G4VGraphicsScene &scene) const
G4Polyhedron *	CreatePolyhedron () const
	G4Cons (__void__ &)
	G4Cons (const G4Cons &rhs)
G4Cons &	operator= (const G4Cons &rhs)
G4double	GetRmin1 () const
G4double	GetRmax1 () const
G4double	GetRmin2 () const
G4double	GetRmax2 () const
G4double	GetDz () const
G4double	GetSPhi () const
G4double	GetDPhi () const
Public Member Functions inherited from G4CSGSolid
	G4CSGSolid (const G4String &pName)
virtual	~G4CSGSolid ()
virtual G4Polyhedron *	GetPolyhedron () const
	G4CSGSolid (__void__ &)
	G4CSGSolid (const G4CSGSolid &rhs)
G4CSGSolid &	operator= (const G4CSGSolid &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 G4VisExtent	GetExtent () 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

Private Types
enum	ESide {   kNull, kRMin, kRMax, kSPhi,   kEPhi, kPZ, kMZ }
enum	ENorm {   kNRMin, kNRMax, kNSPhi, kNEPhi,   kNZ }

Private Member Functions
G4ThreeVectorList *	CreateRotatedVertices (const G4AffineTransform &pTransform) const
void	Initialize ()
void	CheckSPhiAngle (G4double sPhi)
void	CheckDPhiAngle (G4double dPhi)
void	CheckPhiAngles (G4double sPhi, G4double dPhi)
void	InitializeTrigonometry ()
G4ThreeVector	ApproxSurfaceNormal (const G4ThreeVector &p) const

Private Attributes
G4double	kRadTolerance
G4double	kAngTolerance
G4double	fRmin1
G4double	fRmin2
G4double	fRmax1
G4double	fRmax2
G4double	fDz
G4double	fSPhi
G4double	fDPhi
G4double	sinCPhi
G4double	cosCPhi
G4double	cosHDPhiOT
G4double	cosHDPhiIT
G4double	sinSPhi
G4double	cosSPhi
G4double	sinEPhi
G4double	cosEPhi
G4bool	fPhiFullCone
G4double	halfCarTolerance
G4double	halfRadTolerance
G4double	halfAngTolerance

Additional Inherited Members
Protected Member Functions inherited from G4CSGSolid
G4double	GetRadiusInRing (G4double rmin, G4double rmax) 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 inherited from G4CSGSolid
G4double	fCubicVolume
G4double	fSurfaceArea
G4bool	fRebuildPolyhedron
G4Polyhedron *	fpPolyhedron
Protected Attributes inherited from G4VSolid
G4double	kCarTolerance

Detailed Description

Definition at line 83 of file G4Cons.hh.

Member Enumeration Documentation

ENorm

enum G4Cons::ENorm

private

Enumerator
kNRMin
kNRMax
kNSPhi
kNEPhi
kNZ

Definition at line 214 of file G4Cons.hh.

{kNRMin,kNRMax,kNSPhi,kNEPhi,kNZ};

ESide

enum G4Cons::ESide

private

Enumerator
kNull
kRMin
kRMax
kSPhi
kEPhi
kPZ
kMZ

Definition at line 210 of file G4Cons.hh.

{kNull,kRMin,kRMax,kSPhi,kEPhi,kPZ,kMZ};

Constructor & Destructor Documentation

G4Cons() [1/3]

G4Cons::G4Cons	(	const G4String &	pName,
		G4double	pRmin1,
		G4double	pRmax1,
		G4double	pRmin2,
		G4double	pRmax2,
		G4double	pDz,
		G4double	pSPhi,
		G4double	pDPhi
	)

Definition at line 82 of file G4Cons.cc.

  : G4CSGSolid(pName), fRmin1(pRmin1), fRmin2(pRmin2),
    fRmax1(pRmax1), fRmax2(pRmax2), fDz(pDz), fSPhi(0.), fDPhi(0.)
{
  kRadTolerance = G4GeometryTolerance::GetInstance()->GetRadialTolerance();
  kAngTolerance = G4GeometryTolerance::GetInstance()->GetAngularTolerance();

  halfCarTolerance=kCarTolerance*0.5;
  halfRadTolerance=kRadTolerance*0.5;
  halfAngTolerance=kAngTolerance*0.5;

  // Check z-len
  //
  if ( pDz < 0 )
  {
    std::ostringstream message;
    message << "Invalid Z half-length for Solid: " << GetName() << G4endl
            << "        hZ = " << pDz;
    G4Exception("G4Cons::G4Cons()", "GeomSolids0002",
                FatalException, message);
  }

  // Check radii
  //
  if (((pRmin1>=pRmax1) || (pRmin2>=pRmax2) || (pRmin1<0)) && (pRmin2<0))
  {
    std::ostringstream message;
    message << "Invalid values of radii for Solid: " << GetName() << G4endl
            << "        pRmin1 = " << pRmin1 << ", pRmin2 = " << pRmin2
            << ", pRmax1 = " << pRmax1 << ", pRmax2 = " << pRmax2;
    G4Exception("G4Cons::G4Cons()", "GeomSolids0002",
                FatalException, message) ;
  }
  if( (pRmin1 == 0.0) && (pRmin2 > 0.0) ) { fRmin1 = 1e3*kRadTolerance ; }
  if( (pRmin2 == 0.0) && (pRmin1 > 0.0) ) { fRmin2 = 1e3*kRadTolerance ; }

  // Check angles
  //
  CheckPhiAngles(pSPhi, pDPhi);
}

Here is the call graph for this function:

Here is the caller graph for this function:

~G4Cons()

G4Cons::~G4Cons

(

)

Definition at line 146 of file G4Cons.cc.

{
}

G4Cons() [2/3]

G4Cons::G4Cons

(

__void__ &

a

)

Definition at line 132 of file G4Cons.cc.

  : G4CSGSolid(a), kRadTolerance(0.), kAngTolerance(0.),
    fRmin1(0.), fRmin2(0.), fRmax1(0.), fRmax2(0.), fDz(0.),
    fSPhi(0.), fDPhi(0.), sinCPhi(0.), cosCPhi(0.), cosHDPhiOT(0.),
    cosHDPhiIT(0.), sinSPhi(0.), cosSPhi(0.), sinEPhi(0.), cosEPhi(0.),
    fPhiFullCone(false), halfCarTolerance(0.), halfRadTolerance(0.),
    halfAngTolerance(0.)
{
}

G4Cons() [3/3]

G4Cons::G4Cons

(

const G4Cons &

rhs

)

Definition at line 154 of file G4Cons.cc.

  : G4CSGSolid(rhs), kRadTolerance(rhs.kRadTolerance),
    kAngTolerance(rhs.kAngTolerance), fRmin1(rhs.fRmin1), fRmin2(rhs.fRmin2),
    fRmax1(rhs.fRmax1), fRmax2(rhs.fRmax2), fDz(rhs.fDz), fSPhi(rhs.fSPhi),
    fDPhi(rhs.fDPhi), sinCPhi(rhs.sinCPhi), cosCPhi(rhs.cosCPhi),
    cosHDPhiOT(rhs.cosHDPhiOT), cosHDPhiIT(rhs.cosHDPhiIT),
    sinSPhi(rhs.sinSPhi), cosSPhi(rhs.cosSPhi), sinEPhi(rhs.sinEPhi),
    cosEPhi(rhs.cosEPhi), fPhiFullCone(rhs.fPhiFullCone),
    halfCarTolerance(rhs.halfCarTolerance),
    halfRadTolerance(rhs.halfRadTolerance),
    halfAngTolerance(rhs.halfAngTolerance)
{
}

Member Function Documentation

ApproxSurfaceNormal()

G4ThreeVector G4Cons::ApproxSurfaceNormal ( const G4ThreeVector &  p ) const

private

Definition at line 594 of file G4Cons.cc.

{
  ENorm side ;
  G4ThreeVector norm ;
  G4double rho, phi ;
  G4double distZ, distRMin, distRMax, distSPhi, distEPhi, distMin ;
  G4double tanRMin, secRMin, pRMin, widRMin ;
  G4double tanRMax, secRMax, pRMax, widRMax ;

  distZ = std::fabs(std::fabs(p.z()) - fDz) ;
  rho   = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;

  tanRMin  = (fRmin2 - fRmin1)*0.5/fDz ;
  secRMin  = std::sqrt(1 + tanRMin*tanRMin) ;
  pRMin    = rho - p.z()*tanRMin ;
  widRMin  = fRmin2 - fDz*tanRMin ;
  distRMin = std::fabs(pRMin - widRMin)/secRMin ;

  tanRMax  = (fRmax2 - fRmax1)*0.5/fDz ;
  secRMax  = std::sqrt(1+tanRMax*tanRMax) ;
  pRMax    = rho - p.z()*tanRMax ;
  widRMax  = fRmax2 - fDz*tanRMax ;
  distRMax = std::fabs(pRMax - widRMax)/secRMax ;
  
  if (distRMin < distRMax)  // First minimum
  {
    if (distZ < distRMin)
    {
      distMin = distZ ;
      side    = kNZ ;
    }
    else
    {
      distMin = distRMin ;
      side    = kNRMin ;
    }
  }
  else
  {
    if (distZ < distRMax)
    {
      distMin = distZ ;
      side    = kNZ ;
    }
    else
    {
      distMin = distRMax ;
      side    = kNRMax ;
    }
  }
  if ( !fPhiFullCone && rho )  // Protected against (0,0,z) 
  {
    phi = std::atan2(p.y(),p.x()) ;

    if (phi < 0)  { phi += twopi; }

    if (fSPhi < 0)  { distSPhi = std::fabs(phi - (fSPhi + twopi))*rho; }
    else            { distSPhi = std::fabs(phi - fSPhi)*rho; }

    distEPhi = std::fabs(phi - fSPhi - fDPhi)*rho ;

    // Find new minimum

    if (distSPhi < distEPhi)
    {
      if (distSPhi < distMin)  { side = kNSPhi; }
    }
    else 
    {
      if (distEPhi < distMin)  { side = kNEPhi; }
    }
  }    
  switch (side)
  {
    case kNRMin:      // Inner radius
      rho *= secRMin ;
      norm = G4ThreeVector(-p.x()/rho, -p.y()/rho, tanRMin/secRMin) ;
      break ;
    case kNRMax:      // Outer radius
      rho *= secRMax ;
      norm = G4ThreeVector(p.x()/rho, p.y()/rho, -tanRMax/secRMax) ;
      break ;
    case kNZ:         // +/- dz
      if (p.z() > 0)  { norm = G4ThreeVector(0,0,1);  }
      else            { norm = G4ThreeVector(0,0,-1); }
      break ;
    case kNSPhi:
      norm = G4ThreeVector(std::sin(fSPhi), -std::cos(fSPhi), 0) ;
      break ;
    case kNEPhi:
      norm=G4ThreeVector(-std::sin(fSPhi+fDPhi), std::cos(fSPhi+fDPhi), 0) ;
      break ;
    default:          // Should never reach this case...
      DumpInfo();
      G4Exception("G4Cons::ApproxSurfaceNormal()",
                  "GeomSolids1002", JustWarning,
                  "Undefined side for valid surface normal to solid.");
      break ;    
  }
  return norm ;
}

Here is the call graph for this function:

Here is the caller graph for this function:

CalculateExtent()

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

virtual

Implements G4VSolid.

Definition at line 272 of file G4Cons.cc.

{
  if ( !pTransform.IsRotated() && (fDPhi == twopi)
    && (fRmin1 == 0) && (fRmin2 == 0) )
  {
    // Special case handling for unrotated solid cones
    // Compute z/x/y 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 diff1, diff2, delta, maxDiff, newMin, newMax, RMax ;
    G4double xoff1, xoff2, yoff1, yoff2 ;
      
    zoffset = pTransform.NetTranslation().z();
    zMin    = zoffset - fDz ;
    zMax    = zoffset + fDz ;

    if (pVoxelLimit.IsZLimited())
    {
      if( (zMin > pVoxelLimit.GetMaxZExtent() + kCarTolerance) || 
          (zMax < pVoxelLimit.GetMinZExtent() - kCarTolerance)  )
      {
        return false ;
      }
      else
      {
        if ( zMin < pVoxelLimit.GetMinZExtent() )
        {
          zMin = pVoxelLimit.GetMinZExtent() ;
        }
        if ( zMax > pVoxelLimit.GetMaxZExtent() )
        {
          zMax = pVoxelLimit.GetMaxZExtent() ;
        }
      }
    }
    xoffset = pTransform.NetTranslation().x() ;
    RMax    = (fRmax2 >= fRmax1) ?  zMax : zMin  ;                  
    xMax    = xoffset + (fRmax1 + fRmax2)*0.5 + 
              (RMax - zoffset)*(fRmax2 - fRmax1)/(2*fDz) ;
    xMin    = 2*xoffset-xMax ;

    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() ;
    yMax    = yoffset + (fRmax1 + fRmax2)*0.5 + 
              (RMax - zoffset)*(fRmax2 - fRmax1)/(2*fDz) ;
    yMin    = 2*yoffset-yMax ;
    RMax    = yMax - yoffset ;  // = max radius due to Zmax/Zmin cuttings

    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() ;
        }
      }
    }    
    switch (pAxis) // Known to cut cones
    {
      case kXAxis:
        yoff1 = yoffset - yMin ;
        yoff2 = yMax - yoffset ;

        if ((yoff1 >= 0) && (yoff2 >= 0)) // Y limits cross max/min x
        {                                 // => no change
          pMin = xMin ;
          pMax = xMax ;
        }
        else
        {
          // Y limits don't cross max/min x => compute max delta x,
          // hence new mins/maxs
          delta=RMax*RMax-yoff1*yoff1;
          diff1=(delta>0.) ? std::sqrt(delta) : 0.;
          delta=RMax*RMax-yoff2*yoff2;
          diff2=(delta>0.) ? std::sqrt(delta) : 0.;
          maxDiff = (diff1>diff2) ? diff1:diff2 ;
          newMin  = xoffset - maxDiff ;
          newMax  = xoffset + maxDiff ;
          pMin    = ( newMin < xMin ) ? xMin : newMin  ;
          pMax    = ( newMax > xMax) ? xMax : newMax ;
        } 
      break ;

      case kYAxis:
        xoff1 = xoffset - xMin ;
        xoff2 = xMax - xoffset ;

        if ((xoff1 >= 0) && (xoff2 >= 0) ) // X limits cross max/min y
        {                                  // => no change
          pMin = yMin ;
          pMax = yMax ;
        }
        else
        {
          // X limits don't cross max/min y => compute max delta y,
          // hence new mins/maxs
          delta=RMax*RMax-xoff1*xoff1;
          diff1=(delta>0.) ? std::sqrt(delta) : 0.;
          delta=RMax*RMax-xoff2*xoff2;
          diff2=(delta>0.) ? std::sqrt(delta) : 0.;
          maxDiff = (diff1 > diff2) ? diff1:diff2 ;
          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   // Calculate rotated vertex coordinates
  {
    G4int i, noEntries, noBetweenSections4 ;
    G4bool existsAfterClip = false ;
    G4ThreeVectorList* vertices = CreateRotatedVertices(pTransform) ;

    pMin = +kInfinity ;
    pMax = -kInfinity ;

    noEntries          = vertices->size() ;
    noBetweenSections4 = noEntries-4 ;
      
    for ( i = 0 ; i < noEntries ; i += 4 )
    {
      ClipCrossSection(vertices, i, pVoxelLimit, pAxis, pMin, pMax) ;
    }
    for ( i = 0 ; i < noBetweenSections4 ; i += 4 )
    {
      ClipBetweenSections(vertices, i, pVoxelLimit, pAxis, pMin, pMax) ;
    }    
    if ( (pMin != kInfinity) || (pMax != -kInfinity) )
    {
      existsAfterClip = true ;
        
      // Add 2*tolerance to avoid precision troubles

      pMin -= kCarTolerance ;
      pMax += kCarTolerance ;
    }
    else
    {
      // Check for case where completely enveloping clipping volume
      // If point inside then we are confident that the solid completely
      // envelopes the clipping volume. Hence set min/max extents according
      // to clipping volume extents along the specified axis.
       
      G4ThreeVector clipCentre(
      (pVoxelLimit.GetMinXExtent() + pVoxelLimit.GetMaxXExtent())*0.5,
      (pVoxelLimit.GetMinYExtent() + pVoxelLimit.GetMaxYExtent())*0.5,
      (pVoxelLimit.GetMinZExtent() + pVoxelLimit.GetMaxZExtent())*0.5  ) ;
        
      if (Inside(pTransform.Inverse().TransformPoint(clipCentre)) != kOutside)
      {
        existsAfterClip = true ;
        pMin            = pVoxelLimit.GetMinExtent(pAxis) ;
        pMax            = pVoxelLimit.GetMaxExtent(pAxis) ;
      }
    }
    delete vertices ;
    return existsAfterClip ;
  }
}

Here is the call graph for this function:

CheckDPhiAngle()

void G4Cons::CheckDPhiAngle ( G4double  dPhi )

inlineprivate

CheckPhiAngles()

void G4Cons::CheckPhiAngles ( G4double  sPhi, G4double  dPhi  )

inlineprivate

Here is the caller graph for this function:

CheckSPhiAngle()

void G4Cons::CheckSPhiAngle ( G4double  sPhi )

inlineprivate

Clone()

G4VSolid * G4Cons::Clone ( ) const

virtual

Reimplemented from G4VSolid.

Definition at line 2256 of file G4Cons.cc.

{
  return new G4Cons(*this);
}

Here is the call graph for this function:

ComputeDimensions()

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

virtual

Reimplemented from G4VSolid.

Definition at line 260 of file G4Cons.cc.

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

Here is the call graph for this function:

CreatePolyhedron()

G4Polyhedron * G4Cons::CreatePolyhedron ( ) const

virtual

Reimplemented from G4VSolid.

Definition at line 2384 of file G4Cons.cc.

{
  return new G4PolyhedronCons(fRmin1,fRmax1,fRmin2,fRmax2,fDz,fSPhi,fDPhi);
}

CreateRotatedVertices()

G4ThreeVectorList * G4Cons::CreateRotatedVertices ( const G4AffineTransform &  pTransform ) const

private

Definition at line 2161 of file G4Cons.cc.

{
  G4ThreeVectorList* vertices ;
  G4ThreeVector vertex0, vertex1, vertex2, vertex3 ;
  G4double meshAngle, meshRMax1, meshRMax2, crossAngle;
  G4double cosCrossAngle, sinCrossAngle, sAngle ;
  G4double rMaxX1, rMaxX2, rMaxY1, rMaxY2, rMinX1, rMinX2, rMinY1, rMinY2 ;
  G4int crossSection, noCrossSections ;

  // Compute no of cross-sections necessary to mesh cone
    
  noCrossSections = G4int(fDPhi/kMeshAngleDefault) + 1 ;

  if (noCrossSections < kMinMeshSections)
  {
    noCrossSections = kMinMeshSections ;
  }
  else if (noCrossSections > kMaxMeshSections)
  {
    noCrossSections = kMaxMeshSections ;
  }
  meshAngle = fDPhi/(noCrossSections - 1) ;

  meshRMax1 = fRmax1/std::cos(meshAngle*0.5) ;
  meshRMax2 = fRmax2/std::cos(meshAngle*0.5) ;

  // If complete in phi, set start angle such that mesh will be at RMax
  // on the x axis. Will give better extent calculations when not rotated.

  if ( fPhiFullCone && (fSPhi == 0.0) )
  {
    sAngle = -meshAngle*0.5 ;
  }
  else
  {
    sAngle = fSPhi ;
  } 
  vertices = new G4ThreeVectorList();

  if (vertices)
  {
    vertices->reserve(noCrossSections*4) ;
    for (crossSection = 0 ; crossSection < noCrossSections ; crossSection++)
    {
      // Compute coordinates of cross section at section crossSection

      crossAngle    = sAngle + crossSection*meshAngle ;
      cosCrossAngle = std::cos(crossAngle) ;
      sinCrossAngle = std::sin(crossAngle) ;

      rMaxX1 = meshRMax1*cosCrossAngle ;
      rMaxY1 = meshRMax1*sinCrossAngle ;
      rMaxX2 = meshRMax2*cosCrossAngle ;
      rMaxY2 = meshRMax2*sinCrossAngle ;
        
      rMinX1 = fRmin1*cosCrossAngle ;
      rMinY1 = fRmin1*sinCrossAngle ;
      rMinX2 = fRmin2*cosCrossAngle ;
      rMinY2 = fRmin2*sinCrossAngle ;
        
      vertex0 = G4ThreeVector(rMinX1,rMinY1,-fDz) ;
      vertex1 = G4ThreeVector(rMaxX1,rMaxY1,-fDz) ;
      vertex2 = G4ThreeVector(rMaxX2,rMaxY2,+fDz) ;
      vertex3 = G4ThreeVector(rMinX2,rMinY2,+fDz) ;

      vertices->push_back(pTransform.TransformPoint(vertex0)) ;
      vertices->push_back(pTransform.TransformPoint(vertex1)) ;
      vertices->push_back(pTransform.TransformPoint(vertex2)) ;
      vertices->push_back(pTransform.TransformPoint(vertex3)) ;
    }
  }
  else
  {
    DumpInfo();
    G4Exception("G4Cons::CreateRotatedVertices()",
                "GeomSolids0003", FatalException,
                "Error in allocation of vertices. Out of memory !");
  }

  return vertices ;
}

Here is the call graph for this function:

Here is the caller graph for this function:

DescribeYourselfTo()

void G4Cons::DescribeYourselfTo ( G4VGraphicsScene &  scene ) const

virtual

Implements G4VSolid.

Definition at line 2379 of file G4Cons.cc.

{
  scene.AddSolid (*this);
}

Here is the call graph for this function:

DistanceToIn() [1/2]

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

virtual

Implements G4VSolid.

Definition at line 720 of file G4Cons.cc.

{
  G4double snxt = kInfinity ;      // snxt = default return value
  const G4double dRmax = 50*(fRmax1+fRmax2);// 100*(Rmax1+Rmax2)/2.

  G4double tanRMax,secRMax,rMaxAv,rMaxOAv ;  // Data for cones
  G4double tanRMin,secRMin,rMinAv,rMinOAv ;
  G4double rout,rin ;

  G4double tolORMin,tolORMin2,tolIRMin,tolIRMin2 ; // `generous' radii squared
  G4double tolORMax2,tolIRMax,tolIRMax2 ;
  G4double tolODz,tolIDz ;

  G4double Dist,sd,xi,yi,zi,ri=0.,risec,rhoi2,cosPsi ; // Intersection point vars

  G4double t1,t2,t3,b,c,d ;    // Quadratic solver variables 
  G4double nt1,nt2,nt3 ;
  G4double Comp ;

  G4ThreeVector Normal;

  // Cone Precalcs

  tanRMin = (fRmin2 - fRmin1)*0.5/fDz ;
  secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ;
  rMinAv  = (fRmin1 + fRmin2)*0.5 ;

  if (rMinAv > halfRadTolerance)
  {
    rMinOAv = rMinAv - halfRadTolerance ;
  }
  else
  {
    rMinOAv = 0.0 ;
  }  
  tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
  secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
  rMaxAv  = (fRmax1 + fRmax2)*0.5 ;
  rMaxOAv = rMaxAv + halfRadTolerance ;
   
  // Intersection with z-surfaces

  tolIDz = fDz - halfCarTolerance ;
  tolODz = fDz + halfCarTolerance ;

  if (std::fabs(p.z()) >= tolIDz)
  {
    if ( p.z()*v.z() < 0 )    // at +Z going in -Z or visa versa
    {
      sd = (std::fabs(p.z()) - fDz)/std::fabs(v.z()) ; // Z intersect distance

      if( sd < 0.0 )  { sd = 0.0; }                    // negative dist -> zero

      xi   = p.x() + sd*v.x() ;  // Intersection coords
      yi   = p.y() + sd*v.y() ;
      rhoi2 = xi*xi + yi*yi  ;

      // Check validity of intersection
      // Calculate (outer) tolerant radi^2 at intersecion

      if (v.z() > 0)
      {
        tolORMin  = fRmin1 - halfRadTolerance*secRMin ;
        tolIRMin  = fRmin1 + halfRadTolerance*secRMin ;
        tolIRMax  = fRmax1 - halfRadTolerance*secRMin ;
        // tolORMax2 = (fRmax1 + halfRadTolerance*secRMax)*
        //             (fRmax1 + halfRadTolerance*secRMax) ;
      }
      else
      {
        tolORMin  = fRmin2 - halfRadTolerance*secRMin ;
        tolIRMin  = fRmin2 + halfRadTolerance*secRMin ;
        tolIRMax  = fRmax2 - halfRadTolerance*secRMin ;
        // tolORMax2 = (fRmax2 + halfRadTolerance*secRMax)*
        //             (fRmax2 + halfRadTolerance*secRMax) ;
      }
      if ( tolORMin > 0 ) 
      {
        // tolORMin2 = tolORMin*tolORMin ;
        tolIRMin2 = tolIRMin*tolIRMin ;
      }
      else                
      {
        // tolORMin2 = 0.0 ;
        tolIRMin2 = 0.0 ;
      }
      if ( tolIRMax > 0 )  { tolIRMax2 = tolIRMax*tolIRMax; }     
      else                 { tolIRMax2 = 0.0; }
      
      if ( (tolIRMin2 <= rhoi2) && (rhoi2 <= tolIRMax2) )
      {
        if ( !fPhiFullCone && rhoi2 )
        {
          // Psi = angle made with central (average) phi of shape

          cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;

          if (cosPsi >= cosHDPhiIT)  { return sd; }
        }
        else
        {
          return sd;
        }
      }
    }
    else  // On/outside extent, and heading away  -> cannot intersect
    {
      return snxt ;  
    }
  }
    
// ----> Can not intersect z surfaces


// Intersection with outer cone (possible return) and
//                   inner cone (must also check phi)
//
// Intersection point (xi,yi,zi) on line x=p.x+t*v.x etc.
//
// Intersects with x^2+y^2=(a*z+b)^2
//
// where a=tanRMax or tanRMin
//       b=rMaxAv  or rMinAv
//
// (vx^2+vy^2-(a*vz)^2)t^2+2t(pxvx+pyvy-a*vz(a*pz+b))+px^2+py^2-(a*pz+b)^2=0 ;
//     t1                        t2                      t3  
//
//  \--------u-------/       \-----------v----------/ \---------w--------/
//

  t1   = 1.0 - v.z()*v.z() ;
  t2   = p.x()*v.x() + p.y()*v.y() ;
  t3   = p.x()*p.x() + p.y()*p.y() ;
  rin  = tanRMin*p.z() + rMinAv ;
  rout = tanRMax*p.z() + rMaxAv ;

  // Outer Cone Intersection
  // Must be outside/on outer cone for valid intersection

  nt1 = t1 - (tanRMax*v.z())*(tanRMax*v.z()) ;
  nt2 = t2 - tanRMax*v.z()*rout ;
  nt3 = t3 - rout*rout ;

  if (std::fabs(nt1) > kRadTolerance)  // Equation quadratic => 2 roots
  {
    b = nt2/nt1;
    c = nt3/nt1;
    d = b*b-c  ;
    if ( (nt3 > rout*rout*kRadTolerance*kRadTolerance*secRMax*secRMax)
      || (rout < 0) )
    {
      // If outside real cone (should be rho-rout>kRadTolerance*0.5
      // NOT rho^2 etc) saves a std::sqrt() at expense of accuracy

      if (d >= 0)
      {
          
        if ((rout < 0) && (nt3 <= 0))
        {
          // Inside `shadow cone' with -ve radius
          // -> 2nd root could be on real cone

          if (b>0) { sd = c/(-b-std::sqrt(d)); }
          else     { sd = -b + std::sqrt(d);   }
        }
        else
        {
          if ((b <= 0) && (c >= 0)) // both >=0, try smaller root
          {
            sd=c/(-b+std::sqrt(d));
          }
          else
          {
            if ( c <= 0 ) // second >=0
            {
              sd = -b + std::sqrt(d) ;
              if((sd<0) & (sd>-halfRadTolerance)) sd=0;
            }
            else  // both negative, travel away
            {
              return kInfinity ;
            }
          }
        }
        if ( sd >= 0 )  // If 'forwards'. Check z intersection
        {
          if ( sd>dRmax ) // Avoid rounding errors due to precision issues on
          {               // 64 bits systems. Split long distances and recompute
            G4double fTerm = sd-std::fmod(sd,dRmax);
            sd = fTerm + DistanceToIn(p+fTerm*v,v);
          } 
          zi = p.z() + sd*v.z() ;

          if (std::fabs(zi) <= tolODz)
          {
            // Z ok. Check phi intersection if reqd

            if ( fPhiFullCone )  { return sd; }
            else
            {
              xi     = p.x() + sd*v.x() ;
              yi     = p.y() + sd*v.y() ;
              ri     = rMaxAv + zi*tanRMax ;
              cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;

              if ( cosPsi >= cosHDPhiIT )  { return sd; }
            }
          }
        }                // end if (sd>0)
      }
    }
    else
    {
      // Inside outer cone
      // check not inside, and heading through G4Cons (-> 0 to in)

      if ( ( t3  > (rin + halfRadTolerance*secRMin)*
                   (rin + halfRadTolerance*secRMin) )
        && (nt2 < 0) && (d >= 0) && (std::fabs(p.z()) <= tolIDz) )
      {
        // Inside cones, delta r -ve, inside z extent
        // Point is on the Surface => check Direction using  Normal.dot(v)

        xi     = p.x() ;
        yi     = p.y()  ;
        risec  = std::sqrt(xi*xi + yi*yi)*secRMax ;
        Normal = G4ThreeVector(xi/risec,yi/risec,-tanRMax/secRMax) ;
        if ( !fPhiFullCone )
        {
          cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(t3) ;
          if ( cosPsi >= cosHDPhiIT )
          {
            if ( Normal.dot(v) <= 0 )  { return 0.0; }
          }
        }
        else
        {             
          if ( Normal.dot(v) <= 0 )  { return 0.0; }
        }
      }
    }
  }
  else  //  Single root case 
  {
    if ( std::fabs(nt2) > kRadTolerance )
    {
      sd = -0.5*nt3/nt2 ;

      if ( sd < 0 )  { return kInfinity; }   // travel away
      else  // sd >= 0,  If 'forwards'. Check z intersection
      {
        zi = p.z() + sd*v.z() ;

        if ((std::fabs(zi) <= tolODz) && (nt2 < 0))
        {
          // Z ok. Check phi intersection if reqd

          if ( fPhiFullCone )  { return sd; }
          else
          {
            xi     = p.x() + sd*v.x() ;
            yi     = p.y() + sd*v.y() ;
            ri     = rMaxAv + zi*tanRMax ;
            cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;

            if (cosPsi >= cosHDPhiIT)  { return sd; }
          }
        }
      }
    }
    else  //    travel || cone surface from its origin
    {
      sd = kInfinity ;
    }
  }

  // Inner Cone Intersection
  // o Space is divided into 3 areas:
  //   1) Radius greater than real inner cone & imaginary cone & outside
  //      tolerance
  //   2) Radius less than inner or imaginary cone & outside tolarance
  //   3) Within tolerance of real or imaginary cones
  //      - Extra checks needed for 3's intersections
  //        => lots of duplicated code

  if (rMinAv)
  { 
    nt1 = t1 - (tanRMin*v.z())*(tanRMin*v.z()) ;
    nt2 = t2 - tanRMin*v.z()*rin ;
    nt3 = t3 - rin*rin ;
 
    if ( nt1 )
    {
      if ( nt3 > rin*kRadTolerance*secRMin )
      {
        // At radius greater than real & imaginary cones
        // -> 2nd root, with zi check

        b = nt2/nt1 ;
        c = nt3/nt1 ;
        d = b*b-c ;
        if (d >= 0)   // > 0
        {
           if(b>0){sd = c/( -b-std::sqrt(d));}
           else   {sd = -b + std::sqrt(d) ;}

          if ( sd >= 0 )   // > 0
          {
            if ( sd>dRmax ) // Avoid rounding errors due to precision issues on
            {               // 64 bits systems. Split long distance and recompute
              G4double fTerm = sd-std::fmod(sd,dRmax);
              sd = fTerm + DistanceToIn(p+fTerm*v,v);
            } 
            zi = p.z() + sd*v.z() ;

            if ( std::fabs(zi) <= tolODz )
            {
              if ( !fPhiFullCone )
              {
                xi     = p.x() + sd*v.x() ;
                yi     = p.y() + sd*v.y() ;
                ri     = rMinAv + zi*tanRMin ;
                cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;

                if (cosPsi >= cosHDPhiIT)
                { 
                  if ( sd > halfRadTolerance )  { snxt=sd; }
                  else
                  {
                    // Calculate a normal vector in order to check Direction

                    risec  = std::sqrt(xi*xi + yi*yi)*secRMin ;
                    Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin);
                    if ( Normal.dot(v) <= 0 )  { snxt = sd; }
                  } 
                }
              }
              else
              {
                if ( sd > halfRadTolerance )  { return sd; }
                else
                {
                  // Calculate a normal vector in order to check Direction

                  xi     = p.x() + sd*v.x() ;
                  yi     = p.y() + sd*v.y() ;
                  risec  = std::sqrt(xi*xi + yi*yi)*secRMin ;
                  Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ;
                  if ( Normal.dot(v) <= 0 )  { return sd; }
                }
              }
            }
          }
        }
      }
      else  if ( nt3 < -rin*kRadTolerance*secRMin )
      {
        // Within radius of inner cone (real or imaginary)
        // -> Try 2nd root, with checking intersection is with real cone
        // -> If check fails, try 1st root, also checking intersection is
        //    on real cone

        b = nt2/nt1 ;
        c = nt3/nt1 ;
        d = b*b - c ;

        if ( d >= 0 )  // > 0
        {
          if (b>0) { sd = c/(-b-std::sqrt(d)); }
          else     { sd = -b + std::sqrt(d);   }
          zi = p.z() + sd*v.z() ;
          ri = rMinAv + zi*tanRMin ;

          if ( ri > 0 )
          {
            if ( (sd >= 0) && (std::fabs(zi) <= tolODz) )  // sd > 0
            {
              if ( sd>dRmax ) // Avoid rounding errors due to precision issues
              {               // seen on 64 bits systems. Split and recompute
                G4double fTerm = sd-std::fmod(sd,dRmax);
                sd = fTerm + DistanceToIn(p+fTerm*v,v);
              } 
              if ( !fPhiFullCone )
              {
                xi     = p.x() + sd*v.x() ;
                yi     = p.y() + sd*v.y() ;
                cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;

                if (cosPsi >= cosHDPhiOT)
                {
                  if ( sd > halfRadTolerance )  { snxt=sd; }
                  else
                  {
                    // Calculate a normal vector in order to check Direction

                    risec  = std::sqrt(xi*xi + yi*yi)*secRMin ;
                    Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin);
                    if ( Normal.dot(v) <= 0 )  { snxt = sd; } 
                  }
                }
              }
              else
              {
                if( sd > halfRadTolerance )  { return sd; }
                else
                {
                  // Calculate a normal vector in order to check Direction

                  xi     = p.x() + sd*v.x() ;
                  yi     = p.y() + sd*v.y() ;
                  risec  = std::sqrt(xi*xi + yi*yi)*secRMin ;
                  Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ;
                  if ( Normal.dot(v) <= 0 )  { return sd; }
                } 
              }
            }
          }
          else
          {
        if (b>0) { sd = -b - std::sqrt(d);   }
            else     { sd = c/(-b+std::sqrt(d)); }
            zi = p.z() + sd*v.z() ;
            ri = rMinAv + zi*tanRMin ;

            if ( (sd >= 0) && (ri > 0) && (std::fabs(zi) <= tolODz) ) // sd>0
            {
              if ( sd>dRmax ) // Avoid rounding errors due to precision issues
              {               // seen on 64 bits systems. Split and recompute
                G4double fTerm = sd-std::fmod(sd,dRmax);
                sd = fTerm + DistanceToIn(p+fTerm*v,v);
              } 
              if ( !fPhiFullCone )
              {
                xi     = p.x() + sd*v.x() ;
                yi     = p.y() + sd*v.y() ;
                cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;

                if (cosPsi >= cosHDPhiIT)
                {
                  if ( sd > halfRadTolerance )  { snxt=sd; }
                  else
                  {
                    // Calculate a normal vector in order to check Direction

                    risec  = std::sqrt(xi*xi + yi*yi)*secRMin ;
                    Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin);
                    if ( Normal.dot(v) <= 0 )  { snxt = sd; } 
                  }
                }
              }
              else
              {
                if ( sd > halfRadTolerance )  { return sd; }
                else
                {
                  // Calculate a normal vector in order to check Direction

                  xi     = p.x() + sd*v.x() ;
                  yi     = p.y() + sd*v.y() ;
                  risec  = std::sqrt(xi*xi + yi*yi)*secRMin ;
                  Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ;
                  if ( Normal.dot(v) <= 0 )  { return sd; }
                } 
              }
            }
          }
        }
      }
      else
      {
        // Within kRadTol*0.5 of inner cone (real OR imaginary)
        // ----> Check not travelling through (=>0 to in)
        // ----> if not:
        //    -2nd root with validity check

        if ( std::fabs(p.z()) <= tolODz )
        {
          if ( nt2 > 0 )
          {
            // Inside inner real cone, heading outwards, inside z range

            if ( !fPhiFullCone )
            {
              cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(t3) ;

              if (cosPsi >= cosHDPhiIT)  { return 0.0; }
            }
            else  { return 0.0; }
          }
          else
          {
            // Within z extent, but not travelling through
            // -> 2nd root or kInfinity if 1st root on imaginary cone

            b = nt2/nt1 ;
            c = nt3/nt1 ;
            d = b*b - c ;

            if ( d >= 0 )   // > 0
            {
              if (b>0) { sd = -b - std::sqrt(d);   }
              else     { sd = c/(-b+std::sqrt(d)); }
              zi = p.z() + sd*v.z() ;
              ri = rMinAv + zi*tanRMin ;
              
              if ( ri > 0 )   // 2nd root
              {
                if (b>0) { sd = c/(-b-std::sqrt(d)); }
                else     { sd = -b + std::sqrt(d);   }
                
                zi = p.z() + sd*v.z() ;

                if ( (sd >= 0) && (std::fabs(zi) <= tolODz) )  // sd>0
                {
                  if ( sd>dRmax ) // Avoid rounding errors due to precision issue
                  {               // seen on 64 bits systems. Split and recompute
                    G4double fTerm = sd-std::fmod(sd,dRmax);
                    sd = fTerm + DistanceToIn(p+fTerm*v,v);
                  } 
                  if ( !fPhiFullCone )
                  {
                    xi     = p.x() + sd*v.x() ;
                    yi     = p.y() + sd*v.y() ;
                    ri     = rMinAv + zi*tanRMin ;
                    cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ;

                    if ( cosPsi >= cosHDPhiIT )  { snxt = sd; }
                  }
                  else  { return sd; }
                }
              }
              else  { return kInfinity; }
            }
          }
        }
        else   // 2nd root
        {
          b = nt2/nt1 ;
          c = nt3/nt1 ;
          d = b*b - c ;

          if ( d > 0 )
          {  
            if (b>0) { sd = c/(-b-std::sqrt(d)); }
            else     { sd = -b + std::sqrt(d) ;  }
            zi = p.z() + sd*v.z() ;

            if ( (sd >= 0) && (std::fabs(zi) <= tolODz) )  // sd>0
            {
              if ( sd>dRmax ) // Avoid rounding errors due to precision issues
              {               // seen on 64 bits systems. Split and recompute
                G4double fTerm = sd-std::fmod(sd,dRmax);
                sd = fTerm + DistanceToIn(p+fTerm*v,v);
              } 
              if ( !fPhiFullCone )
              {
                xi     = p.x() + sd*v.x();
                yi     = p.y() + sd*v.y();
                ri     = rMinAv + zi*tanRMin ;
                cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri;

                if (cosPsi >= cosHDPhiIT)  { snxt = sd; }
              }
              else  { return sd; }
            }
          }
        }
      }
    }
  }

  // Phi segment intersection
  //
  // o Tolerant of points inside phi planes by up to kCarTolerance*0.5
  //
  // o NOTE: Large duplication of code between sphi & ephi checks
  //         -> only diffs: sphi -> ephi, Comp -> -Comp and half-plane
  //            intersection check <=0 -> >=0
  //         -> Should use some form of loop Construct

  if ( !fPhiFullCone )
  {
    // First phi surface (starting phi)

    Comp    = v.x()*sinSPhi - v.y()*cosSPhi ;
                    
    if ( Comp < 0 )    // Component in outwards normal dirn
    {
      Dist = (p.y()*cosSPhi - p.x()*sinSPhi) ;

      if (Dist < halfCarTolerance)
      {
        sd = Dist/Comp ;

        if ( sd < snxt )
        {
          if ( sd < 0 )  { sd = 0.0; }

          zi = p.z() + sd*v.z() ;

          if ( std::fabs(zi) <= tolODz )
          {
            xi        = p.x() + sd*v.x() ;
            yi        = p.y() + sd*v.y() ;
            rhoi2     = xi*xi + yi*yi ;
            tolORMin2 = (rMinOAv + zi*tanRMin)*(rMinOAv + zi*tanRMin) ;
            tolORMax2 = (rMaxOAv + zi*tanRMax)*(rMaxOAv + zi*tanRMax) ;

            if ( (rhoi2 >= tolORMin2) && (rhoi2 <= tolORMax2) )
            {
              // z and r intersections good - check intersecting with
              // correct half-plane

              if ((yi*cosCPhi - xi*sinCPhi) <= 0 )  { snxt = sd; }
            }
          }
        }
      }
    }

    // Second phi surface (Ending phi)

    Comp    = -(v.x()*sinEPhi - v.y()*cosEPhi) ;
        
    if ( Comp < 0 )   // Component in outwards normal dirn
    {
      Dist = -(p.y()*cosEPhi - p.x()*sinEPhi) ;
      if (Dist < halfCarTolerance)
      {
        sd = Dist/Comp ;

        if ( sd < snxt )
        {
          if ( sd < 0 )  { sd = 0.0; }

          zi = p.z() + sd*v.z() ;

          if (std::fabs(zi) <= tolODz)
          {
            xi        = p.x() + sd*v.x() ;
            yi        = p.y() + sd*v.y() ;
            rhoi2     = xi*xi + yi*yi ;
            tolORMin2 = (rMinOAv + zi*tanRMin)*(rMinOAv + zi*tanRMin) ;
            tolORMax2 = (rMaxOAv + zi*tanRMax)*(rMaxOAv + zi*tanRMax) ;

            if ( (rhoi2 >= tolORMin2) && (rhoi2 <= tolORMax2) )
            {
              // z and r intersections good - check intersecting with
              // correct half-plane

              if ( (yi*cosCPhi - xi*sinCPhi) >= 0.0 )  { snxt = sd; }
            }
          }
        }
      }
    }
  }
  if (snxt < halfCarTolerance)  { snxt = 0.; }

  return snxt ;
}

Here is the call graph for this function:

DistanceToIn() [2/2]

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

virtual

Implements G4VSolid.

Definition at line 1390 of file G4Cons.cc.

{
  G4double safe=0.0, rho, safeR1, safeR2, safeZ, safePhi, cosPsi ;
  G4double tanRMin, secRMin, pRMin ;
  G4double tanRMax, secRMax, pRMax ;

  rho   = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
  safeZ = std::fabs(p.z()) - fDz ;

  if ( fRmin1 || fRmin2 )
  {
    tanRMin = (fRmin2 - fRmin1)*0.5/fDz ;
    secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ;
    pRMin   = tanRMin*p.z() + (fRmin1 + fRmin2)*0.5 ;
    safeR1  = (pRMin - rho)/secRMin ;

    tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
    secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
    pRMax   = tanRMax*p.z() + (fRmax1 + fRmax2)*0.5 ;
    safeR2  = (rho - pRMax)/secRMax ;

    if ( safeR1 > safeR2) { safe = safeR1; }
    else                  { safe = safeR2; }
  }
  else
  {
    tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
    secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
    pRMax   = tanRMax*p.z() + (fRmax1 + fRmax2)*0.5 ;
    safe    = (rho - pRMax)/secRMax ;
  }
  if ( safeZ > safe )  { safe = safeZ; }

  if ( !fPhiFullCone && rho )
  {
    // Psi=angle from central phi to point

    cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/rho ;

    if ( cosPsi < std::cos(fDPhi*0.5) ) // Point lies outside phi range
    {
      if ( (p.y()*cosCPhi - p.x()*sinCPhi) <= 0.0 )
      {
        safePhi = std::fabs(p.x()*std::sin(fSPhi)-p.y()*std::cos(fSPhi));
      }
      else
      {
        safePhi = std::fabs(p.x()*sinEPhi-p.y()*cosEPhi);
      }
      if ( safePhi > safe )  { safe = safePhi; }
    }
  }
  if ( safe < 0.0 )  { safe = 0.0; }

  return safe ;
}

Here is the call graph for this function:

DistanceToOut() [1/2]

G4double G4Cons::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 1452 of file G4Cons.cc.

{
  ESide side = kNull, sider = kNull, sidephi = kNull;

  G4double snxt,srd,sphi,pdist ;

  G4double tanRMax, secRMax, rMaxAv ;  // Data for outer cone
  G4double tanRMin, secRMin, rMinAv ;  // Data for inner cone

  G4double t1, t2, t3, rout, rin, nt1, nt2, nt3 ;
  G4double b, c, d, sr2, sr3 ;

  // Vars for intersection within tolerance

  ESide    sidetol = kNull ;
  G4double slentol = kInfinity ;

  // Vars for phi intersection:

  G4double pDistS, compS, pDistE, compE, sphi2, xi, yi, risec, vphi ;
  G4double zi, ri, deltaRoi2 ;

  // Z plane intersection

  if ( v.z() > 0.0 )
  {
    pdist = fDz - p.z() ;

    if (pdist > halfCarTolerance)
    {
      snxt = pdist/v.z() ;
      side = kPZ ;
    }
    else
    {
      if (calcNorm)
      {
        *n         = G4ThreeVector(0,0,1) ;
        *validNorm = true ;
      }
      return  snxt = 0.0;
    }
  }
  else if ( v.z() < 0.0 )
  {
    pdist = fDz + p.z() ;

    if ( pdist > halfCarTolerance)
    {
      snxt = -pdist/v.z() ;
      side = kMZ ;
    }
    else
    {
      if ( calcNorm )
      {
        *n         = G4ThreeVector(0,0,-1) ;
        *validNorm = true ;
      }
      return snxt = 0.0 ;
    }
  }
  else     // Travel perpendicular to z axis
  {
    snxt = kInfinity ;    
    side = kNull ;
  }

  // Radial Intersections
  //
  // Intersection with outer cone (possible return) and
  //                   inner cone (must also check phi)
  //
  // Intersection point (xi,yi,zi) on line x=p.x+t*v.x etc.
  //
  // Intersects with x^2+y^2=(a*z+b)^2
  //
  // where a=tanRMax or tanRMin
  //       b=rMaxAv  or rMinAv
  //
  // (vx^2+vy^2-(a*vz)^2)t^2+2t(pxvx+pyvy-a*vz(a*pz+b))+px^2+py^2-(a*pz+b)^2=0 ;
  //     t1                        t2                      t3  
  //
  //  \--------u-------/       \-----------v----------/ \---------w--------/

  tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
  secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
  rMaxAv  = (fRmax1 + fRmax2)*0.5 ;


  t1   = 1.0 - v.z()*v.z() ;      // since v normalised
  t2   = p.x()*v.x() + p.y()*v.y() ;
  t3   = p.x()*p.x() + p.y()*p.y() ;
  rout = tanRMax*p.z() + rMaxAv ;

  nt1 = t1 - (tanRMax*v.z())*(tanRMax*v.z()) ;
  nt2 = t2 - tanRMax*v.z()*rout ;
  nt3 = t3 - rout*rout ;

  if (v.z() > 0.0)
  {
    deltaRoi2 = snxt*snxt*t1 + 2*snxt*t2 + t3
                - fRmax2*(fRmax2 + kRadTolerance*secRMax);
  }
  else if ( v.z() < 0.0 )
  {
    deltaRoi2 = snxt*snxt*t1 + 2*snxt*t2 + t3
                - fRmax1*(fRmax1 + kRadTolerance*secRMax);
  }
  else
  {
    deltaRoi2 = 1.0;
  }

  if ( nt1 && (deltaRoi2 > 0.0) )  
  {
    // Equation quadratic => 2 roots : second root must be leaving

    b = nt2/nt1 ;
    c = nt3/nt1 ;
    d = b*b - c ;

    if ( d >= 0 )
    {
      // Check if on outer cone & heading outwards
      // NOTE: Should use rho-rout>-kRadTolerance*0.5
        
      if (nt3 > -halfRadTolerance && nt2 >= 0 )
      {
        if (calcNorm)
        {
          risec      = std::sqrt(t3)*secRMax ;
          *validNorm = true ;
          *n         = G4ThreeVector(p.x()/risec,p.y()/risec,-tanRMax/secRMax);
        }
        return snxt=0 ;
      }
      else
      {
        sider = kRMax  ;
        if (b>0) { srd = -b - std::sqrt(d);    }
        else     { srd = c/(-b+std::sqrt(d)) ; }

        zi    = p.z() + srd*v.z() ;
        ri    = tanRMax*zi + rMaxAv ;
          
        if ((ri >= 0) && (-halfRadTolerance <= srd) && (srd <= halfRadTolerance))
        {
          // An intersection within the tolerance
          //   we will Store it in case it is good -
          // 
          slentol = srd ;
          sidetol = kRMax ;
        }            
        if ( (ri < 0) || (srd < halfRadTolerance) )
        {
          // Safety: if both roots -ve ensure that srd cannot `win'
          //         distance to out

          if (b>0) { sr2 = c/(-b-std::sqrt(d)); }
          else     { sr2 = -b + std::sqrt(d);   }
          zi  = p.z() + sr2*v.z() ;
          ri  = tanRMax*zi + rMaxAv ;

          if ((ri >= 0) && (sr2 > halfRadTolerance))
          {
            srd = sr2;
          }
          else
          {
            srd = kInfinity ;

            if( (-halfRadTolerance <= sr2) && ( sr2 <= halfRadTolerance) )
            {
              // An intersection within the tolerance.
              // Storing it in case it is good.

              slentol = sr2 ;
              sidetol = kRMax ;
            }
          }
        }
      }
    }
    else
    {
      // No intersection with outer cone & not parallel
      // -> already outside, no intersection

      if ( calcNorm )
      {
        risec      = std::sqrt(t3)*secRMax;
        *validNorm = true;
        *n         = G4ThreeVector(p.x()/risec,p.y()/risec,-tanRMax/secRMax);
      }
      return snxt = 0.0 ;
    }
  }
  else if ( nt2 && (deltaRoi2 > 0.0) )
  {
    // Linear case (only one intersection) => point outside outer cone

    if ( calcNorm )
    {
      risec      = std::sqrt(t3)*secRMax;
      *validNorm = true;
      *n         = G4ThreeVector(p.x()/risec,p.y()/risec,-tanRMax/secRMax);
    }
    return snxt = 0.0 ;
  }
  else
  {
    // No intersection -> parallel to outer cone
    // => Z or inner cone intersection

    srd = kInfinity ;
  }

  // Check possible intersection within tolerance

  if ( slentol <= halfCarTolerance )
  {
    // An intersection within the tolerance was found.  
    // We must accept it only if the momentum points outwards.  
    //
    // G4ThreeVector ptTol ;  // The point of the intersection  
    // ptTol= p + slentol*v ;
    // ri=tanRMax*zi+rMaxAv ;
    //
    // Calculate a normal vector,  as below

    xi    = p.x() + slentol*v.x();
    yi    = p.y() + slentol*v.y();
    risec = std::sqrt(xi*xi + yi*yi)*secRMax;
    G4ThreeVector Normal = G4ThreeVector(xi/risec,yi/risec,-tanRMax/secRMax);

    if ( Normal.dot(v) > 0 )    // We will leave the Cone immediatelly
    {
      if ( calcNorm ) 
      {
        *n         = Normal.unit() ;
        *validNorm = true ;
      }
      return snxt = 0.0 ;
    }
    else // On the surface, but not heading out so we ignore this intersection
    {    //                                        (as it is within tolerance).
      slentol = kInfinity ;
    }
  }

  // Inner Cone intersection

  if ( fRmin1 || fRmin2 )
  {
    tanRMin = (fRmin2 - fRmin1)*0.5/fDz ;
    nt1     = t1 - (tanRMin*v.z())*(tanRMin*v.z()) ;

    if ( nt1 )
    {
      secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ;
      rMinAv  = (fRmin1 + fRmin2)*0.5 ;    
      rin     = tanRMin*p.z() + rMinAv ;
      nt2     = t2 - tanRMin*v.z()*rin ;
      nt3     = t3 - rin*rin ;
      
      // Equation quadratic => 2 roots : first root must be leaving

      b = nt2/nt1 ;
      c = nt3/nt1 ;
      d = b*b - c ;

      if ( d >= 0.0 )
      {
        // NOTE: should be rho-rin<kRadTolerance*0.5,
        //       but using squared versions for efficiency

        if (nt3 < kRadTolerance*(rin + kRadTolerance*0.25)) 
        {
          if ( nt2 < 0.0 )
          {
            if (calcNorm)  { *validNorm = false; }
            return          snxt      = 0.0;
          }
        }
        else
        {
          if (b>0) { sr2 = -b - std::sqrt(d);   }
          else     { sr2 = c/(-b+std::sqrt(d)); }
          zi  = p.z() + sr2*v.z() ;
          ri  = tanRMin*zi + rMinAv ;

          if( (ri>=0.0)&&(-halfRadTolerance<=sr2)&&(sr2<=halfRadTolerance) )
          {
            // An intersection within the tolerance
            // storing it in case it is good.

            slentol = sr2 ;
            sidetol = kRMax ;
          }
          if( (ri<0) || (sr2 < halfRadTolerance) )
          {
            if (b>0) { sr3 = c/(-b-std::sqrt(d)); }
            else     { sr3 = -b + std::sqrt(d) ;  }

            // Safety: if both roots -ve ensure that srd cannot `win'
            //         distancetoout

            if  ( sr3 > halfRadTolerance )
            {
              if( sr3 < srd )
              {
                zi = p.z() + sr3*v.z() ;
                ri = tanRMin*zi + rMinAv ;

                if ( ri >= 0.0 )
                {
                  srd=sr3 ;
                  sider=kRMin ;
                }
              } 
            }
            else if ( sr3 > -halfRadTolerance )
            {
              // Intersection in tolerance. Store to check if it's good

              slentol = sr3 ;
              sidetol = kRMin ;
            }
          }
          else if ( (sr2 < srd) && (sr2 > halfCarTolerance) )
          {
            srd   = sr2 ;
            sider = kRMin ;
          }
          else if (sr2 > -halfCarTolerance)
          {
            // Intersection in tolerance. Store to check if it's good

            slentol = sr2 ;
            sidetol = kRMin ;
          }    
          if( slentol <= halfCarTolerance  )
          {
            // An intersection within the tolerance was found. 
            // We must accept it only if  the momentum points outwards. 

            G4ThreeVector Normal ; 
            
            // Calculate a normal vector,  as below

            xi     = p.x() + slentol*v.x() ;
            yi     = p.y() + slentol*v.y() ;
            if( sidetol==kRMax )
            {
              risec  = std::sqrt(xi*xi + yi*yi)*secRMax ;
              Normal = G4ThreeVector(xi/risec,yi/risec,-tanRMax/secRMax) ;
            }
            else
            {
              risec  = std::sqrt(xi*xi + yi*yi)*secRMin ;
              Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ;
            }
            if( Normal.dot(v) > 0 )
            {
              // We will leave the cone immediately

              if( calcNorm ) 
              {
                *n         = Normal.unit() ;
                *validNorm = true ;
              }
              return snxt = 0.0 ;
            }
            else 
            { 
              // On the surface, but not heading out so we ignore this
              // intersection (as it is within tolerance). 

              slentol = kInfinity ;
            }        
          }
        }
      }
    }
  }

  // Linear case => point outside inner cone ---> outer cone intersect
  //
  // Phi Intersection
  
  if ( !fPhiFullCone )
  {
    // add angle calculation with correction 
    // of the difference in domain of atan2 and Sphi

    vphi = std::atan2(v.y(),v.x()) ;

    if ( vphi < fSPhi - halfAngTolerance  )              { vphi += twopi; }
    else if ( vphi > fSPhi + fDPhi + halfAngTolerance )  { vphi -= twopi; }

    if ( p.x() || p.y() )   // Check if on z axis (rho not needed later)
    {
      // pDist -ve when inside

      pDistS = p.x()*sinSPhi - p.y()*cosSPhi ;
      pDistE = -p.x()*sinEPhi + p.y()*cosEPhi ;

      // Comp -ve when in direction of outwards normal

      compS = -sinSPhi*v.x() + cosSPhi*v.y() ;
      compE = sinEPhi*v.x() - cosEPhi*v.y() ;

      sidephi = kNull ;

      if( ( (fDPhi <= pi) && ( (pDistS <= halfCarTolerance)
                            && (pDistE <= halfCarTolerance) ) )
         || ( (fDPhi >  pi) && !((pDistS >  halfCarTolerance)
                              && (pDistE >  halfCarTolerance) ) )  )
      {
        // Inside both phi *full* planes
        if ( compS < 0 )
        {
          sphi = pDistS/compS ;
          if (sphi >= -halfCarTolerance)
          {
            xi = p.x() + sphi*v.x() ;
            yi = p.y() + sphi*v.y() ;

            // Check intersecting with correct half-plane
            // (if not -> no intersect)
            //
            if ( (std::fabs(xi)<=kCarTolerance)
              && (std::fabs(yi)<=kCarTolerance) )
            {
              sidephi= kSPhi;
              if ( ( fSPhi-halfAngTolerance <= vphi )
                && ( fSPhi+fDPhi+halfAngTolerance >=vphi ) )
              {
                sphi = kInfinity;
              }
            }
            else
            if ( (yi*cosCPhi-xi*sinCPhi)>=0 )
            {
              sphi = kInfinity ;
            }
            else
            {
              sidephi = kSPhi ;
              if ( pDistS > -halfCarTolerance )
              {
                sphi = 0.0 ; // Leave by sphi immediately
              }    
            }       
          }
          else
          {
            sphi = kInfinity ;
          }
        }
        else
        {
          sphi = kInfinity ;
        }

        if ( compE < 0 )
        {
          sphi2 = pDistE/compE ;

          // Only check further if < starting phi intersection
          //
          if ( (sphi2 > -halfCarTolerance) && (sphi2 < sphi) )
          {
            xi = p.x() + sphi2*v.x() ;
            yi = p.y() + sphi2*v.y() ;

            // Check intersecting with correct half-plane

            if ( (std::fabs(xi)<=kCarTolerance)
              && (std::fabs(yi)<=kCarTolerance) )
            {
              // Leaving via ending phi

              if(!( (fSPhi-halfAngTolerance <= vphi)
                 && (fSPhi+fDPhi+halfAngTolerance >= vphi) ) )
              {
                sidephi = kEPhi ;
                if ( pDistE <= -halfCarTolerance )  { sphi = sphi2; }
                else                                { sphi = 0.0; }
              }
            }
            else // Check intersecting with correct half-plane
            if ( yi*cosCPhi-xi*sinCPhi >= 0 )
            {
              // Leaving via ending phi

              sidephi = kEPhi ;
              if ( pDistE <= -halfCarTolerance )  { sphi = sphi2; }
              else                                { sphi = 0.0; }
            }
          }
        }
      }
      else
      {
        sphi = kInfinity ;
      }
    }
    else
    {
      // On z axis + travel not || to z axis -> if phi of vector direction
      // within phi of shape, Step limited by rmax, else Step =0

      if ( (fSPhi-halfAngTolerance <= vphi)
        && (vphi <= fSPhi+fDPhi+halfAngTolerance) )
      {
        sphi = kInfinity ;
      }
      else
      {
        sidephi = kSPhi  ;   // arbitrary 
        sphi    = 0.0 ;
      }
    }      
    if ( sphi < snxt )  // Order intersecttions
    {
      snxt=sphi ;
      side=sidephi ;
    }
  }
  if ( srd < snxt )  // Order intersections
  {
    snxt = srd   ;
    side = sider ;
  }
  if (calcNorm)
  {
    switch(side)
    {                     // Note: returned vector not normalised
      case kRMax:         // (divide by frmax for unit vector)
        xi         = p.x() + snxt*v.x() ;
        yi         = p.y() + snxt*v.y() ;
        risec      = std::sqrt(xi*xi + yi*yi)*secRMax ;
        *n         = G4ThreeVector(xi/risec,yi/risec,-tanRMax/secRMax) ;
        *validNorm = true ;
        break ;
      case kRMin:
        *validNorm = false ;  // Rmin is inconvex
        break ;
      case kSPhi:
        if ( fDPhi <= pi )
        {
          *n         = G4ThreeVector(sinSPhi, -cosSPhi, 0);
          *validNorm = true ;
        }
        else
        {
          *validNorm = false ;
        }
        break ;
      case kEPhi:
        if ( fDPhi <= pi )
        {
          *n = G4ThreeVector(-sinEPhi, cosEPhi, 0);
          *validNorm = true ;
        }
        else
        {
          *validNorm = false ;
        }
        break ;
      case kPZ:
        *n         = G4ThreeVector(0,0,1) ;
        *validNorm = true ;
        break ;
      case kMZ:
        *n         = G4ThreeVector(0,0,-1) ;
        *validNorm = true ;
        break ;
      default:
        G4cout << G4endl ;
        DumpInfo();
        std::ostringstream message;
        G4int oldprc = message.precision(16) ;
        message << "Undefined side for valid surface normal to solid."
                << G4endl
                << "Position:"  << G4endl << G4endl
                << "p.x() = "   << p.x()/mm << " mm" << G4endl
                << "p.y() = "   << p.y()/mm << " mm" << G4endl
                << "p.z() = "   << p.z()/mm << " mm" << G4endl << G4endl
                << "pho at z = "   << std::sqrt( p.x()*p.x()+p.y()*p.y() )/mm
                << " mm" << G4endl << G4endl ;
        if( p.x() != 0. || p.y() != 0.)
        {
           message << "point phi = "   << std::atan2(p.y(),p.x())/degree
                   << " degree" << G4endl << G4endl ; 
        }
        message << "Direction:" << G4endl << G4endl
                << "v.x() = "   << v.x() << G4endl
                << "v.y() = "   << v.y() << G4endl
                << "v.z() = "   << v.z() << G4endl<< G4endl
                << "Proposed distance :" << G4endl<< G4endl
                << "snxt = "    << snxt/mm << " mm" << G4endl ;
        message.precision(oldprc) ;
        G4Exception("G4Cons::DistanceToOut(p,v,..)","GeomSolids1002",
                    JustWarning, message) ;
        break ;
    }
  }
  if (snxt < halfCarTolerance)  { snxt = 0.; }

  return snxt ;
}

Here is the call graph for this function:

DistanceToOut() [2/2]

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

virtual

Implements G4VSolid.

Definition at line 2075 of file G4Cons.cc.

{
  G4double safe=0.0, rho, safeR1, safeR2, safeZ, safePhi;
  G4double tanRMin, secRMin, pRMin;
  G4double tanRMax, secRMax, pRMax;

#ifdef G4CSGDEBUG
  if( Inside(p) == kOutside )
  {
    G4int oldprc=G4cout.precision(16) ;
    G4cout << G4endl ;
    DumpInfo();
    G4cout << "Position:"  << G4endl << G4endl ;
    G4cout << "p.x() = "   << p.x()/mm << " mm" << G4endl ;
    G4cout << "p.y() = "   << p.y()/mm << " mm" << G4endl ;
    G4cout << "p.z() = "   << p.z()/mm << " mm" << G4endl << G4endl ;
    G4cout << "pho at z = "   << std::sqrt( p.x()*p.x()+p.y()*p.y() )/mm
           << " mm" << G4endl << G4endl ;
    if( (p.x() != 0.) || (p.x() != 0.) )
    {
      G4cout << "point phi = "   << std::atan2(p.y(),p.x())/degree
             << " degree" << G4endl << G4endl ; 
    }
    G4cout.precision(oldprc) ;
    G4Exception("G4Cons::DistanceToOut(p)", "GeomSolids1002",
                JustWarning, "Point p is outside !?" );
  }
#endif

  rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ;
  safeZ = fDz - std::fabs(p.z()) ;

  if (fRmin1 || fRmin2)
  {
    tanRMin = (fRmin2 - fRmin1)*0.5/fDz ;
    secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ;
    pRMin   = tanRMin*p.z() + (fRmin1 + fRmin2)*0.5 ;
    safeR1  = (rho - pRMin)/secRMin ;
  }
  else
  {
    safeR1 = kInfinity ;
  }

  tanRMax = (fRmax2 - fRmax1)*0.5/fDz ;
  secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ;
  pRMax   = tanRMax*p.z() + (fRmax1+fRmax2)*0.5 ;
  safeR2  = (pRMax - rho)/secRMax ;

  if (safeR1 < safeR2)  { safe = safeR1; }
  else                  { safe = safeR2; }
  if (safeZ < safe)     { safe = safeZ ; }

  // Check if phi divided, Calc distances closest phi plane

  if (!fPhiFullCone)
  {
    // Above/below central phi of G4Cons?

    if ( (p.y()*cosCPhi - p.x()*sinCPhi) <= 0 )
    {
      safePhi = -(p.x()*sinSPhi - p.y()*cosSPhi) ;
    }
    else
    {
      safePhi = (p.x()*sinEPhi - p.y()*cosEPhi) ;
    }
    if (safePhi < safe)  { safe = safePhi; }
  }
  if ( safe < 0 )  { safe = 0; }

  return safe ;
}

Here is the call graph for this function:

GetCubicVolume()

G4double G4Cons::GetCubicVolume ( )

inlinevirtual

Reimplemented from G4VSolid.

GetDeltaPhiAngle()

G4double G4Cons::GetDeltaPhiAngle ( ) const

inline

Here is the caller graph for this function:

GetDPhi()

G4double G4Cons::GetDPhi ( ) const

inline

GetDz()

G4double G4Cons::GetDz ( ) const

inline

GetEntityType()

G4GeometryType G4Cons::GetEntityType ( ) const

virtual

Implements G4VSolid.

Definition at line 2247 of file G4Cons.cc.

{
  return G4String("G4Cons");
}

GetInnerRadiusMinusZ()

G4double G4Cons::GetInnerRadiusMinusZ ( ) const

inline

Here is the caller graph for this function:

GetInnerRadiusPlusZ()

G4double G4Cons::GetInnerRadiusPlusZ ( ) const

inline

Here is the caller graph for this function:

GetOuterRadiusMinusZ()

G4double G4Cons::GetOuterRadiusMinusZ ( ) const

inline

Here is the caller graph for this function:

GetOuterRadiusPlusZ()

G4double G4Cons::GetOuterRadiusPlusZ ( ) const

inline

Here is the caller graph for this function:

GetPointOnSurface()

G4ThreeVector G4Cons::GetPointOnSurface ( ) const

virtual

Reimplemented from G4VSolid.

Definition at line 2292 of file G4Cons.cc.

{   
  // declare working variables
  //
  G4double Aone, Atwo, Athree, Afour, Afive, slin, slout, phi;
  G4double zRand, cosu, sinu, rRand1, rRand2, chose, rone, rtwo, qone, qtwo;
  rone = (fRmax1-fRmax2)/(2.*fDz);
  rtwo = (fRmin1-fRmin2)/(2.*fDz);
  qone=0.; qtwo=0.;
  if(fRmax1!=fRmax2) { qone = fDz*(fRmax1+fRmax2)/(fRmax1-fRmax2); }
  if(fRmin1!=fRmin2) { qtwo = fDz*(fRmin1+fRmin2)/(fRmin1-fRmin2); }
  slin   = std::sqrt(sqr(fRmin1-fRmin2)+sqr(2.*fDz));
  slout  = std::sqrt(sqr(fRmax1-fRmax2)+sqr(2.*fDz));
  Aone   = 0.5*fDPhi*(fRmax2 + fRmax1)*slout;       
  Atwo   = 0.5*fDPhi*(fRmin2 + fRmin1)*slin;
  Athree = 0.5*fDPhi*(fRmax1*fRmax1-fRmin1*fRmin1); 
  Afour  = 0.5*fDPhi*(fRmax2*fRmax2-fRmin2*fRmin2);
  Afive  = fDz*(fRmax1-fRmin1+fRmax2-fRmin2);
  
  phi    = G4RandFlat::shoot(fSPhi,fSPhi+fDPhi);
  cosu   = std::cos(phi);  sinu = std::sin(phi);
  rRand1 = GetRadiusInRing(fRmin1, fRmax1);
  rRand2 = GetRadiusInRing(fRmin2, fRmax2);
  
  if ( (fSPhi == 0.) && fPhiFullCone )  { Afive = 0.; }
  chose  = G4RandFlat::shoot(0.,Aone+Atwo+Athree+Afour+2.*Afive);
 
  if( (chose >= 0.) && (chose < Aone) )
  {
    if(fRmin1 != fRmin2)
    {
      zRand = G4RandFlat::shoot(-1.*fDz,fDz); 
      return G4ThreeVector (rtwo*cosu*(qtwo-zRand),
                            rtwo*sinu*(qtwo-zRand), zRand);
    }
    else
    {
      return G4ThreeVector(fRmin1*cosu, fRmin2*sinu,
                           G4RandFlat::shoot(-1.*fDz,fDz));
    }
  }
  else if( (chose >= Aone) && (chose <= Aone + Atwo) )
  {
    if(fRmax1 != fRmax2)
    {
      zRand = G4RandFlat::shoot(-1.*fDz,fDz); 
      return G4ThreeVector (rone*cosu*(qone-zRand),
                            rone*sinu*(qone-zRand), zRand);
    }    
    else
    {
      return G4ThreeVector(fRmax1*cosu, fRmax2*sinu,
                           G4RandFlat::shoot(-1.*fDz,fDz));
    }
  }
  else if( (chose >= Aone + Atwo) && (chose < Aone + Atwo + Athree) )
  {
    return G4ThreeVector (rRand1*cosu, rRand1*sinu, -1*fDz);
  }
  else if( (chose >= Aone + Atwo + Athree)
        && (chose < Aone + Atwo + Athree + Afour) )
  {
    return G4ThreeVector (rRand2*cosu,rRand2*sinu,fDz);
  }
  else if( (chose >= Aone + Atwo + Athree + Afour)
        && (chose < Aone + Atwo + Athree + Afour + Afive) )
  {
    zRand  = G4RandFlat::shoot(-1.*fDz,fDz);
    rRand1 = G4RandFlat::shoot(fRmin2-((zRand-fDz)/(2.*fDz))*(fRmin1-fRmin2),
                             fRmax2-((zRand-fDz)/(2.*fDz))*(fRmax1-fRmax2)); 
    return G4ThreeVector (rRand1*std::cos(fSPhi),
                          rRand1*std::sin(fSPhi), zRand);
  }
  else
  { 
    zRand  = G4RandFlat::shoot(-1.*fDz,fDz);
    rRand1 = G4RandFlat::shoot(fRmin2-((zRand-fDz)/(2.*fDz))*(fRmin1-fRmin2),
                             fRmax2-((zRand-fDz)/(2.*fDz))*(fRmax1-fRmax2)); 
    return G4ThreeVector (rRand1*std::cos(fSPhi+fDPhi),
                          rRand1*std::sin(fSPhi+fDPhi), zRand);
  }
}

Here is the call graph for this function:

GetRmax1()

G4double G4Cons::GetRmax1 ( ) const

inline

GetRmax2()

G4double G4Cons::GetRmax2 ( ) const

inline

GetRmin1()

G4double G4Cons::GetRmin1 ( ) const

inline

GetRmin2()

G4double G4Cons::GetRmin2 ( ) const

inline

GetSPhi()

G4double G4Cons::GetSPhi ( ) const

inline

GetStartPhiAngle()

G4double G4Cons::GetStartPhiAngle ( ) const

inline

Here is the caller graph for this function:

GetSurfaceArea()

G4double G4Cons::GetSurfaceArea ( )

inlinevirtual

Reimplemented from G4VSolid.

GetZHalfLength()

G4double G4Cons::GetZHalfLength ( ) const

inline

Here is the caller graph for this function:

Initialize()

void G4Cons::Initialize ( )

inlineprivate

InitializeTrigonometry()

void G4Cons::InitializeTrigonometry ( )

inlineprivate

Inside()

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

virtual

Implements G4VSolid.

Definition at line 205 of file G4Cons.cc.

{
  G4double r2, rl, rh, pPhi, tolRMin, tolRMax; // rh2, rl2 ;
  EInside in;

  if (std::fabs(p.z()) > fDz + halfCarTolerance )  { return in = kOutside; }
  else if(std::fabs(p.z()) >= fDz - halfCarTolerance )    { in = kSurface; }
  else                                                    { in = kInside;  }

  r2 = p.x()*p.x() + p.y()*p.y() ;
  rl = 0.5*(fRmin2*(p.z() + fDz) + fRmin1*(fDz - p.z()))/fDz ;
  rh = 0.5*(fRmax2*(p.z()+fDz)+fRmax1*(fDz-p.z()))/fDz;

  // rh2 = rh*rh;

  tolRMin = rl - halfRadTolerance;
  if ( tolRMin < 0 )  { tolRMin = 0; }
  tolRMax = rh + halfRadTolerance;

  if ( (r2<tolRMin*tolRMin) || (r2>tolRMax*tolRMax) ) { return in = kOutside; }

  if (rl) { tolRMin = rl + halfRadTolerance; }
  else    { tolRMin = 0.0; }
  tolRMax = rh - halfRadTolerance;
      
  if (in == kInside) // else it's kSurface already
  {
     if ( (r2 < tolRMin*tolRMin) || (r2 >= tolRMax*tolRMax) ) { in = kSurface; }
  }
  if ( !fPhiFullCone && ((p.x() != 0.0) || (p.y() != 0.0)) )
  {
    pPhi = std::atan2(p.y(),p.x()) ;

    if ( pPhi < fSPhi - halfAngTolerance  )             { pPhi += twopi; }
    else if ( pPhi > fSPhi + fDPhi + halfAngTolerance ) { pPhi -= twopi; }
    
    if ( (pPhi < fSPhi - halfAngTolerance) ||          
         (pPhi > fSPhi + fDPhi + halfAngTolerance) )  { return in = kOutside; }
      
    else if (in == kInside)  // else it's kSurface anyway already
    {
       if ( (pPhi < fSPhi + halfAngTolerance) || 
            (pPhi > fSPhi + fDPhi - halfAngTolerance) )  { in = kSurface; }
    }
  }
  else if ( !fPhiFullCone )  { in = kSurface; }

  return in ;
}

Here is the call graph for this function:

Here is the caller graph for this function:

operator=()

G4Cons & G4Cons::operator=

(

const G4Cons &

rhs

)

Definition at line 172 of file G4Cons.cc.

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

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

   // Copy data
   //
   kRadTolerance = rhs.kRadTolerance;
   kAngTolerance = rhs.kAngTolerance;
   fRmin1 = rhs.fRmin1; fRmin2 = rhs.fRmin2;
   fRmax1 = rhs.fRmax1; fRmax2 = rhs.fRmax2;
   fDz = rhs.fDz; fSPhi = rhs.fSPhi; fDPhi = rhs.fDPhi;
   sinCPhi = rhs.sinCPhi; cosCPhi = rhs.cosCPhi;
   cosHDPhiOT = rhs.cosHDPhiOT; cosHDPhiIT = rhs.cosHDPhiIT;
   sinSPhi = rhs.sinSPhi; cosSPhi = rhs.cosSPhi;
   sinEPhi = rhs.sinEPhi; cosEPhi = rhs.cosEPhi;
   fPhiFullCone = rhs.fPhiFullCone;
   halfCarTolerance = rhs.halfCarTolerance;
   halfRadTolerance = rhs.halfRadTolerance;
   halfAngTolerance = rhs.halfAngTolerance;

   return *this;
}

Here is the call graph for this function:

SetDeltaPhiAngle()

void G4Cons::SetDeltaPhiAngle ( G4double  newDPhi )

inline

Here is the caller graph for this function:

SetInnerRadiusMinusZ()

void G4Cons::SetInnerRadiusMinusZ ( G4double  Rmin1 )

inline

Here is the caller graph for this function:

SetInnerRadiusPlusZ()

void G4Cons::SetInnerRadiusPlusZ ( G4double  Rmin2 )

inline

Here is the caller graph for this function:

SetOuterRadiusMinusZ()

void G4Cons::SetOuterRadiusMinusZ ( G4double  Rmax1 )

inline

Here is the caller graph for this function:

SetOuterRadiusPlusZ()

void G4Cons::SetOuterRadiusPlusZ ( G4double  Rmax2 )

inline

Here is the caller graph for this function:

SetStartPhiAngle()

void G4Cons::SetStartPhiAngle ( G4double  newSPhi, G4bool  trig = true  )

inline

Here is the caller graph for this function:

SetZHalfLength()

void G4Cons::SetZHalfLength ( G4double  newDz )

inline

Here is the caller graph for this function:

StreamInfo()

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

virtual

Reimplemented from G4CSGSolid.

Definition at line 2265 of file G4Cons.cc.

{
  G4int oldprc = os.precision(16);
  os << "-----------------------------------------------------------\n"
     << "    *** Dump for solid - " << GetName() << " ***\n"
     << "    ===================================================\n"
     << " Solid type: G4Cons\n"
     << " Parameters: \n"
     << "   inside  -fDz radius: "  << fRmin1/mm << " mm \n"
     << "   outside -fDz radius: "  << fRmax1/mm << " mm \n"
     << "   inside  +fDz radius: "  << fRmin2/mm << " mm \n"
     << "   outside +fDz radius: "  << fRmax2/mm << " mm \n"
     << "   half length in Z   : "  << fDz/mm << " mm \n"
     << "   starting angle of segment: " << fSPhi/degree << " degrees \n"
     << "   delta angle of segment   : " << fDPhi/degree << " degrees \n"
     << "-----------------------------------------------------------\n";
  os.precision(oldprc);

  return os;
}

Here is the call graph for this function:

SurfaceNormal()

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

virtual

Implements G4VSolid.

Definition at line 490 of file G4Cons.cc.

{
  G4int noSurfaces = 0;
  G4double rho, pPhi;
  G4double distZ, distRMin, distRMax;
  G4double distSPhi = kInfinity, distEPhi = kInfinity;
  G4double tanRMin, secRMin, pRMin, widRMin;
  G4double tanRMax, secRMax, pRMax, widRMax;

  G4ThreeVector norm, sumnorm(0.,0.,0.), nZ = G4ThreeVector(0.,0.,1.);
  G4ThreeVector nR, nr(0.,0.,0.), nPs, nPe;

  distZ = std::fabs(std::fabs(p.z()) - fDz);
  rho   = std::sqrt(p.x()*p.x() + p.y()*p.y());

  tanRMin  = (fRmin2 - fRmin1)*0.5/fDz;
  secRMin  = std::sqrt(1 + tanRMin*tanRMin);
  pRMin    = rho - p.z()*tanRMin;
  widRMin  = fRmin2 - fDz*tanRMin;
  distRMin = std::fabs(pRMin - widRMin)/secRMin;

  tanRMax  = (fRmax2 - fRmax1)*0.5/fDz;
  secRMax  = std::sqrt(1+tanRMax*tanRMax);
  pRMax    = rho - p.z()*tanRMax;
  widRMax  = fRmax2 - fDz*tanRMax;
  distRMax = std::fabs(pRMax - widRMax)/secRMax;

  if (!fPhiFullCone)   // Protected against (0,0,z) 
  {
    if ( rho )
    {
      pPhi = std::atan2(p.y(),p.x());

      if (pPhi  < fSPhi-halfCarTolerance)            { pPhi += twopi; }
      else if (pPhi > fSPhi+fDPhi+halfCarTolerance)  { pPhi -= twopi; }

      distSPhi = std::fabs( pPhi - fSPhi ); 
      distEPhi = std::fabs( pPhi - fSPhi - fDPhi ); 
    }
    else if( !(fRmin1) || !(fRmin2) )
    {
      distSPhi = 0.; 
      distEPhi = 0.; 
    }
    nPs = G4ThreeVector(std::sin(fSPhi), -std::cos(fSPhi), 0);
    nPe = G4ThreeVector(-std::sin(fSPhi+fDPhi), std::cos(fSPhi+fDPhi), 0);
  }
  if ( rho > halfCarTolerance )   
  {
    nR = G4ThreeVector(p.x()/rho/secRMax, p.y()/rho/secRMax, -tanRMax/secRMax);
    if (fRmin1 || fRmin2)
    {
      nr = G4ThreeVector(-p.x()/rho/secRMin,-p.y()/rho/secRMin,tanRMin/secRMin);
    }
  }

  if( distRMax <= halfCarTolerance )
  {
    noSurfaces ++;
    sumnorm += nR;
  }
  if( (fRmin1 || fRmin2) && (distRMin <= halfCarTolerance) )
  {
    noSurfaces ++;
    sumnorm += nr;
  }
  if( !fPhiFullCone )   
  {
    if (distSPhi <= halfAngTolerance)
    {
      noSurfaces ++;
      sumnorm += nPs;
    }
    if (distEPhi <= halfAngTolerance) 
    {
      noSurfaces ++;
      sumnorm += nPe;
    }
  }
  if (distZ <= halfCarTolerance)  
  {
    noSurfaces ++;
    if ( p.z() >= 0.)  { sumnorm += nZ; }
    else               { sumnorm -= nZ; }
  }
  if ( noSurfaces == 0 )
  {
#ifdef G4CSGDEBUG
    G4Exception("G4Cons::SurfaceNormal(p)", "GeomSolids1002",
                JustWarning, "Point p is not on surface !?" );
#endif 
     norm = ApproxSurfaceNormal(p);
  }
  else if ( noSurfaces == 1 )  { norm = sumnorm; }
  else                         { norm = sumnorm.unit(); }

  return norm ;
}

Here is the call graph for this function:

Member Data Documentation

cosCPhi

G4double G4Cons::cosCPhi

private

Definition at line 224 of file G4Cons.hh.

cosEPhi

G4double G4Cons::cosEPhi

private

Definition at line 224 of file G4Cons.hh.

cosHDPhiIT

G4double G4Cons::cosHDPhiIT

private

Definition at line 224 of file G4Cons.hh.

cosHDPhiOT

G4double G4Cons::cosHDPhiOT

private

Definition at line 224 of file G4Cons.hh.

cosSPhi

G4double G4Cons::cosSPhi

private

Definition at line 224 of file G4Cons.hh.

fDPhi

G4double G4Cons::fDPhi

private

Definition at line 220 of file G4Cons.hh.

fDz

G4double G4Cons::fDz

private

Definition at line 220 of file G4Cons.hh.

fPhiFullCone

G4bool G4Cons::fPhiFullCone

private

Definition at line 229 of file G4Cons.hh.

fRmax1

G4double G4Cons::fRmax1

private

Definition at line 220 of file G4Cons.hh.

fRmax2

G4double G4Cons::fRmax2

private

Definition at line 220 of file G4Cons.hh.

fRmin1

G4double G4Cons::fRmin1

private

Definition at line 220 of file G4Cons.hh.

fRmin2

G4double G4Cons::fRmin2

private

Definition at line 220 of file G4Cons.hh.

fSPhi

G4double G4Cons::fSPhi

private

Definition at line 220 of file G4Cons.hh.

halfAngTolerance

G4double G4Cons::halfAngTolerance

private

Definition at line 233 of file G4Cons.hh.

halfCarTolerance

G4double G4Cons::halfCarTolerance

private

Definition at line 233 of file G4Cons.hh.

halfRadTolerance

G4double G4Cons::halfRadTolerance

private

Definition at line 233 of file G4Cons.hh.

kAngTolerance

G4double G4Cons::kAngTolerance

private

Definition at line 216 of file G4Cons.hh.

kRadTolerance

G4double G4Cons::kRadTolerance

private

Definition at line 216 of file G4Cons.hh.

sinCPhi

G4double G4Cons::sinCPhi

private

Definition at line 224 of file G4Cons.hh.

sinEPhi

G4double G4Cons::sinEPhi

private

Definition at line 224 of file G4Cons.hh.

sinSPhi

G4double G4Cons::sinSPhi

private

Definition at line 224 of file G4Cons.hh.

---

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

• G4Cons.hh
• G4Cons.cc