#include <G4Ellipsoid.hh>
Inheritance diagram for G4Ellipsoid:
Definition at line 59 of file G4Ellipsoid.hh.
G4Ellipsoid::G4Ellipsoid | ( | const G4String & | pName, | |
G4double | pxSemiAxis, | |||
G4double | pySemiAxis, | |||
G4double | pzSemiAxis, | |||
G4double | pzBottomCut = 0 , |
|||
G4double | pzTopCut = 0 | |||
) |
Definition at line 64 of file G4Ellipsoid.cc.
References FatalErrorInArgument, G4Exception(), G4GeometryTolerance::GetInstance(), G4VSolid::GetName(), G4GeometryTolerance::GetRadialTolerance(), SetSemiAxis(), and SetZCuts().
Referenced by Clone().
00070 : G4VSolid(pName), fpPolyhedron(0), fCubicVolume(0.), fSurfaceArea(0.), 00071 zBottomCut(0.), zTopCut(0.) 00072 { 00073 // note: for users that want to use the full ellipsoid it is useful 00074 // to include a default for the cuts 00075 00076 kRadTolerance = G4GeometryTolerance::GetInstance()->GetRadialTolerance(); 00077 00078 // Check Semi-Axis 00079 if ( (pxSemiAxis<=0.) || (pySemiAxis<=0.) || (pzSemiAxis<=0.) ) 00080 { 00081 std::ostringstream message; 00082 message << "Invalid semi-axis - " << GetName(); 00083 G4Exception("G4Ellipsoid::G4Ellipsoid()", "GeomSolids0002", 00084 FatalErrorInArgument, message); 00085 } 00086 SetSemiAxis(pxSemiAxis, pySemiAxis, pzSemiAxis); 00087 00088 if ( pzBottomCut == 0 && pzTopCut == 0 ) 00089 { 00090 SetZCuts(-pzSemiAxis, pzSemiAxis); 00091 } 00092 else if ( (pzBottomCut < pzSemiAxis) && (pzTopCut > -pzSemiAxis) 00093 && (pzBottomCut < pzTopCut) ) 00094 { 00095 SetZCuts(pzBottomCut, pzTopCut); 00096 } 00097 else 00098 { 00099 std::ostringstream message; 00100 message << "Invalid z-coordinate for cutting plane - " << GetName(); 00101 G4Exception("G4Ellipsoid::G4Ellipsoid()", "GeomSolids0002", 00102 FatalErrorInArgument, message); 00103 } 00104 }
G4Ellipsoid::~G4Ellipsoid | ( | ) | [virtual] |
G4Ellipsoid::G4Ellipsoid | ( | __void__ & | ) |
Definition at line 111 of file G4Ellipsoid.cc.
00112 : G4VSolid(a), fpPolyhedron(0), kRadTolerance(0.), fCubicVolume(0.), 00113 fSurfaceArea(0.), xSemiAxis(0.), ySemiAxis(0.), zSemiAxis(0.), 00114 semiAxisMax(0.), zBottomCut(0.), zTopCut(0.) 00115 { 00116 }
G4Ellipsoid::G4Ellipsoid | ( | const G4Ellipsoid & | rhs | ) |
Definition at line 130 of file G4Ellipsoid.cc.
00131 : G4VSolid(rhs), 00132 fpPolyhedron(0), kRadTolerance(rhs.kRadTolerance), 00133 fCubicVolume(rhs.fCubicVolume), fSurfaceArea(rhs.fSurfaceArea), 00134 xSemiAxis(rhs.xSemiAxis), ySemiAxis(rhs.ySemiAxis), 00135 zSemiAxis(rhs.zSemiAxis), semiAxisMax(rhs.semiAxisMax), 00136 zBottomCut(rhs.zBottomCut), zTopCut(rhs.zTopCut) 00137 { 00138 }
G4bool G4Ellipsoid::CalculateExtent | ( | const EAxis | pAxis, | |
const G4VoxelLimits & | pVoxelLimit, | |||
const G4AffineTransform & | pTransform, | |||
G4double & | pmin, | |||
G4double & | pmax | |||
) | const [virtual] |
Implements G4VSolid.
Definition at line 170 of file G4Ellipsoid.cc.
References G4VSolid::CalculateClippedPolygonExtent(), CreateRotatedVertices(), 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(), sqr(), and G4AffineTransform::TransformPoint().
00174 { 00175 if (!pTransform.IsRotated()) 00176 { 00177 // Special case handling for unrotated solid ellipsoid 00178 // Compute x/y/z mins and maxs for bounding box respecting limits, 00179 // with early returns if outside limits. Then switch() on pAxis, 00180 // and compute exact x and y limit for x/y case 00181 00182 G4double xoffset,xMin,xMax; 00183 G4double yoffset,yMin,yMax; 00184 G4double zoffset,zMin,zMax; 00185 00186 G4double maxDiff,newMin,newMax; 00187 G4double xoff,yoff; 00188 00189 xoffset=pTransform.NetTranslation().x(); 00190 xMin=xoffset - xSemiAxis; 00191 xMax=xoffset + xSemiAxis; 00192 if (pVoxelLimit.IsXLimited()) 00193 { 00194 if ( (xMin>pVoxelLimit.GetMaxXExtent()+kCarTolerance) 00195 || (xMax<pVoxelLimit.GetMinXExtent()-kCarTolerance) ) 00196 { 00197 return false; 00198 } 00199 else 00200 { 00201 if (xMin<pVoxelLimit.GetMinXExtent()) 00202 { 00203 xMin=pVoxelLimit.GetMinXExtent(); 00204 } 00205 if (xMax>pVoxelLimit.GetMaxXExtent()) 00206 { 00207 xMax=pVoxelLimit.GetMaxXExtent(); 00208 } 00209 } 00210 } 00211 00212 yoffset=pTransform.NetTranslation().y(); 00213 yMin=yoffset - ySemiAxis; 00214 yMax=yoffset + ySemiAxis; 00215 if (pVoxelLimit.IsYLimited()) 00216 { 00217 if ( (yMin>pVoxelLimit.GetMaxYExtent()+kCarTolerance) 00218 || (yMax<pVoxelLimit.GetMinYExtent()-kCarTolerance) ) 00219 { 00220 return false; 00221 } 00222 else 00223 { 00224 if (yMin<pVoxelLimit.GetMinYExtent()) 00225 { 00226 yMin=pVoxelLimit.GetMinYExtent(); 00227 } 00228 if (yMax>pVoxelLimit.GetMaxYExtent()) 00229 { 00230 yMax=pVoxelLimit.GetMaxYExtent(); 00231 } 00232 } 00233 } 00234 00235 zoffset=pTransform.NetTranslation().z(); 00236 zMin=zoffset + (-zSemiAxis > zBottomCut ? -zSemiAxis : zBottomCut); 00237 zMax=zoffset + ( zSemiAxis < zTopCut ? zSemiAxis : zTopCut); 00238 if (pVoxelLimit.IsZLimited()) 00239 { 00240 if ( (zMin>pVoxelLimit.GetMaxZExtent()+kCarTolerance) 00241 || (zMax<pVoxelLimit.GetMinZExtent()-kCarTolerance) ) 00242 { 00243 return false; 00244 } 00245 else 00246 { 00247 if (zMin<pVoxelLimit.GetMinZExtent()) 00248 { 00249 zMin=pVoxelLimit.GetMinZExtent(); 00250 } 00251 if (zMax>pVoxelLimit.GetMaxZExtent()) 00252 { 00253 zMax=pVoxelLimit.GetMaxZExtent(); 00254 } 00255 } 00256 } 00257 00258 // if here, then known to cut bounding box around ellipsoid 00259 // 00260 xoff = (xoffset < xMin) ? (xMin-xoffset) 00261 : (xoffset > xMax) ? (xoffset-xMax) : 0.0; 00262 yoff = (yoffset < yMin) ? (yMin-yoffset) 00263 : (yoffset > yMax) ? (yoffset-yMax) : 0.0; 00264 00265 // detailed calculations 00266 // NOTE: does not use X or Y offsets to adjust Z range, 00267 // and does not use Z offset to adjust X or Y range, 00268 // which is consistent with G4Sphere::CalculateExtent behavior 00269 // 00270 switch (pAxis) 00271 { 00272 case kXAxis: 00273 if (yoff==0.) 00274 { 00275 // YZ limits cross max/min x => no change 00276 // 00277 pMin=xMin; 00278 pMax=xMax; 00279 } 00280 else 00281 { 00282 // YZ limits don't cross max/min x => compute max delta x, 00283 // hence new mins/maxs 00284 // 00285 maxDiff= 1.0-sqr(yoff/ySemiAxis); 00286 if (maxDiff < 0.0) { return false; } 00287 maxDiff= xSemiAxis * std::sqrt(maxDiff); 00288 newMin=xoffset-maxDiff; 00289 newMax=xoffset+maxDiff; 00290 pMin=(newMin<xMin) ? xMin : newMin; 00291 pMax=(newMax>xMax) ? xMax : newMax; 00292 } 00293 break; 00294 case kYAxis: 00295 if (xoff==0.) 00296 { 00297 // XZ limits cross max/min y => no change 00298 // 00299 pMin=yMin; 00300 pMax=yMax; 00301 } 00302 else 00303 { 00304 // XZ limits don't cross max/min y => compute max delta y, 00305 // hence new mins/maxs 00306 // 00307 maxDiff= 1.0-sqr(xoff/xSemiAxis); 00308 if (maxDiff < 0.0) { return false; } 00309 maxDiff= ySemiAxis * std::sqrt(maxDiff); 00310 newMin=yoffset-maxDiff; 00311 newMax=yoffset+maxDiff; 00312 pMin=(newMin<yMin) ? yMin : newMin; 00313 pMax=(newMax>yMax) ? yMax : newMax; 00314 } 00315 break; 00316 case kZAxis: 00317 pMin=zMin; 00318 pMax=zMax; 00319 break; 00320 default: 00321 break; 00322 } 00323 00324 pMin-=kCarTolerance; 00325 pMax+=kCarTolerance; 00326 return true; 00327 } 00328 else // not rotated 00329 { 00330 G4int i,j,noEntries,noBetweenSections; 00331 G4bool existsAfterClip=false; 00332 00333 // Calculate rotated vertex coordinates 00334 00335 G4int noPolygonVertices=0; 00336 G4ThreeVectorList* vertices = 00337 CreateRotatedVertices(pTransform,noPolygonVertices); 00338 00339 pMin=+kInfinity; 00340 pMax=-kInfinity; 00341 00342 noEntries=vertices->size(); // noPolygonVertices*noPhiCrossSections 00343 noBetweenSections=noEntries-noPolygonVertices; 00344 00345 G4ThreeVectorList ThetaPolygon; 00346 for (i=0;i<noEntries;i+=noPolygonVertices) 00347 { 00348 for(j=0;j<(noPolygonVertices/2)-1;j++) 00349 { 00350 ThetaPolygon.push_back((*vertices)[i+j]); 00351 ThetaPolygon.push_back((*vertices)[i+j+1]); 00352 ThetaPolygon.push_back((*vertices)[i+noPolygonVertices-2-j]); 00353 ThetaPolygon.push_back((*vertices)[i+noPolygonVertices-1-j]); 00354 CalculateClippedPolygonExtent(ThetaPolygon,pVoxelLimit,pAxis,pMin,pMax); 00355 ThetaPolygon.clear(); 00356 } 00357 } 00358 for (i=0;i<noBetweenSections;i+=noPolygonVertices) 00359 { 00360 for(j=0;j<noPolygonVertices-1;j++) 00361 { 00362 ThetaPolygon.push_back((*vertices)[i+j]); 00363 ThetaPolygon.push_back((*vertices)[i+j+1]); 00364 ThetaPolygon.push_back((*vertices)[i+noPolygonVertices+j+1]); 00365 ThetaPolygon.push_back((*vertices)[i+noPolygonVertices+j]); 00366 CalculateClippedPolygonExtent(ThetaPolygon,pVoxelLimit,pAxis,pMin,pMax); 00367 ThetaPolygon.clear(); 00368 } 00369 ThetaPolygon.push_back((*vertices)[i+noPolygonVertices-1]); 00370 ThetaPolygon.push_back((*vertices)[i]); 00371 ThetaPolygon.push_back((*vertices)[i+noPolygonVertices]); 00372 ThetaPolygon.push_back((*vertices)[i+2*noPolygonVertices-1]); 00373 CalculateClippedPolygonExtent(ThetaPolygon,pVoxelLimit,pAxis,pMin,pMax); 00374 ThetaPolygon.clear(); 00375 } 00376 if ( (pMin!=kInfinity) || (pMax!=-kInfinity) ) 00377 { 00378 existsAfterClip=true; 00379 00380 // Add 2*tolerance to avoid precision troubles 00381 // 00382 pMin-=kCarTolerance; 00383 pMax+=kCarTolerance; 00384 00385 } 00386 else 00387 { 00388 // Check for case where completely enveloping clipping volume 00389 // If point inside then we are confident that the solid completely 00390 // envelopes the clipping volume. Hence set min/max extents according 00391 // to clipping volume extents along the specified axis. 00392 // 00393 G4ThreeVector 00394 clipCentre((pVoxelLimit.GetMinXExtent()+pVoxelLimit.GetMaxXExtent())*0.5, 00395 (pVoxelLimit.GetMinYExtent()+pVoxelLimit.GetMaxYExtent())*0.5, 00396 (pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5); 00397 00398 if (Inside(pTransform.Inverse().TransformPoint(clipCentre))!=kOutside) 00399 { 00400 existsAfterClip=true; 00401 pMin=pVoxelLimit.GetMinExtent(pAxis); 00402 pMax=pVoxelLimit.GetMaxExtent(pAxis); 00403 } 00404 } 00405 delete vertices; 00406 return existsAfterClip; 00407 } 00408 }
G4VSolid * G4Ellipsoid::Clone | ( | ) | const [virtual] |
Reimplemented from G4VSolid.
Definition at line 947 of file G4Ellipsoid.cc.
References G4Ellipsoid().
00948 { 00949 return new G4Ellipsoid(*this); 00950 }
G4NURBS * G4Ellipsoid::CreateNURBS | ( | ) | const [virtual] |
Reimplemented from G4VSolid.
Definition at line 1063 of file G4Ellipsoid.cc.
01064 { 01065 // Box for now!!! 01066 // 01067 return new G4NURBSbox(semiAxisMax, semiAxisMax, semiAxisMax); 01068 }
G4Polyhedron * G4Ellipsoid::CreatePolyhedron | ( | ) | const [virtual] |
Reimplemented from G4VSolid.
Definition at line 1070 of file G4Ellipsoid.cc.
Referenced by GetPolyhedron().
01071 { 01072 return new G4PolyhedronEllipsoid(xSemiAxis, ySemiAxis, zSemiAxis, 01073 zBottomCut, zTopCut); 01074 }
G4ThreeVectorList * G4Ellipsoid::CreateRotatedVertices | ( | const G4AffineTransform & | pT, | |
G4int & | noPV | |||
) | const [protected] |
Definition at line 822 of file G4Ellipsoid.cc.
References G4VSolid::DumpInfo(), FatalException, G4Exception(), kMeshAngleDefault, G4INCL::Math::pi, and G4AffineTransform::TransformPoint().
Referenced by CalculateExtent().
00824 { 00825 G4ThreeVectorList *vertices; 00826 G4ThreeVector vertex; 00827 G4double meshAnglePhi, meshRMaxFactor, 00828 crossAnglePhi, coscrossAnglePhi, sincrossAnglePhi, sAnglePhi; 00829 G4double meshTheta, crossTheta, startTheta; 00830 G4double rMaxX, rMaxY, rMaxZ, rMaxMax, rx, ry, rz; 00831 G4int crossSectionPhi, noPhiCrossSections, crossSectionTheta, noThetaSections; 00832 00833 // Phi cross sections 00834 // 00835 noPhiCrossSections=G4int (twopi/kMeshAngleDefault)+1; // = 9! 00836 00837 /* 00838 if (noPhiCrossSections<kMinMeshSections) // <3 00839 { 00840 noPhiCrossSections=kMinMeshSections; 00841 } 00842 else if (noPhiCrossSections>kMaxMeshSections) // >37 00843 { 00844 noPhiCrossSections=kMaxMeshSections; 00845 } 00846 */ 00847 meshAnglePhi=twopi/(noPhiCrossSections-1); 00848 00849 // Set start angle such that mesh will be at fRMax 00850 // on the x axis. Will give better extent calculations when not rotated. 00851 00852 sAnglePhi = -meshAnglePhi*0.5; 00853 00854 // Theta cross sections 00855 00856 noThetaSections = G4int(pi/kMeshAngleDefault)+3; // = 7! 00857 00858 /* 00859 if (noThetaSections<kMinMeshSections) // <3 00860 { 00861 noThetaSections=kMinMeshSections; 00862 } 00863 else if (noThetaSections>kMaxMeshSections) // >37 00864 { 00865 noThetaSections=kMaxMeshSections; 00866 } 00867 */ 00868 meshTheta= pi/(noThetaSections-2); 00869 00870 // Set start angle such that mesh will be at fRMax 00871 // on the z axis. Will give better extent calculations when not rotated. 00872 00873 startTheta = -meshTheta*0.5; 00874 00875 meshRMaxFactor = 1.0/std::cos(0.5* 00876 std::sqrt(meshAnglePhi*meshAnglePhi+meshTheta*meshTheta)); 00877 rMaxMax= (xSemiAxis > ySemiAxis ? xSemiAxis : ySemiAxis); 00878 if (zSemiAxis > rMaxMax) rMaxMax= zSemiAxis; 00879 rMaxX= xSemiAxis + rMaxMax*(meshRMaxFactor-1.0); 00880 rMaxY= ySemiAxis + rMaxMax*(meshRMaxFactor-1.0); 00881 rMaxZ= zSemiAxis + rMaxMax*(meshRMaxFactor-1.0); 00882 G4double* cosCrossTheta = new G4double[noThetaSections]; 00883 G4double* sinCrossTheta = new G4double[noThetaSections]; 00884 vertices=new G4ThreeVectorList(noPhiCrossSections*noThetaSections); 00885 if (vertices && cosCrossTheta && sinCrossTheta) 00886 { 00887 for (crossSectionTheta=0; crossSectionTheta<noThetaSections; 00888 crossSectionTheta++) 00889 { 00890 // Compute sine and cosine table (for historical reasons) 00891 // 00892 crossTheta=startTheta+crossSectionTheta*meshTheta; 00893 cosCrossTheta[crossSectionTheta]=std::cos(crossTheta); 00894 sinCrossTheta[crossSectionTheta]=std::sin(crossTheta); 00895 } 00896 for (crossSectionPhi=0; crossSectionPhi<noPhiCrossSections; 00897 crossSectionPhi++) 00898 { 00899 crossAnglePhi=sAnglePhi+crossSectionPhi*meshAnglePhi; 00900 coscrossAnglePhi=std::cos(crossAnglePhi); 00901 sincrossAnglePhi=std::sin(crossAnglePhi); 00902 for (crossSectionTheta=0; crossSectionTheta<noThetaSections; 00903 crossSectionTheta++) 00904 { 00905 // Compute coordinates of cross section at section crossSectionPhi 00906 // 00907 rx= sinCrossTheta[crossSectionTheta]*coscrossAnglePhi*rMaxX; 00908 ry= sinCrossTheta[crossSectionTheta]*sincrossAnglePhi*rMaxY; 00909 rz= cosCrossTheta[crossSectionTheta]*rMaxZ; 00910 if (rz < zBottomCut) 00911 { rz= zBottomCut; } 00912 if (rz > zTopCut) 00913 { rz= zTopCut; } 00914 vertex= G4ThreeVector(rx,ry,rz); 00915 vertices->push_back(pTransform.TransformPoint(vertex)); 00916 } // Theta forward 00917 } // Phi 00918 noPolygonVertices = noThetaSections ; 00919 } 00920 else 00921 { 00922 DumpInfo(); 00923 G4Exception("G4Ellipsoid::CreateRotatedVertices()", 00924 "GeomSolids0003", FatalException, 00925 "Error in allocation of vertices. Out of memory !"); 00926 } 00927 00928 delete[] cosCrossTheta; 00929 delete[] sinCrossTheta; 00930 00931 return vertices; 00932 }
void G4Ellipsoid::DescribeYourselfTo | ( | G4VGraphicsScene & | scene | ) | const [virtual] |
Implements G4VSolid.
Definition at line 1049 of file G4Ellipsoid.cc.
References G4VGraphicsScene::AddSolid().
01050 { 01051 scene.AddSolid(*this); 01052 }
G4double G4Ellipsoid::DistanceToIn | ( | const G4ThreeVector & | p | ) | const [virtual] |
Implements G4VSolid.
Definition at line 590 of file G4Ellipsoid.cc.
00591 { 00592 G4double distR, distZ; 00593 00594 // normal vector: parallel to normal, magnitude 1/(characteristic radius) 00595 // 00596 G4ThreeVector norm(p.x()/(xSemiAxis*xSemiAxis), 00597 p.y()/(ySemiAxis*ySemiAxis), 00598 p.z()/(zSemiAxis*zSemiAxis)); 00599 G4double radius= 1.0/norm.mag(); 00600 00601 // approximate distance to curved surface ( <= actual distance ) 00602 // 00603 distR= (p*norm - 1.0) * radius / 2.0; 00604 00605 // Distance to z-cut plane 00606 // 00607 distZ= zBottomCut - p.z(); 00608 if (distZ < 0.0) 00609 { 00610 distZ = p.z() - zTopCut; 00611 } 00612 00613 // Distance to closest surface from outside 00614 // 00615 if (distZ < 0.0) 00616 { 00617 return (distR < 0.0) ? 0.0 : distR; 00618 } 00619 else if (distR < 0.0) 00620 { 00621 return distZ; 00622 } 00623 else 00624 { 00625 return (distZ < distR) ? distZ : distR; 00626 } 00627 }
G4double G4Ellipsoid::DistanceToIn | ( | const G4ThreeVector & | p, | |
const G4ThreeVector & | v | |||
) | const [virtual] |
Implements G4VSolid.
Definition at line 494 of file G4Ellipsoid.cc.
References Inside(), G4VSolid::kCarTolerance, kOutside, sqr(), and SurfaceNormal().
00496 { 00497 static const G4double halfCarTolerance=kCarTolerance*0.5; 00498 static const G4double halfRadTolerance=kRadTolerance*0.5; 00499 00500 G4double distMin = std::min(xSemiAxis,ySemiAxis); 00501 const G4double dRmax = 100.*std::min(distMin,zSemiAxis); 00502 distMin= kInfinity; 00503 00504 // check to see if Z plane is relevant 00505 if (p.z() <= zBottomCut+halfCarTolerance) 00506 { 00507 if (v.z() <= 0.0) { return distMin; } 00508 G4double distZ = (zBottomCut - p.z()) / v.z(); 00509 00510 if ( (distZ > -halfRadTolerance) && (Inside(p+distZ*v) != kOutside) ) 00511 { 00512 // early exit since can't intercept curved surface if we reach here 00513 if ( std::fabs(distZ) < halfRadTolerance ) { distZ=0.; } 00514 return distMin= distZ; 00515 } 00516 } 00517 if (p.z() >= zTopCut-halfCarTolerance) 00518 { 00519 if (v.z() >= 0.0) { return distMin;} 00520 G4double distZ = (zTopCut - p.z()) / v.z(); 00521 if ( (distZ > -halfRadTolerance) && (Inside(p+distZ*v) != kOutside) ) 00522 { 00523 // early exit since can't intercept curved surface if we reach here 00524 if ( std::fabs(distZ) < halfRadTolerance ) { distZ=0.; } 00525 return distMin= distZ; 00526 } 00527 } 00528 // if fZCut1 <= p.z() <= fZCut2, then must hit curved surface 00529 00530 // now check curved surface intercept 00531 G4double A,B,C; 00532 00533 A= sqr(v.x()/xSemiAxis) + sqr(v.y()/ySemiAxis) + sqr(v.z()/zSemiAxis); 00534 C= sqr(p.x()/xSemiAxis) + sqr(p.y()/ySemiAxis) + sqr(p.z()/zSemiAxis) - 1.0; 00535 B= 2.0 * ( p.x()*v.x()/(xSemiAxis*xSemiAxis) 00536 + p.y()*v.y()/(ySemiAxis*ySemiAxis) 00537 + p.z()*v.z()/(zSemiAxis*zSemiAxis) ); 00538 00539 C= B*B - 4.0*A*C; 00540 if (C > 0.0) 00541 { 00542 G4double distR= (-B - std::sqrt(C)) / (2.0*A); 00543 G4double intZ = p.z()+distR*v.z(); 00544 if ( (distR > halfRadTolerance) 00545 && (intZ >= zBottomCut-halfRadTolerance) 00546 && (intZ <= zTopCut+halfRadTolerance) ) 00547 { 00548 distMin = distR; 00549 } 00550 else if( (distR >- halfRadTolerance) 00551 && (intZ >= zBottomCut-halfRadTolerance) 00552 && (intZ <= zTopCut+halfRadTolerance) ) 00553 { 00554 // p is on the curved surface, DistanceToIn returns 0 or kInfinity: 00555 // DistanceToIn returns 0, if second root is positive (means going inside) 00556 // If second root is negative, DistanceToIn returns kInfinity (outside) 00557 // 00558 distR = (-B + std::sqrt(C) ) / (2.0*A); 00559 if(distR>0.) { distMin=0.; } 00560 } 00561 else 00562 { 00563 distR= (-B + std::sqrt(C)) / (2.0*A); 00564 intZ = p.z()+distR*v.z(); 00565 if ( (distR > halfRadTolerance) 00566 && (intZ >= zBottomCut-halfRadTolerance) 00567 && (intZ <= zTopCut+halfRadTolerance) ) 00568 { 00569 G4ThreeVector norm=SurfaceNormal(p); 00570 if (norm.dot(v)<0.) { distMin = distR; } 00571 } 00572 } 00573 if ( (distMin!=kInfinity) && (distMin>dRmax) ) 00574 { // Avoid rounding errors due to precision issues on 00575 // 64 bits systems. Split long distances and recompute 00576 G4double fTerm = distMin-std::fmod(distMin,dRmax); 00577 distMin = fTerm + DistanceToIn(p+fTerm*v,v); 00578 } 00579 } 00580 00581 if (std::fabs(distMin)<halfRadTolerance) { distMin=0.; } 00582 return distMin; 00583 }
G4double G4Ellipsoid::DistanceToOut | ( | const G4ThreeVector & | p | ) | const [virtual] |
Implements G4VSolid.
Definition at line 756 of file G4Ellipsoid.cc.
References G4VSolid::DumpInfo(), G4endl, G4Exception(), Inside(), JustWarning, and kOutside.
00757 { 00758 G4double distR, distZ; 00759 00760 #ifdef G4SPECSDEBUG 00761 if( Inside(p) == kOutside ) 00762 { 00763 DumpInfo(); 00764 std::ostringstream message; 00765 G4int oldprc = message.precision(16); 00766 message << "Point p is outside !?" << G4endl 00767 << "Position:" << G4endl 00768 << " p.x() = " << p.x()/mm << " mm" << G4endl 00769 << " p.y() = " << p.y()/mm << " mm" << G4endl 00770 << " p.z() = " << p.z()/mm << " mm"; 00771 message.precision(oldprc) ; 00772 G4Exception("G4Ellipsoid::DistanceToOut(p)", "GeomSolids1002", 00773 JustWarning, message); 00774 } 00775 #endif 00776 00777 // Normal vector: parallel to normal, magnitude 1/(characteristic radius) 00778 // 00779 G4ThreeVector norm(p.x()/(xSemiAxis*xSemiAxis), 00780 p.y()/(ySemiAxis*ySemiAxis), 00781 p.z()/(zSemiAxis*zSemiAxis)); 00782 00783 // the following is a safe inlined "radius= min(1.0/norm.mag(),p.mag()) 00784 // 00785 G4double radius= p.mag(); 00786 G4double tmp= norm.mag(); 00787 if ( (tmp > 0.0) && (1.0 < radius*tmp) ) {radius = 1.0/tmp;} 00788 00789 // Approximate distance to curved surface ( <= actual distance ) 00790 // 00791 distR = (1.0 - p*norm) * radius / 2.0; 00792 00793 // Distance to z-cut plane 00794 // 00795 distZ = p.z() - zBottomCut; 00796 if (distZ < 0.0) {distZ= zTopCut - p.z();} 00797 00798 // Distance to closest surface from inside 00799 // 00800 if ( (distZ < 0.0) || (distR < 0.0) ) 00801 { 00802 return 0.0; 00803 } 00804 else 00805 { 00806 return (distZ < distR) ? distZ : distR; 00807 } 00808 }
G4double G4Ellipsoid::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 633 of file G4Ellipsoid.cc.
References G4VSolid::DumpInfo(), G4endl, G4Exception(), JustWarning, and sqr().
00638 { 00639 G4double distMin; 00640 enum surface_e {kPlaneSurf, kCurvedSurf, kNoSurf} surface; 00641 00642 distMin= kInfinity; 00643 surface= kNoSurf; 00644 00645 // check to see if Z plane is relevant 00646 // 00647 if (v.z() < 0.0) 00648 { 00649 G4double distZ = (zBottomCut - p.z()) / v.z(); 00650 if (distZ < 0.0) 00651 { 00652 distZ= 0.0; 00653 if (!calcNorm) {return 0.0;} 00654 } 00655 distMin= distZ; 00656 surface= kPlaneSurf; 00657 } 00658 if (v.z() > 0.0) 00659 { 00660 G4double distZ = (zTopCut - p.z()) / v.z(); 00661 if (distZ < 0.0) 00662 { 00663 distZ= 0.0; 00664 if (!calcNorm) {return 0.0;} 00665 } 00666 distMin= distZ; 00667 surface= kPlaneSurf; 00668 } 00669 00670 // normal vector: parallel to normal, magnitude 1/(characteristic radius) 00671 // 00672 G4ThreeVector nearnorm(p.x()/(xSemiAxis*xSemiAxis), 00673 p.y()/(ySemiAxis*ySemiAxis), 00674 p.z()/(zSemiAxis*zSemiAxis)); 00675 00676 // now check curved surface intercept 00677 // 00678 G4double A,B,C; 00679 00680 A= sqr(v.x()/xSemiAxis) + sqr(v.y()/ySemiAxis) + sqr(v.z()/zSemiAxis); 00681 C= (p * nearnorm) - 1.0; 00682 B= 2.0 * (v * nearnorm); 00683 00684 C= B*B - 4.0*A*C; 00685 if (C > 0.0) 00686 { 00687 G4double distR= (-B + std::sqrt(C) ) / (2.0*A); 00688 if (distR < 0.0) 00689 { 00690 distR= 0.0; 00691 if (!calcNorm) {return 0.0;} 00692 } 00693 if (distR < distMin) 00694 { 00695 distMin= distR; 00696 surface= kCurvedSurf; 00697 } 00698 } 00699 00700 // set normal if requested 00701 // 00702 if (calcNorm) 00703 { 00704 if (surface == kNoSurf) 00705 { 00706 *validNorm = false; 00707 } 00708 else 00709 { 00710 *validNorm = true; 00711 switch (surface) 00712 { 00713 case kPlaneSurf: 00714 *n= G4ThreeVector(0.,0.,(v.z() > 0.0 ? 1. : -1.)); 00715 break; 00716 case kCurvedSurf: 00717 { 00718 G4ThreeVector pexit= p + distMin*v; 00719 G4ThreeVector truenorm(pexit.x()/(xSemiAxis*xSemiAxis), 00720 pexit.y()/(ySemiAxis*ySemiAxis), 00721 pexit.z()/(zSemiAxis*zSemiAxis)); 00722 truenorm *= 1.0/truenorm.mag(); 00723 *n= truenorm; 00724 } break; 00725 default: // Should never reach this case ... 00726 DumpInfo(); 00727 std::ostringstream message; 00728 G4int oldprc = message.precision(16); 00729 message << "Undefined side for valid surface normal to solid." 00730 << G4endl 00731 << "Position:" << G4endl 00732 << " p.x() = " << p.x()/mm << " mm" << G4endl 00733 << " p.y() = " << p.y()/mm << " mm" << G4endl 00734 << " p.z() = " << p.z()/mm << " mm" << G4endl 00735 << "Direction:" << G4endl << G4endl 00736 << " v.x() = " << v.x() << G4endl 00737 << " v.y() = " << v.y() << G4endl 00738 << " v.z() = " << v.z() << G4endl 00739 << "Proposed distance :" << G4endl 00740 << " distMin = " << distMin/mm << " mm"; 00741 message.precision(oldprc); 00742 G4Exception("G4Ellipsoid::DistanceToOut(p,v,..)", 00743 "GeomSolids1002", JustWarning, message); 00744 break; 00745 } 00746 } 00747 } 00748 00749 return distMin; 00750 }
G4double G4Ellipsoid::GetCubicVolume | ( | ) | [inline, virtual] |
Reimplemented from G4VSolid.
Definition at line 85 of file G4Ellipsoid.icc.
References G4INCL::Math::pi, and sqr().
00086 { 00087 if(fCubicVolume != 0 ) {;} 00088 else 00089 { 00090 if ((zTopCut > +zSemiAxis && zBottomCut < -zSemiAxis) 00091 || (zTopCut == 0 && zBottomCut == 0) ) 00092 { 00093 fCubicVolume = (4./3.)*CLHEP::pi*xSemiAxis*ySemiAxis*zSemiAxis; 00094 } 00095 else 00096 { 00097 fCubicVolume = CLHEP::pi*xSemiAxis*ySemiAxis 00098 * ((zTopCut-std::pow(zTopCut,3.)/(3.*sqr(zSemiAxis))) 00099 - (zBottomCut-std::pow(zBottomCut,3.)/(3.*sqr(zSemiAxis)))); 00100 } 00101 } 00102 return fCubicVolume; 00103 }
G4GeometryType G4Ellipsoid::GetEntityType | ( | ) | const [virtual] |
Implements G4VSolid.
Definition at line 938 of file G4Ellipsoid.cc.
00939 { 00940 return G4String("G4Ellipsoid"); 00941 }
G4VisExtent G4Ellipsoid::GetExtent | ( | ) | const [virtual] |
Reimplemented from G4VSolid.
Definition at line 1054 of file G4Ellipsoid.cc.
01055 { 01056 // Define the sides of the box into which the G4Ellipsoid instance would fit. 01057 // 01058 return G4VisExtent (-semiAxisMax, semiAxisMax, 01059 -semiAxisMax, semiAxisMax, 01060 -semiAxisMax, semiAxisMax); 01061 }
G4ThreeVector G4Ellipsoid::GetPointOnSurface | ( | ) | const [virtual] |
Reimplemented from G4VSolid.
Definition at line 981 of file G4Ellipsoid.cc.
References G4InuclParticleNames::alpha, G4INCL::Math::pi, and sqr().
00982 { 00983 G4double aTop, aBottom, aCurved, chose, xRand, yRand, zRand, phi; 00984 G4double cosphi, sinphi, costheta, sintheta, alpha, beta, max1, max2, max3; 00985 00986 max1 = xSemiAxis > ySemiAxis ? xSemiAxis : ySemiAxis; 00987 max1 = max1 > zSemiAxis ? max1 : zSemiAxis; 00988 if (max1 == xSemiAxis) { max2 = ySemiAxis; max3 = zSemiAxis; } 00989 else if (max1 == ySemiAxis) { max2 = xSemiAxis; max3 = zSemiAxis; } 00990 else { max2 = xSemiAxis; max3 = ySemiAxis; } 00991 00992 phi = RandFlat::shoot(0.,twopi); 00993 00994 cosphi = std::cos(phi); sinphi = std::sin(phi); 00995 costheta = RandFlat::shoot(zBottomCut,zTopCut)/zSemiAxis; 00996 sintheta = std::sqrt(1.-sqr(costheta)); 00997 00998 alpha = 1.-sqr(max2/max1); beta = 1.-sqr(max3/max1); 00999 01000 aTop = pi*xSemiAxis*ySemiAxis*(1 - sqr(zTopCut/zSemiAxis)); 01001 aBottom = pi*xSemiAxis*ySemiAxis*(1 - sqr(zBottomCut/zSemiAxis)); 01002 01003 // approximation 01004 // from:" http://www.citr.auckland.ac.nz/techreports/2004/CITR-TR-139.pdf" 01005 aCurved = 4.*pi*max1*max2*(1.-1./6.*(alpha+beta)- 01006 1./120.*(3.*sqr(alpha)+2.*alpha*beta+3.*sqr(beta))); 01007 01008 aCurved *= 0.5*(1.2*zTopCut/zSemiAxis - 1.2*zBottomCut/zSemiAxis); 01009 01010 if( ( zTopCut >= zSemiAxis && zBottomCut <= -1.*zSemiAxis ) 01011 || ( zTopCut == 0 && zBottomCut ==0 ) ) 01012 { 01013 aTop = 0; aBottom = 0; 01014 } 01015 01016 chose = RandFlat::shoot(0.,aTop + aBottom + aCurved); 01017 01018 if(chose < aCurved) 01019 { 01020 xRand = xSemiAxis*sintheta*cosphi; 01021 yRand = ySemiAxis*sintheta*sinphi; 01022 zRand = zSemiAxis*costheta; 01023 return G4ThreeVector (xRand,yRand,zRand); 01024 } 01025 else if(chose >= aCurved && chose < aCurved + aTop) 01026 { 01027 xRand = RandFlat::shoot(-1.,1.)*xSemiAxis 01028 * std::sqrt(1-sqr(zTopCut/zSemiAxis)); 01029 yRand = RandFlat::shoot(-1.,1.)*ySemiAxis 01030 * std::sqrt(1.-sqr(zTopCut/zSemiAxis)-sqr(xRand/xSemiAxis)); 01031 zRand = zTopCut; 01032 return G4ThreeVector (xRand,yRand,zRand); 01033 } 01034 else 01035 { 01036 xRand = RandFlat::shoot(-1.,1.)*xSemiAxis 01037 * std::sqrt(1-sqr(zBottomCut/zSemiAxis)); 01038 yRand = RandFlat::shoot(-1.,1.)*ySemiAxis 01039 * std::sqrt(1.-sqr(zBottomCut/zSemiAxis)-sqr(xRand/xSemiAxis)); 01040 zRand = zBottomCut; 01041 return G4ThreeVector (xRand,yRand,zRand); 01042 } 01043 }
G4Polyhedron * G4Ellipsoid::GetPolyhedron | ( | ) | const [virtual] |
Reimplemented from G4VSolid.
Definition at line 1076 of file G4Ellipsoid.cc.
References CreatePolyhedron(), fpPolyhedron, and G4Polyhedron::GetNumberOfRotationStepsAtTimeOfCreation().
01077 { 01078 if (!fpPolyhedron || 01079 fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() != 01080 fpPolyhedron->GetNumberOfRotationSteps()) 01081 { 01082 delete fpPolyhedron; 01083 fpPolyhedron = CreatePolyhedron(); 01084 } 01085 return fpPolyhedron; 01086 }
Definition at line 39 of file G4Ellipsoid.icc.
Referenced by G4GDMLWriteSolids::EllipsoidWrite(), and G4tgbGeometryDumper::GetSolidParams().
G4double G4Ellipsoid::GetSurfaceArea | ( | ) | [inline, virtual] |
Reimplemented from G4VSolid.
Definition at line 106 of file G4Ellipsoid.icc.
References G4VSolid::GetSurfaceArea().
00107 { 00108 if(fSurfaceArea != 0.) {;} 00109 else { fSurfaceArea = G4VSolid::GetSurfaceArea(); } 00110 return fSurfaceArea; 00111 }
G4double G4Ellipsoid::GetZBottomCut | ( | ) | const [inline] |
Definition at line 47 of file G4Ellipsoid.icc.
Referenced by G4GDMLWriteSolids::EllipsoidWrite(), and G4tgbGeometryDumper::GetSolidParams().
G4double G4Ellipsoid::GetZTopCut | ( | ) | const [inline] |
Definition at line 53 of file G4Ellipsoid.icc.
Referenced by G4GDMLWriteSolids::EllipsoidWrite(), and G4tgbGeometryDumper::GetSolidParams().
EInside G4Ellipsoid::Inside | ( | const G4ThreeVector & | p | ) | const [virtual] |
Implements G4VSolid.
Definition at line 416 of file G4Ellipsoid.cc.
References kInside, kOutside, kSurface, and sqr().
Referenced by CalculateExtent(), DistanceToIn(), and DistanceToOut().
00417 { 00418 G4double rad2oo, // outside surface outer tolerance 00419 rad2oi; // outside surface inner tolerance 00420 EInside in; 00421 00422 static const G4double halfRadTolerance=kRadTolerance*0.5; 00423 00424 // check this side of z cut first, because that's fast 00425 // 00426 if (p.z() < zBottomCut-halfRadTolerance) { return in=kOutside; } 00427 if (p.z() > zTopCut+halfRadTolerance) { return in=kOutside; } 00428 00429 rad2oo= sqr(p.x()/(xSemiAxis+halfRadTolerance)) 00430 + sqr(p.y()/(ySemiAxis+halfRadTolerance)) 00431 + sqr(p.z()/(zSemiAxis+halfRadTolerance)); 00432 00433 if (rad2oo > 1.0) { return in=kOutside; } 00434 00435 rad2oi= sqr(p.x()*(1.0+halfRadTolerance/xSemiAxis)/xSemiAxis) 00436 + sqr(p.y()*(1.0+halfRadTolerance/ySemiAxis)/ySemiAxis) 00437 + sqr(p.z()*(1.0+halfRadTolerance/zSemiAxis)/zSemiAxis); 00438 00439 // Check radial surfaces 00440 // sets `in' (already checked for rad2oo > 1.0) 00441 // 00442 if (rad2oi < 1.0) 00443 { 00444 in = ( (p.z() < zBottomCut+halfRadTolerance) 00445 || (p.z() > zTopCut-halfRadTolerance) ) ? kSurface : kInside; 00446 if ( rad2oi > 1.0-halfRadTolerance ) { in=kSurface; } 00447 } 00448 else 00449 { 00450 in = kSurface; 00451 } 00452 return in; 00453 00454 }
G4Ellipsoid & G4Ellipsoid::operator= | ( | const G4Ellipsoid & | rhs | ) |
Definition at line 144 of file G4Ellipsoid.cc.
References fCubicVolume, fpPolyhedron, fSurfaceArea, kRadTolerance, G4VSolid::operator=(), semiAxisMax, xSemiAxis, ySemiAxis, zBottomCut, zSemiAxis, and zTopCut.
00145 { 00146 // Check assignment to self 00147 // 00148 if (this == &rhs) { return *this; } 00149 00150 // Copy base class data 00151 // 00152 G4VSolid::operator=(rhs); 00153 00154 // Copy data 00155 // 00156 fpPolyhedron = 0; kRadTolerance = rhs.kRadTolerance; 00157 fCubicVolume = rhs.fCubicVolume; fSurfaceArea = rhs.fSurfaceArea; 00158 xSemiAxis = rhs.xSemiAxis; ySemiAxis = rhs.ySemiAxis; 00159 zSemiAxis = rhs.zSemiAxis; semiAxisMax = rhs.semiAxisMax; 00160 zBottomCut = rhs.zBottomCut; zTopCut = rhs.zTopCut; 00161 00162 return *this; 00163 }
Definition at line 59 of file G4Ellipsoid.icc.
Referenced by G4Ellipsoid().
00062 { 00063 xSemiAxis= newxSemiAxis; ySemiAxis= newySemiAxis; zSemiAxis= newzSemiAxis; 00064 semiAxisMax = xSemiAxis > ySemiAxis ? xSemiAxis : ySemiAxis; 00065 if (zSemiAxis > semiAxisMax) { semiAxisMax= zSemiAxis; } 00066 if (zBottomCut < -zSemiAxis) { zBottomCut = -zSemiAxis; } 00067 if (zTopCut > +zSemiAxis) { zTopCut = +zSemiAxis; } 00068 }
Definition at line 71 of file G4Ellipsoid.icc.
Referenced by G4Ellipsoid().
00072 { 00073 if (newzBottomCut < -zSemiAxis) 00074 { zBottomCut = -zSemiAxis; } 00075 else 00076 { zBottomCut = newzBottomCut; } 00077 00078 if (newzTopCut > +zSemiAxis) 00079 { zTopCut = +zSemiAxis; } 00080 else 00081 { zTopCut = newzTopCut; } 00082 }
std::ostream & G4Ellipsoid::StreamInfo | ( | std::ostream & | os | ) | const [virtual] |
Implements G4VSolid.
Definition at line 956 of file G4Ellipsoid.cc.
References G4VSolid::GetName().
00957 { 00958 G4int oldprc = os.precision(16); 00959 os << "-----------------------------------------------------------\n" 00960 << " *** Dump for solid - " << GetName() << " ***\n" 00961 << " ===================================================\n" 00962 << " Solid type: G4Ellipsoid\n" 00963 << " Parameters: \n" 00964 00965 << " semi-axis x: " << xSemiAxis/mm << " mm \n" 00966 << " semi-axis y: " << ySemiAxis/mm << " mm \n" 00967 << " semi-axis z: " << zSemiAxis/mm << " mm \n" 00968 << " max semi-axis: " << semiAxisMax/mm << " mm \n" 00969 << " lower cut plane level z: " << zBottomCut/mm << " mm \n" 00970 << " upper cut plane level z: " << zTopCut/mm << " mm \n" 00971 << "-----------------------------------------------------------\n"; 00972 os.precision(oldprc); 00973 00974 return os; 00975 }
G4ThreeVector G4Ellipsoid::SurfaceNormal | ( | const G4ThreeVector & | p | ) | const [virtual] |
Implements G4VSolid.
Definition at line 460 of file G4Ellipsoid.cc.
Referenced by DistanceToIn().
00461 { 00462 G4double distR, distZBottom, distZTop; 00463 00464 // normal vector with special magnitude: parallel to normal, units 1/length 00465 // norm*p == 1.0 if on surface, >1.0 if outside, <1.0 if inside 00466 // 00467 G4ThreeVector norm(p.x()/(xSemiAxis*xSemiAxis), 00468 p.y()/(ySemiAxis*ySemiAxis), 00469 p.z()/(zSemiAxis*zSemiAxis)); 00470 G4double radius = 1.0/norm.mag(); 00471 00472 // approximate distance to curved surface 00473 // 00474 distR = std::fabs( (p*norm - 1.0) * radius ) / 2.0; 00475 00476 // Distance to z-cut plane 00477 // 00478 distZBottom = std::fabs( p.z() - zBottomCut ); 00479 distZTop = std::fabs( p.z() - zTopCut ); 00480 00481 if ( (distZBottom < distR) || (distZTop < distR) ) 00482 { 00483 return G4ThreeVector(0.,0.,(distZBottom < distZTop) ? -1.0 : 1.0); 00484 } 00485 return ( norm *= radius ); 00486 }
G4Polyhedron* G4Ellipsoid::fpPolyhedron [mutable, protected] |