G4TwistTrapAlphaSide Class Reference

G4TwistTrapAlphaSide describes a twisted boundary surface for a trapezoid. More...

#include <G4TwistTrapAlphaSide.hh>

Inheritance diagram for G4TwistTrapAlphaSide:

Public Member Functions
	G4TwistTrapAlphaSide (const G4String &name, G4double PhiTwist, G4double pDz, G4double pTheta, G4double pPhi, G4double pDy1, G4double pDx1, G4double pDx2, G4double pDy2, G4double pDx3, G4double pDx4, G4double pAlph, G4double AngleSide)
	~G4TwistTrapAlphaSide () override=default
G4ThreeVector	GetNormal (const G4ThreeVector &p, G4bool isGlobal=false) override
G4int	DistanceToSurface (const G4ThreeVector &gp, const G4ThreeVector &gv, G4ThreeVector gxx[], G4double distance[], G4int areacode[], G4bool isvalid[], EValidate validate=kValidateWithTol) override
G4int	DistanceToSurface (const G4ThreeVector &gp, G4ThreeVector gxx[], G4double distance[], G4int areacode[]) override
	G4TwistTrapAlphaSide (__void__ &)
Public Member Functions inherited from G4VTwistSurface
	G4VTwistSurface (const G4String &name)
	G4VTwistSurface (const G4String &name, const G4RotationMatrix &rot, const G4ThreeVector &tlate, G4int handedness, const EAxis axis0, const EAxis axis1, G4double axis0min=-kInfinity, G4double axis1min=-kInfinity, G4double axis0max=kInfinity, G4double axis1max=kInfinity)
virtual	~G4VTwistSurface ()=default
virtual G4int	AmIOnLeftSide (const G4ThreeVector &me, const G4ThreeVector &vec, G4bool withTol=true)
virtual G4double	DistanceToBoundary (G4int areacode, G4ThreeVector &xx, const G4ThreeVector &p)
virtual G4double	DistanceToIn (const G4ThreeVector &gp, const G4ThreeVector &gv, G4ThreeVector &gxxbest)
virtual G4double	DistanceToOut (const G4ThreeVector &gp, const G4ThreeVector &gv, G4ThreeVector &gxxbest)
virtual G4double	DistanceTo (const G4ThreeVector &gp, G4ThreeVector &gxx)
virtual void	GetBoundaryParameters (const G4int &areacode, G4ThreeVector &d, G4ThreeVector &x0, G4int &boundarytype) const
virtual G4ThreeVector	GetBoundaryAtPZ (G4int areacode, const G4ThreeVector &p) const
G4double	DistanceToPlaneWithV (const G4ThreeVector &p, const G4ThreeVector &v, const G4ThreeVector &x0, const G4ThreeVector &n0, G4ThreeVector &xx)
G4double	DistanceToPlane (const G4ThreeVector &p, const G4ThreeVector &x0, const G4ThreeVector &n0, G4ThreeVector &xx)
G4double	DistanceToPlane (const G4ThreeVector &p, const G4ThreeVector &x0, const G4ThreeVector &t1, const G4ThreeVector &t2, G4ThreeVector &xx, G4ThreeVector &n)
G4double	DistanceToLine (const G4ThreeVector &p, const G4ThreeVector &x0, const G4ThreeVector &d, G4ThreeVector &xx)
G4bool	IsAxis0 (G4int areacode) const
G4bool	IsAxis1 (G4int areacode) const
G4bool	IsOutside (G4int areacode) const
G4bool	IsInside (G4int areacode, G4bool testbitmode=false) const
G4bool	IsBoundary (G4int areacode, G4bool testbitmode=false) const
G4bool	IsCorner (G4int areacode, G4bool testbitmode=false) const
G4bool	IsValidNorm () const
G4bool	IsSameBoundary (G4VTwistSurface *surface1, G4int areacode1, G4VTwistSurface *surface2, G4int areacode2) const
G4int	GetAxisType (G4int areacode, G4int whichaxis) const
G4ThreeVector	ComputeGlobalPoint (const G4ThreeVector &lp) const
G4ThreeVector	ComputeLocalPoint (const G4ThreeVector &gp) const
G4ThreeVector	ComputeGlobalDirection (const G4ThreeVector &lp) const
G4ThreeVector	ComputeLocalDirection (const G4ThreeVector &gp) const
void	SetAxis (G4int i, const EAxis axis)
void	SetNeighbours (G4VTwistSurface *ax0min, G4VTwistSurface *ax1min, G4VTwistSurface *ax0max, G4VTwistSurface *ax1max)
G4int	GetNode (G4int i, G4int j, G4int m, G4int n, G4int iside)
G4int	GetFace (G4int i, G4int j, G4int m, G4int n, G4int iside)
G4int	GetEdgeVisibility (G4int i, G4int j, G4int m, G4int n, G4int number, G4int orientation)
const G4String &	GetName () const
void	DebugPrint () const
	G4VTwistSurface (__void__ &)

