G4TwistTubsSide Class Reference

#include <G4TwistTubsSide.hh>

Inheritance diagram for G4TwistTubsSide:

Public Member Functions
	G4TwistTubsSide (const G4String &name, const G4RotationMatrix &rot, const G4ThreeVector &tlate, G4int handedness, const G4double kappa, const EAxis axis0=kXAxis, const EAxis axis1=kZAxis, G4double axis0min=-kInfinity, G4double axis1min=-kInfinity, G4double axis0max=kInfinity, G4double axis1max=kInfinity)
	G4TwistTubsSide (const G4String &name, G4double EndInnerRadius[2], G4double EndOuterRadius[2], G4double DPhi, G4double EndPhi[2], G4double EndZ[2], G4double InnerRadius, G4double OuterRadius, G4double Kappa, G4int handedness)
virtual	~G4TwistTubsSide ()
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[])
G4ThreeVector	ProjectAtPXPZ (const G4ThreeVector &p, G4bool isglobal=false) const
virtual G4ThreeVector	SurfacePoint (G4double, G4double, G4bool isGlobal=false)
virtual G4double	GetBoundaryMin (G4double phi)
virtual G4double	GetBoundaryMax (G4double phi)
virtual G4double	GetSurfaceArea ()
virtual void	GetFacets (G4int m, G4int n, G4double xyz[][3], G4int faces[][4], G4int iside)
	G4TwistTubsSide (__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 axis1, const EAxis axis2, G4double axis0min=-kInfinity, G4double axis1min=-kInfinity, G4double axis0max=kInfinity, G4double axis1max=kInfinity)
virtual	~G4VTwistSurface ()
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)
void	DebugPrint () const
virtual G4String	GetName () const
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 *axis0min, G4VTwistSurface *axis1min, G4VTwistSurface *axis0max, G4VTwistSurface *axis1max)
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)
	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

Definition at line 52 of file G4TwistTubsSide.hh.

Constructor & Destructor Documentation

G4TwistTubsSide::G4TwistTubsSide	(	const G4String &	name,
		const G4RotationMatrix &	rot,
		const G4ThreeVector &	tlate,
		G4int	handedness,
		const G4double	kappa,
		const EAxis	axis0 = kXAxis,
		const EAxis	axis1 = kZAxis,
		G4double	axis0min = -kInfinity,
		G4double	axis1min = -kInfinity,
		G4double	axis0max = kInfinity,
		G4double	axis1max = kInfinity
	)

Definition at line 50 of file G4TwistTubsSide.cc.

References FatalErrorInArgument, G4VTwistSurface::fIsValidNorm, G4Exception(), kXAxis, and kZAxis.

   : G4VTwistSurface(name, rot, tlate, handedness, axis0, axis1,
                axis0min, axis1min, axis0max, axis1max),
     fKappa(kappa)
{
   if (axis0 == kZAxis && axis1 == kXAxis) {
      G4Exception("G4TwistTubsSide::G4TwistTubsSide()", "GeomSolids0002",
                  FatalErrorInArgument, "Should swap axis0 and axis1!");
   }
   fIsValidNorm = false;
   SetCorners();
   SetBoundaries();
}

G4TwistTubsSide::G4TwistTubsSide	(	const G4String &	name,
		G4double	EndInnerRadius[2],
		G4double	EndOuterRadius[2],
		G4double	DPhi,
		G4double	EndPhi[2],
		G4double	EndZ[2],
		G4double	InnerRadius,
		G4double	OuterRadius,
		G4double	Kappa,
		G4int	handedness
	)

Definition at line 74 of file G4TwistTubsSide.cc.

