G4GenericTrap Class Reference

G4GenericTrap is a solid which represents an arbitrary trapezoid with up to 8 vertices standing on two parallel planes perpendicular to the Z axis. Points can be identical in order to create shapes with less than 8 vertices. More...

#include <G4GenericTrap.hh>

Inheritance diagram for G4GenericTrap:

Public Member Functions
	G4GenericTrap (const G4String &name, G4double halfZ, const std::vector< G4TwoVector > &vertices)
	G4GenericTrap (__void__ &)
	G4GenericTrap (const G4GenericTrap &rhs)
G4GenericTrap &	operator= (const G4GenericTrap &rhs)
	~G4GenericTrap () override=default
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 override
G4ThreeVector	SurfaceNormal (const G4ThreeVector &p) const override
G4double	DistanceToIn (const G4ThreeVector &p, const G4ThreeVector &v) const override
G4double	DistanceToIn (const G4ThreeVector &p) const override
G4double	DistanceToOut (const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=false, G4bool *validNorm=nullptr, G4ThreeVector *n=nullptr) const override
G4double	DistanceToOut (const G4ThreeVector &p) const override
void	BoundingLimits (G4ThreeVector &pMin, G4ThreeVector &pMax) const override
G4bool	CalculateExtent (const EAxis pAxis, const G4VoxelLimits &pVoxelLimit, const G4AffineTransform &pTransform, G4double &pmin, G4double &pmax) const override
G4GeometryType	GetEntityType () const override
G4bool	IsFaceted () const override
G4VSolid *	Clone () const override
std::ostream &	StreamInfo (std::ostream &os) const override
G4ThreeVector	GetPointOnSurface () const override
G4double	GetCubicVolume () override
G4double	GetSurfaceArea () override
void	DescribeYourselfTo (G4VGraphicsScene &scene) const override
G4VisExtent	GetExtent () const override
G4Polyhedron *	CreatePolyhedron () const override
G4Polyhedron *	GetPolyhedron () const override
Public Member Functions inherited from G4VSolid
	G4VSolid (const G4String &name)
virtual	~G4VSolid ()
	G4VSolid (const G4VSolid &rhs)
G4VSolid &	operator= (const G4VSolid &rhs)
G4bool	operator== (const G4VSolid &s) const
G4String	GetName () const
void	SetName (const G4String &name)
G4double	GetTolerance () const
virtual void	ComputeDimensions (G4VPVParameterisation *p, const G4int n, const G4VPhysicalVolume *pRep)
virtual G4int	GetNumOfConstituents () const
void	DumpInfo () const
virtual const G4VSolid *	GetConstituentSolid (G4int no) const
virtual G4VSolid *	GetConstituentSolid (G4int no)
virtual const G4DisplacedSolid *	GetDisplacedSolidPtr () const
virtual G4DisplacedSolid *	GetDisplacedSolidPtr ()
	G4VSolid (__void__ &)
G4double	EstimateCubicVolume (G4int nStat, G4double epsilon) const
G4double	EstimateSurfaceArea (G4int nStat, G4double epsilon) const

Additional Inherited Members
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 G4VSolid
G4double	kCarTolerance

Detailed Description

G4GenericTrap is a solid which represents an arbitrary trapezoid with up to 8 vertices standing on two parallel planes perpendicular to the Z axis. Points can be identical in order to create shapes with less than 8 vertices.

Definition at line 83 of file G4GenericTrap.hh.

Constructor & Destructor Documentation

G4GenericTrap() [1/3]

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

Constructs an generic trapezoid, given its vertices.

Parameters

[in]	name	The solid name.
[in]	halfZ	Half length in Z.
[in]	vertices	The (x,y) coordinates of the vertices.

Definition at line 67 of file G4GenericTrap.cc.

  : G4VSolid(name)
{
  halfTolerance = 0.5*kCarTolerance;
  CheckParameters(halfZ, vertices);
  ComputeLateralSurfaces();
  ComputeBoundingBox();
  ComputeScratchLength();
}

Referenced by Clone(), G4GenericTrap(), and operator=().

G4GenericTrap() [2/3]

G4GenericTrap::G4GenericTrap

(

__void__ &

a

)

Fake default constructor for usage restricted to direct object persistency for clients requiring preallocation of memory for persistifiable objects.

Definition at line 82 of file G4GenericTrap.cc.

  : G4VSolid(a)
{
}

G4GenericTrap() [3/3]

G4GenericTrap::G4GenericTrap

(

const G4GenericTrap &

rhs

)

Copy constructor and assignment operator.

Definition at line 91 of file G4GenericTrap.cc.

  : G4VSolid(rhs),
    halfTolerance(rhs.halfTolerance), fScratch(rhs.fScratch),
    fDz(rhs.fDz), fVertices(rhs.fVertices), fIsTwisted(rhs.fIsTwisted),
    fMinBBox(rhs.fMinBBox), fMaxBBox(rhs.fMaxBBox),
    fVisSubdivisions(rhs.fVisSubdivisions),
    fSurfaceArea(rhs.fSurfaceArea), fCubicVolume(rhs.fCubicVolume)
{
  for (auto i = 0; i < 5; ++i) { fTwist[i] = rhs.fTwist[i]; }
  for (auto i = 0; i < 4; ++i) { fDelta[i] = rhs.fDelta[i]; }
  for (auto i = 0; i < 8; ++i) { fPlane[i] = rhs.fPlane[i]; }
  for (auto i = 0; i < 4; ++i) { fSurf[i] = rhs.fSurf[i]; }
  for (auto i = 0; i < 4; ++i) { f4k[i] = rhs.f4k[i]; }
  for (auto i = 0; i < 4; ++i) { fArea[i] = rhs.fArea[i]; }
}

~G4GenericTrap()

G4GenericTrap::~G4GenericTrap ( )

overridedefault

Default Destructor.

Member Function Documentation

BoundingLimits()

void G4GenericTrap::BoundingLimits ( G4ThreeVector & pMin, G4ThreeVector & pMax ) const

overridevirtual

Computes the bounding limits of the solid.

Parameters

[out]	pMin	The minimum bounding limit point.
[out]	pMax	The maximum bounding limit point.

Reimplemented from G4VSolid.

Definition at line 870 of file G4GenericTrap.cc.

{
  pMin = fMinBBox;
  pMax = fMaxBBox;
}

Referenced by CalculateExtent().

CalculateExtent()

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

overridevirtual

Calculates the minimum and maximum extent of the solid, when under the specified transform, and within the specified limits.

Parameters