Additional Inherited Members
Public Types inherited from G4VTwistSurface
enum	EValidate { kDontValidate = 0 , kValidateWithTol = 1 , kValidateWithoutTol = 2 , kUninitialized = 3 }
Static Public Attributes inherited from G4VTwistSurface
static const G4int	sOutside = 0x00000000
static const G4int	sInside = 0x10000000
static const G4int	sBoundary = 0x20000000
static const G4int	sCorner = 0x40000000
static const G4int	sC0Min1Min = 0x40000101
static const G4int	sC0Max1Min = 0x40000201
static const G4int	sC0Max1Max = 0x40000202
static const G4int	sC0Min1Max = 0x40000102
static const G4int	sAxisMin = 0x00000101
static const G4int	sAxisMax = 0x00000202
static const G4int	sAxisX = 0x00000404
static const G4int	sAxisY = 0x00000808
static const G4int	sAxisZ = 0x00000C0C
static const G4int	sAxisRho = 0x00001010
static const G4int	sAxisPhi = 0x00001414
static const G4int	sAxis0 = 0x0000FF00
static const G4int	sAxis1 = 0x000000FF
static const G4int	sSizeMask = 0x00000303
static const G4int	sAxisMask = 0x0000FCFC
static const G4int	sAreaMask = 0XF0000000
Protected Member Functions inherited from G4VTwistSurface
G4VTwistSurface **	GetNeighbours ()
G4int	GetNeighbours (G4int areacode, G4VTwistSurface *surfaces[])
G4ThreeVector	GetCorner (G4int areacode) const
void	GetBoundaryAxis (G4int areacode, EAxis axis[]) const
void	GetBoundaryLimit (G4int areacode, G4double limit[]) const
virtual void	SetBoundary (const G4int &axiscode, const G4ThreeVector &direction, const G4ThreeVector &x0, const G4int &boundarytype)
void	SetCorner (G4int areacode, G4double x, G4double y, G4double z)
Protected Attributes inherited from G4VTwistSurface
EAxis	fAxis [2]
G4double	fAxisMin [2]
G4double	fAxisMax [2]
CurrentStatus	fCurStatWithV
CurrentStatus	fCurStat
G4RotationMatrix	fRot
G4ThreeVector	fTrans
G4int	fHandedness
G4SurfCurNormal	fCurrentNormal
G4bool	fIsValidNorm
G4double	kCarTolerance

Detailed Description

G4TwistTrapAlphaSide describes a twisted boundary surface for a trapezoid.

Definition at line 46 of file G4TwistTrapAlphaSide.hh.

Constructor & Destructor Documentation

G4TwistTrapAlphaSide() [1/2]

G4TwistTrapAlphaSide::G4TwistTrapAlphaSide	(	const G4String &	name,
		G4double	PhiTwist,
		G4double	pDz,
		G4double	pTheta,
		G4double	pPhi,
		G4double	pDy1,
		G4double	pDx1,
		G4double	pDx2,
		G4double	pDy2,
		G4double	pDx3,
		G4double	pDx4,
		G4double	pAlph,
		G4double	AngleSide )

Constructs a trapezoid twisted boundary surface, given its parameters.

Parameters

[in]	name	The surface name.
[in]	PhiTwist	The twist angle.
[in]	pDz	Half z length.
[in]	pTheta	Direction between end planes - polar angle.
[in]	pPhi	Direction between end planes - azimuthal angle.
[in]	pDy1	Half y length at -pDz.
[in]	pDx1	Half x length at -pDz,-pDy.
[in]	pDx2	Half x length at -pDz,+pDy.
[in]	pDy2	Half y length at +pDz.
[in]	pDx3	Half x length at +pDz,-pDy.
[in]	pDx4	Half x length at +pDz,+pDy.
[in]	pAlph	Tilt angle at +pDz.
[in]	AngleSide	Parity.

