G4GenericTrap Class Reference

#include <G4GenericTrap.hh>

Inheritance diagram for G4GenericTrap:

Public Member Functions
	G4GenericTrap (const G4String &name, G4double halfZ, const std::vector< G4TwoVector > &vertices)
	~G4GenericTrap ()
G4double	GetZHalfLength () const
G4int	GetNofVertices () const
G4TwoVector	GetVertex (G4int index) const
const std::vector< G4TwoVector > &	GetVertices () const
G4double	GetTwistAngle (G4int index) const
G4bool	IsTwisted () const
G4int	GetVisSubdivisions () const
void	SetVisSubdivisions (G4int subdiv)
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=false, G4bool *validNorm=0, G4ThreeVector *n=0) const
G4double	DistanceToOut (const G4ThreeVector &p) const
G4bool	CalculateExtent (const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pmin, G4double &pmax) const
G4GeometryType	GetEntityType () const
G4VSolid *	Clone () const
std::ostream &	StreamInfo (std::ostream &os) const
G4ThreeVector	GetPointOnSurface () const
G4double	GetCubicVolume ()
G4double	GetSurfaceArea ()
G4Polyhedron *	GetPolyhedron () const
void	DescribeYourselfTo (G4VGraphicsScene &scene) const
G4VisExtent	GetExtent () const
G4Polyhedron *	CreatePolyhedron () const
G4NURBS *	CreateNURBS () const
	G4GenericTrap (__void__ &)
	G4GenericTrap (const G4GenericTrap &rhs)
G4GenericTrap &	operator= (const G4GenericTrap &rhs)
Protected Attributes
G4Polyhedron *	fpPolyhedron

---

Detailed Description

Definition at line 79 of file G4GenericTrap.hh.

---

Constructor & Destructor Documentation

G4GenericTrap::G4GenericTrap	(	const G4String &	name,
		G4double	halfZ,
		const std::vector< G4TwoVector > &	vertices
	)

Definition at line 68 of file G4GenericTrap.cc.

References FatalErrorInArgument, G4endl, G4Exception(), JustWarning, and G4VSolid::kCarTolerance.

Referenced by Clone().

  : G4VSolid(name),
    fpPolyhedron(0),
    fDz(halfZ),
    fVertices(),
    fIsTwisted(false),
    fTessellatedSolid(0),
    fMinBBoxVector(G4ThreeVector(0,0,0)),
    fMaxBBoxVector(G4ThreeVector(0,0,0)),
    fVisSubdivisions(0),
    fSurfaceArea(0.),
    fCubicVolume(0.)
   
{
  // General constructor
  const G4double min_length=5*1.e-6;
  G4double length = 0.;
  G4int k=0;
  G4String errorDescription = "InvalidSetup in \" ";
  errorDescription += name;
  errorDescription += "\""; 
  
  // Check vertices size

  if ( G4int(vertices.size()) != fgkNofVertices )
  {
    G4Exception("G4GenericTrap::G4GenericTrap()", "GeomSolids0002",
                FatalErrorInArgument, "Number of vertices != 8");
  }            
  
  // Check dZ
  // 
  if (halfZ < kCarTolerance)
  {
     G4Exception("G4GenericTrap::G4GenericTrap()", "GeomSolids0002",
                FatalErrorInArgument, "dZ is too small or negative");
  }           
 
  // Check Ordering and Copy vertices 
  //
  if(CheckOrder(vertices))
  {
    for (G4int i=0; i<fgkNofVertices; ++i) {fVertices.push_back(vertices[i]);}
  }
  else
  { 
    for (G4int i=0; i <4; ++i) {fVertices.push_back(vertices[3-i]);}
    for (G4int i=0; i <4; ++i) {fVertices.push_back(vertices[7-i]);}
  }

   // Check length of segments and Adjust
  // 
  for (G4int j=0; j < 2; j++)
  {
    for (G4int i=1; i<4; ++i)
    {
      k = j*4+i;
      length = (fVertices[k]-fVertices[k-1]).mag();
      if ( ( length < min_length) && ( length > kCarTolerance ) )
      {
        std::ostringstream message;
        message << "Length segment is too small." << G4endl
                << "Distance between " << fVertices[k-1] << " and "
                << fVertices[k] << " is only " << length << " mm !"; 
        G4Exception("G4GenericTrap::G4GenericTrap()", "GeomSolids1001",
                    JustWarning, message, "Vertices will be collapsed.");
        fVertices[k]=fVertices[k-1];
      }
    }
  }

  // Compute Twist
  //
  for( G4int i=0; i<4; i++) { fTwist[i]=0.; }
  fIsTwisted = ComputeIsTwisted();

  // Compute Bounding Box 
  //
  ComputeBBox();
  
  // If not twisted - create tessellated solid 
  // (an alternative implementation for testing)
  //
#ifdef G4TESS_TEST
   if ( !fIsTwisted )  { fTessellatedSolid = CreateTessellatedSolid(); }
#endif
}

G4GenericTrap::~G4GenericTrap

(

)

Definition at line 180 of file G4GenericTrap.cc.

{
  // Destructor
  delete fTessellatedSolid;
}

G4GenericTrap::G4GenericTrap

(

__void__ &

)

Definition at line 159 of file G4GenericTrap.cc.

  : G4VSolid(a),
    fpPolyhedron(0),
    fDz(0.),
    fVertices(),
    fIsTwisted(false),
    fTessellatedSolid(0),
    fMinBBoxVector(G4ThreeVector(0,0,0)),
    fMaxBBoxVector(G4ThreeVector(0,0,0)),
    fVisSubdivisions(0),
    fSurfaceArea(0.),
    fCubicVolume(0.)
{
  // Fake default constructor - sets only member data and allocates memory
  //                            for usage restricted to object persistency.

  for (size_t i=0; i<4; ++i)  { fTwist[i]=0.; }
}

G4GenericTrap::G4GenericTrap

