#include <G4TwistTrapParallelSide.hh>
Inheritance diagram for G4TwistTrapParallelSide:
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__ &) |
Definition at line 51 of file G4TwistTrapParallelSide.hh.
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] |
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 }
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 }