#include <G4TwistBoxSide.hh>
Inheritance diagram for G4TwistBoxSide:
Public Member Functions | |
G4TwistBoxSide (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 | ~G4TwistBoxSide () |
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[]) |
G4TwistBoxSide (__void__ &) |
Definition at line 51 of file G4TwistBoxSide.hh.
G4TwistBoxSide::G4TwistBoxSide | ( | 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 G4TwistBoxSide.cc.
References FatalException, G4VTwistSurface::fAxis, G4VTwistSurface::fAxisMax, G4VTwistSurface::fAxisMin, G4VTwistSurface::fIsValidNorm, G4VTwistSurface::fRot, G4VTwistSurface::fTrans, G4endl, G4Exception(), G4VTwistSurface::GetName(), kYAxis, and kZAxis.
00064 : G4VTwistSurface(name) 00065 { 00066 00067 00068 fAxis[0] = kYAxis; // in local coordinate system 00069 fAxis[1] = kZAxis; 00070 fAxisMin[0] = -kInfinity ; // Y 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 ; // box 00077 fDx3 = pDx3 ; 00078 fDx4 = pDx4 ; // box 00079 00080 // this is an overhead. But the parameter naming scheme fits to the other surfaces. 00081 00082 if ( ! (fDx1 == fDx2 && fDx3 == fDx4 ) ) { 00083 std::ostringstream message; 00084 message << "TwistedTrapBoxSide is not used as a the side of a box: " 00085 << GetName() << G4endl 00086 << " Not a box !"; 00087 G4Exception("G4TwistBoxSide::G4TwistBoxSide()", "GeomSolids0002", 00088 FatalException, message); 00089 } 00090 00091 fDy1 = pDy1 ; 00092 fDy2 = pDy2 ; 00093 00094 fDz = pDz ; 00095 00096 fAlph = pAlph ; 00097 fTAlph = std::tan(fAlph) ; 00098 00099 fTheta = pTheta ; 00100 fPhi = pPhi ; 00101 00102 // precalculate frequently used parameters 00103 fDx4plus2 = fDx4 + fDx2 ; 00104 fDx4minus2 = fDx4 - fDx2 ; 00105 fDx3plus1 = fDx3 + fDx1 ; 00106 fDx3minus1 = fDx3 - fDx1 ; 00107 fDy2plus1 = fDy2 + fDy1 ; 00108 fDy2minus1 = fDy2 - fDy1 ; 00109 00110 fa1md1 = 2*fDx2 - 2*fDx1 ; 00111 fa2md2 = 2*fDx4 - 2*fDx3 ; 00112 00113 00114 fPhiTwist = PhiTwist ; // dphi 00115 fAngleSide = AngleSide ; // 0,90,180,270 deg 00116 00117 fdeltaX = 2 * fDz * std::tan(fTheta) * std::cos(fPhi) ; // dx in surface equation 00118 fdeltaY = 2 * fDz * std::tan(fTheta) * std::sin(fPhi) ; // dy in surface equation 00119 00120 fRot.rotateZ( AngleSide ) ; 00121 00122 fTrans.set(0, 0, 0); // No Translation 00123 fIsValidNorm = false; 00124 00125 SetCorners() ; 00126 SetBoundaries() ; 00127 00128 }
G4TwistBoxSide::~G4TwistBoxSide | ( | ) | [virtual] |
G4TwistBoxSide::G4TwistBoxSide | ( | __void__ & | ) |
Definition at line 134 of file G4TwistBoxSide.cc.
00135 : G4VTwistSurface(a), fTheta(0.), fPhi(0.), fDy1(0.), fDx1(0.), fDx2(0.), 00136 fDy2(0.), fDx3(0.), fDx4(0.), fDz(0.), fAlph(0.), fTAlph(0.), fPhiTwist(0.), 00137 fAngleSide(0.), fdeltaX(0.), fdeltaY(0.), fDx4plus2(0.), fDx4minus2(0.), 00138 fDx3plus1(0.), fDx3minus1(0.), fDy2plus1(0.), fDy2minus1(0.), fa1md1(0.), 00139 fa2md2(0.) 00140 { 00141 }
G4int G4TwistBoxSide::DistanceToSurface | ( | const G4ThreeVector & | gp, | |
G4ThreeVector | gxx[], | |||
G4double | distance[], | |||
G4int | areacode[] | |||
) | [virtual] |
Implements G4VTwistSurface.
Definition at line 673 of file G4TwistBoxSide.cc.
References G4VTwistSurface::ComputeGlobalPoint(), G4VTwistSurface::ComputeLocalPoint(), G4VTwistSurface::DistanceToPlane(), G4VTwistSurface::fCurStat, G4cout, G4endl, G4VSURFACENXX, G4VTwistSurface::kCarTolerance, G4VTwistSurface::kDontValidate, and G4VTwistSurface::sOutside.
00677 { 00678 // to do 00679 00680 static const G4double ctol = 0.5 * kCarTolerance; 00681 00682 fCurStat.ResetfDone(kDontValidate, &gp); 00683 00684 if (fCurStat.IsDone()) { 00685 G4int i; 00686 for (i=0; i<fCurStat.GetNXX(); i++) { 00687 gxx[i] = fCurStat.GetXX(i); 00688 distance[i] = fCurStat.GetDistance(i); 00689 areacode[i] = fCurStat.GetAreacode(i); 00690 } 00691 return fCurStat.GetNXX(); 00692 } else { 00693 // initialize 00694 G4int i; 00695 for (i=0; i<G4VSURFACENXX; i++) { 00696 distance[i] = kInfinity; 00697 areacode[i] = sOutside; 00698 gxx[i].set(kInfinity, kInfinity, kInfinity); 00699 } 00700 } 00701 00702 G4ThreeVector p = ComputeLocalPoint(gp); 00703 G4ThreeVector xx; // intersection point 00704 G4ThreeVector xxonsurface ; // interpolated intersection point 00705 00706 // the surfacenormal at that surface point 00707 G4double phiR = 0 ; // 00708 G4double uR = 0 ; 00709 00710 G4ThreeVector surfacenormal ; 00711 G4double deltaX ; 00712 00713 G4int maxint = 20 ; 00714 00715 for ( G4int i = 1 ; i<maxint ; i++ ) { 00716 00717 xxonsurface = SurfacePoint(phiR,uR) ; 00718 surfacenormal = NormAng(phiR,uR) ; 00719 distance[0] = DistanceToPlane(p, xxonsurface, surfacenormal, xx); // new XX 00720 deltaX = ( xx - xxonsurface ).mag() ; 00721 00722 #ifdef G4TWISTDEBUG 00723 G4cout << "i = " << i << ", distance = " << distance[0] << ", " << deltaX << G4endl ; 00724 G4cout << "X = " << xx << G4endl ; 00725 #endif 00726 00727 // the new point xx is accepted and phi/psi replaced 00728 GetPhiUAtX(xx, phiR, uR) ; 00729 00730 if ( deltaX <= ctol ) { break ; } 00731 00732 } 00733 00734 // check validity of solution ( valid phi,psi ) 00735 00736 G4double halfphi = 0.5*fPhiTwist ; 00737 G4double uMax = GetBoundaryMax(phiR) ; 00738 00739 if ( phiR > halfphi ) phiR = halfphi ; 00740 if ( phiR < -halfphi ) phiR = -halfphi ; 00741 if ( uR > uMax ) uR = uMax ; 00742 if ( uR < -uMax ) uR = -uMax ; 00743 00744 xxonsurface = SurfacePoint(phiR,uR) ; 00745 distance[0] = ( p - xx ).mag() ; 00746 if ( distance[0] <= ctol ) { distance[0] = 0 ; } 00747 00748 // end of validity 00749 00750 #ifdef G4TWISTDEBUG 00751 G4cout << "refined solution " << phiR << " , " << uR << " , " << G4endl ; 00752 G4cout << "distance = " << distance[0] << G4endl ; 00753 G4cout << "X = " << xx << G4endl ; 00754 #endif 00755 00756 G4bool isvalid = true; 00757 gxx[0] = ComputeGlobalPoint(xx); 00758 00759 #ifdef G4TWISTDEBUG 00760 G4cout << "intersection Point found: " << gxx[0] << G4endl ; 00761 G4cout << "distance = " << distance[0] << G4endl ; 00762 #endif 00763 00764 fCurStat.SetCurrentStatus(0, gxx[0], distance[0], areacode[0], 00765 isvalid, 1, kDontValidate, &gp); 00766 return 1; 00767 00768 00769 }
G4int G4TwistBoxSide::DistanceToSurface | ( | const G4ThreeVector & | gp, | |
const G4ThreeVector & | gv, | |||
G4ThreeVector | gxx[], | |||
G4double | distance[], | |||
G4int | areacode[], | |||
G4bool | isvalid[], | |||
EValidate | validate = kValidateWithTol | |||
) | [virtual] |
Definition at line 200 of file G4TwistBoxSide.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.
00207 { 00208 00209 static const G4double ctol = 0.5 * kCarTolerance; 00210 static const G4double pihalf = pi/2 ; 00211 00212 G4bool IsParallel = false ; 00213 G4bool IsConverged = false ; 00214 00215 G4int nxx = 0 ; // number of physical solutions 00216 00217 fCurStatWithV.ResetfDone(validate, &gp, &gv); 00218 00219 if (fCurStatWithV.IsDone()) { 00220 G4int i; 00221 for (i=0; i<fCurStatWithV.GetNXX(); i++) { 00222 gxx[i] = fCurStatWithV.GetXX(i); 00223 distance[i] = fCurStatWithV.GetDistance(i); 00224 areacode[i] = fCurStatWithV.GetAreacode(i); 00225 isvalid[i] = fCurStatWithV.IsValid(i); 00226 } 00227 return fCurStatWithV.GetNXX(); 00228 } else { 00229 00230 // initialize 00231 G4int i; 00232 for (i=0; i<G4VSURFACENXX ; i++) { 00233 distance[i] = kInfinity; 00234 areacode[i] = sOutside; 00235 isvalid[i] = false; 00236 gxx[i].set(kInfinity, kInfinity, kInfinity); 00237 } 00238 } 00239 00240 G4ThreeVector p = ComputeLocalPoint(gp); 00241 G4ThreeVector v = ComputeLocalDirection(gv); 00242 00243 #ifdef G4TWISTDEBUG 00244 G4cout << "Local point p = " << p << G4endl ; 00245 G4cout << "Local direction v = " << v << G4endl ; 00246 #endif 00247 00248 G4double phi=0.,u=0.,q=0.; // parameters 00249 00250 // temporary variables 00251 00252 G4double tmpdist = kInfinity ; 00253 G4ThreeVector tmpxx; 00254 G4int tmpareacode = sOutside ; 00255 G4bool tmpisvalid = false ; 00256 00257 std::vector<Intersection> xbuf ; 00258 Intersection xbuftmp ; 00259 00260 // prepare some variables for the intersection finder 00261 00262 G4double L = 2*fDz ; 00263 00264 G4double phixz = fPhiTwist * ( p.x() * v.z() - p.z() * v.x() ) ; 00265 G4double phiyz = fPhiTwist * ( p.y() * v.z() - p.z() * v.y() ) ; 00266 00267 // special case vz = 0 00268 00269 if ( v.z() == 0. ) { 00270 00271 if ( (std::fabs(p.z()) <= L) && fPhiTwist ) { // intersection possible in z 00272 00273 phi = p.z() * fPhiTwist / L ; // phi is determined by the z-position 00274 00275 q = (2.* fPhiTwist*((v.x() - fTAlph*v.y())*std::cos(phi) 00276 + (fTAlph*v.x() + v.y())*std::sin(phi))); 00277 00278 if (q) { 00279 u = (2*(-(fdeltaY*phi*v.x()) + fPhiTwist*p.y()*v.x() 00280 + fdeltaX*phi*v.y() - fPhiTwist*p.x()*v.y()) 00281 + (fDx4plus2*fPhiTwist + 2*fDx4minus2*phi) 00282 * (v.y()*std::cos(phi) - v.x()*std::sin(phi))) / q; 00283 } 00284 00285 xbuftmp.phi = phi ; 00286 xbuftmp.u = u ; 00287 xbuftmp.areacode = sOutside ; 00288 xbuftmp.distance = kInfinity ; 00289 xbuftmp.isvalid = false ; 00290 00291 xbuf.push_back(xbuftmp) ; // store it to xbuf 00292 00293 } 00294 00295 else { // no intersection possible 00296 00297 distance[0] = kInfinity; 00298 gxx[0].set(kInfinity,kInfinity,kInfinity); 00299 isvalid[0] = false ; 00300 areacode[0] = sOutside ; 00301 fCurStatWithV.SetCurrentStatus(0, gxx[0], distance[0], 00302 areacode[0], isvalid[0], 00303 0, validate, &gp, &gv); 00304 00305 return 0; 00306 00307 00308 } // end std::fabs(p.z() <= L 00309 00310 } // end v.z() == 0 00311 00312 00313 // general solution for non-zero vz 00314 00315 else { 00316 00317 G4double c[8],srd[7],si[7] ; 00318 00319 c[7] = -14400*(-2*phixz + 2*fTAlph*phiyz + fDx4plus2*fPhiTwist*v.z()) ; 00320 c[6] = 28800*(phiyz + 2*fDz*v.x() - (fdeltaX + fDx4minus2)*v.z() + fTAlph*(phixz - 2*fDz*v.y() + fdeltaY*v.z())) ; 00321 c[5] = -1200*(10*phixz - 48*fDz*v.y() + 24*fdeltaY*v.z() + fDx4plus2*fPhiTwist*v.z() - 2*fTAlph*(5*phiyz + 24*fDz*v.x() - 12*fdeltaX*v.z())) ; 00322 c[4] = -2400*(phiyz + 10*fDz*v.x() - 5*fdeltaX*v.z() + fDx4minus2*v.z() + fTAlph*(phixz - 10*fDz*v.y() + 5*fdeltaY*v.z())) ; 00323 c[3] = 24*(2*phixz - 200*fDz*v.y() + 100*fdeltaY*v.z() - fDx4plus2*fPhiTwist*v.z() - 2*fTAlph*(phiyz + 100*fDz*v.x() - 50*fdeltaX*v.z())) ; 00324 c[2] = -16*(7*fTAlph* phixz + 7*phiyz - 6*fDz*v.x() + 6*fDz*fTAlph*v.y() + 3*(fdeltaX + fDx4minus2 - fdeltaY*fTAlph)*v.z()) ; 00325 c[1] = 4*(9*phixz - 9*fTAlph*phiyz - 56*fDz*fTAlph*v.x() - 56*fDz*v.y() + 28*(fdeltaY + fdeltaX*fTAlph)*v.z()) ; 00326 c[0] = 36*(2* fDz*(v.x() - fTAlph*v.y()) - fdeltaX*v.z() + fdeltaY*fTAlph*v.z()) ; 00327 00328 00329 #ifdef G4TWISTDEBUG 00330 G4cout << "coef = " << c[0] << " " 00331 << c[1] << " " 00332 << c[2] << " " 00333 << c[3] << " " 00334 << c[4] << " " 00335 << c[5] << " " 00336 << c[6] << " " 00337 << c[7] << G4endl ; 00338 #endif 00339 00340 G4JTPolynomialSolver trapEq ; 00341 G4int num = trapEq.FindRoots(c,7,srd,si); 00342 00343 00344 for (G4int i = 0 ; i<num ; i++ ) { // loop over all mathematical solutions 00345 if ( (si[i]==0.0) && fPhiTwist ) { // only real solutions 00346 #ifdef G4TWISTDEBUG 00347 G4cout << "Solution " << i << " : " << srd[i] << G4endl ; 00348 #endif 00349 phi = std::fmod(srd[i] , pihalf) ; 00350 00351 u = (2*phiyz + 4*fDz*phi*v.y() - 2*fdeltaY*phi*v.z() 00352 - fDx4plus2*fPhiTwist*v.z()*std::sin(phi) 00353 - 2*fDx4minus2*phi*v.z()*std::sin(phi)) 00354 /(2*fPhiTwist*v.z()*std::cos(phi) 00355 + 2*fPhiTwist*fTAlph*v.z()*std::sin(phi)) ; 00356 00357 xbuftmp.phi = phi ; 00358 xbuftmp.u = u ; 00359 xbuftmp.areacode = sOutside ; 00360 xbuftmp.distance = kInfinity ; 00361 xbuftmp.isvalid = false ; 00362 00363 xbuf.push_back(xbuftmp) ; // store it to xbuf 00364 00365 #ifdef G4TWISTDEBUG 00366 G4cout << "solution " << i << " = " << phi << " , " << u << G4endl ; 00367 #endif 00368 00369 } // end if real solution 00370 } // end loop i 00371 00372 } // end general case 00373 00374 00375 nxx = xbuf.size() ; // save the number of solutions 00376 00377 G4ThreeVector xxonsurface ; // point on surface 00378 G4ThreeVector surfacenormal ; // normal vector 00379 G4double deltaX ; // distance between intersection point and point on surface 00380 G4double theta ; // angle between track and surfacenormal 00381 G4double factor ; // a scaling factor 00382 G4int maxint = 30 ; // number of iterations 00383 00384 00385 for ( size_t k = 0 ; k<xbuf.size() ; k++ ) { 00386 00387 #ifdef G4TWISTDEBUG 00388 G4cout << "Solution " << k << " : " 00389 << "reconstructed phiR = " << xbuf[k].phi 00390 << ", uR = " << xbuf[k].u << G4endl ; 00391 #endif 00392 00393 phi = xbuf[k].phi ; // get the stored values for phi and u 00394 u = xbuf[k].u ; 00395 00396 IsConverged = false ; // no convergence at the beginning 00397 00398 for ( G4int i = 1 ; i<maxint ; i++ ) { 00399 00400 xxonsurface = SurfacePoint(phi,u) ; 00401 surfacenormal = NormAng(phi,u) ; 00402 tmpdist = DistanceToPlaneWithV(p, v, xxonsurface, surfacenormal, tmpxx); 00403 deltaX = ( tmpxx - xxonsurface ).mag() ; 00404 theta = std::fabs(std::acos(v*surfacenormal) - pihalf) ; 00405 if ( theta < 0.001 ) { 00406 factor = 50 ; 00407 IsParallel = true ; 00408 } 00409 else { 00410 factor = 1 ; 00411 } 00412 00413 #ifdef G4TWISTDEBUG 00414 G4cout << "Step i = " << i << ", distance = " << tmpdist << ", " << deltaX << G4endl ; 00415 G4cout << "X = " << tmpxx << G4endl ; 00416 #endif 00417 00418 GetPhiUAtX(tmpxx, phi, u) ; // the new point xx is accepted and phi/u replaced 00419 00420 #ifdef G4TWISTDEBUG 00421 G4cout << "approximated phi = " << phi << ", u = " << u << G4endl ; 00422 #endif 00423 00424 if ( deltaX <= factor*ctol ) { IsConverged = true ; break ; } 00425 00426 } // end iterative loop (i) 00427 00428 00429 // new code 21.09.05 O.Link 00430 if ( std::fabs(tmpdist)<ctol ) tmpdist = 0 ; 00431 00432 #ifdef G4TWISTDEBUG 00433 G4cout << "refined solution " << phi << " , " << u << G4endl ; 00434 G4cout << "distance = " << tmpdist << G4endl ; 00435 G4cout << "local X = " << tmpxx << G4endl ; 00436 #endif 00437 00438 tmpisvalid = false ; // init 00439 00440 if ( IsConverged ) { 00441 00442 if (validate == kValidateWithTol) { 00443 tmpareacode = GetAreaCode(tmpxx); 00444 if (!IsOutside(tmpareacode)) { 00445 if (tmpdist >= 0) tmpisvalid = true; 00446 } 00447 } else if (validate == kValidateWithoutTol) { 00448 tmpareacode = GetAreaCode(tmpxx, false); 00449 if (IsInside(tmpareacode)) { 00450 if (tmpdist >= 0) tmpisvalid = true; 00451 } 00452 } else { // kDontValidate 00453 G4Exception("G4TwistBoxSide::DistanceToSurface()", 00454 "GeomSolids0001", FatalException, 00455 "Feature NOT implemented !"); 00456 } 00457 00458 } 00459 else { 00460 tmpdist = kInfinity; // no convergence after 10 steps 00461 tmpisvalid = false ; // solution is not vaild 00462 } 00463 00464 00465 // store the found values 00466 xbuf[k].xx = tmpxx ; 00467 xbuf[k].distance = tmpdist ; 00468 xbuf[k].areacode = tmpareacode ; 00469 xbuf[k].isvalid = tmpisvalid ; 00470 00471 00472 } // end loop over physical solutions (variable k) 00473 00474 00475 std::sort(xbuf.begin() , xbuf.end(), DistanceSort ) ; // sorting 00476 00477 #ifdef G4TWISTDEBUG 00478 G4cout << G4endl << "list xbuf after sorting : " << G4endl ; 00479 G4cout << G4endl << G4endl ; 00480 #endif 00481 00482 00483 // erase identical intersection (within kCarTolerance) 00484 xbuf.erase( std::unique(xbuf.begin(), xbuf.end() , EqualIntersection ) , xbuf.end() ) ; 00485 00486 00487 // add guesses 00488 00489 G4int nxxtmp = xbuf.size() ; 00490 00491 if ( nxxtmp<2 || IsParallel ) { 00492 00493 // positive end 00494 #ifdef G4TWISTDEBUG 00495 G4cout << "add guess at +z/2 .. " << G4endl ; 00496 #endif 00497 00498 phi = fPhiTwist/2 ; 00499 u = 0 ; 00500 00501 00502 00503 xbuftmp.phi = phi ; 00504 xbuftmp.u = u ; 00505 xbuftmp.areacode = sOutside ; 00506 xbuftmp.distance = kInfinity ; 00507 xbuftmp.isvalid = false ; 00508 00509 xbuf.push_back(xbuftmp) ; // store it to xbuf 00510 00511 00512 #ifdef G4TWISTDEBUG 00513 G4cout << "add guess at -z/2 .. " << G4endl ; 00514 #endif 00515 00516 phi = -fPhiTwist/2 ; 00517 u = 0 ; 00518 00519 xbuftmp.phi = phi ; 00520 xbuftmp.u = u ; 00521 xbuftmp.areacode = sOutside ; 00522 xbuftmp.distance = kInfinity ; 00523 xbuftmp.isvalid = false ; 00524 00525 xbuf.push_back(xbuftmp) ; // store it to xbuf 00526 00527 for ( size_t k = nxxtmp ; k<xbuf.size() ; k++ ) { 00528 00529 #ifdef G4TWISTDEBUG 00530 G4cout << "Solution " << k << " : " 00531 << "reconstructed phiR = " << xbuf[k].phi 00532 << ", uR = " << xbuf[k].u << G4endl ; 00533 #endif 00534 00535 phi = xbuf[k].phi ; // get the stored values for phi and u 00536 u = xbuf[k].u ; 00537 00538 IsConverged = false ; // no convergence at the beginning 00539 00540 for ( G4int i = 1 ; i<maxint ; i++ ) { 00541 00542 xxonsurface = SurfacePoint(phi,u) ; 00543 surfacenormal = NormAng(phi,u) ; 00544 tmpdist = DistanceToPlaneWithV(p, v, xxonsurface, surfacenormal, tmpxx); 00545 deltaX = ( tmpxx - xxonsurface ).mag() ; 00546 theta = std::fabs(std::acos(v*surfacenormal) - pihalf) ; 00547 if ( theta < 0.001 ) { 00548 factor = 50 ; 00549 } 00550 else { 00551 factor = 1 ; 00552 } 00553 00554 #ifdef G4TWISTDEBUG 00555 G4cout << "Step i = " << i << ", distance = " << tmpdist << ", " << deltaX << G4endl ; 00556 G4cout << "X = " << tmpxx << G4endl ; 00557 #endif 00558 00559 GetPhiUAtX(tmpxx, phi, u) ; // the new point xx is accepted and phi/u replaced 00560 00561 #ifdef G4TWISTDEBUG 00562 G4cout << "approximated phi = " << phi << ", u = " << u << G4endl ; 00563 #endif 00564 00565 if ( deltaX <= factor*ctol ) { IsConverged = true ; break ; } 00566 00567 } // end iterative loop (i) 00568 00569 00570 // new code 21.09.05 O.Link 00571 if ( std::fabs(tmpdist)<ctol ) tmpdist = 0 ; 00572 00573 #ifdef G4TWISTDEBUG 00574 G4cout << "refined solution " << phi << " , " << u << G4endl ; 00575 G4cout << "distance = " << tmpdist << G4endl ; 00576 G4cout << "local X = " << tmpxx << G4endl ; 00577 #endif 00578 00579 tmpisvalid = false ; // init 00580 00581 if ( IsConverged ) { 00582 00583 if (validate == kValidateWithTol) { 00584 tmpareacode = GetAreaCode(tmpxx); 00585 if (!IsOutside(tmpareacode)) { 00586 if (tmpdist >= 0) tmpisvalid = true; 00587 } 00588 } else if (validate == kValidateWithoutTol) { 00589 tmpareacode = GetAreaCode(tmpxx, false); 00590 if (IsInside(tmpareacode)) { 00591 if (tmpdist >= 0) tmpisvalid = true; 00592 } 00593 } else { // kDontValidate 00594 G4Exception("G4TwistedBoxSide::DistanceToSurface()", 00595 "GeomSolids0001", FatalException, 00596 "Feature NOT implemented !"); 00597 } 00598 00599 } 00600 else { 00601 tmpdist = kInfinity; // no convergence after 10 steps 00602 tmpisvalid = false ; // solution is not vaild 00603 } 00604 00605 00606 // store the found values 00607 xbuf[k].xx = tmpxx ; 00608 xbuf[k].distance = tmpdist ; 00609 xbuf[k].areacode = tmpareacode ; 00610 xbuf[k].isvalid = tmpisvalid ; 00611 00612 00613 } // end loop over physical solutions 00614 00615 00616 } // end less than 2 solutions 00617 00618 00619 // sort again 00620 std::sort(xbuf.begin() , xbuf.end(), DistanceSort ) ; // sorting 00621 00622 // erase identical intersection (within kCarTolerance) 00623 xbuf.erase( std::unique(xbuf.begin(), xbuf.end() , EqualIntersection ) , xbuf.end() ) ; 00624 00625 #ifdef G4TWISTDEBUG 00626 G4cout << G4endl << "list xbuf after sorting : " << G4endl ; 00627 G4cout << G4endl << G4endl ; 00628 #endif 00629 00630 nxx = xbuf.size() ; // determine number of solutions again. 00631 00632 for ( size_t i = 0 ; i<xbuf.size() ; i++ ) { 00633 00634 distance[i] = xbuf[i].distance; 00635 gxx[i] = ComputeGlobalPoint(xbuf[i].xx); 00636 areacode[i] = xbuf[i].areacode ; 00637 isvalid[i] = xbuf[i].isvalid ; 00638 00639 fCurStatWithV.SetCurrentStatus(i, gxx[i], distance[i], areacode[i], 00640 isvalid[i], nxx, validate, &gp, &gv); 00641 00642 #ifdef G4TWISTDEBUG 00643 G4cout << "element Nr. " << i 00644 << ", local Intersection = " << xbuf[i].xx 00645 << ", distance = " << xbuf[i].distance 00646 << ", u = " << xbuf[i].u 00647 << ", phi = " << xbuf[i].phi 00648 << ", isvalid = " << xbuf[i].isvalid 00649 << G4endl ; 00650 #endif 00651 00652 } // end for( i ) loop 00653 00654 00655 #ifdef G4TWISTDEBUG 00656 G4cout << "G4TwistBoxSide finished " << G4endl ; 00657 G4cout << nxx << " possible physical solutions found" << G4endl ; 00658 for ( G4int k= 0 ; k< nxx ; k++ ) { 00659 G4cout << "global intersection Point found: " << gxx[k] << G4endl ; 00660 G4cout << "distance = " << distance[k] << G4endl ; 00661 G4cout << "isvalid = " << isvalid[k] << G4endl ; 00662 } 00663 #endif 00664 00665 return nxx ; 00666 00667 }
G4ThreeVector G4TwistBoxSide::GetNormal | ( | const G4ThreeVector & | xx, | |
G4bool | isGlobal = false | |||
) | [virtual] |
Implements G4VTwistSurface.
Definition at line 154 of file G4TwistBoxSide.cc.
References G4VTwistSurface::ComputeGlobalDirection(), G4VTwistSurface::ComputeLocalPoint(), G4VTwistSurface::fCurrentNormal, G4cout, G4endl, and G4VTwistSurface::kCarTolerance.
00156 { 00157 // GetNormal returns a normal vector at a surface (or very close 00158 // to surface) point at tmpxx. 00159 // If isGlobal=true, it returns the normal in global coordinate. 00160 // 00161 00162 G4ThreeVector xx; 00163 if (isGlobal) { 00164 xx = ComputeLocalPoint(tmpxx); 00165 if ((xx - fCurrentNormal.p).mag() < 0.5 * kCarTolerance) { 00166 return ComputeGlobalDirection(fCurrentNormal.normal); 00167 } 00168 } else { 00169 xx = tmpxx; 00170 if (xx == fCurrentNormal.p) { 00171 return fCurrentNormal.normal; 00172 } 00173 } 00174 00175 G4double phi ; 00176 G4double u ; 00177 00178 GetPhiUAtX(xx,phi,u) ; // phi,u for point xx close to surface 00179 00180 G4ThreeVector normal = NormAng(phi,u) ; // the normal vector at phi,u 00181 00182 #ifdef G4TWISTDEBUG 00183 G4cout << "normal vector = " << normal << G4endl ; 00184 G4cout << "phi = " << phi << " , u = " << u << G4endl ; 00185 #endif 00186 00187 // normal = normal/normal.mag() ; 00188 00189 if (isGlobal) { 00190 fCurrentNormal.normal = ComputeGlobalDirection(normal.unit()); 00191 } else { 00192 fCurrentNormal.normal = normal.unit(); 00193 } 00194 return fCurrentNormal.normal; 00195 }