(

const G4GenericTrap &

rhs

)

Definition at line 188 of file G4GenericTrap.cc.

References fTessellatedSolid, and fTwist.

  : G4VSolid(rhs),
    fpPolyhedron(0), fDz(rhs.fDz), fVertices(rhs.fVertices),
    fIsTwisted(rhs.fIsTwisted), fTessellatedSolid(0),
    fMinBBoxVector(rhs.fMinBBoxVector), fMaxBBoxVector(rhs.fMaxBBoxVector),
    fVisSubdivisions(rhs.fVisSubdivisions),
    fSurfaceArea(rhs.fSurfaceArea), fCubicVolume(rhs.fCubicVolume) 
{
   for (size_t i=0; i<4; ++i)  { fTwist[i] = rhs.fTwist[i]; }
#ifdef G4TESS_TEST
   if (rhs.fTessellatedSolid && !fIsTwisted )
   { fTessellatedSolid = CreateTessellatedSolid(); } 
#endif
}

---

Member Function Documentation

G4bool G4GenericTrap::CalculateExtent	(	const EAxis	pAxis,
		const G4VoxelLimits &	pVoxelLimit,
		const G4AffineTransform &	pTransform,
		G4double &	pmin,
		G4double &	pmax
	)			const [virtual]

Implements G4VSolid.

Definition at line 1169 of file G4GenericTrap.cc.

References G4TessellatedSolid::CalculateExtent(), G4VSolid::ClipBetweenSections(), G4VSolid::ClipCrossSection(), G4VoxelLimits::GetMaxExtent(), G4VoxelLimits::GetMaxXExtent(), G4VoxelLimits::GetMaxYExtent(), G4VoxelLimits::GetMaxZExtent(), G4VoxelLimits::GetMinExtent(), G4VoxelLimits::GetMinXExtent(), G4VoxelLimits::GetMinYExtent(), G4VoxelLimits::GetMinZExtent(), Inside(), G4AffineTransform::Inverse(), G4AffineTransform::IsRotated(), G4VoxelLimits::IsXLimited(), G4VoxelLimits::IsYLimited(), G4VoxelLimits::IsZLimited(), G4VSolid::kCarTolerance, kOutside, kXAxis, kYAxis, kZAxis, G4AffineTransform::NetTranslation(), and G4AffineTransform::TransformPoint().

{
#ifdef G4TESS_TEST
  if ( fTessellatedSolid )
  {
    return fTessellatedSolid->CalculateExtent(pAxis, pVoxelLimit,
                                              pTransform, pMin, pMax);
  }     
#endif 

  // Computes bounding vectors for a shape
  //
  G4double Dx,Dy;
  G4ThreeVector minVec = GetMinimumBBox();
  G4ThreeVector maxVec = GetMaximumBBox();
  Dx = 0.5*(maxVec.x()- minVec.x());
  Dy = 0.5*(maxVec.y()- minVec.y());

  if (!pTransform.IsRotated())
  {
    // Special case handling for unrotated shapes
    // Compute x/y/z mins and maxs respecting limits, with early returns
    // if outside limits. Then switch() on pAxis
    //
    G4double xoffset,xMin,xMax;
    G4double yoffset,yMin,yMax;
    G4double zoffset,zMin,zMax;

    xoffset=pTransform.NetTranslation().x();
    xMin=xoffset-Dx;
    xMax=xoffset+Dx;
    if (pVoxelLimit.IsXLimited())
    {
      if ( (xMin>pVoxelLimit.GetMaxXExtent()+kCarTolerance)
        || (xMax<pVoxelLimit.GetMinXExtent()-kCarTolerance) )
      {
        return false;
      }
      else
      {
        if (xMin<pVoxelLimit.GetMinXExtent())
        {
          xMin=pVoxelLimit.GetMinXExtent();
        }
        if (xMax>pVoxelLimit.GetMaxXExtent())
        {
          xMax=pVoxelLimit.GetMaxXExtent();
        }
      }
    }

    yoffset=pTransform.NetTranslation().y();
    yMin=yoffset-Dy;
    yMax=yoffset+Dy;
    if (pVoxelLimit.IsYLimited())
    {
      if ( (yMin>pVoxelLimit.GetMaxYExtent()+kCarTolerance)
        || (yMax<pVoxelLimit.GetMinYExtent()-kCarTolerance) )
      {
        return false;
      }
      else
      {
        if (yMin<pVoxelLimit.GetMinYExtent())
        {
          yMin=pVoxelLimit.GetMinYExtent();
        }
        if (yMax>pVoxelLimit.GetMaxYExtent())
        {
          yMax=pVoxelLimit.GetMaxYExtent();
        }
      }
    }

    zoffset=pTransform.NetTranslation().z();
    zMin=zoffset-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();
        }
      }
    }

    switch (pAxis)
    {
      case kXAxis:
        pMin = xMin;
        pMax = xMax;
        break;
      case kYAxis:
        pMin = yMin;
        pMax = yMax;
        break;
      case kZAxis:
        pMin = zMin;
        pMax = zMax;
        break;
      default:
        break;
    }
    pMin-=kCarTolerance;
    pMax+=kCarTolerance;

    return true;
  }
  else
  {
    // General rotated case - create and clip mesh to boundaries

    G4bool existsAfterClip=false;
    G4ThreeVectorList *vertices;

    pMin=+kInfinity;
    pMax=-kInfinity;

    // Calculate rotated vertex coordinates
    //
    vertices=CreateRotatedVertices(pTransform);
    ClipCrossSection(vertices,0,pVoxelLimit,pAxis,pMin,pMax);
    ClipCrossSection(vertices,4,pVoxelLimit,pAxis,pMin,pMax);
    ClipBetweenSections(vertices,0,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;
  }
}

G4VSolid * G4GenericTrap::Clone

(

)

const [virtual]

Reimplemented from G4VSolid.

Definition at line 1396 of file G4GenericTrap.cc.

References G4GenericTrap().

