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__ &)

---

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, and kZAxis.

  : 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__ &

)

Definition at line 109 of file G4TwistTubsSide.cc.

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

---

Member Function Documentation

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::GetBoundaryAtPZ(), G4VTwistSurface::GetBoundaryParameters(), G4VTwistSurface::kCarTolerance, G4VTwistSurface::kDontValidate, kXAxis, kZAxis, G4VTwistSurface::sAxis0, G4VTwistSurface::sAxisMax, G4VTwistSurface::sAxisMin, G4VTwistSurface::sInside, and G4VTwistSurface::sOutside.

{  
   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);
      }
   }
   
   static 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;
}

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

Definition at line 160 of file G4TwistTubsSide.cc.

References G4VTwistSurface::ComputeGlobalPoint(), G4VTwistSurface::ComputeLocalDirection(), G4VTwistSurface::ComputeLocalPoint(), DBL_MIN, G4VTwistSurface::DistanceToPlaneWithV(), FatalException, G4VTwistSurface::fCurStat, G4VTwistSurface::fCurStatWithV, G4endl, G4Exception(), GetNormal(), G4VTwistSurface::IsInside(), G4VTwistSurface::IsOutside(), G4VTwistSurface::kCarTolerance, G4VTwistSurface::kValidateWithoutTol, G4VTwistSurface::kValidateWithTol, G4VTwistSurface::sInside, and G4VTwistSurface::sOutside.

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

G4double G4TwistTubsSide::GetBoundaryMax

(

G4double

phi

)

[inline, virtual]

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

)

[inline, virtual]

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(), and SurfacePoint().

{

  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, and G4VTwistSurface::kCarTolerance.

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

(

)

[inline, virtual]

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, and G4VTwistSurface::fTrans.

{
  // 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	,
		G4double	,
		G4bool	isGlobal = false
	)			[inline, virtual]

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

---