Definition at line 40 of file G4TwistTrapAlphaSide.cc.

  : G4VTwistSurface(name)
{
  fAxis[0]    = kYAxis; // in local coordinate system
  fAxis[1]    = kZAxis;
  fAxisMin[0] = -kInfinity ;  // Y Axis boundary
  fAxisMax[0] = kInfinity ;   //   depends on z !!
  fAxisMin[1] = -pDz ;      // Z Axis boundary
  fAxisMax[1] = pDz ;
  
  fDx1  = pDx1 ;
  fDx2  = pDx2 ;
  fDx3  = pDx3 ;
  fDx4  = pDx4 ;
 
  fDy1   = pDy1 ;
  fDy2   = pDy2 ;
 
  fDz   = pDz ;
 
  fAlph = pAlph  ;
  fTAlph = std::tan(fAlph) ;
 
  fTheta = pTheta ;
  fPhi   = pPhi ;
 
  // precalculate frequently used parameters
  fDx4plus2  = fDx4 + fDx2 ;
  fDx4minus2 = fDx4 - fDx2 ;
  fDx3plus1  = fDx3 + fDx1 ; 
  fDx3minus1 = fDx3 - fDx1 ;
  fDy2plus1  = fDy2 + fDy1 ;
  fDy2minus1 = fDy2 - fDy1 ;
 
  fa1md1 = 2*fDx2 - 2*fDx1  ; 
  fa2md2 = 2*fDx4 - 2*fDx3 ;
 
  fPhiTwist = PhiTwist ;    // dphi
  fAngleSide = AngleSide ;  // 0,90,180,270 deg
 
  fdeltaX = 2 * fDz * std::tan(fTheta) * std::cos(fPhi);
    // dx in surface equation
  fdeltaY = 2 * fDz * std::tan(fTheta) * std::sin(fPhi);
    // dy in surface equation
  
  fRot.rotateZ( AngleSide ) ; 
  
  fTrans.set(0, 0, 0);  // No Translation
  fIsValidNorm = false;
  
  SetCorners() ;
  SetBoundaries() ;
}

~G4TwistTrapAlphaSide()

G4TwistTrapAlphaSide::~G4TwistTrapAlphaSide ( )

overridedefault

Default destructor.

G4TwistTrapAlphaSide() [2/2]

G4TwistTrapAlphaSide::G4TwistTrapAlphaSide

(

__void__ &

a

)

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

Definition at line 112 of file G4TwistTrapAlphaSide.cc.

  : G4VTwistSurface(a)
{
}

Member Function Documentation

DistanceToSurface() [1/2]

G4int G4TwistTrapAlphaSide::DistanceToSurface ( const G4ThreeVector & gp, const G4ThreeVector & gv, G4ThreeVector gxx[], G4double distance[], G4int areacode[], G4bool isvalid[], EValidate validate = kValidateWithTol )

overridevirtual

Returns the distance to surface, given point 'gp' and direction 'gv'.

Parameters

[in]	gp	The point from where computing the distance.
[in]	gv	The direction along which computing the distance.
[out]	gxx	Vector of global points based on number of solutions.
[out]	distance	The distance vector based on number of solutions.
[out]	areacode	The location vector based on number of solutions.
[out]	isvalid	Validity vector based on number of solutions.
[in]	validate	Adopted validation criteria.

Returns

The number of solutions.

Implements G4VTwistSurface.

Definition at line 176 of file G4TwistTrapAlphaSide.cc.