{
  return new G4GenericTrap(*this);
}

G4NURBS * G4GenericTrap::CreateNURBS

(

)

const [virtual]

Reimplemented from G4VSolid.

Definition at line 2219 of file G4GenericTrap.cc.

References G4TessellatedSolid::CreateNURBS().

{
#ifdef G4TESS_TEST
  if ( fTessellatedSolid )
  { 
    return fTessellatedSolid->CreateNURBS();
  }
#endif

  // Computes bounding vectors for the shape
  //
  G4double Dx,Dy;
  G4ThreeVector minVec = GetMinimumBBox();
  G4ThreeVector maxVec = GetMaximumBBox();
  Dx = 0.5*(maxVec.x()- minVec.y());
  Dy = 0.5*(maxVec.y()- minVec.y());

  return new G4NURBSbox (Dx, Dy, fDz);
}    

G4Polyhedron * G4GenericTrap::CreatePolyhedron

(

)

const [virtual]

Reimplemented from G4VSolid.

Definition at line 2123 of file G4GenericTrap.cc.

References G4PolyhedronArbitrary::AddFacet(), G4PolyhedronArbitrary::AddVertex(), G4TessellatedSolid::CreatePolyhedron(), GetTwistAngle(), GetVisSubdivisions(), G4PolyhedronArbitrary::InvertFacets(), and G4PolyhedronArbitrary::SetReferences().

Referenced by GetPolyhedron().

{

#ifdef G4TESS_TEST
  if ( fTessellatedSolid )
  { 
    return fTessellatedSolid->CreatePolyhedron();
  }  
#endif 
  
  // Approximation of Twisted Side
  // Construct extra Points, if Twisted Side
  //
  G4PolyhedronArbitrary* polyhedron;
  size_t nVertices, nFacets;

  G4int subdivisions=0;
  G4int i;
  if(fIsTwisted)
  {
    if ( GetVisSubdivisions()!= 0 )
    {
      subdivisions=GetVisSubdivisions();
    }
    else
    {
      // Estimation of Number of Subdivisions for smooth visualisation
      //
      G4double maxTwist=0.;
      for(i=0; i<4; i++)
      {
        if(GetTwistAngle(i)>maxTwist) { maxTwist=GetTwistAngle(i); }
      }

      // Computes bounding vectors for the shape
      //
      G4double Dx,Dy;
      G4ThreeVector minVec = GetMinimumBBox();
      G4ThreeVector maxVec = GetMaximumBBox();
      Dx = 0.5*(maxVec.x()- minVec.y());
      Dy = 0.5*(maxVec.y()- minVec.y());
      if (Dy > Dx)  { Dx=Dy; }
    
      subdivisions=8*G4int(maxTwist/(Dx*Dx*Dx)*fDz);
      if (subdivisions<4)  { subdivisions=4; }
      if (subdivisions>30) { subdivisions=30; }
    }
  }
  G4int sub4=4*subdivisions;
  nVertices = 8+subdivisions*4;
  nFacets = 6+subdivisions*4;
  G4double cf=1./(subdivisions+1);
  polyhedron = new G4PolyhedronArbitrary (nVertices, nFacets);

  // Add Vertex
  //
  for (i=0;i<4;i++)
  {
    polyhedron->AddVertex(G4ThreeVector(fVertices[i].x(),
                                        fVertices[i].y(),-fDz));
  }
  for( i=0;i<subdivisions;i++)
  {
    for(G4int j=0;j<4;j++)
    {
      G4TwoVector u=fVertices[j]+cf*(i+1)*( fVertices[j+4]-fVertices[j]);
      polyhedron->AddVertex(G4ThreeVector(u.x(),u.y(),-fDz+cf*2*fDz*(i+1)));
    }    
  }
  for (i=4;i<8;i++)
  {
    polyhedron->AddVertex(G4ThreeVector(fVertices[i].x(),
                                        fVertices[i].y(),fDz));
  }

  // Add Facets
  //
  polyhedron->AddFacet(1,4,3,2);  //Z-plane
  for (i=0;i<subdivisions+1;i++)
  {
    G4int is=i*4;
    polyhedron->AddFacet(5+is,8+is,4+is,1+is);
    polyhedron->AddFacet(8+is,7+is,3+is,4+is);
    polyhedron->AddFacet(7+is,6+is,2+is,3+is);
    polyhedron->AddFacet(6+is,5+is,1+is,2+is); 
  }
  polyhedron->AddFacet(5+sub4,6+sub4,7+sub4,8+sub4);  //Z-plane

  polyhedron->SetReferences();
  polyhedron->InvertFacets();

  return (G4Polyhedron*) polyhedron;
}

void G4GenericTrap::DescribeYourselfTo

(

G4VGraphicsScene &

scene

)

const [virtual]

Implements G4VSolid.

Definition at line 2086 of file G4GenericTrap.cc.

References G4VGraphicsScene::AddSolid(), and G4TessellatedSolid::DescribeYourselfTo().

{

#ifdef G4TESS_TEST
  if ( fTessellatedSolid )
  { 
    return fTessellatedSolid->DescribeYourselfTo(scene);
  }
#endif  
  
  scene.AddSolid(*this);
}

G4double G4GenericTrap::DistanceToIn

(

const G4ThreeVector &

p

)

const [virtual]

Implements G4VSolid.

Definition at line 789 of file G4GenericTrap.cc.

References G4TessellatedSolid::DistanceToIn().

{
  // Computes the closest distance from given point to this shape

#ifdef G4TESS_TEST
  if ( fTessellatedSolid )
  {
    return fTessellatedSolid->DistanceToIn(p);
  }  
#endif
 
  G4double safz = std::fabs(p.z())-fDz;
  if(safz<0) { safz=0; }

  G4int iseg;
  G4double safe  = safz;
  G4double safxy = safz;
 
  for (iseg=0; iseg<4; iseg++)
  { 
    safxy = SafetyToFace(p,iseg);
    if (safxy>safe)  { safe=safxy; }
  }

  return safe;
}