[in]	pAxis	The axis along which compute the extent.
[in]	pVoxelLimit	The limiting space dictated by voxels.
[in]	pTransform	The internal transformation applied to the solid.
[out]	pMin	The minimum extent value.
[out]	pMax	The maximum extent value.

Returns

True if the solid is intersected by the extent region.

Implements G4VSolid.

Definition at line 882 of file G4GenericTrap.cc.

{
  G4ThreeVector bmin, bmax;
  G4bool exist;
 
  // Check bounding box (bbox)
  //
  BoundingLimits(bmin,bmax);
  G4BoundingEnvelope bbox(bmin,bmax);
#ifdef G4BBOX_EXTENT
  return bbox.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
#endif
  if (bbox.BoundingBoxVsVoxelLimits(pAxis,pVoxelLimit,pTransform,pMin,pMax))
  {
    return exist = pMin < pMax;
  }
 
  // Set bounding envelope (benv) and calculate extent
  //
  // To build the bounding envelope with plane faces, each lateral face of
  // the trapezoid is subdivided in two triangles. Subdivision is done by
  // duplication of vertices in the bases in a way that the envelope be
  // a convex polyhedron (some faces of the envelope can be degenerated)
  //
  G4double dz = GetZHalfLength();
  G4ThreeVectorList baseA(8), baseB(8);
  for (G4int i = 0; i < 4; ++i)
  {
    G4TwoVector va = GetVertex(i);
    G4TwoVector vb = GetVertex(i+4);
    baseA[2*i].set(va.x(), va.y(),-dz);
    baseB[2*i].set(vb.x(), vb.y(), dz);
  }
  for (G4int i = 0; i < 4; ++i)
  {
    G4int k1 = 2*i, k2 = (2*i + 2)%8;
    G4double ax = (baseA[k2].x() - baseA[k1].x());
    G4double ay = (baseA[k2].y() - baseA[k1].y());
    G4double bx = (baseB[k2].x() - baseB[k1].x());
    G4double by = (baseB[k2].y() - baseB[k1].y());
    G4double znorm = ax*by - ay*bx;
    baseA[k1+1] = (znorm < 0.0) ? baseA[k2] : baseA[k1];
    baseB[k1+1] = (znorm < 0.0) ? baseB[k1] : baseB[k2];
  }
 
  std::vector<const G4ThreeVectorList *> polygons(2);
  polygons[0] = &baseA;
  polygons[1] = &baseB;
 
  G4BoundingEnvelope benv(bmin, bmax, polygons);
  exist = benv.CalculateExtent(pAxis,pVoxelLimit,pTransform,pMin,pMax);
  return exist;
}

Clone()

G4VSolid * G4GenericTrap::Clone ( ) const

overridevirtual

Makes a clone of the object for use in multi-treading.

Returns

A pointer to the new cloned allocated solid.

Reimplemented from G4VSolid.

Definition at line 961 of file G4GenericTrap.cc.

{
  return new G4GenericTrap(*this);
}

CreatePolyhedron()

G4Polyhedron * G4GenericTrap::CreatePolyhedron ( ) const

overridevirtual

Creates a Polyhedron used for Visualisation. It is the caller's responsibility to delete it. A null pointer means "not created".

Reimplemented from G4VSolid.

Definition at line 1531 of file G4GenericTrap.cc.