{
  static const G4double pihalf = pi/2 ;
  const G4double ctol = 0.5 * kCarTolerance;
 
  G4bool IsParallel = false ;
  G4bool IsConverged =  false ;
 
  G4int nxx = 0 ;  // number of physical solutions
 
  fCurStatWithV.ResetfDone(validate, &gp, &gv);
 
  if (fCurStatWithV.IsDone())
  {
    for (G4int i=0; i<fCurStatWithV.GetNXX(); ++i)
    {
      gxx[i] = fCurStatWithV.GetXX(i);
      distance[i] = fCurStatWithV.GetDistance(i);
      areacode[i] = fCurStatWithV.GetAreacode(i);
      isvalid[i]  = fCurStatWithV.IsValid(i);
    }
    return fCurStatWithV.GetNXX();
  }
 
  // initialise
  for (G4int j=0; j<G4VSURFACENXX ; ++j)
  {
    distance[j] = kInfinity;
    areacode[j] = sOutside;
    isvalid[j]  = false;
    gxx[j].set(kInfinity, kInfinity, kInfinity);
  }
 
  G4ThreeVector p = ComputeLocalPoint(gp);
  G4ThreeVector v = ComputeLocalDirection(gv);
  
#ifdef G4TWISTDEBUG
  G4cout << "Local point p = " << p << G4endl ;
  G4cout << "Local direction v = " << v << G4endl ; 
#endif
 
  G4double phi,u ;  // parameters
 
  // temporary variables
 
  G4double      tmpdist = kInfinity ;
  G4ThreeVector tmpxx;
  G4int         tmpareacode = sOutside ;
  G4bool        tmpisvalid  = false ;
 
  std::vector<Intersection> xbuf ;
  Intersection xbuftmp ;
  
  // prepare some variables for the intersection finder
 
  G4double L = 2*fDz ;
 
  G4double phixz = fPhiTwist * ( p.x() * v.z() - p.z() * v.x() ) ;
  G4double phiyz = fPhiTwist * ( p.y() * v.z() - p.z() * v.y() ) ;
 
 
  // special case vz = 0
 
  if ( v.z() == 0. )
  {
    if ( std::fabs(p.z()) <= L )   // intersection possible in z
    {
      phi = p.z() * fPhiTwist / L ;  // phi is determined by the z-position 
      u = (fDy1*(4*(-(fdeltaY*phi*v.x()) + fPhiTwist*p.y()*v.x()
                    + fdeltaX*phi*v.y() - fPhiTwist*p.x()*v.y())
                 + ((fDx3plus1 + fDx4plus2)*fPhiTwist
                 +  2*(fDx3minus1 + fDx4minus2)*phi)
                     *(v.y()*std::cos(phi) - v.x()*std::sin(phi))))
         /(fPhiTwist*(4*fDy1* v.x() - (fa1md1 + 4*fDy1*fTAlph)*v.y())
                    *std::cos(phi) + fPhiTwist*(fa1md1*v.x()
                       + 4*fDy1*(fTAlph*v.x() + v.y()))*std::sin(phi));
      xbuftmp.phi = phi ;
      xbuftmp.u = u ;
      xbuftmp.areacode = sOutside ;
      xbuftmp.distance = kInfinity ;
      xbuftmp.isvalid = false ;
 
      xbuf.push_back(xbuftmp) ;  // store it to xbuf
    }
    else   // no intersection possible
    {
      distance[0] = kInfinity;
      gxx[0].set(kInfinity,kInfinity,kInfinity);
      isvalid[0] = false ;
      areacode[0] = sOutside ;
      fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0],
                                     areacode[0], isvalid[0],
                                     0, validate, &gp, &gv);
      return 0;
    }  // end std::fabs(p.z() <= L 
  } // end v.z() == 0
  else  // general solution for non-zero vz
  {
 
    G4double c[8],srd[7],si[7] ;  
 
    c[7] = 57600*
      fDy1*(fa1md1*phiyz + 
            fDy1*(-4*phixz + 4*fTAlph*phiyz
                 + (fDx3plus1 + fDx4plus2)*fPhiTwist*v.z())) ;
    c[6] = -57600*
      fDy1*(4*fDy1*(phiyz + 2*fDz*v.x() + fTAlph*(phixz - 2*fDz*v.y()))
          - 2*fDy1*(2*fdeltaX + fDx3minus1 + fDx4minus2
                  - 2*fdeltaY*fTAlph)*v.z()
          + fa1md1*(phixz - 2*fDz*v.y() + fdeltaY*v.z()));
    c[5] = 4800*
      fDy1*(fa1md1*(-5*phiyz - 24*fDz*v.x() + 12*fdeltaX*v.z()) + 
            fDy1*(20*phixz - 4*(5*fTAlph*phiyz + 24*fDz*fTAlph*v.x()
                   + 24*fDz*v.y()) + (48*fdeltaY + (fDx3plus1 + fDx4plus2)
                                    *fPhiTwist + 48*fdeltaX*fTAlph)*v.z()));
    c[4] = 4800*
      fDy1*(fa1md1*(phixz - 10*fDz*v.y() + 5*fdeltaY*v.z())
          + 2*fDy1*(2*phiyz + 20*fDz*v.x()
                   + (-10*fdeltaX + fDx3minus1 + fDx4minus2)*v.z()
                   + 2*fTAlph*(phixz - 10*fDz*v.y() + 5*fdeltaY*v.z())));
    c[3] = -96*
      fDy1*(-(fa1md1*(phiyz + 100*fDz*v.x() - 50*fdeltaX*v.z()))
          + fDy1*(4*phixz - 400*fDz*v.y()
                  + (200*fdeltaY - (fDx3plus1 + fDx4plus2)*fPhiTwist)*v.z()
                  - 4*fTAlph*(phiyz + 100*fDz*v.x() - 50*fdeltaX*v.z())));
    c[2] = 32*
      fDy1*(4*fDy1*(7*fTAlph*phixz + 7*phiyz - 6*fDz*v.x() + 6*fDz*fTAlph*v.y())
          + 6*fDy1*(2*fdeltaX+fDx3minus1+fDx4minus2-2*fdeltaY*fTAlph)*v.z()
          + fa1md1*(7*phixz + 6*fDz*v.y() - 3*fdeltaY*v.z()));
    c[1] = -8*
      fDy1*(fa1md1*(-9*phiyz - 56*fDz*v.x() + 28*fdeltaX*v.z())
          + 4*fDy1*(9*phixz - 9*fTAlph*phiyz - 56*fDz*fTAlph*v.x()
                  - 56*fDz*v.y() + 28*(fdeltaY + fdeltaX*fTAlph)*v.z()));
    c[0] = 72*
      fDy1*(fa1md1*(2*fDz*v.y() - fdeltaY*v.z())
          + fDy1*(-8*fDz*v.x() + 8*fDz*fTAlph*v.y()
                + 4*fdeltaX*v.z() - 4*fdeltaY*fTAlph*v.z()));
 
#ifdef G4TWISTDEBUG
    G4cout << "coef = " << c[0] << " " 
           <<  c[1] << " "  
           <<  c[2] << " "  
           <<  c[3] << " "  
           <<  c[4] << " "  
           <<  c[5] << " "  
           <<  c[6] << " "  
           <<  c[7] << G4endl ; 
#endif    
 
    G4JTPolynomialSolver trapEq ;
    G4int num = trapEq.FindRoots(c,7,srd,si);
 
    for (G4int i = 0 ; i<num ; i++ )   // loop over all math solutions
    {  
      if ( si[i]==0.0 )  // only real solutions
      { 
#ifdef G4TWISTDEBUG
        G4cout << "Solution " << i << " : " << srd[i] << G4endl ;
#endif
        phi = std::fmod(srd[i] , pihalf)  ;
        u   = (fDy1*(4*(phiyz + 2*fDz*phi*v.y() - fdeltaY*phi*v.z())
                     - ((fDx3plus1 + fDx4plus2)*fPhiTwist
                       + 2*(fDx3minus1 + fDx4minus2)*phi)*v.z()*std::sin(phi)))
             /(fPhiTwist*v.z()*(4*fDy1*std::cos(phi)
                             + (fa1md1 + 4*fDy1*fTAlph)*std::sin(phi)));        
        xbuftmp.phi = phi ;
        xbuftmp.u = u ;
        xbuftmp.areacode = sOutside ;
        xbuftmp.distance = kInfinity ;
        xbuftmp.isvalid = false ;
        
        xbuf.push_back(xbuftmp) ;  // store it to xbuf
      
#ifdef G4TWISTDEBUG
        G4cout << "solution " << i << " = " << phi << " , " << u  << G4endl ;
#endif
      }  // end if real solution
    }  // end loop i    
  }    // end general case
 
  nxx = (G4int)xbuf.size() ;  // save the number of  solutions
 
  G4ThreeVector xxonsurface  ;    // point on surface
  G4ThreeVector surfacenormal  ;  // normal vector  
  G4double deltaX;  // distance between intersection point and point on surface
  G4double theta;   // angle between track and surfacenormal
  G4double factor;  // a scaling factor
  G4int maxint=30;  // number of iterations
 
  for (auto & k : xbuf)
  {
#ifdef G4TWISTDEBUG
    G4cout << "Solution " << k << " : " 
           << "reconstructed phiR = " << xbuf[k].phi
           << ", uR = " << xbuf[k].u << G4endl ; 
#endif
    
    phi = k.phi ;  // get the stored values for phi and u
    u = k.u ;
 
    IsConverged = false ;   // no convergence at the beginning
    
    for ( G4int i = 1 ; i<maxint ; ++i )
    {
      xxonsurface = SurfacePoint(phi,u) ;
      surfacenormal = NormAng(phi,u) ;
 
      tmpdist = DistanceToPlaneWithV(p, v, xxonsurface, surfacenormal, tmpxx); 
      deltaX = ( tmpxx - xxonsurface ).mag() ; 
      theta = std::fabs(std::acos(v*surfacenormal) - pihalf) ;
      if ( theta < 0.001 )
      { 
        factor = 50 ;
        IsParallel = true ;
      }
      else
      {
        factor = 1 ;
      }
 
#ifdef G4TWISTDEBUG
      G4cout << "Step i = " << i << ", distance = " << tmpdist
             << ", " << deltaX << G4endl ;
      G4cout << "X = " << tmpxx << G4endl ;
#endif
      
      GetPhiUAtX(tmpxx, phi, u) ;
        // the new point xx is accepted and phi/u replaced
      
#ifdef G4TWISTDEBUG
      G4cout << "approximated phi = " << phi << ", u = " << u << G4endl ; 
#endif
      
      if ( deltaX <= factor*ctol ) { IsConverged = true ; break ; }
      
    }  // end iterative loop (i)
    
    if ( std::fabs(tmpdist)<ctol ) { tmpdist = 0 ; }
 
#ifdef G4TWISTDEBUG
    G4cout << "refined solution "  << phi << " , " << u  <<  G4endl ;
    G4cout << "distance = " << tmpdist << G4endl ;
    G4cout << "local X = " << tmpxx << G4endl ;
#endif
    
    tmpisvalid = false ;  // init 
 
    if ( IsConverged )
    {
      if (validate == kValidateWithTol)
      {
        tmpareacode = GetAreaCode(tmpxx);
        if (!IsOutside(tmpareacode))
        {
          if (tmpdist >= 0) { tmpisvalid = true; }
        }
      }
      else if (validate == kValidateWithoutTol)
      {
        tmpareacode = GetAreaCode(tmpxx, false);
        if (IsInside(tmpareacode))
        {
          if (tmpdist >= 0) { tmpisvalid = true; }
        }
      }
      else  // kDontValidate
      {
        G4Exception("G4TwistTrapAlphaSide::DistanceToSurface()",
                    "GeomSolids0001", FatalException,
                    "Feature NOT implemented !");
      }
    } 
    else
    {
      tmpdist = kInfinity;     // no convergence after 10 steps 
      tmpisvalid = false ;     // solution is not vaild
    }  
 
    // store the found values
    // 
    k.xx = tmpxx ;
    k.distance = tmpdist ;
    k.areacode = tmpareacode ;
    k.isvalid = tmpisvalid ;
 
  }  // end loop over physical solutions (variable k)
 
  std::sort(xbuf.begin() , xbuf.end(), DistanceSort ) ;  // sorting
 