References G4VTwistSurface::fAxis, G4VTwistSurface::fAxisMax, G4VTwistSurface::fAxisMin, G4VTwistSurface::fHandedness, G4VTwistSurface::fIsValidNorm, G4VTwistSurface::fRot, G4VTwistSurface::fTrans, kXAxis, kZAxis, CLHEP::HepRotation::rotateZ(), and CLHEP::Hep3Vector::set().

  : G4VTwistSurface(name)
{  
   fHandedness = handedness;   // +z = +ve, -z = -ve
   fAxis[0]    = kXAxis; // in local coordinate system
   fAxis[1]    = kZAxis;
   fAxisMin[0] = InnerRadius;  // Inner-hype radius at z=0
   fAxisMax[0] = OuterRadius;  // Outer-hype radius at z=0
   fAxisMin[1] = EndZ[0];
   fAxisMax[1] = EndZ[1];

   fKappa = Kappa;
   fRot.rotateZ( fHandedness > 0
                 ? -0.5*DPhi
                 :  0.5*DPhi );
   fTrans.set(0, 0, 0);
   fIsValidNorm = false;
   
   SetCorners( EndInnerRadius, EndOuterRadius, EndPhi, EndZ) ;
   SetBoundaries();
}

G4TwistTubsSide::~G4TwistTubsSide ( )

virtual

Definition at line 118 of file G4TwistTubsSide.cc.

{
}

G4TwistTubsSide::G4TwistTubsSide

(

__void__ &

a

)

Definition at line 109 of file G4TwistTubsSide.cc.

  : G4VTwistSurface(a), fKappa(0.)
{
}

Member Function Documentation

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

virtual

Implements G4VTwistSurface.

Definition at line 160 of file G4TwistTubsSide.cc.

References test::a, test::b, test::c, G4VTwistSurface::ComputeGlobalPoint(), G4VTwistSurface::ComputeLocalDirection(), G4VTwistSurface::ComputeLocalPoint(), DBL_MIN, G4VTwistSurface::DistanceToPlaneWithV(), FatalException, G4VTwistSurface::fCurStat, G4VTwistSurface::fCurStatWithV, G4endl, G4Exception(), G4VTwistSurface::CurrentStatus::GetAreacode(), G4VTwistSurface::CurrentStatus::GetDistance(), GetNormal(), G4VTwistSurface::CurrentStatus::GetNXX(), G4VTwistSurface::CurrentStatus::GetXX(), G4VTwistSurface::CurrentStatus::IsDone(), G4VTwistSurface::IsInside(), G4VTwistSurface::IsOutside(), G4VTwistSurface::CurrentStatus::IsValid(), G4VTwistSurface::kCarTolerance, G4VTwistSurface::kValidateWithoutTol, G4VTwistSurface::kValidateWithTol, G4VTwistSurface::CurrentStatus::ResetfDone(), CLHEP::Hep3Vector::set(), G4VTwistSurface::CurrentStatus::SetCurrentStatus(), G4VTwistSurface::sInside, G4VTwistSurface::sOutside, test::v, test::x, CLHEP::Hep3Vector::x(), CLHEP::Hep3Vector::y(), z, and CLHEP::Hep3Vector::z().