{
  // Approximation of Twisted Side
  // Construct extra Points, if Twisted Side
  //
  G4Polyhedron* polyhedron;
  G4int nVertices, nFacets;
 
  G4int subdivisions = 0;
  if (fIsTwisted)
  {
    if (GetVisSubdivisions() != 0)
    {
      subdivisions = GetVisSubdivisions();
    }
    else
    {
      // Estimation of Number of Subdivisions for smooth visualisation
      //
      G4double maxTwist = 0.;
      for(G4int i = 0; i < 4; ++i)
      {
        if (GetTwistAngle(i) > maxTwist) { maxTwist = GetTwistAngle(i); }
      }
 
      // Computes bounding vectors for the shape
      //
      G4double Dx, Dy;
      Dx = 0.5*(fMaxBBox.x() - fMinBBox.y());
      Dy = 0.5*(fMaxBBox.y() - fMinBBox.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 G4Polyhedron(nVertices, nFacets);
 
  // Set vertices
  //
  G4int icur = 0;
  for (G4int i = 0; i < 4; ++i)
  {
    G4ThreeVector v(fVertices[i].x(),fVertices[i].y(),-fDz);
    polyhedron->SetVertex(++icur, v);
  }
  for (G4int i = 0; i < subdivisions; ++i)
  {
    for (G4int j = 0; j < 4; ++j)
    {
      G4TwoVector u = fVertices[j]+cf*(i+1)*(fVertices[j+4]-fVertices[j]);
      G4ThreeVector v(u.x(),u.y(),-fDz+cf*2*fDz*(i+1));
      polyhedron->SetVertex(++icur, v);
    }
  }
  for (G4int i = 4; i < 8; ++i)
  {
    G4ThreeVector v(fVertices[i].x(),fVertices[i].y(),fDz);
    polyhedron->SetVertex(++icur, v);
  }
 
  // Set facets
  //
  icur = 0;
  polyhedron->SetFacet(++icur, 1, 4, 3, 2); // Z-plane
  for (G4int i = 0; i < subdivisions + 1; ++i)
  {
    G4int is = i*4;
    polyhedron->SetFacet(++icur, 5+is, 8+is, 4+is, 1+is);
    polyhedron->SetFacet(++icur, 8+is, 7+is, 3+is, 4+is);
    polyhedron->SetFacet(++icur, 7+is, 6+is, 2+is, 3+is);
    polyhedron->SetFacet(++icur, 6+is, 5+is, 1+is, 2+is);
  }
  polyhedron->SetFacet(++icur, 5+sub4, 6+sub4, 7+sub4, 8+sub4); // Z-plane
 
  polyhedron->SetReferences();
  polyhedron->InvertFacets();
 
  return polyhedron;
}

Referenced by GetPolyhedron().

DescribeYourselfTo()

void G4GenericTrap::DescribeYourselfTo ( G4VGraphicsScene & scene ) const

overridevirtual

Methods for creating graphical representations (i.e. for visualisation).

Implements G4VSolid.

Definition at line 1513 of file G4GenericTrap.cc.

{
  scene.AddSolid(*this);
}

DistanceToIn() [1/2]

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

overridevirtual

Calculates the distance to the nearest surface of a shape from an outside point. The distance can be an underestimate.

Parameters

[in]

p

The point at offset p.

Returns

The safety distance to enter the shape.

Implements G4VSolid.

Definition at line 570 of file G4GenericTrap.cc.

{
  G4double x = p.x(), y = p.y(), z = p.z();
  G4double dist = std::abs(z) - fDz;
  for (auto i = 0; i < 4; ++i)
  {
    if (fTwist[i] == 0.)
    {
      G4double distPlane = fSurf[i].D*x + fSurf[i].E*y + fSurf[i].F*z + fSurf[i].G;
      dist = std::max(distPlane, dist);
    }
    else
    {
      G4double A = fSurf[i].A;
      G4double B = fSurf[i].B;
      G4double C = fSurf[i].C;
      G4double Cz = C*z;
      G4double CzF = Cz + fSurf[i].F;
      G4double gradx = A*z + fSurf[i].D;
      G4double grady = B*z + fSurf[i].E;
      G4double gradz = A*x + B*y + Cz + CzF;
      G4double fun = gradx*x + grady*y + CzF*z + fSurf[i].G;
      G4double grad2 = gradx*gradx + grady*grady + gradz*gradz;
      G4double divisor = std::sqrt(grad2) + std::sqrt(grad2 + f4k[i]*std::abs(fun));
      G4double distSurf = 2.*fun/divisor;
      dist = std::max(distSurf, dist);
    }
  }
  return (dist > 0.) ? dist : 0.; // return safety distance
}

DistanceToIn() [2/2]

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

overridevirtual

Returns the distance along the normalised vector 'v' to the shape, from the point at offset 'p'. If there is no intersection, returns kInfinity. The first intersection resulting from 'leaving' a surface/volume is discarded. Hence, it is tolerant of points on the surface of the shape.

Parameters

[in]	p	The point at offset p.
[in]	v	The normalised direction vector.

Returns

The distance to enter the shape.

Implements G4VSolid.

Definition at line 315 of file G4GenericTrap.cc.

{
  G4double px = p.x(), py = p.y(), pz = p.z();
  G4double vx = v.x(), vy = v.y(), vz = v.z();
 
  // Find intersections with the bounding polyhedron
  //
  if (std::abs(pz) - fDz >= -halfTolerance && pz*vz >= 0) { return kInfinity; }
  G4double invz = (vz == 0) ? kInfinity : -1./vz;
  G4double dz = std::copysign(fDz,invz);
  G4double xin  = (pz - dz)*invz;
  G4double xout = (pz + dz)*invz;
 
  // Check plane faces
  for (auto k = 0; k < 4; ++k)
  {
    if (fTwist[k] != 0) { continue; } // skip twisted faces
    const G4GenericTrapPlane& surf = fPlane[2*k];
    G4double cosa = surf.A*vx + surf.B*vy + surf.C*vz;
    G4double dist = surf.A*px + surf.B*py + surf.C*pz + surf.D;
    if (dist >= -halfTolerance)
    {
      if (cosa >= 0.) { return kInfinity; } // point flies away
      G4double tmp  = -dist/cosa;
      xin = std::max(tmp, xin);
    }
    else
    {
      if (cosa > 0) { xout = std::min(-dist/cosa, xout); }
      if (cosa < 0) { xin = std::max(-dist/cosa, xin); }
    }
  }
 
  // Check planes around twisted faces, each twisted face is bounded by two planes
  G4double tin  = xin;
  G4double tout = xout;
  for (auto i = 0; i < 4; ++i)
  {
    if (fTwist[i] == 0) { continue; } // skip plane faces
 
    // check intersection with 1st bounding plane
    const G4GenericTrapPlane& surf1 = fPlane[2*i];
    G4double cosa = surf1.A*vx + surf1.B*vy + surf1.C*vz;
    G4double dist = surf1.A*px + surf1.B*py + surf1.C*pz + surf1.D;
    if (dist >= -halfTolerance)
    {
      if (cosa >= 0.) { return kInfinity; } // point flies away
      G4double tmp  = -dist/cosa;
      tin = std::max(tmp, tin);
    }
    else
    {
      if (cosa > 0) { tout = std::min(-dist/cosa, tout); }
      if (cosa < 0) { tin = std::max(-dist/cosa, tin); }
    }
 
    // check intersection with 2nd bounding plane
    const G4GenericTrapPlane& surf2 = fPlane[2*i + 1];
    cosa = surf2.A*vx + surf2.B*vy + surf2.C*vz;
    dist = surf2.A*px + surf2.B*py + surf2.C*pz + surf2.D;
    if (dist >= -halfTolerance)
    {
      if (cosa >= 0.) { return kInfinity; } // point flies away
      G4double tmp  = -dist/cosa;
      tin = std::max(tmp, tin);
    }
    else
    {
      if (cosa > 0) { tout = std::min(-dist/cosa, tout); }
      if (cosa < 0) { tin = std::max(-dist/cosa, tin); }
    }
  }
  if (tout - tin <= halfTolerance) { return kInfinity; } // touch or no hit
 
  // if point is outside of the bounding box
  // then move it to the surface of the bounding polyhedron
  G4double t0 = 0., x0 = px, y0 = py, z0 = pz;
  if (x0 < fMinBBox.x() - halfTolerance ||
      y0 < fMinBBox.y() - halfTolerance ||
      z0 < fMinBBox.z() - halfTolerance ||
      x0 > fMaxBBox.x() + halfTolerance ||
      y0 > fMaxBBox.y() + halfTolerance ||
      z0 > fMaxBBox.z() + halfTolerance)
  {
    t0 = tin;
    x0 += vx*t0;
    y0 += vy*t0;
    z0 += vz*t0;
  }
 
  // Check intersections with twisted faces
  //
  G4double ttin[2] = { DBL_MAX, DBL_MAX };
  G4double ttout[2] = { tout, tout };
 
  if (tin - xin < halfTolerance) { ttin[0] = xin; }
  if (xout - tout < halfTolerance) { ttout[0] = xout; }
  G4double tminimal = tin - halfTolerance;
  G4double tmaximal = tout + halfTolerance;
 
  constexpr G4double EPSILON = 1000.*DBL_EPSILON;
  for (auto i = 0; i < 4; ++i)
  {
    if (fTwist[i] == 0) { continue; } // skip plane faces
 
    // twisted face, solve quadratic equation
    G4double ABC  = fSurf[i].A*vx + fSurf[i].B*vy + fSurf[i].C*vz;
    G4double ABCF = fSurf[i].A*x0 + fSurf[i].B*y0 + fSurf[i].C*z0 + fSurf[i].F;
    G4double A = ABC*vz;
    G4double B = 0.5*(fSurf[i].D*vx + fSurf[i].E*vy + ABCF*vz + ABC*z0);
    G4double C = fSurf[i].D*x0 + fSurf[i].E*y0 + ABCF*z0 + fSurf[i].G;
    if (std::abs(A) <= EPSILON)
    {
      // case 1: track is parallel to the surface
      if (std::abs(B) <= EPSILON)
      {
        // check position of the track relative to the surface,
        // for the reason of precision it's better to use (x0,y0,z0) instead of (px,py,pz)
        auto j = (i + 1)%4;
        G4double z = (z0 + fDz);
        G4TwoVector a = fVertices[i] + fDelta[i]*z;
        G4TwoVector b = fVertices[j] + fDelta[j]*z;
        G4double dx = b.x() - a.x();
        G4double dy = b.y() - a.y();
        G4double leng = std::sqrt(dx*dx + dy*dy);
        G4double dist = dx*(y0 - a.y()) - dy*(x0 - a.x());
        if (dist >= -halfTolerance*leng) { return kInfinity; }
        continue;
      }
 
      // case 2: single root
      G4double tmp = t0 - 0.5*C/B;
      // compute normal at the intersection point and check direction
      G4double x = px + vx*tmp;
      G4double y = py + vy*tmp;
      G4double z = pz + vz*tmp;
      const G4GenericTrapSurface& surf = fSurf[i];
      G4double nx = surf.A*z + surf.D;
      G4double ny = surf.B*z + surf.E;
      G4double nz = surf.A*x + surf.B*y + 2.*surf.C*z + surf.F;
 
      if (nx*vx + ny*vy + nz*vz >= 0.) // point is flying to outside
      {
        auto k = (i == 3) ? 0 : i + 1;
        G4double tz = (pz + fDz);
        G4TwoVector a = fVertices[i] + fDelta[i]*tz;
        G4TwoVector b = fVertices[k] + fDelta[k]*tz;
        G4double dx = b.x() - a.x();
        G4double dy = b.y() - a.y();
        G4double leng = std::sqrt(dx*dx + dy*dy);
        G4double dist = dx*(py - a.y()) - dy*(px - a.x());
        if (dist >= -halfTolerance*leng) { return kInfinity; } // point is flies away
 
        if (tmp < tminimal || tmp > tmaximal) { continue; }
        if (std::abs(tmp - ttout[0]) < halfTolerance) { continue; }
        if (tmp < ttout[0])
        {
          ttout[1] = ttout[0];
          ttout[0] = tmp;
        }
        else { ttout[1] = std::min(tmp, ttout[1]); }
      }
      else // point is flying to inside
      {
        if (tmp < tminimal || tmp > tmaximal) { continue; }
        if (std::abs(tmp - ttin[0]) < halfTolerance) { continue; }
        if (tmp < ttin[0])
        {
          ttin[1] = ttin[0];
          ttin[0] = tmp;
        }
        else { ttin[1] = std::min(tmp, ttin[1]); }
      }
      continue;
    }
 
    // case 3: scratch or no intersection
    G4double D = B*B - A*C;
    if (D < 0.25*fScratch*fScratch*A*A)
    {
      if (A > 0) { return kInfinity; }
      continue;
    }
 
    // case 4: two intersection points
    G4double tmp = -B - std::copysign(std::sqrt(D), B);
    G4double t1 = tmp/A + t0;
    G4double t2 = C/tmp + t0;
    G4double tsurfin = std::min(t1, t2);
    G4double tsurfout = std::max(t1, t2);
    if (A < 0) { std::swap(tsurfin, tsurfout); }
 
    if (tsurfin >= tminimal && tsurfin <= tmaximal)
    {
      if (std::abs(tsurfin - ttin[0]) >= halfTolerance)
      {
        if (tsurfin < ttin[0])
        {
          ttin[1] = ttin[0];
          ttin[0] = tsurfin;
        }
        else { ttin[1] = std::min(tsurfin, ttin[1]); }
      }
    }
    if (tsurfout >= tminimal && tsurfout <= tmaximal)
    {
      if (std::abs(tsurfout - ttout[0]) >= halfTolerance)
      {
        if (tsurfout < ttout[0])
        {
          ttout[1] = ttout[0];
          ttout[0] = tsurfout;
        }
        else { ttout[1] = std::min(tsurfout, ttout[1]); }
      }
    }
  }
 
  // Compute distance to In
  //
  if (ttin[0] == DBL_MAX) { return kInfinity; } // no entry point
 
  // single entry point
  if (ttin[1] == DBL_MAX)
  {
    G4double distin = ttin[0];
    G4double distout = (distin >= ttout[0] - halfTolerance) ? ttout[1] : ttout[0];
    G4double dist = (distout <= halfTolerance || distout - distin <= halfTolerance) ? kInfinity : distin;
    return (dist < halfTolerance) ? 0. : dist;
  }
 
  // two entry points
  if (ttin[1] < ttout[0])
  {
    G4double distin = ttin[1], distout = ttout[0];
    G4double dist = (distout <= halfTolerance || distout - distin <= halfTolerance) ? kInfinity : distin;
    return (dist < halfTolerance) ? 0. : dist;
  }
 
  // check 1st pair of in-out points
  G4double distin1 = ttin[0], distout1 = ttout[0];
  G4double dist1 = (distout1 <= halfTolerance || distout1 - distin1 <= halfTolerance) ? kInfinity : distin1;
  if (dist1 != kInfinity) { return (dist1 < halfTolerance) ? 0. : dist1; }
 
  // check 2nd pair of in-out points
  G4double distin2 = ttin[1], distout2 = ttout[1];
  G4double dist2 = (distout2 <= halfTolerance || distout2 - distin2 <= halfTolerance) ? kInfinity : distin2;
  return (dist2 < halfTolerance) ? 0. : dist2;
}

DistanceToOut() [1/2]

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

overridevirtual

Calculates the distance to the nearest surface of a shape from an inside point 'p'. The distance can be an underestimate.

Parameters

[in]

p

The point at offset p.

Returns

The safety distance to exit the shape.

Implements G4VSolid.

Definition at line 835 of file G4GenericTrap.cc.

{
  G4double x = p.x(), y = p.y(), z = p.z();
  G4double dist = std::abs(z) - fDz;
  for (auto i = 0; i < 4; ++i)
  {
    if (fTwist[i] == 0.)
    {
      G4double distPlane = fSurf[i].D*x + fSurf[i].E*y + fSurf[i].F*z + fSurf[i].G;
      dist = std::max(distPlane, dist);
    }
    else
    {
      G4double A = fSurf[i].A;
      G4double B = fSurf[i].B;
      G4double C = fSurf[i].C;
      G4double Cz = C*z;
      G4double CzF = Cz + fSurf[i].F;
      G4double gradx = A*z + fSurf[i].D;
      G4double grady = B*z + fSurf[i].E;
      G4double gradz = A*x + B*y + Cz + CzF;
      G4double fun = gradx*x + grady*y + CzF*z + fSurf[i].G;
      G4double grad2 = gradx*gradx + grady*grady + gradz*gradz;
      G4double divisor = std::sqrt(grad2) + std::sqrt(grad2 + f4k[i]*std::abs(fun));
      G4double distSurf = 2.*fun/divisor;
      dist = std::max(distSurf, dist);
    }
  }
  return (dist < 0.) ? -dist : 0.; // return safety distance
}

DistanceToOut() [2/2]

G4double G4GenericTrap::DistanceToOut ( const G4ThreeVector & p, const G4ThreeVector & v, const G4bool calcNorm = false, G4bool * validNorm = nullptr, G4ThreeVector * n = nullptr ) const

overridevirtual

Returns the distance along the normalised vector 'v' to the shape, from a point at an offset 'p' inside or on the surface of the shape. Intersections with surfaces, when the point is less than Tolerance/2 from a surface must be ignored.

Parameters

[in]	p	The point at offset p.
[in]	v	The normalised direction vector.
[in]	calcNorm	Flag to indicate if to calculate the normal or not.
[out]	validNorm	Flag set to true if the solid lies entirely behind or on the exiting surface. It is set false if the solid does not lie entirely behind or on the exiting surface. 'calcNorm' must be true, otherwise it is unused.
[out]	n	The exiting outwards normal vector (undefined Magnitude). 'calcNorm' must be true, otherwise it is unused.

Returns

The distance to exit the shape.

Implements G4VSolid.

Definition at line 605 of file G4GenericTrap.cc.

{
  G4double px = p.x(), py = p.y(), pz = p.z();
  G4double vx = v.x(), vy = v.y(), vz = v.z();
 
  // Check intersections with plane faces
  //
  if ((std::abs(pz) - fDz) >= -halfTolerance && pz*vz > 0.)
  {
    if (calcNorm)
    {
      *validNorm = true;
      n->set(0., 0., std::copysign(1., pz));
    }
    return 0.;
  }
  G4double tout = (vz == 0) ? DBL_MAX : (std::copysign(fDz, vz) - pz)/vz;
  G4int iface = (vz < 0) ? -4 : -2; // little trick for z-normal: (-4+3)=-1, (-2+3)=+1
 
  for (auto i = 0; i < 4; ++i)
  {
    if (fTwist[i] != 0) { continue; } // skip twisted faces
    const G4GenericTrapPlane& surf = fPlane[2*i];
    G4double cosa = surf.A*vx + surf.B*vy + surf.C*vz;
    if (cosa > 0)
    {
      G4double dist = surf.A*px + surf.B*py + surf.C*pz + surf.D;
      if (dist >= -halfTolerance)
      {
        if (calcNorm)
        {
           *validNorm = true;
           n->set(surf.A, surf.B, surf.C);
        }
        return 0.;
      }
      G4double tmp = -dist/cosa;
      if (tout > tmp) { tout = tmp; iface = i; }
    }
  }
 
  // Check intersections with twisted faces
  //
  constexpr G4double EPSILON = 1000.*DBL_EPSILON;
  for (auto i = 0; i < 4; ++i)
  {
    if (fTwist[i] == 0) { continue; } // skip plane faces
 
    // twisted face, solve quadratic equation
    const G4GenericTrapSurface& surf = fSurf[i];
    G4double ABC  = surf.A*vx + surf.B*vy + surf.C*vz;
    G4double ABCF = surf.A*px + surf.B*py + surf.C*pz + surf.F;
    G4double A = ABC*vz;
    G4double B = 0.5*(surf.D*vx + surf.E*vy + ABCF*vz + ABC*pz);
    G4double C = surf.D*px + surf.E*py + ABCF*pz + surf.G;
 
    if (std::abs(A) <= EPSILON)
    {
      // case 1: track is parallel to the surface
      if (std::abs(B) <= EPSILON) { continue; }
 
      // case 2: single root
      G4double tmp = -0.5*C/B;
      // compute normal at intersection point and check direction
      G4double x = px + vx*tmp;
      G4double y = py + vy*tmp;
      G4double z = pz + vz*tmp;
      G4double nx = surf.A*z + surf.D;
      G4double ny = surf.B*z + surf.E;
      G4double nz = surf.A*x + surf.B*y + 2.*surf.C*z + surf.F;
 
      if (nx*vx + ny*vy + nz*vz > 0.) // point is flying to outside
      {
        auto k = (i + 1)%4;
        G4double tz = (pz + fDz);
        G4TwoVector a = fVertices[i] + fDelta[i]*tz;
        G4TwoVector b = fVertices[k] + fDelta[k]*tz;
        G4double dx = b.x() - a.x();
        G4double dy = b.y() - a.y();
        G4double leng = std::sqrt(dx*dx + dy*dy);
        G4double dist = dx*(py - a.y()) - dy*(px - a.x());
        if (dist >= -halfTolerance*leng) // point is on the surface
        {
          if (calcNorm)
          {
            *validNorm = false;
            G4double normx = surf.A*pz + surf.D;
            G4double normy = surf.B*pz + surf.E;
            G4double normz = surf.A*px + surf.B*py + 2.*surf.C*pz + surf.F;
            G4double inv = 1./std::sqrt(normx*normx + normy*normy + normz*normz);
            n->set(normx*inv, normy*inv, normz*inv);
          }
          return 0.;
        }
        if (tout > tmp) { tout = tmp; iface = i; }
      }
      continue;
    }
 
    // case 3: scratch or no intersection
    G4double D = B*B - A*C;
    if (D < 0.25*fScratch*fScratch*A*A)
    {
      // check position of the point
      auto j = (i + 1)%4;
      G4double tz = pz + fDz;
      G4TwoVector a = fVertices[i] + fDelta[i]*tz;
      G4TwoVector b = fVertices[j] + fDelta[j]*tz;
      G4double dx = b.x() - a.x();
      G4double dy = b.y() - a.y();
      G4double leng = std::sqrt(dx*dx + dy*dy);
      G4double dist = dx*(py - a.y()) - dy*(px - a.x());
      if  (dist <= -halfTolerance*leng) { continue; } // point is inside
      if  (A > 0 || dist > halfTolerance*leng) // convex surface (or point is outside)
      {
        if (calcNorm)
        {
          *validNorm = false;
          G4double nx = surf.A*pz + surf.D;
          G4double ny = surf.B*pz + surf.E;
          G4double nz = surf.A*px + surf.B*py + 2.*surf.C*pz + surf.F;
          G4double inv = 1./std::sqrt(nx*nx + ny*ny + nz*nz);
          n->set(nx*inv, ny*inv, nz*inv);
        }
        return 0.;
      }
      continue;
    }
 
    // case 4: two intersection points
    G4double tmp = -B - std::copysign(std::sqrt(D), B);
    G4double t1 = tmp / A;
    G4double t2 = C / tmp;
    G4double tmin = std::min(t1, t2);
    G4double tmax = std::max(t1, t2);
 
    if (A < 0) // concave profile: tmin(out) -> tmax(in)
    {
      if (std::abs(tmax) < std::abs(tmin)) { continue; } // point flies inside
      if (tmin <= halfTolerance) // point is on external side of the surface
      {
        G4double t = 0.5*(tmin + tmax);
        G4double x = px + vx*t;
        G4double y = py + vy*t;
        G4double z = pz + vz*t;
 
        auto j = (i + 1)%4;
        G4double tz = z + fDz;
        G4TwoVector a = fVertices[i] + fDelta[i]*tz;
        G4TwoVector b = fVertices[j] + fDelta[j]*tz;
        G4double dx = b.x() - a.x();
        G4double dy = b.y() - a.y();
        G4double leng = std::sqrt(dx*dx + dy*dy);
        G4double dist = dx*(y - a.y()) - dy*(x - a.x());
        if  (dist <= -halfTolerance*leng) { continue; } // scratch
        if (calcNorm)
        {
          *validNorm = false;
          G4double nx = surf.A*pz + surf.D;
          G4double ny = surf.B*pz + surf.E;
          G4double nz = surf.A*px + surf.B*py + 2.*surf.C*pz + surf.F;
          G4double inv = 1./std::sqrt(nx*nx + ny*ny + nz*nz);
          n->set(nx*inv, ny*inv, nz*inv);
        }
        return 0.;
      }
      if (tout > tmin) { tout = tmin; iface = i; }
    }
    else // convex profile: tmin(in) -> tmax(out)
    {
      if (tmax < halfTolerance) // point is on the surface
      {
        if (calcNorm)
        {
          *validNorm = false;
          G4double nx = surf.A*pz + surf.D;
          G4double ny = surf.B*pz + surf.E;
          G4double nz = surf.A*px + surf.B*py + 2.*surf.C*pz + surf.F;
          G4double inv = 1./std::sqrt(nx*nx + ny*ny + nz*nz);
          n->set(nx*inv, ny*inv, nz*inv);
        }
        return 0.;
      }
      if (tout > tmax) { tout = tmax; iface = i; }
    }
  }
 
  // Compute normal, if required, and return distance to out
  //
  if (tout < halfTolerance) { tout = 0.; }
  if (calcNorm)
  {
    if (iface < 0)
    {
      *validNorm = true;
      n->set(0, 0, iface + 3); // little trick: (-4+3)=-1, (-2+3)=+1
    }
    else
    {
      const G4GenericTrapSurface& surf = fSurf[iface];
      if (fTwist[iface] == 0)
      {
        *validNorm = true;
        n->set(surf.D, surf.E, surf.F);
      }
      else
      {
        *validNorm = false;
        G4double x = px + vx*tout;
        G4double y = py + vy*tout;
        G4double z = pz + vz*tout;
        G4double nx = surf.A*z + surf.D;
        G4double ny = surf.B*z + surf.E;
        G4double nz = surf.A*x + surf.B*y + 2.*surf.C*z + surf.F;
        G4double inv = 1./std::sqrt(nx*nx + ny*ny + nz*nz);
        n->set(nx*inv, ny*inv, nz*inv);
      }
    }
  }
  return tout;
}

GetCubicVolume()

G4double G4GenericTrap::GetCubicVolume ( )

overridevirtual

Returning an estimation of the solid volume (capacity) and surface area, in internal units.

Reimplemented from G4VSolid.

Definition at line 1136 of file G4GenericTrap.cc.

{
  if (fCubicVolume == 0)
  {
    G4AutoLock l(&gentrapMutex);
    // diagonals
    G4TwoVector A = fVertices[3] - fVertices[1];
    G4TwoVector B = fVertices[2] - fVertices[0];
    G4TwoVector C = fVertices[7] - fVertices[5];
    G4TwoVector D = fVertices[6] - fVertices[4];
 
    // kross products
    G4double AB = A.x()*B.y() - A.y()*B.x();
    G4double CD = C.x()*D.y() - C.y()*D.x();
    G4double AD = A.x()*D.y() - A.y()*D.x();
    G4double CB = C.x()*B.y() - C.y()*B.x();
 
    fCubicVolume = ((AB + CD)/3. + (AD + CB)/6.)*fDz;
    l.unlock();
  }
  return fCubicVolume;
}

GetEntityType()

G4GeometryType G4GenericTrap::GetEntityType ( ) const

overridevirtual

Returns the type ID, "G4GenericTrap" of the solid.

Implements G4VSolid.

Definition at line 943 of file G4GenericTrap.cc.

{
  return { "G4GenericTrap" };
}

Referenced by StreamInfo().

GetExtent()

G4VisExtent G4GenericTrap::GetExtent ( ) const

overridevirtual

Provides extent (bounding box) as possible hint to the graphics view.

Reimplemented from G4VSolid.

Definition at line 1522 of file G4GenericTrap.cc.

{
  return { fMinBBox.x(), fMaxBBox.x(),
           fMinBBox.y(), fMaxBBox.y(),
           fMinBBox.z(), fMaxBBox.z() };
}

GetNofVertices()

G4int G4GenericTrap::GetNofVertices ( ) const

inline

GetPointOnSurface()

G4ThreeVector G4GenericTrap::GetPointOnSurface ( ) const

overridevirtual

Returns a random point located and uniformly distributed on the surface of the solid.

Reimplemented from G4VSolid.

Definition at line 992 of file G4GenericTrap.cc.

{
  if (fArea[0] + fArea[1] + fArea[2] + fArea[3] == 0.)
  {
    G4AutoLock l(&lateralareaMutex);
    fArea[0] = GetLateralFaceArea(0);
    fArea[1] = GetLateralFaceArea(1);
    fArea[2] = GetLateralFaceArea(2);
    fArea[3] = GetLateralFaceArea(3);
    l.unlock();
  }
 
  constexpr G4int iface[6][4] =
    { {0,1,2,3}, {0,4,5,1}, {1,5,6,2}, {2,6,7,3}, {3,7,4,0}, {4,5,6,7} };
 
  G4bool isTwisted[6] = {false};
  for (auto i = 0; i < 4; ++i)
  {
    isTwisted[i + 1] = (fTwist[i] != 0.0);
  }
 
  G4double ssurf[6];
  G4TwoVector A = fVertices[3] - fVertices[1];
  G4TwoVector B = fVertices[2] - fVertices[0];
  G4TwoVector C = fVertices[7] - fVertices[5];
  G4TwoVector D = fVertices[6] - fVertices[4];
  ssurf[0] = (A.x()*B.y() - A.y()*B.x())*0.5; // -fDz face
  ssurf[1] = ssurf[0] + fArea[0];
  ssurf[2] = ssurf[1] + fArea[1];
  ssurf[3] = ssurf[2] + fArea[2];
  ssurf[4] = ssurf[3] + fArea[3];
  ssurf[5] = ssurf[4] + (C.x()*D.y() - C.y()*D.x())*.5; // +fDz face
 
  G4double select = ssurf[5]*G4QuickRand();
  G4int k;
  for (k = 0; k < 5; ++k)
  {
    if (select <= ssurf[k]) { break; }
  }
 
  G4int i0 = iface[k][0];
  G4int i1 = iface[k][1];
  G4int i2 = iface[k][2];
  G4int i3 = iface[k][3];
  G4ThreeVector pp[4];
  pp[0].set(fVertices[i0].x(), fVertices[i0].y(), ((k == 5) ?  fDz : -fDz));
  pp[1].set(fVertices[i1].x(), fVertices[i1].y(), ((k == 0) ? -fDz :  fDz));
  pp[2].set(fVertices[i2].x(), fVertices[i2].y(), ((k == 0) ? -fDz :  fDz));
  pp[3].set(fVertices[i3].x(), fVertices[i3].y(), ((k == 5) ?  fDz : -fDz));
 
  G4ThreeVector point;
  if (isTwisted[k])
  { // twisted face, rejection sampling
    G4double maxlength = std::max((pp[2] - pp[1]).mag(), (pp[3] - pp[0]).mag());
    point = (pp[0] + pp[1] + pp[2] + pp[3])*0.25;
    for (auto i = 0; i < 10000; ++i)
    {
      G4double u = G4QuickRand();
      G4ThreeVector p0 = pp[0] + (pp[1] - pp[0])*u;
      G4ThreeVector p1 = pp[3] + (pp[2] - pp[3])*u;
      G4double v = G4QuickRand()*(maxlength/(p1 - p0).mag());
      if (v >= 1.) { continue; }
      point = p0 + (p1 - p0)*v;
      break;
    }
  }
  else
  { // plane face
    G4double u = G4QuickRand();
    G4double v = G4QuickRand();
    if (u + v > 1.) { u = 1. - u; v = 1. - v; }
    G4double ss = (((pp[2] - pp[0]).cross(pp[3] - pp[0])).mag())*0.5;
    G4int j = (select > ssurf[k] - ss) ? 3 : 1;
    point = pp[0] + (pp[2] - pp[0])*u + (pp[j] - pp[0])*v;
  }
  return point;
}

GetPolyhedron()

G4Polyhedron * G4GenericTrap::GetPolyhedron ( ) const

overridevirtual

Smart access function - creates on request and stores for future access. A null pointer means "not available".

Reimplemented from G4VSolid.

Definition at line 1493 of file G4GenericTrap.cc.

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

GetSurfaceArea()

G4double G4GenericTrap::GetSurfaceArea ( )

overridevirtual

Returns an estimation of the solid surface area in internal units. This method may be overloaded by derived classes to compute the exact geometrical quantity for solids where this is possible, or anyway to cache the computed value. Note: the computed value is NOT cached.

Reimplemented from G4VSolid.

Definition at line 1163 of file G4GenericTrap.cc.

{
  if (fSurfaceArea == 0)
  {
    G4AutoLock l(&gentrapMutex);
    G4TwoVector A = fVertices[3] - fVertices[1];
    G4TwoVector B = fVertices[2] - fVertices[0];
    G4TwoVector C = fVertices[7] - fVertices[5];
    G4TwoVector D = fVertices[6] - fVertices[4];
    G4double S_bot = (A.x()*B.y() - A.y()*B.x())*0.5;
    G4double S_top = (C.x()*D.y() - C.y()*D.x())*0.5;
    fArea[0] = GetLateralFaceArea(0);
    fArea[1] = GetLateralFaceArea(1);
    fArea[2] = GetLateralFaceArea(2);
    fArea[3] = GetLateralFaceArea(3);
    fSurfaceArea = S_bot + S_top + fArea[0] + fArea[1] + fArea[2] + fArea[3];
    l.unlock();
  }
  return fSurfaceArea;
}

GetTwistAngle()

G4double G4GenericTrap::GetTwistAngle ( G4int index ) const

inline

Referenced by CreatePolyhedron().

GetVertex()

G4TwoVector G4GenericTrap::GetVertex ( G4int index ) const

inline

Referenced by CalculateExtent().

GetVertices()

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

inline

Referenced by G4GDMLWriteSolids::GenTrapWrite().

GetVisSubdivisions()

G4int G4GenericTrap::GetVisSubdivisions ( ) const

inline

Referenced by CreatePolyhedron().

GetZHalfLength()

G4double G4GenericTrap::GetZHalfLength ( ) const

inline

Accessors and modifiers.

Referenced by CalculateExtent(), and G4GDMLWriteSolids::GenTrapWrite().

Inside()

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

overridevirtual

Concrete implementations of the expected query interfaces for solids, as defined in the base class G4VSolid.

Implements G4VSolid.

Definition at line 144 of file G4GenericTrap.cc.

{
  G4double px = p.x(), py = p.y(), pz = p.z();
 
  // intersect edges by z = pz plane
  G4TwoVector pp[4];
  G4double z = (pz + fDz);
  for (auto i = 0; i < 4; ++i)
  {
    pp[i] = fVertices[i] + fDelta[i]*z;
  }
 
  // estimate distance to the solid
  G4double dist = std::abs(pz) - fDz;
  for (auto i = 0; i < 4; ++i)
  {
    if (fTwist[i] == 0.)
    {
      G4double dd = fSurf[i].D*px + fSurf[i].E*py + fSurf[i].F*pz + fSurf[i].G;
      dist = std::max(dd, dist);
    }
    else
    {
      auto j = (i + 1)%4;
      G4TwoVector a = pp[i];
      G4TwoVector b = pp[j];
      G4double dx = b.x() - a.x();
      G4double dy = b.y() - a.y();
      G4double leng = std::sqrt(dx*dx + dy*dy);
      G4double dd = (dx*(py - a.y()) - dy*(px - a.x()))/leng;
      dist = std::max(dd, dist);
    }
  }
  return (dist > halfTolerance) ? kOutside :
    ((dist > -halfTolerance) ? kSurface : kInside);
}

IsFaceted()

G4bool G4GenericTrap::IsFaceted ( ) const

overridevirtual

Returns true if the solid has only planar faces; false if twisted.

Reimplemented from G4VSolid.

Definition at line 952 of file G4GenericTrap.cc.

{
  return (!fIsTwisted);
}

IsTwisted()

G4bool G4GenericTrap::IsTwisted ( ) const

inline

operator=()

G4GenericTrap & G4GenericTrap::operator=

(

const G4GenericTrap &

rhs

)

Definition at line 111 of file G4GenericTrap.cc.

{
  // Check assignment to self
  if (this == &rhs)  { return *this; }
 
  // Copy base class data
  G4VSolid::operator=(rhs);
 
  // Copy data
  halfTolerance = rhs.halfTolerance; fScratch = rhs.fScratch;
  fDz = rhs.fDz; fVertices = rhs.fVertices; fIsTwisted = rhs.fIsTwisted;
  fMinBBox = rhs.fMinBBox; fMaxBBox = rhs.fMaxBBox;
  fVisSubdivisions = rhs.fVisSubdivisions;
  fSurfaceArea = rhs.fSurfaceArea; fCubicVolume = rhs.fCubicVolume;
 
  for (auto i = 0; i < 8; ++i) { fVertices[i] = rhs.fVertices[i]; }
  for (auto i = 0; i < 5; ++i) { fTwist[i] = rhs.fTwist[i]; }
  for (auto i = 0; i < 4; ++i) { fDelta[i] = rhs.fDelta[i]; }
  for (auto i = 0; i < 8; ++i) { fPlane[i] = rhs.fPlane[i]; }
  for (auto i = 0; i < 4; ++i) { fSurf[i] = rhs.fSurf[i]; }
  for (auto i = 0; i < 4; ++i) { f4k[i] = rhs.f4k[i]; }
  for (auto i = 0; i < 4; ++i) { fArea[i] = rhs.fArea[i]; }
 
  fRebuildPolyhedron = false;
  delete fpPolyhedron; fpPolyhedron = nullptr;
 
  return *this;
}

SetVisSubdivisions()

void G4GenericTrap::SetVisSubdivisions ( G4int subdiv )

inline

StreamInfo()

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

overridevirtual

Streams the object contents to an output stream.

Implements G4VSolid.

Definition at line 970 of file G4GenericTrap.cc.

{
  G4long oldprc = os.precision(16);
  os << "-----------------------------------------------------------\n"
     << "    *** Dump for solid - " << GetName() << " ***\n"
     << "    ===================================================\n"
     << "Solid geometry type: " << GetEntityType() << "\n"
     << "   half length Z: " << fDz/mm << "\n"
     << "   list of vertices:\n";
  for (G4int i = 0; i < 8; ++i)
  {
    os << "    #" << i << " " << fVertices[i] << "\n";
  }
  os << "-----------------------------------------------------------\n";
  os.precision(oldprc);
  return os;
}

SurfaceNormal()

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

overridevirtual

Returns the outwards pointing unit normal of the shape for the surface closest to the point at offset 'p'.

Parameters

[in]

p

The point at offset p.

Returns

The outwards pointing unit normal.

Implements G4VSolid.

Definition at line 185 of file G4GenericTrap.cc.

{
  G4double halfToleranceSquared = halfTolerance*halfTolerance;
  G4double px = p.x(), py = p.y(), pz = p.z();
  G4double nx = 0, ny = 0, nz = 0;
 
  // intersect edges by z = pz plane
  G4TwoVector pp[4];
  G4double tz = (pz + fDz);
  for (auto i = 0; i < 4; ++i)
  {
    pp[i] = fVertices[i] + fDelta[i]*tz;
  }
 
  // check bottom and top faces
  G4double dz = std::abs(pz) - fDz;
  nz = std::copysign(G4double(std::abs(dz) <= halfTolerance), pz);
 
  // check lateral faces
  for (auto i = 0; i < 4; ++i)
  {
    if (fTwist[i] == 0.)
    {
      G4double dd = fSurf[i].D*px + fSurf[i].E*py + fSurf[i].F*pz + fSurf[i].G;
      if (std::abs(dd) <= halfTolerance)
      {
        nx += fSurf[i].D;
        ny += fSurf[i].E;
        nz += fSurf[i].F;
      }
    }
    else
    {
      auto j = (i + 1)%4;
      G4TwoVector a = pp[i];
      G4TwoVector b = pp[j];
      G4double dx = b.x() - a.x();
      G4double dy = b.y() - a.y();
      G4double ll = dx*dx + dy*dy;
      G4double dd = dx*(py - a.y()) - dy*(px - a.x());
      if (dd*dd <= halfToleranceSquared*ll)
      {
        G4double x = fSurf[i].A*pz + fSurf[i].D;
        G4double y = fSurf[i].B*pz + fSurf[i].E;
        G4double z = fSurf[i].A*px + fSurf[i].B*py + 2.*fSurf[i].C*pz + fSurf[i].F;
        G4double inv = 1./std::sqrt(x*x + y*y + z*z);
        nx += x*inv;
        ny += y*inv;
        nz += z*inv;
      }
    }
  }
 
  // return normal
  G4double mag2 = nx*nx + ny*ny + nz*nz;
  if (mag2 == 0.)
  {
    return ApproxSurfaceNormal(p); // point is not on the surface
  }
  G4double mag = std::sqrt(mag2);
  G4double inv = (mag == 1.) ? 1. : 1./mag;
  return { nx*inv, ny*inv, nz*inv };
}

---

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

• G4GenericTrap.hh
• G4GenericTrap.cc