#ifdef G4TWISTDEBUG
  G4cout << G4endl << "list xbuf after sorting : " << G4endl ;
  G4cout << G4endl << G4endl ;
#endif
 
  // erase identical intersection (within kCarTolerance) 
  //
  xbuf.erase( std::unique(xbuf.begin(), xbuf.end() , EqualIntersection ),
              xbuf.end() );
 
 
  // add guesses
  //
  auto nxxtmp = (G4int)xbuf.size() ;
 
  if ( nxxtmp<2 || IsParallel  )  // positive end
  {
 
#ifdef G4TWISTDEBUG
    G4cout << "add guess at +z/2 .. " << G4endl ;
#endif
 
    phi = fPhiTwist/2 ;
    u   = 0 ;
    
    xbuftmp.phi = phi ;
    xbuftmp.u = u ;
    xbuftmp.areacode = sOutside ;
    xbuftmp.distance = kInfinity ;
    xbuftmp.isvalid = false ;
    
    xbuf.push_back(xbuftmp) ;  // store it to xbuf
 
#ifdef G4TWISTDEBUG
    G4cout << "add guess at -z/2 .. " << G4endl ;
#endif
 
    phi = -fPhiTwist/2 ;
    u   = 0 ;
    
    xbuftmp.phi = phi ;
    xbuftmp.u = u ;
    xbuftmp.areacode = sOutside ;
    xbuftmp.distance = kInfinity ;
    xbuftmp.isvalid = false ;
    
    xbuf.push_back(xbuftmp) ;  // store it to xbuf
 
    for ( std::size_t k = nxxtmp ; k<xbuf.size() ; ++k )
    {
 
#ifdef G4TWISTDEBUG
      G4cout << "Solution " << k << " : " 
             << "reconstructed phiR = " << xbuf[k].phi
             << ", uR = " << xbuf[k].u << G4endl ; 
#endif
      
      phi = xbuf[k].phi ;  // get the stored values for phi and u
      u   = xbuf[k].u ;
 
      IsConverged = false ;   // no convergence at the beginning
      
      for ( G4int i = 1 ; i<maxint ; ++i )
      {
        xxonsurface = SurfacePoint(phi,u) ;
        surfacenormal = NormAng(phi,u) ;
        tmpdist = DistanceToPlaneWithV(p, v, xxonsurface, surfacenormal, tmpxx);
        deltaX = ( tmpxx - xxonsurface ).mag();
        theta = std::fabs(std::acos(v*surfacenormal) - pihalf);
        if ( theta < 0.001 )
        { 
          factor = 50 ;    
        }
        else
        {
          factor = 1 ;
        }
        
#ifdef G4TWISTDEBUG
        G4cout << "Step i = " << i << ", distance = " << tmpdist
               << ", " << deltaX << G4endl
               << "X = " << tmpxx << G4endl ;
#endif
 
        GetPhiUAtX(tmpxx, phi, u) ;
          // the new point xx is accepted and phi/u replaced
      
#ifdef G4TWISTDEBUG
        G4cout << "approximated phi = " << phi << ", u = " << u << G4endl ; 
#endif
      
        if ( deltaX <= factor*ctol ) { IsConverged = true ; break ; }
      
      }  // end iterative loop (i)
    
      if ( std::fabs(tmpdist)<ctol )  { tmpdist = 0; }
 
#ifdef G4TWISTDEBUG
      G4cout << "refined solution "  << phi << " , " << u  <<  G4endl ;
      G4cout << "distance = " << tmpdist << G4endl ;
      G4cout << "local X = " << tmpxx << G4endl ;
#endif
 
      tmpisvalid = false ;  // init 
 
      if ( IsConverged )
      {
        if (validate == kValidateWithTol)
        {
          tmpareacode = GetAreaCode(tmpxx);
          if (!IsOutside(tmpareacode))
          {
            if (tmpdist >= 0)  { tmpisvalid = true; }
          }
        }
        else if (validate == kValidateWithoutTol)
        {
          tmpareacode = GetAreaCode(tmpxx, false);
          if (IsInside(tmpareacode))
          {
            if (tmpdist >= 0)  { tmpisvalid = true; }
          }
        }
        else  // kDontValidate
        {
          G4Exception("G4TwistedBoxSide::DistanceToSurface()",
                      "GeomSolids0001", FatalException,
                      "Feature NOT implemented !");
        }
      } 
      else
      {
        tmpdist = kInfinity;     // no convergence after 10 steps 
        tmpisvalid = false ;     // solution is not vaild
      }  
 
      // store the found values
      //
      xbuf[k].xx = tmpxx ;
      xbuf[k].distance = tmpdist ;
      xbuf[k].areacode = tmpareacode ;
      xbuf[k].isvalid = tmpisvalid ;
 
    }  // end loop over physical solutions 
  }  // end less than 2 solutions
 
  // sort again
  std::sort(xbuf.begin() , xbuf.end(), DistanceSort ) ;  // sorting
 
  // erase identical intersection (within kCarTolerance) 
  xbuf.erase( std::unique(xbuf.begin(), xbuf.end() , EqualIntersection ) ,
              xbuf.end() );
 