G4double G4GenericTrap::DistanceToIn	(	const G4ThreeVector &	p,
		const G4ThreeVector &	v
	)			const [virtual]

Implements G4VSolid.

Definition at line 719 of file G4GenericTrap.cc.

References G4TessellatedSolid::DistanceToIn(), Inside(), G4VSolid::kCarTolerance, kOutside, and CLHEP::detail::n.

{
#ifdef G4TESS_TEST
  if ( fTessellatedSolid )
  { 
    return fTessellatedSolid->DistanceToIn(p, v);
  }  
#endif
    
  static const G4double halfCarTolerance=kCarTolerance*0.5;

  G4double dist[5];
  G4ThreeVector n;

  // Check lateral faces
  //
  G4int i;
  for (i=0; i<4; i++)
  {
    dist[i]=DistToPlane(p, v, i);  
  }

  // Check Z planes
  //
  dist[4]=kInfinity;
  if (std::fabs(p.z())>fDz-halfCarTolerance)
  {
    if (v.z())
    {
      G4ThreeVector pt;
      if (p.z()>0)
      {
        dist[4] = (fDz-p.z())/v.z();
      }
      else
      {   
        dist[4] = (-fDz-p.z())/v.z();
      }   
      if (dist[4]<-halfCarTolerance)
      {
        dist[4]=kInfinity;
      }
      else
      {
        if(dist[4]<halfCarTolerance)
        {
          if(p.z()>0)  { n=G4ThreeVector(0,0,1); }
          else         { n=G4ThreeVector(0,0,-1); }
          if (n.dot(v)<0) { dist[4]=0.; }
          else            { dist[4]=kInfinity; }
        }
        pt=p+dist[4]*v;
        if (Inside(pt)==kOutside)  { dist[4]=kInfinity; }
      }
    }
  }   
  G4double distmin = dist[0];
  for (i=1;i<5;i++)
  {
    if (dist[i] < distmin)  { distmin = dist[i]; }
  }

  if (distmin<halfCarTolerance)  { distmin=0.; }

  return distmin;
}    

G4double G4GenericTrap::DistanceToOut

(

const G4ThreeVector &

p

)

const [virtual]

Implements G4VSolid.

Definition at line 1142 of file G4GenericTrap.cc.

References G4TessellatedSolid::DistanceToOut().

{

#ifdef G4TESS_TEST
  if ( fTessellatedSolid )
  { 
    return fTessellatedSolid->DistanceToOut(p);
  }  
#endif

  G4double safz = fDz-std::fabs(p.z());
  if (safz<0) { safz = 0; }

  G4double safe  = safz;
  G4double safxy = safz;
 
  for (G4int iseg=0; iseg<4; iseg++)
  { 
    safxy = std::fabs(SafetyToFace(p,iseg));
    if (safxy < safe)  { safe = safxy; }
  }

  return safe;
}    

G4double G4GenericTrap::DistanceToOut	(	const G4ThreeVector &	p,
		const G4ThreeVector &	v,
		const G4bool	calcNorm = false,
		G4bool *	validNorm = 0,
		G4ThreeVector *	n = 0
	)			const [virtual]

Implements G4VSolid.

Definition at line 895 of file G4GenericTrap.cc.

References G4TessellatedSolid::DistanceToOut(), G4VSolid::DumpInfo(), G4endl, G4Exception(), JustWarning, G4VSolid::kCarTolerance, and kOutside.

