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 "G4TwistTrapParallelSide.hh"
00045 #include "G4PhysicalConstants.hh"
00046 #include "G4JTPolynomialSolver.hh"
00047
00048
00049
00050
00051 G4TwistTrapParallelSide::G4TwistTrapParallelSide(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 )
00065 : G4VTwistSurface(name)
00066 {
00067
00068 fAxis[0] = kXAxis;
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 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
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 ;
00103 fAngleSide = AngleSide ;
00104
00105 fdeltaX = 2 * fDz * std::tan(fTheta) * std::cos(fPhi) ;
00106 fdeltaY = 2 * fDz * std::tan(fTheta) * std::sin(fPhi) ;
00107
00108 fRot.rotateZ( AngleSide ) ;
00109
00110 fTrans.set(0, 0, 0);
00111 fIsValidNorm = false;
00112
00113 SetCorners() ;
00114 SetBoundaries() ;
00115
00116 }
00117
00118
00119
00120
00121
00122 G4TwistTrapParallelSide::G4TwistTrapParallelSide( __void__& a )
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 }
00130
00131
00132
00133
00134
00135 G4TwistTrapParallelSide::~G4TwistTrapParallelSide()
00136 {
00137 }
00138
00139
00140
00141
00142 G4ThreeVector G4TwistTrapParallelSide::GetNormal(const G4ThreeVector &tmpxx,
00143 G4bool isGlobal)
00144 {
00145
00146
00147
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) ;
00167
00168 G4ThreeVector normal = NormAng(phi,u) ;
00169
00170 #ifdef G4TWISTDEBUG
00171 G4cout << "normal vector = " << normal << G4endl ;
00172 G4cout << "phi = " << phi << " , u = " << u << G4endl ;
00173 #endif
00174
00175
00176
00177 if (isGlobal) {
00178 fCurrentNormal.normal = ComputeGlobalDirection(normal.unit());
00179 } else {
00180 fCurrentNormal.normal = normal.unit();
00181 }
00182 return fCurrentNormal.normal;
00183 }
00184
00185
00186
00187
00188 G4int G4TwistTrapParallelSide::DistanceToSurface(const G4ThreeVector &gp,
00189 const G4ThreeVector &gv,
00190 G4ThreeVector gxx[],
00191 G4double distance[],
00192 G4int areacode[],
00193 G4bool isvalid[],
00194 EValidate validate)
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 ;
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
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 ;
00237
00238
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
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
00256
00257 if ( v.z() == 0. ) {
00258
00259 if ( std::fabs(p.z()) <= L ) {
00260
00261 phi = p.z() * fPhiTwist / L ;
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) ;
00275
00276 }
00277
00278 else {
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 }
00292
00293 }
00294
00295
00296
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++ ) {
00330 if ( si[i]==0.0 ) {
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) ;
00345
00346 #ifdef G4TWISTDEBUG
00347 G4cout << "solution " << i << " = " << phi << " , " << u << G4endl ;
00348 #endif
00349
00350 }
00351 }
00352
00353 }
00354
00355
00356 nxx = xbuf.size() ;
00357
00358 G4ThreeVector xxonsurface ;
00359 G4ThreeVector surfacenormal ;
00360 G4double deltaX ;
00361 G4double theta ;
00362 G4double factor ;
00363 G4int maxint = 30 ;
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 ;
00375 u = xbuf[k].u ;
00376
00377 IsConverged = false ;
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) ;
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 }
00408
00409
00410
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 ;
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 {
00435 G4Exception("G4TwistTrapParallelSide::DistanceToSurface()",
00436 "GeomSolids0001", FatalException,
00437 "Feature NOT implemented !");
00438 }
00439
00440 }
00441 else {
00442 tmpdist = kInfinity;
00443 tmpisvalid = false ;
00444 }
00445
00446
00447
00448 xbuf[k].xx = tmpxx ;
00449 xbuf[k].distance = tmpdist ;
00450 xbuf[k].areacode = tmpareacode ;
00451 xbuf[k].isvalid = tmpisvalid ;
00452
00453
00454 }
00455
00456
00457 std::sort(xbuf.begin() , xbuf.end(), DistanceSort ) ;
00458
00459 #ifdef G4TWISTDEBUG
00460 G4cout << G4endl << "list xbuf after sorting : " << G4endl ;
00461 G4cout << G4endl << G4endl ;
00462 #endif
00463
00464
00465
00466 xbuf.erase( std::unique(xbuf.begin(), xbuf.end() , EqualIntersection ) , xbuf.end() ) ;
00467
00468
00469
00470
00471 G4int nxxtmp = xbuf.size() ;
00472
00473 if ( nxxtmp<2 || IsParallel ) {
00474
00475
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) ;
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) ;
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 ;
00516 u = xbuf[k].u ;
00517
00518 IsConverged = false ;
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) ;
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 }
00548
00549
00550
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 ;
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 {
00575 G4Exception("G4TwistedBoxSide::DistanceToSurface()",
00576 "GeomSolids0001", FatalException,
00577 "Feature NOT implemented !");
00578 }
00579
00580 }
00581 else {
00582 tmpdist = kInfinity;
00583 tmpisvalid = false ;
00584 }
00585
00586
00587
00588 xbuf[k].xx = tmpxx ;
00589 xbuf[k].distance = tmpdist ;
00590 xbuf[k].areacode = tmpareacode ;
00591 xbuf[k].isvalid = tmpisvalid ;
00592
00593
00594 }
00595
00596
00597 }
00598
00599
00600
00601 std::sort(xbuf.begin() , xbuf.end(), DistanceSort ) ;
00602
00603
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() ;
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 }
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 }
00649
00650
00651
00652
00653
00654
00655 G4int G4TwistTrapParallelSide::DistanceToSurface(const G4ThreeVector &gp,
00656 G4ThreeVector gxx[],
00657 G4double distance[],
00658 G4int areacode[])
00659 {
00660
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
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;
00686 G4ThreeVector xxonsurface ;
00687
00688
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);
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
00710 GetPhiUAtX(xx, phiR, uR) ;
00711
00712 if ( deltaX <= ctol ) { break ; }
00713
00714 }
00715
00716
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
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 }
00751
00752
00753
00754
00755
00756 G4int G4TwistTrapParallelSide::GetAreaCode(const G4ThreeVector &xx,
00757 G4bool withTol)
00758 {
00759
00760
00761
00762 static const G4double ctol = 0.5 * kCarTolerance;
00763
00764 G4double phi ;
00765 G4double yprime ;
00766 GetPhiUAtX(xx, phi,yprime ) ;
00767
00768 G4double fXAxisMax = GetBoundaryMax(phi) ;
00769 G4double fXAxisMin = GetBoundaryMin(phi) ;
00770
00771 #ifdef G4TWISTDEBUG
00772 G4cout << "GetAreaCode: phi = " << phi << G4endl ;
00773 G4cout << "GetAreaCode: yprime = " << yprime << G4endl ;
00774 G4cout << "Intervall is " << fXAxisMin << " to " << fXAxisMax << G4endl ;
00775 #endif
00776
00777 G4int areacode = sInside;
00778
00779 if (fAxis[0] == kXAxis && fAxis[1] == kZAxis) {
00780
00781 G4int zaxis = 1;
00782
00783 if (withTol) {
00784
00785 G4bool isoutside = false;
00786
00787
00788
00789 if (yprime < fXAxisMin + ctol) {
00790 areacode |= (sAxis0 & (sAxisX | sAxisMin)) | sBoundary;
00791 if (yprime <= fXAxisMin - ctol) isoutside = true;
00792
00793 } else if (yprime > fXAxisMax - ctol) {
00794 areacode |= (sAxis0 & (sAxisX | sAxisMax)) | sBoundary;
00795 if (yprime >= fXAxisMax + ctol) isoutside = true;
00796 }
00797
00798
00799
00800 if (xx.z() < fAxisMin[zaxis] + ctol) {
00801 areacode |= (sAxis1 & (sAxisZ | sAxisMin));
00802
00803 if (areacode & sBoundary) areacode |= sCorner;
00804 else areacode |= sBoundary;
00805 if (xx.z() <= fAxisMin[zaxis] - ctol) isoutside = true;
00806
00807 } else if (xx.z() > fAxisMax[zaxis] - ctol) {
00808 areacode |= (sAxis1 & (sAxisZ | sAxisMax));
00809
00810 if (areacode & sBoundary) areacode |= sCorner;
00811 else areacode |= sBoundary;
00812 if (xx.z() >= fAxisMax[zaxis] + ctol) isoutside = true;
00813 }
00814
00815
00816
00817
00818 if (isoutside) {
00819 G4int tmpareacode = areacode & (~sInside);
00820 areacode = tmpareacode;
00821 } else if ((areacode & sBoundary) != sBoundary) {
00822 areacode |= (sAxis0 & sAxisX) | (sAxis1 & sAxisZ);
00823 }
00824
00825 } else {
00826
00827
00828
00829 if (yprime < fXAxisMin ) {
00830 areacode |= (sAxis0 & (sAxisX | sAxisMin)) | sBoundary;
00831 } else if (yprime > fXAxisMax) {
00832 areacode |= (sAxis0 & (sAxisX | sAxisMax)) | sBoundary;
00833 }
00834
00835
00836
00837 if (xx.z() < fAxisMin[zaxis]) {
00838 areacode |= (sAxis1 & (sAxisZ | sAxisMin));
00839 if (areacode & sBoundary) areacode |= sCorner;
00840 else areacode |= sBoundary;
00841
00842 } else if (xx.z() > fAxisMax[zaxis]) {
00843 areacode |= (sAxis1 & (sAxisZ | sAxisMax)) ;
00844 if (areacode & sBoundary) areacode |= sCorner;
00845 else areacode |= sBoundary;
00846 }
00847
00848 if ((areacode & sBoundary) != sBoundary) {
00849 areacode |= (sAxis0 & sAxisX) | (sAxis1 & sAxisZ);
00850 }
00851 }
00852 return areacode;
00853 } else {
00854 G4Exception("G4TwistTrapParallelSide::GetAreaCode()",
00855 "GeomSolids0001", FatalException,
00856 "Feature NOT implemented !");
00857 }
00858 return areacode;
00859 }
00860
00861
00862
00863
00864 void G4TwistTrapParallelSide::SetCorners()
00865 {
00866
00867
00868
00869 if (fAxis[0] == kXAxis && fAxis[1] == kZAxis) {
00870
00871 G4double x, y, z;
00872
00873
00874
00875 x = -fdeltaX/2. + (-fDx2 + fDy1*fTAlph)*std::cos(fPhiTwist/2.) + fDy1*std::sin(fPhiTwist/2.) ;
00876 y = -fdeltaY/2. + fDy1*std::cos(fPhiTwist/2.) + (fDx2 - fDy1*fTAlph)*std::sin(fPhiTwist/2.) ;
00877 z = -fDz ;
00878
00879 SetCorner(sC0Min1Min, x, y, z);
00880
00881
00882
00883 x = -fdeltaX/2. + (fDx2 + fDy1*fTAlph)*std::cos(fPhiTwist/2.) + fDy1*std::sin(fPhiTwist/2.) ;
00884 y = -fdeltaY/2. + fDy1*std::cos(fPhiTwist/2.) - (fDx2 + fDy1*fTAlph)*std::sin(fPhiTwist/2.) ;
00885 z = -fDz;
00886
00887 SetCorner(sC0Max1Min, x, y, z);
00888
00889
00890 x = fdeltaX/2. + (fDx4 + fDy2*fTAlph)*std::cos(fPhiTwist/2.) - fDy2*std::sin(fPhiTwist/2.) ;
00891 y = fdeltaY/2. + fDy2*std::cos(fPhiTwist/2.) + (fDx4 + fDy2*fTAlph)*std::sin(fPhiTwist/2.) ;
00892 z = fDz ;
00893
00894 SetCorner(sC0Max1Max, x, y, z);
00895
00896
00897 x = fdeltaX/2. + (-fDx4 + fDy2*fTAlph)*std::cos(fPhiTwist/2.) - fDy2*std::sin(fPhiTwist/2.) ;
00898 y = fdeltaY/2. + fDy2*std::cos(fPhiTwist/2.) + (-fDx4 + fDy2*fTAlph)*std::sin(fPhiTwist/2.) ;
00899 z = fDz ;
00900
00901 SetCorner(sC0Min1Max, x, y, z);
00902
00903 } else {
00904
00905 G4Exception("G4TwistTrapParallelSide::SetCorners()",
00906 "GeomSolids0001", FatalException,
00907 "Method NOT implemented !");
00908 }
00909 }
00910
00911
00912
00913
00914 void G4TwistTrapParallelSide::SetBoundaries()
00915 {
00916
00917
00918
00919 G4ThreeVector direction;
00920
00921 if (fAxis[0] == kXAxis && fAxis[1] == kZAxis) {
00922
00923
00924 direction = GetCorner(sC0Min1Max) - GetCorner(sC0Min1Min);
00925 direction = direction.unit();
00926 SetBoundary(sAxis0 & (sAxisX | sAxisMin), direction,
00927 GetCorner(sC0Min1Min), sAxisZ) ;
00928
00929
00930 direction = GetCorner(sC0Max1Max) - GetCorner(sC0Max1Min);
00931 direction = direction.unit();
00932 SetBoundary(sAxis0 & (sAxisX | sAxisMax), direction,
00933 GetCorner(sC0Max1Min), sAxisZ);
00934
00935
00936 direction = GetCorner(sC0Max1Min) - GetCorner(sC0Min1Min);
00937 direction = direction.unit();
00938 SetBoundary(sAxis1 & (sAxisZ | sAxisMin), direction,
00939 GetCorner(sC0Min1Min), sAxisX);
00940
00941
00942 direction = GetCorner(sC0Max1Max) - GetCorner(sC0Min1Max);
00943 direction = direction.unit();
00944 SetBoundary(sAxis1 & (sAxisZ | sAxisMax), direction,
00945 GetCorner(sC0Min1Max), sAxisX);
00946
00947 } else {
00948
00949 G4Exception("G4TwistTrapParallelSide::SetCorners()",
00950 "GeomSolids0001", FatalException,
00951 "Feature NOT implemented !");
00952 }
00953
00954 }
00955
00956
00957
00958
00959 void
00960 G4TwistTrapParallelSide::GetPhiUAtX( G4ThreeVector p, G4double &phi, G4double &u)
00961 {
00962
00963
00964
00965
00966
00967 phi = p.z()/(2*fDz)*fPhiTwist ;
00968
00969 u = ((-(fdeltaX*phi) + fPhiTwist*p.x())* std::cos(phi) + (-(fdeltaY*phi) + fPhiTwist*p.y())*std::sin(phi))/fPhiTwist ;
00970
00971 }
00972
00973
00974
00975
00976 G4ThreeVector G4TwistTrapParallelSide::ProjectPoint(const G4ThreeVector &p,
00977 G4bool isglobal)
00978 {
00979
00980 G4ThreeVector tmpp;
00981 if (isglobal) {
00982 tmpp = fRot.inverse()*p - fTrans;
00983 } else {
00984 tmpp = p;
00985 }
00986
00987 G4double phi ;
00988 G4double u ;
00989
00990 GetPhiUAtX( tmpp, phi, u ) ;
00991
00992 G4ThreeVector xx = SurfacePoint(phi,u) ;
00993
00994 if (isglobal) {
00995 return (fRot * xx + fTrans);
00996 } else {
00997 return xx;
00998 }
00999 }
01000
01001
01002
01003
01004 void G4TwistTrapParallelSide::GetFacets( G4int k, G4int n, G4double xyz[][3],
01005 G4int faces[][4], G4int iside )
01006 {
01007
01008 G4double phi ;
01009 G4double z, u ;
01010 G4ThreeVector p ;
01011
01012 G4int nnode ;
01013 G4int nface ;
01014
01015 G4double umin, umax ;
01016
01017
01018
01019 G4int i,j ;
01020
01021 for ( i = 0 ; i<n ; i++ ) {
01022
01023 z = -fDz+i*(2.*fDz)/(n-1) ;
01024 phi = z*fPhiTwist/(2*fDz) ;
01025 umin = GetBoundaryMin(phi) ;
01026 umax = GetBoundaryMax(phi) ;
01027
01028 for ( j = 0 ; j<k ; j++ ) {
01029
01030 nnode = GetNode(i,j,k,n,iside) ;
01031 u = umax - j*(umax-umin)/(k-1) ;
01032 p = SurfacePoint(phi,u,true) ;
01033
01034 xyz[nnode][0] = p.x() ;
01035 xyz[nnode][1] = p.y() ;
01036 xyz[nnode][2] = p.z() ;
01037
01038 if ( i<n-1 && j<k-1 ) {
01039
01040 nface = GetFace(i,j,k,n,iside) ;
01041 faces[nface][0] = GetEdgeVisibility(i,j,k,n,0,-1) * (GetNode(i ,j ,k,n,iside)+1) ;
01042 faces[nface][1] = GetEdgeVisibility(i,j,k,n,1,-1) * (GetNode(i ,j+1,k,n,iside)+1) ;
01043 faces[nface][2] = GetEdgeVisibility(i,j,k,n,2,-1) * (GetNode(i+1,j+1,k,n,iside)+1) ;
01044 faces[nface][3] = GetEdgeVisibility(i,j,k,n,3,-1) * (GetNode(i+1,j ,k,n,iside)+1) ;
01045
01046 }
01047 }
01048 }
01049 }