G4TwistTrapParallelSide Class Reference

#include <G4TwistTrapParallelSide.hh>

Inheritance diagram for G4TwistTrapParallelSide:

G4VTwistSurface

Public Member Functions

 G4TwistTrapParallelSide (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 ~G4TwistTrapParallelSide ()
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[])
 G4TwistTrapParallelSide (__void__ &)

Detailed Description

Definition at line 51 of file G4TwistTrapParallelSide.hh.


Constructor & Destructor Documentation

G4TwistTrapParallelSide::G4TwistTrapParallelSide ( 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 51 of file G4TwistTrapParallelSide.cc.

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

00065   : G4VTwistSurface(name)
00066 {  
00067   
00068   fAxis[0]    = kXAxis; // in local coordinate system
00069   fAxis[1]    = kZAxis;
00070   fAxisMin[0] = -kInfinity ;  // X Axis boundary
00071   fAxisMax[0] = kInfinity ;   //   depends on z !!
00072   fAxisMin[1] = -pDz ;      // Z Axis boundary
00073   fAxisMax[1] = pDz ;
00074   
00075   fDx1  = pDx1 ;
00076   fDx2  = pDx2 ;
00077   fDx3  = pDx3 ;
00078   fDx4  = pDx4 ;
00079 
00080   fDy1   = pDy1 ;
00081   fDy2   = pDy2 ;
00082 
00083   fDz   = pDz ;
00084 
00085   fAlph = pAlph  ;
00086   fTAlph = std::tan(fAlph) ;
00087 
00088   fTheta = pTheta ;
00089   fPhi   = pPhi ;
00090 
00091  // precalculate frequently used parameters
00092   fDx4plus2  = fDx4 + fDx2 ;
00093   fDx4minus2 = fDx4 - fDx2 ;
00094   fDx3plus1  = fDx3 + fDx1 ; 
00095   fDx3minus1 = fDx3 - fDx1 ;
00096   fDy2plus1  = fDy2 + fDy1 ;
00097   fDy2minus1 = fDy2 - fDy1 ;
00098 
00099   fa1md1 = 2*fDx2 - 2*fDx1  ; 
00100   fa2md2 = 2*fDx4 - 2*fDx3 ;
00101 
00102   fPhiTwist = PhiTwist ;     // dphi
00103   fAngleSide = AngleSide ;  // 0,90,180,270 deg
00104 
00105   fdeltaX = 2 * fDz * std::tan(fTheta) * std::cos(fPhi)  ;  // dx in surface equation
00106   fdeltaY = 2 * fDz * std::tan(fTheta) * std::sin(fPhi)  ;  // dy in surface equation
00107   
00108   fRot.rotateZ( AngleSide ) ; 
00109   
00110   fTrans.set(0, 0, 0);  // No Translation
00111   fIsValidNorm = false;
00112   
00113   SetCorners() ;
00114   SetBoundaries() ;
00115 
00116 }

G4TwistTrapParallelSide::~G4TwistTrapParallelSide (  )  [virtual]

Definition at line 135 of file G4TwistTrapParallelSide.cc.

00136 {
00137 }

G4TwistTrapParallelSide::G4TwistTrapParallelSide ( __void__ &   ) 

Definition at line 122 of file G4TwistTrapParallelSide.cc.

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


Member Function Documentation

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

Implements G4VTwistSurface.

Definition at line 655 of file G4TwistTrapParallelSide.cc.

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

00659 {  
00660   // to do
00661 
00662   static const G4double ctol = 0.5 * kCarTolerance;
00663 
00664   fCurStat.ResetfDone(kDontValidate, &gp);
00665 
00666    if (fCurStat.IsDone()) {
00667       G4int i;
00668       for (i=0; i<fCurStat.GetNXX(); i++) {
00669          gxx[i] = fCurStat.GetXX(i);
00670          distance[i] = fCurStat.GetDistance(i);
00671          areacode[i] = fCurStat.GetAreacode(i);
00672       }
00673       return fCurStat.GetNXX();
00674    } else {
00675       // initialize
00676       G4int i;
00677       for (i=0; i<G4VSURFACENXX; i++) {
00678          distance[i] = kInfinity;
00679          areacode[i] = sOutside;
00680          gxx[i].set(kInfinity, kInfinity, kInfinity);
00681       }
00682    }
00683    
00684    G4ThreeVector p = ComputeLocalPoint(gp);
00685    G4ThreeVector xx;  // intersection point
00686    G4ThreeVector xxonsurface ; // interpolated intersection point 
00687 
00688    // the surfacenormal at that surface point
00689    G4double phiR = 0  ; // 
00690    G4double uR = 0 ;
00691 
00692    G4ThreeVector surfacenormal ; 
00693    G4double deltaX ;
00694    
00695    G4int maxint = 20 ;
00696 
00697    for ( G4int i = 1 ; i<maxint ; i++ ) {
00698 
00699      xxonsurface = SurfacePoint(phiR,uR) ;
00700      surfacenormal = NormAng(phiR,uR) ;
00701      distance[0] = DistanceToPlane(p, xxonsurface, surfacenormal, xx); // new XX
00702      deltaX = ( xx - xxonsurface ).mag() ; 
00703 
00704 #ifdef G4TWISTDEBUG
00705      G4cout << "i = " << i << ", distance = " << distance[0] << ", " << deltaX << G4endl ;
00706      G4cout << "X = " << xx << G4endl ;
00707 #endif
00708 
00709      // the new point xx is accepted and phi/psi replaced
00710      GetPhiUAtX(xx, phiR, uR) ;
00711      
00712      if ( deltaX <= ctol ) { break ; }
00713 
00714    }
00715 
00716    // check validity of solution ( valid phi,psi ) 
00717 
00718    G4double halfphi = 0.5*fPhiTwist ;
00719    G4double uMax = GetBoundaryMax(phiR) ;
00720    G4double uMin = GetBoundaryMin(phiR) ;
00721 
00722    if (  phiR > halfphi ) phiR =  halfphi ;
00723    if ( phiR < -halfphi ) phiR = -halfphi ;
00724    if ( uR > uMax ) uR = uMax ;
00725    if ( uR < uMin ) uR = uMin ;
00726 
00727    xxonsurface = SurfacePoint(phiR,uR) ;
00728    distance[0] = (  p - xx ).mag() ;
00729    if ( distance[0] <= ctol ) { distance[0] = 0 ; } 
00730 
00731    // end of validity 
00732 
00733 #ifdef G4TWISTDEBUG
00734    G4cout << "refined solution "  << phiR << " , " << uR << " , " <<  G4endl ;
00735    G4cout << "distance = " << distance[0] << G4endl ;
00736    G4cout << "X = " << xx << G4endl ;
00737 #endif
00738 
00739    G4bool isvalid = true;
00740    gxx[0]      = ComputeGlobalPoint(xx);
00741    
00742 #ifdef G4TWISTDEBUG
00743    G4cout << "intersection Point found: " << gxx[0] << G4endl ;
00744    G4cout << "distance = " << distance[0] << G4endl ;
00745 #endif
00746 
00747    fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0],
00748                             isvalid, 1, kDontValidate, &gp);
00749    return 1;
00750 }

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