{
   // Coordinate system:
   //
   //    The coordinate system is so chosen that the intersection of
   //    the twisted surface with the z=0 plane coincides with the
   //    x-axis. 
   //    Rotation matrix from this coordinate system (local system)
   //    to global system is saved in fRot field.
   //    So the (global) particle position and (global) velocity vectors, 
   //    p and v, should be rotated fRot.inverse() in order to convert
   //    to local vectors.
   //
   // Equation of a twisted surface:
   //
   //    x(rho(z=0), z) = rho(z=0)
   //    y(rho(z=0), z) = rho(z=0)*K*z
   //    z(rho(z=0), z) = z
   //    with
   //       K = std::tan(fPhiTwist/2)/fZHalfLen
   //
   // Equation of a line:
   //
   //    gxx = p + t*v
   //    with
   //       p = fRot.inverse()*gp  
   //       v = fRot.inverse()*gv
   //
   // Solution for intersection:
   //
   //    Required time for crossing is given by solving the
   //    following quadratic equation:
   //
   //       a*t^2 + b*t + c = 0
   //
   //    where
   //
   //       a = K*v_x*v_z
   //       b = K*(v_x*p_z + v_z*p_x) - v_y
   //       c = K*p_x*p_z - p_y
   //
   //    Out of the possible two solutions you must choose
   //    the one that gives a positive rho(z=0).
   //
   //
      
   fCurStatWithV.ResetfDone(validate, &gp, &gv);

   if (fCurStatWithV.IsDone()) {
      G4int i;
      for (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 {
      // initialize
      G4int i;
      for (i=0; i<2; i++) {
         distance[i] = kInfinity;
         areacode[i] = sOutside;
         isvalid[i]  = false;
         gxx[i].set(kInfinity, kInfinity, kInfinity);
      }
   }

   G4ThreeVector p = ComputeLocalPoint(gp);
   G4ThreeVector v = ComputeLocalDirection(gv);
   G4ThreeVector xx[2]; 


   // 
   // special case! 
   // p is origin or
   //

   G4double absvz = std::fabs(v.z());

   if ((absvz < DBL_MIN) && (std::fabs(p.x() * v.y() - p.y() * v.x()) < DBL_MIN)) {
      // no intersection

      isvalid[0] = false;
      fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
                                isvalid[0], 0, validate, &gp, &gv);
      return 0;
   } 
   
   // 
   // special case end
   //

   
   G4double a = fKappa * v.x() * v.z();
   G4double b = fKappa * (v.x() * p.z() + v.z() * p.x()) - v.y();
   G4double c = fKappa * p.x() * p.z() - p.y();
   G4double D = b * b - 4 * a * c;             // discriminant
   G4int vout = 0;

   if (std::fabs(a) < DBL_MIN) {
      if (std::fabs(b) > DBL_MIN) { 

         // single solution

         distance[0] = - c / b;
         xx[0]       = p + distance[0]*v;
         gxx[0]      = ComputeGlobalPoint(xx[0]);

         if (validate == kValidateWithTol) {
            areacode[0] = GetAreaCode(xx[0]);
            if (!IsOutside(areacode[0])) {
               if (distance[0] >= 0) isvalid[0] = true;
            }
         } else if (validate == kValidateWithoutTol) {
            areacode[0] = GetAreaCode(xx[0], false);
            if (IsInside(areacode[0])) {
               if (distance[0] >= 0) isvalid[0] = true;
            }
         } else { // kDontValidate                       
            // we must omit x(rho,z) = rho(z=0) < 0
            if (xx[0].x() > 0) {
               areacode[0] = sInside;
               if (distance[0] >= 0) isvalid[0] = true;
            } else {
               distance[0] = kInfinity;
               fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0],
                                              areacode[0], isvalid[0],
                                              0, validate, &gp, &gv);
               return vout;
            } 
         }
                 
         fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
                                        isvalid[0], 1, validate, &gp, &gv);
         vout = 1;

      } else {
         // if a=b=0 , v.y=0 and (v.x=0 && p.x=0) or (v.z=0 && p.z=0) .
         //    if v.x=0 && p.x=0, no intersection unless p is on z-axis
         //    (in that case, v is paralell to surface). 
         //    if v.z=0 && p.z=0, no intersection unless p is on x-axis
         //    (in that case, v is paralell to surface). 
         // return distance = infinity.

         fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
                                        isvalid[0], 0, validate, &gp, &gv);
      }
      
   } else if (D > DBL_MIN) {   

      // double solutions

      D = std::sqrt(D);
      G4double      factor = 0.5/a;
      G4double      tmpdist[2] = {kInfinity, kInfinity};
      G4ThreeVector tmpxx[2];
      G4int         tmpareacode[2] = {sOutside, sOutside};
      G4bool        tmpisvalid[2]  = {false, false};
      G4int i;

      for (i=0; i<2; i++) {
         G4double bminusD = - b - D;

         // protection against round off error  
         //G4double protection = 1.0e-6;
         G4double protection = 0;
         if ( b * D < 0 && std::fabs(bminusD / D) < protection ) {
            G4double acovbb = (a*c)/(b*b);
            tmpdist[i] = - c/b * ( 1 - acovbb * (1 + 2*acovbb));
         } else { 
            tmpdist[i] = factor * bminusD;
         }

         D = -D;
         tmpxx[i] = p + tmpdist[i]*v;
         
         if (validate == kValidateWithTol) {
            tmpareacode[i] = GetAreaCode(tmpxx[i]);
            if (!IsOutside(tmpareacode[i])) {
               if (tmpdist[i] >= 0) tmpisvalid[i] = true;
               continue;
            }
         } else if (validate == kValidateWithoutTol) {
            tmpareacode[i] = GetAreaCode(tmpxx[i], false);
            if (IsInside(tmpareacode[i])) {
               if (tmpdist[i] >= 0) tmpisvalid[i] = true;
               continue;
            }
         } else { // kDontValidate
            // we must choose x(rho,z) = rho(z=0) > 0
            if (tmpxx[i].x() > 0) {
               tmpareacode[i] = sInside;
               if (tmpdist[i] >= 0) tmpisvalid[i] = true;
               continue;
            } else {
               tmpdist[i] = kInfinity;
               continue;
            }                     
         }
      }
      
      if (tmpdist[0] <= tmpdist[1]) {
         distance[0] = tmpdist[0];
         distance[1] = tmpdist[1];
         xx[0]       = tmpxx[0];
         xx[1]       = tmpxx[1];
         gxx[0]      = ComputeGlobalPoint(tmpxx[0]);
         gxx[1]      = ComputeGlobalPoint(tmpxx[1]);
         areacode[0] = tmpareacode[0];
         areacode[1] = tmpareacode[1];
         isvalid[0]  = tmpisvalid[0];
         isvalid[1]  = tmpisvalid[1];
      } else {
         distance[0] = tmpdist[1];
         distance[1] = tmpdist[0];
         xx[0]       = tmpxx[1];
         xx[1]       = tmpxx[0];
         gxx[0]      = ComputeGlobalPoint(tmpxx[1]);
         gxx[1]      = ComputeGlobalPoint(tmpxx[0]);
         areacode[0] = tmpareacode[1];
         areacode[1] = tmpareacode[0];
         isvalid[0]  = tmpisvalid[1];
         isvalid[1]  = tmpisvalid[0];
      }
         
      fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
                                     isvalid[0], 2, validate, &gp, &gv);
      fCurStatWithV.SetCurrentStatus(1, gxx[1], distance[1], areacode[1],
                                     isvalid[1], 2, validate, &gp, &gv);

      // protection against roundoff error

      for (G4int k=0; k<2; k++) {
         if (!isvalid[k]) continue;

         G4ThreeVector xxonsurface(xx[k].x(), fKappa * std::fabs(xx[k].x())
                                              * xx[k].z() , xx[k].z());
         G4double      deltaY  =  (xx[k] - xxonsurface).mag();

         if ( deltaY > 0.5*kCarTolerance ) {

           G4int maxcount = 10;
           G4int l;
           G4double      lastdeltaY = deltaY; 
           for (l=0; l<maxcount; l++) {
             G4ThreeVector surfacenormal = GetNormal(xxonsurface); 
             distance[k] = DistanceToPlaneWithV(p, v, xxonsurface,
                                                surfacenormal, xx[k]);
             deltaY      = (xx[k] - xxonsurface).mag();
             if (deltaY > lastdeltaY) {
               
             }
             gxx[k]      = ComputeGlobalPoint(xx[k]);

               if (deltaY <= 0.5*kCarTolerance) {
            
                 break;
               }
               xxonsurface.set(xx[k].x(),
                               fKappa * std::fabs(xx[k].x()) * xx[k].z(),
                               xx[k].z());
            }
            if (l == maxcount) {
               std::ostringstream message;
               message << "Exceeded maxloop count!" << G4endl
                      << "        maxloop count " << maxcount;
               G4Exception("G4TwistTubsFlatSide::DistanceToSurface(p,v)",
                           "GeomSolids0003",  FatalException, message);
            }
         }

      }
      vout = 2;
   } else {
      // if D<0, no solution
      // if D=0, just grazing the surfaces, return kInfinity

      fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
                                     isvalid[0], 0, validate, &gp, &gv);
   }

   return vout;
}

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