{
#ifdef G4TESS_TEST
  if ( fTessellatedSolid )
  { 
    return fTessellatedSolid->DistanceToOut(p, v, calcNorm, validNorm, n);
  }
#endif

  static const G4double halfCarTolerance=kCarTolerance*0.5;

  G4double distmin;
  G4bool lateral_cross = false;
  ESide side = kUndefined;
 
  if (calcNorm)  { *validNorm=true; } // All normals are valid

  if (v.z() < 0)
  {
    distmin=(-fDz-p.z())/v.z();
    if (calcNorm) { side=kMZ; *n=G4ThreeVector(0,0,-1); }
  }
  else
  {
    if (v.z() > 0)
    { 
      distmin = (fDz-p.z())/v.z(); 
      if (calcNorm) { side=kPZ; *n=G4ThreeVector(0,0,1); } 
    }
    else            { distmin = kInfinity; }
  }      

  G4double dz2 =0.5/fDz;
  G4double xa,xb,xc,xd;
  G4double ya,yb,yc,yd;

  for (G4int ipl=0; ipl<4; ipl++)
  {
    G4int j = (ipl+1)%4;
    xa=fVertices[ipl].x();
    ya=fVertices[ipl].y();
    xb=fVertices[ipl+4].x();
    yb=fVertices[ipl+4].y();
    xc=fVertices[j].x();
    yc=fVertices[j].y();
    xd=fVertices[4+j].x();
    yd=fVertices[4+j].y();

    if ( ((std::fabs(xb-xd)+std::fabs(yb-yd))<halfCarTolerance)
      || ((std::fabs(xa-xc)+std::fabs(ya-yc))<halfCarTolerance) )
    {
      G4double q=DistToTriangle(p,v,ipl) ;
      if ( (q>=0) && (q<distmin) )
      {
        distmin=q;
        lateral_cross=true;
        side=ESide(ipl+1);
      }
      continue;
    }
    G4double tx1 =dz2*(xb-xa);
    G4double ty1 =dz2*(yb-ya);
    G4double tx2 =dz2*(xd-xc);
    G4double ty2 =dz2*(yd-yc);
    G4double dzp =fDz+p.z();
    G4double xs1 =xa+tx1*dzp;
    G4double ys1 =ya+ty1*dzp;
    G4double xs2 =xc+tx2*dzp;
    G4double ys2 =yc+ty2*dzp;
    G4double dxs =xs2-xs1;
    G4double dys =ys2-ys1;
    G4double dtx =tx2-tx1;
    G4double dty =ty2-ty1;
    G4double a = (dtx*v.y()-dty*v.x()+(tx1*ty2-tx2*ty1)*v.z())*v.z();
    G4double b = dxs*v.y()-dys*v.x()+(dtx*p.y()-dty*p.x()+ty2*xs1-ty1*xs2
               + tx1*ys2-tx2*ys1)*v.z();
    G4double c=dxs*p.y()-dys*p.x()+xs1*ys2-xs2*ys1;
    G4double q=kInfinity;

    if (std::fabs(a) < kCarTolerance)
    {           
      if (std::fabs(b) < kCarTolerance) { continue; }
      q=-c/b;
         
      // Check for Point on the Surface
      //
      if ((q > -halfCarTolerance) && (q < distmin))
      {
        if (q < halfCarTolerance)
        {
          if (NormalToPlane(p,ipl).dot(v)<0.)  { continue; }
        }
        distmin =q;
        lateral_cross=true;
        side=ESide(ipl+1);
      }   
      continue;
    }
    G4double d=b*b-4*a*c;
    if (d >= 0.)
    {
      if (a > 0)  { q=0.5*(-b-std::sqrt(d))/a; }
      else        { q=0.5*(-b+std::sqrt(d))/a; }

      // Check for Point on the Surface
      //
      if (q > -halfCarTolerance )
      {
        if (q < distmin)
        {
          if(q < halfCarTolerance)
          {
            if (NormalToPlane(p,ipl).dot(v)<0.)  // Check second root
            {
              if (a > 0)  { q=0.5*(-b+std::sqrt(d))/a; }
              else        { q=0.5*(-b-std::sqrt(d))/a; }
              if (( q > halfCarTolerance) && (q < distmin))
              {
                distmin=q;
                lateral_cross = true;
                side=ESide(ipl+1);
              }
              continue;
            }
          }
          distmin = q;
          lateral_cross = true;
          side=ESide(ipl+1);
        }
      }
      else
      {
        if (a > 0)  { q=0.5*(-b+std::sqrt(d))/a; }
        else        { q=0.5*(-b-std::sqrt(d))/a; }

        // Check for Point on the Surface
        //
        if ((q > -halfCarTolerance) && (q < distmin))
        {
          if (q < halfCarTolerance)
          {
            if (NormalToPlane(p,ipl).dot(v)<0.)  // Check second root
            {
              if (a > 0)  { q=0.5*(-b-std::sqrt(d))/a; }
              else        { q=0.5*(-b+std::sqrt(d))/a; }
              if ( ( q > halfCarTolerance) && (q < distmin) )
              {
                distmin=q;
                lateral_cross = true;
                side=ESide(ipl+1);
              }
              continue;
            }  
          }
          distmin =q;
          lateral_cross = true;
          side=ESide(ipl+1);
        }   
      }
    }
  }
  if (!lateral_cross)  // Make sure that track crosses the top or bottom
  {
    if (distmin >= kInfinity)  { distmin=kCarTolerance; }
    G4ThreeVector pt=p+distmin*v;
    
    // Check if propagated point is in the polygon
    //
    G4int i=0;
    if (v.z()>0.) { i=4; }
    std::vector<G4TwoVector> xy;
    for ( G4int j=0; j<4; j++)  { xy.push_back(fVertices[i+j]); }

    // Check Inside
    //
    if (InsidePolygone(pt,xy)==kOutside)
    { 
      if(calcNorm)
      {
        if (v.z()>0) {side= kPZ; *n = G4ThreeVector(0,0,1);}
        else         { side=kMZ; *n = G4ThreeVector(0,0,-1);}
      } 
      return 0.;
    }
    else
    {
      if(v.z()>0) {side=kPZ;}
      else        {side=kMZ;}
    }
  }

  if (calcNorm)
  {
    G4ThreeVector pt=p+v*distmin;     
    switch (side)
    {
      case kXY0:
        *n=NormalToPlane(pt,0);
        break;
      case kXY1:
        *n=NormalToPlane(pt,1);
        break;
      case kXY2:
        *n=NormalToPlane(pt,2);
        break;
      case kXY3:
        *n=NormalToPlane(pt,3);
        break;
      case kPZ:
        *n=G4ThreeVector(0,0,1);
        break;
      case kMZ:
        *n=G4ThreeVector(0,0,-1);
        break;
      default:
        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("G4GenericTrap::DistanceToOut(p,v,..)",
                    "GeomSolids1002", JustWarning, message);
        break;
     }
  }
 
  if (distmin<halfCarTolerance)  { distmin=0.; }
 
  return distmin;
}    

G4double G4GenericTrap::GetCubicVolume

(

)

[virtual]

Reimplemented from G4VSolid.

Definition at line 1514 of file G4GenericTrap.cc.

References G4VSolid::GetCubicVolume().

{
  if(fCubicVolume != 0.) {;}
  else   { fCubicVolume = G4VSolid::GetCubicVolume(); }
  return fCubicVolume;
}

G4GeometryType G4GenericTrap::GetEntityType

(

)

const [virtual]

Implements G4VSolid.

Definition at line 1389 of file G4GenericTrap.cc.

Referenced by StreamInfo().

{
  return G4String("G4GenericTrap");
}

G4VisExtent G4GenericTrap::GetExtent

(

)

const [virtual]

Reimplemented from G4VSolid.

Definition at line 2101 of file G4GenericTrap.cc.

References G4TessellatedSolid::GetExtent().

{
  // Computes bounding vectors for the shape

#ifdef G4TESS_TEST
  if ( fTessellatedSolid )
  { 
    return fTessellatedSolid->GetExtent();
  } 
#endif
   
  G4double Dx,Dy;
  G4ThreeVector minVec = GetMinimumBBox();
  G4ThreeVector maxVec = GetMaximumBBox();
  Dx = 0.5*(maxVec.x()- minVec.x());
  Dy = 0.5*(maxVec.y()- minVec.y());

  return G4VisExtent (-Dx, Dx, -Dy, Dy, -fDz, fDz); 
}    