#ifdef G4TWISTDEBUG
  G4cout << G4endl << "list xbuf after sorting : " << G4endl ;
  G4cout << G4endl << G4endl ;
#endif
 
  nxx = (G4int)xbuf.size() ;   // determine number of solutions again.
 
  for ( G4int i = 0 ; i<(G4int)xbuf.size() ; ++i )
  {
    distance[i] = xbuf[i].distance;
    gxx[i]      = ComputeGlobalPoint(xbuf[i].xx);
    areacode[i] = xbuf[i].areacode ;
    isvalid[i]  = xbuf[i].isvalid ;
    
    fCurStatWithV.SetCurrentStatus(i, gxx[i], distance[i], areacode[i],
                                     isvalid[i], nxx, validate, &gp, &gv);
#ifdef G4TWISTDEBUG
    G4cout << "element Nr. " << i 
           << ", local Intersection = " << xbuf[i].xx 
           << ", distance = " << xbuf[i].distance 
           << ", u = " << xbuf[i].u 
           << ", phi = " << xbuf[i].phi 
           << ", isvalid = " << xbuf[i].isvalid 
           << G4endl ;
#endif
 
  }  // end for( i ) loop
 
#ifdef G4TWISTDEBUG
  G4cout << "G4TwistTrapAlphaSide finished " << G4endl ;
  G4cout << nxx << " possible physical solutions found" << G4endl ;
  for ( G4int k= 0 ; k< nxx ; k++ )
  {
    G4cout << "global intersection Point found: " << gxx[k] << G4endl ;
    G4cout << "distance = " << distance[k] << G4endl ;
    G4cout << "isvalid = " << isvalid[k] << G4endl ;
  }
