#include <G4Cons.hh>
Inheritance diagram for G4Cons:
Definition at line 74 of file G4Cons.hh.
G4Cons::G4Cons | ( | const G4String & | pName, | |
G4double | pRmin1, | |||
G4double | pRmax1, | |||
G4double | pRmin2, | |||
G4double | pRmax2, | |||
G4double | pDz, | |||
G4double | pSPhi, | |||
G4double | pDPhi | |||
) |
Definition at line 79 of file G4Cons.cc.
References FatalException, G4endl, G4Exception(), G4GeometryTolerance::GetAngularTolerance(), G4GeometryTolerance::GetInstance(), G4VSolid::GetName(), and G4GeometryTolerance::GetRadialTolerance().
Referenced by Clone().
00084 : G4CSGSolid(pName), fRmin1(pRmin1), fRmin2(pRmin2), 00085 fRmax1(pRmax1), fRmax2(pRmax2), fDz(pDz), fSPhi(0.), fDPhi(0.) 00086 { 00087 kRadTolerance = G4GeometryTolerance::GetInstance()->GetRadialTolerance(); 00088 kAngTolerance = G4GeometryTolerance::GetInstance()->GetAngularTolerance(); 00089 00090 // Check z-len 00091 // 00092 if ( pDz < 0 ) 00093 { 00094 std::ostringstream message; 00095 message << "Invalid Z half-length for Solid: " << GetName() << G4endl 00096 << " hZ = " << pDz; 00097 G4Exception("G4Cons::G4Cons()", "GeomSolids0002", 00098 FatalException, message); 00099 } 00100 00101 // Check radii 00102 // 00103 if (((pRmin1>=pRmax1) || (pRmin2>=pRmax2) || (pRmin1<0)) && (pRmin2<0)) 00104 { 00105 std::ostringstream message; 00106 message << "Invalid values of radii for Solid: " << GetName() << G4endl 00107 << " pRmin1 = " << pRmin1 << ", pRmin2 = " << pRmin2 00108 << ", pRmax1 = " << pRmax1 << ", pRmax2 = " << pRmax2; 00109 G4Exception("G4Cons::G4Cons()", "GeomSolids0002", 00110 FatalException, message) ; 00111 } 00112 if( (pRmin1 == 0.0) && (pRmin2 > 0.0) ) { fRmin1 = 1e3*kRadTolerance ; } 00113 if( (pRmin2 == 0.0) && (pRmin1 > 0.0) ) { fRmin2 = 1e3*kRadTolerance ; } 00114 00115 // Check angles 00116 // 00117 CheckPhiAngles(pSPhi, pDPhi); 00118 }
G4Cons::G4Cons | ( | __void__ & | ) |
Definition at line 125 of file G4Cons.cc.
00126 : G4CSGSolid(a), kRadTolerance(0.), kAngTolerance(0.), 00127 fRmin1(0.), fRmin2(0.), fRmax1(0.), fRmax2(0.), fDz(0.), 00128 fSPhi(0.), fDPhi(0.), sinCPhi(0.), cosCPhi(0.), cosHDPhiOT(0.), 00129 cosHDPhiIT(0.), sinSPhi(0.), cosSPhi(0.), sinEPhi(0.), cosEPhi(0.), 00130 fPhiFullCone(false) 00131 { 00132 }
G4Cons::G4Cons | ( | const G4Cons & | rhs | ) |
Definition at line 146 of file G4Cons.cc.
00147 : G4CSGSolid(rhs), kRadTolerance(rhs.kRadTolerance), 00148 kAngTolerance(rhs.kAngTolerance), fRmin1(rhs.fRmin1), fRmin2(rhs.fRmin2), 00149 fRmax1(rhs.fRmax1), fRmax2(rhs.fRmax2), fDz(rhs.fDz), fSPhi(rhs.fSPhi), 00150 fDPhi(rhs.fDPhi), sinCPhi(rhs.sinCPhi), cosCPhi(rhs.cosCPhi), 00151 cosHDPhiOT(rhs.cosHDPhiOT), cosHDPhiIT(rhs.cosHDPhiIT), 00152 sinSPhi(rhs.sinSPhi), cosSPhi(rhs.cosSPhi), sinEPhi(rhs.sinEPhi), 00153 cosEPhi(rhs.cosEPhi), fPhiFullCone(rhs.fPhiFullCone) 00154 { 00155 }
G4bool G4Cons::CalculateExtent | ( | const EAxis | pAxis, | |
const G4VoxelLimits & | pVoxelLimit, | |||
const G4AffineTransform & | pTransform, | |||
G4double & | pmin, | |||
G4double & | pmax | |||
) | const [virtual] |
Implements G4VSolid.
Definition at line 261 of file G4Cons.cc.
References G4VSolid::ClipBetweenSections(), G4VSolid::ClipCrossSection(), G4VoxelLimits::GetMaxExtent(), G4VoxelLimits::GetMaxXExtent(), G4VoxelLimits::GetMaxYExtent(), G4VoxelLimits::GetMaxZExtent(), G4VoxelLimits::GetMinExtent(), G4VoxelLimits::GetMinXExtent(), G4VoxelLimits::GetMinYExtent(), G4VoxelLimits::GetMinZExtent(), Inside(), G4AffineTransform::Inverse(), G4AffineTransform::IsRotated(), G4VoxelLimits::IsXLimited(), G4VoxelLimits::IsYLimited(), G4VoxelLimits::IsZLimited(), G4VSolid::kCarTolerance, kOutside, kXAxis, kYAxis, kZAxis, G4AffineTransform::NetTranslation(), and G4AffineTransform::TransformPoint().
00266 { 00267 if ( !pTransform.IsRotated() && (fDPhi == twopi) 00268 && (fRmin1 == 0) && (fRmin2 == 0) ) 00269 { 00270 // Special case handling for unrotated solid cones 00271 // Compute z/x/y mins and maxs for bounding box respecting limits, 00272 // with early returns if outside limits. Then switch() on pAxis, 00273 // and compute exact x and y limit for x/y case 00274 00275 G4double xoffset, xMin, xMax ; 00276 G4double yoffset, yMin, yMax ; 00277 G4double zoffset, zMin, zMax ; 00278 00279 G4double diff1, diff2, delta, maxDiff, newMin, newMax, RMax ; 00280 G4double xoff1, xoff2, yoff1, yoff2 ; 00281 00282 zoffset = pTransform.NetTranslation().z(); 00283 zMin = zoffset - fDz ; 00284 zMax = zoffset + fDz ; 00285 00286 if (pVoxelLimit.IsZLimited()) 00287 { 00288 if( (zMin > pVoxelLimit.GetMaxZExtent() + kCarTolerance) || 00289 (zMax < pVoxelLimit.GetMinZExtent() - kCarTolerance) ) 00290 { 00291 return false ; 00292 } 00293 else 00294 { 00295 if ( zMin < pVoxelLimit.GetMinZExtent() ) 00296 { 00297 zMin = pVoxelLimit.GetMinZExtent() ; 00298 } 00299 if ( zMax > pVoxelLimit.GetMaxZExtent() ) 00300 { 00301 zMax = pVoxelLimit.GetMaxZExtent() ; 00302 } 00303 } 00304 } 00305 xoffset = pTransform.NetTranslation().x() ; 00306 RMax = (fRmax2 >= fRmax1) ? zMax : zMin ; 00307 xMax = xoffset + (fRmax1 + fRmax2)*0.5 + 00308 (RMax - zoffset)*(fRmax2 - fRmax1)/(2*fDz) ; 00309 xMin = 2*xoffset-xMax ; 00310 00311 if (pVoxelLimit.IsXLimited()) 00312 { 00313 if ( (xMin > pVoxelLimit.GetMaxXExtent() + kCarTolerance) || 00314 (xMax < pVoxelLimit.GetMinXExtent() - kCarTolerance) ) 00315 { 00316 return false ; 00317 } 00318 else 00319 { 00320 if ( xMin < pVoxelLimit.GetMinXExtent() ) 00321 { 00322 xMin = pVoxelLimit.GetMinXExtent() ; 00323 } 00324 if ( xMax > pVoxelLimit.GetMaxXExtent() ) 00325 { 00326 xMax=pVoxelLimit.GetMaxXExtent() ; 00327 } 00328 } 00329 } 00330 yoffset = pTransform.NetTranslation().y() ; 00331 yMax = yoffset + (fRmax1 + fRmax2)*0.5 + 00332 (RMax - zoffset)*(fRmax2 - fRmax1)/(2*fDz) ; 00333 yMin = 2*yoffset-yMax ; 00334 RMax = yMax - yoffset ; // = max radius due to Zmax/Zmin cuttings 00335 00336 if (pVoxelLimit.IsYLimited()) 00337 { 00338 if ( (yMin > pVoxelLimit.GetMaxYExtent() + kCarTolerance) || 00339 (yMax < pVoxelLimit.GetMinYExtent() - kCarTolerance) ) 00340 { 00341 return false ; 00342 } 00343 else 00344 { 00345 if ( yMin < pVoxelLimit.GetMinYExtent() ) 00346 { 00347 yMin = pVoxelLimit.GetMinYExtent() ; 00348 } 00349 if ( yMax > pVoxelLimit.GetMaxYExtent() ) 00350 { 00351 yMax = pVoxelLimit.GetMaxYExtent() ; 00352 } 00353 } 00354 } 00355 switch (pAxis) // Known to cut cones 00356 { 00357 case kXAxis: 00358 yoff1 = yoffset - yMin ; 00359 yoff2 = yMax - yoffset ; 00360 00361 if ((yoff1 >= 0) && (yoff2 >= 0)) // Y limits cross max/min x 00362 { // => no change 00363 pMin = xMin ; 00364 pMax = xMax ; 00365 } 00366 else 00367 { 00368 // Y limits don't cross max/min x => compute max delta x, 00369 // hence new mins/maxs 00370 delta=RMax*RMax-yoff1*yoff1; 00371 diff1=(delta>0.) ? std::sqrt(delta) : 0.; 00372 delta=RMax*RMax-yoff2*yoff2; 00373 diff2=(delta>0.) ? std::sqrt(delta) : 0.; 00374 maxDiff = (diff1>diff2) ? diff1:diff2 ; 00375 newMin = xoffset - maxDiff ; 00376 newMax = xoffset + maxDiff ; 00377 pMin = ( newMin < xMin ) ? xMin : newMin ; 00378 pMax = ( newMax > xMax) ? xMax : newMax ; 00379 } 00380 break ; 00381 00382 case kYAxis: 00383 xoff1 = xoffset - xMin ; 00384 xoff2 = xMax - xoffset ; 00385 00386 if ((xoff1 >= 0) && (xoff2 >= 0) ) // X limits cross max/min y 00387 { // => no change 00388 pMin = yMin ; 00389 pMax = yMax ; 00390 } 00391 else 00392 { 00393 // X limits don't cross max/min y => compute max delta y, 00394 // hence new mins/maxs 00395 delta=RMax*RMax-xoff1*xoff1; 00396 diff1=(delta>0.) ? std::sqrt(delta) : 0.; 00397 delta=RMax*RMax-xoff2*xoff2; 00398 diff2=(delta>0.) ? std::sqrt(delta) : 0.; 00399 maxDiff = (diff1 > diff2) ? diff1:diff2 ; 00400 newMin = yoffset - maxDiff ; 00401 newMax = yoffset + maxDiff ; 00402 pMin = (newMin < yMin) ? yMin : newMin ; 00403 pMax = (newMax > yMax) ? yMax : newMax ; 00404 } 00405 break ; 00406 00407 case kZAxis: 00408 pMin = zMin ; 00409 pMax = zMax ; 00410 break ; 00411 00412 default: 00413 break ; 00414 } 00415 pMin -= kCarTolerance ; 00416 pMax += kCarTolerance ; 00417 00418 return true ; 00419 } 00420 else // Calculate rotated vertex coordinates 00421 { 00422 G4int i, noEntries, noBetweenSections4 ; 00423 G4bool existsAfterClip = false ; 00424 G4ThreeVectorList* vertices = CreateRotatedVertices(pTransform) ; 00425 00426 pMin = +kInfinity ; 00427 pMax = -kInfinity ; 00428 00429 noEntries = vertices->size() ; 00430 noBetweenSections4 = noEntries-4 ; 00431 00432 for ( i = 0 ; i < noEntries ; i += 4 ) 00433 { 00434 ClipCrossSection(vertices, i, pVoxelLimit, pAxis, pMin, pMax) ; 00435 } 00436 for ( i = 0 ; i < noBetweenSections4 ; i += 4 ) 00437 { 00438 ClipBetweenSections(vertices, i, pVoxelLimit, pAxis, pMin, pMax) ; 00439 } 00440 if ( (pMin != kInfinity) || (pMax != -kInfinity) ) 00441 { 00442 existsAfterClip = true ; 00443 00444 // Add 2*tolerance to avoid precision troubles 00445 00446 pMin -= kCarTolerance ; 00447 pMax += kCarTolerance ; 00448 } 00449 else 00450 { 00451 // Check for case where completely enveloping clipping volume 00452 // If point inside then we are confident that the solid completely 00453 // envelopes the clipping volume. Hence set min/max extents according 00454 // to clipping volume extents along the specified axis. 00455 00456 G4ThreeVector clipCentre( 00457 (pVoxelLimit.GetMinXExtent() + pVoxelLimit.GetMaxXExtent())*0.5, 00458 (pVoxelLimit.GetMinYExtent() + pVoxelLimit.GetMaxYExtent())*0.5, 00459 (pVoxelLimit.GetMinZExtent() + pVoxelLimit.GetMaxZExtent())*0.5 ) ; 00460 00461 if (Inside(pTransform.Inverse().TransformPoint(clipCentre)) != kOutside) 00462 { 00463 existsAfterClip = true ; 00464 pMin = pVoxelLimit.GetMinExtent(pAxis) ; 00465 pMax = pVoxelLimit.GetMaxExtent(pAxis) ; 00466 } 00467 } 00468 delete vertices ; 00469 return existsAfterClip ; 00470 } 00471 }
G4VSolid * G4Cons::Clone | ( | ) | const [virtual] |
void G4Cons::ComputeDimensions | ( | G4VPVParameterisation * | p, | |
const G4int | n, | |||
const G4VPhysicalVolume * | pRep | |||
) | [virtual] |
Reimplemented from G4VSolid.
Definition at line 249 of file G4Cons.cc.
References G4VPVParameterisation::ComputeDimensions().
00252 { 00253 p->ComputeDimensions(*this,n,pRep) ; 00254 }
G4NURBS * G4Cons::CreateNURBS | ( | ) | const [virtual] |
Reimplemented from G4VSolid.
Definition at line 2387 of file G4Cons.cc.
02388 { 02389 G4double RMax = (fRmax2 >= fRmax1) ? fRmax2 : fRmax1 ; 02390 return new G4NURBSbox (RMax, RMax, fDz); // Box for now!!! 02391 }
G4Polyhedron * G4Cons::CreatePolyhedron | ( | ) | const [virtual] |
Reimplemented from G4VSolid.
Definition at line 2382 of file G4Cons.cc.
02383 { 02384 return new G4PolyhedronCons(fRmin1,fRmax1,fRmin2,fRmax2,fDz,fSPhi,fDPhi); 02385 }
void G4Cons::DescribeYourselfTo | ( | G4VGraphicsScene & | scene | ) | const [virtual] |
Implements G4VSolid.
Definition at line 2377 of file G4Cons.cc.
References G4VGraphicsScene::AddSolid().
02378 { 02379 scene.AddSolid (*this); 02380 }
G4double G4Cons::DistanceToIn | ( | const G4ThreeVector & | p | ) | const [virtual] |
Implements G4VSolid.
Definition at line 1384 of file G4Cons.cc.
01385 { 01386 G4double safe=0.0, rho, safeR1, safeR2, safeZ, safePhi, cosPsi ; 01387 G4double tanRMin, secRMin, pRMin ; 01388 G4double tanRMax, secRMax, pRMax ; 01389 01390 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ; 01391 safeZ = std::fabs(p.z()) - fDz ; 01392 01393 if ( fRmin1 || fRmin2 ) 01394 { 01395 tanRMin = (fRmin2 - fRmin1)*0.5/fDz ; 01396 secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ; 01397 pRMin = tanRMin*p.z() + (fRmin1 + fRmin2)*0.5 ; 01398 safeR1 = (pRMin - rho)/secRMin ; 01399 01400 tanRMax = (fRmax2 - fRmax1)*0.5/fDz ; 01401 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ; 01402 pRMax = tanRMax*p.z() + (fRmax1 + fRmax2)*0.5 ; 01403 safeR2 = (rho - pRMax)/secRMax ; 01404 01405 if ( safeR1 > safeR2) { safe = safeR1; } 01406 else { safe = safeR2; } 01407 } 01408 else 01409 { 01410 tanRMax = (fRmax2 - fRmax1)*0.5/fDz ; 01411 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ; 01412 pRMax = tanRMax*p.z() + (fRmax1 + fRmax2)*0.5 ; 01413 safe = (rho - pRMax)/secRMax ; 01414 } 01415 if ( safeZ > safe ) { safe = safeZ; } 01416 01417 if ( !fPhiFullCone && rho ) 01418 { 01419 // Psi=angle from central phi to point 01420 01421 cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/rho ; 01422 01423 if ( cosPsi < std::cos(fDPhi*0.5) ) // Point lies outside phi range 01424 { 01425 if ( (p.y()*cosCPhi - p.x()*sinCPhi) <= 0.0 ) 01426 { 01427 safePhi = std::fabs(p.x()*std::sin(fSPhi)-p.y()*std::cos(fSPhi)); 01428 } 01429 else 01430 { 01431 safePhi = std::fabs(p.x()*sinEPhi-p.y()*cosEPhi); 01432 } 01433 if ( safePhi > safe ) { safe = safePhi; } 01434 } 01435 } 01436 if ( safe < 0.0 ) { safe = 0.0; } 01437 01438 return safe ; 01439 }
G4double G4Cons::DistanceToIn | ( | const G4ThreeVector & | p, | |
const G4ThreeVector & | v | |||
) | const [virtual] |
Implements G4VSolid.
Definition at line 712 of file G4Cons.cc.
References G4VSolid::kCarTolerance.
00714 { 00715 G4double snxt = kInfinity ; // snxt = default return value 00716 const G4double dRmax = 50*(fRmax1+fRmax2);// 100*(Rmax1+Rmax2)/2. 00717 static const G4double halfCarTolerance=kCarTolerance*0.5; 00718 static const G4double halfRadTolerance=kRadTolerance*0.5; 00719 00720 G4double tanRMax,secRMax,rMaxAv,rMaxOAv ; // Data for cones 00721 G4double tanRMin,secRMin,rMinAv,rMinOAv ; 00722 G4double rout,rin ; 00723 00724 G4double tolORMin,tolORMin2,tolIRMin,tolIRMin2 ; // `generous' radii squared 00725 G4double tolORMax2,tolIRMax,tolIRMax2 ; 00726 G4double tolODz,tolIDz ; 00727 00728 G4double Dist,sd,xi,yi,zi,ri=0.,risec,rhoi2,cosPsi ; // Intersection point vars 00729 00730 G4double t1,t2,t3,b,c,d ; // Quadratic solver variables 00731 G4double nt1,nt2,nt3 ; 00732 G4double Comp ; 00733 00734 G4ThreeVector Normal; 00735 00736 // Cone Precalcs 00737 00738 tanRMin = (fRmin2 - fRmin1)*0.5/fDz ; 00739 secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ; 00740 rMinAv = (fRmin1 + fRmin2)*0.5 ; 00741 00742 if (rMinAv > halfRadTolerance) 00743 { 00744 rMinOAv = rMinAv - halfRadTolerance ; 00745 } 00746 else 00747 { 00748 rMinOAv = 0.0 ; 00749 } 00750 tanRMax = (fRmax2 - fRmax1)*0.5/fDz ; 00751 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ; 00752 rMaxAv = (fRmax1 + fRmax2)*0.5 ; 00753 rMaxOAv = rMaxAv + halfRadTolerance ; 00754 00755 // Intersection with z-surfaces 00756 00757 tolIDz = fDz - halfCarTolerance ; 00758 tolODz = fDz + halfCarTolerance ; 00759 00760 if (std::fabs(p.z()) >= tolIDz) 00761 { 00762 if ( p.z()*v.z() < 0 ) // at +Z going in -Z or visa versa 00763 { 00764 sd = (std::fabs(p.z()) - fDz)/std::fabs(v.z()) ; // Z intersect distance 00765 00766 if( sd < 0.0 ) { sd = 0.0; } // negative dist -> zero 00767 00768 xi = p.x() + sd*v.x() ; // Intersection coords 00769 yi = p.y() + sd*v.y() ; 00770 rhoi2 = xi*xi + yi*yi ; 00771 00772 // Check validity of intersection 00773 // Calculate (outer) tolerant radi^2 at intersecion 00774 00775 if (v.z() > 0) 00776 { 00777 tolORMin = fRmin1 - halfRadTolerance*secRMin ; 00778 tolIRMin = fRmin1 + halfRadTolerance*secRMin ; 00779 tolIRMax = fRmax1 - halfRadTolerance*secRMin ; 00780 tolORMax2 = (fRmax1 + halfRadTolerance*secRMax)* 00781 (fRmax1 + halfRadTolerance*secRMax) ; 00782 } 00783 else 00784 { 00785 tolORMin = fRmin2 - halfRadTolerance*secRMin ; 00786 tolIRMin = fRmin2 + halfRadTolerance*secRMin ; 00787 tolIRMax = fRmax2 - halfRadTolerance*secRMin ; 00788 tolORMax2 = (fRmax2 + halfRadTolerance*secRMax)* 00789 (fRmax2 + halfRadTolerance*secRMax) ; 00790 } 00791 if ( tolORMin > 0 ) 00792 { 00793 tolORMin2 = tolORMin*tolORMin ; 00794 tolIRMin2 = tolIRMin*tolIRMin ; 00795 } 00796 else 00797 { 00798 tolORMin2 = 0.0 ; 00799 tolIRMin2 = 0.0 ; 00800 } 00801 if ( tolIRMax > 0 ) { tolIRMax2 = tolIRMax*tolIRMax; } 00802 else { tolIRMax2 = 0.0; } 00803 00804 if ( (tolIRMin2 <= rhoi2) && (rhoi2 <= tolIRMax2) ) 00805 { 00806 if ( !fPhiFullCone && rhoi2 ) 00807 { 00808 // Psi = angle made with central (average) phi of shape 00809 00810 cosPsi = (xi*cosCPhi + yi*sinCPhi)/std::sqrt(rhoi2) ; 00811 00812 if (cosPsi >= cosHDPhiIT) { return sd; } 00813 } 00814 else 00815 { 00816 return sd; 00817 } 00818 } 00819 } 00820 else // On/outside extent, and heading away -> cannot intersect 00821 { 00822 return snxt ; 00823 } 00824 } 00825 00826 // ----> Can not intersect z surfaces 00827 00828 00829 // Intersection with outer cone (possible return) and 00830 // inner cone (must also check phi) 00831 // 00832 // Intersection point (xi,yi,zi) on line x=p.x+t*v.x etc. 00833 // 00834 // Intersects with x^2+y^2=(a*z+b)^2 00835 // 00836 // where a=tanRMax or tanRMin 00837 // b=rMaxAv or rMinAv 00838 // 00839 // (vx^2+vy^2-(a*vz)^2)t^2+2t(pxvx+pyvy-a*vz(a*pz+b))+px^2+py^2-(a*pz+b)^2=0 ; 00840 // t1 t2 t3 00841 // 00842 // \--------u-------/ \-----------v----------/ \---------w--------/ 00843 // 00844 00845 t1 = 1.0 - v.z()*v.z() ; 00846 t2 = p.x()*v.x() + p.y()*v.y() ; 00847 t3 = p.x()*p.x() + p.y()*p.y() ; 00848 rin = tanRMin*p.z() + rMinAv ; 00849 rout = tanRMax*p.z() + rMaxAv ; 00850 00851 // Outer Cone Intersection 00852 // Must be outside/on outer cone for valid intersection 00853 00854 nt1 = t1 - (tanRMax*v.z())*(tanRMax*v.z()) ; 00855 nt2 = t2 - tanRMax*v.z()*rout ; 00856 nt3 = t3 - rout*rout ; 00857 00858 if (std::fabs(nt1) > kRadTolerance) // Equation quadratic => 2 roots 00859 { 00860 b = nt2/nt1; 00861 c = nt3/nt1; 00862 d = b*b-c ; 00863 if ( (nt3 > rout*rout*kRadTolerance*kRadTolerance*secRMax*secRMax) 00864 || (rout < 0) ) 00865 { 00866 // If outside real cone (should be rho-rout>kRadTolerance*0.5 00867 // NOT rho^2 etc) saves a std::sqrt() at expense of accuracy 00868 00869 if (d >= 0) 00870 { 00871 00872 if ((rout < 0) && (nt3 <= 0)) 00873 { 00874 // Inside `shadow cone' with -ve radius 00875 // -> 2nd root could be on real cone 00876 00877 if (b>0) { sd = c/(-b-std::sqrt(d)); } 00878 else { sd = -b + std::sqrt(d); } 00879 } 00880 else 00881 { 00882 if ((b <= 0) && (c >= 0)) // both >=0, try smaller root 00883 { 00884 sd=c/(-b+std::sqrt(d)); 00885 } 00886 else 00887 { 00888 if ( c <= 0 ) // second >=0 00889 { 00890 sd = -b + std::sqrt(d) ; 00891 if((sd<0) & (sd>-halfRadTolerance)) sd=0; 00892 } 00893 else // both negative, travel away 00894 { 00895 return kInfinity ; 00896 } 00897 } 00898 } 00899 if ( sd > 0 ) // If 'forwards'. Check z intersection 00900 { 00901 if ( sd>dRmax ) // Avoid rounding errors due to precision issues on 00902 { // 64 bits systems. Split long distances and recompute 00903 G4double fTerm = sd-std::fmod(sd,dRmax); 00904 sd = fTerm + DistanceToIn(p+fTerm*v,v); 00905 } 00906 zi = p.z() + sd*v.z() ; 00907 00908 if (std::fabs(zi) <= tolODz) 00909 { 00910 // Z ok. Check phi intersection if reqd 00911 00912 if ( fPhiFullCone ) { return sd; } 00913 else 00914 { 00915 xi = p.x() + sd*v.x() ; 00916 yi = p.y() + sd*v.y() ; 00917 ri = rMaxAv + zi*tanRMax ; 00918 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ; 00919 00920 if ( cosPsi >= cosHDPhiIT ) { return sd; } 00921 } 00922 } 00923 } // end if (sd>0) 00924 } 00925 } 00926 else 00927 { 00928 // Inside outer cone 00929 // check not inside, and heading through G4Cons (-> 0 to in) 00930 00931 if ( ( t3 > (rin + halfRadTolerance*secRMin)* 00932 (rin + halfRadTolerance*secRMin) ) 00933 && (nt2 < 0) && (d >= 0) && (std::fabs(p.z()) <= tolIDz) ) 00934 { 00935 // Inside cones, delta r -ve, inside z extent 00936 // Point is on the Surface => check Direction using Normal.dot(v) 00937 00938 xi = p.x() ; 00939 yi = p.y() ; 00940 risec = std::sqrt(xi*xi + yi*yi)*secRMax ; 00941 Normal = G4ThreeVector(xi/risec,yi/risec,-tanRMax/secRMax) ; 00942 if ( !fPhiFullCone ) 00943 { 00944 cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(t3) ; 00945 if ( cosPsi >= cosHDPhiIT ) 00946 { 00947 if ( Normal.dot(v) <= 0 ) { return 0.0; } 00948 } 00949 } 00950 else 00951 { 00952 if ( Normal.dot(v) <= 0 ) { return 0.0; } 00953 } 00954 } 00955 } 00956 } 00957 else // Single root case 00958 { 00959 if ( std::fabs(nt2) > kRadTolerance ) 00960 { 00961 sd = -0.5*nt3/nt2 ; 00962 00963 if ( sd < 0 ) { return kInfinity; } // travel away 00964 else // sd >= 0, If 'forwards'. Check z intersection 00965 { 00966 zi = p.z() + sd*v.z() ; 00967 00968 if ((std::fabs(zi) <= tolODz) && (nt2 < 0)) 00969 { 00970 // Z ok. Check phi intersection if reqd 00971 00972 if ( fPhiFullCone ) { return sd; } 00973 else 00974 { 00975 xi = p.x() + sd*v.x() ; 00976 yi = p.y() + sd*v.y() ; 00977 ri = rMaxAv + zi*tanRMax ; 00978 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ; 00979 00980 if (cosPsi >= cosHDPhiIT) { return sd; } 00981 } 00982 } 00983 } 00984 } 00985 else // travel || cone surface from its origin 00986 { 00987 sd = kInfinity ; 00988 } 00989 } 00990 00991 // Inner Cone Intersection 00992 // o Space is divided into 3 areas: 00993 // 1) Radius greater than real inner cone & imaginary cone & outside 00994 // tolerance 00995 // 2) Radius less than inner or imaginary cone & outside tolarance 00996 // 3) Within tolerance of real or imaginary cones 00997 // - Extra checks needed for 3's intersections 00998 // => lots of duplicated code 00999 01000 if (rMinAv) 01001 { 01002 nt1 = t1 - (tanRMin*v.z())*(tanRMin*v.z()) ; 01003 nt2 = t2 - tanRMin*v.z()*rin ; 01004 nt3 = t3 - rin*rin ; 01005 01006 if ( nt1 ) 01007 { 01008 if ( nt3 > rin*kRadTolerance*secRMin ) 01009 { 01010 // At radius greater than real & imaginary cones 01011 // -> 2nd root, with zi check 01012 01013 b = nt2/nt1 ; 01014 c = nt3/nt1 ; 01015 d = b*b-c ; 01016 if (d >= 0) // > 0 01017 { 01018 if(b>0){sd = c/( -b-std::sqrt(d));} 01019 else {sd = -b + std::sqrt(d) ;} 01020 01021 if ( sd >= 0 ) // > 0 01022 { 01023 if ( sd>dRmax ) // Avoid rounding errors due to precision issues on 01024 { // 64 bits systems. Split long distance and recompute 01025 G4double fTerm = sd-std::fmod(sd,dRmax); 01026 sd = fTerm + DistanceToIn(p+fTerm*v,v); 01027 } 01028 zi = p.z() + sd*v.z() ; 01029 01030 if ( std::fabs(zi) <= tolODz ) 01031 { 01032 if ( !fPhiFullCone ) 01033 { 01034 xi = p.x() + sd*v.x() ; 01035 yi = p.y() + sd*v.y() ; 01036 ri = rMinAv + zi*tanRMin ; 01037 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ; 01038 01039 if (cosPsi >= cosHDPhiIT) 01040 { 01041 if ( sd > halfRadTolerance ) { snxt=sd; } 01042 else 01043 { 01044 // Calculate a normal vector in order to check Direction 01045 01046 risec = std::sqrt(xi*xi + yi*yi)*secRMin ; 01047 Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin); 01048 if ( Normal.dot(v) <= 0 ) { snxt = sd; } 01049 } 01050 } 01051 } 01052 else 01053 { 01054 if ( sd > halfRadTolerance ) { return sd; } 01055 else 01056 { 01057 // Calculate a normal vector in order to check Direction 01058 01059 xi = p.x() + sd*v.x() ; 01060 yi = p.y() + sd*v.y() ; 01061 risec = std::sqrt(xi*xi + yi*yi)*secRMin ; 01062 Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ; 01063 if ( Normal.dot(v) <= 0 ) { return sd; } 01064 } 01065 } 01066 } 01067 } 01068 } 01069 } 01070 else if ( nt3 < -rin*kRadTolerance*secRMin ) 01071 { 01072 // Within radius of inner cone (real or imaginary) 01073 // -> Try 2nd root, with checking intersection is with real cone 01074 // -> If check fails, try 1st root, also checking intersection is 01075 // on real cone 01076 01077 b = nt2/nt1 ; 01078 c = nt3/nt1 ; 01079 d = b*b - c ; 01080 01081 if ( d >= 0 ) // > 0 01082 { 01083 if (b>0) { sd = c/(-b-std::sqrt(d)); } 01084 else { sd = -b + std::sqrt(d); } 01085 zi = p.z() + sd*v.z() ; 01086 ri = rMinAv + zi*tanRMin ; 01087 01088 if ( ri > 0 ) 01089 { 01090 if ( (sd >= 0) && (std::fabs(zi) <= tolODz) ) // sd > 0 01091 { 01092 if ( sd>dRmax ) // Avoid rounding errors due to precision issues 01093 { // seen on 64 bits systems. Split and recompute 01094 G4double fTerm = sd-std::fmod(sd,dRmax); 01095 sd = fTerm + DistanceToIn(p+fTerm*v,v); 01096 } 01097 if ( !fPhiFullCone ) 01098 { 01099 xi = p.x() + sd*v.x() ; 01100 yi = p.y() + sd*v.y() ; 01101 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ; 01102 01103 if (cosPsi >= cosHDPhiOT) 01104 { 01105 if ( sd > halfRadTolerance ) { snxt=sd; } 01106 else 01107 { 01108 // Calculate a normal vector in order to check Direction 01109 01110 risec = std::sqrt(xi*xi + yi*yi)*secRMin ; 01111 Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin); 01112 if ( Normal.dot(v) <= 0 ) { snxt = sd; } 01113 } 01114 } 01115 } 01116 else 01117 { 01118 if( sd > halfRadTolerance ) { return sd; } 01119 else 01120 { 01121 // Calculate a normal vector in order to check Direction 01122 01123 xi = p.x() + sd*v.x() ; 01124 yi = p.y() + sd*v.y() ; 01125 risec = std::sqrt(xi*xi + yi*yi)*secRMin ; 01126 Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ; 01127 if ( Normal.dot(v) <= 0 ) { return sd; } 01128 } 01129 } 01130 } 01131 } 01132 else 01133 { 01134 if (b>0) { sd = -b - std::sqrt(d); } 01135 else { sd = c/(-b+std::sqrt(d)); } 01136 zi = p.z() + sd*v.z() ; 01137 ri = rMinAv + zi*tanRMin ; 01138 01139 if ( (sd >= 0) && (ri > 0) && (std::fabs(zi) <= tolODz) ) // sd>0 01140 { 01141 if ( sd>dRmax ) // Avoid rounding errors due to precision issues 01142 { // seen on 64 bits systems. Split and recompute 01143 G4double fTerm = sd-std::fmod(sd,dRmax); 01144 sd = fTerm + DistanceToIn(p+fTerm*v,v); 01145 } 01146 if ( !fPhiFullCone ) 01147 { 01148 xi = p.x() + sd*v.x() ; 01149 yi = p.y() + sd*v.y() ; 01150 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ; 01151 01152 if (cosPsi >= cosHDPhiIT) 01153 { 01154 if ( sd > halfRadTolerance ) { snxt=sd; } 01155 else 01156 { 01157 // Calculate a normal vector in order to check Direction 01158 01159 risec = std::sqrt(xi*xi + yi*yi)*secRMin ; 01160 Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin); 01161 if ( Normal.dot(v) <= 0 ) { snxt = sd; } 01162 } 01163 } 01164 } 01165 else 01166 { 01167 if ( sd > halfRadTolerance ) { return sd; } 01168 else 01169 { 01170 // Calculate a normal vector in order to check Direction 01171 01172 xi = p.x() + sd*v.x() ; 01173 yi = p.y() + sd*v.y() ; 01174 risec = std::sqrt(xi*xi + yi*yi)*secRMin ; 01175 Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ; 01176 if ( Normal.dot(v) <= 0 ) { return sd; } 01177 } 01178 } 01179 } 01180 } 01181 } 01182 } 01183 else 01184 { 01185 // Within kRadTol*0.5 of inner cone (real OR imaginary) 01186 // ----> Check not travelling through (=>0 to in) 01187 // ----> if not: 01188 // -2nd root with validity check 01189 01190 if ( std::fabs(p.z()) <= tolODz ) 01191 { 01192 if ( nt2 > 0 ) 01193 { 01194 // Inside inner real cone, heading outwards, inside z range 01195 01196 if ( !fPhiFullCone ) 01197 { 01198 cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/std::sqrt(t3) ; 01199 01200 if (cosPsi >= cosHDPhiIT) { return 0.0; } 01201 } 01202 else { return 0.0; } 01203 } 01204 else 01205 { 01206 // Within z extent, but not travelling through 01207 // -> 2nd root or kInfinity if 1st root on imaginary cone 01208 01209 b = nt2/nt1 ; 01210 c = nt3/nt1 ; 01211 d = b*b - c ; 01212 01213 if ( d >= 0 ) // > 0 01214 { 01215 if (b>0) { sd = -b - std::sqrt(d); } 01216 else { sd = c/(-b+std::sqrt(d)); } 01217 zi = p.z() + sd*v.z() ; 01218 ri = rMinAv + zi*tanRMin ; 01219 01220 if ( ri > 0 ) // 2nd root 01221 { 01222 if (b>0) { sd = c/(-b-std::sqrt(d)); } 01223 else { sd = -b + std::sqrt(d); } 01224 01225 zi = p.z() + sd*v.z() ; 01226 01227 if ( (sd >= 0) && (std::fabs(zi) <= tolODz) ) // sd>0 01228 { 01229 if ( sd>dRmax ) // Avoid rounding errors due to precision issue 01230 { // seen on 64 bits systems. Split and recompute 01231 G4double fTerm = sd-std::fmod(sd,dRmax); 01232 sd = fTerm + DistanceToIn(p+fTerm*v,v); 01233 } 01234 if ( !fPhiFullCone ) 01235 { 01236 xi = p.x() + sd*v.x() ; 01237 yi = p.y() + sd*v.y() ; 01238 ri = rMinAv + zi*tanRMin ; 01239 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri ; 01240 01241 if ( cosPsi >= cosHDPhiIT ) { snxt = sd; } 01242 } 01243 else { return sd; } 01244 } 01245 } 01246 else { return kInfinity; } 01247 } 01248 } 01249 } 01250 else // 2nd root 01251 { 01252 b = nt2/nt1 ; 01253 c = nt3/nt1 ; 01254 d = b*b - c ; 01255 01256 if ( d > 0 ) 01257 { 01258 if (b>0) { sd = c/(-b-std::sqrt(d)); } 01259 else { sd = -b + std::sqrt(d) ; } 01260 zi = p.z() + sd*v.z() ; 01261 01262 if ( (sd >= 0) && (std::fabs(zi) <= tolODz) ) // sd>0 01263 { 01264 if ( sd>dRmax ) // Avoid rounding errors due to precision issues 01265 { // seen on 64 bits systems. Split and recompute 01266 G4double fTerm = sd-std::fmod(sd,dRmax); 01267 sd = fTerm + DistanceToIn(p+fTerm*v,v); 01268 } 01269 if ( !fPhiFullCone ) 01270 { 01271 xi = p.x() + sd*v.x(); 01272 yi = p.y() + sd*v.y(); 01273 ri = rMinAv + zi*tanRMin ; 01274 cosPsi = (xi*cosCPhi + yi*sinCPhi)/ri; 01275 01276 if (cosPsi >= cosHDPhiIT) { snxt = sd; } 01277 } 01278 else { return sd; } 01279 } 01280 } 01281 } 01282 } 01283 } 01284 } 01285 01286 // Phi segment intersection 01287 // 01288 // o Tolerant of points inside phi planes by up to kCarTolerance*0.5 01289 // 01290 // o NOTE: Large duplication of code between sphi & ephi checks 01291 // -> only diffs: sphi -> ephi, Comp -> -Comp and half-plane 01292 // intersection check <=0 -> >=0 01293 // -> Should use some form of loop Construct 01294 01295 if ( !fPhiFullCone ) 01296 { 01297 // First phi surface (starting phi) 01298 01299 Comp = v.x()*sinSPhi - v.y()*cosSPhi ; 01300 01301 if ( Comp < 0 ) // Component in outwards normal dirn 01302 { 01303 Dist = (p.y()*cosSPhi - p.x()*sinSPhi) ; 01304 01305 if (Dist < halfCarTolerance) 01306 { 01307 sd = Dist/Comp ; 01308 01309 if ( sd < snxt ) 01310 { 01311 if ( sd < 0 ) { sd = 0.0; } 01312 01313 zi = p.z() + sd*v.z() ; 01314 01315 if ( std::fabs(zi) <= tolODz ) 01316 { 01317 xi = p.x() + sd*v.x() ; 01318 yi = p.y() + sd*v.y() ; 01319 rhoi2 = xi*xi + yi*yi ; 01320 tolORMin2 = (rMinOAv + zi*tanRMin)*(rMinOAv + zi*tanRMin) ; 01321 tolORMax2 = (rMaxOAv + zi*tanRMax)*(rMaxOAv + zi*tanRMax) ; 01322 01323 if ( (rhoi2 >= tolORMin2) && (rhoi2 <= tolORMax2) ) 01324 { 01325 // z and r intersections good - check intersecting with 01326 // correct half-plane 01327 01328 if ((yi*cosCPhi - xi*sinCPhi) <= 0 ) { snxt = sd; } 01329 } 01330 } 01331 } 01332 } 01333 } 01334 01335 // Second phi surface (Ending phi) 01336 01337 Comp = -(v.x()*sinEPhi - v.y()*cosEPhi) ; 01338 01339 if ( Comp < 0 ) // Component in outwards normal dirn 01340 { 01341 Dist = -(p.y()*cosEPhi - p.x()*sinEPhi) ; 01342 if (Dist < halfCarTolerance) 01343 { 01344 sd = Dist/Comp ; 01345 01346 if ( sd < snxt ) 01347 { 01348 if ( sd < 0 ) { sd = 0.0; } 01349 01350 zi = p.z() + sd*v.z() ; 01351 01352 if (std::fabs(zi) <= tolODz) 01353 { 01354 xi = p.x() + sd*v.x() ; 01355 yi = p.y() + sd*v.y() ; 01356 rhoi2 = xi*xi + yi*yi ; 01357 tolORMin2 = (rMinOAv + zi*tanRMin)*(rMinOAv + zi*tanRMin) ; 01358 tolORMax2 = (rMaxOAv + zi*tanRMax)*(rMaxOAv + zi*tanRMax) ; 01359 01360 if ( (rhoi2 >= tolORMin2) && (rhoi2 <= tolORMax2) ) 01361 { 01362 // z and r intersections good - check intersecting with 01363 // correct half-plane 01364 01365 if ( (yi*cosCPhi - xi*sinCPhi) >= 0.0 ) { snxt = sd; } 01366 } 01367 } 01368 } 01369 } 01370 } 01371 } 01372 if (snxt < halfCarTolerance) { snxt = 0.; } 01373 01374 return snxt ; 01375 }
G4double G4Cons::DistanceToOut | ( | const G4ThreeVector & | p | ) | const [virtual] |
Implements G4VSolid.
Definition at line 2073 of file G4Cons.cc.
References G4VSolid::DumpInfo(), G4cout, G4endl, G4Exception(), Inside(), JustWarning, and kOutside.
02074 { 02075 G4double safe=0.0, rho, safeR1, safeR2, safeZ, safePhi; 02076 G4double tanRMin, secRMin, pRMin; 02077 G4double tanRMax, secRMax, pRMax; 02078 02079 #ifdef G4CSGDEBUG 02080 if( Inside(p) == kOutside ) 02081 { 02082 G4int oldprc=G4cout.precision(16) ; 02083 G4cout << G4endl ; 02084 DumpInfo(); 02085 G4cout << "Position:" << G4endl << G4endl ; 02086 G4cout << "p.x() = " << p.x()/mm << " mm" << G4endl ; 02087 G4cout << "p.y() = " << p.y()/mm << " mm" << G4endl ; 02088 G4cout << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl ; 02089 G4cout << "pho at z = " << std::sqrt( p.x()*p.x()+p.y()*p.y() )/mm 02090 << " mm" << G4endl << G4endl ; 02091 if( (p.x() != 0.) || (p.x() != 0.) ) 02092 { 02093 G4cout << "point phi = " << std::atan2(p.y(),p.x())/degree 02094 << " degree" << G4endl << G4endl ; 02095 } 02096 G4cout.precision(oldprc) ; 02097 G4Exception("G4Cons::DistanceToOut(p)", "GeomSolids1002", 02098 JustWarning, "Point p is outside !?" ); 02099 } 02100 #endif 02101 02102 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ; 02103 safeZ = fDz - std::fabs(p.z()) ; 02104 02105 if (fRmin1 || fRmin2) 02106 { 02107 tanRMin = (fRmin2 - fRmin1)*0.5/fDz ; 02108 secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ; 02109 pRMin = tanRMin*p.z() + (fRmin1 + fRmin2)*0.5 ; 02110 safeR1 = (rho - pRMin)/secRMin ; 02111 } 02112 else 02113 { 02114 safeR1 = kInfinity ; 02115 } 02116 02117 tanRMax = (fRmax2 - fRmax1)*0.5/fDz ; 02118 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ; 02119 pRMax = tanRMax*p.z() + (fRmax1+fRmax2)*0.5 ; 02120 safeR2 = (pRMax - rho)/secRMax ; 02121 02122 if (safeR1 < safeR2) { safe = safeR1; } 02123 else { safe = safeR2; } 02124 if (safeZ < safe) { safe = safeZ ; } 02125 02126 // Check if phi divided, Calc distances closest phi plane 02127 02128 if (!fPhiFullCone) 02129 { 02130 // Above/below central phi of G4Cons? 02131 02132 if ( (p.y()*cosCPhi - p.x()*sinCPhi) <= 0 ) 02133 { 02134 safePhi = -(p.x()*sinSPhi - p.y()*cosSPhi) ; 02135 } 02136 else 02137 { 02138 safePhi = (p.x()*sinEPhi - p.y()*cosEPhi) ; 02139 } 02140 if (safePhi < safe) { safe = safePhi; } 02141 } 02142 if ( safe < 0 ) { safe = 0; } 02143 02144 return safe ; 02145 }
G4double G4Cons::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 1446 of file G4Cons.cc.
References G4VSolid::DumpInfo(), G4cout, G4endl, G4Exception(), JustWarning, G4VSolid::kCarTolerance, and G4INCL::Math::pi.
01451 { 01452 ESide side = kNull, sider = kNull, sidephi = kNull; 01453 01454 static const G4double halfCarTolerance=kCarTolerance*0.5; 01455 static const G4double halfRadTolerance=kRadTolerance*0.5; 01456 static const G4double halfAngTolerance=kAngTolerance*0.5; 01457 01458 G4double snxt,srd,sphi,pdist ; 01459 01460 G4double tanRMax, secRMax, rMaxAv ; // Data for outer cone 01461 G4double tanRMin, secRMin, rMinAv ; // Data for inner cone 01462 01463 G4double t1, t2, t3, rout, rin, nt1, nt2, nt3 ; 01464 G4double b, c, d, sr2, sr3 ; 01465 01466 // Vars for intersection within tolerance 01467 01468 ESide sidetol = kNull ; 01469 G4double slentol = kInfinity ; 01470 01471 // Vars for phi intersection: 01472 01473 G4double pDistS, compS, pDistE, compE, sphi2, xi, yi, risec, vphi ; 01474 G4double zi, ri, deltaRoi2 ; 01475 01476 // Z plane intersection 01477 01478 if ( v.z() > 0.0 ) 01479 { 01480 pdist = fDz - p.z() ; 01481 01482 if (pdist > halfCarTolerance) 01483 { 01484 snxt = pdist/v.z() ; 01485 side = kPZ ; 01486 } 01487 else 01488 { 01489 if (calcNorm) 01490 { 01491 *n = G4ThreeVector(0,0,1) ; 01492 *validNorm = true ; 01493 } 01494 return snxt = 0.0; 01495 } 01496 } 01497 else if ( v.z() < 0.0 ) 01498 { 01499 pdist = fDz + p.z() ; 01500 01501 if ( pdist > halfCarTolerance) 01502 { 01503 snxt = -pdist/v.z() ; 01504 side = kMZ ; 01505 } 01506 else 01507 { 01508 if ( calcNorm ) 01509 { 01510 *n = G4ThreeVector(0,0,-1) ; 01511 *validNorm = true ; 01512 } 01513 return snxt = 0.0 ; 01514 } 01515 } 01516 else // Travel perpendicular to z axis 01517 { 01518 snxt = kInfinity ; 01519 side = kNull ; 01520 } 01521 01522 // Radial Intersections 01523 // 01524 // Intersection with outer cone (possible return) and 01525 // inner cone (must also check phi) 01526 // 01527 // Intersection point (xi,yi,zi) on line x=p.x+t*v.x etc. 01528 // 01529 // Intersects with x^2+y^2=(a*z+b)^2 01530 // 01531 // where a=tanRMax or tanRMin 01532 // b=rMaxAv or rMinAv 01533 // 01534 // (vx^2+vy^2-(a*vz)^2)t^2+2t(pxvx+pyvy-a*vz(a*pz+b))+px^2+py^2-(a*pz+b)^2=0 ; 01535 // t1 t2 t3 01536 // 01537 // \--------u-------/ \-----------v----------/ \---------w--------/ 01538 01539 tanRMax = (fRmax2 - fRmax1)*0.5/fDz ; 01540 secRMax = std::sqrt(1.0 + tanRMax*tanRMax) ; 01541 rMaxAv = (fRmax1 + fRmax2)*0.5 ; 01542 01543 01544 t1 = 1.0 - v.z()*v.z() ; // since v normalised 01545 t2 = p.x()*v.x() + p.y()*v.y() ; 01546 t3 = p.x()*p.x() + p.y()*p.y() ; 01547 rout = tanRMax*p.z() + rMaxAv ; 01548 01549 nt1 = t1 - (tanRMax*v.z())*(tanRMax*v.z()) ; 01550 nt2 = t2 - tanRMax*v.z()*rout ; 01551 nt3 = t3 - rout*rout ; 01552 01553 if (v.z() > 0.0) 01554 { 01555 deltaRoi2 = snxt*snxt*t1 + 2*snxt*t2 + t3 01556 - fRmax2*(fRmax2 + kRadTolerance*secRMax); 01557 } 01558 else if ( v.z() < 0.0 ) 01559 { 01560 deltaRoi2 = snxt*snxt*t1 + 2*snxt*t2 + t3 01561 - fRmax1*(fRmax1 + kRadTolerance*secRMax); 01562 } 01563 else 01564 { 01565 deltaRoi2 = 1.0; 01566 } 01567 01568 if ( nt1 && (deltaRoi2 > 0.0) ) 01569 { 01570 // Equation quadratic => 2 roots : second root must be leaving 01571 01572 b = nt2/nt1 ; 01573 c = nt3/nt1 ; 01574 d = b*b - c ; 01575 01576 if ( d >= 0 ) 01577 { 01578 // Check if on outer cone & heading outwards 01579 // NOTE: Should use rho-rout>-kRadTolerance*0.5 01580 01581 if (nt3 > -halfRadTolerance && nt2 >= 0 ) 01582 { 01583 if (calcNorm) 01584 { 01585 risec = std::sqrt(t3)*secRMax ; 01586 *validNorm = true ; 01587 *n = G4ThreeVector(p.x()/risec,p.y()/risec,-tanRMax/secRMax); 01588 } 01589 return snxt=0 ; 01590 } 01591 else 01592 { 01593 sider = kRMax ; 01594 if (b>0) { srd = -b - std::sqrt(d); } 01595 else { srd = c/(-b+std::sqrt(d)) ; } 01596 01597 zi = p.z() + srd*v.z() ; 01598 ri = tanRMax*zi + rMaxAv ; 01599 01600 if ((ri >= 0) && (-halfRadTolerance <= srd) && (srd <= halfRadTolerance)) 01601 { 01602 // An intersection within the tolerance 01603 // we will Store it in case it is good - 01604 // 01605 slentol = srd ; 01606 sidetol = kRMax ; 01607 } 01608 if ( (ri < 0) || (srd < halfRadTolerance) ) 01609 { 01610 // Safety: if both roots -ve ensure that srd cannot `win' 01611 // distance to out 01612 01613 if (b>0) { sr2 = c/(-b-std::sqrt(d)); } 01614 else { sr2 = -b + std::sqrt(d); } 01615 zi = p.z() + sr2*v.z() ; 01616 ri = tanRMax*zi + rMaxAv ; 01617 01618 if ((ri >= 0) && (sr2 > halfRadTolerance)) 01619 { 01620 srd = sr2; 01621 } 01622 else 01623 { 01624 srd = kInfinity ; 01625 01626 if( (-halfRadTolerance <= sr2) && ( sr2 <= halfRadTolerance) ) 01627 { 01628 // An intersection within the tolerance. 01629 // Storing it in case it is good. 01630 01631 slentol = sr2 ; 01632 sidetol = kRMax ; 01633 } 01634 } 01635 } 01636 } 01637 } 01638 else 01639 { 01640 // No intersection with outer cone & not parallel 01641 // -> already outside, no intersection 01642 01643 if ( calcNorm ) 01644 { 01645 risec = std::sqrt(t3)*secRMax; 01646 *validNorm = true; 01647 *n = G4ThreeVector(p.x()/risec,p.y()/risec,-tanRMax/secRMax); 01648 } 01649 return snxt = 0.0 ; 01650 } 01651 } 01652 else if ( nt2 && (deltaRoi2 > 0.0) ) 01653 { 01654 // Linear case (only one intersection) => point outside outer cone 01655 01656 if ( calcNorm ) 01657 { 01658 risec = std::sqrt(t3)*secRMax; 01659 *validNorm = true; 01660 *n = G4ThreeVector(p.x()/risec,p.y()/risec,-tanRMax/secRMax); 01661 } 01662 return snxt = 0.0 ; 01663 } 01664 else 01665 { 01666 // No intersection -> parallel to outer cone 01667 // => Z or inner cone intersection 01668 01669 srd = kInfinity ; 01670 } 01671 01672 // Check possible intersection within tolerance 01673 01674 if ( slentol <= halfCarTolerance ) 01675 { 01676 // An intersection within the tolerance was found. 01677 // We must accept it only if the momentum points outwards. 01678 // 01679 // G4ThreeVector ptTol ; // The point of the intersection 01680 // ptTol= p + slentol*v ; 01681 // ri=tanRMax*zi+rMaxAv ; 01682 // 01683 // Calculate a normal vector, as below 01684 01685 xi = p.x() + slentol*v.x(); 01686 yi = p.y() + slentol*v.y(); 01687 risec = std::sqrt(xi*xi + yi*yi)*secRMax; 01688 G4ThreeVector Normal = G4ThreeVector(xi/risec,yi/risec,-tanRMax/secRMax); 01689 01690 if ( Normal.dot(v) > 0 ) // We will leave the Cone immediatelly 01691 { 01692 if ( calcNorm ) 01693 { 01694 *n = Normal.unit() ; 01695 *validNorm = true ; 01696 } 01697 return snxt = 0.0 ; 01698 } 01699 else // On the surface, but not heading out so we ignore this intersection 01700 { // (as it is within tolerance). 01701 slentol = kInfinity ; 01702 } 01703 } 01704 01705 // Inner Cone intersection 01706 01707 if ( fRmin1 || fRmin2 ) 01708 { 01709 tanRMin = (fRmin2 - fRmin1)*0.5/fDz ; 01710 nt1 = t1 - (tanRMin*v.z())*(tanRMin*v.z()) ; 01711 01712 if ( nt1 ) 01713 { 01714 secRMin = std::sqrt(1.0 + tanRMin*tanRMin) ; 01715 rMinAv = (fRmin1 + fRmin2)*0.5 ; 01716 rin = tanRMin*p.z() + rMinAv ; 01717 nt2 = t2 - tanRMin*v.z()*rin ; 01718 nt3 = t3 - rin*rin ; 01719 01720 // Equation quadratic => 2 roots : first root must be leaving 01721 01722 b = nt2/nt1 ; 01723 c = nt3/nt1 ; 01724 d = b*b - c ; 01725 01726 if ( d >= 0.0 ) 01727 { 01728 // NOTE: should be rho-rin<kRadTolerance*0.5, 01729 // but using squared versions for efficiency 01730 01731 if (nt3 < kRadTolerance*(rin + kRadTolerance*0.25)) 01732 { 01733 if ( nt2 < 0.0 ) 01734 { 01735 if (calcNorm) { *validNorm = false; } 01736 return snxt = 0.0; 01737 } 01738 } 01739 else 01740 { 01741 if (b>0) { sr2 = -b - std::sqrt(d); } 01742 else { sr2 = c/(-b+std::sqrt(d)); } 01743 zi = p.z() + sr2*v.z() ; 01744 ri = tanRMin*zi + rMinAv ; 01745 01746 if( (ri>=0.0)&&(-halfRadTolerance<=sr2)&&(sr2<=halfRadTolerance) ) 01747 { 01748 // An intersection within the tolerance 01749 // storing it in case it is good. 01750 01751 slentol = sr2 ; 01752 sidetol = kRMax ; 01753 } 01754 if( (ri<0) || (sr2 < halfRadTolerance) ) 01755 { 01756 if (b>0) { sr3 = c/(-b-std::sqrt(d)); } 01757 else { sr3 = -b + std::sqrt(d) ; } 01758 01759 // Safety: if both roots -ve ensure that srd cannot `win' 01760 // distancetoout 01761 01762 if ( sr3 > halfRadTolerance ) 01763 { 01764 if( sr3 < srd ) 01765 { 01766 zi = p.z() + sr3*v.z() ; 01767 ri = tanRMin*zi + rMinAv ; 01768 01769 if ( ri >= 0.0 ) 01770 { 01771 srd=sr3 ; 01772 sider=kRMin ; 01773 } 01774 } 01775 } 01776 else if ( sr3 > -halfRadTolerance ) 01777 { 01778 // Intersection in tolerance. Store to check if it's good 01779 01780 slentol = sr3 ; 01781 sidetol = kRMin ; 01782 } 01783 } 01784 else if ( (sr2 < srd) && (sr2 > halfCarTolerance) ) 01785 { 01786 srd = sr2 ; 01787 sider = kRMin ; 01788 } 01789 else if (sr2 > -halfCarTolerance) 01790 { 01791 // Intersection in tolerance. Store to check if it's good 01792 01793 slentol = sr2 ; 01794 sidetol = kRMin ; 01795 } 01796 if( slentol <= halfCarTolerance ) 01797 { 01798 // An intersection within the tolerance was found. 01799 // We must accept it only if the momentum points outwards. 01800 01801 G4ThreeVector Normal ; 01802 01803 // Calculate a normal vector, as below 01804 01805 xi = p.x() + slentol*v.x() ; 01806 yi = p.y() + slentol*v.y() ; 01807 if( sidetol==kRMax ) 01808 { 01809 risec = std::sqrt(xi*xi + yi*yi)*secRMax ; 01810 Normal = G4ThreeVector(xi/risec,yi/risec,-tanRMax/secRMax) ; 01811 } 01812 else 01813 { 01814 risec = std::sqrt(xi*xi + yi*yi)*secRMin ; 01815 Normal = G4ThreeVector(-xi/risec,-yi/risec,tanRMin/secRMin) ; 01816 } 01817 if( Normal.dot(v) > 0 ) 01818 { 01819 // We will leave the cone immediately 01820 01821 if( calcNorm ) 01822 { 01823 *n = Normal.unit() ; 01824 *validNorm = true ; 01825 } 01826 return snxt = 0.0 ; 01827 } 01828 else 01829 { 01830 // On the surface, but not heading out so we ignore this 01831 // intersection (as it is within tolerance). 01832 01833 slentol = kInfinity ; 01834 } 01835 } 01836 } 01837 } 01838 } 01839 } 01840 01841 // Linear case => point outside inner cone ---> outer cone intersect 01842 // 01843 // Phi Intersection 01844 01845 if ( !fPhiFullCone ) 01846 { 01847 // add angle calculation with correction 01848 // of the difference in domain of atan2 and Sphi 01849 01850 vphi = std::atan2(v.y(),v.x()) ; 01851 01852 if ( vphi < fSPhi - halfAngTolerance ) { vphi += twopi; } 01853 else if ( vphi > fSPhi + fDPhi + halfAngTolerance ) { vphi -= twopi; } 01854 01855 if ( p.x() || p.y() ) // Check if on z axis (rho not needed later) 01856 { 01857 // pDist -ve when inside 01858 01859 pDistS = p.x()*sinSPhi - p.y()*cosSPhi ; 01860 pDistE = -p.x()*sinEPhi + p.y()*cosEPhi ; 01861 01862 // Comp -ve when in direction of outwards normal 01863 01864 compS = -sinSPhi*v.x() + cosSPhi*v.y() ; 01865 compE = sinEPhi*v.x() - cosEPhi*v.y() ; 01866 01867 sidephi = kNull ; 01868 01869 if( ( (fDPhi <= pi) && ( (pDistS <= halfCarTolerance) 01870 && (pDistE <= halfCarTolerance) ) ) 01871 || ( (fDPhi > pi) && !((pDistS > halfCarTolerance) 01872 && (pDistE > halfCarTolerance) ) ) ) 01873 { 01874 // Inside both phi *full* planes 01875 if ( compS < 0 ) 01876 { 01877 sphi = pDistS/compS ; 01878 if (sphi >= -halfCarTolerance) 01879 { 01880 xi = p.x() + sphi*v.x() ; 01881 yi = p.y() + sphi*v.y() ; 01882 01883 // Check intersecting with correct half-plane 01884 // (if not -> no intersect) 01885 // 01886 if ( (std::fabs(xi)<=kCarTolerance) 01887 && (std::fabs(yi)<=kCarTolerance) ) 01888 { 01889 sidephi= kSPhi; 01890 if ( ( fSPhi-halfAngTolerance <= vphi ) 01891 && ( fSPhi+fDPhi+halfAngTolerance >=vphi ) ) 01892 { 01893 sphi = kInfinity; 01894 } 01895 } 01896 else 01897 if ( (yi*cosCPhi-xi*sinCPhi)>=0 ) 01898 { 01899 sphi = kInfinity ; 01900 } 01901 else 01902 { 01903 sidephi = kSPhi ; 01904 if ( pDistS > -halfCarTolerance ) 01905 { 01906 sphi = 0.0 ; // Leave by sphi immediately 01907 } 01908 } 01909 } 01910 else 01911 { 01912 sphi = kInfinity ; 01913 } 01914 } 01915 else 01916 { 01917 sphi = kInfinity ; 01918 } 01919 01920 if ( compE < 0 ) 01921 { 01922 sphi2 = pDistE/compE ; 01923 01924 // Only check further if < starting phi intersection 01925 // 01926 if ( (sphi2 > -halfCarTolerance) && (sphi2 < sphi) ) 01927 { 01928 xi = p.x() + sphi2*v.x() ; 01929 yi = p.y() + sphi2*v.y() ; 01930 01931 // Check intersecting with correct half-plane 01932 01933 if ( (std::fabs(xi)<=kCarTolerance) 01934 && (std::fabs(yi)<=kCarTolerance) ) 01935 { 01936 // Leaving via ending phi 01937 01938 if(!( (fSPhi-halfAngTolerance <= vphi) 01939 && (fSPhi+fDPhi+halfAngTolerance >= vphi) ) ) 01940 { 01941 sidephi = kEPhi ; 01942 if ( pDistE <= -halfCarTolerance ) { sphi = sphi2; } 01943 else { sphi = 0.0; } 01944 } 01945 } 01946 else // Check intersecting with correct half-plane 01947 if ( yi*cosCPhi-xi*sinCPhi >= 0 ) 01948 { 01949 // Leaving via ending phi 01950 01951 sidephi = kEPhi ; 01952 if ( pDistE <= -halfCarTolerance ) { sphi = sphi2; } 01953 else { sphi = 0.0; } 01954 } 01955 } 01956 } 01957 } 01958 else 01959 { 01960 sphi = kInfinity ; 01961 } 01962 } 01963 else 01964 { 01965 // On z axis + travel not || to z axis -> if phi of vector direction 01966 // within phi of shape, Step limited by rmax, else Step =0 01967 01968 if ( (fSPhi-halfAngTolerance <= vphi) 01969 && (vphi <= fSPhi+fDPhi+halfAngTolerance) ) 01970 { 01971 sphi = kInfinity ; 01972 } 01973 else 01974 { 01975 sidephi = kSPhi ; // arbitrary 01976 sphi = 0.0 ; 01977 } 01978 } 01979 if ( sphi < snxt ) // Order intersecttions 01980 { 01981 snxt=sphi ; 01982 side=sidephi ; 01983 } 01984 } 01985 if ( srd < snxt ) // Order intersections 01986 { 01987 snxt = srd ; 01988 side = sider ; 01989 } 01990 if (calcNorm) 01991 { 01992 switch(side) 01993 { // Note: returned vector not normalised 01994 case kRMax: // (divide by frmax for unit vector) 01995 xi = p.x() + snxt*v.x() ; 01996 yi = p.y() + snxt*v.y() ; 01997 risec = std::sqrt(xi*xi + yi*yi)*secRMax ; 01998 *n = G4ThreeVector(xi/risec,yi/risec,-tanRMax/secRMax) ; 01999 *validNorm = true ; 02000 break ; 02001 case kRMin: 02002 *validNorm = false ; // Rmin is inconvex 02003 break ; 02004 case kSPhi: 02005 if ( fDPhi <= pi ) 02006 { 02007 *n = G4ThreeVector(sinSPhi, -cosSPhi, 0); 02008 *validNorm = true ; 02009 } 02010 else 02011 { 02012 *validNorm = false ; 02013 } 02014 break ; 02015 case kEPhi: 02016 if ( fDPhi <= pi ) 02017 { 02018 *n = G4ThreeVector(-sinEPhi, cosEPhi, 0); 02019 *validNorm = true ; 02020 } 02021 else 02022 { 02023 *validNorm = false ; 02024 } 02025 break ; 02026 case kPZ: 02027 *n = G4ThreeVector(0,0,1) ; 02028 *validNorm = true ; 02029 break ; 02030 case kMZ: 02031 *n = G4ThreeVector(0,0,-1) ; 02032 *validNorm = true ; 02033 break ; 02034 default: 02035 G4cout << G4endl ; 02036 DumpInfo(); 02037 std::ostringstream message; 02038 G4int oldprc = message.precision(16) ; 02039 message << "Undefined side for valid surface normal to solid." 02040 << G4endl 02041 << "Position:" << G4endl << G4endl 02042 << "p.x() = " << p.x()/mm << " mm" << G4endl 02043 << "p.y() = " << p.y()/mm << " mm" << G4endl 02044 << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl 02045 << "pho at z = " << std::sqrt( p.x()*p.x()+p.y()*p.y() )/mm 02046 << " mm" << G4endl << G4endl ; 02047 if( p.x() != 0. || p.x() != 0.) 02048 { 02049 message << "point phi = " << std::atan2(p.y(),p.x())/degree 02050 << " degree" << G4endl << G4endl ; 02051 } 02052 message << "Direction:" << G4endl << G4endl 02053 << "v.x() = " << v.x() << G4endl 02054 << "v.y() = " << v.y() << G4endl 02055 << "v.z() = " << v.z() << G4endl<< G4endl 02056 << "Proposed distance :" << G4endl<< G4endl 02057 << "snxt = " << snxt/mm << " mm" << G4endl ; 02058 message.precision(oldprc) ; 02059 G4Exception("G4Cons::DistanceToOut(p,v,..)","GeomSolids1002", 02060 JustWarning, message) ; 02061 break ; 02062 } 02063 } 02064 if (snxt < halfCarTolerance) { snxt = 0.; } 02065 02066 return snxt ; 02067 }
G4double G4Cons::GetCubicVolume | ( | ) | [inline, virtual] |
Reimplemented from G4VSolid.
Definition at line 254 of file G4Cons.icc.
References G4CSGSolid::fCubicVolume.
00255 { 00256 if(fCubicVolume != 0.) {;} 00257 else 00258 { 00259 G4double Rmean, rMean, deltaR, deltar; 00260 00261 Rmean = 0.5*(fRmax1+fRmax2); 00262 deltaR = fRmax1-fRmax2; 00263 00264 rMean = 0.5*(fRmin1+fRmin2); 00265 deltar = fRmin1-fRmin2; 00266 fCubicVolume = fDPhi*fDz*(Rmean*Rmean-rMean*rMean 00267 +(deltaR*deltaR-deltar*deltar)/12); 00268 } 00269 return fCubicVolume; 00270 }
G4double G4Cons::GetDeltaPhiAngle | ( | ) | const [inline] |
Definition at line 74 of file G4Cons.icc.
Referenced by G4HepRepSceneHandler::AddSolid(), G4HepRepFileSceneHandler::AddSolid(), G4ParameterisationConsZ::ComputeDimensions(), G4ParameterisationConsRho::ComputeDimensions(), G4GDMLWriteParamvol::Cone_dimensionsWrite(), G4GDMLWriteSolids::ConeWrite(), G4ParameterisationConsPhi::G4ParameterisationConsPhi(), G4VParameterisationCons::G4VParameterisationCons(), GetDPhi(), G4ParameterisationConsPhi::GetMaxParameter(), and G4tgbGeometryDumper::GetSolidParams().
G4double G4Cons::GetDPhi | ( | ) | const [inline] |
Definition at line 248 of file G4Cons.icc.
References GetDeltaPhiAngle().
00249 { 00250 return GetDeltaPhiAngle(); 00251 }
G4double G4Cons::GetDz | ( | ) | const [inline] |
Definition at line 236 of file G4Cons.icc.
References GetZHalfLength().
00237 { 00238 return GetZHalfLength(); 00239 }
G4GeometryType G4Cons::GetEntityType | ( | ) | const [virtual] |
G4double G4Cons::GetInnerRadiusMinusZ | ( | ) | const [inline] |
Definition at line 38 of file G4Cons.icc.
Referenced by G4HepRepSceneHandler::AddSolid(), G4HepRepFileSceneHandler::AddSolid(), G4ParameterisationConsZ::ComputeDimensions(), G4ParameterisationConsPhi::ComputeDimensions(), G4ParameterisationConsRho::ComputeDimensions(), G4GDMLWriteParamvol::Cone_dimensionsWrite(), G4GDMLWriteSolids::ConeWrite(), G4ParameterisationConsRho::G4ParameterisationConsRho(), G4VParameterisationCons::G4VParameterisationCons(), G4ParameterisationConsRho::GetMaxParameter(), GetRmin1(), and G4tgbGeometryDumper::GetSolidParams().
G4double G4Cons::GetInnerRadiusPlusZ | ( | ) | const [inline] |
Definition at line 50 of file G4Cons.icc.
Referenced by G4HepRepSceneHandler::AddSolid(), G4HepRepFileSceneHandler::AddSolid(), G4ParameterisationConsZ::ComputeDimensions(), G4ParameterisationConsPhi::ComputeDimensions(), G4ParameterisationConsRho::ComputeDimensions(), G4GDMLWriteParamvol::Cone_dimensionsWrite(), G4GDMLWriteSolids::ConeWrite(), G4ParameterisationConsRho::G4ParameterisationConsRho(), G4VParameterisationCons::G4VParameterisationCons(), GetRmin2(), and G4tgbGeometryDumper::GetSolidParams().
G4double G4Cons::GetOuterRadiusMinusZ | ( | ) | const [inline] |
Definition at line 44 of file G4Cons.icc.
Referenced by G4HepRepSceneHandler::AddSolid(), G4HepRepFileSceneHandler::AddSolid(), G4ParameterisationConsZ::ComputeDimensions(), G4ParameterisationConsPhi::ComputeDimensions(), G4GDMLWriteParamvol::Cone_dimensionsWrite(), G4GDMLWriteSolids::ConeWrite(), G4ParameterisationConsRho::G4ParameterisationConsRho(), G4VParameterisationCons::G4VParameterisationCons(), G4ParameterisationConsRho::GetMaxParameter(), GetRmax1(), and G4tgbGeometryDumper::GetSolidParams().
G4double G4Cons::GetOuterRadiusPlusZ | ( | ) | const [inline] |
Definition at line 56 of file G4Cons.icc.
Referenced by G4HepRepSceneHandler::AddSolid(), G4HepRepFileSceneHandler::AddSolid(), G4ParameterisationConsZ::ComputeDimensions(), G4ParameterisationConsPhi::ComputeDimensions(), G4ParameterisationConsRho::ComputeDimensions(), G4GDMLWriteParamvol::Cone_dimensionsWrite(), G4GDMLWriteSolids::ConeWrite(), G4VParameterisationCons::G4VParameterisationCons(), GetRmax2(), and G4tgbGeometryDumper::GetSolidParams().
G4ThreeVector G4Cons::GetPointOnSurface | ( | ) | const [virtual] |
Reimplemented from G4VSolid.
Definition at line 2290 of file G4Cons.cc.
References G4CSGSolid::GetRadiusInRing(), and sqr().
02291 { 02292 // declare working variables 02293 // 02294 G4double Aone, Atwo, Athree, Afour, Afive, slin, slout, phi; 02295 G4double zRand, cosu, sinu, rRand1, rRand2, chose, rone, rtwo, qone, qtwo; 02296 rone = (fRmax1-fRmax2)/(2.*fDz); 02297 rtwo = (fRmin1-fRmin2)/(2.*fDz); 02298 qone=0.; qtwo=0.; 02299 if(fRmax1!=fRmax2) { qone = fDz*(fRmax1+fRmax2)/(fRmax1-fRmax2); } 02300 if(fRmin1!=fRmin2) { qtwo = fDz*(fRmin1+fRmin2)/(fRmin1-fRmin2); } 02301 slin = std::sqrt(sqr(fRmin1-fRmin2)+sqr(2.*fDz)); 02302 slout = std::sqrt(sqr(fRmax1-fRmax2)+sqr(2.*fDz)); 02303 Aone = 0.5*fDPhi*(fRmax2 + fRmax1)*slout; 02304 Atwo = 0.5*fDPhi*(fRmin2 + fRmin1)*slin; 02305 Athree = 0.5*fDPhi*(fRmax1*fRmax1-fRmin1*fRmin1); 02306 Afour = 0.5*fDPhi*(fRmax2*fRmax2-fRmin2*fRmin2); 02307 Afive = fDz*(fRmax1-fRmin1+fRmax2-fRmin2); 02308 02309 phi = RandFlat::shoot(fSPhi,fSPhi+fDPhi); 02310 cosu = std::cos(phi); sinu = std::sin(phi); 02311 rRand1 = GetRadiusInRing(fRmin1, fRmin2); 02312 rRand2 = GetRadiusInRing(fRmax1, fRmax2); 02313 02314 if ( (fSPhi == 0.) && fPhiFullCone ) { Afive = 0.; } 02315 chose = RandFlat::shoot(0.,Aone+Atwo+Athree+Afour+2.*Afive); 02316 02317 if( (chose >= 0.) && (chose < Aone) ) 02318 { 02319 if(fRmin1 != fRmin2) 02320 { 02321 zRand = RandFlat::shoot(-1.*fDz,fDz); 02322 return G4ThreeVector (rtwo*cosu*(qtwo-zRand), 02323 rtwo*sinu*(qtwo-zRand), zRand); 02324 } 02325 else 02326 { 02327 return G4ThreeVector(fRmin1*cosu, fRmin2*sinu, 02328 RandFlat::shoot(-1.*fDz,fDz)); 02329 } 02330 } 02331 else if( (chose >= Aone) && (chose <= Aone + Atwo) ) 02332 { 02333 if(fRmax1 != fRmax2) 02334 { 02335 zRand = RandFlat::shoot(-1.*fDz,fDz); 02336 return G4ThreeVector (rone*cosu*(qone-zRand), 02337 rone*sinu*(qone-zRand), zRand); 02338 } 02339 else 02340 { 02341 return G4ThreeVector(fRmax1*cosu, fRmax2*sinu, 02342 RandFlat::shoot(-1.*fDz,fDz)); 02343 } 02344 } 02345 else if( (chose >= Aone + Atwo) && (chose < Aone + Atwo + Athree) ) 02346 { 02347 return G4ThreeVector (rRand1*cosu, rRand1*sinu, -1*fDz); 02348 } 02349 else if( (chose >= Aone + Atwo + Athree) 02350 && (chose < Aone + Atwo + Athree + Afour) ) 02351 { 02352 return G4ThreeVector (rRand2*cosu,rRand2*sinu,fDz); 02353 } 02354 else if( (chose >= Aone + Atwo + Athree + Afour) 02355 && (chose < Aone + Atwo + Athree + Afour + Afive) ) 02356 { 02357 zRand = RandFlat::shoot(-1.*fDz,fDz); 02358 rRand1 = RandFlat::shoot(fRmin2-((zRand-fDz)/(2.*fDz))*(fRmin1-fRmin2), 02359 fRmax2-((zRand-fDz)/(2.*fDz))*(fRmax1-fRmax2)); 02360 return G4ThreeVector (rRand1*std::cos(fSPhi), 02361 rRand1*std::sin(fSPhi), zRand); 02362 } 02363 else 02364 { 02365 zRand = RandFlat::shoot(-1.*fDz,fDz); 02366 rRand1 = RandFlat::shoot(fRmin2-((zRand-fDz)/(2.*fDz))*(fRmin1-fRmin2), 02367 fRmax2-((zRand-fDz)/(2.*fDz))*(fRmax1-fRmax2)); 02368 return G4ThreeVector (rRand1*std::cos(fSPhi+fDPhi), 02369 rRand1*std::sin(fSPhi+fDPhi), zRand); 02370 } 02371 }
G4double G4Cons::GetRmax1 | ( | ) | const [inline] |
Definition at line 218 of file G4Cons.icc.
References GetOuterRadiusMinusZ().
00219 { 00220 return GetOuterRadiusMinusZ(); 00221 }
G4double G4Cons::GetRmax2 | ( | ) | const [inline] |
Definition at line 230 of file G4Cons.icc.
References GetOuterRadiusPlusZ().
00231 { 00232 return GetOuterRadiusPlusZ(); 00233 }
G4double G4Cons::GetRmin1 | ( | ) | const [inline] |
Definition at line 212 of file G4Cons.icc.
References GetInnerRadiusMinusZ().
00213 { 00214 return GetInnerRadiusMinusZ(); 00215 }
G4double G4Cons::GetRmin2 | ( | ) | const [inline] |
Definition at line 224 of file G4Cons.icc.
References GetInnerRadiusPlusZ().
00225 { 00226 return GetInnerRadiusPlusZ(); 00227 }
G4double G4Cons::GetSPhi | ( | ) | const [inline] |
Definition at line 242 of file G4Cons.icc.
References GetStartPhiAngle().
00243 { 00244 return GetStartPhiAngle(); 00245 }
G4double G4Cons::GetStartPhiAngle | ( | ) | const [inline] |
Definition at line 68 of file G4Cons.icc.
Referenced by G4ParameterisationConsZ::ComputeDimensions(), G4ParameterisationConsPhi::ComputeDimensions(), G4ParameterisationConsRho::ComputeDimensions(), G4GDMLWriteParamvol::Cone_dimensionsWrite(), G4GDMLWriteSolids::ConeWrite(), G4VParameterisationCons::G4VParameterisationCons(), G4tgbGeometryDumper::GetSolidParams(), and GetSPhi().
G4double G4Cons::GetSurfaceArea | ( | ) | [inline, virtual] |
Reimplemented from G4VSolid.
Definition at line 273 of file G4Cons.icc.
References G4CSGSolid::fSurfaceArea.
00274 { 00275 if(fSurfaceArea != 0.) {;} 00276 else 00277 { 00278 G4double mmin, mmax, dmin, dmax; 00279 00280 mmin= (fRmin1+fRmin2)*0.5; 00281 mmax= (fRmax1+fRmax2)*0.5; 00282 dmin= (fRmin2-fRmin1); 00283 dmax= (fRmax2-fRmax1); 00284 00285 fSurfaceArea = fDPhi*( mmin * std::sqrt(dmin*dmin+4*fDz*fDz) 00286 + mmax * std::sqrt(dmax*dmax+4*fDz*fDz) 00287 + 0.5*(fRmax1*fRmax1-fRmin1*fRmin1 00288 +fRmax2*fRmax2-fRmin2*fRmin2 )); 00289 if(!fPhiFullCone) 00290 { 00291 fSurfaceArea = fSurfaceArea+4*fDz*(mmax-mmin); 00292 } 00293 } 00294 return fSurfaceArea; 00295 }
G4double G4Cons::GetZHalfLength | ( | ) | const [inline] |
Definition at line 62 of file G4Cons.icc.
Referenced by G4HepRepSceneHandler::AddSolid(), G4HepRepFileSceneHandler::AddSolid(), G4ParameterisationConsZ::ComputeDimensions(), G4ParameterisationConsPhi::ComputeDimensions(), G4ParameterisationConsRho::ComputeDimensions(), G4ParameterisationConsZ::ComputeTransformation(), G4GDMLWriteParamvol::Cone_dimensionsWrite(), G4GDMLWriteSolids::ConeWrite(), G4ParameterisationConsZ::G4ParameterisationConsZ(), G4VParameterisationCons::G4VParameterisationCons(), GetDz(), G4ParameterisationConsZ::GetMaxParameter(), and G4tgbGeometryDumper::GetSolidParams().
EInside G4Cons::Inside | ( | const G4ThreeVector & | p | ) | const [virtual] |
Implements G4VSolid.
Definition at line 191 of file G4Cons.cc.
References G4VSolid::kCarTolerance, kInside, kOutside, and kSurface.
Referenced by CalculateExtent(), and DistanceToOut().
00192 { 00193 G4double r2, rl, rh, pPhi, tolRMin, tolRMax; // rh2, rl2 ; 00194 EInside in; 00195 static const G4double halfCarTolerance=kCarTolerance*0.5; 00196 static const G4double halfRadTolerance=kRadTolerance*0.5; 00197 static const G4double halfAngTolerance=kAngTolerance*0.5; 00198 00199 if (std::fabs(p.z()) > fDz + halfCarTolerance ) { return in = kOutside; } 00200 else if(std::fabs(p.z()) >= fDz - halfCarTolerance ) { in = kSurface; } 00201 else { in = kInside; } 00202 00203 r2 = p.x()*p.x() + p.y()*p.y() ; 00204 rl = 0.5*(fRmin2*(p.z() + fDz) + fRmin1*(fDz - p.z()))/fDz ; 00205 rh = 0.5*(fRmax2*(p.z()+fDz)+fRmax1*(fDz-p.z()))/fDz; 00206 00207 // rh2 = rh*rh; 00208 00209 tolRMin = rl - halfRadTolerance; 00210 if ( tolRMin < 0 ) { tolRMin = 0; } 00211 tolRMax = rh + halfRadTolerance; 00212 00213 if ( (r2<tolRMin*tolRMin) || (r2>tolRMax*tolRMax) ) { return in = kOutside; } 00214 00215 if (rl) { tolRMin = rl + halfRadTolerance; } 00216 else { tolRMin = 0.0; } 00217 tolRMax = rh - halfRadTolerance; 00218 00219 if (in == kInside) // else it's kSurface already 00220 { 00221 if ( (r2 < tolRMin*tolRMin) || (r2 >= tolRMax*tolRMax) ) { in = kSurface; } 00222 } 00223 if ( !fPhiFullCone && ((p.x() != 0.0) || (p.y() != 0.0)) ) 00224 { 00225 pPhi = std::atan2(p.y(),p.x()) ; 00226 00227 if ( pPhi < fSPhi - halfAngTolerance ) { pPhi += twopi; } 00228 else if ( pPhi > fSPhi + fDPhi + halfAngTolerance ) { pPhi -= twopi; } 00229 00230 if ( (pPhi < fSPhi - halfAngTolerance) || 00231 (pPhi > fSPhi + fDPhi + halfAngTolerance) ) { return in = kOutside; } 00232 00233 else if (in == kInside) // else it's kSurface anyway already 00234 { 00235 if ( (pPhi < fSPhi + halfAngTolerance) || 00236 (pPhi > fSPhi + fDPhi - halfAngTolerance) ) { in = kSurface; } 00237 } 00238 } 00239 else if ( !fPhiFullCone ) { in = kSurface; } 00240 00241 return in ; 00242 }
Definition at line 161 of file G4Cons.cc.
References cosCPhi, cosEPhi, cosHDPhiIT, cosHDPhiOT, cosSPhi, fDPhi, fDz, fPhiFullCone, fRmax1, fRmax2, fRmin1, fRmin2, fSPhi, kAngTolerance, kRadTolerance, G4CSGSolid::operator=(), sinCPhi, sinEPhi, and sinSPhi.
00162 { 00163 // Check assignment to self 00164 // 00165 if (this == &rhs) { return *this; } 00166 00167 // Copy base class data 00168 // 00169 G4CSGSolid::operator=(rhs); 00170 00171 // Copy data 00172 // 00173 kRadTolerance = rhs.kRadTolerance; 00174 kAngTolerance = rhs.kAngTolerance; 00175 fRmin1 = rhs.fRmin1; fRmin2 = rhs.fRmin2; 00176 fRmax1 = rhs.fRmax1; fRmax2 = rhs.fRmax2; 00177 fDz = rhs.fDz; fSPhi = rhs.fSPhi; fDPhi = rhs.fDPhi; 00178 sinCPhi = rhs.sinCPhi; cosCPhi = rhs.cosCPhi; 00179 cosHDPhiOT = rhs.cosHDPhiOT; cosHDPhiIT = rhs.cosHDPhiIT; 00180 sinSPhi = rhs.sinSPhi; cosSPhi = rhs.cosSPhi; 00181 sinEPhi = rhs.sinEPhi; cosEPhi = rhs.cosEPhi; 00182 fPhiFullCone = rhs.fPhiFullCone; 00183 00184 return *this; 00185 }
void G4Cons::SetDeltaPhiAngle | ( | G4double | newDPhi | ) | [inline] |
Definition at line 203 of file G4Cons.icc.
Referenced by G4ParameterisationConsZ::ComputeDimensions(), G4ParameterisationConsPhi::ComputeDimensions(), and G4ParameterisationConsRho::ComputeDimensions().
void G4Cons::SetInnerRadiusMinusZ | ( | G4double | Rmin1 | ) | [inline] |
Definition at line 157 of file G4Cons.icc.
Referenced by G4ParameterisationConsZ::ComputeDimensions(), G4ParameterisationConsPhi::ComputeDimensions(), and G4ParameterisationConsRho::ComputeDimensions().
void G4Cons::SetInnerRadiusPlusZ | ( | G4double | Rmin2 | ) | [inline] |
Definition at line 171 of file G4Cons.icc.
Referenced by G4ParameterisationConsZ::ComputeDimensions(), G4ParameterisationConsPhi::ComputeDimensions(), and G4ParameterisationConsRho::ComputeDimensions().
void G4Cons::SetOuterRadiusMinusZ | ( | G4double | Rmax1 | ) | [inline] |
Definition at line 164 of file G4Cons.icc.
Referenced by G4ParameterisationConsZ::ComputeDimensions(), G4ParameterisationConsPhi::ComputeDimensions(), and G4ParameterisationConsRho::ComputeDimensions().
void G4Cons::SetOuterRadiusPlusZ | ( | G4double | Rmax2 | ) | [inline] |
Definition at line 178 of file G4Cons.icc.
Referenced by G4ParameterisationConsZ::ComputeDimensions(), G4ParameterisationConsPhi::ComputeDimensions(), and G4ParameterisationConsRho::ComputeDimensions().
Definition at line 192 of file G4Cons.icc.
Referenced by G4ParameterisationConsZ::ComputeDimensions(), G4ParameterisationConsPhi::ComputeDimensions(), and G4ParameterisationConsRho::ComputeDimensions().
00193 { 00194 // Flag 'compute' can be used to explicitely avoid recomputation of 00195 // trigonometry in case SetDeltaPhiAngle() is invoked afterwards 00196 00197 CheckSPhiAngle(newSPhi); 00198 fPhiFullCone = false; 00199 if (compute) { InitializeTrigonometry(); } 00200 Initialize(); 00201 }
void G4Cons::SetZHalfLength | ( | G4double | newDz | ) | [inline] |
Definition at line 185 of file G4Cons.icc.
Referenced by G4ParameterisationConsZ::ComputeDimensions(), G4ParameterisationConsPhi::ComputeDimensions(), and G4ParameterisationConsRho::ComputeDimensions().
std::ostream & G4Cons::StreamInfo | ( | std::ostream & | os | ) | const [virtual] |
Reimplemented from G4CSGSolid.
Definition at line 2263 of file G4Cons.cc.
References G4VSolid::GetName().
02264 { 02265 G4int oldprc = os.precision(16); 02266 os << "-----------------------------------------------------------\n" 02267 << " *** Dump for solid - " << GetName() << " ***\n" 02268 << " ===================================================\n" 02269 << " Solid type: G4Cons\n" 02270 << " Parameters: \n" 02271 << " inside -fDz radius: " << fRmin1/mm << " mm \n" 02272 << " outside -fDz radius: " << fRmax1/mm << " mm \n" 02273 << " inside +fDz radius: " << fRmin2/mm << " mm \n" 02274 << " outside +fDz radius: " << fRmax2/mm << " mm \n" 02275 << " half length in Z : " << fDz/mm << " mm \n" 02276 << " starting angle of segment: " << fSPhi/degree << " degrees \n" 02277 << " delta angle of segment : " << fDPhi/degree << " degrees \n" 02278 << "-----------------------------------------------------------\n"; 02279 os.precision(oldprc); 02280 02281 return os; 02282 }
G4ThreeVector G4Cons::SurfaceNormal | ( | const G4ThreeVector & | p | ) | const [virtual] |
Implements G4VSolid.
Definition at line 479 of file G4Cons.cc.
References G4Exception(), JustWarning, and G4VSolid::kCarTolerance.
00480 { 00481 G4int noSurfaces = 0; 00482 G4double rho, pPhi; 00483 G4double distZ, distRMin, distRMax; 00484 G4double distSPhi = kInfinity, distEPhi = kInfinity; 00485 G4double tanRMin, secRMin, pRMin, widRMin; 00486 G4double tanRMax, secRMax, pRMax, widRMax; 00487 00488 static const G4double delta = 0.5*kCarTolerance; 00489 static const G4double dAngle = 0.5*kAngTolerance; 00490 00491 G4ThreeVector norm, sumnorm(0.,0.,0.), nZ = G4ThreeVector(0.,0.,1.); 00492 G4ThreeVector nR, nr(0.,0.,0.), nPs, nPe; 00493 00494 distZ = std::fabs(std::fabs(p.z()) - fDz); 00495 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()); 00496 00497 tanRMin = (fRmin2 - fRmin1)*0.5/fDz; 00498 secRMin = std::sqrt(1 + tanRMin*tanRMin); 00499 pRMin = rho - p.z()*tanRMin; 00500 widRMin = fRmin2 - fDz*tanRMin; 00501 distRMin = std::fabs(pRMin - widRMin)/secRMin; 00502 00503 tanRMax = (fRmax2 - fRmax1)*0.5/fDz; 00504 secRMax = std::sqrt(1+tanRMax*tanRMax); 00505 pRMax = rho - p.z()*tanRMax; 00506 widRMax = fRmax2 - fDz*tanRMax; 00507 distRMax = std::fabs(pRMax - widRMax)/secRMax; 00508 00509 if (!fPhiFullCone) // Protected against (0,0,z) 00510 { 00511 if ( rho ) 00512 { 00513 pPhi = std::atan2(p.y(),p.x()); 00514 00515 if (pPhi < fSPhi-delta) { pPhi += twopi; } 00516 else if (pPhi > fSPhi+fDPhi+delta) { pPhi -= twopi; } 00517 00518 distSPhi = std::fabs( pPhi - fSPhi ); 00519 distEPhi = std::fabs( pPhi - fSPhi - fDPhi ); 00520 } 00521 else if( !(fRmin1) || !(fRmin2) ) 00522 { 00523 distSPhi = 0.; 00524 distEPhi = 0.; 00525 } 00526 nPs = G4ThreeVector(std::sin(fSPhi), -std::cos(fSPhi), 0); 00527 nPe = G4ThreeVector(-std::sin(fSPhi+fDPhi), std::cos(fSPhi+fDPhi), 0); 00528 } 00529 if ( rho > delta ) 00530 { 00531 nR = G4ThreeVector(p.x()/rho/secRMax, p.y()/rho/secRMax, -tanRMax/secRMax); 00532 if (fRmin1 || fRmin2) 00533 { 00534 nr = G4ThreeVector(-p.x()/rho/secRMin,-p.y()/rho/secRMin,tanRMin/secRMin); 00535 } 00536 } 00537 00538 if( distRMax <= delta ) 00539 { 00540 noSurfaces ++; 00541 sumnorm += nR; 00542 } 00543 if( (fRmin1 || fRmin2) && (distRMin <= delta) ) 00544 { 00545 noSurfaces ++; 00546 sumnorm += nr; 00547 } 00548 if( !fPhiFullCone ) 00549 { 00550 if (distSPhi <= dAngle) 00551 { 00552 noSurfaces ++; 00553 sumnorm += nPs; 00554 } 00555 if (distEPhi <= dAngle) 00556 { 00557 noSurfaces ++; 00558 sumnorm += nPe; 00559 } 00560 } 00561 if (distZ <= delta) 00562 { 00563 noSurfaces ++; 00564 if ( p.z() >= 0.) { sumnorm += nZ; } 00565 else { sumnorm -= nZ; } 00566 } 00567 if ( noSurfaces == 0 ) 00568 { 00569 #ifdef G4CSGDEBUG 00570 G4Exception("G4Cons::SurfaceNormal(p)", "GeomSolids1002", 00571 JustWarning, "Point p is not on surface !?" ); 00572 #endif 00573 norm = ApproxSurfaceNormal(p); 00574 } 00575 else if ( noSurfaces == 1 ) { norm = sumnorm; } 00576 else { norm = sumnorm.unit(); } 00577 00578 return norm ; 00579 }