G4int G4GenericTrap::GetNofVertices

(

)

const [inline]

Definition at line 46 of file G4GenericTrap.icc.

{
  return fVertices.size();
}  

G4ThreeVector G4GenericTrap::GetPointOnSurface

(

)

const [virtual]

Reimplemented from G4VSolid.

Definition at line 1426 of file G4GenericTrap.cc.

References G4UniformRand, and G4TessellatedSolid::GetPointOnSurface().

{

#ifdef G4TESS_TEST
  if ( fTessellatedSolid )
  { 
    return fTessellatedSolid->GetPointOnSurface();
  }  
#endif

  G4ThreeVector point;
  G4TwoVector u,v,w;
  G4double rand,area,chose,cf,lambda0,lambda1,alfa,beta,zp;
  G4int ipl,j;
 
  std::vector<G4ThreeVector> vertices;
  for (G4int i=0; i<4;i++)
  {
    vertices.push_back(G4ThreeVector(fVertices[i].x(),fVertices[i].y(),-fDz));
  }
  for (G4int i=4; i<8;i++)
  {
    vertices.push_back(G4ThreeVector(fVertices[i].x(),fVertices[i].y(),fDz));
  }

  // Surface Area of Planes(only estimation for twisted)
  //
  G4double Surface0=GetFaceSurfaceArea(vertices[0],vertices[1],
                                       vertices[2],vertices[3]);//-fDz plane 
  G4double Surface1=GetFaceSurfaceArea(vertices[0],vertices[1],
                                       vertices[5],vertices[4]);// Lat plane
  G4double Surface2=GetFaceSurfaceArea(vertices[3],vertices[0],
                                       vertices[4],vertices[7]);// Lat plane
  G4double Surface3=GetFaceSurfaceArea(vertices[2],vertices[3],
                                       vertices[7],vertices[6]);// Lat plane
  G4double Surface4=GetFaceSurfaceArea(vertices[2],vertices[1],
                                       vertices[5],vertices[6]);// Lat plane
  G4double Surface5=GetFaceSurfaceArea(vertices[4],vertices[5],
                                       vertices[6],vertices[7]);// fDz plane
  rand = G4UniformRand();
  area = Surface0+Surface1+Surface2+Surface3+Surface4+Surface5;
  chose = rand*area;

  if ( ( chose < Surface0)
    || ( chose > (Surface0+Surface1+Surface2+Surface3+Surface4)) )
  {                                        // fDz or -fDz Plane
    ipl = G4int(G4UniformRand()*4);
    j = (ipl+1)%4;
    if(chose < Surface0)
    { 
      zp = -fDz;
      u = fVertices[ipl]; v = fVertices[j];
      w = fVertices[(ipl+3)%4]; 
    }
    else
    {
      zp = fDz;
      u = fVertices[ipl+4]; v = fVertices[j+4];
      w = fVertices[(ipl+3)%4+4]; 
    }
    alfa = G4UniformRand();
    beta = G4UniformRand();
    lambda1=alfa*beta;
    lambda0=alfa-lambda1;
    v = v-u;
    w = w-u;
    v = u+lambda0*v+lambda1*w;
  }
  else                                     // Lateral Plane Twisted or Not
  {
    if (chose < Surface0+Surface1) { ipl=0; }
    else if (chose < Surface0+Surface1+Surface2) { ipl=1; }
    else if (chose < Surface0+Surface1+Surface2+Surface3) { ipl=2; }
    else { ipl=3; }
    j = (ipl+1)%4;
    zp = -fDz+G4UniformRand()*2*fDz;
    cf = 0.5*(fDz-zp)/fDz;
    u = fVertices[ipl+4]+cf*( fVertices[ipl]-fVertices[ipl+4]);
    v = fVertices[j+4]+cf*(fVertices[j]-fVertices[j+4]); 
    v = u+(v-u)*G4UniformRand();
  }
  point=G4ThreeVector(v.x(),v.y(),zp);

  return point;
}

G4Polyhedron * G4GenericTrap::GetPolyhedron

(

)

const [virtual]

Reimplemented from G4VSolid.

Definition at line 2064 of file G4GenericTrap.cc.

References CreatePolyhedron(), fpPolyhedron, G4Polyhedron::GetNumberOfRotationStepsAtTimeOfCreation(), and G4TessellatedSolid::GetPolyhedron().

{

#ifdef G4TESS_TEST
  if ( fTessellatedSolid )
  { 
    return fTessellatedSolid->GetPolyhedron();
  }
#endif  
  
  if ( (!fpPolyhedron)
    || (fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() !=
        fpPolyhedron->GetNumberOfRotationSteps()) )
  {
    delete fpPolyhedron;
    fpPolyhedron = CreatePolyhedron();
  }
  return fpPolyhedron;
}    

G4double G4GenericTrap::GetSurfaceArea

(

)

[virtual]

Reimplemented from G4VSolid.

Definition at line 1523 of file G4GenericTrap.cc.

References G4VSolid::GetSurfaceArea().

