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
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057 #include "G4Torus.hh"
00058
00059 #include "G4VoxelLimits.hh"
00060 #include "G4AffineTransform.hh"
00061 #include "G4GeometryTolerance.hh"
00062 #include "G4JTPolynomialSolver.hh"
00063
00064 #include "G4VPVParameterisation.hh"
00065
00066 #include "meshdefs.hh"
00067
00068 #include "Randomize.hh"
00069
00070 #include "G4VGraphicsScene.hh"
00071 #include "G4Polyhedron.hh"
00072 #include "G4NURBS.hh"
00073 #include "G4NURBStube.hh"
00074 #include "G4NURBScylinder.hh"
00075 #include "G4NURBStubesector.hh"
00076
00077 using namespace CLHEP;
00078
00080
00081
00082
00083
00084 G4Torus::G4Torus( const G4String &pName,
00085 G4double pRmin,
00086 G4double pRmax,
00087 G4double pRtor,
00088 G4double pSPhi,
00089 G4double pDPhi)
00090 : G4CSGSolid(pName)
00091 {
00092 SetAllParameters(pRmin, pRmax, pRtor, pSPhi, pDPhi);
00093 }
00094
00096
00097
00098
00099 void
00100 G4Torus::SetAllParameters( G4double pRmin,
00101 G4double pRmax,
00102 G4double pRtor,
00103 G4double pSPhi,
00104 G4double pDPhi )
00105 {
00106 const G4double fEpsilon = 4.e-11;
00107
00108 fCubicVolume = 0.;
00109 fSurfaceArea = 0.;
00110 fpPolyhedron = 0;
00111
00112 kRadTolerance = G4GeometryTolerance::GetInstance()->GetRadialTolerance();
00113 kAngTolerance = G4GeometryTolerance::GetInstance()->GetAngularTolerance();
00114
00115 if ( pRtor >= pRmax+1.e3*kCarTolerance )
00116 {
00117 fRtor = pRtor ;
00118 }
00119 else
00120 {
00121 std::ostringstream message;
00122 message << "Invalid swept radius for Solid: " << GetName() << G4endl
00123 << " pRtor = " << pRtor << ", pRmax = " << pRmax;
00124 G4Exception("G4Torus::SetAllParameters()",
00125 "GeomSolids0002", FatalException, message);
00126 }
00127
00128
00129
00130 if ( pRmin < pRmax - 1.e2*kCarTolerance && pRmin >= 0 )
00131 {
00132 if (pRmin >= 1.e2*kCarTolerance) { fRmin = pRmin ; }
00133 else { fRmin = 0.0 ; }
00134 fRmax = pRmax ;
00135 }
00136 else
00137 {
00138 std::ostringstream message;
00139 message << "Invalid values of radii for Solid: " << GetName() << G4endl
00140 << " pRmin = " << pRmin << ", pRmax = " << pRmax;
00141 G4Exception("G4Torus::SetAllParameters()",
00142 "GeomSolids0002", FatalException, message);
00143 }
00144
00145
00146
00147 fRminTolerance = (fRmin)
00148 ? 0.5*std::max( kRadTolerance, fEpsilon*(fRtor-fRmin )) : 0;
00149 fRmaxTolerance = 0.5*std::max( kRadTolerance, fEpsilon*(fRtor+fRmax) );
00150
00151
00152
00153 if ( pDPhi >= twopi ) { fDPhi = twopi ; }
00154 else
00155 {
00156 if (pDPhi > 0) { fDPhi = pDPhi ; }
00157 else
00158 {
00159 std::ostringstream message;
00160 message << "Invalid Z delta-Phi for Solid: " << GetName() << G4endl
00161 << " pDPhi = " << pDPhi;
00162 G4Exception("G4Torus::SetAllParameters()",
00163 "GeomSolids0002", FatalException, message);
00164 }
00165 }
00166
00167
00168
00169 fSPhi = pSPhi;
00170
00171 if (fSPhi < 0) { fSPhi = twopi-std::fmod(std::fabs(fSPhi),twopi) ; }
00172 else { fSPhi = std::fmod(fSPhi,twopi) ; }
00173
00174 if (fSPhi+fDPhi > twopi) { fSPhi-=twopi ; }
00175 }
00176
00178
00179
00180
00181
00182 G4Torus::G4Torus( __void__& a )
00183 : G4CSGSolid(a), fRmin(0.), fRmax(0.), fRtor(0.), fSPhi(0.),
00184 fDPhi(0.), fRminTolerance(0.), fRmaxTolerance(0. ),
00185 kRadTolerance(0.), kAngTolerance(0.)
00186 {
00187 }
00188
00190
00191
00192
00193 G4Torus::~G4Torus()
00194 {}
00195
00197
00198
00199
00200 G4Torus::G4Torus(const G4Torus& rhs)
00201 : G4CSGSolid(rhs), fRmin(rhs.fRmin),fRmax(rhs.fRmax),
00202 fRtor(rhs.fRtor),fSPhi(rhs.fSPhi),fDPhi(rhs.fDPhi),
00203 fRminTolerance(rhs.fRminTolerance), fRmaxTolerance(rhs.fRmaxTolerance),
00204 kRadTolerance(rhs.kRadTolerance), kAngTolerance(rhs.kAngTolerance)
00205 {
00206 }
00207
00209
00210
00211
00212 G4Torus& G4Torus::operator = (const G4Torus& rhs)
00213 {
00214
00215
00216 if (this == &rhs) { return *this; }
00217
00218
00219
00220 G4CSGSolid::operator=(rhs);
00221
00222
00223
00224 fRmin = rhs.fRmin; fRmax = rhs.fRmax;
00225 fRtor = rhs.fRtor; fSPhi = rhs.fSPhi; fDPhi = rhs.fDPhi;
00226 fRminTolerance = rhs.fRminTolerance; fRmaxTolerance = rhs.fRmaxTolerance;
00227 kRadTolerance = rhs.kRadTolerance; kAngTolerance = rhs.kAngTolerance;
00228
00229 return *this;
00230 }
00231
00233
00234
00235
00236
00237 void G4Torus::ComputeDimensions( G4VPVParameterisation* p,
00238 const G4int n,
00239 const G4VPhysicalVolume* pRep )
00240 {
00241 p->ComputeDimensions(*this,n,pRep);
00242 }
00243
00244
00245
00247
00248
00249
00250
00251 void G4Torus::TorusRootsJT( const G4ThreeVector& p,
00252 const G4ThreeVector& v,
00253 G4double r,
00254 std::vector<G4double>& roots ) const
00255 {
00256
00257 G4int i, num ;
00258 G4double c[5], srd[4], si[4] ;
00259
00260 G4double Rtor2 = fRtor*fRtor, r2 = r*r ;
00261
00262 G4double pDotV = p.x()*v.x() + p.y()*v.y() + p.z()*v.z() ;
00263 G4double pRad2 = p.x()*p.x() + p.y()*p.y() + p.z()*p.z() ;
00264
00265 c[0] = 1.0 ;
00266 c[1] = 4*pDotV ;
00267 c[2] = 2*(pRad2 + 2*pDotV*pDotV - Rtor2 - r2 + 2*Rtor2*v.z()*v.z()) ;
00268 c[3] = 4*(pDotV*(pRad2 - Rtor2 - r2) + 2*Rtor2*p.z()*v.z()) ;
00269 c[4] = pRad2*pRad2 - 2*pRad2*(Rtor2+r2)
00270 + 4*Rtor2*p.z()*p.z() + (Rtor2-r2)*(Rtor2-r2) ;
00271
00272 G4JTPolynomialSolver torusEq;
00273
00274 num = torusEq.FindRoots( c, 4, srd, si );
00275
00276 for ( i = 0; i < num; i++ )
00277 {
00278 if( si[i] == 0. ) { roots.push_back(srd[i]) ; }
00279 }
00280
00281 std::sort(roots.begin() , roots.end() ) ;
00282 }
00283
00285
00286
00287
00288
00289
00290
00291 G4double G4Torus::SolveNumericJT( const G4ThreeVector& p,
00292 const G4ThreeVector& v,
00293 G4double r,
00294 G4bool IsDistanceToIn ) const
00295 {
00296 G4double bigdist = 10*mm ;
00297 G4double tmin = kInfinity ;
00298 G4double t, scal ;
00299 static const G4double halfCarTolerance = 0.5*kCarTolerance;
00300 static const G4double halfAngTolerance = 0.5*kAngTolerance;
00301
00302
00303
00304
00305 std::vector<G4double> roots ;
00306 std::vector<G4double> rootsrefined ;
00307 TorusRootsJT(p,v,r,roots) ;
00308
00309 G4ThreeVector ptmp ;
00310
00311
00312
00313 for ( size_t k = 0 ; k<roots.size() ; k++ )
00314 {
00315 t = roots[k] ;
00316
00317 if ( t < -halfCarTolerance ) { continue ; }
00318
00319 if ( t > bigdist && t<kInfinity )
00320 {
00321 ptmp = p + t*v ;
00322 TorusRootsJT(ptmp,v,r,rootsrefined) ;
00323 if ( rootsrefined.size()==roots.size() )
00324 {
00325 t = t + rootsrefined[k] ;
00326 }
00327 }
00328
00329 ptmp = p + t*v ;
00330
00331 G4double theta = std::atan2(ptmp.y(),ptmp.x());
00332
00333 if ( fSPhi >= 0 )
00334 {
00335 if ( theta < - halfAngTolerance ) { theta += twopi; }
00336 if ( (std::fabs(theta) < halfAngTolerance)
00337 && (std::fabs(fSPhi + fDPhi - twopi) < halfAngTolerance) )
00338 {
00339 theta += twopi ;
00340 }
00341 }
00342 if ((fSPhi <= -pi )&&(theta>halfAngTolerance)) { theta = theta-twopi; }
00343
00344
00345
00346
00347 if ( (theta - fSPhi >= - halfAngTolerance)
00348 && (theta - (fSPhi + fDPhi) <= halfAngTolerance) )
00349 {
00350
00351
00352
00353 if ( IsDistanceToIn == true )
00354 {
00355 if (std::fabs(t) < halfCarTolerance )
00356 {
00357
00358
00359
00360 scal = v* G4ThreeVector( p.x()*(1-fRtor/std::sqrt(p.x()*p.x()
00361 + p.y()*p.y())),
00362 p.y()*(1-fRtor/std::sqrt(p.x()*p.x()
00363 + p.y()*p.y())),
00364 p.z() );
00365
00366
00367
00368 if ( r == GetRmin() ) { scal = -scal ; }
00369 if ( scal < 0 ) { return 0.0 ; }
00370 }
00371 }
00372
00373
00374
00375
00376 if ( IsDistanceToIn == false )
00377 {
00378 if (std::fabs(t) < halfCarTolerance )
00379 {
00380
00381
00382 scal = v* G4ThreeVector( p.x()*(1-fRtor/std::sqrt(p.x()*p.x()
00383 + p.y()*p.y())),
00384 p.y()*(1-fRtor/std::sqrt(p.x()*p.x()
00385 + p.y()*p.y())),
00386 p.z() );
00387
00388
00389
00390 if ( r == GetRmin() ) { scal = -scal ; }
00391 if ( scal > 0 ) { return 0.0 ; }
00392 }
00393 }
00394
00395
00396
00397 if( t > halfCarTolerance )
00398 {
00399 tmin = t ;
00400 return tmin ;
00401 }
00402 }
00403 }
00404
00405 return tmin;
00406 }
00407
00409
00410
00411
00412 G4bool G4Torus::CalculateExtent( const EAxis pAxis,
00413 const G4VoxelLimits& pVoxelLimit,
00414 const G4AffineTransform& pTransform,
00415 G4double& pMin, G4double& pMax) const
00416 {
00417 if ((!pTransform.IsRotated()) && (fDPhi==twopi) && (fRmin==0))
00418 {
00419
00420
00421
00422
00423
00424 G4double xoffset,xMin,xMax;
00425 G4double yoffset,yMin,yMax;
00426 G4double zoffset,zMin,zMax;
00427
00428 G4double RTorus,delta,diff1,diff2,maxDiff,newMin,newMax;
00429 G4double xoff1,xoff2,yoff1,yoff2;
00430
00431 xoffset = pTransform.NetTranslation().x();
00432 xMin = xoffset - fRmax - fRtor ;
00433 xMax = xoffset + fRmax + fRtor ;
00434
00435 if (pVoxelLimit.IsXLimited())
00436 {
00437 if ( (xMin > pVoxelLimit.GetMaxXExtent()+kCarTolerance)
00438 || (xMax < pVoxelLimit.GetMinXExtent()-kCarTolerance) )
00439 return false ;
00440 else
00441 {
00442 if (xMin < pVoxelLimit.GetMinXExtent())
00443 {
00444 xMin = pVoxelLimit.GetMinXExtent() ;
00445 }
00446 if (xMax > pVoxelLimit.GetMaxXExtent())
00447 {
00448 xMax = pVoxelLimit.GetMaxXExtent() ;
00449 }
00450 }
00451 }
00452 yoffset = pTransform.NetTranslation().y();
00453 yMin = yoffset - fRmax - fRtor ;
00454 yMax = yoffset + fRmax + fRtor ;
00455
00456 if (pVoxelLimit.IsYLimited())
00457 {
00458 if ( (yMin > pVoxelLimit.GetMaxYExtent()+kCarTolerance)
00459 || (yMax < pVoxelLimit.GetMinYExtent()-kCarTolerance) )
00460 {
00461 return false ;
00462 }
00463 else
00464 {
00465 if (yMin < pVoxelLimit.GetMinYExtent() )
00466 {
00467 yMin = pVoxelLimit.GetMinYExtent() ;
00468 }
00469 if (yMax > pVoxelLimit.GetMaxYExtent() )
00470 {
00471 yMax = pVoxelLimit.GetMaxYExtent() ;
00472 }
00473 }
00474 }
00475 zoffset = pTransform.NetTranslation().z() ;
00476 zMin = zoffset - fRmax ;
00477 zMax = zoffset + fRmax ;
00478
00479 if (pVoxelLimit.IsZLimited())
00480 {
00481 if ( (zMin > pVoxelLimit.GetMaxZExtent()+kCarTolerance)
00482 || (zMax < pVoxelLimit.GetMinZExtent()-kCarTolerance) )
00483 {
00484 return false ;
00485 }
00486 else
00487 {
00488 if (zMin < pVoxelLimit.GetMinZExtent() )
00489 {
00490 zMin = pVoxelLimit.GetMinZExtent() ;
00491 }
00492 if (zMax > pVoxelLimit.GetMaxZExtent() )
00493 {
00494 zMax = pVoxelLimit.GetMaxZExtent() ;
00495 }
00496 }
00497 }
00498
00499
00500
00501 switch (pAxis)
00502 {
00503 case kXAxis:
00504 yoff1=yoffset-yMin;
00505 yoff2=yMax-yoffset;
00506 if ( yoff1 >= 0 && yoff2 >= 0 )
00507 {
00508
00509
00510 pMin = xMin ;
00511 pMax = xMax ;
00512 }
00513 else
00514 {
00515
00516
00517
00518
00519 RTorus=fRmax+fRtor;
00520 delta = RTorus*RTorus - yoff1*yoff1;
00521 diff1 = (delta>0.) ? std::sqrt(delta) : 0.;
00522 delta = RTorus*RTorus - yoff2*yoff2;
00523 diff2 = (delta>0.) ? std::sqrt(delta) : 0.;
00524 maxDiff = (diff1 > diff2) ? diff1:diff2 ;
00525 newMin = xoffset - maxDiff ;
00526 newMax = xoffset + maxDiff ;
00527 pMin = (newMin < xMin) ? xMin : newMin ;
00528 pMax = (newMax > xMax) ? xMax : newMax ;
00529 }
00530 break;
00531
00532 case kYAxis:
00533 xoff1 = xoffset - xMin ;
00534 xoff2 = xMax - xoffset ;
00535 if (xoff1 >= 0 && xoff2 >= 0 )
00536 {
00537
00538
00539 pMin = yMin ;
00540 pMax = yMax ;
00541 }
00542 else
00543 {
00544
00545
00546
00547 RTorus=fRmax+fRtor;
00548 delta = RTorus*RTorus - xoff1*xoff1;
00549 diff1 = (delta>0.) ? std::sqrt(delta) : 0.;
00550 delta = RTorus*RTorus - xoff2*xoff2;
00551 diff2 = (delta>0.) ? std::sqrt(delta) : 0.;
00552 maxDiff = (diff1 > diff2) ? diff1 : diff2 ;
00553 newMin = yoffset - maxDiff ;
00554 newMax = yoffset + maxDiff ;
00555 pMin = (newMin < yMin) ? yMin : newMin ;
00556 pMax = (newMax > yMax) ? yMax : newMax ;
00557 }
00558 break;
00559
00560 case kZAxis:
00561 pMin=zMin;
00562 pMax=zMax;
00563 break;
00564
00565 default:
00566 break;
00567 }
00568 pMin -= kCarTolerance ;
00569 pMax += kCarTolerance ;
00570
00571 return true;
00572 }
00573 else
00574 {
00575 G4int i, noEntries, noBetweenSections4 ;
00576 G4bool existsAfterClip = false ;
00577
00578
00579
00580 G4ThreeVectorList *vertices ;
00581 G4int noPolygonVertices ;
00582 vertices = CreateRotatedVertices(pTransform,noPolygonVertices) ;
00583
00584 pMin = +kInfinity ;
00585 pMax = -kInfinity ;
00586
00587 noEntries = vertices->size() ;
00588 noBetweenSections4 = noEntries - noPolygonVertices ;
00589
00590 for (i=0;i<noEntries;i+=noPolygonVertices)
00591 {
00592 ClipCrossSection(vertices,i,pVoxelLimit,pAxis,pMin,pMax);
00593 }
00594 for (i=0;i<noBetweenSections4;i+=noPolygonVertices)
00595 {
00596 ClipBetweenSections(vertices,i,pVoxelLimit,pAxis,pMin,pMax);
00597 }
00598 if (pMin!=kInfinity||pMax!=-kInfinity)
00599 {
00600 existsAfterClip = true ;
00601 pMin -= kCarTolerance ;
00602 pMax += kCarTolerance ;
00603 }
00604 else
00605 {
00606
00607
00608
00609
00610
00611 G4ThreeVector clipCentre(
00612 (pVoxelLimit.GetMinXExtent()+pVoxelLimit.GetMaxXExtent())*0.5,
00613 (pVoxelLimit.GetMinYExtent()+pVoxelLimit.GetMaxYExtent())*0.5,
00614 (pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5 ) ;
00615
00616 if (Inside(pTransform.Inverse().TransformPoint(clipCentre)) != kOutside )
00617 {
00618 existsAfterClip = true ;
00619 pMin = pVoxelLimit.GetMinExtent(pAxis) ;
00620 pMax = pVoxelLimit.GetMaxExtent(pAxis) ;
00621 }
00622 }
00623 delete vertices;
00624 return existsAfterClip;
00625 }
00626 }
00627
00629
00630
00631
00632 EInside G4Torus::Inside( const G4ThreeVector& p ) const
00633 {
00634 G4double r2, pt2, pPhi, tolRMin, tolRMax ;
00635
00636 EInside in = kOutside ;
00637 static const G4double halfAngTolerance = 0.5*kAngTolerance;
00638
00639
00640 r2 = p.x()*p.x() + p.y()*p.y() ;
00641 pt2 = r2 + p.z()*p.z() + fRtor*fRtor - 2*fRtor*std::sqrt(r2) ;
00642
00643 if (fRmin) tolRMin = fRmin + fRminTolerance ;
00644 else tolRMin = 0 ;
00645
00646 tolRMax = fRmax - fRmaxTolerance;
00647
00648 if (pt2 >= tolRMin*tolRMin && pt2 <= tolRMax*tolRMax )
00649 {
00650 if ( fDPhi == twopi || pt2 == 0 )
00651 {
00652 in = kInside ;
00653 }
00654 else
00655 {
00656
00657
00658
00659 pPhi = std::atan2(p.y(),p.x()) ;
00660
00661 if ( pPhi < -halfAngTolerance ) { pPhi += twopi ; }
00662 if ( fSPhi >= 0 )
00663 {
00664 if ( (std::fabs(pPhi) < halfAngTolerance)
00665 && (std::fabs(fSPhi + fDPhi - twopi) < halfAngTolerance) )
00666 {
00667 pPhi += twopi ;
00668 }
00669 if ( (pPhi >= fSPhi + halfAngTolerance)
00670 && (pPhi <= fSPhi + fDPhi - halfAngTolerance) )
00671 {
00672 in = kInside ;
00673 }
00674 else if ( (pPhi >= fSPhi - halfAngTolerance)
00675 && (pPhi <= fSPhi + fDPhi + halfAngTolerance) )
00676 {
00677 in = kSurface ;
00678 }
00679 }
00680 else
00681 {
00682 if ( (pPhi <= fSPhi + twopi - halfAngTolerance)
00683 && (pPhi >= fSPhi + fDPhi + halfAngTolerance) ) {;}
00684 else
00685 {
00686 in = kSurface ;
00687 }
00688 }
00689 }
00690 }
00691 else
00692 {
00693 tolRMin = fRmin - fRminTolerance ;
00694 tolRMax = fRmax + fRmaxTolerance ;
00695
00696 if (tolRMin < 0 ) { tolRMin = 0 ; }
00697
00698 if ( (pt2 >= tolRMin*tolRMin) && (pt2 <= tolRMax*tolRMax) )
00699 {
00700 if ( (fDPhi == twopi) || (pt2 == 0) )
00701 {
00702 in = kSurface ;
00703 }
00704 else
00705 {
00706 pPhi = std::atan2(p.y(),p.x()) ;
00707
00708 if ( pPhi < -halfAngTolerance ) { pPhi += twopi ; }
00709 if ( fSPhi >= 0 )
00710 {
00711 if ( (std::fabs(pPhi) < halfAngTolerance)
00712 && (std::fabs(fSPhi + fDPhi - twopi) < halfAngTolerance) )
00713 {
00714 pPhi += twopi ;
00715 }
00716 if ( (pPhi >= fSPhi - halfAngTolerance)
00717 && (pPhi <= fSPhi + fDPhi + halfAngTolerance) )
00718 {
00719 in = kSurface;
00720 }
00721 }
00722 else
00723 {
00724 if ( (pPhi <= fSPhi + twopi - halfAngTolerance)
00725 && (pPhi >= fSPhi + fDPhi + halfAngTolerance) ) {;}
00726 else
00727 {
00728 in = kSurface ;
00729 }
00730 }
00731 }
00732 }
00733 }
00734 return in ;
00735 }
00736
00738
00739
00740
00741
00742
00743 G4ThreeVector G4Torus::SurfaceNormal( const G4ThreeVector& p ) const
00744 {
00745 G4int noSurfaces = 0;
00746 G4double rho2, rho, pt2, pt, pPhi;
00747 G4double distRMin = kInfinity;
00748 G4double distSPhi = kInfinity, distEPhi = kInfinity;
00749
00750
00751
00752 const G4double delta = std::max(10.0*kCarTolerance,
00753 1.0e-8*(fRtor+fRmax));
00754 const G4double dAngle = 10.0*kAngTolerance;
00755
00756 G4ThreeVector nR, nPs, nPe;
00757 G4ThreeVector norm, sumnorm(0.,0.,0.);
00758
00759 rho2 = p.x()*p.x() + p.y()*p.y();
00760 rho = std::sqrt(rho2);
00761 pt2 = rho2+p.z()*p.z() +fRtor * (fRtor-2*rho);
00762 pt2 = std::max(pt2, 0.0);
00763 pt = std::sqrt(pt2) ;
00764
00765 G4double distRMax = std::fabs(pt - fRmax);
00766 if(fRmin) distRMin = std::fabs(pt - fRmin);
00767
00768 if( rho > delta && pt != 0.0 )
00769 {
00770 G4double redFactor= (rho-fRtor)/rho;
00771 nR = G4ThreeVector( p.x()*redFactor,
00772 p.y()*redFactor,
00773 p.z() );
00774 nR *= 1.0/pt;
00775 }
00776
00777 if ( fDPhi < twopi )
00778 {
00779 if ( rho )
00780 {
00781 pPhi = std::atan2(p.y(),p.x());
00782
00783 if(pPhi < fSPhi-delta) { pPhi += twopi; }
00784 else if(pPhi > fSPhi+fDPhi+delta) { pPhi -= twopi; }
00785
00786 distSPhi = std::fabs( pPhi - fSPhi );
00787 distEPhi = std::fabs(pPhi-fSPhi-fDPhi);
00788 }
00789 nPs = G4ThreeVector(std::sin(fSPhi),-std::cos(fSPhi),0);
00790 nPe = G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0);
00791 }
00792 if( distRMax <= delta )
00793 {
00794 noSurfaces ++;
00795 sumnorm += nR;
00796 }
00797 else if( fRmin && (distRMin <= delta) )
00798 {
00799 noSurfaces ++;
00800 sumnorm -= nR;
00801 }
00802
00803
00804
00805
00806 if( (fDPhi < twopi) && (fRmin-delta <= pt) && (pt <= (fRmax+delta)) )
00807 {
00808 if (distSPhi <= dAngle)
00809 {
00810 noSurfaces ++;
00811 sumnorm += nPs;
00812 }
00813 if (distEPhi <= dAngle)
00814 {
00815 noSurfaces ++;
00816 sumnorm += nPe;
00817 }
00818 }
00819 if ( noSurfaces == 0 )
00820 {
00821 #ifdef G4CSGDEBUG
00822 G4ExceptionDescription ed;
00823 ed.precision(16);
00824
00825 EInside inIt= Inside( p );
00826
00827 if( inIt != kSurface )
00828 {
00829 ed << " ERROR> Surface Normal was called for Torus,"
00830 << " with point not on surface." << G4endl;
00831 }
00832 else
00833 {
00834 ed << " ERROR> Surface Normal has not found a surface, "
00835 << " despite the point being on the surface. " <<G4endl;
00836 }
00837
00838 if( inIt != kInside)
00839 {
00840 ed << " Safety (Dist To In) = " << DistanceToIn(p) << G4endl;
00841 }
00842 if( inIt != kOutside)
00843 {
00844 ed << " Safety (Dist to Out) = " << DistanceToOut(p) << G4endl;
00845 }
00846 ed << " Coordinates of point : " << p << G4endl;
00847 ed << " Parameters of solid : " << G4endl << *this << G4endl;
00848
00849 if( inIt == kSurface )
00850 {
00851 G4Exception("G4Torus::SurfaceNormal(p)", "GeomSolids1002",
00852 JustWarning, ed,
00853 "Failing to find normal, even though point is on surface!" );
00854 }
00855 else
00856 {
00857 static const char* NameInside[3]= { "Inside", "Surface", "Outside" };
00858 ed << " The point is " << NameInside[inIt] << " the solid. "<< G4endl;
00859 G4Exception("G4Torus::SurfaceNormal(p)", "GeomSolids1002",
00860 JustWarning, ed, "Point p is not on surface !?" );
00861 }
00862 #endif
00863 norm = ApproxSurfaceNormal(p);
00864 }
00865 else if ( noSurfaces == 1 ) { norm = sumnorm; }
00866 else { norm = sumnorm.unit(); }
00867
00868
00869
00870 return norm ;
00871 }
00872
00874
00875
00876
00877
00878 G4ThreeVector G4Torus::ApproxSurfaceNormal( const G4ThreeVector& p ) const
00879 {
00880 ENorm side ;
00881 G4ThreeVector norm;
00882 G4double rho2,rho,pt2,pt,phi;
00883 G4double distRMin,distRMax,distSPhi,distEPhi,distMin;
00884
00885 rho2 = p.x()*p.x() + p.y()*p.y();
00886 rho = std::sqrt(rho2) ;
00887 pt2 = std::fabs(rho2+p.z()*p.z() +fRtor*fRtor - 2*fRtor*rho) ;
00888 pt = std::sqrt(pt2) ;
00889
00890 #ifdef G4CSGDEBUG
00891 G4cout << " G4Torus::ApproximateSurfaceNormal called for point " << p
00892 << G4endl;
00893 #endif
00894
00895 distRMax = std::fabs(pt - fRmax) ;
00896
00897 if(fRmin)
00898 {
00899 distRMin = std::fabs(pt - fRmin) ;
00900
00901 if (distRMin < distRMax)
00902 {
00903 distMin = distRMin ;
00904 side = kNRMin ;
00905 }
00906 else
00907 {
00908 distMin = distRMax ;
00909 side = kNRMax ;
00910 }
00911 }
00912 else
00913 {
00914 distMin = distRMax ;
00915 side = kNRMax ;
00916 }
00917 if ( (fDPhi < twopi) && rho )
00918 {
00919 phi = std::atan2(p.y(),p.x()) ;
00920
00921 if (phi < 0) { phi += twopi ; }
00922
00923 if (fSPhi < 0 ) { distSPhi = std::fabs(phi-(fSPhi+twopi))*rho ; }
00924 else { distSPhi = std::fabs(phi-fSPhi)*rho ; }
00925
00926 distEPhi = std::fabs(phi - fSPhi - fDPhi)*rho ;
00927
00928 if (distSPhi < distEPhi)
00929 {
00930 if (distSPhi<distMin) side = kNSPhi ;
00931 }
00932 else
00933 {
00934 if (distEPhi < distMin) { side = kNEPhi ; }
00935 }
00936 }
00937 switch (side)
00938 {
00939 case kNRMin:
00940 norm = G4ThreeVector( -p.x()*(1-fRtor/rho)/pt,
00941 -p.y()*(1-fRtor/rho)/pt,
00942 -p.z()/pt ) ;
00943 break ;
00944 case kNRMax:
00945 norm = G4ThreeVector( p.x()*(1-fRtor/rho)/pt,
00946 p.y()*(1-fRtor/rho)/pt,
00947 p.z()/pt ) ;
00948 break;
00949 case kNSPhi:
00950 norm = G4ThreeVector(std::sin(fSPhi),-std::cos(fSPhi),0) ;
00951 break;
00952 case kNEPhi:
00953 norm = G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0) ;
00954 break;
00955 default:
00956 DumpInfo();
00957 G4Exception("G4Torus::ApproxSurfaceNormal()",
00958 "GeomSolids1002", JustWarning,
00959 "Undefined side for valid surface normal to solid.");
00960 break ;
00961 }
00962 return norm ;
00963 }
00964
00966
00967
00968
00969
00970
00971
00972
00973
00974
00975
00976
00977
00978
00979
00980
00981
00982
00983
00984
00985
00986
00987 G4double G4Torus::DistanceToIn( const G4ThreeVector& p,
00988 const G4ThreeVector& v ) const
00989 {
00990
00991 G4double snxt=kInfinity, sphi=kInfinity;
00992
00993 G4double sd[4] ;
00994
00995
00996
00997
00998 G4bool seg;
00999 G4double hDPhi;
01000 G4double cPhi,sinCPhi=0.,cosCPhi=0.;
01001
01002 G4double tolORMin2;
01003 G4double tolORMax2;
01004
01005 static const G4double halfCarTolerance = 0.5*kCarTolerance;
01006
01007 G4double Dist,xi,yi,zi,rhoi2,it2;
01008
01009 G4double Comp;
01010 G4double cosSPhi,sinSPhi;
01011 G4double ePhi,cosEPhi,sinEPhi;
01012
01013
01014
01015 if ( fDPhi < twopi )
01016 {
01017 seg = true ;
01018 hDPhi = 0.5*fDPhi ;
01019 cPhi = fSPhi + hDPhi ;
01020 sinCPhi = std::sin(cPhi) ;
01021 cosCPhi = std::cos(cPhi) ;
01022 }
01023 else
01024 {
01025 seg = false ;
01026 }
01027
01028 if (fRmin > fRminTolerance)
01029 {
01030 tolORMin2 = (fRmin - fRminTolerance)*(fRmin - fRminTolerance) ;
01031 }
01032 else
01033 {
01034 tolORMin2 = 0 ;
01035 }
01036 tolORMax2 = (fRmax + fRmaxTolerance)*(fRmax + fRmaxTolerance) ;
01037
01038
01039
01040 G4double Rtor2 = fRtor*fRtor ;
01041
01042 snxt = SolveNumericJT(p,v,fRmax,true);
01043
01044 if (fRmin)
01045 {
01046 sd[0] = SolveNumericJT(p,v,fRmin,true);
01047 if ( sd[0] < snxt ) { snxt = sd[0] ; }
01048 }
01049
01050
01051
01052
01053
01054
01055
01056
01057
01058
01059
01060 if (seg)
01061 {
01062 sinSPhi = std::sin(fSPhi) ;
01063 cosSPhi = std::cos(fSPhi) ;
01064 Comp = v.x()*sinSPhi - v.y()*cosSPhi ;
01065
01066 if (Comp < 0 )
01067 {
01068 Dist = (p.y()*cosSPhi - p.x()*sinSPhi) ;
01069
01070 if (Dist < halfCarTolerance)
01071 {
01072 sphi = Dist/Comp ;
01073 if (sphi < snxt)
01074 {
01075 if ( sphi < 0 ) { sphi = 0 ; }
01076
01077 xi = p.x() + sphi*v.x() ;
01078 yi = p.y() + sphi*v.y() ;
01079 zi = p.z() + sphi*v.z() ;
01080 rhoi2 = xi*xi + yi*yi ;
01081 it2 = std::fabs(rhoi2 + zi*zi + Rtor2 - 2*fRtor*std::sqrt(rhoi2)) ;
01082
01083 if ( it2 >= tolORMin2 && it2 <= tolORMax2 )
01084 {
01085
01086
01087
01088 if ((yi*cosCPhi-xi*sinCPhi)<=0) { snxt=sphi; }
01089 }
01090 }
01091 }
01092 }
01093 ePhi=fSPhi+fDPhi;
01094 sinEPhi=std::sin(ePhi);
01095 cosEPhi=std::cos(ePhi);
01096 Comp=-(v.x()*sinEPhi-v.y()*cosEPhi);
01097
01098 if ( Comp < 0 )
01099 {
01100 Dist = -(p.y()*cosEPhi - p.x()*sinEPhi) ;
01101
01102 if (Dist < halfCarTolerance )
01103 {
01104 sphi = Dist/Comp ;
01105
01106 if (sphi < snxt )
01107 {
01108 if (sphi < 0 ) { sphi = 0 ; }
01109
01110 xi = p.x() + sphi*v.x() ;
01111 yi = p.y() + sphi*v.y() ;
01112 zi = p.z() + sphi*v.z() ;
01113 rhoi2 = xi*xi + yi*yi ;
01114 it2 = std::fabs(rhoi2 + zi*zi + Rtor2 - 2*fRtor*std::sqrt(rhoi2)) ;
01115
01116 if (it2 >= tolORMin2 && it2 <= tolORMax2)
01117 {
01118
01119
01120
01121 if ((yi*cosCPhi-xi*sinCPhi)>=0) { snxt=sphi; }
01122 }
01123 }
01124 }
01125 }
01126 }
01127 if(snxt < halfCarTolerance) { snxt = 0.0 ; }
01128
01129 return snxt ;
01130 }
01131
01133
01134
01135
01136
01137
01138
01139 G4double G4Torus::DistanceToIn( const G4ThreeVector& p ) const
01140 {
01141 G4double safe=0.0, safe1, safe2 ;
01142 G4double phiC, cosPhiC, sinPhiC, safePhi, ePhi, cosPsi ;
01143 G4double rho2, rho, pt2, pt ;
01144
01145 rho2 = p.x()*p.x() + p.y()*p.y() ;
01146 rho = std::sqrt(rho2) ;
01147 pt2 = std::fabs(rho2 + p.z()*p.z() + fRtor*fRtor - 2*fRtor*rho) ;
01148 pt = std::sqrt(pt2) ;
01149
01150 safe1 = fRmin - pt ;
01151 safe2 = pt - fRmax ;
01152
01153 if (safe1 > safe2) { safe = safe1; }
01154 else { safe = safe2; }
01155
01156 if ( fDPhi < twopi && rho )
01157 {
01158 phiC = fSPhi + fDPhi*0.5 ;
01159 cosPhiC = std::cos(phiC) ;
01160 sinPhiC = std::sin(phiC) ;
01161 cosPsi = (p.x()*cosPhiC + p.y()*sinPhiC)/rho ;
01162
01163 if (cosPsi < std::cos(fDPhi*0.5) )
01164 {
01165 if ((p.y()*cosPhiC - p.x()*sinPhiC) <= 0 )
01166 {
01167 safePhi = std::fabs(p.x()*std::sin(fSPhi) - p.y()*std::cos(fSPhi)) ;
01168 }
01169 else
01170 {
01171 ePhi = fSPhi + fDPhi ;
01172 safePhi = std::fabs(p.x()*std::sin(ePhi) - p.y()*std::cos(ePhi)) ;
01173 }
01174 if (safePhi > safe) { safe = safePhi ; }
01175 }
01176 }
01177 if (safe < 0 ) { safe = 0 ; }
01178 return safe;
01179 }
01180
01182
01183
01184
01185
01186
01187 G4double G4Torus::DistanceToOut( const G4ThreeVector& p,
01188 const G4ThreeVector& v,
01189 const G4bool calcNorm,
01190 G4bool *validNorm,
01191 G4ThreeVector *n ) const
01192 {
01193 ESide side = kNull, sidephi = kNull ;
01194 G4double snxt = kInfinity, sphi, sd[4] ;
01195
01196 static const G4double halfCarTolerance = 0.5*kCarTolerance;
01197 static const G4double halfAngTolerance = 0.5*kAngTolerance;
01198
01199
01200
01201 G4double sinSPhi, cosSPhi, ePhi, sinEPhi, cosEPhi;
01202 G4double cPhi, sinCPhi, cosCPhi ;
01203 G4double pDistS, compS, pDistE, compE, sphi2, xi, yi, zi, vphi ;
01204
01205
01206
01208
01209 #if 1
01210
01211
01212
01213
01214
01215 G4double rho2 = p.x()*p.x()+p.y()*p.y();
01216 G4double rho = std::sqrt(rho2) ;
01217
01218 G4double pt2 = rho2 + p.z()*p.z() + fRtor * (fRtor - 2.0*rho);
01219
01220
01221 if( pt2 < 0.0)
01222 {
01223 pt2= std::fabs( pt2 );
01224 }
01225
01226 G4double pt = std::sqrt(pt2) ;
01227
01228 G4double pDotV = p.x()*v.x() + p.y()*v.y() + p.z()*v.z() ;
01229
01230 G4double tolRMax = fRmax - fRmaxTolerance ;
01231
01232 G4double vDotNmax = pDotV - fRtor*(v.x()*p.x() + v.y()*p.y())/rho ;
01233 G4double pDotxyNmax = (1 - fRtor/rho) ;
01234
01235 if( (pt2 > tolRMax*tolRMax) && (vDotNmax >= 0) )
01236 {
01237
01238
01239
01240
01241 if ( calcNorm && (pDotxyNmax >= -2.*fRmaxTolerance) )
01242 {
01243 *n = G4ThreeVector( p.x()*(1 - fRtor/rho)/pt,
01244 p.y()*(1 - fRtor/rho)/pt,
01245 p.z()/pt ) ;
01246 *validNorm = true ;
01247 }
01248
01249 return snxt = 0 ;
01250 }
01251
01252 snxt = SolveNumericJT(p,v,fRmax,false);
01253 side = kRMax ;
01254
01255
01256
01257 if ( fRmin )
01258 {
01259 G4double tolRMin = fRmin + fRminTolerance ;
01260
01261 if ( (pt2 < tolRMin*tolRMin) && (vDotNmax < 0) )
01262 {
01263 if (calcNorm) { *validNorm = false ; }
01264 return snxt = 0 ;
01265 }
01266
01267 sd[0] = SolveNumericJT(p,v,fRmin,false);
01268 if ( sd[0] < snxt )
01269 {
01270 snxt = sd[0] ;
01271 side = kRMin ;
01272 }
01273 }
01274
01275 #else
01276
01277
01278
01279
01280 snxt = SolveNumericJT(p,v,fRmax,false);
01281 side = kRMax ;
01282
01283 if ( fRmin )
01284 {
01285 sd[0] = SolveNumericJT(p,v,fRmin,false);
01286 if ( sd[0] < snxt )
01287 {
01288 snxt = sd[0] ;
01289 side = kRMin ;
01290 }
01291 }
01292
01293 if ( calcNorm && (snxt == 0.0) )
01294 {
01295 *validNorm = false ;
01296 return snxt ;
01297 }
01298
01299 #endif
01300
01301 if (fDPhi < twopi)
01302 {
01303 sinSPhi = std::sin(fSPhi) ;
01304 cosSPhi = std::cos(fSPhi) ;
01305 ePhi = fSPhi + fDPhi ;
01306 sinEPhi = std::sin(ePhi) ;
01307 cosEPhi = std::cos(ePhi) ;
01308 cPhi = fSPhi + fDPhi*0.5 ;
01309 sinCPhi = std::sin(cPhi) ;
01310 cosCPhi = std::cos(cPhi) ;
01311
01312
01313
01314
01315 vphi = std::atan2(v.y(),v.x()) ;
01316
01317 if ( vphi < fSPhi - halfAngTolerance ) { vphi += twopi; }
01318 else if ( vphi > ePhi + halfAngTolerance ) { vphi -= twopi; }
01319
01320 if ( p.x() || p.y() )
01321 {
01322 pDistS = p.x()*sinSPhi - p.y()*cosSPhi ;
01323 pDistE = -p.x()*sinEPhi + p.y()*cosEPhi ;
01324
01325
01326
01327 compS = -sinSPhi*v.x() + cosSPhi*v.y() ;
01328 compE = sinEPhi*v.x() - cosEPhi*v.y() ;
01329 sidephi = kNull ;
01330
01331 if( ( (fDPhi <= pi) && ( (pDistS <= halfCarTolerance)
01332 && (pDistE <= halfCarTolerance) ) )
01333 || ( (fDPhi > pi) && !((pDistS > halfCarTolerance)
01334 && (pDistE > halfCarTolerance) ) ) )
01335 {
01336
01337
01338 if ( compS < 0 )
01339 {
01340 sphi = pDistS/compS ;
01341
01342 if (sphi >= -halfCarTolerance)
01343 {
01344 xi = p.x() + sphi*v.x() ;
01345 yi = p.y() + sphi*v.y() ;
01346
01347
01348
01349
01350 if ( (std::fabs(xi)<=kCarTolerance)
01351 && (std::fabs(yi)<=kCarTolerance) )
01352 {
01353 sidephi = kSPhi;
01354 if ( ((fSPhi-halfAngTolerance)<=vphi)
01355 && ((ePhi+halfAngTolerance)>=vphi) )
01356 {
01357 sphi = kInfinity;
01358 }
01359 }
01360 else if ( yi*cosCPhi-xi*sinCPhi >=0 )
01361 {
01362 sphi = kInfinity ;
01363 }
01364 else
01365 {
01366 sidephi = kSPhi ;
01367 }
01368 }
01369 else
01370 {
01371 sphi = kInfinity ;
01372 }
01373 }
01374 else
01375 {
01376 sphi = kInfinity ;
01377 }
01378
01379 if ( compE < 0 )
01380 {
01381 sphi2 = pDistE/compE ;
01382
01383
01384
01385 if ( (sphi2 > -kCarTolerance) && (sphi2 < sphi) )
01386 {
01387 xi = p.x() + sphi2*v.x() ;
01388 yi = p.y() + sphi2*v.y() ;
01389
01390 if ( (std::fabs(xi)<=kCarTolerance)
01391 && (std::fabs(yi)<=kCarTolerance) )
01392 {
01393
01394
01395 if( !( (fSPhi-halfAngTolerance <= vphi)
01396 && (ePhi+halfAngTolerance >= vphi) ) )
01397 {
01398 sidephi = kEPhi ;
01399 sphi = sphi2;
01400 }
01401 }
01402 else
01403 {
01404 if ( (yi*cosCPhi-xi*sinCPhi) >= 0)
01405 {
01406
01407
01408 sidephi = kEPhi ;
01409 sphi = sphi2;
01410
01411 }
01412 }
01413 }
01414 }
01415 }
01416 else
01417 {
01418 sphi = kInfinity ;
01419 }
01420 }
01421 else
01422 {
01423
01424
01425
01426 vphi = std::atan2(v.y(),v.x());
01427
01428 if ( ( fSPhi-halfAngTolerance <= vphi ) &&
01429 ( vphi <= ( ePhi+halfAngTolerance ) ) )
01430 {
01431 sphi = kInfinity;
01432 }
01433 else
01434 {
01435 sidephi = kSPhi ;
01436 sphi=0;
01437 }
01438 }
01439
01440
01441
01442 if (sphi<snxt)
01443 {
01444 snxt=sphi;
01445 side=sidephi;
01446 }
01447 }
01448
01449 G4double rhoi2,rhoi,it2,it,iDotxyNmax ;
01450
01451
01452
01453 if (calcNorm)
01454 {
01455 switch(side)
01456 {
01457 case kRMax:
01458 xi = p.x() + snxt*v.x() ;
01459 yi =p.y() + snxt*v.y() ;
01460 zi = p.z() + snxt*v.z() ;
01461 rhoi2 = xi*xi + yi*yi ;
01462 rhoi = std::sqrt(rhoi2) ;
01463 it2 = std::fabs(rhoi2 + zi*zi + fRtor*fRtor - 2*fRtor*rhoi) ;
01464 it = std::sqrt(it2) ;
01465 iDotxyNmax = (1-fRtor/rhoi) ;
01466 if(iDotxyNmax >= -2.*fRmaxTolerance)
01467 {
01468 *n = G4ThreeVector( xi*(1-fRtor/rhoi)/it,
01469 yi*(1-fRtor/rhoi)/it,
01470 zi/it ) ;
01471 *validNorm = true ;
01472 }
01473 else
01474 {
01475 *validNorm = false ;
01476 }
01477 break ;
01478
01479 case kRMin:
01480 *validNorm = false ;
01481 break;
01482
01483 case kSPhi:
01484 if (fDPhi <= pi )
01485 {
01486 *n=G4ThreeVector(std::sin(fSPhi),-std::cos(fSPhi),0);
01487 *validNorm=true;
01488 }
01489 else
01490 {
01491 *validNorm = false ;
01492 }
01493 break ;
01494
01495 case kEPhi:
01496 if (fDPhi <= pi)
01497 {
01498 *n=G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0);
01499 *validNorm=true;
01500 }
01501 else
01502 {
01503 *validNorm = false ;
01504 }
01505 break;
01506
01507 default:
01508
01509
01510
01511 G4cout << G4endl;
01512 DumpInfo();
01513 std::ostringstream message;
01514 G4int oldprc = message.precision(16);
01515 message << "Undefined side for valid surface normal to solid."
01516 << G4endl
01517 << "Position:" << G4endl << G4endl
01518 << "p.x() = " << p.x()/mm << " mm" << G4endl
01519 << "p.y() = " << p.y()/mm << " mm" << G4endl
01520 << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl
01521 << "Direction:" << G4endl << G4endl
01522 << "v.x() = " << v.x() << G4endl
01523 << "v.y() = " << v.y() << G4endl
01524 << "v.z() = " << v.z() << G4endl << G4endl
01525 << "Proposed distance :" << G4endl << G4endl
01526 << "snxt = " << snxt/mm << " mm" << G4endl;
01527 message.precision(oldprc);
01528 G4Exception("G4Torus::DistanceToOut(p,v,..)",
01529 "GeomSolids1002",JustWarning, message);
01530 break;
01531 }
01532 }
01533 if ( snxt<halfCarTolerance ) { snxt=0 ; }
01534
01535 return snxt;
01536 }
01537
01539
01540
01541
01542 G4double G4Torus::DistanceToOut( const G4ThreeVector& p ) const
01543 {
01544 G4double safe=0.0,safeR1,safeR2;
01545 G4double rho2,rho,pt2,pt ;
01546 G4double safePhi,phiC,cosPhiC,sinPhiC,ePhi;
01547 rho2 = p.x()*p.x() + p.y()*p.y() ;
01548 rho = std::sqrt(rho2) ;
01549 pt2 = std::fabs(rho2 + p.z()*p.z() + fRtor*fRtor - 2*fRtor*rho) ;
01550 pt = std::sqrt(pt2) ;
01551
01552 #ifdef G4CSGDEBUG
01553 if( Inside(p) == kOutside )
01554 {
01555 G4int oldprc = G4cout.precision(16) ;
01556 G4cout << G4endl ;
01557 DumpInfo();
01558 G4cout << "Position:" << G4endl << G4endl ;
01559 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ;
01560 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ;
01561 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ;
01562 G4cout.precision(oldprc);
01563 G4Exception("G4Torus::DistanceToOut(p)", "GeomSolids1002",
01564 JustWarning, "Point p is outside !?" );
01565 }
01566 #endif
01567
01568 if (fRmin)
01569 {
01570 safeR1 = pt - fRmin ;
01571 safeR2 = fRmax - pt ;
01572
01573 if (safeR1 < safeR2) { safe = safeR1 ; }
01574 else { safe = safeR2 ; }
01575 }
01576 else
01577 {
01578 safe = fRmax - pt ;
01579 }
01580
01581
01582
01583 if (fDPhi<twopi)
01584 {
01585 phiC = fSPhi + fDPhi*0.5 ;
01586 cosPhiC = std::cos(phiC) ;
01587 sinPhiC = std::sin(phiC) ;
01588
01589 if ((p.y()*cosPhiC-p.x()*sinPhiC)<=0)
01590 {
01591 safePhi = -(p.x()*std::sin(fSPhi) - p.y()*std::cos(fSPhi)) ;
01592 }
01593 else
01594 {
01595 ePhi = fSPhi + fDPhi ;
01596 safePhi = (p.x()*std::sin(ePhi) - p.y()*std::cos(ePhi)) ;
01597 }
01598 if (safePhi < safe) { safe = safePhi ; }
01599 }
01600 if (safe < 0) { safe = 0 ; }
01601 return safe ;
01602 }
01603
01605
01606
01607
01608
01609
01610
01611
01612
01613
01614
01615 G4ThreeVectorList*
01616 G4Torus::CreateRotatedVertices( const G4AffineTransform& pTransform,
01617 G4int& noPolygonVertices ) const
01618 {
01619 G4ThreeVectorList *vertices;
01620 G4ThreeVector vertex0,vertex1,vertex2,vertex3;
01621 G4double meshAngle,meshRMax,crossAngle,cosCrossAngle,sinCrossAngle,sAngle;
01622 G4double rMaxX,rMaxY,rMinX,rMinY;
01623 G4int crossSection,noCrossSections;
01624
01625
01626
01627 noCrossSections = G4int (fDPhi/kMeshAngleDefault) + 1 ;
01628
01629 if (noCrossSections < kMinMeshSections)
01630 {
01631 noCrossSections = kMinMeshSections ;
01632 }
01633 else if (noCrossSections>kMaxMeshSections)
01634 {
01635 noCrossSections=kMaxMeshSections;
01636 }
01637 meshAngle = fDPhi/(noCrossSections - 1) ;
01638 meshRMax = (fRtor + fRmax)/std::cos(meshAngle*0.5) ;
01639
01640
01641
01642
01643 if ( (fDPhi == twopi) && (fSPhi == 0) )
01644 {
01645 sAngle = -meshAngle*0.5 ;
01646 }
01647 else
01648 {
01649 sAngle = fSPhi ;
01650 }
01651 vertices = new G4ThreeVectorList();
01652
01653 if (vertices)
01654 {
01655 vertices->reserve(noCrossSections*4) ;
01656 for (crossSection=0;crossSection<noCrossSections;crossSection++)
01657 {
01658
01659
01660 crossAngle=sAngle+crossSection*meshAngle;
01661 cosCrossAngle=std::cos(crossAngle);
01662 sinCrossAngle=std::sin(crossAngle);
01663
01664 rMaxX=meshRMax*cosCrossAngle;
01665 rMaxY=meshRMax*sinCrossAngle;
01666 rMinX=(fRtor-fRmax)*cosCrossAngle;
01667 rMinY=(fRtor-fRmax)*sinCrossAngle;
01668 vertex0=G4ThreeVector(rMinX,rMinY,-fRmax);
01669 vertex1=G4ThreeVector(rMaxX,rMaxY,-fRmax);
01670 vertex2=G4ThreeVector(rMaxX,rMaxY,+fRmax);
01671 vertex3=G4ThreeVector(rMinX,rMinY,+fRmax);
01672
01673 vertices->push_back(pTransform.TransformPoint(vertex0));
01674 vertices->push_back(pTransform.TransformPoint(vertex1));
01675 vertices->push_back(pTransform.TransformPoint(vertex2));
01676 vertices->push_back(pTransform.TransformPoint(vertex3));
01677 }
01678 noPolygonVertices = 4 ;
01679 }
01680 else
01681 {
01682 DumpInfo();
01683 G4Exception("G4Torus::CreateRotatedVertices()",
01684 "GeomSolids0003", FatalException,
01685 "Error in allocation of vertices. Out of memory !");
01686 }
01687 return vertices;
01688 }
01689
01691
01692
01693
01694 G4GeometryType G4Torus::GetEntityType() const
01695 {
01696 return G4String("G4Torus");
01697 }
01698
01700
01701
01702
01703 G4VSolid* G4Torus::Clone() const
01704 {
01705 return new G4Torus(*this);
01706 }
01707
01709
01710
01711
01712 std::ostream& G4Torus::StreamInfo( std::ostream& os ) const
01713 {
01714 G4int oldprc = os.precision(16);
01715 os << "-----------------------------------------------------------\n"
01716 << " *** Dump for solid - " << GetName() << " ***\n"
01717 << " ===================================================\n"
01718 << " Solid type: G4Torus\n"
01719 << " Parameters: \n"
01720 << " inner radius: " << fRmin/mm << " mm \n"
01721 << " outer radius: " << fRmax/mm << " mm \n"
01722 << " swept radius: " << fRtor/mm << " mm \n"
01723 << " starting phi: " << fSPhi/degree << " degrees \n"
01724 << " delta phi : " << fDPhi/degree << " degrees \n"
01725 << "-----------------------------------------------------------\n";
01726 os.precision(oldprc);
01727
01728 return os;
01729 }
01730
01732
01733
01734
01735 G4ThreeVector G4Torus::GetPointOnSurface() const
01736 {
01737 G4double cosu, sinu,cosv, sinv, aOut, aIn, aSide, chose, phi, theta, rRand;
01738
01739 phi = RandFlat::shoot(fSPhi,fSPhi+fDPhi);
01740 theta = RandFlat::shoot(0.,twopi);
01741
01742 cosu = std::cos(phi); sinu = std::sin(phi);
01743 cosv = std::cos(theta); sinv = std::sin(theta);
01744
01745
01746
01747 aOut = (fDPhi)*twopi*fRtor*fRmax;
01748 aIn = (fDPhi)*twopi*fRtor*fRmin;
01749 aSide = pi*(fRmax*fRmax-fRmin*fRmin);
01750
01751 if ((fSPhi == 0) && (fDPhi == twopi)){ aSide = 0; }
01752 chose = RandFlat::shoot(0.,aOut + aIn + 2.*aSide);
01753
01754 if(chose < aOut)
01755 {
01756 return G4ThreeVector ((fRtor+fRmax*cosv)*cosu,
01757 (fRtor+fRmax*cosv)*sinu, fRmax*sinv);
01758 }
01759 else if( (chose >= aOut) && (chose < aOut + aIn) )
01760 {
01761 return G4ThreeVector ((fRtor+fRmin*cosv)*cosu,
01762 (fRtor+fRmin*cosv)*sinu, fRmin*sinv);
01763 }
01764 else if( (chose >= aOut + aIn) && (chose < aOut + aIn + aSide) )
01765 {
01766 rRand = GetRadiusInRing(fRmin,fRmax);
01767 return G4ThreeVector ((fRtor+rRand*cosv)*std::cos(fSPhi),
01768 (fRtor+rRand*cosv)*std::sin(fSPhi), rRand*sinv);
01769 }
01770 else
01771 {
01772 rRand = GetRadiusInRing(fRmin,fRmax);
01773 return G4ThreeVector ((fRtor+rRand*cosv)*std::cos(fSPhi+fDPhi),
01774 (fRtor+rRand*cosv)*std::sin(fSPhi+fDPhi),
01775 rRand*sinv);
01776 }
01777 }
01778
01780
01781
01782
01783 void G4Torus::DescribeYourselfTo ( G4VGraphicsScene& scene ) const
01784 {
01785 scene.AddSolid (*this);
01786 }
01787
01788 G4Polyhedron* G4Torus::CreatePolyhedron () const
01789 {
01790 return new G4PolyhedronTorus (fRmin, fRmax, fRtor, fSPhi, fDPhi);
01791 }
01792
01793 G4NURBS* G4Torus::CreateNURBS () const
01794 {
01795 G4NURBS* pNURBS;
01796 if (fRmin != 0)
01797 {
01798 if (fDPhi >= twopi)
01799 {
01800 pNURBS = new G4NURBStube(fRmin, fRmax, fRtor);
01801 }
01802 else
01803 {
01804 pNURBS = new G4NURBStubesector(fRmin, fRmax, fRtor, fSPhi, fSPhi + fDPhi);
01805 }
01806 }
01807 else
01808 {
01809 if (fDPhi >= twopi)
01810 {
01811 pNURBS = new G4NURBScylinder (fRmax, fRtor);
01812 }
01813 else
01814 {
01815 const G4double epsilon = 1.e-4;
01816 pNURBS = new G4NURBStubesector (epsilon, fRmax, fRtor,
01817 fSPhi, fSPhi + fDPhi);
01818 }
01819 }
01820 return pNURBS;
01821 }