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