Definition at line 188 of file G4TwistTrapParallelSide.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.

00195 {
00196 
00197   static const G4double ctol = 0.5 * kCarTolerance;
00198   static const G4double pihalf = pi/2 ;
00199 
00200   G4bool IsParallel = false ;
00201   G4bool IsConverged =  false ;
00202 
00203   G4int nxx = 0 ;  // number of physical solutions
00204 
00205   fCurStatWithV.ResetfDone(validate, &gp, &gv);
00206 
00207   if (fCurStatWithV.IsDone()) {
00208     G4int i;
00209     for (i=0; i<fCurStatWithV.GetNXX(); i++) {
00210       gxx[i] = fCurStatWithV.GetXX(i);
00211       distance[i] = fCurStatWithV.GetDistance(i);
00212       areacode[i] = fCurStatWithV.GetAreacode(i);
00213       isvalid[i]  = fCurStatWithV.IsValid(i);
00214     }
00215     return fCurStatWithV.GetNXX();
00216   } else {
00217    
00218    // initialize
00219     G4int i;
00220     for (i=0; i<G4VSURFACENXX ; i++) {
00221       distance[i] = kInfinity;
00222       areacode[i] = sOutside;
00223       isvalid[i]  = false;
00224       gxx[i].set(kInfinity, kInfinity, kInfinity);
00225     }
00226   }
00227 
00228   G4ThreeVector p = ComputeLocalPoint(gp);
00229   G4ThreeVector v = ComputeLocalDirection(gv);
00230   
00231 #ifdef G4TWISTDEBUG
00232   G4cout << "Local point p = " << p << G4endl ;
00233   G4cout << "Local direction v = " << v << G4endl ; 
00234 #endif
00235 
00236   G4double phi,u ;  // parameters
00237 
00238   // temporary variables
00239 
00240   G4double      tmpdist = kInfinity ;
00241   G4ThreeVector tmpxx;
00242   G4int         tmpareacode = sOutside ;
00243   G4bool        tmpisvalid  = false ;
00244 
00245   std::vector<Intersection> xbuf ;
00246   Intersection xbuftmp ;
00247   
00248   // prepare some variables for the intersection finder
00249 
00250   G4double L = 2*fDz ;
00251 
00252   G4double phixz = fPhiTwist * ( p.x() * v.z() - p.z() * v.x() ) ;
00253   G4double phiyz = fPhiTwist * ( p.y() * v.z() - p.z() * v.y() ) ;
00254 
00255   // special case vz = 0
00256 
00257   if ( v.z() == 0. ) {         
00258 
00259     if ( std::fabs(p.z()) <= L ) {     // intersection possible in z
00260       
00261       phi = p.z() * fPhiTwist / L ;  // phi is determined by the z-position 
00262 
00263       u = (2*(fdeltaY*phi*v.x() - fPhiTwist*p.y()*v.x() - fdeltaX*phi*v.y()
00264               + fPhiTwist*p.x()*v.y()) + (fDy2plus1*fPhiTwist
00265               + 2*fDy2minus1*phi)*(v.x()*std::cos(phi) + v.y()*std::sin(phi)))
00266         / (2.* fPhiTwist*(v.y()*std::cos(phi) - v.x()*std::sin(phi)));
00267 
00268       xbuftmp.phi = phi ;
00269       xbuftmp.u = u ;
00270       xbuftmp.areacode = sOutside ;
00271       xbuftmp.distance = kInfinity ;
00272       xbuftmp.isvalid = false ;
00273 
00274       xbuf.push_back(xbuftmp) ;  // store it to xbuf
00275 
00276     }
00277 
00278     else {                        // no intersection possible
00279 
00280       distance[0] = kInfinity;
00281       gxx[0].set(kInfinity,kInfinity,kInfinity);
00282       isvalid[0] = false ;
00283       areacode[0] = sOutside ;
00284       fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0],
00285                                      areacode[0], isvalid[0],
00286                                      0, validate, &gp, &gv);
00287       
00288       return 0;
00289 
00290 
00291     }  // end std::fabs(p.z() <= L 
00292 
00293   } // end v.z() == 0
00294   
00295 
00296   // general solution for non-zero vz
00297 
00298   else {
00299 
00300     G4double c[9],srd[8],si[8] ;  
00301 
00302     c[8] = -3600*(-2*phiyz + fDy2plus1*fPhiTwist*v.z()) ;
00303     c[7] = -7200*(phixz - 2*fDz*v.y() + (fdeltaY + fDy2minus1)*v.z()) ;
00304     c[6] = 120*(-52*phiyz - 120*fDz*v.x() + 60*fdeltaX*v.z() + 11*fDy2plus1*fPhiTwist*v.z()) ;
00305     c[5] = 240*(16*phixz - 52*fDz*v.y() + 26*fdeltaY*v.z() + 11*fDy2minus1*v.z()) ;
00306     c[4] = 12*(127*phiyz + 640*fDz*v.x() - 320*fdeltaX*v.z() + 4*fDy2plus1*fPhiTwist*v.z()) ;
00307     c[3] = -404*phixz + 3048*fDz*v.y() - 1524*fdeltaY*v.z() + 96*fDy2minus1*v.z() ;
00308     c[2] = -72*phiyz + 404*(-2*fDz*v.x() + fdeltaX*v.z()) ;
00309     c[1] = 12*(phixz - 12*fDz*v.y() + 6*fdeltaY*v.z()) ;
00310     c[0] = 24*fDz*v.x() - 12*fdeltaX*v.z() ;
00311 
00312 
00313 #ifdef G4TWISTDEBUG
00314     G4cout << "coef = " << c[0] << " " 
00315            <<  c[1] << " "  
00316            <<  c[2] << " "  
00317            <<  c[3] << " "  
00318            <<  c[4] << " "  
00319            <<  c[5] << " "  
00320            <<  c[6] << " "  
00321            <<  c[7] << " "  
00322            <<  c[8] << G4endl ;
00323 #endif    
00324 
00325     G4JTPolynomialSolver trapEq ;
00326     G4int num = trapEq.FindRoots(c,8,srd,si);
00327   
00328 
00329     for (G4int i = 0 ; i<num ; i++ ) {  // loop over all mathematical solutions
00330       if ( si[i]==0.0 ) {  // only real solutions
00331 #ifdef G4TWISTDEBUG
00332         G4cout << "Solution " << i << " : " << srd[i] << G4endl ;
00333 #endif
00334         phi = std::fmod(srd[i] , pihalf)  ;
00335 
00336         u = (1/std::cos(phi)*(2*phixz + 4*fDz*phi*v.x() - 2*fdeltaX*phi*v.z() + (fDy2plus1*fPhiTwist + 2*fDy2minus1*phi)*v.z()* std::sin(phi)))/(2.*fPhiTwist*v.z()) ;
00337 
00338         xbuftmp.phi = phi ;
00339         xbuftmp.u = u ;
00340         xbuftmp.areacode = sOutside ;
00341         xbuftmp.distance = kInfinity ;
00342         xbuftmp.isvalid = false ;
00343         
00344         xbuf.push_back(xbuftmp) ;  // store it to xbuf
00345       
00346 #ifdef G4TWISTDEBUG
00347         G4cout << "solution " << i << " = " << phi << " , " << u  << G4endl ;
00348 #endif
00349 
00350       }  // end if real solution
00351     }  // end loop i
00352     
00353   }    // end general case
00354 
00355 
00356   nxx = xbuf.size() ;  // save the number of  solutions
00357 
00358   G4ThreeVector xxonsurface  ;       // point on surface
00359   G4ThreeVector surfacenormal  ;     // normal vector  
00360   G4double deltaX  ;                 // distance between intersection point and point on surface
00361   G4double theta  ;                  // angle between track and surfacenormal
00362   G4double factor ;                  // a scaling factor
00363   G4int maxint = 30 ;                // number of iterations
00364 
00365 
00366   for ( size_t k = 0 ; k<xbuf.size() ; k++ ) {
00367 
00368 #ifdef G4TWISTDEBUG
00369     G4cout << "Solution " << k << " : " 
00370            << "reconstructed phiR = " << xbuf[k].phi
00371            << ", uR = " << xbuf[k].u << G4endl ; 
00372 #endif
00373     
00374     phi = xbuf[k].phi ;  // get the stored values for phi and u
00375     u = xbuf[k].u ;
00376 
00377     IsConverged = false ;   // no convergence at the beginning
00378     
00379     for ( G4int i = 1 ; i<maxint ; i++ ) {
00380       
00381       xxonsurface = SurfacePoint(phi,u) ;
00382       surfacenormal = NormAng(phi,u) ;
00383       tmpdist = DistanceToPlaneWithV(p, v, xxonsurface, surfacenormal, tmpxx); 
00384       deltaX = ( tmpxx - xxonsurface ).mag() ; 
00385       theta = std::fabs(std::acos(v*surfacenormal) - pihalf) ;
00386       if ( theta < 0.001 ) { 
00387         factor = 50 ;
00388         IsParallel = true ;
00389       }
00390       else {
00391         factor = 1 ;
00392       }
00393 
00394 #ifdef G4TWISTDEBUG
00395       G4cout << "Step i = " << i << ", distance = " << tmpdist << ", " << deltaX << G4endl ;
00396       G4cout << "X = " << tmpxx << G4endl ;
00397 #endif
00398       
00399       GetPhiUAtX(tmpxx, phi, u) ; // the new point xx is accepted and phi/u replaced
00400       
00401 #ifdef G4TWISTDEBUG
00402       G4cout << "approximated phi = " << phi << ", u = " << u << G4endl ; 
00403 #endif
00404       
00405       if ( deltaX <= factor*ctol ) { IsConverged = true ; break ; }
00406       
00407     }  // end iterative loop (i)
00408     
00409 
00410     // new code  21.09.05 O.Link
00411     if ( std::fabs(tmpdist)<ctol ) tmpdist = 0 ; 
00412 
00413 
00414 #ifdef G4TWISTDEBUG
00415     G4cout << "refined solution "  << phi << " , " << u  <<  G4endl ;
00416     G4cout << "distance = " << tmpdist << G4endl ;
00417     G4cout << "local X = " << tmpxx << G4endl ;
00418 #endif
00419     
00420     tmpisvalid = false ;  // init 
00421 
00422     if ( IsConverged ) {
00423       
00424       if (validate == kValidateWithTol) {
00425         tmpareacode = GetAreaCode(tmpxx);
00426         if (!IsOutside(tmpareacode)) {
00427           if (tmpdist >= 0) tmpisvalid = true;
00428         }
00429       } else if (validate == kValidateWithoutTol) {
00430         tmpareacode = GetAreaCode(tmpxx, false);
00431         if (IsInside(tmpareacode)) {
00432           if (tmpdist >= 0) tmpisvalid = true;
00433         }
00434       } else { // kDontValidate
00435         G4Exception("G4TwistTrapParallelSide::DistanceToSurface()",
00436                     "GeomSolids0001", FatalException,
00437                     "Feature NOT implemented !");
00438       }
00439 
00440     } 
00441     else {
00442       tmpdist = kInfinity;     // no convergence after 10 steps 
00443       tmpisvalid = false ;     // solution is not vaild
00444     }  
00445 
00446 
00447     // store the found values 
00448     xbuf[k].xx = tmpxx ;
00449     xbuf[k].distance = tmpdist ;
00450     xbuf[k].areacode = tmpareacode ;
00451     xbuf[k].isvalid = tmpisvalid ;
00452 
00453 
00454   }  // end loop over physical solutions (variable k)
00455 
00456 
00457   std::sort(xbuf.begin() , xbuf.end(), DistanceSort ) ;  // sorting
00458 
00459 #ifdef G4TWISTDEBUG
00460   G4cout << G4endl << "list xbuf after sorting : " << G4endl ;
00461   G4cout << G4endl << G4endl ;
00462 #endif
00463 
00464 
00465   // erase identical intersection (within kCarTolerance) 
00466   xbuf.erase( std::unique(xbuf.begin(), xbuf.end() , EqualIntersection ) , xbuf.end() ) ;
00467 
00468 
00469   // add guesses
00470 
00471   G4int nxxtmp = xbuf.size() ;
00472 
00473   if ( nxxtmp<2 || IsParallel  ) {
00474 
00475     // positive end
00476 #ifdef G4TWISTDEBUG
00477     G4cout << "add guess at +z/2 .. " << G4endl ;
00478 #endif
00479 
00480     phi = fPhiTwist/2 ;
00481     u   =  0 ;
00482     
00483     xbuftmp.phi = phi ;
00484     xbuftmp.u = u ;
00485     xbuftmp.areacode = sOutside ;
00486     xbuftmp.distance = kInfinity ;
00487     xbuftmp.isvalid = false ;
00488     
00489     xbuf.push_back(xbuftmp) ;  // store it to xbuf
00490 
00491 
00492 #ifdef G4TWISTDEBUG
00493     G4cout << "add guess at -z/2 .. " << G4endl ;
00494 #endif
00495 
00496     phi = -fPhiTwist/2 ;
00497     u   = 0 ;
00498 
00499     xbuftmp.phi = phi ;
00500     xbuftmp.u = u ;
00501     xbuftmp.areacode = sOutside ;
00502     xbuftmp.distance = kInfinity ;
00503     xbuftmp.isvalid = false ;
00504     
00505     xbuf.push_back(xbuftmp) ;  // store it to xbuf
00506 
00507     for ( size_t k = nxxtmp ; k<xbuf.size() ; k++ ) {
00508 
00509 #ifdef G4TWISTDEBUG
00510       G4cout << "Solution " << k << " : " 
00511              << "reconstructed phiR = " << xbuf[k].phi
00512              << ", uR = " << xbuf[k].u << G4endl ; 
00513 #endif
00514       
00515       phi = xbuf[k].phi ;  // get the stored values for phi and u
00516       u   = xbuf[k].u ;
00517 
00518       IsConverged = false ;   // no convergence at the beginning
00519       
00520       for ( G4int i = 1 ; i<maxint ; i++ ) {
00521         
00522         xxonsurface = SurfacePoint(phi,u) ;
00523         surfacenormal = NormAng(phi,u) ;
00524         tmpdist = DistanceToPlaneWithV(p, v, xxonsurface, surfacenormal, tmpxx); 
00525         deltaX = ( tmpxx - xxonsurface ).mag() ; 
00526         theta = std::fabs(std::acos(v*surfacenormal) - pihalf) ;
00527         if ( theta < 0.001 ) { 
00528           factor = 50 ;    
00529         }
00530         else {
00531           factor = 1 ;
00532         }
00533         
00534 #ifdef G4TWISTDEBUG
00535         G4cout << "Step i = " << i << ", distance = " << tmpdist << ", " << deltaX << G4endl ;
00536         G4cout << "X = " << tmpxx << G4endl ;
00537 #endif
00538 
00539         GetPhiUAtX(tmpxx, phi, u) ; // the new point xx is accepted and phi/u replaced
00540       
00541 #ifdef G4TWISTDEBUG
00542         G4cout << "approximated phi = " << phi << ", u = " << u << G4endl ; 
00543 #endif
00544       
00545         if ( deltaX <= factor*ctol ) { IsConverged = true ; break ; }
00546       
00547       }  // end iterative loop (i)
00548     
00549 
00550       // new code  21.09.05 O.Link
00551       if ( std::fabs(tmpdist)<ctol ) tmpdist = 0 ; 
00552 
00553 
00554 #ifdef G4TWISTDEBUG
00555       G4cout << "refined solution "  << phi << " , " << u  <<  G4endl ;
00556       G4cout << "distance = " << tmpdist << G4endl ;
00557       G4cout << "local X = " << tmpxx << G4endl ;
00558 #endif
00559 
00560       tmpisvalid = false ;  // init 
00561 
00562       if ( IsConverged ) {
00563 
00564         if (validate == kValidateWithTol) {
00565           tmpareacode = GetAreaCode(tmpxx);
00566           if (!IsOutside(tmpareacode)) {
00567             if (tmpdist >= 0) tmpisvalid = true;
00568           }
00569         } else if (validate == kValidateWithoutTol) {
00570           tmpareacode = GetAreaCode(tmpxx, false);
00571           if (IsInside(tmpareacode)) {
00572             if (tmpdist >= 0) tmpisvalid = true;
00573           }
00574         } else { // kDontValidate
00575           G4Exception("G4TwistedBoxSide::DistanceToSurface()",
00576                       "GeomSolids0001", FatalException,
00577                       "Feature NOT implemented !");
00578         }
00579         
00580       } 
00581       else {
00582         tmpdist = kInfinity;     // no convergence after 10 steps 
00583         tmpisvalid = false ;     // solution is not vaild
00584       }  
00585         
00586         
00587       // store the found values 
00588       xbuf[k].xx = tmpxx ;
00589       xbuf[k].distance = tmpdist ;
00590       xbuf[k].areacode = tmpareacode ;
00591       xbuf[k].isvalid = tmpisvalid ;
00592 
00593 
00594     }  // end loop over physical solutions 
00595 
00596 
00597   }  // end less than 2 solutions
00598 
00599 
00600   // sort again
00601   std::sort(xbuf.begin() , xbuf.end(), DistanceSort ) ;  // sorting
00602 
00603   // erase identical intersection (within kCarTolerance) 
00604   xbuf.erase( std::unique(xbuf.begin(), xbuf.end() , EqualIntersection ) , xbuf.end() ) ;
00605 
00606 #ifdef G4TWISTDEBUG
00607   G4cout << G4endl << "list xbuf after sorting : " << G4endl ;
00608   G4cout << G4endl << G4endl ;
00609 #endif
00610 
00611   nxx = xbuf.size() ;   // determine number of solutions again.
00612 
00613   for ( size_t i = 0 ; i<xbuf.size() ; i++ ) {
00614     
00615     distance[i] = xbuf[i].distance;
00616     gxx[i]      = ComputeGlobalPoint(xbuf[i].xx);
00617     areacode[i] = xbuf[i].areacode ;
00618     isvalid[i]  = xbuf[i].isvalid ;
00619     
00620     fCurStatWithV.SetCurrentStatus(i, gxx[i], distance[i], areacode[i],
00621                                      isvalid[i], nxx, validate, &gp, &gv);
00622 
00623 #ifdef G4TWISTDEBUG
00624     G4cout << "element Nr. " << i 
00625            << ", local Intersection = " << xbuf[i].xx 
00626            << ", distance = " << xbuf[i].distance 
00627            << ", u = " << xbuf[i].u 
00628            << ", phi = " << xbuf[i].phi 
00629            << ", isvalid = " << xbuf[i].isvalid 
00630            << G4endl ;
00631 #endif
00632 
00633   }  // end for( i ) loop
00634 
00635     
00636 #ifdef G4TWISTDEBUG
00637   G4cout << "G4TwistTrapParallelSide finished " << G4endl ;
00638   G4cout << nxx << " possible physical solutions found" << G4endl ;
00639   for ( G4int k= 0 ; k< nxx ; k++ ) {
00640     G4cout << "global intersection Point found: " << gxx[k] << G4endl ;
00641     G4cout << "distance = " << distance[k] << G4endl ;
00642     G4cout << "isvalid = " << isvalid[k] << G4endl ;
00643   }
00644 #endif
00645 
00646   return nxx ;
00647     
00648 }

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

