#include <G4Sphere.hh>
Inheritance diagram for G4Sphere:
Definition at line 77 of file G4Sphere.hh.
G4Sphere::G4Sphere | ( | const G4String & | pName, | |
G4double | pRmin, | |||
G4double | pRmax, | |||
G4double | pSPhi, | |||
G4double | pDPhi, | |||
G4double | pSTheta, | |||
G4double | pDTheta | |||
) |
Definition at line 90 of file G4Sphere.cc.
References FatalException, G4endl, G4Exception(), G4GeometryTolerance::GetAngularTolerance(), G4GeometryTolerance::GetInstance(), G4VSolid::GetName(), and G4GeometryTolerance::GetRadialTolerance().
Referenced by Clone().
00094 : G4CSGSolid(pName), fEpsilon(2.e-11), 00095 fFullPhiSphere(true), fFullThetaSphere(true) 00096 { 00097 kAngTolerance = G4GeometryTolerance::GetInstance()->GetAngularTolerance(); 00098 00099 // Check radii and set radial tolerances 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 // Check angles 00115 00116 CheckPhiAngles(pSPhi, pDPhi); 00117 CheckThetaAngles(pSTheta, pDTheta); 00118 }
G4Sphere::~G4Sphere | ( | ) |
G4Sphere::G4Sphere | ( | __void__ & | ) |
Definition at line 125 of file G4Sphere.cc.
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 }
G4Sphere::G4Sphere | ( | const G4Sphere & | rhs | ) |
Definition at line 149 of file G4Sphere.cc.
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 }
G4bool G4Sphere::CalculateExtent | ( | const EAxis | pAxis, | |
const G4VoxelLimits & | pVoxelLimit, | |||
const G4AffineTransform & | pTransform, | |||
G4double & | pmin, | |||
G4double & | pmax | |||
) | const [virtual] |
Implements G4VSolid.
Definition at line 220 of file G4Sphere.cc.
References G4VSolid::CalculateClippedPolygonExtent(), G4VoxelLimits::GetMaxExtent(), G4VoxelLimits::GetMaxXExtent(), G4VoxelLimits::GetMaxYExtent(), G4VoxelLimits::GetMaxZExtent(), G4VoxelLimits::GetMinExtent(), G4VoxelLimits::GetMinXExtent(), G4VoxelLimits::GetMinYExtent(), G4VoxelLimits::GetMinZExtent(), Inside(), G4AffineTransform::Inverse(), G4VoxelLimits::IsXLimited(), G4VoxelLimits::IsYLimited(), G4VoxelLimits::IsZLimited(), G4VSolid::kCarTolerance, kOutside, kXAxis, kYAxis, kZAxis, G4AffineTransform::NetTranslation(), and G4AffineTransform::TransformPoint().
00224 { 00225 if ( fFullSphere ) 00226 { 00227 // Special case handling for solid spheres-shells 00228 // (rotation doesn't influence). 00229 // Compute x/y/z mins and maxs for bounding box respecting limits, 00230 // with early returns if outside limits. Then switch() on pAxis, 00231 // and compute exact x and y limit for x/y case 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 // Known to cut sphere 00310 00311 switch (pAxis) 00312 { 00313 case kXAxis: 00314 yoff1=yoffset-yMin; 00315 yoff2=yMax-yoffset; 00316 if ((yoff1>=0) && (yoff2>=0)) 00317 { 00318 // Y limits cross max/min x => no change 00319 // 00320 pMin=xMin; 00321 pMax=xMax; 00322 } 00323 else 00324 { 00325 // Y limits don't cross max/min x => compute max delta x, 00326 // hence new mins/maxs 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 // X limits cross max/min y => no change 00345 // 00346 pMin=yMin; 00347 pMax=yMax; 00348 } 00349 else 00350 { 00351 // X limits don't cross max/min y => compute max delta y, 00352 // hence new mins/maxs 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 // Transformed cutted sphere 00378 { 00379 G4int i,j,noEntries,noBetweenSections; 00380 G4bool existsAfterClip=false; 00381 00382 // Calculate rotated vertex coordinates 00383 00384 G4ThreeVectorList* vertices; 00385 G4int noPolygonVertices ; 00386 vertices=CreateRotatedVertices(pTransform,noPolygonVertices); 00387 00388 pMin=+kInfinity; 00389 pMax=-kInfinity; 00390 00391 noEntries=vertices->size(); // noPolygonVertices*noPhiCrossSections 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 // Add 2*tolerance to avoid precision troubles 00431 // 00432 pMin-=kCarTolerance; 00433 pMax+=kCarTolerance; 00434 } 00435 else 00436 { 00437 // Check for case where completely enveloping clipping volume 00438 // If point inside then we are confident that the solid completely 00439 // envelopes the clipping volume. Hence set min/max extents according 00440 // to clipping volume extents along the specified axis. 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 }
G4VSolid * G4Sphere::Clone | ( | ) | const [virtual] |
Reimplemented from G4VSolid.
Definition at line 3009 of file G4Sphere.cc.
References G4Sphere().
03010 { 03011 return new G4Sphere(*this); 03012 }
void G4Sphere::ComputeDimensions | ( | G4VPVParameterisation * | p, | |
const G4int | n, | |||
const G4VPhysicalVolume * | pRep | |||
) | [virtual] |
Reimplemented from G4VSolid.
Definition at line 209 of file G4Sphere.cc.
References G4VPVParameterisation::ComputeDimensions().
00212 { 00213 p->ComputeDimensions(*this,n,pRep); 00214 }
G4NURBS * G4Sphere::CreateNURBS | ( | ) | const [virtual] |
Reimplemented from G4VSolid.
Definition at line 3187 of file G4Sphere.cc.
03188 { 03189 return new G4NURBSbox (fRmax, fRmax, fRmax); // Box for now!!! 03190 }
G4Polyhedron * G4Sphere::CreatePolyhedron | ( | ) | const [virtual] |
Reimplemented from G4VSolid.
Definition at line 3182 of file G4Sphere.cc.
03183 { 03184 return new G4PolyhedronSphere (fRmin, fRmax, fSPhi, fDPhi, fSTheta, fDTheta); 03185 }
void G4Sphere::DescribeYourselfTo | ( | G4VGraphicsScene & | scene | ) | const [virtual] |
Implements G4VSolid.
Definition at line 3177 of file G4Sphere.cc.
References G4VGraphicsScene::AddSolid().
03178 { 03179 scene.AddSolid (*this); 03180 }
G4double G4Sphere::DistanceToIn | ( | const G4ThreeVector & | p | ) | const [virtual] |
Implements G4VSolid.
Definition at line 1801 of file G4Sphere.cc.
References G4INCL::Math::pi.
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 // Distance to r shells 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 // Distance to phi extent 01834 // 01835 if (!fFullPhiSphere && rho) 01836 { 01837 // Psi=angle from central phi to point 01838 // 01839 cosPsi=(p.x()*cosCPhi+p.y()*sinCPhi)/rho; 01840 if (cosPsi<std::cos(hDPhi)) 01841 { 01842 // Point lies outside phi range 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 // Distance to Theta extent 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) // WHY ??????????? 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 }
G4double G4Sphere::DistanceToIn | ( | const G4ThreeVector & | p, | |
const G4ThreeVector & | v | |||
) | const [virtual] |
Implements G4VSolid.
Definition at line 867 of file G4Sphere.cc.
References G4VSolid::kCarTolerance, and G4INCL::Math::pi.
00869 { 00870 G4double snxt = kInfinity ; // snxt = default return value 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 // Intersection point 00889 // 00890 G4double xi, yi, zi, rhoi, rhoi2, radi2, iTheta ; 00891 00892 // Phi intersection 00893 // 00894 G4double Comp ; 00895 00896 // Phi precalcs 00897 // 00898 G4double Dist, cosPsi ; 00899 00900 // Theta precalcs 00901 // 00902 G4double dist2STheta, dist2ETheta ; 00903 G4double t1, t2, b, c, d2, d, sd = kInfinity ; 00904 00905 // General Precalcs 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 // Theta precalcs 00915 // 00916 if (!fFullThetaSphere) 00917 { 00918 tolSTheta = fSTheta - halfAngTolerance ; 00919 tolETheta = eTheta + halfAngTolerance ; 00920 } 00921 00922 // Outer spherical shell intersection 00923 // - Only if outside tolerant fRmax 00924 // - Check for if inside and outer G4Sphere heading through solid (-> 0) 00925 // - No intersect -> no intersection with G4Sphere 00926 // 00927 // Shell eqn: x^2+y^2+z^2=RSPH^2 00928 // 00929 // => (px+svx)^2+(py+svy)^2+(pz+svz)^2=R^2 00930 // 00931 // => (px^2+py^2+pz^2) +2sd(pxvx+pyvy+pzvz)+sd^2(vx^2+vy^2+vz^2)=R^2 00932 // => rad2 +2sd(pDotV3d) +sd^2 =R^2 00933 // 00934 // => sd=-pDotV3d+-std::sqrt(pDotV3d^2-(rad2-R^2)) 00935 00936 c = rad2 - fRmax*fRmax ; 00937 00938 if (c > fRmaxTolerance*fRmax) 00939 { 00940 // If outside tolerant boundary of outer G4Sphere 00941 // [should be std::sqrt(rad2)-fRmax > halfRmaxTolerance] 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 ) // Avoid rounding errors due to precision issues seen on 00952 { // 64 bits systems. Split long distances and recompute 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) // Check phi intersection 00961 { 00962 cosPsi = (xi*cosCPhi + yi*sinCPhi)/rhoi ; 00963 00964 if (cosPsi >= cosHDPhiOT) 00965 { 00966 if (!fFullThetaSphere) // Check theta intersection 00967 { 00968 zi = p.z() + sd*v.z() ; 00969 00970 // rhoi & zi can never both be 0 00971 // (=>intersect at origin =>fRmax=0) 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) // Check theta intersection 00988 { 00989 zi = p.z() + sd*v.z() ; 00990 00991 // rhoi & zi can never both be 0 00992 // (=>intersect at origin => fRmax=0 !) 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 // No intersection with G4Sphere 01008 { 01009 return snxt=kInfinity; 01010 } 01011 } 01012 else 01013 { 01014 // Inside outer radius 01015 // check not inside, and heading through G4Sphere (-> 0 to in) 01016 01017 d2 = pDotV3d*pDotV3d - c ; 01018 01019 if ( (rad2 > tolIRMax2) 01020 && ( (d2 >= fRmaxTolerance*fRmax) && (pDotV3d < 0) ) ) 01021 { 01022 if (!fFullPhiSphere) 01023 { 01024 // Use inner phi tolerant boundary -> if on tolerant 01025 // phi boundaries, phi intersect code handles leaving/entering checks 01026 01027 cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(rho2) ; 01028 01029 if (cosPsi>=cosHDPhiIT) 01030 { 01031 // inside radii, delta r -ve, inside phi 01032 01033 if ( !fFullThetaSphere ) 01034 { 01035 if ( (pTheta >= tolSTheta + kAngTolerance) 01036 && (pTheta <= tolETheta - kAngTolerance) ) 01037 { 01038 return snxt=0; 01039 } 01040 } 01041 else // strictly inside Theta in both cases 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 // strictly inside Theta in both cases 01058 { 01059 return snxt=0; 01060 } 01061 } 01062 } 01063 } 01064 01065 // Inner spherical shell intersection 01066 // - Always farthest root, because would have passed through outer 01067 // surface first. 01068 // - Tolerant check if travelling through solid 01069 01070 if (fRmin) 01071 { 01072 c = rad2 - fRmin*fRmin ; 01073 d2 = pDotV3d*pDotV3d - c ; 01074 01075 // Within tolerance inner radius of inner G4Sphere 01076 // Check for immediate entry/already inside and travelling outwards 01077 01078 if ( (c > -halfRminTolerance) && (rad2 < tolIRMin2) 01079 && ( (d2 < fRmin*kCarTolerance) || (pDotV3d >= 0) ) ) 01080 { 01081 if ( !fFullPhiSphere ) 01082 { 01083 // Use inner phi tolerant boundary -> if on tolerant 01084 // phi boundaries, phi intersect code handles leaving/entering checks 01085 01086 cosPsi = (p.x()*cosCPhi+p.y()*sinCPhi)/std::sqrt(rho2) ; 01087 if (cosPsi >= cosHDPhiIT) 01088 { 01089 // inside radii, delta r -ve, inside phi 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 // Not special tolerant case 01122 { 01123 if (d2 >= 0) 01124 { 01125 sd = -pDotV3d + std::sqrt(d2) ; 01126 if ( sd >= halfRminTolerance ) // It was >= 0 ?? 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 ) // Check phi intersection 01133 { 01134 cosPsi = (xi*cosCPhi + yi*sinCPhi)/rhoi ; 01135 01136 if (cosPsi >= cosHDPhiOT) 01137 { 01138 if ( !fFullThetaSphere ) // Check theta intersection 01139 { 01140 zi = p.z() + sd*v.z() ; 01141 01142 // rhoi & zi can never both be 0 01143 // (=>intersect at origin =>fRmax=0) 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 ) // Check theta intersection 01160 { 01161 zi = p.z() + sd*v.z() ; 01162 01163 // rhoi & zi can never both be 0 01164 // (=>intersect at origin => fRmax=0 !) 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 // Phi segment intersection 01183 // 01184 // o Tolerant of points inside phi planes by up to kCarTolerance*0.5 01185 // 01186 // o NOTE: Large duplication of code between sphi & ephi checks 01187 // -> only diffs: sphi -> ephi, Comp -> -Comp and half-plane 01188 // intersection check <=0 -> >=0 01189 // -> Should use some form of loop Construct 01190 // 01191 if ( !fFullPhiSphere ) 01192 { 01193 // First phi surface ('S'tarting phi) 01194 // Comp = Component in outwards normal dirn 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 // Check theta intersection 01230 // rhoi & zi can never both be 0 01231 // (=>intersect at origin =>fRmax=0) 01232 // 01233 if ( !fFullThetaSphere ) 01234 { 01235 iTheta = std::atan2(std::sqrt(rhoi2),zi) ; 01236 if ( (iTheta >= tolSTheta) && (iTheta <= tolETheta) ) 01237 { 01238 // r and theta intersections good 01239 // - check intersecting with correct half-plane 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 // Second phi surface ('E'nding phi) 01257 // Component in outwards normal dirn 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 // Check theta intersection 01292 // rhoi & zi can never both be 0 01293 // (=>intersect at origin =>fRmax=0) 01294 // 01295 if ( !fFullThetaSphere ) 01296 { 01297 iTheta = std::atan2(std::sqrt(rhoi2),zi) ; 01298 if ( (iTheta >= tolSTheta) && (iTheta <= tolETheta) ) 01299 { 01300 // r and theta intersections good 01301 // - check intersecting with correct half-plane 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 // Theta segment intersection 01320 01321 if ( !fFullThetaSphere ) 01322 { 01323 01324 // Intersection with theta surfaces 01325 // Known failure cases: 01326 // o Inside tolerance of stheta surface, skim 01327 // ~parallel to cone and Hit & enter etheta surface [& visa versa] 01328 // 01329 // To solve: Check 2nd root of etheta surface in addition to stheta 01330 // 01331 // o start/end theta is exactly pi/2 01332 // Intersections with cones 01333 // 01334 // Cone equation: x^2+y^2=z^2tan^2(t) 01335 // 01336 // => (px+svx)^2+(py+svy)^2=(pz+svz)^2tan^2(t) 01337 // 01338 // => (px^2+py^2-pz^2tan^2(t))+2sd(pxvx+pyvy-pzvztan^2(t)) 01339 // + sd^2(vx^2+vy^2-vz^2tan^2(t)) = 0 01340 // 01341 // => sd^2(1-vz^2(1+tan^2(t))+2sd(pdotv2d-pzvztan^2(t))+(rho2-pz^2tan^2(t))=0 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 // Inside (theta<stheta-tol) stheta cone 01362 // First root of stheta cone, second if first root -ve 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 ; // First root 01376 zi = p.z() + sd*v.z(); 01377 01378 if ( (sd < 0) || (zi*(fSTheta - halfpi) > 0) ) 01379 { 01380 sd = -b+d; // Second root 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 ) // Check phi intersection 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 // Possible intersection with ETheta cone. 01411 // Second >= 0 root should be considered 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 ; // Second root 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) // Check phi intersection 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 // dist2ETheta<-kRadTolerance*0.5 && dist2STheta>0) 01461 // Inside (theta > etheta+tol) e-theta cone 01462 // First root of etheta cone, second if first root 'imaginary' 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 ; // First root 01476 zi = p.z() + sd*v.z(); 01477 01478 if ( (sd < 0) || (zi*(eTheta - halfpi) > 0) ) 01479 { 01480 sd = -b + d ; // second root 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) // Check phi intersection 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 // Possible intersection with STheta cone. 01512 // Second >= 0 root should be considered 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 ; // Second root 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) // Check phi intersection 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 // In tolerance of stheta 01563 // If entering through solid [r,phi] => 0 to in 01564 // else try 2nd root 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) // Check phi intersection 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 // Not entering immediately/travelling through 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 { // ^^^^^^^^^^^^^^^^^^^^^ shouldn't it be >=0 instead ? 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 ) // Check phi intersection 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 // In tolerance of etheta 01631 // If entering through solid [r,phi] => 0 to in 01632 // else try 2nd root 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) // Check phi intersection 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 // Not entering immediately/travelling through 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) // Check phi intersection 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 // stheta+tol<theta<etheta-tol 01704 // For BOTH stheta & etheta check 2nd root for validity [r,phi] 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 ; // second root 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) // Check phi intersection 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; // second root 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) // Check phi intersection 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 }
G4double G4Sphere::DistanceToOut | ( | const G4ThreeVector & | p | ) | const [virtual] |
Implements G4VSolid.
Definition at line 2781 of file G4Sphere.cc.
References G4VSolid::DumpInfo(), G4cout, G4endl, G4Exception(), Inside(), JustWarning, kOutside, and G4INCL::Math::pi.
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 // Distance to r shells 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 // Distance to phi extent 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 // Distance to Theta extent 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 }
G4double G4Sphere::DistanceToOut | ( | const G4ThreeVector & | p, | |
const G4ThreeVector & | v, | |||
const G4bool | calcNorm = G4bool(false) , |
|||
G4bool * | validNorm = 0 , |
|||
G4ThreeVector * | n = 0 | |||
) | const [virtual] |
Implements G4VSolid.
Definition at line 1897 of file G4Sphere.cc.
References G4VSolid::DumpInfo(), G4cout, G4endl, G4Exception(), JustWarning, G4VSolid::kCarTolerance, and G4INCL::Math::pi.
01902 { 01903 G4double snxt = kInfinity; // snxt is default return value 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 // Variables for phi intersection: 01917 01918 G4double pDistS,compS,pDistE,compE,sphi2,vphi; 01919 01920 G4double rho2,rad2,pDotV2d,pDotV3d; 01921 01922 G4double xi,yi,zi; // Intersection point 01923 01924 // Theta precals 01925 // 01926 G4double rhoSecTheta; 01927 G4double dist2STheta, dist2ETheta, distTheta; 01928 G4double d2,sd; 01929 01930 // General Precalcs 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 // Radial Intersections from G4Sphere::DistanceToIn 01939 // 01940 // Outer spherical shell intersection 01941 // - Only if outside tolerant fRmax 01942 // - Check for if inside and outer G4Sphere heading through solid (-> 0) 01943 // - No intersect -> no intersection with G4Sphere 01944 // 01945 // Shell eqn: x^2+y^2+z^2=RSPH^2 01946 // 01947 // => (px+svx)^2+(py+svy)^2+(pz+svz)^2=R^2 01948 // 01949 // => (px^2+py^2+pz^2) +2sd(pxvx+pyvy+pzvz)+sd^2(vx^2+vy^2+vz^2)=R^2 01950 // => rad2 +2sd(pDotV3d) +sd^2 =R^2 01951 // 01952 // => sd=-pDotV3d+-std::sqrt(pDotV3d^2-(rad2-R^2)) 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 // Within tolerant Outer radius 01961 // 01962 // The test is 01963 // rad - fRmax < 0.5*kRadTolerance 01964 // => rad < fRmax + 0.5*kRadTol 01965 // => rad2 < (fRmax + 0.5*kRadTol)^2 01966 // => rad2 < fRmax^2 + 2.*0.5*fRmax*kRadTol + 0.25*kRadTol*kRadTol 01967 // => rad2 - fRmax^2 <~ fRmax*kRadTol 01968 01969 d2 = pDotV3d*pDotV3d - c; 01970 01971 if( (c >- fRmaxTolerance*fRmax) // on tolerant surface 01972 && ((pDotV3d >=0) || (d2 < 0)) ) // leaving outside from Rmax 01973 // not re-entering 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); // second root since inside Rmax 01985 side = kRMax ; 01986 } 01987 } 01988 01989 // Inner spherical shell intersection: 01990 // Always first >=0 root, because would have passed 01991 // from outside of Rmin surface . 01992 01993 if (fRmin) 01994 { 01995 c = rad2 - fRmin*fRmin; 01996 d2 = pDotV3d*pDotV3d - c; 01997 01998 if (c >- fRminTolerance*fRmin) // 2.0 * (0.5*kRadTolerance) * fRmin 01999 { 02000 if ( (c < fRminTolerance*fRmin) // leaving from Rmin 02001 && (d2 >= fRminTolerance*fRmin) && (pDotV3d < 0) ) 02002 { 02003 if(calcNorm) { *validNorm = false; } // Rmin surface is concave 02004 return snxt = 0 ; 02005 } 02006 else 02007 { 02008 if ( d2 >= 0. ) 02009 { 02010 sd = -pDotV3d-std::sqrt(d2); 02011 02012 if ( sd >= 0. ) // Always intersect Rmin first 02013 { 02014 snxt = sd ; 02015 side = kRMin ; 02016 } 02017 } 02018 } 02019 } 02020 } 02021 } 02022 02023 // Theta segment intersection 02024 02025 if ( !fFullThetaSphere ) 02026 { 02027 // Intersection with theta surfaces 02028 // 02029 // Known failure cases: 02030 // o Inside tolerance of stheta surface, skim 02031 // ~parallel to cone and Hit & enter etheta surface [& visa versa] 02032 // 02033 // To solve: Check 2nd root of etheta surface in addition to stheta 02034 // 02035 // o start/end theta is exactly pi/2 02036 // 02037 // Intersections with cones 02038 // 02039 // Cone equation: x^2+y^2=z^2tan^2(t) 02040 // 02041 // => (px+svx)^2+(py+svy)^2=(pz+svz)^2tan^2(t) 02042 // 02043 // => (px^2+py^2-pz^2tan^2(t))+2sd(pxvx+pyvy-pzvztan^2(t)) 02044 // + sd^2(vx^2+vy^2-vz^2tan^2(t)) = 0 02045 // 02046 // => sd^2(1-vz^2(1+tan^2(t))+2sd(pdotv2d-pzvztan^2(t))+(rho2-pz^2tan^2(t))=0 02047 // 02048 02049 if(fSTheta) // intersection with first cons 02050 { 02051 if( std::fabs(tanSTheta) > 5./kAngTolerance ) // kons is plane z=0 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 // kons is not plane 02069 { 02070 t1 = 1-v.z()*v.z()*(1+tanSTheta2); 02071 t2 = pDotV2d-p.z()*v.z()*tanSTheta2; // ~vDotN if p on cons 02072 dist2STheta = rho2-p.z()*p.z()*tanSTheta2; // t3 02073 02074 distTheta = std::sqrt(rho2)-p.z()*tanSTheta; 02075 02076 if( std::fabs(t1) < halfAngTolerance ) // 1st order equation, 02077 { // v parallel to kons 02078 if( v.z() > 0. ) 02079 { 02080 if(std::fabs(distTheta) < halfRmaxTolerance) // p on surface 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 } // 2nd order equation, 1st root of fSTheta cone, 02109 else // 2nd if 1st root -ve 02110 { 02111 if( std::fabs(distTheta) < halfRmaxTolerance ) 02112 { 02113 if( (fSTheta > halfpi) && (t2 >= 0.) ) // leave 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.) ) // leave 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; // First root 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 ; // 2nd root 02152 } 02153 if( (sd > halfRmaxTolerance) && (p.z() + sd*v.z() <= 0.) ) 02154 { 02155 stheta = sd; 02156 sidetheta = kSTheta; 02157 } 02158 } 02159 else // sTheta < pi/2, concave surface, no normal 02160 { 02161 sd = -b - d; // First root 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 ; // 2nd root 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) // intersection with second cons 02179 { 02180 if( std::fabs(tanETheta) > 5./kAngTolerance ) // kons is plane z=0 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 // kons is not plane 02203 { 02204 t1 = 1-v.z()*v.z()*(1+tanETheta2); 02205 t2 = pDotV2d-p.z()*v.z()*tanETheta2; // ~vDotN if p on cons 02206 dist2ETheta = rho2-p.z()*p.z()*tanETheta2; // t3 02207 02208 distTheta = std::sqrt(rho2)-p.z()*tanETheta; 02209 02210 if( std::fabs(t1) < halfAngTolerance ) // 1st order equation, 02211 { // v parallel to kons 02212 if( v.z() < 0. ) 02213 { 02214 if(std::fabs(distTheta) < halfRmaxTolerance) // p on surface 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 } // 2nd order equation, 1st root of fSTheta cone 02247 else // 2nd if 1st root -ve 02248 { 02249 if ( std::fabs(distTheta) < halfRmaxTolerance ) 02250 { 02251 if( (eTheta < halfpi) && (t2 >= 0.) ) // leave 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.) ) // leave 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; // First root 02285 02286 if( ((std::fabs(sd) < halfRmaxTolerance) && (t2 < 0.)) 02287 || (sd < 0.) ) 02288 { 02289 sd = -b + d ; // 2nd root 02290 } 02291 if( sd > halfRmaxTolerance ) 02292 { 02293 if( sd < stheta ) 02294 { 02295 stheta = sd; 02296 sidetheta = kETheta; 02297 } 02298 } 02299 } 02300 else // sTheta+fDTheta > pi/2, concave surface, no normal 02301 { 02302 sd = -b - d; // First root 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 ; // 2nd root 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 } // end theta intersections 02324 02325 // Phi Intersection 02326 02327 if ( !fFullPhiSphere ) 02328 { 02329 if ( p.x() || p.y() ) // Check if on z axis (rho not needed later) 02330 { 02331 // pDist -ve when inside 02332 02333 pDistS=p.x()*sinSPhi-p.y()*cosSPhi; 02334 pDistE=-p.x()*sinEPhi+p.y()*cosEPhi; 02335 02336 // Comp -ve when in direction of outwards normal 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 // Inside both phi *full* planes 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 // Check intersection with correct half-plane (if not -> no intersect) 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; } // Leave by sphi 02372 } 02373 } 02374 else { sphi = kInfinity; } 02375 02376 if ( compE < 0 ) 02377 { 02378 sphi2=pDistE/compE ; 02379 if (sphi2 < sphi) // Only check further if < starting phi intersection 02380 { 02381 xi = p.x()+sphi2*v.x() ; 02382 yi = p.y()+sphi2*v.y() ; 02383 02384 // Check intersection with correct half-plane 02385 // 02386 if ((std::fabs(xi)<=kCarTolerance) && (std::fabs(yi)<=kCarTolerance)) 02387 { 02388 // Leaving via ending phi 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) // Leaving via ending phi 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)) // Outside both *full* phi planes 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 // if towards both >=0 then once inside (after error) 02433 // will remain inside 02434 02435 if ( (compS >= 0) && (compE >= 0) ) { sphi = kInfinity; } 02436 else { sphi = 0; } 02437 } 02438 } 02439 else if ( (pDistS > 0) && (pDistE < 0) ) 02440 { 02441 // Outside full starting plane, inside full ending plane 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 // Check intersection in correct half-plane 02452 // (if not -> not leaving phi extent) 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 // Leaving via Ending phi 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 // Check intersection in correct half-plane 02490 // (if not -> remain in extent) 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 // otherwise leaving via Ending phi 02507 { 02508 sidephi = kEPhi ; 02509 } 02510 } 02511 else sphi=kInfinity; 02512 } 02513 else // leaving immediately by starting phi 02514 { 02515 sidephi = kSPhi ; 02516 sphi = 0 ; 02517 } 02518 } 02519 } 02520 else 02521 { 02522 // Must be pDistS < 0 && pDistE > 0 02523 // Inside full starting plane, outside full ending plane 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 // Check intersection in correct half-plane 02534 // (if not -> not leaving phi extent) 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 // Leaving via Starting phi 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 // Check intersection in correct half-plane 02572 // (if not -> remain in extent) 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 // otherwise leaving via Starting phi 02589 { 02590 sidephi = kSPhi ; 02591 } 02592 } 02593 else 02594 { 02595 sphi = kInfinity ; 02596 } 02597 } 02598 else // leaving immediately by ending 02599 { 02600 sidephi = kEPhi ; 02601 sphi = 0 ; 02602 } 02603 } 02604 } 02605 } 02606 else 02607 { 02608 // On z axis + travel not || to z axis -> if phi of vector direction 02609 // within phi of shape, Step limited by rmax, else Step =0 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 ; // arbitrary 02621 sphi = 0 ; 02622 } 02623 } 02624 else // travel along z - no phi intersection 02625 { 02626 sphi = kInfinity ; 02627 } 02628 } 02629 if ( sphi < snxt ) // Order intersecttions 02630 { 02631 snxt = sphi ; 02632 side = sidephi ; 02633 } 02634 } 02635 if (stheta < snxt ) // Order intersections 02636 { 02637 snxt = stheta ; 02638 side = sidetheta ; 02639 } 02640 02641 if (calcNorm) // Output switch operator 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; // Rmin is concave 02655 break; 02656 02657 case kSPhi: 02658 if ( fDPhi <= pi ) // Normal to Phi- 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 ) // Normal to Phi+ 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; } // Concave STheta cone 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; } // Concave ETheta cone 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 }
G4double G4Sphere::GetCubicVolume | ( | ) | [inline, virtual] |
Reimplemented from G4VSolid.
Definition at line 309 of file G4Sphere.icc.
References G4CSGSolid::fCubicVolume.
00310 { 00311 if(fCubicVolume != 0.) {;} 00312 else { fCubicVolume = fDPhi*(std::cos(fSTheta)-std::cos(fSTheta+fDTheta))* 00313 (fRmax*fRmax*fRmax-fRmin*fRmin*fRmin)/3.; } 00314 return fCubicVolume; 00315 }
G4double G4Sphere::GetDeltaPhiAngle | ( | ) | const [inline] |
Definition at line 62 of file G4Sphere.icc.
Referenced by GetDPhi(), G4tgbGeometryDumper::GetSolidParams(), G4PSSphereSurfaceFlux::ProcessHits(), G4PSSphereSurfaceCurrent::ProcessHits(), G4GDMLWriteParamvol::Sphere_dimensionsWrite(), and G4GDMLWriteSolids::SphereWrite().
G4double G4Sphere::GetDeltaThetaAngle | ( | ) | const [inline] |
Definition at line 73 of file G4Sphere.icc.
Referenced by GetDTheta(), G4tgbGeometryDumper::GetSolidParams(), G4PSSphereSurfaceFlux::ProcessHits(), G4PSSphereSurfaceCurrent::ProcessHits(), G4GDMLWriteParamvol::Sphere_dimensionsWrite(), and G4GDMLWriteSolids::SphereWrite().
G4double G4Sphere::GetDPhi | ( | ) | const [inline] |
Definition at line 291 of file G4Sphere.icc.
References GetDeltaPhiAngle().
00292 { 00293 return GetDeltaPhiAngle(); 00294 }
G4double G4Sphere::GetDTheta | ( | ) | const [inline] |
Definition at line 303 of file G4Sphere.icc.
References GetDeltaThetaAngle().
00304 { 00305 return GetDeltaThetaAngle(); 00306 }
G4GeometryType G4Sphere::GetEntityType | ( | ) | const [virtual] |
Implements G4VSolid.
Definition at line 3000 of file G4Sphere.cc.
03001 { 03002 return G4String("G4Sphere"); 03003 }
G4VisExtent G4Sphere::GetExtent | ( | ) | const [virtual] |
Reimplemented from G4VSolid.
Definition at line 3171 of file G4Sphere.cc.
03172 { 03173 return G4VisExtent(-fRmax, fRmax,-fRmax, fRmax,-fRmax, fRmax ); 03174 }
G4double G4Sphere::GetInnerRadius | ( | ) | const [inline] |
G4double G4Sphere::GetInsideRadius | ( | ) | const [inline] |
Definition at line 38 of file G4Sphere.icc.
Referenced by GetRmin(), G4PSSphereSurfaceFlux::IsSelectedSurface(), G4PSSphereSurfaceCurrent::IsSelectedSurface(), G4PSSphereSurfaceFlux::ProcessHits(), G4PSSphereSurfaceCurrent::ProcessHits(), G4GDMLWriteParamvol::Sphere_dimensionsWrite(), and G4GDMLWriteSolids::SphereWrite().
G4double G4Sphere::GetOuterRadius | ( | ) | const [inline] |
Definition at line 50 of file G4Sphere.icc.
Referenced by GetRmax(), G4tgbGeometryDumper::GetSolidParams(), G4GDMLWriteParamvol::Sphere_dimensionsWrite(), and G4GDMLWriteSolids::SphereWrite().
G4ThreeVector G4Sphere::GetPointOnSurface | ( | ) | const [virtual] |
Reimplemented from G4VSolid.
Definition at line 3042 of file G4Sphere.cc.
References G4CSGSolid::GetRadiusInRing(), G4INCL::Math::pi, and sqr().
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 }
G4double G4Sphere::GetRmax | ( | ) | const [inline] |
Definition at line 279 of file G4Sphere.icc.
References GetOuterRadius().
00280 { 00281 return GetOuterRadius(); 00282 }
G4double G4Sphere::GetRmin | ( | ) | const [inline] |
Definition at line 273 of file G4Sphere.icc.
References GetInsideRadius().
00274 { 00275 return GetInsideRadius(); 00276 }
G4double G4Sphere::GetSPhi | ( | ) | const [inline] |
Definition at line 285 of file G4Sphere.icc.
References GetStartPhiAngle().
00286 { 00287 return GetStartPhiAngle(); 00288 }
G4double G4Sphere::GetStartPhiAngle | ( | ) | const [inline] |
Definition at line 56 of file G4Sphere.icc.
Referenced by G4tgbGeometryDumper::GetSolidParams(), GetSPhi(), G4GDMLWriteParamvol::Sphere_dimensionsWrite(), and G4GDMLWriteSolids::SphereWrite().
G4double G4Sphere::GetStartThetaAngle | ( | ) | const [inline] |
Definition at line 68 of file G4Sphere.icc.
Referenced by G4tgbGeometryDumper::GetSolidParams(), GetSTheta(), G4PSSphereSurfaceFlux::ProcessHits(), G4PSSphereSurfaceCurrent::ProcessHits(), G4GDMLWriteParamvol::Sphere_dimensionsWrite(), and G4GDMLWriteSolids::SphereWrite().
G4double G4Sphere::GetSTheta | ( | ) | const [inline] |
Definition at line 297 of file G4Sphere.icc.
References GetStartThetaAngle().
00298 { 00299 return GetStartThetaAngle(); 00300 }
G4double G4Sphere::GetSurfaceArea | ( | ) | [virtual] |
Reimplemented from G4VSolid.
Definition at line 3124 of file G4Sphere.cc.
References G4CSGSolid::fSurfaceArea, and G4INCL::Math::pi.
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 }
EInside G4Sphere::Inside | ( | const G4ThreeVector & | p | ) | const [virtual] |
Implements G4VSolid.
Definition at line 465 of file G4Sphere.cc.
References kInside, kOutside, and kSurface.
Referenced by CalculateExtent(), and DistanceToOut().
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 // Check radial surfaces. Sets 'in' 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; // outside case 00491 tolRMin = std::max(fRmin-halfRminTolerance, 0.); // outside case 00492 if ( (rad2 <= tolRMax*tolRMax) && (rad2 >= tolRMin*tolRMin) ) 00493 { 00494 in = kSurface; 00495 } 00496 else 00497 { 00498 return in = kOutside; 00499 } 00500 } 00501 00502 // Phi boundaries : Do not check if it has no phi boundary! 00503 00504 if ( !fFullPhiSphere && rho2 ) // [fDPhi < twopi] and [p.x or p.y] 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) // else it's kSurface anyway already 00515 { 00516 if ( (pPhi < fSPhi + halfAngTolerance) 00517 || (pPhi > ePhi - halfAngTolerance) ) { in = kSurface; } 00518 } 00519 } 00520 00521 // Theta bondaries 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 }
Definition at line 173 of file G4Sphere.cc.
References cosCPhi, cosEPhi, cosETheta, cosHDPhiIT, cosHDPhiOT, cosSPhi, cosSTheta, cPhi, ePhi, eTheta, fDPhi, fDTheta, fEpsilon, fFullPhiSphere, fFullSphere, fFullThetaSphere, fRmax, fRmaxTolerance, fRmin, fRminTolerance, fSPhi, fSTheta, hDPhi, kAngTolerance, kRadTolerance, G4CSGSolid::operator=(), sinCPhi, sinEPhi, sinETheta, sinSPhi, sinSTheta, tanETheta, tanETheta2, tanSTheta, and tanSTheta2.
00174 { 00175 // Check assignment to self 00176 // 00177 if (this == &rhs) { return *this; } 00178 00179 // Copy base class data 00180 // 00181 G4CSGSolid::operator=(rhs); 00182 00183 // Copy data 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 }
void G4Sphere::SetDeltaPhiAngle | ( | G4double | newDphi | ) | [inline] |
void G4Sphere::SetDeltaThetaAngle | ( | G4double | newDTheta | ) | [inline] |
void G4Sphere::SetInnerRadius | ( | G4double | newRMin | ) | [inline] |
Definition at line 224 of file G4Sphere.icc.
References SetInsideRadius().
00225 { 00226 SetInsideRadius(newRmin); 00227 }
void G4Sphere::SetInsideRadius | ( | G4double | newRmin | ) | [inline] |
Definition at line 216 of file G4Sphere.icc.
Referenced by SetInnerRadius().
00217 { 00218 fRmin= newRmin; 00219 fRminTolerance = (fRmin) ? std::max( kRadTolerance, fEpsilon*fRmin ) : 0; 00220 Initialize(); 00221 }
void G4Sphere::SetOuterRadius | ( | G4double | newRmax | ) | [inline] |
Definition at line 230 of file G4Sphere.icc.
00231 { 00232 fRmax= newRmax; 00233 fRmaxTolerance = std::max( kRadTolerance, fEpsilon*fRmax ); 00234 Initialize(); 00235 }
Definition at line 238 of file G4Sphere.icc.
00239 { 00240 // Flag 'compute' can be used to explicitely avoid recomputation of 00241 // trigonometry in case SetDeltaPhiAngle() is invoked afterwards 00242 00243 CheckSPhiAngle(newSPhi); 00244 fFullPhiSphere = false; 00245 if (compute) { InitializePhiTrigonometry(); } 00246 Initialize(); 00247 }
void G4Sphere::SetStartThetaAngle | ( | G4double | newSTheta | ) | [inline] |
std::ostream & G4Sphere::StreamInfo | ( | std::ostream & | os | ) | const [virtual] |
Reimplemented from G4CSGSolid.
Definition at line 3018 of file G4Sphere.cc.
References G4VSolid::GetName().
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 }
G4ThreeVector G4Sphere::SurfaceNormal | ( | const G4ThreeVector & | p | ) | const [virtual] |
Implements G4VSolid.
Definition at line 562 of file G4Sphere.cc.
References G4Exception(), JustWarning, G4VSolid::kCarTolerance, and G4INCL::Math::pi.
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 }