virtual

Implements G4VTwistSurface.

Definition at line 454 of file G4TwistTubsSide.cc.

References G4VTwistSurface::AmIOnLeftSide(), G4VTwistSurface::ComputeGlobalPoint(), G4VTwistSurface::ComputeLocalPoint(), G4VTwistSurface::DistanceToBoundary(), G4VTwistSurface::DistanceToLine(), G4VTwistSurface::DistanceToPlane(), G4VTwistSurface::fAxis, G4VTwistSurface::fCurStat, G4VTwistSurface::fCurStatWithV, G4VTwistSurface::CurrentStatus::GetAreacode(), G4VTwistSurface::GetBoundaryAtPZ(), G4VTwistSurface::GetBoundaryParameters(), G4VTwistSurface::CurrentStatus::GetDistance(), G4VTwistSurface::CurrentStatus::GetNXX(), CLHEP::Hep3Vector::getRho(), G4VTwistSurface::CurrentStatus::GetXX(), G4VTwistSurface::CurrentStatus::IsDone(), G4VTwistSurface::kCarTolerance, G4VTwistSurface::kDontValidate, kXAxis, kZAxis, G4VTwistSurface::CurrentStatus::ResetfDone(), G4VTwistSurface::sAxis0, G4VTwistSurface::sAxisMax, G4VTwistSurface::sAxisMin, CLHEP::Hep3Vector::set(), G4VTwistSurface::CurrentStatus::SetCurrentStatus(), G4VTwistSurface::sInside, G4VTwistSurface::sOutside, mcscore::test(), CLHEP::Hep3Vector::x(), CLHEP::Hep3Vector::y(), z, and CLHEP::Hep3Vector::z().