#endif
 
  return nxx ;
}

DistanceToSurface() [2/2]

G4int G4TwistTrapAlphaSide::DistanceToSurface ( const G4ThreeVector & gp, G4ThreeVector gxx[], G4double distance[], G4int areacode[] )

overridevirtual

Returns the safety distance to surface, given point 'gp'.

Parameters

[in]	gp	The point from where computing the safety distance.
[out]	gxx	Vector of global points based on number of solutions.
[out]	distance	The distance vector based on number of solutions.
[out]	areacode	The location vector based on number of solutions.

Returns

The number of solutions.

Implements G4VTwistSurface.

Definition at line 671 of file G4TwistTrapAlphaSide.cc.

{  
   const G4double ctol = 0.5 * kCarTolerance;
 
   fCurStat.ResetfDone(kDontValidate, &gp);
 
   if (fCurStat.IsDone())
   {
      for (G4int i=0; i<fCurStat.GetNXX(); ++i)
      {
         gxx[i] = fCurStat.GetXX(i);
         distance[i] = fCurStat.GetDistance(i);
         areacode[i] = fCurStat.GetAreacode(i);
      }
      return fCurStat.GetNXX();
   }
 
   // initialize
   for (G4int i=0; i<G4VSURFACENXX; ++i)
   {
      distance[i] = kInfinity;
      areacode[i] = sOutside;
      gxx[i].set(kInfinity, kInfinity, kInfinity);
   }
 
   G4ThreeVector p = ComputeLocalPoint(gp);
   G4ThreeVector xx;  // intersection point
   G4ThreeVector xxonsurface ; // interpolated intersection point 
 
   // the surfacenormal at that surface point
   //
   G4double phiR = 0  ;
   G4double uR = 0 ;
 
   G4ThreeVector surfacenormal ; 
   G4double deltaX, uMax ;
   G4double halfphi = 0.5*fPhiTwist ;
   
   for ( G4int i = 1 ; i<20 ; ++i )
   {
     xxonsurface = SurfacePoint(phiR,uR) ;
     surfacenormal = NormAng(phiR,uR) ;
     distance[0] = DistanceToPlane(p,xxonsurface,surfacenormal,xx); // new XX
     deltaX = ( xx - xxonsurface ).mag() ; 
 
#ifdef G4TWISTDEBUG
     G4cout << "i = " << i << ", distance = " << distance[0]
            << ", " << deltaX << G4endl
            << "X = " << xx << G4endl ;
#endif
 
     // the new point xx is accepted and phi/psi replaced
     //
     GetPhiUAtX(xx, phiR, uR) ;
     
     if ( deltaX <= ctol ) { break ; }
   }
 
   // check validity of solution ( valid phi,psi ) 
 
   uMax = GetBoundaryMax(phiR) ;
 
   if (  phiR > halfphi ) { phiR =  halfphi ; }
   if ( phiR < -halfphi ) { phiR = -halfphi ; }
   if ( uR > uMax )  { uR = uMax ;  }
   if ( uR < -uMax ) { uR = -uMax ; }
 
   xxonsurface = SurfacePoint(phiR,uR) ;
   distance[0] = (  p - xx ).mag() ;
   if ( distance[0] <= ctol ) { distance[0] = 0 ; } 
 
   // end of validity 
 
#ifdef G4TWISTDEBUG
   G4cout << "refined solution "  << phiR << " , " << uR << " , " << G4endl ;
   G4cout << "distance = " << distance[0] << G4endl ;
   G4cout << "X = " << xx << G4endl ;
#endif
 
   G4bool isvalid = true;
   gxx[0]      = ComputeGlobalPoint(xx);
   
#ifdef G4TWISTDEBUG
   G4cout << "intersection Point found: " << gxx[0] << G4endl ;
   G4cout << "distance = " << distance[0] << G4endl ;
#endif
 
   fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
                            isvalid, 1, kDontValidate, &gp);
   return 1;
}