{
  if(fSurfaceArea != 0.) {;}
  else
  {
    std::vector<G4ThreeVector> vertices;
    for (G4int i=0; i<4;i++)
    {
      vertices.push_back(G4ThreeVector(fVertices[i].x(),fVertices[i].y(),-fDz));
    }
    for (G4int i=4; i<8;i++)
    {
      vertices.push_back(G4ThreeVector(fVertices[i].x(),fVertices[i].y(),fDz));
    }

    // Surface Area of Planes(only estimation for twisted)
    //
    G4double fSurface0=GetFaceSurfaceArea(vertices[0],vertices[1],
                                          vertices[2],vertices[3]);//-fDz plane
    G4double fSurface1=GetFaceSurfaceArea(vertices[0],vertices[1],
                                          vertices[5],vertices[4]);// Lat plane
    G4double fSurface2=GetFaceSurfaceArea(vertices[3],vertices[0],
                                          vertices[4],vertices[7]);// Lat plane 
    G4double fSurface3=GetFaceSurfaceArea(vertices[2],vertices[3],
                                          vertices[7],vertices[6]);// Lat plane
    G4double fSurface4=GetFaceSurfaceArea(vertices[2],vertices[1],
                                          vertices[5],vertices[6]);// Lat plane
    G4double fSurface5=GetFaceSurfaceArea(vertices[4],vertices[5],
                                          vertices[6],vertices[7]);// fDz plane 

    // Total Surface Area
    //
    if(!fIsTwisted)
    {
      fSurfaceArea = fSurface0+fSurface1+fSurface2
                   + fSurface3+fSurface4+fSurface5;
    }
    else
    {
      fSurfaceArea = G4VSolid::GetSurfaceArea();
    }
  }
  return fSurfaceArea;
}

G4double G4GenericTrap::GetTwistAngle

(

G4int

index

)

const [inline]

Definition at line 76 of file G4GenericTrap.icc.

References FatalException, and G4Exception().

Referenced by CreatePolyhedron(), and SurfaceNormal().

{
  if ( (index<0) || (index >= G4int(fVertices.size())) )
  {
    G4Exception ("G4GenericTrap::GetTwistAngle()", "GeomSolids0003",
                 FatalException, "Index outside range.");
    return 0.;
  }
  return fTwist[index];
}

G4TwoVector G4GenericTrap::GetVertex

(

G4int

index

)

const [inline]

Definition at line 54 of file G4GenericTrap.icc.

References FatalException, and G4Exception().

{
  if ( index<0 || index >= G4int(fVertices.size()) )
  {
    G4Exception ("G4GenericTrap::GetVertex()", "GeomSolids0003",
                 FatalException, "Index outside range.");
    return G4TwoVector();
  }
  return fVertices[index];
}  

const std::vector< G4TwoVector > & G4GenericTrap::GetVertices

(

)

const [inline]

Definition at line 68 of file G4GenericTrap.icc.

Referenced by G4GDMLWriteSolids::GenTrapWrite().

{
  return fVertices;
}  

G4int G4GenericTrap::GetVisSubdivisions

(

)

const [inline]

Definition at line 114 of file G4GenericTrap.icc.

Referenced by CreatePolyhedron().

{
  return fVisSubdivisions;
}

G4double G4GenericTrap::GetZHalfLength

(

)

const [inline]

Definition at line 38 of file G4GenericTrap.icc.

Referenced by G4GDMLWriteSolids::GenTrapWrite().

{
  return fDz;
}

EInside G4GenericTrap::Inside

(

const G4ThreeVector &

p

)

const [virtual]

Implements G4VSolid.

Definition at line 304 of file G4GenericTrap.cc.

References G4TessellatedSolid::Inside(), G4VSolid::kCarTolerance, kInside, kOutside, and kSurface.

Referenced by CalculateExtent(), and DistanceToIn().

{
  // Test if point is inside this shape

#ifdef G4TESS_TEST
   if ( fTessellatedSolid )
   { 
     return fTessellatedSolid->Inside(p);
   }
#endif  

  static const G4double halfCarTolerance=kCarTolerance*0.5;
  EInside innew=kOutside;
  std::vector<G4TwoVector> xy;
 
  if (std::fabs(p.z()) <= fDz+halfCarTolerance)  // First check Z range
  {
    // Compute intersection between Z plane containing point and the shape
    //
    G4double cf = 0.5*(fDz-p.z())/fDz;
    for (G4int i=0; i<4; i++)
    {
      xy.push_back(fVertices[i+4]+cf*( fVertices[i]-fVertices[i+4]));
    }

    innew=InsidePolygone(p,xy);

    if( (innew==kInside) || (innew==kSurface) )
    { 
      if(std::fabs(p.z()) > fDz-halfCarTolerance)  { innew=kSurface; }
    }
  }
  return innew;    
} 

G4bool G4GenericTrap::IsTwisted

(

)

const [inline]

Definition at line 90 of file G4GenericTrap.icc.

{
  return fIsTwisted;
}

G4GenericTrap & G4GenericTrap::operator=

(

const G4GenericTrap &

rhs

)

Definition at line 205 of file G4GenericTrap.cc.

References fCubicVolume, fDz, fIsTwisted, fMaxBBoxVector, fMinBBoxVector, fpPolyhedron, fSurfaceArea, fTessellatedSolid, fTwist, fVertices, fVisSubdivisions, and G4VSolid::operator=().

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

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

   // Copy data
   //
   fpPolyhedron = 0; fDz = rhs.fDz; fVertices = rhs.fVertices;
   fIsTwisted = rhs.fIsTwisted; fTessellatedSolid = 0;
   fMinBBoxVector = rhs.fMinBBoxVector; fMaxBBoxVector = rhs.fMaxBBoxVector;
   fVisSubdivisions = rhs.fVisSubdivisions;
   fSurfaceArea = rhs.fSurfaceArea; fCubicVolume = rhs.fCubicVolume;

   for (size_t i=0; i<4; ++i)  { fTwist[i] = rhs.fTwist[i]; }
#ifdef G4TESS_TEST
   if (rhs.fTessellatedSolid && !fIsTwisted )
   { delete fTessellatedSolid; fTessellatedSolid = CreateTessellatedSolid(); } 
#endif

   return *this;
}

void G4GenericTrap::SetVisSubdivisions

(

G4int

subdiv

)

[inline]

Definition at line 122 of file G4GenericTrap.icc.

{
  fVisSubdivisions=subdiv;
}

std::ostream & G4GenericTrap::StreamInfo

(

std::ostream &

os

)

const [virtual]

Implements G4VSolid.

Definition at line 1403 of file G4GenericTrap.cc.

References G4endl, GetEntityType(), and G4VSolid::GetName().