{  
   fCurStat.ResetfDone(kDontValidate, &gp);
   G4int i = 0;
   if (fCurStat.IsDone()) {
      for (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 (i=0; i<2; i++) {
         distance[i] = kInfinity;
         areacode[i] = sOutside;
         gxx[i].set(kInfinity, kInfinity, kInfinity);
      }
   }
   
   const G4double halftol = 0.5 * kCarTolerance; 

   G4ThreeVector  p       = ComputeLocalPoint(gp);
   G4ThreeVector  xx;
   G4int          parity  = (fKappa >= 0 ? 1 : -1);
 
   // 
   // special case! 
   // If p is on surface, or
   // p is on z-axis, 
   // return here immediatery.
   //
   
   G4ThreeVector  lastgxx[2];
   for (i=0; i<2; i++) {
      lastgxx[i] = fCurStatWithV.GetXX(i);
   } 

   if  ((gp - lastgxx[0]).mag() < halftol
     || (gp - lastgxx[1]).mag() < halftol) { 
      // last winner, or last poststep point is on the surface.
      xx = p;
      distance[0] = 0;
      gxx[0] = gp;

      G4bool isvalid = true;
      fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
                             isvalid, 1, kDontValidate, &gp);
      return 1;
   }
          
   if (p.getRho() == 0) { 
      // p is on z-axis. Namely, p is on twisted surface (invalid area).
      // We must return here, however, returning distance to x-minimum
      // boundary is better than return 0-distance.
      //
      G4bool isvalid = true;
      if (fAxis[0] == kXAxis && fAxis[1] == kZAxis) {
         distance[0] = DistanceToBoundary(sAxis0 & sAxisMin, xx, p);
         areacode[0] = sInside;
      } else {
         distance[0] = 0;
         xx.set(0., 0., 0.);
      }
      gxx[0] = ComputeGlobalPoint(xx);
      fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
                                isvalid, 0, kDontValidate, &gp);
      return 1;
   } 

   // 
   // special case end
   //

   // set corner points of quadrangle try area ...

   G4ThreeVector A;  // foot of normal from p to boundary of sAxis0 & sAxisMin
   G4ThreeVector C;  // foot of normal from p to boundary of sAxis0 & sAxisMax
   G4ThreeVector B;       // point on boundary sAxis0 & sAxisMax at z = A.z()
   G4ThreeVector D;       // point on boundary sAxis0 & sAxisMin at z = C.z()

   // G4double      distToA; // distance from p to A
   DistanceToBoundary(sAxis0 & sAxisMin, A, p);
   // G4double      distToC; // distance from p to C 
   DistanceToBoundary(sAxis0 & sAxisMax, C, p);
   
   // is p.z between a.z and c.z?
   // p.z must be bracketed a.z and c.z.
   if (A.z() > C.z()) {
      if (p.z() > A.z()) {
         A = GetBoundaryAtPZ(sAxis0 & sAxisMin, p);
      } else if (p.z() < C.z()) {
         C = GetBoundaryAtPZ(sAxis0 & sAxisMax, p);
      }
   } else {
      if (p.z() > C.z()) {
         C = GetBoundaryAtPZ(sAxis0 & sAxisMax, p);
      } else if (p.z() < A.z()) {
         A = GetBoundaryAtPZ(sAxis0 & sAxisMin, p);
      }
   }

   G4ThreeVector  d[2];     // direction vectors of boundary
   G4ThreeVector  x0[2];    // foot of normal from line to p 
   G4int          btype[2]; // boundary type

   for (i=0; i<2; i++) {
      if (i == 0) {
         GetBoundaryParameters((sAxis0 & sAxisMax), d[i], x0[i], btype[i]);
         B = x0[i] + ((A.z() - x0[i].z()) / d[i].z()) * d[i]; 
         // x0 + t*d , d is direction unit vector.
      } else {
         GetBoundaryParameters((sAxis0 & sAxisMin), d[i], x0[i], btype[i]);
         D = x0[i] + ((C.z() - x0[i].z()) / d[i].z()) * d[i]; 
      }
   }

   // In order to set correct diagonal, swap A and D, C and B if needed.  
   G4ThreeVector pt(p.x(), p.y(), 0.);
   G4double      rc = std::fabs(p.x());
   G4ThreeVector surfacevector(rc, rc * fKappa * p.z(), 0.); 
   G4int         pside = AmIOnLeftSide(pt, surfacevector); 
   G4double      test  = (A.z() - C.z()) * parity * pside;  

   if (test == 0) {
      if (pside == 0) {
         // p is on surface.
         xx = p;
         distance[0] = 0;
         gxx[0] = gp;

         G4bool isvalid = true;
         fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
                                isvalid, 1, kDontValidate, &gp);
         return 1;
      } else {
         // A.z = C.z(). return distance to line.
         d[0] = C - A;
         distance[0] = DistanceToLine(p, A, d[0], xx);
         areacode[0] = sInside;
         gxx[0] = ComputeGlobalPoint(xx);
         G4bool isvalid = true;
         fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
                                isvalid, 1, kDontValidate, &gp);
         return 1;
      } 

   } else if (test < 0) {

      // wrong diagonal. vector AC is crossing the surface!  
      // swap A and D, C and B
      G4ThreeVector tmp;
      tmp = A;
      A = D;
      D = tmp;
      tmp = C;
      C = B;
      B = tmp; 

   } else {
      // correct diagonal. nothing to do.  
   }

   // Now, we chose correct diaglnal.
   // First try. divide quadrangle into double triangle by diagonal and 
   // calculate distance to both surfaces.

   G4ThreeVector xxacb;   // foot of normal from plane ACB to p
   G4ThreeVector nacb;    // normal of plane ACD
   G4ThreeVector xxcad;   // foot of normal from plane CAD to p
   G4ThreeVector ncad;    // normal of plane CAD
   G4ThreeVector AB(A.x(), A.y(), 0);
   G4ThreeVector DC(C.x(), C.y(), 0);

   G4double distToACB = G4VTwistSurface::DistanceToPlane(p, A, C-A, AB, xxacb, nacb) * parity;
   G4double distToCAD = G4VTwistSurface::DistanceToPlane(p, C, C-A, DC, xxcad, ncad) * parity;

   // if calculated distance = 0, return  

   if (std::fabs(distToACB) <= halftol || std::fabs(distToCAD) <= halftol) {
      xx = (std::fabs(distToACB) < std::fabs(distToCAD) ? xxacb : xxcad); 
      areacode[0] = sInside;
      gxx[0] = ComputeGlobalPoint(xx);
      distance[0] = 0;
      G4bool isvalid = true;
      fCurStat.SetCurrentStatus(0, gxx[0], distance[0] , areacode[0],
                                isvalid, 1, kDontValidate, &gp);
      return 1;
   }
   
   if (distToACB * distToCAD > 0 && distToACB < 0) {
      // both distToACB and distToCAD are negative.
      // divide quadrangle into double triangle by diagonal
      G4ThreeVector normal;
      distance[0] = DistanceToPlane(p, A, B, C, D, parity, xx, normal);
   } else {
      if (distToACB * distToCAD > 0) {
         // both distToACB and distToCAD are positive.
         // Take smaller one.
         if (distToACB <= distToCAD) {
            distance[0] = distToACB;
            xx   = xxacb;
         } else {
            distance[0] = distToCAD;
            xx   = xxcad;
         }
      } else {
         // distToACB * distToCAD is negative.
         // take positive one
         if (distToACB > 0) {
            distance[0] = distToACB;
            xx   = xxacb;
         } else {
            distance[0] = distToCAD;
            xx   = xxcad;
         }
      }
      
   }
   areacode[0] = sInside;
   gxx[0]      = ComputeGlobalPoint(xx);
   G4bool isvalid = true;
   fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
                             isvalid, 1, kDontValidate, &gp);
   return 1;
}

