G4TwistTrapAlphaSide Class Reference

#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)
virtual	~G4TwistTrapAlphaSide ()
virtual G4ThreeVector	GetNormal (const G4ThreeVector &xx, G4bool isGlobal=false)
virtual G4int	DistanceToSurface (const G4ThreeVector &gp, const G4ThreeVector &gv, G4ThreeVector gxx[], G4double distance[], G4int areacode[], G4bool isvalid[], EValidate validate=kValidateWithTol)
virtual G4int	DistanceToSurface (const G4ThreeVector &gp, G4ThreeVector gxx[], G4double distance[], G4int areacode[])
	G4TwistTrapAlphaSide (__void__ &)

---

Detailed Description

Definition at line 51 of file G4TwistTrapAlphaSide.hh.

---

Constructor & Destructor Documentation

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
	)

Definition at line 52 of file G4TwistTrapAlphaSide.cc.

References G4VTwistSurface::fAxis, G4VTwistSurface::fAxisMax, G4VTwistSurface::fAxisMin, G4VTwistSurface::fIsValidNorm, G4VTwistSurface::fRot, G4VTwistSurface::fTrans, kYAxis, and kZAxis.

  : 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

(

)

[virtual]

Definition at line 137 of file G4TwistTrapAlphaSide.cc.

{
}

G4TwistTrapAlphaSide::G4TwistTrapAlphaSide

(

__void__ &

)

Definition at line 124 of file G4TwistTrapAlphaSide.cc.

  : G4VTwistSurface(a), fTheta(0.), fPhi(0.), fDy1(0.), fDx1(0.), fDx2(0.), 
    fDy2(0.), fDx3(0.), fDx4(0.), fDz(0.), fAlph(0.), fTAlph(0.), fPhiTwist(0.), 
    fAngleSide(0.), fDx4plus2(0.), fDx4minus2(0.), fDx3plus1(0.), fDx3minus1(0.), 
    fDy2plus1(0.), fDy2minus1(0.), fa1md1(0.), fa2md2(0.), fdeltaX(0.),
    fdeltaY(0.)
{
}

---

Member Function Documentation

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

Implements G4VTwistSurface.

Definition at line 696 of file G4TwistTrapAlphaSide.cc.

References G4VTwistSurface::ComputeGlobalPoint(), G4VTwistSurface::ComputeLocalPoint(), G4VTwistSurface::DistanceToPlane(), G4VTwistSurface::fCurStat, G4cout, G4endl, G4VSURFACENXX, G4VTwistSurface::kCarTolerance, G4VTwistSurface::kDontValidate, and G4VTwistSurface::sOutside.

{  
  static const G4double ctol = 0.5 * kCarTolerance;

  fCurStat.ResetfDone(kDontValidate, &gp);

   if (fCurStat.IsDone())
   {
      for (register int i=0; i<fCurStat.GetNXX(); i++)
      {
         gxx[i] = fCurStat.GetXX(i);
         distance[i] = fCurStat.GetDistance(i);
         areacode[i] = fCurStat.GetAreacode(i);
      }
      return fCurStat.GetNXX();
   }
   else  // initialize
   {
      for (register int 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 ( register int 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;
}

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

Definition at line 200 of file G4TwistTrapAlphaSide.cc.

References Intersection::areacode, G4VTwistSurface::ComputeGlobalPoint(), G4VTwistSurface::ComputeLocalDirection(), G4VTwistSurface::ComputeLocalPoint(), Intersection::distance, DistanceSort(), G4VTwistSurface::DistanceToPlaneWithV(), EqualIntersection(), FatalException, G4VTwistSurface::fCurStatWithV, G4JTPolynomialSolver::FindRoots(), G4cout, G4endl, G4Exception(), G4VSURFACENXX, G4VTwistSurface::IsInside(), G4VTwistSurface::IsOutside(), Intersection::isvalid, G4VTwistSurface::kCarTolerance, G4VTwistSurface::kValidateWithoutTol, G4VTwistSurface::kValidateWithTol, Intersection::phi, G4INCL::Math::pi, G4VTwistSurface::sOutside, and Intersection::u.

{
  static const G4double ctol = 0.5 * kCarTolerance;
  static const G4double pihalf = pi/2 ;

  G4bool IsParallel = false ;
  G4bool IsConverged =  false ;

  G4int nxx = 0 ;  // number of physical solutions

  fCurStatWithV.ResetfDone(validate, &gp, &gv);

  if (fCurStatWithV.IsDone())
  {
    for (register int 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();
  }
  else  // initialise
  {
    for (register int 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 (register int 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 = 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 ( register size_t k = 0 ; 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 ( register int 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
    // 
    xbuf[k].xx = tmpxx ;
    xbuf[k].distance = tmpdist ;
    xbuf[k].areacode = tmpareacode ;
    xbuf[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
  //
  G4int nxxtmp = 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 ( register 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 ( register int 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 = xbuf.size() ;   // determine number of solutions again.

  for ( register size_t i = 0 ; i<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 ;
}

G4ThreeVector G4TwistTrapAlphaSide::GetNormal	(	const G4ThreeVector &	xx,
		G4bool	isGlobal = false
	)			[virtual]

Implements G4VTwistSurface.

Definition at line 146 of file G4TwistTrapAlphaSide.cc.

References G4VTwistSurface::ComputeGlobalDirection(), G4VTwistSurface::ComputeLocalPoint(), G4VTwistSurface::fCurrentNormal, G4cout, G4endl, and G4VTwistSurface::kCarTolerance.

{
   // 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

---