00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042 #include <cmath>
00043
00044 #include "G4TwistBoxSide.hh"
00045 #include "G4PhysicalConstants.hh"
00046 #include "G4JTPolynomialSolver.hh"
00047
00048
00049
00050
00051 G4TwistBoxSide::G4TwistBoxSide(const G4String &name,
00052 G4double PhiTwist,
00053 G4double pDz,
00054 G4double pTheta,
00055 G4double pPhi,
00056 G4double pDy1,
00057 G4double pDx1,
00058 G4double pDx2,
00059 G4double pDy2,
00060 G4double pDx3,
00061 G4double pDx4,
00062 G4double pAlph,
00063 G4double AngleSide
00064 ) : G4VTwistSurface(name)
00065 {
00066
00067
00068 fAxis[0] = kYAxis;
00069 fAxis[1] = kZAxis;
00070 fAxisMin[0] = -kInfinity ;
00071 fAxisMax[0] = kInfinity ;
00072 fAxisMin[1] = -pDz ;
00073 fAxisMax[1] = pDz ;
00074
00075 fDx1 = pDx1 ;
00076 fDx2 = pDx2 ;
00077 fDx3 = pDx3 ;
00078 fDx4 = pDx4 ;
00079
00080
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
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 ;
00115 fAngleSide = AngleSide ;
00116
00117 fdeltaX = 2 * fDz * std::tan(fTheta) * std::cos(fPhi) ;
00118 fdeltaY = 2 * fDz * std::tan(fTheta) * std::sin(fPhi) ;
00119
00120 fRot.rotateZ( AngleSide ) ;
00121
00122 fTrans.set(0, 0, 0);
00123 fIsValidNorm = false;
00124
00125 SetCorners() ;
00126 SetBoundaries() ;
00127
00128 }
00129
00130
00131
00132
00133
00134 G4TwistBoxSide::G4TwistBoxSide( __void__& a )
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 }
00142
00143
00144
00145
00146
00147 G4TwistBoxSide::~G4TwistBoxSide()
00148 {
00149 }
00150
00151
00152
00153
00154 G4ThreeVector G4TwistBoxSide::GetNormal(const G4ThreeVector &tmpxx,
00155 G4bool isGlobal)
00156 {
00157
00158
00159
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) ;
00179
00180 G4ThreeVector normal = NormAng(phi,u) ;
00181
00182 #ifdef G4TWISTDEBUG
00183 G4cout << "normal vector = " << normal << G4endl ;
00184 G4cout << "phi = " << phi << " , u = " << u << G4endl ;
00185 #endif
00186
00187
00188
00189 if (isGlobal) {
00190 fCurrentNormal.normal = ComputeGlobalDirection(normal.unit());
00191 } else {
00192 fCurrentNormal.normal = normal.unit();
00193 }
00194 return fCurrentNormal.normal;
00195 }
00196
00197
00198
00199
00200 G4int G4TwistBoxSide::DistanceToSurface(const G4ThreeVector &gp,
00201 const G4ThreeVector &gv,
00202 G4ThreeVector gxx[],
00203 G4double distance[],
00204 G4int areacode[],
00205 G4bool isvalid[],
00206 EValidate validate)
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 ;
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
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.;
00249
00250
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
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
00268
00269 if ( v.z() == 0. ) {
00270
00271 if ( (std::fabs(p.z()) <= L) && fPhiTwist ) {
00272
00273 phi = p.z() * fPhiTwist / L ;
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) ;
00292
00293 }
00294
00295 else {
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 }
00309
00310 }
00311
00312
00313
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++ ) {
00345 if ( (si[i]==0.0) && fPhiTwist ) {
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) ;
00364
00365 #ifdef G4TWISTDEBUG
00366 G4cout << "solution " << i << " = " << phi << " , " << u << G4endl ;
00367 #endif
00368
00369 }
00370 }
00371
00372 }
00373
00374
00375 nxx = xbuf.size() ;
00376
00377 G4ThreeVector xxonsurface ;
00378 G4ThreeVector surfacenormal ;
00379 G4double deltaX ;
00380 G4double theta ;
00381 G4double factor ;
00382 G4int maxint = 30 ;
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 ;
00394 u = xbuf[k].u ;
00395
00396 IsConverged = false ;
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) ;
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 }
00427
00428
00429
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 ;
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 {
00453 G4Exception("G4TwistBoxSide::DistanceToSurface()",
00454 "GeomSolids0001", FatalException,
00455 "Feature NOT implemented !");
00456 }
00457
00458 }
00459 else {
00460 tmpdist = kInfinity;
00461 tmpisvalid = false ;
00462 }
00463
00464
00465
00466 xbuf[k].xx = tmpxx ;
00467 xbuf[k].distance = tmpdist ;
00468 xbuf[k].areacode = tmpareacode ;
00469 xbuf[k].isvalid = tmpisvalid ;
00470
00471
00472 }
00473
00474
00475 std::sort(xbuf.begin() , xbuf.end(), DistanceSort ) ;
00476
00477 #ifdef G4TWISTDEBUG
00478 G4cout << G4endl << "list xbuf after sorting : " << G4endl ;
00479 G4cout << G4endl << G4endl ;
00480 #endif
00481
00482
00483
00484 xbuf.erase( std::unique(xbuf.begin(), xbuf.end() , EqualIntersection ) , xbuf.end() ) ;
00485
00486
00487
00488
00489 G4int nxxtmp = xbuf.size() ;
00490
00491 if ( nxxtmp<2 || IsParallel ) {
00492
00493
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) ;
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) ;
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 ;
00536 u = xbuf[k].u ;
00537
00538 IsConverged = false ;
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) ;
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 }
00568
00569
00570
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 ;
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 {
00594 G4Exception("G4TwistedBoxSide::DistanceToSurface()",
00595 "GeomSolids0001", FatalException,
00596 "Feature NOT implemented !");
00597 }
00598
00599 }
00600 else {
00601 tmpdist = kInfinity;
00602 tmpisvalid = false ;
00603 }
00604
00605
00606
00607 xbuf[k].xx = tmpxx ;
00608 xbuf[k].distance = tmpdist ;
00609 xbuf[k].areacode = tmpareacode ;
00610 xbuf[k].isvalid = tmpisvalid ;
00611
00612
00613 }
00614
00615
00616 }
00617
00618
00619
00620 std::sort(xbuf.begin() , xbuf.end(), DistanceSort ) ;
00621
00622
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() ;
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 }
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 }
00668
00669
00670
00671
00672
00673 G4int G4TwistBoxSide::DistanceToSurface(const G4ThreeVector &gp,
00674 G4ThreeVector gxx[],
00675 G4double distance[],
00676 G4int areacode[])
00677 {
00678
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
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;
00704 G4ThreeVector xxonsurface ;
00705
00706
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);
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
00728 GetPhiUAtX(xx, phiR, uR) ;
00729
00730 if ( deltaX <= ctol ) { break ; }
00731
00732 }
00733
00734
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
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 }
00770
00771
00772
00773
00774
00775 G4int G4TwistBoxSide::GetAreaCode(const G4ThreeVector &xx,
00776 G4bool withTol)
00777 {
00778
00779
00780
00781 static const G4double ctol = 0.5 * kCarTolerance;
00782
00783 G4double phi ;
00784 G4double yprime ;
00785 GetPhiUAtX(xx, phi,yprime ) ;
00786
00787 G4double fYAxisMax = GetBoundaryMax(phi) ;
00788 G4double fYAxisMin = - fYAxisMax ;
00789
00790 #ifdef G4TWISTDEBUG
00791 G4cout << "GetAreaCode: phi = " << phi << G4endl ;
00792 G4cout << "GetAreaCode: yprime = " << yprime << G4endl ;
00793 G4cout << "Intervall is " << fYAxisMin << " to " << fYAxisMax << G4endl ;
00794 #endif
00795
00796 G4int areacode = sInside;
00797
00798 if (fAxis[0] == kYAxis && fAxis[1] == kZAxis) {
00799
00800 G4int zaxis = 1;
00801
00802 if (withTol) {
00803
00804 G4bool isoutside = false;
00805
00806
00807
00808 if (yprime < fYAxisMin + ctol) {
00809 areacode |= (sAxis0 & (sAxisY | sAxisMin)) | sBoundary;
00810 if (yprime <= fYAxisMin - ctol) isoutside = true;
00811
00812 } else if (yprime > fYAxisMax - ctol) {
00813 areacode |= (sAxis0 & (sAxisY | sAxisMax)) | sBoundary;
00814 if (yprime >= fYAxisMax + ctol) isoutside = true;
00815 }
00816
00817
00818
00819 if (xx.z() < fAxisMin[zaxis] + ctol) {
00820 areacode |= (sAxis1 & (sAxisZ | sAxisMin));
00821
00822 if (areacode & sBoundary) areacode |= sCorner;
00823 else areacode |= sBoundary;
00824 if (xx.z() <= fAxisMin[zaxis] - ctol) isoutside = true;
00825
00826 } else if (xx.z() > fAxisMax[zaxis] - ctol) {
00827 areacode |= (sAxis1 & (sAxisZ | sAxisMax));
00828
00829 if (areacode & sBoundary) areacode |= sCorner;
00830 else areacode |= sBoundary;
00831 if (xx.z() >= fAxisMax[zaxis] + ctol) isoutside = true;
00832 }
00833
00834
00835
00836
00837 if (isoutside) {
00838 G4int tmpareacode = areacode & (~sInside);
00839 areacode = tmpareacode;
00840 } else if ((areacode & sBoundary) != sBoundary) {
00841 areacode |= (sAxis0 & sAxisY) | (sAxis1 & sAxisZ);
00842 }
00843
00844 } else {
00845
00846
00847
00848 if (yprime < fYAxisMin ) {
00849 areacode |= (sAxis0 & (sAxisY | sAxisMin)) | sBoundary;
00850 } else if (yprime > fYAxisMax) {
00851 areacode |= (sAxis0 & (sAxisY | sAxisMax)) | sBoundary;
00852 }
00853
00854
00855
00856 if (xx.z() < fAxisMin[zaxis]) {
00857 areacode |= (sAxis1 & (sAxisZ | sAxisMin));
00858 if (areacode & sBoundary) areacode |= sCorner;
00859 else areacode |= sBoundary;
00860
00861 } else if (xx.z() > fAxisMax[zaxis]) {
00862 areacode |= (sAxis1 & (sAxisZ | sAxisMax)) ;
00863 if (areacode & sBoundary) areacode |= sCorner;
00864 else areacode |= sBoundary;
00865 }
00866
00867 if ((areacode & sBoundary) != sBoundary) {
00868 areacode |= (sAxis0 & sAxisY) | (sAxis1 & sAxisZ);
00869 }
00870 }
00871 return areacode;
00872 } else {
00873 G4Exception("G4TwistBoxSide::GetAreaCode()",
00874 "GeomSolids0001", FatalException,
00875 "Feature NOT implemented !");
00876 }
00877 return areacode;
00878 }
00879
00880
00881
00882
00883 void G4TwistBoxSide::SetCorners()
00884 {
00885
00886
00887
00888 if (fAxis[0] == kYAxis && fAxis[1] == kZAxis) {
00889
00890 G4double x, y, z;
00891
00892
00893
00894 x = -fdeltaX/2. + (fDx2 - fDy1*fTAlph)*std::cos(fPhiTwist/2.) - fDy1*std::sin(fPhiTwist/2.) ;
00895 y = -fdeltaY/2. - fDy1*std::cos(fPhiTwist/2.) + (-fDx2 + fDy1*fTAlph)*std::sin(fPhiTwist/2.) ;
00896 z = -fDz ;
00897
00898 SetCorner(sC0Min1Min, x, y, z);
00899
00900
00901
00902 x = -fdeltaX/2. + (fDx2 + fDy1*fTAlph)*std::cos(fPhiTwist/2.) + fDy1*std::sin(fPhiTwist/2.) ;
00903 y = -fdeltaY/2. + fDy1*std::cos(fPhiTwist/2.) - (fDx2 + fDy1*fTAlph)*std::sin(fPhiTwist/2.) ;
00904 z = -fDz ;
00905
00906 SetCorner(sC0Max1Min, x, y, z);
00907
00908
00909
00910 x = fdeltaX/2. + (fDx4 + fDy2*fTAlph)*std::cos(fPhiTwist/2.) - fDy2*std::sin(fPhiTwist/2.) ;
00911 y = fdeltaY/2. + fDy2*std::cos(fPhiTwist/2.) + (fDx4 + fDy2*fTAlph)*std::sin(fPhiTwist/2.) ;
00912 z = fDz ;
00913
00914 SetCorner(sC0Max1Max, x, y, z);
00915
00916
00917
00918 x = fdeltaX/2. + (fDx4 - fDy2*fTAlph)*std::cos(fPhiTwist/2.) + fDy2*std::sin(fPhiTwist/2.) ;
00919 y = fdeltaY/2. - fDy2*std::cos(fPhiTwist/2.) + (fDx4 - fDy2*fTAlph)*std::sin(fPhiTwist/2.) ;
00920 z = fDz ;
00921
00922 SetCorner(sC0Min1Max, x, y, z);
00923
00924 } else {
00925
00926 G4Exception("G4TwistBoxSide::SetCorners()",
00927 "GeomSolids0001", FatalException,
00928 "Method NOT implemented !");
00929 }
00930 }
00931
00932
00933
00934
00935 void G4TwistBoxSide::SetBoundaries()
00936 {
00937
00938
00939
00940 G4ThreeVector direction;
00941
00942 if (fAxis[0] == kYAxis && fAxis[1] == kZAxis) {
00943
00944
00945 direction = GetCorner(sC0Min1Max) - GetCorner(sC0Min1Min);
00946 direction = direction.unit();
00947 SetBoundary(sAxis0 & (sAxisY | sAxisMin), direction,
00948 GetCorner(sC0Min1Min), sAxisZ) ;
00949
00950
00951 direction = GetCorner(sC0Max1Max) - GetCorner(sC0Max1Min);
00952 direction = direction.unit();
00953 SetBoundary(sAxis0 & (sAxisY | sAxisMax), direction,
00954 GetCorner(sC0Max1Min), sAxisZ);
00955
00956
00957 direction = GetCorner(sC0Max1Min) - GetCorner(sC0Min1Min);
00958 direction = direction.unit();
00959 SetBoundary(sAxis1 & (sAxisZ | sAxisMin), direction,
00960 GetCorner(sC0Min1Min), sAxisY);
00961
00962
00963 direction = GetCorner(sC0Max1Max) - GetCorner(sC0Min1Max);
00964 direction = direction.unit();
00965 SetBoundary(sAxis1 & (sAxisZ | sAxisMax), direction,
00966 GetCorner(sC0Min1Max), sAxisY);
00967
00968 } else {
00969
00970 G4Exception("G4TwistBoxSide::SetCorners()",
00971 "GeomSolids0001", FatalException,
00972 "Feature NOT implemented !");
00973 }
00974 }
00975
00976
00977 void G4TwistBoxSide::GetPhiUAtX( G4ThreeVector p, G4double &phi, G4double &u)
00978 {
00979
00980
00981
00982
00983
00984 phi = p.z()/(2*fDz)*fPhiTwist ;
00985
00986 u = -(fTAlph*(fDx4plus2*fPhiTwist + 2*fDx4minus2*phi) + 2*(fdeltaY*phi + fdeltaX*fTAlph*phi - fPhiTwist*(fTAlph*p.x() + p.y()))* std::cos(phi) + 2*(-(fdeltaX*phi) + fdeltaY*fTAlph*phi + fPhiTwist*(p.x() - fTAlph*p.y()))*std::sin(phi))/(2.*(fPhiTwist + fPhiTwist*fTAlph*fTAlph)) ;
00987
00988
00989 }
00990
00991
00992 G4ThreeVector G4TwistBoxSide::ProjectPoint(const G4ThreeVector &p,
00993 G4bool isglobal)
00994 {
00995
00996 G4ThreeVector tmpp;
00997 if (isglobal) {
00998 tmpp = fRot.inverse()*p - fTrans;
00999 } else {
01000 tmpp = p;
01001 }
01002
01003 G4double phi ;
01004 G4double u ;
01005
01006 GetPhiUAtX( tmpp, phi, u ) ;
01007
01008 G4ThreeVector xx = SurfacePoint(phi,u) ;
01009
01010 if (isglobal) {
01011 return (fRot * xx + fTrans);
01012 } else {
01013 return xx;
01014 }
01015 }
01016
01017 void G4TwistBoxSide::GetFacets( G4int k, G4int n, G4double xyz[][3], G4int faces[][4], G4int iside )
01018 {
01019
01020 G4double phi ;
01021 G4double b ;
01022
01023 G4double z, u ;
01024 G4ThreeVector p ;
01025
01026 G4int nnode ;
01027 G4int nface ;
01028
01029
01030
01031 G4int i,j ;
01032
01033 for ( i = 0 ; i<n ; i++ ) {
01034
01035 z = -fDz+i*(2.*fDz)/(n-1) ;
01036 phi = z*fPhiTwist/(2*fDz) ;
01037 b = GetValueB(phi) ;
01038
01039 for ( j = 0 ; j<k ; j++ ) {
01040
01041 nnode = GetNode(i,j,k,n,iside) ;
01042 u = -b/2 +j*b/(k-1) ;
01043 p = SurfacePoint(phi,u,true) ;
01044
01045 xyz[nnode][0] = p.x() ;
01046 xyz[nnode][1] = p.y() ;
01047 xyz[nnode][2] = p.z() ;
01048
01049 if ( i<n-1 && j<k-1 ) {
01050
01051 nface = GetFace(i,j,k,n,iside) ;
01052 faces[nface][0] = GetEdgeVisibility(i,j,k,n,0,-1) * (GetNode(i ,j ,k,n,iside)+1) ;
01053 faces[nface][1] = GetEdgeVisibility(i,j,k,n,1,-1) * (GetNode(i ,j+1,k,n,iside)+1) ;
01054 faces[nface][2] = GetEdgeVisibility(i,j,k,n,2,-1) * (GetNode(i+1,j+1,k,n,iside)+1) ;
01055 faces[nface][3] = GetEdgeVisibility(i,j,k,n,3,-1) * (GetNode(i+1,j ,k,n,iside)+1) ;
01056
01057 }
01058 }
01059 }
01060 }