G4double G4TwistTubsSide::GetBoundaryMax ( G4double  phi )

inlinevirtual

Implements G4VTwistSurface.

Definition at line 181 of file G4TwistTubsSide.hh.

References G4VTwistSurface::fAxisMax.

Referenced by GetFacets().

{
  return  fAxisMax[0] ;  // outer radius at z = 0
}

G4double G4TwistTubsSide::GetBoundaryMin ( G4double  phi )

inlinevirtual

Implements G4VTwistSurface.

Definition at line 175 of file G4TwistTubsSide.hh.

References G4VTwistSurface::fAxisMin.

Referenced by GetFacets().

{
  return  fAxisMin[0] ;  // inner radius at z = 0
}

void G4TwistTubsSide::GetFacets ( G4int  m, G4int  n, G4double  xyz[][3], G4int  faces[][4], G4int  iside  )

virtual

Implements G4VTwistSurface.

Definition at line 953 of file G4TwistTubsSide.cc.

References G4VTwistSurface::fAxisMax, G4VTwistSurface::fAxisMin, G4VTwistSurface::fHandedness, GetBoundaryMax(), GetBoundaryMin(), G4VTwistSurface::GetEdgeVisibility(), G4VTwistSurface::GetFace(), G4VTwistSurface::GetNode(), n, SurfacePoint(), test::x, CLHEP::Hep3Vector::x(), CLHEP::Hep3Vector::y(), z, and CLHEP::Hep3Vector::z().