{
  G4int oldprc = os.precision(16);
  os << "-----------------------------------------------------------\n"
     << "    *** Dump for solid - " << GetName() << " *** \n"
     << "    =================================================== \n"
     << " Solid geometry type: " << GetEntityType()  << G4endl
     << "   half length Z: " << fDz/mm << " mm \n"
     << "   list of vertices:\n";
     
  for ( G4int i=0; i<fgkNofVertices; ++i )
  {
    os << std::setw(5) << "#" << i 
       << "   vx = " << fVertices[i].x()/mm << " mm" 
       << "   vy = " << fVertices[i].y()/mm << " mm" << G4endl;
  }
  os.precision(oldprc);

  return os;
} 

G4ThreeVector G4GenericTrap::SurfaceNormal

(

const G4ThreeVector &

p

)

const [virtual]

Implements G4VSolid.

Definition at line 341 of file G4GenericTrap.cc.

References G4Exception(), GetTwistAngle(), JustWarning, G4VSolid::kCarTolerance, and G4TessellatedSolid::SurfaceNormal().

{
  // Calculate side nearest to p, and return normal
  // If two sides are equidistant, sum of the Normal is returned

#ifdef G4TESS_TEST
  if ( fTessellatedSolid )
  {
    return fTessellatedSolid->SurfaceNormal(p);
  }  
#endif   

  static const G4double halfCarTolerance=kCarTolerance*0.5;
 
  G4ThreeVector lnorm, sumnorm(0.,0.,0.), apprnorm(0.,0.,1.),
                p0, p1, p2, r1, r2, r3, r4;
  G4int noSurfaces = 0; 
  G4double distxy,distz;
  G4bool zPlusSide=false;

  distz = fDz-std::fabs(p.z());
  if (distz < halfCarTolerance)
  {
    if(p.z()>0)
    {
      zPlusSide=true;
      sumnorm=G4ThreeVector(0,0,1);
    }
    else
    {
      sumnorm=G4ThreeVector(0,0,-1);
    }
    noSurfaces ++;
  } 

  // Check lateral planes
  //
  std:: vector<G4TwoVector> vertices;  
  G4double cf = 0.5*(fDz-p.z())/fDz;
  for (G4int i=0; i<4; i++)
  {
    vertices.push_back(fVertices[i+4]+cf*(fVertices[i]-fVertices[i+4]));
  }

  // Compute distance for lateral planes
  //
  for (G4int q=0; q<4; q++)
  {
    p0=G4ThreeVector(vertices[q].x(),vertices[q].y(),p.z());
    if(zPlusSide)
    {
      p1=G4ThreeVector(fVertices[q].x(),fVertices[q].y(),-fDz);
    }
    else
    {
      p1=G4ThreeVector(fVertices[q+4].x(),fVertices[q+4].y(),fDz); 
    }
    p2=G4ThreeVector(vertices[(q+1)%4].x(),vertices[(q+1)%4].y(),p.z());

    // Collapsed vertices
    //
    if ( (p2-p0).mag2() < kCarTolerance )
    {
      if ( std::fabs(p.z()+fDz) > kCarTolerance )
      {
        p2=G4ThreeVector(fVertices[(q+1)%4].x(),fVertices[(q+1)%4].y(),-fDz);
      }
      else
      {
        p2=G4ThreeVector(fVertices[(q+1)%4+4].x(),fVertices[(q+1)%4+4].y(),fDz);
      }
    }
    lnorm = (p1-p0).cross(p2-p0);
    lnorm = lnorm.unit();
    if(zPlusSide)  { lnorm=-lnorm; }

    // Adjust Normal for Twisted Surface
    //
    if ( (fIsTwisted) && (GetTwistAngle(q)!=0) )
    {
      G4double normP=(p2-p0).mag();
      if(normP)
      {
        G4double proj=(p-p0).dot(p2-p0)/normP;
        if(proj<0)     { proj=0; }
        if(proj>normP) { proj=normP; }
        G4int j=(q+1)%4;
        r1=G4ThreeVector(fVertices[q+4].x(),fVertices[q+4].y(),fDz);
        r2=G4ThreeVector(fVertices[j+4].x(),fVertices[j+4].y(),fDz);
        r3=G4ThreeVector(fVertices[q].x(),fVertices[q].y(),-fDz);
        r4=G4ThreeVector(fVertices[j].x(),fVertices[j].y(),-fDz);
        r1=r1+proj*(r2-r1)/normP;
        r3=r3+proj*(r4-r3)/normP;
        r2=r1-r3;
        r4=r2.cross(p2-p0); r4=r4.unit();
        lnorm=r4;
      }
    }   // End if fIsTwisted

    distxy=std::fabs((p0-p).dot(lnorm));
    if ( distxy<halfCarTolerance )
    {
      noSurfaces ++;

      // Negative sign for Normal is taken for Outside Normal
      //
      sumnorm=sumnorm+lnorm;
    }

    // For ApproxSurfaceNormal
    //
    if (distxy<distz)
    {
      distz=distxy;
      apprnorm=lnorm;
    }      
  }  // End for loop

  // Calculate final Normal, add Normal in the Corners and Touching Sides
  //
  if ( noSurfaces == 0 )
  {
    G4Exception("G4GenericTrap::SurfaceNormal(p)", "GeomSolids1002",
                JustWarning, "Point p is not on surface !?" );
    sumnorm=apprnorm;
    // Add Approximative Surface Normal Calculation?
  }
  else if ( noSurfaces == 1 )  { sumnorm = sumnorm; }
  else                         { sumnorm = sumnorm.unit(); }

  return sumnorm ; 
}    

---

Field Documentation

G4Polyhedron* G4GenericTrap::fpPolyhedron [mutable, protected]

Definition at line 192 of file G4GenericTrap.hh.

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

---

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

• G4GenericTrap.hh
• G4GenericTrap.icc
• G4GenericTrap.cc

---