GetNormal()

G4ThreeVector G4TwistTrapAlphaSide::GetNormal ( const G4ThreeVector & p, G4bool isGlobal = false )

overridevirtual

Returns a normal vector at a surface (or very close to the surface) point at 'p'.

Parameters

[in]	p	The point where computing the normal.
[in]	isGlobal	If true, it returns the normal in global coordinates.

Returns

The normal vector.

Implements G4VTwistSurface.

Definition at line 122 of file G4TwistTrapAlphaSide.cc.

{
   // GetNormal returns a normal vector at a surface (or very close
   // to surface) point at tmpxx.
   // If isGlobal=true, it returns the normal in global coordinate.
   //
 
   G4ThreeVector xx;
   if (isGlobal)
   {
      xx = ComputeLocalPoint(tmpxx);
      if ((xx - fCurrentNormal.p).mag() < 0.5 * kCarTolerance)
      {
         return ComputeGlobalDirection(fCurrentNormal.normal);
      }
   }
   else
   {
      xx = tmpxx;
      if (xx == fCurrentNormal.p)
      {
         return fCurrentNormal.normal;
      }
   }
 
   G4double phi ;
   G4double u ;
 
   GetPhiUAtX(xx,phi,u) ;   // phi,u for point xx close to surface 
 
   G4ThreeVector normal =  NormAng(phi,u) ;  // the normal vector at phi,u
 
#ifdef G4TWISTDEBUG
   G4cout  << "normal vector = " << normal << G4endl ;
   G4cout << "phi = " << phi << " , u = " << u << G4endl ;
#endif
 
   if (isGlobal)
   {
      fCurrentNormal.normal = ComputeGlobalDirection(normal.unit());
   }
   else
   {
      fCurrentNormal.normal = normal.unit();
   }
 
   return fCurrentNormal.normal;
}

---

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

• G4TwistTrapAlphaSide.hh
• G4TwistTrapAlphaSide.cc