{

  G4double z ;     // the two parameters for the surface equation
  G4double x,xmin,xmax ;

  G4ThreeVector p ;  // a point on the surface, given by (z,u)

  G4int nnode ;
  G4int nface ;

  // calculate the (n-1)*(k-1) vertices

  G4int i,j ;

  for ( i = 0 ; i<n ; i++ )
  {

    z = fAxisMin[1] + i*(fAxisMax[1]-fAxisMin[1])/(n-1) ;

    for ( j = 0 ; j<k ; j++ ) {

      nnode = GetNode(i,j,k,n,iside) ;

      xmin = GetBoundaryMin(z) ; 
      xmax = GetBoundaryMax(z) ;

      if (fHandedness < 0) { 
        x = xmin + j*(xmax-xmin)/(k-1) ;
      } else {               
        x = xmax - j*(xmax-xmin)/(k-1) ;
      }

      p = SurfacePoint(x,z,true) ;  // surface point in global coord.system

      xyz[nnode][0] = p.x() ;
      xyz[nnode][1] = p.y() ;
      xyz[nnode][2] = p.z() ;

      if ( i<n-1 && j<k-1 ) {   // clock wise filling
        
        nface = GetFace(i,j,k,n,iside) ;

    faces[nface][0] = GetEdgeVisibility(i,j,k,n,0,1) * ( GetNode(i  ,j  ,k,n,iside)+1) ;  
    faces[nface][1] = GetEdgeVisibility(i,j,k,n,1,1) * ( GetNode(i+1,j  ,k,n,iside)+1) ;
    faces[nface][2] = GetEdgeVisibility(i,j,k,n,2,1) * ( GetNode(i+1,j+1,k,n,iside)+1) ;
    faces[nface][3] = GetEdgeVisibility(i,j,k,n,3,1) * ( GetNode(i  ,j+1,k,n,iside)+1) ;

      }
    }
  }
}