Implements G4VTwistSurface.

Definition at line 142 of file G4TwistTrapParallelSide.cc.

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

00144 {
00145    // GetNormal returns a normal vector at a surface (or very close
00146    // to surface) point at tmpxx.
00147    // If isGlobal=true, it returns the normal in global coordinate.
00148    //
00149 
00150    G4ThreeVector xx;
00151    if (isGlobal) {
00152       xx = ComputeLocalPoint(tmpxx);
00153       if ((xx - fCurrentNormal.p).mag() < 0.5 * kCarTolerance) {
00154          return ComputeGlobalDirection(fCurrentNormal.normal);
00155       }
00156    } else {
00157       xx = tmpxx;
00158       if (xx == fCurrentNormal.p) {
00159          return fCurrentNormal.normal;
00160       }
00161    }
00162 
00163    G4double phi ;
00164    G4double u ;
00165 
00166    GetPhiUAtX(xx,phi,u) ;   // phi,u for point xx close to surface 
00167 
00168    G4ThreeVector normal =  NormAng(phi,u) ;  // the normal vector at phi,u
00169 
00170 #ifdef G4TWISTDEBUG
00171    G4cout  << "normal vector = " << normal << G4endl ;
00172    G4cout << "phi = " << phi << " , u = " << u << G4endl ;
00173 #endif
00174 
00175    //    normal = normal/normal.mag() ;
00176 
00177    if (isGlobal) {
00178       fCurrentNormal.normal = ComputeGlobalDirection(normal.unit());
00179    } else {
00180       fCurrentNormal.normal = normal.unit();
00181    }
00182    return fCurrentNormal.normal;
00183 }


The documentation for this class was generated from the following files:
Generated on Mon May 27 17:53:35 2013 for Geant4 by  doxygen 1.4.7