#include <G4CutTubs.hh>
Inheritance diagram for G4CutTubs:
Definition at line 52 of file G4CutTubs.hh.
G4CutTubs::G4CutTubs | ( | const G4String & | pName, | |
G4double | pRMin, | |||
G4double | pRMax, | |||
G4double | pDz, | |||
G4double | pSPhi, | |||
G4double | pDPhi, | |||
G4ThreeVector | pLowNorm, | |||
G4ThreeVector | pHighNorm | |||
) |
Definition at line 66 of file G4CutTubs.cc.
References FatalException, G4endl, G4Exception(), G4GeometryTolerance::GetAngularTolerance(), G4GeometryTolerance::GetInstance(), G4VSolid::GetName(), G4GeometryTolerance::GetRadialTolerance(), IsCrossingCutPlanes(), JustWarning, G4Tubs::kAngTolerance, and G4Tubs::kRadTolerance.
Referenced by Clone().
00071 : G4Tubs(pName, pRMin, pRMax, pDz, pSPhi, pDPhi), 00072 fPhiFullCutTube(true) 00073 { 00074 kRadTolerance = G4GeometryTolerance::GetInstance()->GetRadialTolerance(); 00075 kAngTolerance = G4GeometryTolerance::GetInstance()->GetAngularTolerance(); 00076 00077 // Check on Cutted Planes Normals 00078 // If there is NO CUT, propose to use G4Tubs instead 00079 // 00080 if(pDPhi<twopi) { fPhiFullCutTube=false; } 00081 00082 if ( ( !pLowNorm.x()) && ( !pLowNorm.y()) 00083 && ( !pHighNorm.x()) && (!pHighNorm.y()) ) 00084 { 00085 std::ostringstream message; 00086 message << "Inexisting Low/High Normal to Z plane or Parallel to Z." 00087 << G4endl 00088 << "Normals to Z plane are (" << pLowNorm <<" and " 00089 << pHighNorm << ") in solid: " << GetName(); 00090 G4Exception("G4CutTubs::G4CutTubs()", "GeomSolids1001", 00091 JustWarning, message, "Should use G4Tubs!"); 00092 } 00093 00094 // If Normal is (0,0,0),means parallel to R, give it value of (0,0,+/-1) 00095 // 00096 if (pLowNorm.mag2() == 0.) { pLowNorm.setZ(-1.); } 00097 if (pHighNorm.mag2()== 0.) { pHighNorm.setZ(1.); } 00098 00099 // Given Normals to Cut Planes have to be an unit vectors. 00100 // Normalize if it is needed. 00101 // 00102 if (pLowNorm.mag2() != 1.) { pLowNorm = pLowNorm.unit(); } 00103 if (pHighNorm.mag2()!= 1.) { pHighNorm = pHighNorm.unit(); } 00104 00105 // Normals to cutted planes have to point outside Solid 00106 // 00107 if( (pLowNorm.mag2() != 0.) && (pHighNorm.mag2()!= 0. ) ) 00108 { 00109 if( ( pLowNorm.z()>= 0. ) || ( pHighNorm.z() <= 0.)) 00110 { 00111 std::ostringstream message; 00112 message << "Invalid Low or High Normal to Z plane; " 00113 "has to point outside Solid." << G4endl 00114 << "Invalid Norm to Z plane (" << pLowNorm << " or " 00115 << pHighNorm << ") in solid: " << GetName(); 00116 G4Exception("G4CutTubs::G4CutTubs()", "GeomSolids0002", 00117 FatalException, message); 00118 } 00119 } 00120 fLowNorm = pLowNorm; 00121 fHighNorm = pHighNorm; 00122 00123 // Check Intersection of Cutted planes. They MUST NOT Intersect 00124 // 00125 if(IsCrossingCutPlanes()) 00126 { 00127 std::ostringstream message; 00128 message << "Invalid Low or High Normal to Z plane; " 00129 << "Crossing Cutted Planes." << G4endl 00130 << "Invalid Norm to Z plane (" << pLowNorm << " and " 00131 << pHighNorm << ") in solid: " << GetName(); 00132 G4Exception("G4CutTubs::G4CutTubs()", "GeomSolids0002", 00133 FatalException, message); 00134 } 00135 }
G4CutTubs::~G4CutTubs | ( | ) |
G4CutTubs::G4CutTubs | ( | __void__ & | ) |
Definition at line 142 of file G4CutTubs.cc.
00143 : G4Tubs(a), 00144 fLowNorm(G4ThreeVector()), fHighNorm(G4ThreeVector()), 00145 fPhiFullCutTube(false) 00146 { 00147 }
G4CutTubs::G4CutTubs | ( | const G4CutTubs & | rhs | ) |
Definition at line 161 of file G4CutTubs.cc.
00162 : G4Tubs(rhs), 00163 fLowNorm(rhs.fLowNorm), fHighNorm(rhs.fHighNorm), 00164 fPhiFullCutTube(rhs.fPhiFullCutTube) 00165 { 00166 }
G4ThreeVector G4CutTubs::ApproxSurfaceNormal | ( | const G4ThreeVector & | p | ) | const [protected, virtual] |
Reimplemented from G4Tubs.
Definition at line 623 of file G4CutTubs.cc.
References G4VSolid::DumpInfo(), G4Tubs::fDPhi, G4Tubs::fDz, G4Tubs::fRMax, G4Tubs::fRMin, G4Tubs::fSPhi, G4Exception(), JustWarning, G4Tubs::kNEPhi, G4Tubs::kNRMax, G4Tubs::kNRMin, G4Tubs::kNSPhi, and G4Tubs::kNZ.
Referenced by SurfaceNormal().
00624 { 00625 ENorm side ; 00626 G4ThreeVector norm ; 00627 G4double rho, phi ; 00628 G4double distZLow,distZHigh,distZ; 00629 G4double distRMin, distRMax, distSPhi, distEPhi, distMin ; 00630 G4ThreeVector vZ=G4ThreeVector(0,0,fDz); 00631 00632 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ; 00633 00634 distRMin = std::fabs(rho - fRMin) ; 00635 distRMax = std::fabs(rho - fRMax) ; 00636 00637 //dist to Low Cut 00638 // 00639 distZLow =std::fabs((p+vZ).dot(fLowNorm)); 00640 00641 //dist to High Cut 00642 // 00643 distZHigh = std::fabs((p-vZ).dot(fHighNorm)); 00644 distZ=std::min(distZLow,distZHigh); 00645 00646 if (distRMin < distRMax) // First minimum 00647 { 00648 if ( distZ < distRMin ) 00649 { 00650 distMin = distZ ; 00651 side = kNZ ; 00652 } 00653 else 00654 { 00655 distMin = distRMin ; 00656 side = kNRMin ; 00657 } 00658 } 00659 else 00660 { 00661 if ( distZ < distRMax ) 00662 { 00663 distMin = distZ ; 00664 side = kNZ ; 00665 } 00666 else 00667 { 00668 distMin = distRMax ; 00669 side = kNRMax ; 00670 } 00671 } 00672 if (!fPhiFullCutTube && rho ) // Protected against (0,0,z) 00673 { 00674 phi = std::atan2(p.y(),p.x()) ; 00675 00676 if ( phi < 0 ) { phi += twopi; } 00677 00678 if ( fSPhi < 0 ) 00679 { 00680 distSPhi = std::fabs(phi - (fSPhi + twopi))*rho ; 00681 } 00682 else 00683 { 00684 distSPhi = std::fabs(phi - fSPhi)*rho ; 00685 } 00686 distEPhi = std::fabs(phi - fSPhi - fDPhi)*rho ; 00687 00688 if (distSPhi < distEPhi) // Find new minimum 00689 { 00690 if ( distSPhi < distMin ) 00691 { 00692 side = kNSPhi ; 00693 } 00694 } 00695 else 00696 { 00697 if ( distEPhi < distMin ) 00698 { 00699 side = kNEPhi ; 00700 } 00701 } 00702 } 00703 switch ( side ) 00704 { 00705 case kNRMin : // Inner radius 00706 { 00707 norm = G4ThreeVector(-p.x()/rho, -p.y()/rho, 0) ; 00708 break ; 00709 } 00710 case kNRMax : // Outer radius 00711 { 00712 norm = G4ThreeVector(p.x()/rho, p.y()/rho, 0) ; 00713 break ; 00714 } 00715 case kNZ : // + or - dz 00716 { 00717 if ( distZHigh > distZLow ) { norm = fHighNorm ; } 00718 else { norm = fLowNorm; } 00719 break ; 00720 } 00721 case kNSPhi: 00722 { 00723 norm = G4ThreeVector(std::sin(fSPhi), -std::cos(fSPhi), 0) ; 00724 break ; 00725 } 00726 case kNEPhi: 00727 { 00728 norm = G4ThreeVector(-std::sin(fSPhi+fDPhi), std::cos(fSPhi+fDPhi), 0) ; 00729 break; 00730 } 00731 default: // Should never reach this case ... 00732 { 00733 DumpInfo(); 00734 G4Exception("G4CutTubs::ApproxSurfaceNormal()", 00735 "GeomSolids1002", JustWarning, 00736 "Undefined side for valid surface normal to solid."); 00737 break ; 00738 } 00739 } 00740 return norm; 00741 }
G4bool G4CutTubs::CalculateExtent | ( | const EAxis | pAxis, | |
const G4VoxelLimits & | pVoxelLimit, | |||
const G4AffineTransform & | pTransform, | |||
G4double & | pmin, | |||
G4double & | pmax | |||
) | const [virtual] |
Reimplemented from G4Tubs.
Definition at line 194 of file G4CutTubs.cc.
References G4VSolid::ClipBetweenSections(), G4VSolid::ClipCrossSection(), CreateRotatedVertices(), G4Tubs::fDPhi, G4Tubs::fRMax, G4Tubs::fRMin, G4VoxelLimits::GetMaxExtent(), GetMaxMinZ(), 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().
00199 { 00200 if ( (!pTransform.IsRotated()) && (fDPhi == twopi) && (fRMin == 0) ) 00201 { 00202 // Special case handling for unrotated solid tubes 00203 // Compute x/y/z mins and maxs fro bounding box respecting limits, 00204 // with early returns if outside limits. Then switch() on pAxis, 00205 // and compute exact x and y limit for x/y case 00206 00207 G4double xoffset, xMin, xMax; 00208 G4double yoffset, yMin, yMax; 00209 G4double zoffset, zMin, zMax; 00210 00211 G4double diff1, diff2, maxDiff, newMin, newMax; 00212 G4double xoff1, xoff2, yoff1, yoff2, delta; 00213 00214 xoffset = pTransform.NetTranslation().x(); 00215 xMin = xoffset - fRMax; 00216 xMax = xoffset + fRMax; 00217 00218 if (pVoxelLimit.IsXLimited()) 00219 { 00220 if ( (xMin > pVoxelLimit.GetMaxXExtent()) 00221 || (xMax < pVoxelLimit.GetMinXExtent()) ) 00222 { 00223 return false; 00224 } 00225 else 00226 { 00227 if (xMin < pVoxelLimit.GetMinXExtent()) 00228 { 00229 xMin = pVoxelLimit.GetMinXExtent(); 00230 } 00231 if (xMax > pVoxelLimit.GetMaxXExtent()) 00232 { 00233 xMax = pVoxelLimit.GetMaxXExtent(); 00234 } 00235 } 00236 } 00237 yoffset = pTransform.NetTranslation().y(); 00238 yMin = yoffset - fRMax; 00239 yMax = yoffset + fRMax; 00240 00241 if ( pVoxelLimit.IsYLimited() ) 00242 { 00243 if ( (yMin > pVoxelLimit.GetMaxYExtent()) 00244 || (yMax < pVoxelLimit.GetMinYExtent()) ) 00245 { 00246 return false; 00247 } 00248 else 00249 { 00250 if (yMin < pVoxelLimit.GetMinYExtent()) 00251 { 00252 yMin = pVoxelLimit.GetMinYExtent(); 00253 } 00254 if (yMax > pVoxelLimit.GetMaxYExtent()) 00255 { 00256 yMax=pVoxelLimit.GetMaxYExtent(); 00257 } 00258 } 00259 } 00260 zoffset = pTransform.NetTranslation().z(); 00261 GetMaxMinZ(zMin,zMax); 00262 zMin += zoffset; 00263 zMax += zoffset; 00264 00265 if ( pVoxelLimit.IsZLimited() ) 00266 { 00267 if ( (zMin > pVoxelLimit.GetMaxZExtent()) 00268 || (zMax < pVoxelLimit.GetMinZExtent()) ) 00269 { 00270 return false; 00271 } 00272 else 00273 { 00274 if (zMin < pVoxelLimit.GetMinZExtent()) 00275 { 00276 zMin = pVoxelLimit.GetMinZExtent(); 00277 } 00278 if (zMax > pVoxelLimit.GetMaxZExtent()) 00279 { 00280 zMax = pVoxelLimit.GetMaxZExtent(); 00281 } 00282 } 00283 } 00284 switch ( pAxis ) // Known to cut cylinder 00285 { 00286 case kXAxis : 00287 { 00288 yoff1 = yoffset - yMin; 00289 yoff2 = yMax - yoffset; 00290 00291 if ( (yoff1 >= 0) && (yoff2 >= 0) ) // Y limits cross max/min x 00292 { // => no change 00293 pMin = xMin; 00294 pMax = xMax; 00295 } 00296 else 00297 { 00298 // Y limits don't cross max/min x => compute max delta x, 00299 // hence new mins/maxs 00300 00301 delta = fRMax*fRMax - yoff1*yoff1; 00302 diff1 = (delta>0.) ? std::sqrt(delta) : 0.; 00303 delta = fRMax*fRMax - yoff2*yoff2; 00304 diff2 = (delta>0.) ? std::sqrt(delta) : 0.; 00305 maxDiff = (diff1 > diff2) ? diff1:diff2; 00306 newMin = xoffset - maxDiff; 00307 newMax = xoffset + maxDiff; 00308 pMin = (newMin < xMin) ? xMin : newMin; 00309 pMax = (newMax > xMax) ? xMax : newMax; 00310 } 00311 break; 00312 } 00313 case kYAxis : 00314 { 00315 xoff1 = xoffset - xMin; 00316 xoff2 = xMax - xoffset; 00317 00318 if ( (xoff1 >= 0) && (xoff2 >= 0) ) // X limits cross max/min y 00319 { // => no change 00320 pMin = yMin; 00321 pMax = yMax; 00322 } 00323 else 00324 { 00325 // X limits don't cross max/min y => compute max delta y, 00326 // hence new mins/maxs 00327 00328 delta = fRMax*fRMax - xoff1*xoff1; 00329 diff1 = (delta>0.) ? std::sqrt(delta) : 0.; 00330 delta = fRMax*fRMax - xoff2*xoff2; 00331 diff2 = (delta>0.) ? std::sqrt(delta) : 0.; 00332 maxDiff = (diff1 > diff2) ? diff1 : diff2; 00333 newMin = yoffset - maxDiff; 00334 newMax = yoffset + maxDiff; 00335 pMin = (newMin < yMin) ? yMin : newMin; 00336 pMax = (newMax > yMax) ? yMax : newMax; 00337 } 00338 break; 00339 } 00340 case kZAxis: 00341 { 00342 pMin = zMin; 00343 pMax = zMax; 00344 break; 00345 } 00346 default: 00347 break; 00348 } 00349 pMin -= kCarTolerance; 00350 pMax += kCarTolerance; 00351 return true; 00352 } 00353 else // Calculate rotated vertex coordinates 00354 { 00355 G4int i, noEntries, noBetweenSections4; 00356 G4bool existsAfterClip = false; 00357 G4ThreeVectorList* vertices = CreateRotatedVertices(pTransform); 00358 00359 pMin = kInfinity; 00360 pMax = -kInfinity; 00361 00362 noEntries = vertices->size(); 00363 noBetweenSections4 = noEntries - 4; 00364 00365 for ( i = 0 ; i < noEntries ; i += 4 ) 00366 { 00367 ClipCrossSection(vertices, i, pVoxelLimit, pAxis, pMin, pMax); 00368 } 00369 for ( i = 0 ; i < noBetweenSections4 ; i += 4 ) 00370 { 00371 ClipBetweenSections(vertices, i, pVoxelLimit, pAxis, pMin, pMax); 00372 } 00373 if ( (pMin != kInfinity) || (pMax != -kInfinity) ) 00374 { 00375 existsAfterClip = true; 00376 pMin -= kCarTolerance; // Add 2*tolerance to avoid precision troubles 00377 pMax += kCarTolerance; 00378 } 00379 else 00380 { 00381 // Check for case where completely enveloping clipping volume 00382 // If point inside then we are confident that the solid completely 00383 // envelopes the clipping volume. Hence set min/max extents according 00384 // to clipping volume extents along the specified axis. 00385 00386 G4ThreeVector clipCentre( 00387 (pVoxelLimit.GetMinXExtent()+pVoxelLimit.GetMaxXExtent())*0.5, 00388 (pVoxelLimit.GetMinYExtent()+pVoxelLimit.GetMaxYExtent())*0.5, 00389 (pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5 ); 00390 00391 if ( Inside(pTransform.Inverse().TransformPoint(clipCentre)) != kOutside ) 00392 { 00393 existsAfterClip = true; 00394 pMin = pVoxelLimit.GetMinExtent(pAxis); 00395 pMax = pVoxelLimit.GetMaxExtent(pAxis); 00396 } 00397 } 00398 delete vertices; 00399 return existsAfterClip; 00400 } 00401 }
G4VSolid * G4CutTubs::Clone | ( | ) | const [virtual] |
Reimplemented from G4Tubs.
Definition at line 1862 of file G4CutTubs.cc.
References G4CutTubs().
01863 { 01864 return new G4CutTubs(*this); 01865 }
G4Polyhedron * G4CutTubs::CreatePolyhedron | ( | ) | const [virtual] |
Reimplemented from G4Tubs.
Definition at line 1975 of file G4CutTubs.cc.
References G4Tubs::CreatePolyhedron(), G4Tubs::fDz, GetCutZ(), G4VSolid::kCarTolerance, CLHEP::detail::n, and G4InuclParticleNames::nn.
01976 { 01977 typedef G4double G4double3[3]; 01978 typedef G4int G4int4[4]; 01979 01980 G4Polyhedron *ph = new G4Polyhedron; 01981 G4Polyhedron *ph1 = G4Tubs::CreatePolyhedron(); 01982 G4int nn=ph1->GetNoVertices(); 01983 G4int nf=ph1->GetNoFacets(); 01984 G4double3* xyz = new G4double3[nn]; // number of nodes 01985 G4int4* faces = new G4int4[nf] ; // number of faces 01986 01987 for(G4int i=0;i<nn;i++) 01988 { 01989 xyz[i][0]=ph1->GetVertex(i+1).x(); 01990 xyz[i][1]=ph1->GetVertex(i+1).y(); 01991 G4double tmpZ=ph1->GetVertex(i+1).z(); 01992 if(tmpZ>=fDz-kCarTolerance) 01993 { 01994 xyz[i][2]=GetCutZ(G4ThreeVector(xyz[i][0],xyz[i][1],fDz)); 01995 } 01996 else if(tmpZ<=-fDz+kCarTolerance) 01997 { 01998 xyz[i][2]=GetCutZ(G4ThreeVector(xyz[i][0],xyz[i][1],-fDz)); 01999 } 02000 else 02001 { 02002 xyz[i][2]=tmpZ; 02003 } 02004 } 02005 G4int iNodes[4]; 02006 G4int *iEdge=0; 02007 G4int n; 02008 for(G4int i=0;i<nf;i++) 02009 { 02010 ph1->GetFacet(i+1,n,iNodes,iEdge); 02011 for(G4int k=0;k<n;k++) 02012 { 02013 faces[i][k]=iNodes[k]; 02014 } 02015 for(G4int k=n;k<4;k++) 02016 { 02017 faces[i][k]=0; 02018 } 02019 } 02020 ph->createPolyhedron(nn,nf,xyz,faces); 02021 02022 delete [] xyz; 02023 delete [] faces; 02024 delete ph1; 02025 02026 return ph; 02027 }
G4ThreeVectorList * G4CutTubs::CreateRotatedVertices | ( | const G4AffineTransform & | pTransform | ) | const [protected] |
Reimplemented from G4Tubs.
Definition at line 1767 of file G4CutTubs.cc.
References G4VSolid::DumpInfo(), FatalException, G4Tubs::fDPhi, G4Tubs::fDz, G4Tubs::fRMax, G4Tubs::fRMin, G4Tubs::fSPhi, G4Exception(), GetCutZ(), G4VSolid::kCarTolerance, kMaxMeshSections, kMeshAngleDefault, kMinMeshSections, and G4AffineTransform::TransformPoint().
Referenced by CalculateExtent().
01768 { 01769 G4ThreeVectorList* vertices ; 01770 G4ThreeVector vertex0, vertex1, vertex2, vertex3 ; 01771 G4double meshAngle, meshRMax, crossAngle, 01772 cosCrossAngle, sinCrossAngle, sAngle; 01773 G4double rMaxX, rMaxY, rMinX, rMinY, meshRMin ; 01774 G4int crossSection, noCrossSections; 01775 01776 // Compute no of cross-sections necessary to mesh tube 01777 // 01778 noCrossSections = G4int(fDPhi/kMeshAngleDefault) + 1 ; 01779 01780 if ( noCrossSections < kMinMeshSections ) 01781 { 01782 noCrossSections = kMinMeshSections ; 01783 } 01784 else if (noCrossSections>kMaxMeshSections) 01785 { 01786 noCrossSections = kMaxMeshSections ; 01787 } 01788 // noCrossSections = 4 ; 01789 01790 meshAngle = fDPhi/(noCrossSections - 1) ; 01791 // meshAngle = fDPhi/(noCrossSections) ; 01792 01793 meshRMax = (fRMax+100*kCarTolerance)/std::cos(meshAngle*0.5) ; 01794 meshRMin = fRMin - 100*kCarTolerance ; 01795 01796 // If complete in phi, set start angle such that mesh will be at fRMax 01797 // on the x axis. Will give better extent calculations when not rotated. 01798 01799 if (fPhiFullCutTube && (fSPhi == 0) ) { sAngle = -meshAngle*0.5 ; } 01800 else { sAngle = fSPhi ; } 01801 01802 vertices = new G4ThreeVectorList(); 01803 01804 if ( vertices ) 01805 { 01806 vertices->reserve(noCrossSections*4); 01807 for (crossSection = 0 ; crossSection < noCrossSections ; crossSection++ ) 01808 { 01809 // Compute coordinates of cross section at section crossSection 01810 01811 crossAngle = sAngle + crossSection*meshAngle ; 01812 cosCrossAngle = std::cos(crossAngle) ; 01813 sinCrossAngle = std::sin(crossAngle) ; 01814 01815 rMaxX = meshRMax*cosCrossAngle ; 01816 rMaxY = meshRMax*sinCrossAngle ; 01817 01818 if(meshRMin <= 0.0) 01819 { 01820 rMinX = 0.0 ; 01821 rMinY = 0.0 ; 01822 } 01823 else 01824 { 01825 rMinX = meshRMin*cosCrossAngle ; 01826 rMinY = meshRMin*sinCrossAngle ; 01827 } 01828 vertex0 = G4ThreeVector(rMinX,rMinY,GetCutZ(G4ThreeVector(rMinX,rMinY,-fDz))) ; 01829 vertex1 = G4ThreeVector(rMaxX,rMaxY,GetCutZ(G4ThreeVector(rMaxX,rMaxY,-fDz))) ; 01830 vertex2 = G4ThreeVector(rMaxX,rMaxY,GetCutZ(G4ThreeVector(rMaxX,rMaxY,+fDz))) ; 01831 vertex3 = G4ThreeVector(rMinX,rMinY,GetCutZ(G4ThreeVector(rMinX,rMinY,+fDz))) ; 01832 01833 vertices->push_back(pTransform.TransformPoint(vertex0)) ; 01834 vertices->push_back(pTransform.TransformPoint(vertex1)) ; 01835 vertices->push_back(pTransform.TransformPoint(vertex2)) ; 01836 vertices->push_back(pTransform.TransformPoint(vertex3)) ; 01837 } 01838 } 01839 else 01840 { 01841 DumpInfo(); 01842 G4Exception("G4CutTubs::CreateRotatedVertices()", 01843 "GeomSolids0003", FatalException, 01844 "Error in allocation of vertices. Out of memory !"); 01845 } 01846 return vertices ; 01847 }
void G4CutTubs::DescribeYourselfTo | ( | G4VGraphicsScene & | scene | ) | const [virtual] |
Reimplemented from G4Tubs.
Definition at line 1970 of file G4CutTubs.cc.
References G4VGraphicsScene::AddSolid().
01971 { 01972 scene.AddSolid (*this) ; 01973 }
G4double G4CutTubs::DistanceToIn | ( | const G4ThreeVector & | p | ) | const [virtual] |
Reimplemented from G4Tubs.
Definition at line 1221 of file G4CutTubs.cc.
References G4Tubs::cosCPhi, G4Tubs::cosEPhi, G4Tubs::cosSPhi, G4Tubs::fDPhi, G4Tubs::fDz, G4Tubs::fRMax, G4Tubs::fRMin, G4Tubs::sinCPhi, G4Tubs::sinEPhi, and G4Tubs::sinSPhi.
01222 { 01223 G4double safRMin,safRMax,safZLow,safZHigh,safePhi,safe,rho,cosPsi; 01224 G4ThreeVector vZ=G4ThreeVector(0,0,fDz); 01225 01226 // Distance to R 01227 // 01228 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ; 01229 01230 safRMin = fRMin- rho ; 01231 safRMax = rho - fRMax ; 01232 01233 // Distances to ZCut(Low/High) 01234 01235 // Dist to Low Cut 01236 // 01237 safZLow = (p+vZ).dot(fLowNorm); 01238 01239 // Dist to High Cut 01240 // 01241 safZHigh = (p-vZ).dot(fHighNorm); 01242 01243 safe = std::max(safZLow,safZHigh); 01244 01245 if ( safRMin > safe ) { safe = safRMin; } 01246 if ( safRMax> safe ) { safe = safRMax; } 01247 01248 // Distance to Phi 01249 // 01250 if ( (!fPhiFullCutTube) && (rho) ) 01251 { 01252 // Psi=angle from central phi to point 01253 // 01254 cosPsi = (p.x()*cosCPhi + p.y()*sinCPhi)/rho ; 01255 01256 if ( cosPsi < std::cos(fDPhi*0.5) ) 01257 { 01258 // Point lies outside phi range 01259 01260 if ( (p.y()*cosCPhi - p.x()*sinCPhi) <= 0 ) 01261 { 01262 safePhi = std::fabs(p.x()*sinSPhi - p.y()*cosSPhi) ; 01263 } 01264 else 01265 { 01266 safePhi = std::fabs(p.x()*sinEPhi - p.y()*cosEPhi) ; 01267 } 01268 if ( safePhi > safe ) { safe = safePhi; } 01269 } 01270 } 01271 if ( safe < 0 ) { safe = 0; } 01272 01273 return safe ; 01274 }
G4double G4CutTubs::DistanceToIn | ( | const G4ThreeVector & | p, | |
const G4ThreeVector & | v | |||
) | const [virtual] |
Reimplemented from G4Tubs.
Definition at line 765 of file G4CutTubs.cc.
References G4Tubs::cosCPhi, G4Tubs::cosEPhi, G4Tubs::cosHDPhiIT, G4Tubs::cosSPhi, G4Tubs::fDz, G4Tubs::fRMax, G4Tubs::fRMin, GetCutZ(), G4VSolid::kCarTolerance, G4Tubs::kRadTolerance, G4Tubs::sinCPhi, G4Tubs::sinEPhi, and G4Tubs::sinSPhi.
00767 { 00768 G4double snxt = kInfinity ; // snxt = default return value 00769 G4double tolORMin2, tolIRMax2 ; // 'generous' radii squared 00770 G4double tolORMax2, tolIRMin2; 00771 const G4double dRmax = 100.*fRMax; 00772 G4ThreeVector vZ=G4ThreeVector(0,0,fDz); 00773 00774 static const G4double halfCarTolerance = 0.5*kCarTolerance; 00775 static const G4double halfRadTolerance = 0.5*kRadTolerance; 00776 00777 // Intersection point variables 00778 // 00779 G4double Dist, sd=0, xi, yi, zi, rho2, inum, iden, cosPsi, Comp,calf ; 00780 G4double t1, t2, t3, b, c, d ; // Quadratic solver variables 00781 G4double distZLow,distZHigh; 00782 // Calculate tolerant rmin and rmax 00783 00784 if (fRMin > kRadTolerance) 00785 { 00786 tolORMin2 = (fRMin - halfRadTolerance)*(fRMin - halfRadTolerance) ; 00787 tolIRMin2 = (fRMin + halfRadTolerance)*(fRMin + halfRadTolerance) ; 00788 } 00789 else 00790 { 00791 tolORMin2 = 0.0 ; 00792 tolIRMin2 = 0.0 ; 00793 } 00794 tolORMax2 = (fRMax + halfRadTolerance)*(fRMax + halfRadTolerance) ; 00795 tolIRMax2 = (fRMax - halfRadTolerance)*(fRMax - halfRadTolerance) ; 00796 00797 // Intersection with ZCut surfaces 00798 00799 // dist to Low Cut 00800 // 00801 distZLow =(p+vZ).dot(fLowNorm); 00802 00803 // dist to High Cut 00804 // 00805 distZHigh = (p-vZ).dot(fHighNorm); 00806 00807 if ( distZLow >= -halfCarTolerance ) 00808 { 00809 calf = v.dot(fLowNorm); 00810 if (calf<0) 00811 { 00812 sd = -distZLow/calf; 00813 if(sd < 0.0) { sd = 0.0; } 00814 00815 xi = p.x() + sd*v.x() ; // Intersection coords 00816 yi = p.y() + sd*v.y() ; 00817 rho2 = xi*xi + yi*yi ; 00818 00819 // Check validity of intersection 00820 00821 if ((tolIRMin2 <= rho2) && (rho2 <= tolIRMax2)) 00822 { 00823 if (!fPhiFullCutTube && rho2) 00824 { 00825 // Psi = angle made with central (average) phi of shape 00826 // 00827 inum = xi*cosCPhi + yi*sinCPhi ; 00828 iden = std::sqrt(rho2) ; 00829 cosPsi = inum/iden ; 00830 if (cosPsi >= cosHDPhiIT) { return sd ; } 00831 } 00832 else 00833 { 00834 return sd ; 00835 } 00836 } 00837 } 00838 else 00839 { 00840 if ( sd<halfCarTolerance ) 00841 { 00842 if(calf>=0) { sd=kInfinity; } 00843 return sd ; // On/outside extent, and heading away 00844 } // -> cannot intersect 00845 } 00846 } 00847 00848 if(distZHigh >= -halfCarTolerance ) 00849 { 00850 calf = v.dot(fHighNorm); 00851 if (calf<0) 00852 { 00853 sd = -distZHigh/calf; 00854 00855 if(sd < 0.0) { sd = 0.0; } 00856 00857 xi = p.x() + sd*v.x() ; // Intersection coords 00858 yi = p.y() + sd*v.y() ; 00859 rho2 = xi*xi + yi*yi ; 00860 00861 // Check validity of intersection 00862 00863 if ((tolIRMin2 <= rho2) && (rho2 <= tolIRMax2)) 00864 { 00865 if (!fPhiFullCutTube && rho2) 00866 { 00867 // Psi = angle made with central (average) phi of shape 00868 // 00869 inum = xi*cosCPhi + yi*sinCPhi ; 00870 iden = std::sqrt(rho2) ; 00871 cosPsi = inum/iden ; 00872 if (cosPsi >= cosHDPhiIT) { return sd ; } 00873 } 00874 else 00875 { 00876 return sd ; 00877 } 00878 } 00879 } 00880 else 00881 { 00882 if ( sd<halfCarTolerance ) 00883 { 00884 if(calf>=0) { sd=kInfinity; } 00885 return sd ; // On/outside extent, and heading away 00886 } // -> cannot intersect 00887 } 00888 } 00889 00890 // -> Can not intersect z surfaces 00891 // 00892 // Intersection with rmax (possible return) and rmin (must also check phi) 00893 // 00894 // Intersection point (xi,yi,zi) on line x=p.x+t*v.x etc. 00895 // 00896 // Intersects with x^2+y^2=R^2 00897 // 00898 // Hence (v.x^2+v.y^2)t^2+ 2t(p.x*v.x+p.y*v.y)+p.x^2+p.y^2-R^2=0 00899 // t1 t2 t3 00900 00901 t1 = 1.0 - v.z()*v.z() ; 00902 t2 = p.x()*v.x() + p.y()*v.y() ; 00903 t3 = p.x()*p.x() + p.y()*p.y() ; 00904 if ( t1 > 0 ) // Check not || to z axis 00905 { 00906 b = t2/t1 ; 00907 c = t3 - fRMax*fRMax ; 00908 00909 if ((t3 >= tolORMax2) && (t2<0)) // This also handles the tangent case 00910 { 00911 // Try outer cylinder intersection, c=(t3-fRMax*fRMax)/t1; 00912 00913 c /= t1 ; 00914 d = b*b - c ; 00915 00916 if (d >= 0) // If real root 00917 { 00918 sd = c/(-b+std::sqrt(d)); 00919 if (sd >= 0) // If 'forwards' 00920 { 00921 if ( sd>dRmax ) // Avoid rounding errors due to precision issues on 00922 { // 64 bits systems. Split long distances and recompute 00923 G4double fTerm = sd-std::fmod(sd,dRmax); 00924 sd = fTerm + DistanceToIn(p+fTerm*v,v); 00925 } 00926 // Check z intersection 00927 // 00928 zi = p.z() + sd*v.z() ; 00929 xi = p.x() + sd*v.x() ; 00930 yi = p.y() + sd*v.y() ; 00931 if ((-xi*fLowNorm.x()-yi*fLowNorm.y() 00932 -(zi+fDz)*fLowNorm.z())>-halfCarTolerance) 00933 { 00934 if ((-xi*fHighNorm.x()-yi*fHighNorm.y() 00935 +(fDz-zi)*fHighNorm.z())>-halfCarTolerance) 00936 { 00937 // Z ok. Check phi intersection if reqd 00938 // 00939 if (fPhiFullCutTube) 00940 { 00941 return sd ; 00942 } 00943 else 00944 { 00945 xi = p.x() + sd*v.x() ; 00946 yi = p.y() + sd*v.y() ; 00947 cosPsi = (xi*cosCPhi + yi*sinCPhi)/fRMax ; 00948 if (cosPsi >= cosHDPhiIT) { return sd ; } 00949 } 00950 } // end if std::fabs(zi) 00951 } 00952 } // end if (sd>=0) 00953 } // end if (d>=0) 00954 } // end if (r>=fRMax) 00955 else 00956 { 00957 // Inside outer radius : 00958 // check not inside, and heading through tubs (-> 0 to in) 00959 if ((t3 > tolIRMin2) && (t2 < 0) 00960 && (std::fabs(p.z()) <= std::fabs(GetCutZ(p))-halfCarTolerance )) 00961 { 00962 // Inside both radii, delta r -ve, inside z extent 00963 00964 if (!fPhiFullCutTube) 00965 { 00966 inum = p.x()*cosCPhi + p.y()*sinCPhi ; 00967 iden = std::sqrt(t3) ; 00968 cosPsi = inum/iden ; 00969 if (cosPsi >= cosHDPhiIT) 00970 { 00971 // In the old version, the small negative tangent for the point 00972 // on surface was not taken in account, and returning 0.0 ... 00973 // New version: check the tangent for the point on surface and 00974 // if no intersection, return kInfinity, if intersection instead 00975 // return sd. 00976 // 00977 c = t3-fRMax*fRMax; 00978 if ( c<=0.0 ) 00979 { 00980 return 0.0; 00981 } 00982 else 00983 { 00984 c = c/t1 ; 00985 d = b*b-c; 00986 if ( d>=0.0 ) 00987 { 00988 snxt = c/(-b+std::sqrt(d)); // using safe solution 00989 // for quadratic equation 00990 if ( snxt < halfCarTolerance ) { snxt=0; } 00991 return snxt ; 00992 } 00993 else 00994 { 00995 return kInfinity; 00996 } 00997 } 00998 } 00999 } 01000 else 01001 { 01002 // In the old version, the small negative tangent for the point 01003 // on surface was not taken in account, and returning 0.0 ... 01004 // New version: check the tangent for the point on surface and 01005 // if no intersection, return kInfinity, if intersection instead 01006 // return sd. 01007 // 01008 c = t3 - fRMax*fRMax; 01009 if ( c<=0.0 ) 01010 { 01011 return 0.0; 01012 } 01013 else 01014 { 01015 c = c/t1 ; 01016 d = b*b-c; 01017 if ( d>=0.0 ) 01018 { 01019 snxt= c/(-b+std::sqrt(d)); // using safe solution 01020 // for quadratic equation 01021 if ( snxt < halfCarTolerance ) { snxt=0; } 01022 return snxt ; 01023 } 01024 else 01025 { 01026 return kInfinity; 01027 } 01028 } 01029 } // end if (!fPhiFullCutTube) 01030 } // end if (t3>tolIRMin2) 01031 } // end if (Inside Outer Radius) 01032 01033 if ( fRMin ) // Try inner cylinder intersection 01034 { 01035 c = (t3 - fRMin*fRMin)/t1 ; 01036 d = b*b - c ; 01037 if ( d >= 0.0 ) // If real root 01038 { 01039 // Always want 2nd root - we are outside and know rmax Hit was bad 01040 // - If on surface of rmin also need farthest root 01041 01042 sd =( b > 0. )? c/(-b - std::sqrt(d)) : (-b + std::sqrt(d)); 01043 if (sd >= -10*halfCarTolerance) // check forwards 01044 { 01045 // Check z intersection 01046 // 01047 if (sd < 0.0) { sd = 0.0; } 01048 if (sd>dRmax) // Avoid rounding errors due to precision issues seen 01049 { // 64 bits systems. Split long distances and recompute 01050 G4double fTerm = sd-std::fmod(sd,dRmax); 01051 sd = fTerm + DistanceToIn(p+fTerm*v,v); 01052 } 01053 zi = p.z() + sd*v.z() ; 01054 xi = p.x() + sd*v.x() ; 01055 yi = p.y() + sd*v.y() ; 01056 if ((-xi*fLowNorm.x()-yi*fLowNorm.y() 01057 -(zi+fDz)*fLowNorm.z())>-halfCarTolerance) 01058 { 01059 if ((-xi*fHighNorm.x()-yi*fHighNorm.y() 01060 +(fDz-zi)*fHighNorm.z())>-halfCarTolerance) 01061 { 01062 // Z ok. Check phi 01063 // 01064 if ( fPhiFullCutTube ) 01065 { 01066 return sd ; 01067 } 01068 else 01069 { 01070 cosPsi = (xi*cosCPhi + yi*sinCPhi)/fRMin ; 01071 if (cosPsi >= cosHDPhiIT) 01072 { 01073 // Good inner radius isect 01074 // - but earlier phi isect still possible 01075 // 01076 snxt = sd ; 01077 } 01078 } 01079 } // end if std::fabs(zi) 01080 } 01081 } // end if (sd>=0) 01082 } // end if (d>=0) 01083 } // end if (fRMin) 01084 } 01085 01086 // Phi segment intersection 01087 // 01088 // o Tolerant of points inside phi planes by up to kCarTolerance*0.5 01089 // 01090 // o NOTE: Large duplication of code between sphi & ephi checks 01091 // -> only diffs: sphi -> ephi, Comp -> -Comp and half-plane 01092 // intersection check <=0 -> >=0 01093 // -> use some form of loop Construct ? 01094 // 01095 if ( !fPhiFullCutTube ) 01096 { 01097 // First phi surface (Starting phi) 01098 // 01099 Comp = v.x()*sinSPhi - v.y()*cosSPhi ; 01100 01101 if ( Comp < 0 ) // Component in outwards normal dirn 01102 { 01103 Dist = (p.y()*cosSPhi - p.x()*sinSPhi) ; 01104 01105 if ( Dist < halfCarTolerance ) 01106 { 01107 sd = Dist/Comp ; 01108 01109 if (sd < snxt) 01110 { 01111 if ( sd < 0 ) { sd = 0.0; } 01112 zi = p.z() + sd*v.z() ; 01113 xi = p.x() + sd*v.x() ; 01114 yi = p.y() + sd*v.y() ; 01115 if ((-xi*fLowNorm.x()-yi*fLowNorm.y() 01116 -(zi+fDz)*fLowNorm.z())>-halfCarTolerance) 01117 { 01118 if ((-xi*fHighNorm.x()-yi*fHighNorm.y() 01119 +(fDz-zi)*fHighNorm.z())>-halfCarTolerance) 01120 { 01121 rho2 = xi*xi + yi*yi ; 01122 if ( ( (rho2 >= tolIRMin2) && (rho2 <= tolIRMax2) ) 01123 || ( (rho2 > tolORMin2) && (rho2 < tolIRMin2) 01124 && ( v.y()*cosSPhi - v.x()*sinSPhi > 0 ) 01125 && ( v.x()*cosSPhi + v.y()*sinSPhi >= 0 ) ) 01126 || ( (rho2 > tolIRMax2) && (rho2 < tolORMax2) 01127 && (v.y()*cosSPhi - v.x()*sinSPhi > 0) 01128 && (v.x()*cosSPhi + v.y()*sinSPhi < 0) ) ) 01129 { 01130 // z and r intersections good 01131 // - check intersecting with correct half-plane 01132 // 01133 if ((yi*cosCPhi-xi*sinCPhi) <= halfCarTolerance) { snxt = sd; } 01134 } 01135 } //two Z conditions 01136 } 01137 } 01138 } 01139 } 01140 01141 // Second phi surface (Ending phi) 01142 // 01143 Comp = -(v.x()*sinEPhi - v.y()*cosEPhi) ; 01144 01145 if (Comp < 0 ) // Component in outwards normal dirn 01146 { 01147 Dist = -(p.y()*cosEPhi - p.x()*sinEPhi) ; 01148 01149 if ( Dist < halfCarTolerance ) 01150 { 01151 sd = Dist/Comp ; 01152 01153 if (sd < snxt) 01154 { 01155 if ( sd < 0 ) { sd = 0; } 01156 zi = p.z() + sd*v.z() ; 01157 xi = p.x() + sd*v.x() ; 01158 yi = p.y() + sd*v.y() ; 01159 if ((-xi*fLowNorm.x()-yi*fLowNorm.y() 01160 -(zi+fDz)*fLowNorm.z())>-halfCarTolerance) 01161 { 01162 if ((-xi*fHighNorm.x()-yi*fHighNorm.y() 01163 +(fDz-zi)*fHighNorm.z())>-halfCarTolerance) 01164 { 01165 xi = p.x() + sd*v.x() ; 01166 yi = p.y() + sd*v.y() ; 01167 rho2 = xi*xi + yi*yi ; 01168 if ( ( (rho2 >= tolIRMin2) && (rho2 <= tolIRMax2) ) 01169 || ( (rho2 > tolORMin2) && (rho2 < tolIRMin2) 01170 && (v.x()*sinEPhi - v.y()*cosEPhi > 0) 01171 && (v.x()*cosEPhi + v.y()*sinEPhi >= 0) ) 01172 || ( (rho2 > tolIRMax2) && (rho2 < tolORMax2) 01173 && (v.x()*sinEPhi - v.y()*cosEPhi > 0) 01174 && (v.x()*cosEPhi + v.y()*sinEPhi < 0) ) ) 01175 { 01176 // z and r intersections good 01177 // - check intersecting with correct half-plane 01178 // 01179 if ( (yi*cosCPhi-xi*sinCPhi) >= -halfCarTolerance ) 01180 { 01181 snxt = sd; 01182 } 01183 } //?? >=-halfCarTolerance 01184 } 01185 } // two Z conditions 01186 } 01187 } 01188 } // Comp < 0 01189 } // !fPhiFullTube 01190 if ( snxt<halfCarTolerance ) { snxt=0; } 01191 01192 return snxt ; 01193 }
G4double G4CutTubs::DistanceToOut | ( | const G4ThreeVector & | p | ) | const [virtual] |
Reimplemented from G4Tubs.
Definition at line 1714 of file G4CutTubs.cc.
References G4Tubs::cosCPhi, G4Tubs::cosEPhi, G4Tubs::cosSPhi, G4Tubs::fDz, G4Tubs::fRMax, G4Tubs::fRMin, G4Tubs::sinCPhi, G4Tubs::sinEPhi, and G4Tubs::sinSPhi.
01715 { 01716 G4double safRMin,safRMax,safZLow,safZHigh,safePhi,safe,rho; 01717 G4ThreeVector vZ=G4ThreeVector(0,0,fDz); 01718 01719 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()) ; // Distance to R 01720 01721 safRMin = rho - fRMin ; 01722 safRMax = fRMax - rho ; 01723 01724 // Distances to ZCut(Low/High) 01725 01726 // Dist to Low Cut 01727 // 01728 safZLow = std::fabs((p+vZ).dot(fLowNorm)); 01729 01730 // Dist to High Cut 01731 // 01732 safZHigh = std::fabs((p-vZ).dot(fHighNorm)); 01733 safe = std::min(safZLow,safZHigh); 01734 01735 if ( safRMin < safe ) { safe = safRMin; } 01736 if ( safRMax< safe ) { safe = safRMax; } 01737 01738 // Check if phi divided, Calc distances closest phi plane 01739 // 01740 if ( !fPhiFullCutTube ) 01741 { 01742 if ( p.y()*cosCPhi-p.x()*sinCPhi <= 0 ) 01743 { 01744 safePhi = -(p.x()*sinSPhi - p.y()*cosSPhi) ; 01745 } 01746 else 01747 { 01748 safePhi = (p.x()*sinEPhi - p.y()*cosEPhi) ; 01749 } 01750 if (safePhi < safe) { safe = safePhi ; } 01751 } 01752 if ( safe < 0 ) { safe = 0; } 01753 01754 return safe ; 01755 }
G4double G4CutTubs::DistanceToOut | ( | const G4ThreeVector & | p, | |
const G4ThreeVector & | v, | |||
const G4bool | calcNorm = G4bool(false) , |
|||
G4bool * | validNorm = 0 , |
|||
G4ThreeVector * | n = 0 | |||
) | const [virtual] |
Reimplemented from G4Tubs.
Definition at line 1281 of file G4CutTubs.cc.
References G4Tubs::cosCPhi, G4Tubs::cosEPhi, G4Tubs::cosSPhi, G4VSolid::DumpInfo(), G4Tubs::fDPhi, G4Tubs::fDz, G4Tubs::fRMax, G4Tubs::fRMin, G4Tubs::fSPhi, G4cout, G4endl, G4Exception(), JustWarning, G4Tubs::kAngTolerance, G4VSolid::kCarTolerance, G4Tubs::kEPhi, G4Tubs::kMZ, G4Tubs::kNull, G4Tubs::kPZ, G4Tubs::kRadTolerance, G4Tubs::kRMax, G4Tubs::kRMin, G4Tubs::kSPhi, G4INCL::Math::pi, G4Tubs::sinCPhi, G4Tubs::sinEPhi, and G4Tubs::sinSPhi.
01286 { 01287 ESide side=kNull , sider=kNull, sidephi=kNull ; 01288 G4double snxt=kInfinity, srd=kInfinity,sz=kInfinity, sphi=kInfinity ; 01289 G4double deltaR, t1, t2, t3, b, c, d2, roMin2 ; 01290 G4double distZLow,distZHigh,calfH,calfL; 01291 G4ThreeVector vZ=G4ThreeVector(0,0,fDz); 01292 static const G4double halfCarTolerance = kCarTolerance*0.5; 01293 static const G4double halfAngTolerance = kAngTolerance*0.5; 01294 01295 // Vars for phi intersection: 01296 // 01297 G4double pDistS, compS, pDistE, compE, sphi2, xi, yi, vphi, roi2 ; 01298 01299 // Z plane intersection 01300 // Distances to ZCut(Low/High) 01301 01302 // dist to Low Cut 01303 // 01304 distZLow =(p+vZ).dot(fLowNorm); 01305 01306 // dist to High Cut 01307 // 01308 distZHigh = (p-vZ).dot(fHighNorm); 01309 01310 calfH = v.dot(fHighNorm); 01311 calfL = v.dot(fLowNorm); 01312 01313 if (calfH > 0 ) 01314 { 01315 if ( distZHigh < halfCarTolerance ) 01316 { 01317 snxt = -distZHigh/calfH ; 01318 side = kPZ ; 01319 } 01320 else 01321 { 01322 if (calcNorm) 01323 { 01324 *n = G4ThreeVector(0,0,1) ; 01325 *validNorm = true ; 01326 } 01327 return snxt = 0 ; 01328 } 01329 } 01330 if ( calfL>0) 01331 { 01332 01333 if ( distZLow < halfCarTolerance ) 01334 { 01335 sz = -distZLow/calfL ; 01336 if(sz<snxt){ 01337 snxt=sz; 01338 side = kMZ ; 01339 } 01340 01341 } 01342 else 01343 { 01344 if (calcNorm) 01345 { 01346 *n = G4ThreeVector(0,0,-1) ; 01347 *validNorm = true ; 01348 } 01349 return snxt = 0.0 ; 01350 } 01351 } 01352 if((calfH<=0)&&(calfL<=0)) 01353 { 01354 snxt = kInfinity ; // Travel perpendicular to z axis 01355 side = kNull; 01356 } 01357 // Radial Intersections 01358 // 01359 // Find intersection with cylinders at rmax/rmin 01360 // Intersection point (xi,yi,zi) on line x=p.x+t*v.x etc. 01361 // 01362 // Intersects with x^2+y^2=R^2 01363 // 01364 // Hence (v.x^2+v.y^2)t^2+ 2t(p.x*v.x+p.y*v.y)+p.x^2+p.y^2-R^2=0 01365 // 01366 // t1 t2 t3 01367 01368 t1 = 1.0 - v.z()*v.z() ; // since v normalised 01369 t2 = p.x()*v.x() + p.y()*v.y() ; 01370 t3 = p.x()*p.x() + p.y()*p.y() ; 01371 01372 if ( snxt > 10*(fDz+fRMax) ) { roi2 = 2*fRMax*fRMax; } 01373 else { roi2 = snxt*snxt*t1 + 2*snxt*t2 + t3; } // radius^2 on +-fDz 01374 01375 if ( t1 > 0 ) // Check not parallel 01376 { 01377 // Calculate srd, r exit distance 01378 01379 if ( (t2 >= 0.0) && (roi2 > fRMax*(fRMax + kRadTolerance)) ) 01380 { 01381 // Delta r not negative => leaving via rmax 01382 01383 deltaR = t3 - fRMax*fRMax ; 01384 01385 // NOTE: Should use rho-fRMax<-kRadTolerance*0.5 01386 // - avoid sqrt for efficiency 01387 01388 if ( deltaR < -kRadTolerance*fRMax ) 01389 { 01390 b = t2/t1 ; 01391 c = deltaR/t1 ; 01392 d2 = b*b-c; 01393 if( d2 >= 0 ) { srd = c/( -b - std::sqrt(d2)); } 01394 else { srd = 0.; } 01395 sider = kRMax ; 01396 } 01397 else 01398 { 01399 // On tolerant boundary & heading outwards (or perpendicular to) 01400 // outer radial surface -> leaving immediately 01401 01402 if ( calcNorm ) 01403 { 01404 *n = G4ThreeVector(p.x()/fRMax,p.y()/fRMax,0) ; 01405 *validNorm = true ; 01406 } 01407 return snxt = 0 ; // Leaving by rmax immediately 01408 } 01409 } 01410 else if ( t2 < 0. ) // i.e. t2 < 0; Possible rmin intersection 01411 { 01412 roMin2 = t3 - t2*t2/t1 ; // min ro2 of the plane of movement 01413 01414 if ( fRMin && (roMin2 < fRMin*(fRMin - kRadTolerance)) ) 01415 { 01416 deltaR = t3 - fRMin*fRMin ; 01417 b = t2/t1 ; 01418 c = deltaR/t1 ; 01419 d2 = b*b - c ; 01420 01421 if ( d2 >= 0 ) // Leaving via rmin 01422 { 01423 // NOTE: SHould use rho-rmin>kRadTolerance*0.5 01424 // - avoid sqrt for efficiency 01425 01426 if (deltaR > kRadTolerance*fRMin) 01427 { 01428 srd = c/(-b+std::sqrt(d2)); 01429 sider = kRMin ; 01430 } 01431 else 01432 { 01433 if ( calcNorm ) { *validNorm = false; } // Concave side 01434 return snxt = 0.0; 01435 } 01436 } 01437 else // No rmin intersect -> must be rmax intersect 01438 { 01439 deltaR = t3 - fRMax*fRMax ; 01440 c = deltaR/t1 ; 01441 d2 = b*b-c; 01442 if( d2 >=0. ) 01443 { 01444 srd = -b + std::sqrt(d2) ; 01445 sider = kRMax ; 01446 } 01447 else // Case: On the border+t2<kRadTolerance 01448 // (v is perpendicular to the surface) 01449 { 01450 if (calcNorm) 01451 { 01452 *n = G4ThreeVector(p.x()/fRMax,p.y()/fRMax,0) ; 01453 *validNorm = true ; 01454 } 01455 return snxt = 0.0; 01456 } 01457 } 01458 } 01459 else if ( roi2 > fRMax*(fRMax + kRadTolerance) ) 01460 // No rmin intersect -> must be rmax intersect 01461 { 01462 deltaR = t3 - fRMax*fRMax ; 01463 b = t2/t1 ; 01464 c = deltaR/t1; 01465 d2 = b*b-c; 01466 if( d2 >= 0 ) 01467 { 01468 srd = -b + std::sqrt(d2) ; 01469 sider = kRMax ; 01470 } 01471 else // Case: On the border+t2<kRadTolerance 01472 // (v is perpendicular to the surface) 01473 { 01474 if (calcNorm) 01475 { 01476 *n = G4ThreeVector(p.x()/fRMax,p.y()/fRMax,0) ; 01477 *validNorm = true ; 01478 } 01479 return snxt = 0.0; 01480 } 01481 } 01482 } 01483 // Phi Intersection 01484 01485 if ( !fPhiFullCutTube ) 01486 { 01487 // add angle calculation with correction 01488 // of the difference in domain of atan2 and Sphi 01489 // 01490 vphi = std::atan2(v.y(),v.x()) ; 01491 01492 if ( vphi < fSPhi - halfAngTolerance ) { vphi += twopi; } 01493 else if ( vphi > fSPhi + fDPhi + halfAngTolerance ) { vphi -= twopi; } 01494 01495 01496 if ( p.x() || p.y() ) // Check if on z axis (rho not needed later) 01497 { 01498 // pDist -ve when inside 01499 01500 pDistS = p.x()*sinSPhi - p.y()*cosSPhi ; 01501 pDistE = -p.x()*sinEPhi + p.y()*cosEPhi ; 01502 01503 // Comp -ve when in direction of outwards normal 01504 01505 compS = -sinSPhi*v.x() + cosSPhi*v.y() ; 01506 compE = sinEPhi*v.x() - cosEPhi*v.y() ; 01507 01508 sidephi = kNull; 01509 01510 if( ( (fDPhi <= pi) && ( (pDistS <= halfCarTolerance) 01511 && (pDistE <= halfCarTolerance) ) ) 01512 || ( (fDPhi > pi) && !((pDistS > halfCarTolerance) 01513 && (pDistE > halfCarTolerance) ) ) ) 01514 { 01515 // Inside both phi *full* planes 01516 01517 if ( compS < 0 ) 01518 { 01519 sphi = pDistS/compS ; 01520 01521 if (sphi >= -halfCarTolerance) 01522 { 01523 xi = p.x() + sphi*v.x() ; 01524 yi = p.y() + sphi*v.y() ; 01525 01526 // Check intersecting with correct half-plane 01527 // (if not -> no intersect) 01528 // 01529 if( (std::fabs(xi)<=kCarTolerance) 01530 && (std::fabs(yi)<=kCarTolerance) ) 01531 { 01532 sidephi = kSPhi; 01533 if (((fSPhi-halfAngTolerance)<=vphi) 01534 &&((fSPhi+fDPhi+halfAngTolerance)>=vphi)) 01535 { 01536 sphi = kInfinity; 01537 } 01538 } 01539 else if ( yi*cosCPhi-xi*sinCPhi >=0 ) 01540 { 01541 sphi = kInfinity ; 01542 } 01543 else 01544 { 01545 sidephi = kSPhi ; 01546 if ( pDistS > -halfCarTolerance ) 01547 { 01548 sphi = 0.0 ; // Leave by sphi immediately 01549 } 01550 } 01551 } 01552 else 01553 { 01554 sphi = kInfinity ; 01555 } 01556 } 01557 else 01558 { 01559 sphi = kInfinity ; 01560 } 01561 01562 if ( compE < 0 ) 01563 { 01564 sphi2 = pDistE/compE ; 01565 01566 // Only check further if < starting phi intersection 01567 // 01568 if ( (sphi2 > -halfCarTolerance) && (sphi2 < sphi) ) 01569 { 01570 xi = p.x() + sphi2*v.x() ; 01571 yi = p.y() + sphi2*v.y() ; 01572 01573 if ((std::fabs(xi)<=kCarTolerance)&&(std::fabs(yi)<=kCarTolerance)) 01574 { 01575 // Leaving via ending phi 01576 // 01577 if( !((fSPhi-halfAngTolerance <= vphi) 01578 &&(fSPhi+fDPhi+halfAngTolerance >= vphi)) ) 01579 { 01580 sidephi = kEPhi ; 01581 if ( pDistE <= -halfCarTolerance ) { sphi = sphi2 ; } 01582 else { sphi = 0.0 ; } 01583 } 01584 } 01585 else // Check intersecting with correct half-plane 01586 01587 if ( (yi*cosCPhi-xi*sinCPhi) >= 0) 01588 { 01589 // Leaving via ending phi 01590 // 01591 sidephi = kEPhi ; 01592 if ( pDistE <= -halfCarTolerance ) { sphi = sphi2 ; } 01593 else { sphi = 0.0 ; } 01594 } 01595 } 01596 } 01597 } 01598 else 01599 { 01600 sphi = kInfinity ; 01601 } 01602 } 01603 else 01604 { 01605 // On z axis + travel not || to z axis -> if phi of vector direction 01606 // within phi of shape, Step limited by rmax, else Step =0 01607 01608 if ( (fSPhi - halfAngTolerance <= vphi) 01609 && (vphi <= fSPhi + fDPhi + halfAngTolerance ) ) 01610 { 01611 sphi = kInfinity ; 01612 } 01613 else 01614 { 01615 sidephi = kSPhi ; // arbitrary 01616 sphi = 0.0 ; 01617 } 01618 } 01619 if (sphi < snxt) // Order intersecttions 01620 { 01621 snxt = sphi ; 01622 side = sidephi ; 01623 } 01624 } 01625 if (srd < snxt) // Order intersections 01626 { 01627 snxt = srd ; 01628 side = sider ; 01629 } 01630 } 01631 if (calcNorm) 01632 { 01633 switch(side) 01634 { 01635 case kRMax: 01636 // Note: returned vector not normalised 01637 // (divide by fRMax for unit vector) 01638 // 01639 xi = p.x() + snxt*v.x() ; 01640 yi = p.y() + snxt*v.y() ; 01641 *n = G4ThreeVector(xi/fRMax,yi/fRMax,0) ; 01642 *validNorm = true ; 01643 break ; 01644 01645 case kRMin: 01646 *validNorm = false ; // Rmin is inconvex 01647 break ; 01648 01649 case kSPhi: 01650 if ( fDPhi <= pi ) 01651 { 01652 *n = G4ThreeVector(sinSPhi,-cosSPhi,0) ; 01653 *validNorm = true ; 01654 } 01655 else 01656 { 01657 *validNorm = false ; 01658 } 01659 break ; 01660 01661 case kEPhi: 01662 if (fDPhi <= pi) 01663 { 01664 *n = G4ThreeVector(-sinEPhi,cosEPhi,0) ; 01665 *validNorm = true ; 01666 } 01667 else 01668 { 01669 *validNorm = false ; 01670 } 01671 break ; 01672 01673 case kPZ: 01674 *n = fHighNorm ; 01675 *validNorm = true ; 01676 break ; 01677 01678 case kMZ: 01679 *n = fLowNorm ; 01680 *validNorm = true ; 01681 break ; 01682 01683 default: 01684 G4cout << G4endl ; 01685 DumpInfo(); 01686 std::ostringstream message; 01687 G4int oldprc = message.precision(16); 01688 message << "Undefined side for valid surface normal to solid." 01689 << G4endl 01690 << "Position:" << G4endl << G4endl 01691 << "p.x() = " << p.x()/mm << " mm" << G4endl 01692 << "p.y() = " << p.y()/mm << " mm" << G4endl 01693 << "p.z() = " << p.z()/mm << " mm" << G4endl << G4endl 01694 << "Direction:" << G4endl << G4endl 01695 << "v.x() = " << v.x() << G4endl 01696 << "v.y() = " << v.y() << G4endl 01697 << "v.z() = " << v.z() << G4endl << G4endl 01698 << "Proposed distance :" << G4endl << G4endl 01699 << "snxt = " << snxt/mm << " mm" << G4endl ; 01700 message.precision(oldprc) ; 01701 G4Exception("G4CutTubs::DistanceToOut(p,v,..)", "GeomSolids1002", 01702 JustWarning, message); 01703 break ; 01704 } 01705 } 01706 if ( snxt<halfCarTolerance ) { snxt=0 ; } 01707 return snxt ; 01708 }
G4double G4CutTubs::GetCubicVolume | ( | ) | [inline, virtual] |
Reimplemented from G4Tubs.
Definition at line 38 of file G4CutTubs.icc.
References G4VSolid::GetCubicVolume().
00039 { 00040 return G4VSolid::GetCubicVolume(); 00041 }
G4double G4CutTubs::GetCutZ | ( | const G4ThreeVector & | p | ) | const [protected] |
Definition at line 2058 of file G4CutTubs.cc.
References G4Tubs::fDz.
Referenced by CreatePolyhedron(), CreateRotatedVertices(), DistanceToIn(), GetMaxMinZ(), GetPointOnSurface(), and IsCrossingCutPlanes().
02059 { 02060 G4double newz = p.z(); // p.z() should be either +fDz or -fDz 02061 if (p.z()<0) 02062 { 02063 if(fLowNorm.z()!=0.) 02064 { 02065 newz = -fDz-(p.x()*fLowNorm.x()+p.y()*fLowNorm.y())/fLowNorm.z(); 02066 } 02067 } 02068 else 02069 { 02070 if(fHighNorm.z()!=0.) 02071 { 02072 newz = fDz-(p.x()*fHighNorm.x()+p.y()*fHighNorm.y())/fHighNorm.z(); 02073 } 02074 } 02075 return newz; 02076 }
G4GeometryType G4CutTubs::GetEntityType | ( | ) | const [virtual] |
Reimplemented from G4Tubs.
Definition at line 1853 of file G4CutTubs.cc.
01854 { 01855 return G4String("G4CutTubs"); 01856 }
G4ThreeVector G4CutTubs::GetHighNorm | ( | ) | const [inline] |
G4ThreeVector G4CutTubs::GetLowNorm | ( | ) | const [inline] |
Definition at line 2082 of file G4CutTubs.cc.
References G4Tubs::fDPhi, G4Tubs::fDz, G4Tubs::fRMax, G4Tubs::fRMin, G4Tubs::fSPhi, GetCutZ(), and G4INCL::Math::pi.
Referenced by CalculateExtent().
02084 { 02085 G4double phiLow = std::atan2(fLowNorm.y(),fLowNorm.x()); 02086 G4double phiHigh= std::atan2(fHighNorm.y(),fHighNorm.x()); 02087 02088 G4double xc=0, yc=0,z1; 02089 G4double z[8]; 02090 G4bool in_range_low = false; 02091 G4bool in_range_hi = false; 02092 02093 G4int i; 02094 for (i=0; i<2; i++) 02095 { 02096 if (phiLow<0) { phiLow+=twopi; } 02097 G4double ddp = phiLow-fSPhi; 02098 if (ddp<0) { ddp += twopi; } 02099 if (ddp <= fDPhi) 02100 { 02101 xc = fRMin*std::cos(phiLow); 02102 yc = fRMin*std::sin(phiLow); 02103 z1 = GetCutZ(G4ThreeVector(xc, yc, -fDz)); 02104 xc = fRMax*std::cos(phiLow); 02105 yc = fRMax*std::sin(phiLow); 02106 z1 = std::min(z1, GetCutZ(G4ThreeVector(xc, yc, -fDz))); 02107 if (in_range_low) { zmin = std::min(zmin, z1); } 02108 else { zmin = z1; } 02109 in_range_low = true; 02110 } 02111 phiLow += pi; 02112 if (phiLow>twopi) { phiLow-=twopi; } 02113 } 02114 for (i=0; i<2; i++) 02115 { 02116 if (phiHigh<0) { phiHigh+=twopi; } 02117 G4double ddp = phiHigh-fSPhi; 02118 if (ddp<0) { ddp += twopi; } 02119 if (ddp <= fDPhi) 02120 { 02121 xc = fRMin*std::cos(phiHigh); 02122 yc = fRMin*std::sin(phiHigh); 02123 z1 = GetCutZ(G4ThreeVector(xc, yc, fDz)); 02124 xc = fRMax*std::cos(phiHigh); 02125 yc = fRMax*std::sin(phiHigh); 02126 z1 = std::min(z1, GetCutZ(G4ThreeVector(xc, yc, fDz))); 02127 if (in_range_hi) { zmax = std::min(zmax, z1); } 02128 else { zmax = z1; } 02129 in_range_hi = true; 02130 } 02131 phiHigh += pi; 02132 if (phiLow>twopi) { phiHigh-=twopi; } 02133 } 02134 02135 xc = fRMin*std::cos(fSPhi); 02136 yc = fRMin*std::sin(fSPhi); 02137 z[0] = GetCutZ(G4ThreeVector(xc, yc, -fDz)); 02138 z[4] = GetCutZ(G4ThreeVector(xc, yc, fDz)); 02139 02140 xc = fRMin*std::cos(fSPhi+fDPhi); 02141 yc = fRMin*std::sin(fSPhi+fDPhi); 02142 z[1] = GetCutZ(G4ThreeVector(xc, yc, -fDz)); 02143 z[5] = GetCutZ(G4ThreeVector(xc, yc, fDz)); 02144 02145 xc = fRMax*std::cos(fSPhi); 02146 yc = fRMax*std::sin(fSPhi); 02147 z[2] = GetCutZ(G4ThreeVector(xc, yc, -fDz)); 02148 z[6] = GetCutZ(G4ThreeVector(xc, yc, fDz)); 02149 02150 xc = fRMax*std::cos(fSPhi+fDPhi); 02151 yc = fRMax*std::sin(fSPhi+fDPhi); 02152 z[3] = GetCutZ(G4ThreeVector(xc, yc, -fDz)); 02153 z[7] = GetCutZ(G4ThreeVector(xc, yc, fDz)); 02154 02155 // Find min/max 02156 02157 z1=z[0]; 02158 for (i = 1; i < 4; i++) 02159 { 02160 if(z[i] < z[i-1])z1=z[i]; 02161 } 02162 02163 if (in_range_low) 02164 { 02165 zmin = std::min(zmin, z1); 02166 } 02167 else 02168 { 02169 zmin = z1; 02170 } 02171 z1=z[4]; 02172 for (i = 1; i < 4; i++) 02173 { 02174 if(z[4+i] > z[4+i-1]) { z1=z[4+i]; } 02175 } 02176 02177 if (in_range_hi) { zmax = std::max(zmax, z1); } 02178 else { zmax = z1; } 02179 }
G4ThreeVector G4CutTubs::GetPointOnSurface | ( | ) | const [virtual] |
Reimplemented from G4Tubs.
Definition at line 1896 of file G4CutTubs.cc.
References G4Tubs::fDPhi, G4Tubs::fDz, G4Tubs::fRMax, G4Tubs::fRMin, G4Tubs::fSPhi, GetCutZ(), and G4CSGSolid::GetRadiusInRing().
01897 { 01898 G4double xRand, yRand, zRand, phi, cosphi, sinphi, chose, 01899 aOne, aTwo, aThr, aFou; 01900 G4double rRand; 01901 01902 aOne = 2.*fDz*fDPhi*fRMax; 01903 aTwo = 2.*fDz*fDPhi*fRMin; 01904 aThr = 0.5*fDPhi*(fRMax*fRMax-fRMin*fRMin); 01905 aFou = 2.*fDz*(fRMax-fRMin); 01906 01907 phi = RandFlat::shoot(fSPhi, fSPhi+fDPhi); 01908 cosphi = std::cos(phi); 01909 sinphi = std::sin(phi); 01910 01911 rRand = GetRadiusInRing(fRMin,fRMax); 01912 01913 if( (fSPhi == 0) && (fDPhi == twopi) ) { aFou = 0; } 01914 01915 chose = RandFlat::shoot(0.,aOne+aTwo+2.*aThr+2.*aFou); 01916 01917 if( (chose >=0) && (chose < aOne) ) 01918 { 01919 xRand = fRMax*cosphi; 01920 yRand = fRMax*sinphi; 01921 zRand = RandFlat::shoot(GetCutZ(G4ThreeVector(xRand,yRand,-fDz)), 01922 GetCutZ(G4ThreeVector(xRand,yRand,fDz))); 01923 return G4ThreeVector (xRand, yRand, zRand); 01924 } 01925 else if( (chose >= aOne) && (chose < aOne + aTwo) ) 01926 { 01927 xRand = fRMin*cosphi; 01928 yRand = fRMin*sinphi; 01929 zRand = RandFlat::shoot(GetCutZ(G4ThreeVector(xRand,yRand,-fDz)), 01930 GetCutZ(G4ThreeVector(xRand,yRand,fDz))); 01931 return G4ThreeVector (xRand, yRand, zRand); 01932 } 01933 else if( (chose >= aOne + aTwo) && (chose < aOne + aTwo + aThr) ) 01934 { 01935 xRand = rRand*cosphi; 01936 yRand = rRand*sinphi; 01937 zRand = GetCutZ(G4ThreeVector(xRand,yRand,fDz)); 01938 return G4ThreeVector (xRand, yRand, zRand); 01939 } 01940 else if( (chose >= aOne + aTwo + aThr) && (chose < aOne + aTwo + 2.*aThr) ) 01941 { 01942 xRand = rRand*cosphi; 01943 yRand = rRand*sinphi; 01944 zRand = GetCutZ(G4ThreeVector(xRand,yRand,-fDz)); 01945 return G4ThreeVector (xRand, yRand, zRand); 01946 } 01947 else if( (chose >= aOne + aTwo + 2.*aThr) 01948 && (chose < aOne + aTwo + 2.*aThr + aFou) ) 01949 { 01950 xRand = rRand*std::cos(fSPhi); 01951 yRand = rRand*std::sin(fSPhi); 01952 zRand = RandFlat::shoot(GetCutZ(G4ThreeVector(xRand,yRand,-fDz)), 01953 GetCutZ(G4ThreeVector(xRand,yRand,fDz))); 01954 return G4ThreeVector (xRand, yRand, zRand); 01955 } 01956 else 01957 { 01958 xRand = rRand*std::cos(fSPhi+fDPhi); 01959 yRand = rRand*std::sin(fSPhi+fDPhi); 01960 zRand = RandFlat::shoot(GetCutZ(G4ThreeVector(xRand,yRand,-fDz)), 01961 GetCutZ(G4ThreeVector(xRand,yRand,fDz))); 01962 return G4ThreeVector (xRand, yRand, zRand); 01963 } 01964 }
G4double G4CutTubs::GetSurfaceArea | ( | ) | [inline, virtual] |
Reimplemented from G4Tubs.
Definition at line 43 of file G4CutTubs.icc.
References G4VSolid::GetSurfaceArea().
00044 { 00045 return G4VSolid::GetSurfaceArea(); 00046 }
EInside G4CutTubs::Inside | ( | const G4ThreeVector & | p | ) | const [virtual] |
Reimplemented from G4Tubs.
Definition at line 407 of file G4CutTubs.cc.
References G4Tubs::fDPhi, G4Tubs::fDz, G4Tubs::fRMax, G4Tubs::fRMin, G4Tubs::fSPhi, G4Tubs::kAngTolerance, G4VSolid::kCarTolerance, kInside, kOutside, G4Tubs::kRadTolerance, and kSurface.
Referenced by CalculateExtent().
00408 { 00409 G4double zinLow,zinHigh,r2,pPhi=0.; 00410 G4double tolRMin,tolRMax; 00411 G4ThreeVector vZ=G4ThreeVector(0,0,fDz); 00412 EInside in = kInside; 00413 static const G4double halfCarTolerance=kCarTolerance*0.5; 00414 static const G4double halfRadTolerance=kRadTolerance*0.5; 00415 static const G4double halfAngTolerance=kAngTolerance*0.5; 00416 00417 // Check if point is contained in the cut plane in -/+ Z 00418 00419 // Check the lower cut plane 00420 // 00421 zinLow =(p+vZ).dot(fLowNorm); 00422 if (zinLow>halfCarTolerance) { return kOutside; } 00423 00424 // Check the higher cut plane 00425 // 00426 zinHigh = (p-vZ).dot(fHighNorm); 00427 if (zinHigh>halfCarTolerance) { return kOutside; } 00428 00429 // Check radius 00430 // 00431 r2 = p.x()*p.x() + p.y()*p.y() ; 00432 00433 // First check 'generous' boundaries R+tolerance 00434 // 00435 tolRMin = fRMin - halfRadTolerance ; 00436 tolRMax = fRMax + halfRadTolerance ; 00437 if ( tolRMin < 0 ) { tolRMin = 0; } 00438 00439 if ( ((r2 < tolRMin*tolRMin) || (r2 > tolRMax*tolRMax)) 00440 && (r2 >=halfRadTolerance*halfRadTolerance) ) { return kOutside; } 00441 00442 // Check Phi 00443 // 00444 if(!fPhiFullCutTube) 00445 { 00446 // Try outer tolerant phi boundaries only 00447 00448 if ( (tolRMin==0) && (std::fabs(p.x())<=halfCarTolerance) 00449 && (std::fabs(p.y())<=halfCarTolerance) ) 00450 { 00451 return kSurface; 00452 } 00453 00454 pPhi = std::atan2(p.y(),p.x()) ; 00455 00456 if ( pPhi < -halfAngTolerance) { pPhi += twopi; } // 0<=pPhi<2pi 00457 if ( fSPhi >= 0 ) 00458 { 00459 if ( (std::fabs(pPhi) < halfAngTolerance) 00460 && (std::fabs(fSPhi + fDPhi - twopi) < halfAngTolerance) ) 00461 { 00462 pPhi += twopi ; // 0 <= pPhi < 2pi 00463 } 00464 if ( (pPhi <= fSPhi - halfAngTolerance) 00465 || (pPhi >= fSPhi + fDPhi + halfAngTolerance) ) 00466 { 00467 in = kOutside ; 00468 } 00469 else if ( (pPhi <= fSPhi + halfAngTolerance) 00470 || (pPhi >= fSPhi + fDPhi - halfAngTolerance) ) 00471 { 00472 in=kSurface; 00473 } 00474 } 00475 else // fSPhi < 0 00476 { 00477 if ( (pPhi <= fSPhi + twopi - halfAngTolerance) 00478 && (pPhi >= fSPhi + fDPhi + halfAngTolerance) ) 00479 { 00480 in = kOutside; 00481 } 00482 else 00483 { 00484 in = kSurface ; 00485 } 00486 } 00487 } 00488 00489 // Check on the Surface for Z 00490 // 00491 if ((zinLow>=-halfCarTolerance) 00492 || (zinHigh>=-halfCarTolerance)) 00493 { 00494 in=kSurface; 00495 } 00496 00497 // Check on the Surface for R 00498 // 00499 if (fRMin) { tolRMin = fRMin + halfRadTolerance ; } 00500 else { tolRMin = 0 ; } 00501 tolRMax = fRMax - halfRadTolerance ; 00502 if ( ((r2 <= tolRMin*tolRMin) || (r2 >= tolRMax*tolRMax))&& 00503 (r2 >=halfRadTolerance*halfRadTolerance) ) 00504 { 00505 return kSurface; 00506 } 00507 00508 return in; 00509 }
G4bool G4CutTubs::IsCrossingCutPlanes | ( | ) | const [protected] |
Definition at line 2035 of file G4CutTubs.cc.
References G4Tubs::fDz, G4Tubs::fRMax, and GetCutZ().
Referenced by G4CutTubs().
02036 { 02037 G4double zXLow1,zXLow2,zYLow1,zYLow2; 02038 G4double zXHigh1,zXHigh2,zYHigh1,zYHigh2; 02039 02040 zXLow1 = GetCutZ(G4ThreeVector(-fRMax, 0,-fDz)); 02041 zXLow2 = GetCutZ(G4ThreeVector( fRMax, 0,-fDz)); 02042 zYLow1 = GetCutZ(G4ThreeVector( 0,-fRMax,-fDz)); 02043 zYLow2 = GetCutZ(G4ThreeVector( 0, fRMax,-fDz)); 02044 zXHigh1 = GetCutZ(G4ThreeVector(-fRMax, 0, fDz)); 02045 zXHigh2 = GetCutZ(G4ThreeVector( fRMax, 0, fDz)); 02046 zYHigh1 = GetCutZ(G4ThreeVector( 0,-fRMax, fDz)); 02047 zYHigh2 = GetCutZ(G4ThreeVector( 0, fRMax, fDz)); 02048 if ( (zXLow1>zXHigh1) ||(zXLow2>zXHigh2) 02049 || (zYLow1>zYHigh1) ||(zYLow2>zYHigh2)) { return true; } 02050 02051 return false; 02052 }
Definition at line 172 of file G4CutTubs.cc.
References fHighNorm, fLowNorm, fPhiFullCutTube, and G4Tubs::operator=().
00173 { 00174 // Check assignment to self 00175 // 00176 if (this == &rhs) { return *this; } 00177 00178 // Copy base class data 00179 // 00180 G4Tubs::operator=(rhs); 00181 00182 // Copy data 00183 // 00184 fLowNorm = rhs.fLowNorm; fHighNorm = rhs.fHighNorm; 00185 fPhiFullCutTube = rhs.fPhiFullCutTube; 00186 00187 return *this; 00188 }
std::ostream & G4CutTubs::StreamInfo | ( | std::ostream & | os | ) | const [virtual] |
Reimplemented from G4Tubs.
Definition at line 1871 of file G4CutTubs.cc.
References G4Tubs::fDPhi, G4Tubs::fDz, G4Tubs::fRMax, G4Tubs::fRMin, G4Tubs::fSPhi, and G4VSolid::GetName().
01872 { 01873 G4int oldprc = os.precision(16); 01874 os << "-----------------------------------------------------------\n" 01875 << " *** Dump for solid - " << GetName() << " ***\n" 01876 << " ===================================================\n" 01877 << " Solid type: G4CutTubs\n" 01878 << " Parameters: \n" 01879 << " inner radius : " << fRMin/mm << " mm \n" 01880 << " outer radius : " << fRMax/mm << " mm \n" 01881 << " half length Z: " << fDz/mm << " mm \n" 01882 << " starting phi : " << fSPhi/degree << " degrees \n" 01883 << " delta phi : " << fDPhi/degree << " degrees \n" 01884 << " low Norm : " << fLowNorm << " \n" 01885 << " high Norm : " <<fHighNorm << " \n" 01886 << "-----------------------------------------------------------\n"; 01887 os.precision(oldprc); 01888 01889 return os; 01890 }
G4ThreeVector G4CutTubs::SurfaceNormal | ( | const G4ThreeVector & | p | ) | const [virtual] |
Reimplemented from G4Tubs.
Definition at line 517 of file G4CutTubs.cc.
References ApproxSurfaceNormal(), G4Tubs::fDPhi, G4Tubs::fDz, G4Tubs::fRMax, G4Tubs::fRMin, G4Tubs::fSPhi, G4cout, G4endl, G4Exception(), JustWarning, G4Tubs::kAngTolerance, and G4VSolid::kCarTolerance.
00518 { 00519 G4int noSurfaces = 0; 00520 G4double rho, pPhi; 00521 G4double distZLow,distZHigh, distRMin, distRMax; 00522 G4double distSPhi = kInfinity, distEPhi = kInfinity; 00523 G4ThreeVector vZ=G4ThreeVector(0,0,fDz); 00524 00525 static const G4double halfCarTolerance = 0.5*kCarTolerance; 00526 static const G4double halfAngTolerance = 0.5*kAngTolerance; 00527 00528 G4ThreeVector norm, sumnorm(0.,0.,0.); 00529 G4ThreeVector nZ = G4ThreeVector(0, 0, 1.0); 00530 G4ThreeVector nR, nPs, nPe; 00531 00532 rho = std::sqrt(p.x()*p.x() + p.y()*p.y()); 00533 00534 distRMin = std::fabs(rho - fRMin); 00535 distRMax = std::fabs(rho - fRMax); 00536 00537 // dist to Low Cut 00538 // 00539 distZLow =std::fabs((p+vZ).dot(fLowNorm)); 00540 00541 // dist to High Cut 00542 // 00543 distZHigh = std::fabs((p-vZ).dot(fHighNorm)); 00544 00545 if (!fPhiFullCutTube) // Protected against (0,0,z) 00546 { 00547 if ( rho > halfCarTolerance ) 00548 { 00549 pPhi = std::atan2(p.y(),p.x()); 00550 00551 if(pPhi < fSPhi- halfCarTolerance) { pPhi += twopi; } 00552 else if(pPhi > fSPhi+fDPhi+ halfCarTolerance) { pPhi -= twopi; } 00553 00554 distSPhi = std::fabs(pPhi - fSPhi); 00555 distEPhi = std::fabs(pPhi - fSPhi - fDPhi); 00556 } 00557 else if( !fRMin ) 00558 { 00559 distSPhi = 0.; 00560 distEPhi = 0.; 00561 } 00562 nPs = G4ThreeVector(std::sin(fSPhi),-std::cos(fSPhi),0); 00563 nPe = G4ThreeVector(-std::sin(fSPhi+fDPhi),std::cos(fSPhi+fDPhi),0); 00564 } 00565 if ( rho > halfCarTolerance ) { nR = G4ThreeVector(p.x()/rho,p.y()/rho,0); } 00566 00567 if( distRMax <= halfCarTolerance ) 00568 { 00569 noSurfaces ++; 00570 sumnorm += nR; 00571 } 00572 if( fRMin && (distRMin <= halfCarTolerance) ) 00573 { 00574 noSurfaces ++; 00575 sumnorm -= nR; 00576 } 00577 if( fDPhi < twopi ) 00578 { 00579 if (distSPhi <= halfAngTolerance) 00580 { 00581 noSurfaces ++; 00582 sumnorm += nPs; 00583 } 00584 if (distEPhi <= halfAngTolerance) 00585 { 00586 noSurfaces ++; 00587 sumnorm += nPe; 00588 } 00589 } 00590 if (distZLow <= halfCarTolerance) 00591 { 00592 noSurfaces ++; 00593 sumnorm += fLowNorm; 00594 } 00595 if (distZHigh <= halfCarTolerance) 00596 { 00597 noSurfaces ++; 00598 sumnorm += fHighNorm; 00599 } 00600 if ( noSurfaces == 0 ) 00601 { 00602 #ifdef G4CSGDEBUG 00603 G4Exception("G4CutTubs::SurfaceNormal(p)", "GeomSolids1002", 00604 JustWarning, "Point p is not on surface !?" ); 00605 G4int oldprc = G4cout.precision(20); 00606 G4cout<< "G4CutTubs::SN ( "<<p.x()<<", "<<p.y()<<", "<<p.z()<<" ); " 00607 << G4endl << G4endl; 00608 G4cout.precision(oldprc) ; 00609 #endif 00610 norm = ApproxSurfaceNormal(p); 00611 } 00612 else if ( noSurfaces == 1 ) { norm = sumnorm; } 00613 else { norm = sumnorm.unit(); } 00614 00615 return norm; 00616 }