G4ThreeVector G4TwistTubsSide::GetNormal ( const G4ThreeVector &  xx, G4bool  isGlobal = false  )

virtual

Implements G4VTwistSurface.

Definition at line 125 of file G4TwistTubsSide.cc.

References G4VTwistSurface::ComputeGlobalDirection(), G4VTwistSurface::ComputeLocalPoint(), G4VTwistSurface::fCurrentNormal, G4VTwistSurface::fHandedness, G4VTwistSurface::kCarTolerance, G4VTwistSurface::G4SurfCurNormal::normal, G4VTwistSurface::G4SurfCurNormal::p, CLHEP::Hep3Vector::unit(), CLHEP::Hep3Vector::x(), and CLHEP::Hep3Vector::z().

Referenced by DistanceToSurface().

{
   // 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;
      }
   }
   
   G4ThreeVector er(1, fKappa * xx.z(), 0);
   G4ThreeVector ez(0, fKappa * xx.x(), 1);
   G4ThreeVector normal = fHandedness*(er.cross(ez));

   if (isGlobal) {
      fCurrentNormal.normal = ComputeGlobalDirection(normal.unit());
   } else {
      fCurrentNormal.normal = normal.unit();
   }
   return fCurrentNormal.normal;
}

G4double G4TwistTubsSide::GetSurfaceArea ( )

inlinevirtual

Implements G4VTwistSurface.

Definition at line 187 of file G4TwistTubsSide.hh.

References G4VTwistSurface::fAxisMax, and G4VTwistSurface::fAxisMin.

{
  // approximation only
  return ( fAxisMax[0] - fAxisMin[0] ) * ( fAxisMax[1] - fAxisMin[1] ) ;
}

G4ThreeVector G4TwistTubsSide::ProjectAtPXPZ ( const G4ThreeVector &  p, G4bool  isglobal = false  ) const

inline

Definition at line 149 of file G4TwistTubsSide.hh.

References G4VTwistSurface::fRot, G4VTwistSurface::fTrans, CLHEP::HepRotation::inverse(), CLHEP::Hep3Vector::x(), and CLHEP::Hep3Vector::z().

{
  // Get Rho at p.z() on Hyperbolic Surface.
  G4ThreeVector tmpp;
  if (isglobal) {
     tmpp = fRot.inverse()*p - fTrans;
  } else {
     tmpp = p;
  }
  G4ThreeVector xx(p.x(), p.x() * fKappa * p.z(), p.z());
  if (isglobal) { return (fRot * xx + fTrans); }
  return xx;
}

G4ThreeVector G4TwistTubsSide::SurfacePoint ( G4double  x, G4double  z, G4bool  isGlobal = false  )

inlinevirtual

Implements G4VTwistSurface.

Definition at line 166 of file G4TwistTubsSide.hh.

References G4VTwistSurface::fRot, and G4VTwistSurface::fTrans.

Referenced by GetFacets().

{
  G4ThreeVector SurfPoint( x , x * fKappa * z , z ) ;

  if (isGlobal) { return (fRot * SurfPoint + fTrans); }
  return SurfPoint;
}

---

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

• G4TwistTubsSide.hh
• G4TwistTubsSide.cc