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 "G4Sphere.hh"
00058
00059 #include "G4VoxelLimits.hh"
00060 #include "G4AffineTransform.hh"
00061 #include "G4GeometryTolerance.hh"
00062
00063 #include "G4VPVParameterisation.hh"
00064
00065 #include "Randomize.hh"
00066
00067 #include "meshdefs.hh"
00068
00069 #include "G4VGraphicsScene.hh"
00070 #include "G4VisExtent.hh"
00071 #include "G4Polyhedron.hh"
00072 #include "G4NURBS.hh"
00073 #include "G4NURBSbox.hh"
00074
00075 using namespace CLHEP;
00076
00077
00078
00079 enum ESide {kNull,kRMin,kRMax,kSPhi,kEPhi,kSTheta,kETheta};
00080
00081
00082
00083 enum ENorm {kNRMin,kNRMax,kNSPhi,kNEPhi,kNSTheta,kNETheta};
00084
00086
00087
00088
00089
00090 G4Sphere::G4Sphere( const G4String& pName,
00091 G4double pRmin, G4double pRmax,
00092 G4double pSPhi, G4double pDPhi,
00093 G4double pSTheta, G4double pDTheta )
00094 : G4CSGSolid(pName), fEpsilon(2.e-11),
00095 fFullPhiSphere(true), fFullThetaSphere(true)
00096 {
00097 kAngTolerance = G4GeometryTolerance::GetInstance()->GetAngularTolerance();
00098
00099
00100
00101 kRadTolerance = G4GeometryTolerance::GetInstance()->GetRadialTolerance();
00102 if ( (pRmin >= pRmax) || (pRmax < 1.1*kRadTolerance) || (pRmin < 0) )
00103 {
00104 std::ostringstream message;
00105 message << "Invalid radii for Solid: " << GetName() << G4endl
00106 << " pRmin = " << pRmin << ", pRmax = " << pRmax;
00107 G4Exception("G4Sphere::G4Sphere()", "GeomSolids0002",
00108 FatalException, message);
00109 }
00110 fRmin=pRmin; fRmax=pRmax;
00111 fRminTolerance = (fRmin) ? std::max( kRadTolerance, fEpsilon*fRmin ) : 0;
00112 fRmaxTolerance = std::max( kRadTolerance, fEpsilon*fRmax );
00113
00114
00115
00116 CheckPhiAngles(pSPhi, pDPhi);
00117 CheckThetaAngles(pSTheta, pDTheta);
00118 }
00119
00121
00122
00123
00124
00125 G4Sphere::G4Sphere( __void__& a )
00126 : G4CSGSolid(a), fRminTolerance(0.), fRmaxTolerance(0.),
00127 kAngTolerance(0.), kRadTolerance(0.), fEpsilon(0.),
00128 fRmin(0.), fRmax(0.), fSPhi(0.), fDPhi(0.), fSTheta(0.),
00129 fDTheta(0.), sinCPhi(0.), cosCPhi(0.), cosHDPhiOT(0.), cosHDPhiIT(0.),
00130 sinSPhi(0.), cosSPhi(0.), sinEPhi(0.), cosEPhi(0.), hDPhi(0.), cPhi(0.),
00131 ePhi(0.), sinSTheta(0.), cosSTheta(0.), sinETheta(0.), cosETheta(0.),
00132 tanSTheta(0.), tanSTheta2(0.), tanETheta(0.), tanETheta2(0.), eTheta(0.),
00133 fFullPhiSphere(false), fFullThetaSphere(false), fFullSphere(true)
00134 {
00135 }
00136
00138
00139
00140
00141 G4Sphere::~G4Sphere()
00142 {
00143 }
00144
00146
00147
00148
00149 G4Sphere::G4Sphere(const G4Sphere& rhs)
00150 : G4CSGSolid(rhs), fRminTolerance(rhs.fRminTolerance),
00151 fRmaxTolerance(rhs.fRmaxTolerance), kAngTolerance(rhs.kAngTolerance),
00152 kRadTolerance(rhs.kRadTolerance), fEpsilon(rhs.fEpsilon),
00153 fRmin(rhs.fRmin), fRmax(rhs.fRmax), fSPhi(rhs.fSPhi), fDPhi(rhs.fDPhi),
00154 fSTheta(rhs.fSTheta), fDTheta(rhs.fDTheta),
00155 sinCPhi(rhs.sinCPhi), cosCPhi(rhs.cosCPhi),
00156 cosHDPhiOT(rhs.cosHDPhiOT), cosHDPhiIT(rhs.cosHDPhiIT),
00157 sinSPhi(rhs.sinSPhi), cosSPhi(rhs.cosSPhi),
00158 sinEPhi(rhs.sinEPhi), cosEPhi(rhs.cosEPhi),
00159 hDPhi(rhs.hDPhi), cPhi(rhs.cPhi), ePhi(rhs.ePhi),
00160 sinSTheta(rhs.sinSTheta), cosSTheta(rhs.cosSTheta),
00161 sinETheta(rhs.sinETheta), cosETheta(rhs.cosETheta),
00162 tanSTheta(rhs.tanSTheta), tanSTheta2(rhs.tanSTheta2),
00163 tanETheta(rhs.tanETheta), tanETheta2(rhs.tanETheta2), eTheta(rhs.eTheta),
00164 fFullPhiSphere(rhs.fFullPhiSphere), fFullThetaSphere(rhs.fFullThetaSphere),
00165 fFullSphere(rhs.fFullSphere)
00166 {
00167 }
00168
00170
00171
00172
00173 G4Sphere& G4Sphere::operator = (const G4Sphere& rhs)
00174 {
00175
00176
00177 if (this == &rhs) { return *this; }
00178
00179
00180
00181 G4CSGSolid::operator=(rhs);
00182
00183
00184
00185 fRminTolerance = rhs.fRminTolerance; fRmaxTolerance = rhs.fRmaxTolerance;
00186 kAngTolerance = rhs.kAngTolerance; kRadTolerance = rhs.kRadTolerance;
00187 fEpsilon = rhs.fEpsilon; fRmin = rhs.fRmin; fRmax = rhs.fRmax;
00188 fSPhi = rhs.fSPhi; fDPhi = rhs.fDPhi; fSTheta = rhs.fSTheta;
00189 fDTheta = rhs.fDTheta; sinCPhi = rhs.sinCPhi; cosCPhi = rhs.cosCPhi;
00190 cosHDPhiOT = rhs.cosHDPhiOT; cosHDPhiIT = rhs.cosHDPhiIT;
00191 sinSPhi = rhs.sinSPhi; cosSPhi = rhs.cosSPhi;
00192 sinEPhi = rhs.sinEPhi; cosEPhi = rhs.cosEPhi;
00193 hDPhi = rhs.hDPhi; cPhi = rhs.cPhi; ePhi = rhs.ePhi;
00194 sinSTheta = rhs.sinSTheta; cosSTheta = rhs.cosSTheta;
00195 sinETheta = rhs.sinETheta; cosETheta = rhs.cosETheta;
00196 tanSTheta = rhs.tanSTheta; tanSTheta2 = rhs.tanSTheta2;
00197 tanETheta = rhs.tanETheta; tanETheta2 = rhs.tanETheta2;
00198 eTheta = rhs.eTheta; fFullPhiSphere = rhs.fFullPhiSphere;
00199 fFullThetaSphere = rhs.fFullThetaSphere; fFullSphere = rhs.fFullSphere;
00200
00201 return *this;
00202 }
00203
00205
00206
00207
00208
00209 void G4Sphere::ComputeDimensions( G4VPVParameterisation* p,
00210 const G4int n,
00211 const G4VPhysicalVolume* pRep)
00212 {
00213 p->ComputeDimensions(*this,n,pRep);
00214 }
00215
00217
00218
00219
00220 G4bool G4Sphere::CalculateExtent( const EAxis pAxis,
00221 const G4VoxelLimits& pVoxelLimit,
00222 const G4AffineTransform& pTransform,
00223 G4double& pMin, G4double& pMax ) const
00224 {
00225 if ( fFullSphere )
00226 {
00227
00228
00229
00230
00231
00232
00233 G4double xoffset,xMin,xMax;
00234 G4double yoffset,yMin,yMax;
00235 G4double zoffset,zMin,zMax;
00236
00237 G4double diff1,diff2,delta,maxDiff,newMin,newMax;
00238 G4double xoff1,xoff2,yoff1,yoff2;
00239
00240 xoffset=pTransform.NetTranslation().x();
00241 xMin=xoffset-fRmax;
00242 xMax=xoffset+fRmax;
00243 if (pVoxelLimit.IsXLimited())
00244 {
00245 if ( (xMin>pVoxelLimit.GetMaxXExtent()+kCarTolerance)
00246 || (xMax<pVoxelLimit.GetMinXExtent()-kCarTolerance) )
00247 {
00248 return false;
00249 }
00250 else
00251 {
00252 if (xMin<pVoxelLimit.GetMinXExtent())
00253 {
00254 xMin=pVoxelLimit.GetMinXExtent();
00255 }
00256 if (xMax>pVoxelLimit.GetMaxXExtent())
00257 {
00258 xMax=pVoxelLimit.GetMaxXExtent();
00259 }
00260 }
00261 }
00262
00263 yoffset=pTransform.NetTranslation().y();
00264 yMin=yoffset-fRmax;
00265 yMax=yoffset+fRmax;
00266 if (pVoxelLimit.IsYLimited())
00267 {
00268 if ( (yMin>pVoxelLimit.GetMaxYExtent()+kCarTolerance)
00269 || (yMax<pVoxelLimit.GetMinYExtent()-kCarTolerance) )
00270 {
00271 return false;
00272 }
00273 else
00274 {
00275 if (yMin<pVoxelLimit.GetMinYExtent())
00276 {
00277 yMin=pVoxelLimit.GetMinYExtent();
00278 }
00279 if (yMax>pVoxelLimit.GetMaxYExtent())
00280 {
00281 yMax=pVoxelLimit.GetMaxYExtent();
00282 }
00283 }
00284 }
00285
00286 zoffset=pTransform.NetTranslation().z();
00287 zMin=zoffset-fRmax;
00288 zMax=zoffset+fRmax;
00289 if (pVoxelLimit.IsZLimited())
00290 {
00291 if ( (zMin>pVoxelLimit.GetMaxZExtent()+kCarTolerance)
00292 || (zMax<pVoxelLimit.GetMinZExtent()-kCarTolerance) )
00293 {
00294 return false;
00295 }
00296 else
00297 {
00298 if (zMin<pVoxelLimit.GetMinZExtent())
00299 {
00300 zMin=pVoxelLimit.GetMinZExtent();
00301 }
00302 if (zMax>pVoxelLimit.GetMaxZExtent())
00303 {
00304 zMax=pVoxelLimit.GetMaxZExtent();
00305 }
00306 }
00307 }
00308
00309
00310
00311 switch (pAxis)
00312 {
00313 case kXAxis:
00314 yoff1=yoffset-yMin;
00315 yoff2=yMax-yoffset;
00316 if ((yoff1>=0) && (yoff2>=0))
00317 {
00318
00319
00320 pMin=xMin;
00321 pMax=xMax;
00322 }
00323 else
00324 {
00325
00326
00327
00328 delta=fRmax*fRmax-yoff1*yoff1;
00329 diff1=(delta>0.) ? std::sqrt(delta) : 0.;
00330 delta=fRmax*fRmax-yoff2*yoff2;
00331 diff2=(delta>0.) ? std::sqrt(delta) : 0.;
00332 maxDiff=(diff1>diff2) ? diff1:diff2;
00333 newMin=xoffset-maxDiff;
00334 newMax=xoffset+maxDiff;
00335 pMin=(newMin<xMin) ? xMin : newMin;
00336 pMax=(newMax>xMax) ? xMax : newMax;
00337 }
00338 break;
00339 case kYAxis:
00340 xoff1=xoffset-xMin;
00341 xoff2=xMax-xoffset;
00342 if ((xoff1>=0) && (xoff2>=0))
00343 {
00344
00345
00346 pMin=yMin;
00347 pMax=yMax;
00348 }
00349 else
00350 {
00351
00352
00353
00354 delta=fRmax*fRmax-xoff1*xoff1;
00355 diff1=(delta>0.) ? std::sqrt(delta) : 0.;
00356 delta=fRmax*fRmax-xoff2*xoff2;
00357 diff2=(delta>0.) ? std::sqrt(delta) : 0.;
00358 maxDiff=(diff1>diff2) ? diff1:diff2;
00359 newMin=yoffset-maxDiff;
00360 newMax=yoffset+maxDiff;
00361 pMin=(newMin<yMin) ? yMin : newMin;
00362 pMax=(newMax>yMax) ? yMax : newMax;
00363 }
00364 break;
00365 case kZAxis:
00366 pMin=zMin;
00367 pMax=zMax;
00368 break;
00369 default:
00370 break;
00371 }
00372 pMin-=kCarTolerance;
00373 pMax+=kCarTolerance;
00374
00375 return true;
00376 }
00377 else
00378 {
00379 G4int i,j,noEntries,noBetweenSections;
00380 G4bool existsAfterClip=false;
00381
00382
00383
00384 G4ThreeVectorList* vertices;
00385 G4int noPolygonVertices ;
00386 vertices=CreateRotatedVertices(pTransform,noPolygonVertices);
00387
00388 pMin=+kInfinity;
00389 pMax=-kInfinity;
00390
00391 noEntries=vertices->size();
00392 noBetweenSections=noEntries-noPolygonVertices;
00393
00394 G4ThreeVectorList ThetaPolygon ;
00395 for (i=0;i<noEntries;i+=noPolygonVertices)
00396 {
00397 for(j=0;j<(noPolygonVertices/2)-1;j++)
00398 {
00399 ThetaPolygon.push_back((*vertices)[i+j]) ;
00400 ThetaPolygon.push_back((*vertices)[i+j+1]) ;
00401 ThetaPolygon.push_back((*vertices)[i+noPolygonVertices-2-j]) ;
00402 ThetaPolygon.push_back((*vertices)[i+noPolygonVertices-1-j]) ;
00403 CalculateClippedPolygonExtent(ThetaPolygon,pVoxelLimit,pAxis,pMin,pMax);
00404 ThetaPolygon.clear() ;
00405 }
00406 }
00407 for (i=0;i<noBetweenSections;i+=noPolygonVertices)
00408 {
00409 for(j=0;j<noPolygonVertices-1;j++)
00410 {
00411 ThetaPolygon.push_back((*vertices)[i+j]) ;
00412 ThetaPolygon.push_back((*vertices)[i+j+1]) ;
00413 ThetaPolygon.push_back((*vertices)[i+noPolygonVertices+j+1]) ;
00414 ThetaPolygon.push_back((*vertices)[i+noPolygonVertices+j]) ;
00415 CalculateClippedPolygonExtent(ThetaPolygon,pVoxelLimit,pAxis,pMin,pMax);
00416 ThetaPolygon.clear() ;
00417 }
00418 ThetaPolygon.push_back((*vertices)[i+noPolygonVertices-1]) ;
00419 ThetaPolygon.push_back((*vertices)[i]) ;
00420 ThetaPolygon.push_back((*vertices)[i+noPolygonVertices]) ;
00421 ThetaPolygon.push_back((*vertices)[i+2*noPolygonVertices-1]) ;
00422 CalculateClippedPolygonExtent(ThetaPolygon,pVoxelLimit,pAxis,pMin,pMax);
00423 ThetaPolygon.clear() ;
00424 }
00425
00426 if ((pMin!=kInfinity) || (pMax!=-kInfinity))
00427 {
00428 existsAfterClip=true;
00429
00430
00431
00432 pMin-=kCarTolerance;
00433 pMax+=kCarTolerance;
00434 }
00435 else
00436 {
00437
00438
00439
00440
00441
00442 G4ThreeVector clipCentre(
00443 (pVoxelLimit.GetMinXExtent()+pVoxelLimit.GetMaxXExtent())*0.5,
00444 (pVoxelLimit.GetMinYExtent()+pVoxelLimit.GetMaxYExtent())*0.5,
00445 (pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5);
00446
00447 if (Inside(pTransform.Inverse().TransformPoint(clipCentre))!=kOutside)
00448 {
00449 existsAfterClip=true;
00450 pMin=pVoxelLimit.GetMinExtent(pAxis);
00451 pMax=pVoxelLimit.GetMaxExtent(pAxis);
00452 }
00453 }
00454 delete vertices;
00455 return existsAfterClip;
00456 }
00457 }
00458
00460
00461
00462
00463
00464
00465 EInside G4Sphere::Inside( const G4ThreeVector& p ) const
00466 {
00467 G4double rho,rho2,rad2,tolRMin,tolRMax;
00468 G4double pPhi,pTheta;
00469 EInside in = kOutside;
00470 static const G4double halfAngTolerance = kAngTolerance*0.5;
00471 const G4double halfRmaxTolerance = fRmaxTolerance*0.5;
00472 const G4double halfRminTolerance = fRminTolerance*0.5;
00473 const G4double Rmax_minus = fRmax - halfRmaxTolerance;
00474 const G4double Rmin_plus = (fRmin > 0) ? fRmin+halfRminTolerance : 0;
00475
00476 rho2 = p.x()*p.x() + p.y()*p.y() ;
00477 rad2 = rho2 + p.z()*p.z() ;
00478
00479
00480
00481 tolRMin = Rmin_plus;
00482 tolRMax = Rmax_minus;
00483
00484 if ( (rad2 <= Rmax_minus*Rmax_minus) && (rad2 >= Rmin_plus*Rmin_plus) )
00485 {
00486 in = kInside;
00487 }
00488 else
00489 {
00490 tolRMax = fRmax + halfRmaxTolerance;
00491 tolRMin = std::max(fRmin-halfRminTolerance, 0.);
00492 if ( (rad2 <= tolRMax*tolRMax) && (rad2 >= tolRMin*tolRMin) )
00493 {
00494 in = kSurface;
00495 }
00496 else
00497 {
00498 return in = kOutside;
00499 }
00500 }
00501
00502
00503
00504 if ( !fFullPhiSphere && rho2 )
00505 {
00506 pPhi = std::atan2(p.y(),p.x()) ;
00507
00508 if ( pPhi < fSPhi - halfAngTolerance ) { pPhi += twopi; }
00509 else if ( pPhi > ePhi + halfAngTolerance ) { pPhi -= twopi; }
00510
00511 if ( (pPhi < fSPhi - halfAngTolerance)
00512 || (pPhi > ePhi + halfAngTolerance) ) { return in = kOutside; }
00513
00514 else if (in == kInside)
00515 {
00516 if ( (pPhi < fSPhi + halfAngTolerance)
00517 || (pPhi > ePhi - halfAngTolerance) ) { in = kSurface; }
00518 }
00519 }
00520
00521
00522
00523 if ( (rho2 || p.z()) && (!fFullThetaSphere) )
00524 {
00525 rho = std::sqrt(rho2);
00526 pTheta = std::atan2(rho,p.z());
00527
00528 if ( in == kInside )
00529 {
00530 if ( (pTheta < fSTheta + halfAngTolerance)
00531 || (pTheta > eTheta - halfAngTolerance) )
00532 {
00533 if ( (pTheta >= fSTheta - halfAngTolerance)
00534 && (pTheta <= eTheta + halfAngTolerance) )
00535 {
00536 in = kSurface;
00537 }
00538 else
00539 {
00540 in = kOutside;
00541 }
00542 }
00543 }
00544 else
00545 {
00546 if ( (pTheta < fSTheta - halfAngTolerance)
00547 || (pTheta > eTheta + halfAngTolerance) )
00548 {
00549 in = kOutside;
00550 }
00551 }
00552 }
00553 return in;
00554 }
00555
00557
00558
00559
00560
00561
00562 G4ThreeVector G4Sphere::SurfaceNormal( const G4ThreeVector& p ) const
00563 {
00564 G4int noSurfaces = 0;
00565 G4double rho, rho2, radius, pTheta, pPhi=0.;
00566 G4double distRMin = kInfinity;
00567 G4double distSPhi = kInfinity, distEPhi = kInfinity;
00568 G4double distSTheta = kInfinity, distETheta = kInfinity;
00569 G4ThreeVector nR, nPs, nPe, nTs, nTe, nZ(0.,0.,1.);
00570 G4ThreeVector norm, sumnorm(0.,0.,0.);
00571
00572 static const G4double halfCarTolerance = 0.5*kCarTolerance;
00573 static const G4double halfAngTolerance = 0.5*kAngTolerance;
00574
00575 rho2 = p.x()*p.x()+p.y()*p.y();
00576 radius = std::sqrt(rho2+p.z()*p.z());
00577 rho = std::sqrt(rho2);
00578
00579 G4double distRMax = std::fabs(radius-fRmax);
00580 if (fRmin) distRMin = std::fabs(radius-fRmin);
00581
00582 if ( rho && !fFullSphere )
00583 {
00584 pPhi = std::atan2(p.y(),p.x());
00585
00586 if (pPhi < fSPhi-halfAngTolerance) { pPhi += twopi; }
00587 else if (pPhi > ePhi+halfAngTolerance) { pPhi -= twopi; }
00588 }
00589 if ( !fFullPhiSphere )
00590 {
00591 if ( rho )
00592 {
00593 distSPhi = std::fabs( pPhi-fSPhi );
00594 distEPhi = std::fabs( pPhi-ePhi );
00595 }
00596 else if( !fRmin )
00597 {
00598 distSPhi = 0.;
00599 distEPhi = 0.;
00600 }
00601 nPs = G4ThreeVector(sinSPhi,-cosSPhi,0);
00602 nPe = G4ThreeVector(-sinEPhi,cosEPhi,0);
00603 }
00604 if ( !fFullThetaSphere )
00605 {
00606 if ( rho )
00607 {
00608 pTheta = std::atan2(rho,p.z());
00609 distSTheta = std::fabs(pTheta-fSTheta);
00610 distETheta = std::fabs(pTheta-eTheta);
00611
00612 nTs = G4ThreeVector(-cosSTheta*p.x()/rho,
00613 -cosSTheta*p.y()/rho,
00614 sinSTheta );
00615
00616 nTe = G4ThreeVector( cosETheta*p.x()/rho,
00617 cosETheta*p.y()/rho,
00618 -sinETheta );
00619 }
00620 else if( !fRmin )
00621 {
00622 if ( fSTheta )
00623 {
00624 distSTheta = 0.;
00625 nTs = G4ThreeVector(0.,0.,-1.);
00626 }
00627 if ( eTheta < pi )
00628 {
00629 distETheta = 0.;
00630 nTe = G4ThreeVector(0.,0.,1.);
00631 }
00632 }
00633 }
00634 if( radius ) { nR = G4ThreeVector(p.x()/radius,p.y()/radius,p.z()/radius); }
00635
00636 if( distRMax <= halfCarTolerance )
00637 {
00638 noSurfaces ++;
00639 sumnorm += nR;
00640 }
00641 if( fRmin && (distRMin <= halfCarTolerance) )
00642 {
00643 noSurfaces ++;
00644 sumnorm -= nR;
00645 }
00646 if( !fFullPhiSphere )
00647 {
00648 if (distSPhi <= halfAngTolerance)
00649 {
00650 noSurfaces ++;
00651 sumnorm += nPs;
00652 }
00653 if (distEPhi <= halfAngTolerance)
00654 {
00655 noSurfaces ++;
00656 sumnorm += nPe;
00657 }
00658 }
00659 if ( !fFullThetaSphere )
00660 {
00661 if ((distSTheta <= halfAngTolerance) && (fSTheta > 0.))
00662 {
00663 noSurfaces ++;
00664 if ((radius <= halfCarTolerance) && fFullPhiSphere) { sumnorm += nZ; }
00665 else { sumnorm += nTs; }
00666 }
00667 if ((distETheta <= halfAngTolerance) && (eTheta < pi))
00668 {
00669 noSurfaces ++;
00670 if ((radius <= halfCarTolerance) && fFullPhiSphere) { sumnorm -= nZ; }
00671 else { sumnorm += nTe; }
00672 if(sumnorm.z() == 0.) { sumnorm += nZ; }
00673 }
00674 }
00675 if ( noSurfaces == 0 )
00676 {
00677 #ifdef G4CSGDEBUG
00678 G4Exception("G4Sphere::SurfaceNormal(p)", "GeomSolids1002",
00679 JustWarning, "Point p is not on surface !?" );
00680 #endif
00681 norm = ApproxSurfaceNormal(p);
00682 }
00683 else if ( noSurfaces == 1 ) { norm = sumnorm; }
00684 else { norm = sumnorm.unit(); }
00685 return norm;
00686 }
00687
00688
00690
00691
00692
00693
00694 G4ThreeVector G4Sphere::ApproxSurfaceNormal( const G4ThreeVector& p ) const
00695 {
00696 ENorm side;
00697 G4ThreeVector norm;
00698 G4double rho,rho2,radius,pPhi,pTheta;
00699 G4double distRMin,distRMax,distSPhi,distEPhi,
00700 distSTheta,distETheta,distMin;
00701
00702 rho2=p.x()*p.x()+p.y()*p.y();
00703 radius=std::sqrt(rho2+p.z()*p.z());
00704 rho=std::sqrt(rho2);
00705
00706
00707
00708
00709
00710 distRMax=std::fabs(radius-fRmax);
00711 if (fRmin)
00712 {
00713 distRMin=std::fabs(radius-fRmin);
00714
00715 if (distRMin<distRMax)
00716 {
00717 distMin=distRMin;
00718 side=kNRMin;
00719 }
00720 else
00721 {
00722 distMin=distRMax;
00723 side=kNRMax;
00724 }
00725 }
00726 else
00727 {
00728 distMin=distRMax;
00729 side=kNRMax;
00730 }
00731
00732
00733
00734
00735
00736
00737 pPhi = std::atan2(p.y(),p.x());
00738 if (pPhi<0) { pPhi += twopi; }
00739
00740 if (!fFullPhiSphere && rho)
00741 {
00742 if (fSPhi<0)
00743 {
00744 distSPhi=std::fabs(pPhi-(fSPhi+twopi))*rho;
00745 }
00746 else
00747 {
00748 distSPhi=std::fabs(pPhi-fSPhi)*rho;
00749 }
00750
00751 distEPhi=std::fabs(pPhi-fSPhi-fDPhi)*rho;
00752
00753
00754
00755 if (distSPhi<distEPhi)
00756 {
00757 if (distSPhi<distMin)
00758 {
00759 distMin=distSPhi;
00760 side=kNSPhi;
00761 }
00762 }
00763 else
00764 {
00765 if (distEPhi<distMin)
00766 {
00767 distMin=distEPhi;
00768 side=kNEPhi;
00769 }
00770 }
00771 }
00772
00773
00774
00775
00776
00777 if (!fFullThetaSphere && radius)
00778 {
00779 pTheta=std::atan2(rho,p.z());
00780 distSTheta=std::fabs(pTheta-fSTheta)*radius;
00781 distETheta=std::fabs(pTheta-fSTheta-fDTheta)*radius;
00782
00783
00784
00785 if (distSTheta<distETheta)
00786 {
00787 if (distSTheta<distMin)
00788 {
00789 distMin = distSTheta ;
00790 side = kNSTheta ;
00791 }
00792 }
00793 else
00794 {
00795 if (distETheta<distMin)
00796 {
00797 distMin = distETheta ;
00798 side = kNETheta ;
00799 }
00800 }
00801 }
00802
00803 switch (side)
00804 {
00805 case kNRMin:
00806 norm=G4ThreeVector(-p.x()/radius,-p.y()/radius,-p.z()/radius);
00807 break;
00808 case kNRMax:
00809 norm=G4ThreeVector(p.x()/radius,p.y()/radius,p.z()/radius);
00810 break;
00811 case kNSPhi:
00812 norm=G4ThreeVector(sinSPhi,-cosSPhi,0);
00813 break;
00814 case kNEPhi:
00815 norm=G4ThreeVector(-sinEPhi,cosEPhi,0);
00816 break;
00817 case kNSTheta:
00818 norm=G4ThreeVector(-cosSTheta*std::cos(pPhi),
00819 -cosSTheta*std::sin(pPhi),
00820 sinSTheta );
00821 break;
00822 case kNETheta:
00823 norm=G4ThreeVector( cosETheta*std::cos(pPhi),
00824 cosETheta*std::sin(pPhi),
00825 -sinETheta );
00826 break;
00827 default:
00828 DumpInfo();
00829 G4Exception("G4Sphere::ApproxSurfaceNormal()",
00830 "GeomSolids1002", JustWarning,
00831 "Undefined side for valid surface normal to solid.");
00832 break;
00833 }
00834
00835 return norm;
00836 }
00837
00839
00840
00841
00842
00843
00844
00845
00846
00847
00848
00849
00850
00851
00852
00853
00854
00855
00856
00857
00858
00859
00860
00861
00862
00863
00864
00865
00866
00867 G4double G4Sphere::DistanceToIn( const G4ThreeVector& p,
00868 const G4ThreeVector& v ) const
00869 {
00870 G4double snxt = kInfinity ;
00871 G4double rho2, rad2, pDotV2d, pDotV3d, pTheta ;
00872 G4double tolSTheta=0., tolETheta=0. ;
00873 const G4double dRmax = 100.*fRmax;
00874
00875 static const G4double halfCarTolerance = kCarTolerance*0.5;
00876 static const G4double halfAngTolerance = kAngTolerance*0.5;
00877 const G4double halfRmaxTolerance = fRmaxTolerance*0.5;
00878 const G4double halfRminTolerance = fRminTolerance*0.5;
00879 const G4double tolORMin2 = (fRmin>halfRminTolerance)
00880 ? (fRmin-halfRminTolerance)*(fRmin-halfRminTolerance) : 0;
00881 const G4double tolIRMin2 =
00882 (fRmin+halfRminTolerance)*(fRmin+halfRminTolerance);
00883 const G4double tolORMax2 =
00884 (fRmax+halfRmaxTolerance)*(fRmax+halfRmaxTolerance);
00885 const G4double tolIRMax2 =
00886 (fRmax-halfRmaxTolerance)*(fRmax-halfRmaxTolerance);
00887
00888
00889
00890 G4double xi, yi, zi, rhoi, rhoi2, radi2, iTheta ;
00891
00892
00893
00894 G4double Comp ;
00895
00896
00897
00898 G4double Dist, cosPsi ;
00899
00900
00901
00902 G4double dist2STheta, dist2ETheta ;
00903 G4double t1, t2, b, c, d2, d, sd = kInfinity ;
00904
00905
00906
00907 rho2 = p.x()*p.x() + p.y()*p.y() ;
00908 rad2 = rho2 + p.z()*p.z() ;
00909 pTheta = std::atan2(std::sqrt(rho2),p.z()) ;
00910
00911 pDotV2d = p.x()*v.x() + p.y()*v.y() ;
00912 pDotV3d = pDotV2d + p.z()*v.z() ;
00913
00914
00915
00916 if (!fFullThetaSphere)
00917 {
00918 tolSTheta = fSTheta - halfAngTolerance ;
00919 tolETheta = eTheta + halfAngTolerance ;
00920 }
00921
00922
00923
00924
00925
00926
00927
00928
00929
00930
00931
00932
00933
00934
00935
00936 c = rad2 - fRmax*fRmax ;
00937
00938 if (c > fRmaxTolerance*fRmax)
00939 {
00940
00941
00942
00943 d2 = pDotV3d*pDotV3d - c ;
00944
00945 if ( d2 >= 0 )
00946 {
00947 sd = -pDotV3d - std::sqrt(d2) ;
00948
00949 if (sd >= 0 )
00950 {
00951 if ( sd>dRmax )
00952 {
00953 G4double fTerm = sd-std::fmod(sd,dRmax);
00954 sd = fTerm + DistanceToIn(p+fTerm*v,v);
00955 }
00956 xi = p.x() + sd*v.x() ;
00957 yi = p.y() + sd*v.y() ;
00958 rhoi = std::sqrt(xi*xi + yi*yi) ;
00959
00960 if (!fFullPhiSphere && rhoi)
00961 {
00962 cosPsi = (xi*cosCPhi + yi*sinCPhi)/rhoi ;
00963
00964 if (cosPsi >= cosHDPhiOT)
00965 {
00966 if (!fFullThetaSphere)
00967 {
00968 zi = p.z() + sd*v.z() ;
00969
00970
00971
00972
00973 iTheta = std::atan2(rhoi,zi) ;
00974 if ( (iTheta >= tolSTheta) && (iTheta <= tolETheta) )
00975 {
00976 return snxt = sd ;
00977 }
00978 }
00979 else
00980 {
00981 return snxt=sd;
00982 }
00983 }
00984 }
00985 else
00986 {
00987 if (!fFullThetaSphere)
00988 {
00989 zi = p.z() + sd*v.z() ;
00990
00991
00992
00993
00994 iTheta = std::atan2(rhoi,zi) ;
00995 if ( (iTheta >= tolSTheta) && (iTheta <= tolETheta) )
00996 {
00997 return snxt=sd;
00998 }
00999 }
01000 else
01001 {
01002 return snxt = sd;
01003 }
01004 }
01005 }
01006 }
01007 else
01008 {
01009 return snxt=kInfinity;
01010 }
01011 }
01012 else
01013 {
01014
01015
01016
01017 d2 = pDotV3d*pDotV3d - c ;
01018
01019 if ( (rad2 > tolIRMax2)
01020 && ( (d2 >= fRmaxTolerance*fRmax) && (pDotV3d < 0) ) )
01021 {
01022 if (!fFullPhiSphere)
01023 {
01024
01025
01026
01027 cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(rho2) ;
01028
01029 if (cosPsi>=cosHDPhiIT)
01030 {
01031
01032
01033 if ( !fFullThetaSphere )
01034 {
01035 if ( (pTheta >= tolSTheta + kAngTolerance)
01036 && (pTheta <= tolETheta - kAngTolerance) )
01037 {
01038 return snxt=0;
01039 }
01040 }
01041 else
01042 {
01043 return snxt=0;
01044 }
01045 }
01046 }
01047 else
01048 {
01049 if ( !fFullThetaSphere )
01050 {
01051 if ( (pTheta >= tolSTheta + kAngTolerance)
01052 && (pTheta <= tolETheta - kAngTolerance) )
01053 {
01054 return snxt=0;
01055 }
01056 }
01057 else
01058 {
01059 return snxt=0;
01060 }
01061 }
01062 }
01063 }
01064
01065
01066
01067
01068
01069
01070 if (fRmin)
01071 {
01072 c = rad2 - fRmin*fRmin ;
01073 d2 = pDotV3d*pDotV3d - c ;
01074
01075
01076
01077
01078 if ( (c > -halfRminTolerance) && (rad2 < tolIRMin2)
01079 && ( (d2 < fRmin*kCarTolerance) || (pDotV3d >= 0) ) )
01080 {
01081 if ( !fFullPhiSphere )
01082 {
01083
01084
01085
01086 cosPsi = (p.x()*cosCPhi+p.y()*sinCPhi)/std::sqrt(rho2) ;
01087 if (cosPsi >= cosHDPhiIT)
01088 {
01089
01090
01091 if ( !fFullThetaSphere )
01092 {
01093 if ( (pTheta >= tolSTheta + kAngTolerance)
01094 && (pTheta <= tolETheta - kAngTolerance) )
01095 {
01096 return snxt=0;
01097 }
01098 }
01099 else
01100 {
01101 return snxt = 0 ;
01102 }
01103 }
01104 }
01105 else
01106 {
01107 if ( !fFullThetaSphere )
01108 {
01109 if ( (pTheta >= tolSTheta + kAngTolerance)
01110 && (pTheta <= tolETheta - kAngTolerance) )
01111 {
01112 return snxt = 0 ;
01113 }
01114 }
01115 else
01116 {
01117 return snxt=0;
01118 }
01119 }
01120 }
01121 else
01122 {
01123 if (d2 >= 0)
01124 {
01125 sd = -pDotV3d + std::sqrt(d2) ;
01126 if ( sd >= halfRminTolerance )
01127 {
01128 xi = p.x() + sd*v.x() ;
01129 yi = p.y() + sd*v.y() ;
01130 rhoi = std::sqrt(xi*xi+yi*yi) ;
01131
01132 if ( !fFullPhiSphere && rhoi )
01133 {
01134 cosPsi = (xi*cosCPhi + yi*sinCPhi)/rhoi ;
01135
01136 if (cosPsi >= cosHDPhiOT)
01137 {
01138 if ( !fFullThetaSphere )
01139 {
01140 zi = p.z() + sd*v.z() ;
01141
01142
01143
01144
01145 iTheta = std::atan2(rhoi,zi) ;
01146 if ( (iTheta >= tolSTheta) && (iTheta<=tolETheta) )
01147 {
01148 snxt = sd;
01149 }
01150 }
01151 else
01152 {
01153 snxt=sd;
01154 }
01155 }
01156 }
01157 else
01158 {
01159 if ( !fFullThetaSphere )
01160 {
01161 zi = p.z() + sd*v.z() ;
01162
01163
01164
01165
01166 iTheta = std::atan2(rhoi,zi) ;
01167 if ( (iTheta >= tolSTheta) && (iTheta <= tolETheta) )
01168 {
01169 snxt = sd;
01170 }
01171 }
01172 else
01173 {
01174 snxt = sd;
01175 }
01176 }
01177 }
01178 }
01179 }
01180 }
01181
01182
01183
01184
01185
01186
01187
01188
01189
01190
01191 if ( !fFullPhiSphere )
01192 {
01193
01194
01195
01196 Comp = v.x()*sinSPhi - v.y()*cosSPhi ;
01197
01198 if ( Comp < 0 )
01199 {
01200 Dist = p.y()*cosSPhi - p.x()*sinSPhi ;
01201
01202 if (Dist < halfCarTolerance)
01203 {
01204 sd = Dist/Comp ;
01205
01206 if (sd < snxt)
01207 {
01208 if ( sd > 0 )
01209 {
01210 xi = p.x() + sd*v.x() ;
01211 yi = p.y() + sd*v.y() ;
01212 zi = p.z() + sd*v.z() ;
01213 rhoi2 = xi*xi + yi*yi ;
01214 radi2 = rhoi2 + zi*zi ;
01215 }
01216 else
01217 {
01218 sd = 0 ;
01219 xi = p.x() ;
01220 yi = p.y() ;
01221 zi = p.z() ;
01222 rhoi2 = rho2 ;
01223 radi2 = rad2 ;
01224 }
01225 if ( (radi2 <= tolORMax2)
01226 && (radi2 >= tolORMin2)
01227 && ((yi*cosCPhi-xi*sinCPhi) <= 0) )
01228 {
01229
01230
01231
01232
01233 if ( !fFullThetaSphere )
01234 {
01235 iTheta = std::atan2(std::sqrt(rhoi2),zi) ;
01236 if ( (iTheta >= tolSTheta) && (iTheta <= tolETheta) )
01237 {
01238
01239
01240
01241 if ((yi*cosCPhi-xi*sinCPhi) <= 0)
01242 {
01243 snxt = sd;
01244 }
01245 }
01246 }
01247 else
01248 {
01249 snxt = sd;
01250 }
01251 }
01252 }
01253 }
01254 }
01255
01256
01257
01258
01259 Comp = -( v.x()*sinEPhi-v.y()*cosEPhi ) ;
01260
01261 if (Comp < 0)
01262 {
01263 Dist = -(p.y()*cosEPhi-p.x()*sinEPhi) ;
01264 if ( Dist < halfCarTolerance )
01265 {
01266 sd = Dist/Comp ;
01267
01268 if ( sd < snxt )
01269 {
01270 if (sd > 0)
01271 {
01272 xi = p.x() + sd*v.x() ;
01273 yi = p.y() + sd*v.y() ;
01274 zi = p.z() + sd*v.z() ;
01275 rhoi2 = xi*xi + yi*yi ;
01276 radi2 = rhoi2 + zi*zi ;
01277 }
01278 else
01279 {
01280 sd = 0 ;
01281 xi = p.x() ;
01282 yi = p.y() ;
01283 zi = p.z() ;
01284 rhoi2 = rho2 ;
01285 radi2 = rad2 ;
01286 }
01287 if ( (radi2 <= tolORMax2)
01288 && (radi2 >= tolORMin2)
01289 && ((yi*cosCPhi-xi*sinCPhi) >= 0) )
01290 {
01291
01292
01293
01294
01295 if ( !fFullThetaSphere )
01296 {
01297 iTheta = std::atan2(std::sqrt(rhoi2),zi) ;
01298 if ( (iTheta >= tolSTheta) && (iTheta <= tolETheta) )
01299 {
01300
01301
01302
01303 if ((yi*cosCPhi-xi*sinCPhi) >= 0)
01304 {
01305 snxt = sd;
01306 }
01307 }
01308 }
01309 else
01310 {
01311 snxt = sd;
01312 }
01313 }
01314 }
01315 }
01316 }
01317 }
01318
01319
01320
01321 if ( !fFullThetaSphere )
01322 {
01323
01324
01325
01326
01327
01328
01329
01330
01331
01332
01333
01334
01335
01336
01337
01338
01339
01340
01341
01342
01343 if (fSTheta)
01344 {
01345 dist2STheta = rho2 - p.z()*p.z()*tanSTheta2 ;
01346 }
01347 else
01348 {
01349 dist2STheta = kInfinity ;
01350 }
01351 if ( eTheta < pi )
01352 {
01353 dist2ETheta=rho2-p.z()*p.z()*tanETheta2;
01354 }
01355 else
01356 {
01357 dist2ETheta=kInfinity;
01358 }
01359 if ( pTheta < tolSTheta )
01360 {
01361
01362
01363
01364 t1 = 1 - v.z()*v.z()*(1 + tanSTheta2) ;
01365 t2 = pDotV2d - p.z()*v.z()*tanSTheta2 ;
01366 if (t1)
01367 {
01368 b = t2/t1 ;
01369 c = dist2STheta/t1 ;
01370 d2 = b*b - c ;
01371
01372 if ( d2 >= 0 )
01373 {
01374 d = std::sqrt(d2) ;
01375 sd = -b - d ;
01376 zi = p.z() + sd*v.z();
01377
01378 if ( (sd < 0) || (zi*(fSTheta - halfpi) > 0) )
01379 {
01380 sd = -b+d;
01381 }
01382 if ((sd >= 0) && (sd < snxt))
01383 {
01384 xi = p.x() + sd*v.x();
01385 yi = p.y() + sd*v.y();
01386 zi = p.z() + sd*v.z();
01387 rhoi2 = xi*xi + yi*yi;
01388 radi2 = rhoi2 + zi*zi;
01389 if ( (radi2 <= tolORMax2)
01390 && (radi2 >= tolORMin2)
01391 && (zi*(fSTheta - halfpi) <= 0) )
01392 {
01393 if ( !fFullPhiSphere && rhoi2 )
01394 {
01395 cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
01396 if (cosPsi >= cosHDPhiOT)
01397 {
01398 snxt = sd;
01399 }
01400 }
01401 else
01402 {
01403 snxt = sd;
01404 }
01405 }
01406 }
01407 }
01408 }
01409
01410
01411
01412
01413 if ( eTheta < pi )
01414 {
01415 t1 = 1 - v.z()*v.z()*(1 + tanETheta2) ;
01416 t2 = pDotV2d - p.z()*v.z()*tanETheta2 ;
01417 if (t1)
01418 {
01419 b = t2/t1 ;
01420 c = dist2ETheta/t1 ;
01421 d2 = b*b - c ;
01422
01423 if (d2 >= 0)
01424 {
01425 d = std::sqrt(d2) ;
01426 sd = -b + d ;
01427
01428 if ( (sd >= 0) && (sd < snxt) )
01429 {
01430 xi = p.x() + sd*v.x() ;
01431 yi = p.y() + sd*v.y() ;
01432 zi = p.z() + sd*v.z() ;
01433 rhoi2 = xi*xi + yi*yi ;
01434 radi2 = rhoi2 + zi*zi ;
01435
01436 if ( (radi2 <= tolORMax2)
01437 && (radi2 >= tolORMin2)
01438 && (zi*(eTheta - halfpi) <= 0) )
01439 {
01440 if (!fFullPhiSphere && rhoi2)
01441 {
01442 cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
01443 if (cosPsi >= cosHDPhiOT)
01444 {
01445 snxt = sd;
01446 }
01447 }
01448 else
01449 {
01450 snxt = sd;
01451 }
01452 }
01453 }
01454 }
01455 }
01456 }
01457 }
01458 else if ( pTheta > tolETheta )
01459 {
01460
01461
01462
01463
01464 t1 = 1 - v.z()*v.z()*(1 + tanETheta2) ;
01465 t2 = pDotV2d - p.z()*v.z()*tanETheta2 ;
01466 if (t1)
01467 {
01468 b = t2/t1 ;
01469 c = dist2ETheta/t1 ;
01470 d2 = b*b - c ;
01471
01472 if (d2 >= 0)
01473 {
01474 d = std::sqrt(d2) ;
01475 sd = -b - d ;
01476 zi = p.z() + sd*v.z();
01477
01478 if ( (sd < 0) || (zi*(eTheta - halfpi) > 0) )
01479 {
01480 sd = -b + d ;
01481 }
01482 if ( (sd >= 0) && (sd < snxt) )
01483 {
01484 xi = p.x() + sd*v.x() ;
01485 yi = p.y() + sd*v.y() ;
01486 zi = p.z() + sd*v.z() ;
01487 rhoi2 = xi*xi + yi*yi ;
01488 radi2 = rhoi2 + zi*zi ;
01489
01490 if ( (radi2 <= tolORMax2)
01491 && (radi2 >= tolORMin2)
01492 && (zi*(eTheta - halfpi) <= 0) )
01493 {
01494 if (!fFullPhiSphere && rhoi2)
01495 {
01496 cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
01497 if (cosPsi >= cosHDPhiOT)
01498 {
01499 snxt = sd;
01500 }
01501 }
01502 else
01503 {
01504 snxt = sd;
01505 }
01506 }
01507 }
01508 }
01509 }
01510
01511
01512
01513
01514 if ( fSTheta )
01515 {
01516 t1 = 1 - v.z()*v.z()*(1 + tanSTheta2) ;
01517 t2 = pDotV2d - p.z()*v.z()*tanSTheta2 ;
01518 if (t1)
01519 {
01520 b = t2/t1 ;
01521 c = dist2STheta/t1 ;
01522 d2 = b*b - c ;
01523
01524 if (d2 >= 0)
01525 {
01526 d = std::sqrt(d2) ;
01527 sd = -b + d ;
01528
01529 if ( (sd >= 0) && (sd < snxt) )
01530 {
01531 xi = p.x() + sd*v.x() ;
01532 yi = p.y() + sd*v.y() ;
01533 zi = p.z() + sd*v.z() ;
01534 rhoi2 = xi*xi + yi*yi ;
01535 radi2 = rhoi2 + zi*zi ;
01536
01537 if ( (radi2 <= tolORMax2)
01538 && (radi2 >= tolORMin2)
01539 && (zi*(fSTheta - halfpi) <= 0) )
01540 {
01541 if (!fFullPhiSphere && rhoi2)
01542 {
01543 cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
01544 if (cosPsi >= cosHDPhiOT)
01545 {
01546 snxt = sd;
01547 }
01548 }
01549 else
01550 {
01551 snxt = sd;
01552 }
01553 }
01554 }
01555 }
01556 }
01557 }
01558 }
01559 else if ( (pTheta < tolSTheta + kAngTolerance)
01560 && (fSTheta > halfAngTolerance) )
01561 {
01562
01563
01564
01565
01566 t2 = pDotV2d - p.z()*v.z()*tanSTheta2 ;
01567 if ( (t2>=0 && tolIRMin2<rad2 && rad2<tolIRMax2 && fSTheta<halfpi)
01568 || (t2<0 && tolIRMin2<rad2 && rad2<tolIRMax2 && fSTheta>halfpi)
01569 || (v.z()<0 && tolIRMin2<rad2 && rad2<tolIRMax2 && fSTheta==halfpi) )
01570 {
01571 if (!fFullPhiSphere && rho2)
01572 {
01573 cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(rho2) ;
01574 if (cosPsi >= cosHDPhiIT)
01575 {
01576 return 0 ;
01577 }
01578 }
01579 else
01580 {
01581 return 0 ;
01582 }
01583 }
01584
01585
01586
01587 t1 = 1 - v.z()*v.z()*(1 + tanSTheta2) ;
01588 if (t1)
01589 {
01590 b = t2/t1 ;
01591 c = dist2STheta/t1 ;
01592 d2 = b*b - c ;
01593
01594 if (d2 >= 0)
01595 {
01596 d = std::sqrt(d2) ;
01597 sd = -b + d ;
01598 if ( (sd >= halfCarTolerance) && (sd < snxt) && (fSTheta < halfpi) )
01599 {
01600 xi = p.x() + sd*v.x() ;
01601 yi = p.y() + sd*v.y() ;
01602 zi = p.z() + sd*v.z() ;
01603 rhoi2 = xi*xi + yi*yi ;
01604 radi2 = rhoi2 + zi*zi ;
01605
01606 if ( (radi2 <= tolORMax2)
01607 && (radi2 >= tolORMin2)
01608 && (zi*(fSTheta - halfpi) <= 0) )
01609 {
01610 if ( !fFullPhiSphere && rhoi2 )
01611 {
01612 cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
01613 if ( cosPsi >= cosHDPhiOT )
01614 {
01615 snxt = sd;
01616 }
01617 }
01618 else
01619 {
01620 snxt = sd;
01621 }
01622 }
01623 }
01624 }
01625 }
01626 }
01627 else if ((pTheta > tolETheta-kAngTolerance) && (eTheta < pi-kAngTolerance))
01628 {
01629
01630
01631
01632
01633
01634 t2 = pDotV2d - p.z()*v.z()*tanETheta2 ;
01635
01636 if ( ((t2<0) && (eTheta < halfpi)
01637 && (tolIRMin2 < rad2) && (rad2 < tolIRMax2))
01638 || ((t2>=0) && (eTheta > halfpi)
01639 && (tolIRMin2 < rad2) && (rad2 < tolIRMax2))
01640 || ((v.z()>0) && (eTheta == halfpi)
01641 && (tolIRMin2 < rad2) && (rad2 < tolIRMax2)) )
01642 {
01643 if (!fFullPhiSphere && rho2)
01644 {
01645 cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(rho2) ;
01646 if (cosPsi >= cosHDPhiIT)
01647 {
01648 return 0 ;
01649 }
01650 }
01651 else
01652 {
01653 return 0 ;
01654 }
01655 }
01656
01657
01658
01659 t1 = 1 - v.z()*v.z()*(1 + tanETheta2) ;
01660 if (t1)
01661 {
01662 b = t2/t1 ;
01663 c = dist2ETheta/t1 ;
01664 d2 = b*b - c ;
01665
01666 if (d2 >= 0)
01667 {
01668 d = std::sqrt(d2) ;
01669 sd = -b + d ;
01670
01671 if ( (sd >= halfCarTolerance)
01672 && (sd < snxt) && (eTheta > halfpi) )
01673 {
01674 xi = p.x() + sd*v.x() ;
01675 yi = p.y() + sd*v.y() ;
01676 zi = p.z() + sd*v.z() ;
01677 rhoi2 = xi*xi + yi*yi ;
01678 radi2 = rhoi2 + zi*zi ;
01679
01680 if ( (radi2 <= tolORMax2)
01681 && (radi2 >= tolORMin2)
01682 && (zi*(eTheta - halfpi) <= 0) )
01683 {
01684 if (!fFullPhiSphere && rhoi2)
01685 {
01686 cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
01687 if (cosPsi >= cosHDPhiOT)
01688 {
01689 snxt = sd;
01690 }
01691 }
01692 else
01693 {
01694 snxt = sd;
01695 }
01696 }
01697 }
01698 }
01699 }
01700 }
01701 else
01702 {
01703
01704
01705
01706 t1 = 1 - v.z()*v.z()*(1 + tanSTheta2) ;
01707 t2 = pDotV2d - p.z()*v.z()*tanSTheta2 ;
01708 if (t1)
01709 {
01710 b = t2/t1;
01711 c = dist2STheta/t1 ;
01712 d2 = b*b - c ;
01713
01714 if (d2 >= 0)
01715 {
01716 d = std::sqrt(d2) ;
01717 sd = -b + d ;
01718
01719 if ((sd >= 0) && (sd < snxt))
01720 {
01721 xi = p.x() + sd*v.x() ;
01722 yi = p.y() + sd*v.y() ;
01723 zi = p.z() + sd*v.z() ;
01724 rhoi2 = xi*xi + yi*yi ;
01725 radi2 = rhoi2 + zi*zi ;
01726
01727 if ( (radi2 <= tolORMax2)
01728 && (radi2 >= tolORMin2)
01729 && (zi*(fSTheta - halfpi) <= 0) )
01730 {
01731 if (!fFullPhiSphere && rhoi2)
01732 {
01733 cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
01734 if (cosPsi >= cosHDPhiOT)
01735 {
01736 snxt = sd;
01737 }
01738 }
01739 else
01740 {
01741 snxt = sd;
01742 }
01743 }
01744 }
01745 }
01746 }
01747 t1 = 1 - v.z()*v.z()*(1 + tanETheta2) ;
01748 t2 = pDotV2d - p.z()*v.z()*tanETheta2 ;
01749 if (t1)
01750 {
01751 b = t2/t1 ;
01752 c = dist2ETheta/t1 ;
01753 d2 = b*b - c ;
01754
01755 if (d2 >= 0)
01756 {
01757 d = std::sqrt(d2) ;
01758 sd = -b + d;
01759
01760 if ((sd >= 0) && (sd < snxt))
01761 {
01762 xi = p.x() + sd*v.x() ;
01763 yi = p.y() + sd*v.y() ;
01764 zi = p.z() + sd*v.z() ;
01765 rhoi2 = xi*xi + yi*yi ;
01766 radi2 = rhoi2 + zi*zi ;
01767
01768 if ( (radi2 <= tolORMax2)
01769 && (radi2 >= tolORMin2)
01770 && (zi*(eTheta - halfpi) <= 0) )
01771 {
01772 if (!fFullPhiSphere && rhoi2)
01773 {
01774 cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ;
01775 if ( cosPsi >= cosHDPhiOT )
01776 {
01777 snxt = sd;
01778 }
01779 }
01780 else
01781 {
01782 snxt = sd;
01783 }
01784 }
01785 }
01786 }
01787 }
01788 }
01789 }
01790 return snxt;
01791 }
01792
01794
01795
01796
01797
01798
01799
01800
01801 G4double G4Sphere::DistanceToIn( const G4ThreeVector& p ) const
01802 {
01803 G4double safe=0.0,safeRMin,safeRMax,safePhi,safeTheta;
01804 G4double rho2,rds,rho;
01805 G4double cosPsi;
01806 G4double pTheta,dTheta1,dTheta2;
01807 rho2=p.x()*p.x()+p.y()*p.y();
01808 rds=std::sqrt(rho2+p.z()*p.z());
01809 rho=std::sqrt(rho2);
01810
01811
01812
01813
01814 if (fRmin)
01815 {
01816 safeRMin=fRmin-rds;
01817 safeRMax=rds-fRmax;
01818 if (safeRMin>safeRMax)
01819 {
01820 safe=safeRMin;
01821 }
01822 else
01823 {
01824 safe=safeRMax;
01825 }
01826 }
01827 else
01828 {
01829 safe=rds-fRmax;
01830 }
01831
01832
01833
01834
01835 if (!fFullPhiSphere && rho)
01836 {
01837
01838
01839 cosPsi=(p.x()*cosCPhi+p.y()*sinCPhi)/rho;
01840 if (cosPsi<std::cos(hDPhi))
01841 {
01842
01843
01844 if ((p.y()*cosCPhi-p.x()*sinCPhi)<=0)
01845 {
01846 safePhi=std::fabs(p.x()*sinSPhi-p.y()*cosSPhi);
01847 }
01848 else
01849 {
01850 safePhi=std::fabs(p.x()*sinEPhi-p.y()*cosEPhi);
01851 }
01852 if (safePhi>safe) { safe=safePhi; }
01853 }
01854 }
01855
01856
01857
01858 if ((rds!=0.0) && (!fFullThetaSphere))
01859 {
01860 pTheta=std::acos(p.z()/rds);
01861 if (pTheta<0) { pTheta+=pi; }
01862 dTheta1=fSTheta-pTheta;
01863 dTheta2=pTheta-eTheta;
01864 if (dTheta1>dTheta2)
01865 {
01866 if (dTheta1>=0)
01867 {
01868 safeTheta=rds*std::sin(dTheta1);
01869 if (safe<=safeTheta)
01870 {
01871 safe=safeTheta;
01872 }
01873 }
01874 }
01875 else
01876 {
01877 if (dTheta2>=0)
01878 {
01879 safeTheta=rds*std::sin(dTheta2);
01880 if (safe<=safeTheta)
01881 {
01882 safe=safeTheta;
01883 }
01884 }
01885 }
01886 }
01887
01888 if (safe<0) { safe=0; }
01889 return safe;
01890 }
01891
01893
01894
01895
01896
01897 G4double G4Sphere::DistanceToOut( const G4ThreeVector& p,
01898 const G4ThreeVector& v,
01899 const G4bool calcNorm,
01900 G4bool *validNorm,
01901 G4ThreeVector *n ) const
01902 {
01903 G4double snxt = kInfinity;
01904 G4double sphi= kInfinity,stheta= kInfinity;
01905 ESide side=kNull,sidephi=kNull,sidetheta=kNull;
01906
01907 static const G4double halfCarTolerance = kCarTolerance*0.5;
01908 static const G4double halfAngTolerance = kAngTolerance*0.5;
01909 const G4double halfRmaxTolerance = fRmaxTolerance*0.5;
01910 const G4double halfRminTolerance = fRminTolerance*0.5;
01911 const G4double Rmax_plus = fRmax + halfRmaxTolerance;
01912 const G4double Rmin_minus = (fRmin) ? fRmin-halfRminTolerance : 0;
01913 G4double t1,t2;
01914 G4double b,c,d;
01915
01916
01917
01918 G4double pDistS,compS,pDistE,compE,sphi2,vphi;
01919
01920 G4double rho2,rad2,pDotV2d,pDotV3d;
01921
01922 G4double xi,yi,zi;
01923
01924
01925
01926 G4double rhoSecTheta;
01927 G4double dist2STheta, dist2ETheta, distTheta;
01928 G4double d2,sd;
01929
01930
01931
01932 rho2 = p.x()*p.x()+p.y()*p.y();
01933 rad2 = rho2+p.z()*p.z();
01934
01935 pDotV2d = p.x()*v.x()+p.y()*v.y();
01936 pDotV3d = pDotV2d+p.z()*v.z();
01937
01938
01939
01940
01941
01942
01943
01944
01945
01946
01947
01948
01949
01950
01951
01952
01953
01954 if( (rad2 <= Rmax_plus*Rmax_plus) && (rad2 >= Rmin_minus*Rmin_minus) )
01955 {
01956 c = rad2 - fRmax*fRmax;
01957
01958 if (c < fRmaxTolerance*fRmax)
01959 {
01960
01961
01962
01963
01964
01965
01966
01967
01968
01969 d2 = pDotV3d*pDotV3d - c;
01970
01971 if( (c >- fRmaxTolerance*fRmax)
01972 && ((pDotV3d >=0) || (d2 < 0)) )
01973
01974 {
01975 if(calcNorm)
01976 {
01977 *validNorm = true ;
01978 *n = G4ThreeVector(p.x()/fRmax,p.y()/fRmax,p.z()/fRmax) ;
01979 }
01980 return snxt = 0;
01981 }
01982 else
01983 {
01984 snxt = -pDotV3d+std::sqrt(d2);
01985 side = kRMax ;
01986 }
01987 }
01988
01989
01990
01991
01992
01993 if (fRmin)
01994 {
01995 c = rad2 - fRmin*fRmin;
01996 d2 = pDotV3d*pDotV3d - c;
01997
01998 if (c >- fRminTolerance*fRmin)
01999 {
02000 if ( (c < fRminTolerance*fRmin)
02001 && (d2 >= fRminTolerance*fRmin) && (pDotV3d < 0) )
02002 {
02003 if(calcNorm) { *validNorm = false; }
02004 return snxt = 0 ;
02005 }
02006 else
02007 {
02008 if ( d2 >= 0. )
02009 {
02010 sd = -pDotV3d-std::sqrt(d2);
02011
02012 if ( sd >= 0. )
02013 {
02014 snxt = sd ;
02015 side = kRMin ;
02016 }
02017 }
02018 }
02019 }
02020 }
02021 }
02022
02023
02024
02025 if ( !fFullThetaSphere )
02026 {
02027
02028
02029
02030
02031
02032
02033
02034
02035
02036
02037
02038
02039
02040
02041
02042
02043
02044
02045
02046
02047
02048
02049 if(fSTheta)
02050 {
02051 if( std::fabs(tanSTheta) > 5./kAngTolerance )
02052 {
02053 if( v.z() > 0. )
02054 {
02055 if ( std::fabs( p.z() ) <= halfRmaxTolerance )
02056 {
02057 if(calcNorm)
02058 {
02059 *validNorm = true;
02060 *n = G4ThreeVector(0.,0.,1.);
02061 }
02062 return snxt = 0 ;
02063 }
02064 stheta = -p.z()/v.z();
02065 sidetheta = kSTheta;
02066 }
02067 }
02068 else
02069 {
02070 t1 = 1-v.z()*v.z()*(1+tanSTheta2);
02071 t2 = pDotV2d-p.z()*v.z()*tanSTheta2;
02072 dist2STheta = rho2-p.z()*p.z()*tanSTheta2;
02073
02074 distTheta = std::sqrt(rho2)-p.z()*tanSTheta;
02075
02076 if( std::fabs(t1) < halfAngTolerance )
02077 {
02078 if( v.z() > 0. )
02079 {
02080 if(std::fabs(distTheta) < halfRmaxTolerance)
02081 {
02082 if( (fSTheta < halfpi) && (p.z() > 0.) )
02083 {
02084 if( calcNorm ) { *validNorm = false; }
02085 return snxt = 0.;
02086 }
02087 else if( (fSTheta > halfpi) && (p.z() <= 0) )
02088 {
02089 if( calcNorm )
02090 {
02091 *validNorm = true;
02092 if (rho2)
02093 {
02094 rhoSecTheta = std::sqrt(rho2*(1+tanSTheta2));
02095
02096 *n = G4ThreeVector( p.x()/rhoSecTheta,
02097 p.y()/rhoSecTheta,
02098 std::sin(fSTheta) );
02099 }
02100 else *n = G4ThreeVector(0.,0.,1.);
02101 }
02102 return snxt = 0.;
02103 }
02104 }
02105 stheta = -0.5*dist2STheta/t2;
02106 sidetheta = kSTheta;
02107 }
02108 }
02109 else
02110 {
02111 if( std::fabs(distTheta) < halfRmaxTolerance )
02112 {
02113 if( (fSTheta > halfpi) && (t2 >= 0.) )
02114 {
02115 if( calcNorm )
02116 {
02117 *validNorm = true;
02118 if (rho2)
02119 {
02120 rhoSecTheta = std::sqrt(rho2*(1+tanSTheta2));
02121
02122 *n = G4ThreeVector( p.x()/rhoSecTheta,
02123 p.y()/rhoSecTheta,
02124 std::sin(fSTheta) );
02125 }
02126 else { *n = G4ThreeVector(0.,0.,1.); }
02127 }
02128 return snxt = 0.;
02129 }
02130 else if( (fSTheta < halfpi) && (t2 < 0.) && (p.z() >=0.) )
02131 {
02132 if( calcNorm ) { *validNorm = false; }
02133 return snxt = 0.;
02134 }
02135 }
02136 b = t2/t1;
02137 c = dist2STheta/t1;
02138 d2 = b*b - c ;
02139
02140 if ( d2 >= 0. )
02141 {
02142 d = std::sqrt(d2);
02143
02144 if( fSTheta > halfpi )
02145 {
02146 sd = -b - d;
02147
02148 if ( ((std::fabs(s) < halfRmaxTolerance) && (t2 < 0.))
02149 || (sd < 0.) || ( (sd > 0.) && (p.z() + sd*v.z() > 0.) ) )
02150 {
02151 sd = -b + d ;
02152 }
02153 if( (sd > halfRmaxTolerance) && (p.z() + sd*v.z() <= 0.) )
02154 {
02155 stheta = sd;
02156 sidetheta = kSTheta;
02157 }
02158 }
02159 else
02160 {
02161 sd = -b - d;
02162
02163 if ( ( (std::fabs(sd) < halfRmaxTolerance) && (t2 >= 0.) )
02164 || (sd < 0.) || ( (sd > 0.) && (p.z() + sd*v.z() < 0.) ) )
02165 {
02166 sd = -b + d ;
02167 }
02168 if( (sd > halfRmaxTolerance) && (p.z() + sd*v.z() >= 0.) )
02169 {
02170 stheta = sd;
02171 sidetheta = kSTheta;
02172 }
02173 }
02174 }
02175 }
02176 }
02177 }
02178 if (eTheta < pi)
02179 {
02180 if( std::fabs(tanETheta) > 5./kAngTolerance )
02181 {
02182 if( v.z() < 0. )
02183 {
02184 if ( std::fabs( p.z() ) <= halfRmaxTolerance )
02185 {
02186 if(calcNorm)
02187 {
02188 *validNorm = true;
02189 *n = G4ThreeVector(0.,0.,-1.);
02190 }
02191 return snxt = 0 ;
02192 }
02193 sd = -p.z()/v.z();
02194
02195 if( sd < stheta )
02196 {
02197 stheta = sd;
02198 sidetheta = kETheta;
02199 }
02200 }
02201 }
02202 else
02203 {
02204 t1 = 1-v.z()*v.z()*(1+tanETheta2);
02205 t2 = pDotV2d-p.z()*v.z()*tanETheta2;
02206 dist2ETheta = rho2-p.z()*p.z()*tanETheta2;
02207
02208 distTheta = std::sqrt(rho2)-p.z()*tanETheta;
02209
02210 if( std::fabs(t1) < halfAngTolerance )
02211 {
02212 if( v.z() < 0. )
02213 {
02214 if(std::fabs(distTheta) < halfRmaxTolerance)
02215 {
02216 if( (eTheta > halfpi) && (p.z() < 0.) )
02217 {
02218 if( calcNorm ) { *validNorm = false; }
02219 return snxt = 0.;
02220 }
02221 else if ( (eTheta < halfpi) && (p.z() >= 0) )
02222 {
02223 if( calcNorm )
02224 {
02225 *validNorm = true;
02226 if (rho2)
02227 {
02228 rhoSecTheta = std::sqrt(rho2*(1+tanETheta2));
02229 *n = G4ThreeVector( p.x()/rhoSecTheta,
02230 p.y()/rhoSecTheta,
02231 -sinETheta );
02232 }
02233 else { *n = G4ThreeVector(0.,0.,-1.); }
02234 }
02235 return snxt = 0.;
02236 }
02237 }
02238 sd = -0.5*dist2ETheta/t2;
02239
02240 if( sd < stheta )
02241 {
02242 stheta = sd;
02243 sidetheta = kETheta;
02244 }
02245 }
02246 }
02247 else
02248 {
02249 if ( std::fabs(distTheta) < halfRmaxTolerance )
02250 {
02251 if( (eTheta < halfpi) && (t2 >= 0.) )
02252 {
02253 if( calcNorm )
02254 {
02255 *validNorm = true;
02256 if (rho2)
02257 {
02258 rhoSecTheta = std::sqrt(rho2*(1+tanETheta2));
02259 *n = G4ThreeVector( p.x()/rhoSecTheta,
02260 p.y()/rhoSecTheta,
02261 -sinETheta );
02262 }
02263 else *n = G4ThreeVector(0.,0.,-1.);
02264 }
02265 return snxt = 0.;
02266 }
02267 else if ( (eTheta > halfpi)
02268 && (t2 < 0.) && (p.z() <=0.) )
02269 {
02270 if( calcNorm ) { *validNorm = false; }
02271 return snxt = 0.;
02272 }
02273 }
02274 b = t2/t1;
02275 c = dist2ETheta/t1;
02276 d2 = b*b - c ;
02277
02278 if ( d2 >= 0. )
02279 {
02280 d = std::sqrt(d2);
02281
02282 if( eTheta < halfpi )
02283 {
02284 sd = -b - d;
02285
02286 if( ((std::fabs(sd) < halfRmaxTolerance) && (t2 < 0.))
02287 || (sd < 0.) )
02288 {
02289 sd = -b + d ;
02290 }
02291 if( sd > halfRmaxTolerance )
02292 {
02293 if( sd < stheta )
02294 {
02295 stheta = sd;
02296 sidetheta = kETheta;
02297 }
02298 }
02299 }
02300 else
02301 {
02302 sd = -b - d;
02303
02304 if ( ((std::fabs(sd) < halfRmaxTolerance) && (t2 >= 0.))
02305 || (sd < 0.) || ( (sd > 0.) && (p.z() + sd*v.z() > 0.) ) )
02306 {
02307 sd = -b + d ;
02308 }
02309 if( (sd > halfRmaxTolerance) && (p.z() + sd*v.z() <= 0.) )
02310 {
02311 if( sd < stheta )
02312 {
02313 stheta = sd;
02314 sidetheta = kETheta;
02315 }
02316 }
02317 }
02318 }
02319 }
02320 }
02321 }
02322
02323 }
02324
02325
02326
02327 if ( !fFullPhiSphere )
02328 {
02329 if ( p.x() || p.y() )
02330 {
02331
02332
02333 pDistS=p.x()*sinSPhi-p.y()*cosSPhi;
02334 pDistE=-p.x()*sinEPhi+p.y()*cosEPhi;
02335
02336
02337
02338 compS = -sinSPhi*v.x()+cosSPhi*v.y() ;
02339 compE = sinEPhi*v.x()-cosEPhi*v.y() ;
02340 sidephi = kNull ;
02341
02342 if ( (pDistS <= 0) && (pDistE <= 0) )
02343 {
02344
02345
02346 if ( compS < 0 )
02347 {
02348 sphi = pDistS/compS ;
02349 xi = p.x()+sphi*v.x() ;
02350 yi = p.y()+sphi*v.y() ;
02351
02352
02353
02354 if( (std::fabs(xi)<=kCarTolerance) && (std::fabs(yi)<=kCarTolerance) )
02355 {
02356 vphi = std::atan2(v.y(),v.x());
02357 sidephi = kSPhi;
02358 if ( ( (fSPhi-halfAngTolerance) <= vphi)
02359 && ( (ePhi+halfAngTolerance) >= vphi) )
02360 {
02361 sphi = kInfinity;
02362 }
02363 }
02364 else if ( ( yi*cosCPhi - xi*sinCPhi ) >= 0 )
02365 {
02366 sphi=kInfinity;
02367 }
02368 else
02369 {
02370 sidephi = kSPhi ;
02371 if ( pDistS > -halfCarTolerance) { sphi = 0; }
02372 }
02373 }
02374 else { sphi = kInfinity; }
02375
02376 if ( compE < 0 )
02377 {
02378 sphi2=pDistE/compE ;
02379 if (sphi2 < sphi)
02380 {
02381 xi = p.x()+sphi2*v.x() ;
02382 yi = p.y()+sphi2*v.y() ;
02383
02384
02385
02386 if ((std::fabs(xi)<=kCarTolerance) && (std::fabs(yi)<=kCarTolerance))
02387 {
02388
02389
02390 vphi = std::atan2(v.y(),v.x()) ;
02391
02392 if( !((fSPhi-halfAngTolerance <= vphi)
02393 &&(fSPhi+fDPhi+halfAngTolerance >= vphi)) )
02394 {
02395 sidephi = kEPhi;
02396 if ( pDistE <= -halfCarTolerance ) { sphi = sphi2; }
02397 else { sphi = 0.0; }
02398 }
02399 }
02400 else if ((yi*cosCPhi-xi*sinCPhi)>=0)
02401 {
02402 sidephi = kEPhi ;
02403 if ( pDistE <= -halfCarTolerance )
02404 {
02405 sphi=sphi2;
02406 }
02407 else
02408 {
02409 sphi = 0 ;
02410 }
02411 }
02412 }
02413 }
02414 }
02415 else if ((pDistS >= 0) && (pDistE >= 0))
02416 {
02417 if ( pDistS <= pDistE )
02418 {
02419 sidephi = kSPhi ;
02420 }
02421 else
02422 {
02423 sidephi = kEPhi ;
02424 }
02425 if ( fDPhi > pi )
02426 {
02427 if ( (compS < 0) && (compE < 0) ) { sphi = 0; }
02428 else { sphi = kInfinity; }
02429 }
02430 else
02431 {
02432
02433
02434
02435 if ( (compS >= 0) && (compE >= 0) ) { sphi = kInfinity; }
02436 else { sphi = 0; }
02437 }
02438 }
02439 else if ( (pDistS > 0) && (pDistE < 0) )
02440 {
02441
02442
02443 if ( fDPhi > pi )
02444 {
02445 if ( compE < 0 )
02446 {
02447 sphi = pDistE/compE ;
02448 xi = p.x() + sphi*v.x() ;
02449 yi = p.y() + sphi*v.y() ;
02450
02451
02452
02453
02454 if( (std::fabs(xi)<=kCarTolerance)&&(std::fabs(yi)<=kCarTolerance) )
02455 {
02456 vphi = std::atan2(v.y(),v.x());
02457 sidephi = kSPhi;
02458 if ( ( (fSPhi-halfAngTolerance) <= vphi)
02459 && ( (ePhi+halfAngTolerance) >= vphi) )
02460 {
02461 sphi = kInfinity;
02462 }
02463 }
02464 else if ( ( yi*cosCPhi - xi*sinCPhi ) <= 0 )
02465 {
02466 sphi = kInfinity ;
02467 }
02468 else
02469 {
02470 sidephi = kEPhi ;
02471 if ( pDistE > -halfCarTolerance ) { sphi = 0.; }
02472 }
02473 }
02474 else
02475 {
02476 sphi = kInfinity ;
02477 }
02478 }
02479 else
02480 {
02481 if ( compS >= 0 )
02482 {
02483 if ( compE < 0 )
02484 {
02485 sphi = pDistE/compE ;
02486 xi = p.x() + sphi*v.x() ;
02487 yi = p.y() + sphi*v.y() ;
02488
02489
02490
02491
02492 if( (std::fabs(xi)<=kCarTolerance)&&(std::fabs(yi)<=kCarTolerance) )
02493 {
02494 vphi = std::atan2(v.y(),v.x());
02495 sidephi = kSPhi;
02496 if ( ( (fSPhi-halfAngTolerance) <= vphi)
02497 && ( (ePhi+halfAngTolerance) >= vphi) )
02498 {
02499 sphi = kInfinity;
02500 }
02501 }
02502 else if ( ( yi*cosCPhi - xi*sinCPhi) <= 0 )
02503 {
02504 sphi=kInfinity;
02505 }
02506 else
02507 {
02508 sidephi = kEPhi ;
02509 }
02510 }
02511 else sphi=kInfinity;
02512 }
02513 else
02514 {
02515 sidephi = kSPhi ;
02516 sphi = 0 ;
02517 }
02518 }
02519 }
02520 else
02521 {
02522
02523
02524
02525 if ( fDPhi > pi )
02526 {
02527 if ( compS < 0 )
02528 {
02529 sphi=pDistS/compS;
02530 xi=p.x()+sphi*v.x();
02531 yi=p.y()+sphi*v.y();
02532
02533
02534
02535
02536 if( (std::fabs(xi)<=kCarTolerance)&&(std::fabs(yi)<=kCarTolerance) )
02537 {
02538 vphi = std::atan2(v.y(),v.x()) ;
02539 sidephi = kSPhi;
02540 if ( ( (fSPhi-halfAngTolerance) <= vphi)
02541 && ( (ePhi+halfAngTolerance) >= vphi) )
02542 {
02543 sphi = kInfinity;
02544 }
02545 }
02546 else if ( ( yi*cosCPhi - xi*sinCPhi ) >= 0 )
02547 {
02548 sphi = kInfinity ;
02549 }
02550 else
02551 {
02552 sidephi = kSPhi ;
02553 if ( pDistS > -halfCarTolerance ) { sphi = 0; }
02554 }
02555 }
02556 else
02557 {
02558 sphi = kInfinity ;
02559 }
02560 }
02561 else
02562 {
02563 if ( compE >= 0 )
02564 {
02565 if ( compS < 0 )
02566 {
02567 sphi = pDistS/compS ;
02568 xi = p.x()+sphi*v.x() ;
02569 yi = p.y()+sphi*v.y() ;
02570
02571
02572
02573
02574 if((std::fabs(xi)<=kCarTolerance) && (std::fabs(yi)<=kCarTolerance))
02575 {
02576 vphi = std::atan2(v.y(),v.x()) ;
02577 sidephi = kSPhi;
02578 if ( ( (fSPhi-halfAngTolerance) <= vphi)
02579 && ( (ePhi+halfAngTolerance) >= vphi) )
02580 {
02581 sphi = kInfinity;
02582 }
02583 }
02584 else if ( ( yi*cosCPhi - xi*sinCPhi ) >= 0 )
02585 {
02586 sphi = kInfinity ;
02587 }
02588 else
02589 {
02590 sidephi = kSPhi ;
02591 }
02592 }
02593 else
02594 {
02595 sphi = kInfinity ;
02596 }
02597 }
02598 else
02599 {
02600 sidephi = kEPhi ;
02601 sphi = 0 ;
02602 }
02603 }
02604 }
02605 }
02606 else
02607 {
02608
02609
02610
02611 if ( v.x() || v.y() )
02612 {
02613 vphi = std::atan2(v.y(),v.x()) ;
02614 if ((fSPhi-halfAngTolerance < vphi) && (vphi < ePhi+halfAngTolerance))
02615 {
02616 sphi = kInfinity;
02617 }
02618 else
02619 {
02620 sidephi = kSPhi ;
02621 sphi = 0 ;
02622 }
02623 }
02624 else
02625 {
02626 sphi = kInfinity ;
02627 }
02628 }
02629 if ( sphi < snxt )
02630 {
02631 snxt = sphi ;
02632 side = sidephi ;
02633 }
02634 }
02635 if (stheta < snxt )
02636 {
02637 snxt = stheta ;
02638 side = sidetheta ;
02639 }
02640
02641 if (calcNorm)
02642 {
02643 switch( side )
02644 {
02645 case kRMax:
02646 xi=p.x()+snxt*v.x();
02647 yi=p.y()+snxt*v.y();
02648 zi=p.z()+snxt*v.z();
02649 *n=G4ThreeVector(xi/fRmax,yi/fRmax,zi/fRmax);
02650 *validNorm=true;
02651 break;
02652
02653 case kRMin:
02654 *validNorm=false;
02655 break;
02656
02657 case kSPhi:
02658 if ( fDPhi <= pi )
02659 {
02660 *n=G4ThreeVector(sinSPhi,-cosSPhi,0);
02661 *validNorm=true;
02662 }
02663 else { *validNorm=false; }
02664 break ;
02665
02666 case kEPhi:
02667 if ( fDPhi <= pi )
02668 {
02669 *n=G4ThreeVector(-sinEPhi,cosEPhi,0);
02670 *validNorm=true;
02671 }
02672 else { *validNorm=false; }
02673 break;
02674
02675 case kSTheta:
02676 if( fSTheta == halfpi )
02677 {
02678 *n=G4ThreeVector(0.,0.,1.);
02679 *validNorm=true;
02680 }
02681 else if ( fSTheta > halfpi )
02682 {
02683 xi = p.x() + snxt*v.x();
02684 yi = p.y() + snxt*v.y();
02685 rho2=xi*xi+yi*yi;
02686 if (rho2)
02687 {
02688 rhoSecTheta = std::sqrt(rho2*(1+tanSTheta2));
02689 *n = G4ThreeVector( xi/rhoSecTheta, yi/rhoSecTheta,
02690 -tanSTheta/std::sqrt(1+tanSTheta2));
02691 }
02692 else
02693 {
02694 *n = G4ThreeVector(0.,0.,1.);
02695 }
02696 *validNorm=true;
02697 }
02698 else { *validNorm=false; }
02699 break;
02700
02701 case kETheta:
02702 if( eTheta == halfpi )
02703 {
02704 *n = G4ThreeVector(0.,0.,-1.);
02705 *validNorm = true;
02706 }
02707 else if ( eTheta < halfpi )
02708 {
02709 xi=p.x()+snxt*v.x();
02710 yi=p.y()+snxt*v.y();
02711 rho2=xi*xi+yi*yi;
02712 if (rho2)
02713 {
02714 rhoSecTheta = std::sqrt(rho2*(1+tanETheta2));
02715 *n = G4ThreeVector( xi/rhoSecTheta, yi/rhoSecTheta,
02716 -tanETheta/std::sqrt(1+tanETheta2) );
02717 }
02718 else
02719 {
02720 *n = G4ThreeVector(0.,0.,-1.);
02721 }
02722 *validNorm=true;
02723 }
02724 else { *validNorm=false; }
02725 break;
02726
02727 default:
02728 G4cout << G4endl;
02729 DumpInfo();
02730 std::ostringstream message;
02731 G4int oldprc = message.precision(16);
02732 message << "Undefined side for valid surface normal to solid."
02733 << G4endl
02734 << "Position:" << G4endl << G4endl
02735 << "p.x() = " << p.x()/mm << " mm" << G4endl
02736 << "p.y() = " << p.y()/mm << " mm" << G4endl
02737 << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl
02738 << "Direction:" << G4endl << G4endl
02739 << "v.x() = " << v.x() << G4endl
02740 << "v.y() = " << v.y() << G4endl
02741 << "v.z() = " << v.z() << G4endl << G4endl
02742 << "Proposed distance :" << G4endl << G4endl
02743 << "snxt = " << snxt/mm << " mm" << G4endl;
02744 message.precision(oldprc);
02745 G4Exception("G4Sphere::DistanceToOut(p,v,..)",
02746 "GeomSolids1002", JustWarning, message);
02747 break;
02748 }
02749 }
02750 if (snxt == kInfinity)
02751 {
02752 G4cout << G4endl;
02753 DumpInfo();
02754 std::ostringstream message;
02755 G4int oldprc = message.precision(16);
02756 message << "Logic error: snxt = kInfinity ???" << G4endl
02757 << "Position:" << G4endl << G4endl
02758 << "p.x() = " << p.x()/mm << " mm" << G4endl
02759 << "p.y() = " << p.y()/mm << " mm" << G4endl
02760 << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl
02761 << "Rp = "<< std::sqrt( p.x()*p.x()+p.y()*p.y()+p.z()*p.z() )/mm
02762 << " mm" << G4endl << G4endl
02763 << "Direction:" << G4endl << G4endl
02764 << "v.x() = " << v.x() << G4endl
02765 << "v.y() = " << v.y() << G4endl
02766 << "v.z() = " << v.z() << G4endl << G4endl
02767 << "Proposed distance :" << G4endl << G4endl
02768 << "snxt = " << snxt/mm << " mm" << G4endl;
02769 message.precision(oldprc);
02770 G4Exception("G4Sphere::DistanceToOut(p,v,..)",
02771 "GeomSolids1002", JustWarning, message);
02772 }
02773
02774 return snxt;
02775 }
02776
02778
02779
02780
02781 G4double G4Sphere::DistanceToOut( const G4ThreeVector& p ) const
02782 {
02783 G4double safe=0.0,safeRMin,safeRMax,safePhi,safeTheta;
02784 G4double rho2,rds,rho;
02785 G4double pTheta,dTheta1,dTheta2;
02786 rho2=p.x()*p.x()+p.y()*p.y();
02787 rds=std::sqrt(rho2+p.z()*p.z());
02788 rho=std::sqrt(rho2);
02789
02790 #ifdef G4CSGDEBUG
02791 if( Inside(p) == kOutside )
02792 {
02793 G4int old_prc = G4cout.precision(16);
02794 G4cout << G4endl;
02795 DumpInfo();
02796 G4cout << "Position:" << G4endl << G4endl ;
02797 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ;
02798 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ;
02799 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ;
02800 G4cout.precision(old_prc) ;
02801 G4Exception("G4Sphere::DistanceToOut(p)",
02802 "GeomSolids1002", JustWarning, "Point p is outside !?" );
02803 }
02804 #endif
02805
02806
02807
02808
02809 if (fRmin)
02810 {
02811 safeRMin=rds-fRmin;
02812 safeRMax=fRmax-rds;
02813 if (safeRMin<safeRMax)
02814 {
02815 safe=safeRMin;
02816 }
02817 else
02818 {
02819 safe=safeRMax;
02820 }
02821 }
02822 else
02823 {
02824 safe=fRmax-rds;
02825 }
02826
02827
02828
02829
02830 if (!fFullPhiSphere && rho)
02831 {
02832 if ((p.y()*cosCPhi-p.x()*sinCPhi)<=0)
02833 {
02834 safePhi=-(p.x()*sinSPhi-p.y()*cosSPhi);
02835 }
02836 else
02837 {
02838 safePhi=(p.x()*sinEPhi-p.y()*cosEPhi);
02839 }
02840 if (safePhi<safe) { safe=safePhi; }
02841 }
02842
02843
02844
02845
02846 if (rds)
02847 {
02848 pTheta=std::acos(p.z()/rds);
02849 if (pTheta<0) { pTheta+=pi; }
02850 dTheta1=pTheta-fSTheta;
02851 dTheta2=eTheta-pTheta;
02852 if (dTheta1<dTheta2) { safeTheta=rds*std::sin(dTheta1); }
02853 else { safeTheta=rds*std::sin(dTheta2); }
02854 if (safe>safeTheta) { safe=safeTheta; }
02855 }
02856
02857 if (safe<0) { safe=0; }
02858 return safe;
02859 }
02860
02862
02863
02864
02865
02866
02867
02868
02869
02870
02871
02872 G4ThreeVectorList*
02873 G4Sphere::CreateRotatedVertices( const G4AffineTransform& pTransform,
02874 G4int& noPolygonVertices ) const
02875 {
02876 G4ThreeVectorList *vertices;
02877 G4ThreeVector vertex;
02878 G4double meshAnglePhi,meshRMax,crossAnglePhi,
02879 coscrossAnglePhi,sincrossAnglePhi,sAnglePhi;
02880 G4double meshTheta,crossTheta,startTheta;
02881 G4double rMaxX,rMaxY,rMinX,rMinY,rMinZ,rMaxZ;
02882 G4int crossSectionPhi,noPhiCrossSections,crossSectionTheta,noThetaSections;
02883
02884
02885
02886 noPhiCrossSections = G4int(fDPhi/kMeshAngleDefault)+1;
02887
02888 if (noPhiCrossSections<kMinMeshSections)
02889 {
02890 noPhiCrossSections=kMinMeshSections;
02891 }
02892 else if (noPhiCrossSections>kMaxMeshSections)
02893 {
02894 noPhiCrossSections=kMaxMeshSections;
02895 }
02896 meshAnglePhi=fDPhi/(noPhiCrossSections-1);
02897
02898
02899
02900
02901 if (fFullPhiSphere)
02902 {
02903 sAnglePhi = -meshAnglePhi*0.5;
02904 }
02905 else
02906 {
02907 sAnglePhi=fSPhi;
02908 }
02909
02910
02911
02912 noThetaSections = G4int(fDTheta/kMeshAngleDefault)+1;
02913
02914 if (noThetaSections<kMinMeshSections)
02915 {
02916 noThetaSections=kMinMeshSections;
02917 }
02918 else if (noThetaSections>kMaxMeshSections)
02919 {
02920 noThetaSections=kMaxMeshSections;
02921 }
02922 meshTheta=fDTheta/(noThetaSections-1);
02923
02924
02925
02926
02927 if (fFullThetaSphere)
02928 {
02929 startTheta = -meshTheta*0.5;
02930 }
02931 else
02932 {
02933 startTheta=fSTheta;
02934 }
02935
02936 meshRMax = (meshAnglePhi >= meshTheta) ?
02937 fRmax/std::cos(meshAnglePhi*0.5) : fRmax/std::cos(meshTheta*0.5);
02938 G4double* cosCrossTheta = new G4double[noThetaSections];
02939 G4double* sinCrossTheta = new G4double[noThetaSections];
02940 vertices=new G4ThreeVectorList();
02941 if (vertices && cosCrossTheta && sinCrossTheta)
02942 {
02943 vertices->reserve(noPhiCrossSections*(noThetaSections*2));
02944 for (crossSectionPhi=0;
02945 crossSectionPhi<noPhiCrossSections; crossSectionPhi++)
02946 {
02947 crossAnglePhi=sAnglePhi+crossSectionPhi*meshAnglePhi;
02948 coscrossAnglePhi=std::cos(crossAnglePhi);
02949 sincrossAnglePhi=std::sin(crossAnglePhi);
02950 for (crossSectionTheta=0;
02951 crossSectionTheta<noThetaSections;crossSectionTheta++)
02952 {
02953
02954
02955 crossTheta=startTheta+crossSectionTheta*meshTheta;
02956 cosCrossTheta[crossSectionTheta]=std::cos(crossTheta);
02957 sinCrossTheta[crossSectionTheta]=std::sin(crossTheta);
02958
02959 rMinX=fRmin*sinCrossTheta[crossSectionTheta]*coscrossAnglePhi;
02960 rMinY=fRmin*sinCrossTheta[crossSectionTheta]*sincrossAnglePhi;
02961 rMinZ=fRmin*cosCrossTheta[crossSectionTheta];
02962
02963 vertex=G4ThreeVector(rMinX,rMinY,rMinZ);
02964 vertices->push_back(pTransform.TransformPoint(vertex));
02965
02966 }
02967
02968 for (crossSectionTheta=noThetaSections-1;
02969 crossSectionTheta>=0; crossSectionTheta--)
02970 {
02971 rMaxX=meshRMax*sinCrossTheta[crossSectionTheta]*coscrossAnglePhi;
02972 rMaxY=meshRMax*sinCrossTheta[crossSectionTheta]*sincrossAnglePhi;
02973 rMaxZ=meshRMax*cosCrossTheta[crossSectionTheta];
02974
02975 vertex=G4ThreeVector(rMaxX,rMaxY,rMaxZ);
02976 vertices->push_back(pTransform.TransformPoint(vertex));
02977
02978 }
02979 }
02980 noPolygonVertices = noThetaSections*2 ;
02981 }
02982 else
02983 {
02984 DumpInfo();
02985 G4Exception("G4Sphere::CreateRotatedVertices()",
02986 "GeomSolids0003", FatalException,
02987 "Error in allocation of vertices. Out of memory !");
02988 }
02989
02990 delete [] cosCrossTheta;
02991 delete [] sinCrossTheta;
02992
02993 return vertices;
02994 }
02995
02997
02998
02999
03000 G4GeometryType G4Sphere::GetEntityType() const
03001 {
03002 return G4String("G4Sphere");
03003 }
03004
03006
03007
03008
03009 G4VSolid* G4Sphere::Clone() const
03010 {
03011 return new G4Sphere(*this);
03012 }
03013
03015
03016
03017
03018 std::ostream& G4Sphere::StreamInfo( std::ostream& os ) const
03019 {
03020 G4int oldprc = os.precision(16);
03021 os << "-----------------------------------------------------------\n"
03022 << " *** Dump for solid - " << GetName() << " ***\n"
03023 << " ===================================================\n"
03024 << " Solid type: G4Sphere\n"
03025 << " Parameters: \n"
03026 << " inner radius: " << fRmin/mm << " mm \n"
03027 << " outer radius: " << fRmax/mm << " mm \n"
03028 << " starting phi of segment : " << fSPhi/degree << " degrees \n"
03029 << " delta phi of segment : " << fDPhi/degree << " degrees \n"
03030 << " starting theta of segment: " << fSTheta/degree << " degrees \n"
03031 << " delta theta of segment : " << fDTheta/degree << " degrees \n"
03032 << "-----------------------------------------------------------\n";
03033 os.precision(oldprc);
03034
03035 return os;
03036 }
03037
03039
03040
03041
03042 G4ThreeVector G4Sphere::GetPointOnSurface() const
03043 {
03044 G4double zRand, aOne, aTwo, aThr, aFou, aFiv, chose, phi, sinphi, cosphi;
03045 G4double height1, height2, slant1, slant2, costheta, sintheta, rRand;
03046
03047 height1 = (fRmax-fRmin)*cosSTheta;
03048 height2 = (fRmax-fRmin)*cosETheta;
03049 slant1 = std::sqrt(sqr((fRmax - fRmin)*sinSTheta) + height1*height1);
03050 slant2 = std::sqrt(sqr((fRmax - fRmin)*sinETheta) + height2*height2);
03051 rRand = GetRadiusInRing(fRmin,fRmax);
03052
03053 aOne = fRmax*fRmax*fDPhi*(cosSTheta-cosETheta);
03054 aTwo = fRmin*fRmin*fDPhi*(cosSTheta-cosETheta);
03055 aThr = fDPhi*((fRmax + fRmin)*sinSTheta)*slant1;
03056 aFou = fDPhi*((fRmax + fRmin)*sinETheta)*slant2;
03057 aFiv = 0.5*fDTheta*(fRmax*fRmax-fRmin*fRmin);
03058
03059 phi = RandFlat::shoot(fSPhi, ePhi);
03060 cosphi = std::cos(phi);
03061 sinphi = std::sin(phi);
03062 costheta = RandFlat::shoot(cosETheta,cosSTheta);
03063 sintheta = std::sqrt(1.-sqr(costheta));
03064
03065 if(fFullPhiSphere) { aFiv = 0; }
03066 if(fSTheta == 0) { aThr=0; }
03067 if(eTheta == pi) { aFou = 0; }
03068 if(fSTheta == halfpi) { aThr = pi*(fRmax*fRmax-fRmin*fRmin); }
03069 if(eTheta == halfpi) { aFou = pi*(fRmax*fRmax-fRmin*fRmin); }
03070
03071 chose = RandFlat::shoot(0.,aOne+aTwo+aThr+aFou+2.*aFiv);
03072 if( (chose>=0.) && (chose<aOne) )
03073 {
03074 return G4ThreeVector(fRmax*sintheta*cosphi,
03075 fRmax*sintheta*sinphi, fRmax*costheta);
03076 }
03077 else if( (chose>=aOne) && (chose<aOne+aTwo) )
03078 {
03079 return G4ThreeVector(fRmin*sintheta*cosphi,
03080 fRmin*sintheta*sinphi, fRmin*costheta);
03081 }
03082 else if( (chose>=aOne+aTwo) && (chose<aOne+aTwo+aThr) )
03083 {
03084 if (fSTheta != halfpi)
03085 {
03086 zRand = RandFlat::shoot(fRmin*cosSTheta,fRmax*cosSTheta);
03087 return G4ThreeVector(tanSTheta*zRand*cosphi,
03088 tanSTheta*zRand*sinphi,zRand);
03089 }
03090 else
03091 {
03092 return G4ThreeVector(rRand*cosphi, rRand*sinphi, 0.);
03093 }
03094 }
03095 else if( (chose>=aOne+aTwo+aThr) && (chose<aOne+aTwo+aThr+aFou) )
03096 {
03097 if(eTheta != halfpi)
03098 {
03099 zRand = RandFlat::shoot(fRmin*cosETheta, fRmax*cosETheta);
03100 return G4ThreeVector (tanETheta*zRand*cosphi,
03101 tanETheta*zRand*sinphi,zRand);
03102 }
03103 else
03104 {
03105 return G4ThreeVector(rRand*cosphi, rRand*sinphi, 0.);
03106 }
03107 }
03108 else if( (chose>=aOne+aTwo+aThr+aFou) && (chose<aOne+aTwo+aThr+aFou+aFiv) )
03109 {
03110 return G4ThreeVector(rRand*sintheta*cosSPhi,
03111 rRand*sintheta*sinSPhi,rRand*costheta);
03112 }
03113 else
03114 {
03115 return G4ThreeVector(rRand*sintheta*cosEPhi,
03116 rRand*sintheta*sinEPhi,rRand*costheta);
03117 }
03118 }
03119
03121
03122
03123
03124 G4double G4Sphere::GetSurfaceArea()
03125 {
03126 if(fSurfaceArea != 0.) {;}
03127 else
03128 {
03129 G4double Rsq=fRmax*fRmax;
03130 G4double rsq=fRmin*fRmin;
03131
03132 fSurfaceArea = fDPhi*(rsq+Rsq)*(cosSTheta - cosETheta);
03133 if(!fFullPhiSphere)
03134 {
03135 fSurfaceArea = fSurfaceArea + fDTheta*(Rsq-rsq);
03136 }
03137 if(fSTheta >0)
03138 {
03139 G4double acos1=std::acos( std::pow(sinSTheta,2) * std::cos(fDPhi)
03140 + std::pow(cosSTheta,2));
03141 if(fDPhi>pi)
03142 {
03143 fSurfaceArea = fSurfaceArea + 0.5*(Rsq-rsq)*(twopi-acos1);
03144 }
03145 else
03146 {
03147 fSurfaceArea = fSurfaceArea + 0.5*(Rsq-rsq)*acos1;
03148 }
03149 }
03150 if(eTheta < pi)
03151 {
03152 G4double acos2=std::acos( std::pow(sinETheta,2) * std::cos(fDPhi)
03153 + std::pow(cosETheta,2));
03154 if(fDPhi>pi)
03155 {
03156 fSurfaceArea = fSurfaceArea + 0.5*(Rsq-rsq)*(twopi-acos2);
03157 }
03158 else
03159 {
03160 fSurfaceArea = fSurfaceArea + 0.5*(Rsq-rsq)*acos2;
03161 }
03162 }
03163 }
03164 return fSurfaceArea;
03165 }
03166
03168
03169
03170
03171 G4VisExtent G4Sphere::GetExtent() const
03172 {
03173 return G4VisExtent(-fRmax, fRmax,-fRmax, fRmax,-fRmax, fRmax );
03174 }
03175
03176
03177 void G4Sphere::DescribeYourselfTo ( G4VGraphicsScene& scene ) const
03178 {
03179 scene.AddSolid (*this);
03180 }
03181
03182 G4Polyhedron* G4Sphere::CreatePolyhedron () const
03183 {
03184 return new G4PolyhedronSphere (fRmin, fRmax, fSPhi, fDPhi, fSTheta, fDTheta);
03185 }
03186
03187 G4NURBS* G4Sphere::CreateNURBS () const
03188 {
03189 return new G4NURBSbox (fRmax, fRmax, fRmax);
03190 }