#include <G4EllipticalCone.hh>
Inheritance diagram for G4EllipticalCone:
Definition at line 84 of file G4EllipticalCone.hh.
G4EllipticalCone::G4EllipticalCone | ( | const G4String & | pName, | |
G4double | pxSemiAxis, | |||
G4double | pySemiAxis, | |||
G4double | zMax, | |||
G4double | pzTopCut | |||
) |
Definition at line 68 of file G4EllipticalCone.cc.
References FatalErrorInArgument, G4Exception(), G4GeometryTolerance::GetInstance(), G4VSolid::GetName(), G4GeometryTolerance::GetRadialTolerance(), SetSemiAxis(), and SetZCut().
Referenced by Clone().
00073 : G4VSolid(pName), fpPolyhedron(0), fCubicVolume(0.), fSurfaceArea(0.), 00074 zTopCut(0.) 00075 { 00076 00077 kRadTolerance = G4GeometryTolerance::GetInstance()->GetRadialTolerance(); 00078 00079 // Check Semi-Axis & Z-cut 00080 // 00081 if ( (pxSemiAxis <= 0.) || (pySemiAxis <= 0.) || (pzMax <= 0.) ) 00082 { 00083 std::ostringstream message; 00084 message << "Invalid semi-axis or height - " << GetName(); 00085 G4Exception("G4EllipticalCone::G4EllipticalCone()", "GeomSolids0002", 00086 FatalErrorInArgument, message); 00087 } 00088 if ( pzTopCut <= 0 ) 00089 { 00090 std::ostringstream message; 00091 message << "Invalid z-coordinate for cutting plane - " << GetName(); 00092 G4Exception("G4EllipticalCone::G4EllipticalCone()", "InvalidSetup", 00093 FatalErrorInArgument, message); 00094 } 00095 00096 SetSemiAxis( pxSemiAxis, pySemiAxis, pzMax ); 00097 SetZCut(pzTopCut); 00098 }
G4EllipticalCone::~G4EllipticalCone | ( | ) | [virtual] |
G4EllipticalCone::G4EllipticalCone | ( | __void__ & | ) |
Definition at line 105 of file G4EllipticalCone.cc.
00106 : G4VSolid(a), fpPolyhedron(0), kRadTolerance(0.), fCubicVolume(0.), 00107 fSurfaceArea(0.), xSemiAxis(0.), ySemiAxis(0.), zheight(0.), 00108 semiAxisMax(0.), zTopCut(0.) 00109 { 00110 }
G4EllipticalCone::G4EllipticalCone | ( | const G4EllipticalCone & | rhs | ) |
Definition at line 124 of file G4EllipticalCone.cc.
00125 : G4VSolid(rhs), 00126 fpPolyhedron(0), kRadTolerance(rhs.kRadTolerance), 00127 fCubicVolume(rhs.fCubicVolume), fSurfaceArea(rhs.fSurfaceArea), 00128 xSemiAxis(rhs.xSemiAxis), ySemiAxis(rhs.ySemiAxis), zheight(rhs.zheight), 00129 semiAxisMax(rhs.semiAxisMax), zTopCut(rhs.zTopCut) 00130 { 00131 }
G4bool G4EllipticalCone::CalculateExtent | ( | const EAxis | pAxis, | |
const G4VoxelLimits & | pVoxelLimit, | |||
const G4AffineTransform & | pTransform, | |||
G4double & | pmin, | |||
G4double & | pmax | |||
) | const [virtual] |
Implements G4VSolid.
Definition at line 162 of file G4EllipticalCone.cc.
References G4SolidExtentList::AddSurface(), G4ClippablePolygon::AddVertexInOrder(), G4AffineTransform::ApplyPointTransform(), G4ClippablePolygon::ClearAllVertices(), G4SolidExtentList::GetExtent(), kMaxMeshSections, G4ClippablePolygon::PartialClip(), G4ClippablePolygon::SetNormal(), and G4AffineTransform::TransformAxis().
00166 { 00167 G4SolidExtentList extentList( axis, voxelLimit ); 00168 00169 // 00170 // We are going to divide up our elliptical face into small pieces 00171 // 00172 00173 // 00174 // Choose phi size of our segment(s) based on constants as 00175 // defined in meshdefs.hh 00176 // 00177 G4int numPhi = kMaxMeshSections; 00178 G4double sigPhi = twopi/numPhi; 00179 00180 // 00181 // We have to be careful to keep our segments completely outside 00182 // of the elliptical surface. To do so we imagine we have 00183 // a simple (unit radius) circular cross section (as in G4Tubs) 00184 // and then "stretch" the dimensions as necessary to fit the ellipse. 00185 // 00186 G4double rFudge = 1.0/std::cos(0.5*sigPhi); 00187 G4double dxFudgeBot = xSemiAxis*2.*zheight*rFudge, 00188 dyFudgeBot = ySemiAxis*2.*zheight*rFudge; 00189 G4double dxFudgeTop = xSemiAxis*(zheight-zTopCut)*rFudge, 00190 dyFudgeTop = ySemiAxis*(zheight-zTopCut)*rFudge; 00191 00192 // 00193 // As we work around the elliptical surface, we build 00194 // a "phi" segment on the way, and keep track of two 00195 // additional polygons for the two ends. 00196 // 00197 G4ClippablePolygon endPoly1, endPoly2, phiPoly; 00198 00199 G4double phi = 0, 00200 cosPhi = std::cos(phi), 00201 sinPhi = std::sin(phi); 00202 G4ThreeVector v0( dxFudgeTop*cosPhi, dyFudgeTop*sinPhi, +zTopCut ), 00203 v1( dxFudgeBot*cosPhi, dyFudgeBot*sinPhi, -zTopCut ), 00204 w0, w1; 00205 transform.ApplyPointTransform( v0 ); 00206 transform.ApplyPointTransform( v1 ); 00207 do 00208 { 00209 phi += sigPhi; 00210 if (numPhi == 1) phi = 0; // Try to avoid roundoff 00211 cosPhi = std::cos(phi), 00212 sinPhi = std::sin(phi); 00213 00214 w0 = G4ThreeVector( dxFudgeTop*cosPhi, dyFudgeTop*sinPhi, +zTopCut ); 00215 w1 = G4ThreeVector( dxFudgeBot*cosPhi, dyFudgeBot*sinPhi, -zTopCut ); 00216 transform.ApplyPointTransform( w0 ); 00217 transform.ApplyPointTransform( w1 ); 00218 00219 // 00220 // Add a point to our z ends 00221 // 00222 endPoly1.AddVertexInOrder( v0 ); 00223 endPoly2.AddVertexInOrder( v1 ); 00224 00225 // 00226 // Build phi polygon 00227 // 00228 phiPoly.ClearAllVertices(); 00229 00230 phiPoly.AddVertexInOrder( v0 ); 00231 phiPoly.AddVertexInOrder( v1 ); 00232 phiPoly.AddVertexInOrder( w1 ); 00233 phiPoly.AddVertexInOrder( w0 ); 00234 00235 if (phiPoly.PartialClip( voxelLimit, axis )) 00236 { 00237 // 00238 // Get unit normal 00239 // 00240 phiPoly.SetNormal( (v1-v0).cross(w0-v0).unit() ); 00241 00242 extentList.AddSurface( phiPoly ); 00243 } 00244 00245 // 00246 // Next vertex 00247 // 00248 v0 = w0; 00249 v1 = w1; 00250 } while( --numPhi > 0 ); 00251 00252 // 00253 // Process the end pieces 00254 // 00255 if (endPoly1.PartialClip( voxelLimit, axis )) 00256 { 00257 static const G4ThreeVector normal(0,0,+1); 00258 endPoly1.SetNormal( transform.TransformAxis(normal) ); 00259 extentList.AddSurface( endPoly1 ); 00260 } 00261 00262 if (endPoly2.PartialClip( voxelLimit, axis )) 00263 { 00264 static const G4ThreeVector normal(0,0,-1); 00265 endPoly2.SetNormal( transform.TransformAxis(normal) ); 00266 extentList.AddSurface( endPoly2 ); 00267 } 00268 00269 // 00270 // Return min/max value 00271 // 00272 return extentList.GetExtent( min, max ); 00273 }
G4VSolid * G4EllipticalCone::Clone | ( | ) | const [virtual] |
Reimplemented from G4VSolid.
Definition at line 960 of file G4EllipticalCone.cc.
References G4EllipticalCone().
00961 { 00962 return new G4EllipticalCone(*this); 00963 }
G4NURBS * G4EllipticalCone::CreateNURBS | ( | ) | const [virtual] |
Reimplemented from G4VSolid.
Definition at line 1070 of file G4EllipticalCone.cc.
01071 { 01072 // Box for now!!! 01073 // 01074 return new G4NURBSbox(xSemiAxis, ySemiAxis,zheight); 01075 }
G4Polyhedron * G4EllipticalCone::CreatePolyhedron | ( | ) | const [virtual] |
Reimplemented from G4VSolid.
Definition at line 1077 of file G4EllipticalCone.cc.
Referenced by GetPolyhedron().
01078 { 01079 return new G4PolyhedronEllipticalCone(xSemiAxis, ySemiAxis, zheight, zTopCut); 01080 }
G4ThreeVectorList* G4EllipticalCone::CreateRotatedVertices | ( | const G4AffineTransform & | pT, | |
G4int & | noPV | |||
) | const [protected] |
void G4EllipticalCone::DescribeYourselfTo | ( | G4VGraphicsScene & | scene | ) | const [virtual] |
Implements G4VSolid.
Definition at line 1052 of file G4EllipticalCone.cc.
References G4VGraphicsScene::AddSolid().
01053 { 01054 scene.AddSolid(*this); 01055 }
G4double G4EllipticalCone::DistanceToIn | ( | const G4ThreeVector & | p | ) | const [virtual] |
Implements G4VSolid.
Definition at line 681 of file G4EllipticalCone.cc.
References G4VSolid::kCarTolerance, and sqr().
00682 { 00683 G4double distR, distR2, distZ, maxDim; 00684 G4double distRad; 00685 00686 // check if the point lies either below z=-zTopCut in bottom elliptical 00687 // region or on top within cut elliptical region 00688 // 00689 if( (p.z() <= -zTopCut) && (sqr(p.x()/xSemiAxis) + sqr(p.y()/ySemiAxis) 00690 <= sqr(zTopCut + zheight + 0.5*kCarTolerance )) ) 00691 { 00692 //return distZ = std::fabs(zTopCut - p.z()); 00693 return distZ = std::fabs(zTopCut + p.z()); 00694 } 00695 00696 if( (p.z() >= zTopCut) && (sqr(p.x()/xSemiAxis)+sqr(p.y()/ySemiAxis) 00697 <= sqr(zheight - zTopCut + kCarTolerance/2.0 )) ) 00698 { 00699 return distZ = std::fabs(p.z() - zTopCut); 00700 } 00701 00702 // below we use the following approximation: we take the largest of the 00703 // axes and find the shortest distance to the circular (cut) cone of that 00704 // radius. 00705 // 00706 maxDim = xSemiAxis >= ySemiAxis ? xSemiAxis:ySemiAxis; 00707 distRad = std::sqrt(p.x()*p.x()+p.y()*p.y()); 00708 00709 if( p.z() > maxDim*distRad + zTopCut*(1.+maxDim)-sqr(maxDim)*zheight ) 00710 { 00711 distR2 = sqr(p.z() - zTopCut) + sqr(distRad - maxDim*(zheight - zTopCut)); 00712 return std::sqrt( distR2 ); 00713 } 00714 00715 if( distRad > maxDim*( zheight - p.z() ) ) 00716 { 00717 if( p.z() > maxDim*distRad - (zTopCut*(1.+maxDim)+sqr(maxDim)*zheight) ) 00718 { 00719 G4double zVal = (p.z()-maxDim*(distRad-maxDim*zheight))/(1.+sqr(maxDim)); 00720 G4double rVal = maxDim*(zheight - zVal); 00721 return distR = std::sqrt(sqr(p.z() - zVal) + sqr(distRad - rVal)); 00722 } 00723 } 00724 00725 if( distRad <= maxDim*(zheight - p.z()) ) 00726 { 00727 distR2 = sqr(distRad - maxDim*(zheight + zTopCut)) + sqr(p.z() + zTopCut); 00728 return std::sqrt( distR2 ); 00729 } 00730 00731 return distR = 0; 00732 }
G4double G4EllipticalCone::DistanceToIn | ( | const G4ThreeVector & | p, | |
const G4ThreeVector & | v | |||
) | const [virtual] |
Implements G4VSolid.
Definition at line 430 of file G4EllipticalCone.cc.
References G4VSolid::kCarTolerance, G4InuclParticleNames::lambda, and sqr().
00432 { 00433 static const G4double halfTol = 0.5*kCarTolerance; 00434 00435 G4double distMin = kInfinity; 00436 00437 // code from EllipticalTube 00438 00439 G4double sigz = p.z()+zTopCut; 00440 00441 // 00442 // Check z = -dz planer surface 00443 // 00444 00445 if (sigz < halfTol) 00446 { 00447 // 00448 // We are "behind" the shape in z, and so can 00449 // potentially hit the rear face. Correct direction? 00450 // 00451 if (v.z() <= 0) 00452 { 00453 // 00454 // As long as we are far enough away, we know we 00455 // can't intersect 00456 // 00457 if (sigz < 0) return kInfinity; 00458 00459 // 00460 // Otherwise, we don't intersect unless we are 00461 // on the surface of the ellipse 00462 // 00463 00464 if ( sqr(p.x()/( xSemiAxis - halfTol )) 00465 + sqr(p.y()/( ySemiAxis - halfTol )) <= sqr( zheight+zTopCut ) ) 00466 return kInfinity; 00467 00468 } 00469 else 00470 { 00471 // 00472 // How far? 00473 // 00474 G4double q = -sigz/v.z(); 00475 00476 // 00477 // Where does that place us? 00478 // 00479 G4double xi = p.x() + q*v.x(), 00480 yi = p.y() + q*v.y(); 00481 00482 // 00483 // Is this on the surface (within ellipse)? 00484 // 00485 if ( sqr(xi/xSemiAxis) + sqr(yi/ySemiAxis) <= sqr( zheight + zTopCut ) ) 00486 { 00487 // 00488 // Yup. Return q, unless we are on the surface 00489 // 00490 return (sigz < -halfTol) ? q : 0; 00491 } 00492 else if (xi/(xSemiAxis*xSemiAxis)*v.x() 00493 + yi/(ySemiAxis*ySemiAxis)*v.y() >= 0) 00494 { 00495 // 00496 // Else, if we are traveling outwards, we know 00497 // we must miss 00498 // 00499 // return kInfinity; 00500 } 00501 } 00502 } 00503 00504 // 00505 // Check z = +dz planer surface 00506 // 00507 sigz = p.z() - zTopCut; 00508 00509 if (sigz > -halfTol) 00510 { 00511 if (v.z() >= 0) 00512 { 00513 00514 if (sigz > 0) return kInfinity; 00515 00516 if ( sqr(p.x()/( xSemiAxis - halfTol )) 00517 + sqr(p.y()/( ySemiAxis - halfTol )) <= sqr( zheight-zTopCut ) ) 00518 return kInfinity; 00519 00520 } 00521 else { 00522 G4double q = -sigz/v.z(); 00523 00524 G4double xi = p.x() + q*v.x(), 00525 yi = p.y() + q*v.y(); 00526 00527 if ( sqr(xi/xSemiAxis) + sqr(yi/ySemiAxis) <= sqr( zheight - zTopCut ) ) 00528 { 00529 return (sigz > -halfTol) ? q : 0; 00530 } 00531 else if (xi/(xSemiAxis*xSemiAxis)*v.x() 00532 + yi/(ySemiAxis*ySemiAxis)*v.y() >= 0) 00533 { 00534 // return kInfinity; 00535 } 00536 } 00537 } 00538 00539 00540 #if 0 00541 00542 // check to see if Z plane is relevant 00543 // 00544 if (p.z() < -zTopCut - 0.5*kCarTolerance) 00545 { 00546 if (v.z() <= 0.0) 00547 return distMin; 00548 00549 G4double lambda = (-zTopCut - p.z())/v.z(); 00550 00551 if ( sqr((lambda*v.x()+p.x())/xSemiAxis) + 00552 sqr((lambda*v.y()+p.y())/ySemiAxis) <= 00553 sqr(zTopCut + zheight + 0.5*kRadTolerance) ) 00554 { 00555 return distMin = std::fabs(lambda); 00556 } 00557 } 00558 00559 if (p.z() > zTopCut+0.5*kCarTolerance) 00560 { 00561 if (v.z() >= 0.0) 00562 { return distMin; } 00563 00564 G4double lambda = (zTopCut - p.z()) / v.z(); 00565 00566 if ( sqr((lambda*v.x() + p.x())/xSemiAxis) + 00567 sqr((lambda*v.y() + p.y())/ySemiAxis) <= 00568 sqr(zheight - zTopCut + 0.5*kRadTolerance) ) 00569 { 00570 return distMin = std::fabs(lambda); 00571 } 00572 } 00573 00574 if (p.z() > zTopCut - halfTol 00575 && p.z() < zTopCut + halfTol ) 00576 { 00577 if (v.z() > 0.) 00578 { return kInfinity; } 00579 00580 return distMin = 0.; 00581 } 00582 00583 if (p.z() < -zTopCut + halfTol 00584 && p.z() > -zTopCut - halfTol) 00585 { 00586 if (v.z() < 0.) 00587 { return distMin = kInfinity; } 00588 00589 return distMin = 0.; 00590 } 00591 00592 #endif 00593 00594 // if we are here then it either intersects or grazes the curved surface 00595 // or it does not intersect at all 00596 // 00597 G4double A = sqr(v.x()/xSemiAxis) + sqr(v.y()/ySemiAxis) - sqr(v.z()); 00598 G4double B = 2*(v.x()*p.x()/sqr(xSemiAxis) + 00599 v.y()*p.y()/sqr(ySemiAxis) + v.z()*(zheight-p.z())); 00600 G4double C = sqr(p.x()/xSemiAxis) + sqr(p.y()/ySemiAxis) - 00601 sqr(zheight - p.z()); 00602 00603 G4double discr = B*B - 4.*A*C; 00604 00605 // if the discriminant is negative it never hits the curved object 00606 // 00607 if ( discr < -halfTol ) 00608 { return distMin; } 00609 00610 // case below is when it hits or grazes the surface 00611 // 00612 if ( (discr >= - halfTol ) && (discr < halfTol ) ) 00613 { 00614 return distMin = std::fabs(-B/(2.*A)); 00615 } 00616 00617 G4double plus = (-B+std::sqrt(discr))/(2.*A); 00618 G4double minus = (-B-std::sqrt(discr))/(2.*A); 00619 00620 // Special case::Point on Surface, Check norm.dot(v) 00621 00622 if ( ( std::fabs(plus) < halfTol )||( std::fabs(minus) < halfTol ) ) 00623 { 00624 G4ThreeVector truenorm(p.x()/(xSemiAxis*xSemiAxis), 00625 p.y()/(ySemiAxis*ySemiAxis), 00626 -( p.z() - zheight )); 00627 if ( truenorm*v >= 0) // going outside the solid from surface 00628 { 00629 return kInfinity; 00630 } 00631 else 00632 { 00633 return 0; 00634 } 00635 } 00636 00637 // G4double lambda = std::fabs(plus) < std::fabs(minus) ? plus : minus; 00638 G4double lambda = 0; 00639 00640 if ( minus > halfTol && minus < distMin ) 00641 { 00642 lambda = minus ; 00643 // check normal vector n * v < 0 00644 G4ThreeVector pin = p + lambda*v; 00645 if(std::fabs(pin.z())<zTopCut+0.5*kCarTolerance) 00646 { 00647 G4ThreeVector truenorm(pin.x()/(xSemiAxis*xSemiAxis), 00648 pin.y()/(ySemiAxis*ySemiAxis), 00649 - ( pin.z() - zheight )); 00650 if ( truenorm*v < 0) 00651 { // yes, going inside the solid 00652 distMin = lambda; 00653 } 00654 } 00655 } 00656 if ( plus > halfTol && plus < distMin ) 00657 { 00658 lambda = plus ; 00659 // check normal vector n * v < 0 00660 G4ThreeVector pin = p + lambda*v; 00661 if(std::fabs(pin.z())<zTopCut+0.5*kCarTolerance) 00662 { 00663 G4ThreeVector truenorm(pin.x()/(xSemiAxis*xSemiAxis), 00664 pin.y()/(ySemiAxis*ySemiAxis), 00665 - ( pin.z() - zheight ) ); 00666 if ( truenorm*v < 0) 00667 { // yes, going inside the solid 00668 distMin = lambda; 00669 } 00670 } 00671 } 00672 if (distMin < halfTol) distMin=0.; 00673 return distMin ; 00674 }
G4double G4EllipticalCone::DistanceToOut | ( | const G4ThreeVector & | p | ) | const [virtual] |
Implements G4VSolid.
Definition at line 901 of file G4EllipticalCone.cc.
References G4VSolid::DumpInfo(), G4endl, G4Exception(), Inside(), JustWarning, kOutside, and sqr().
00902 { 00903 G4double rds,roo,roo1, distR, distZ, distMin=0.; 00904 G4double minAxis = xSemiAxis < ySemiAxis ? xSemiAxis : ySemiAxis; 00905 00906 #ifdef G4SPECSDEBUG 00907 if( Inside(p) == kOutside ) 00908 { 00909 DumpInfo(); 00910 std::ostringstream message; 00911 G4int oldprc = message.precision(16); 00912 message << "Point p is outside !?" << G4endl 00913 << "Position:" << G4endl 00914 << " p.x() = " << p.x()/mm << " mm" << G4endl 00915 << " p.y() = " << p.y()/mm << " mm" << G4endl 00916 << " p.z() = " << p.z()/mm << " mm"; 00917 message.precision(oldprc) ; 00918 G4Exception("G4Ellipsoid::DistanceToOut(p)", "GeomSolids1002", 00919 JustWarning, message); 00920 } 00921 #endif 00922 00923 // since we have made the above warning, below we are working assuming p 00924 // is inside check how close it is to the circular cone with radius equal 00925 // to the smaller of the axes 00926 // 00927 if( sqr(p.x()/minAxis)+sqr(p.y()/minAxis) < sqr(zheight - p.z()) ) 00928 { 00929 rds = std::sqrt(sqr(p.x()) + sqr(p.y())); 00930 roo = minAxis*(zheight-p.z()); // radius of cone at z= p.z() 00931 roo1 = minAxis*(zheight-zTopCut); // radius of cone at z=+zTopCut 00932 00933 distZ=zTopCut - std::fabs(p.z()) ; 00934 distR=(roo-rds)/(std::sqrt(1+sqr(minAxis))); 00935 00936 if(rds>roo1) 00937 { 00938 distMin=(zTopCut-p.z())*(roo-rds)/(roo-roo1); 00939 distMin=std::min(distMin,distR); 00940 } 00941 distMin=std::min(distR,distZ); 00942 } 00943 00944 return distMin; 00945 }
G4double G4EllipticalCone::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 739 of file G4EllipticalCone.cc.
References G4VSolid::DumpInfo(), G4endl, G4Exception(), JustWarning, G4VSolid::kCarTolerance, G4InuclParticleNames::lambda, and sqr().
00744 { 00745 G4double distMin, lambda; 00746 enum surface_e {kPlaneSurf, kCurvedSurf, kNoSurf} surface; 00747 00748 distMin = kInfinity; 00749 surface = kNoSurf; 00750 00751 if (v.z() < 0.0) 00752 { 00753 lambda = (-p.z() - zTopCut)/v.z(); 00754 00755 if ( (sqr((p.x() + lambda*v.x())/xSemiAxis) + 00756 sqr((p.y() + lambda*v.y())/ySemiAxis)) < 00757 sqr(zheight + zTopCut + 0.5*kCarTolerance) ) 00758 { 00759 distMin = std::fabs(lambda); 00760 00761 if (!calcNorm) { return distMin; } 00762 } 00763 distMin = std::fabs(lambda); 00764 surface = kPlaneSurf; 00765 } 00766 00767 if (v.z() > 0.0) 00768 { 00769 lambda = (zTopCut - p.z()) / v.z(); 00770 00771 if ( (sqr((p.x() + lambda*v.x())/xSemiAxis) 00772 + sqr((p.y() + lambda*v.y())/ySemiAxis) ) 00773 < (sqr(zheight - zTopCut + 0.5*kCarTolerance)) ) 00774 { 00775 distMin = std::fabs(lambda); 00776 if (!calcNorm) { return distMin; } 00777 } 00778 distMin = std::fabs(lambda); 00779 surface = kPlaneSurf; 00780 } 00781 00782 // if we are here then it either intersects or grazes the 00783 // curved surface... 00784 // 00785 G4double A = sqr(v.x()/xSemiAxis) + sqr(v.y()/ySemiAxis) - sqr(v.z()); 00786 G4double B = 2.*(v.x()*p.x()/sqr(xSemiAxis) + 00787 v.y()*p.y()/sqr(ySemiAxis) + v.z()*(zheight-p.z())); 00788 G4double C = sqr(p.x()/xSemiAxis) + sqr(p.y()/ySemiAxis) 00789 - sqr(zheight - p.z()); 00790 00791 G4double discr = B*B - 4.*A*C; 00792 00793 if ( discr >= - 0.5*kCarTolerance && discr < 0.5*kCarTolerance ) 00794 { 00795 if(!calcNorm) { return distMin = std::fabs(-B/(2.*A)); } 00796 } 00797 00798 else if ( discr > 0.5*kCarTolerance ) 00799 { 00800 G4double plus = (-B+std::sqrt(discr))/(2.*A); 00801 G4double minus = (-B-std::sqrt(discr))/(2.*A); 00802 00803 if ( plus > 0.5*kCarTolerance && minus > 0.5*kCarTolerance ) 00804 { 00805 // take the shorter distance 00806 // 00807 lambda = std::fabs(plus) < std::fabs(minus) ? plus : minus; 00808 } 00809 else 00810 { 00811 // at least one solution is close to zero or negative 00812 // so, take small positive solution or zero 00813 // 00814 lambda = plus > -0.5*kCarTolerance ? plus : 0; 00815 } 00816 00817 if ( std::fabs(lambda) < distMin ) 00818 { 00819 if( std::fabs(lambda) > 0.5*kCarTolerance) 00820 { 00821 distMin = std::fabs(lambda); 00822 surface = kCurvedSurf; 00823 } 00824 else // Point is On the Surface, Check Normal 00825 { 00826 G4ThreeVector truenorm(p.x()/(xSemiAxis*xSemiAxis), 00827 p.y()/(ySemiAxis*ySemiAxis), 00828 -( p.z() - zheight )); 00829 if( truenorm.dot(v) > 0 ) 00830 { 00831 distMin = 0.0; 00832 surface = kCurvedSurf; 00833 } 00834 } 00835 } 00836 } 00837 00838 // set normal if requested 00839 // 00840 if (calcNorm) 00841 { 00842 if (surface == kNoSurf) 00843 { 00844 *validNorm = false; 00845 } 00846 else 00847 { 00848 *validNorm = true; 00849 switch (surface) 00850 { 00851 case kPlaneSurf: 00852 { 00853 *n = G4ThreeVector(0.,0.,(v.z() > 0.0 ? 1. : -1.)); 00854 } 00855 break; 00856 00857 case kCurvedSurf: 00858 { 00859 G4ThreeVector pexit = p + distMin*v; 00860 G4ThreeVector truenorm( pexit.x()/(xSemiAxis*xSemiAxis), 00861 pexit.y()/(ySemiAxis*ySemiAxis), 00862 -( pexit.z() - zheight ) ); 00863 truenorm /= truenorm.mag(); 00864 *n= truenorm; 00865 } 00866 break; 00867 00868 default: // Should never reach this case ... 00869 DumpInfo(); 00870 std::ostringstream message; 00871 G4int oldprc = message.precision(16); 00872 message << "Undefined side for valid surface normal to solid." 00873 << G4endl 00874 << "Position:" << G4endl 00875 << " p.x() = " << p.x()/mm << " mm" << G4endl 00876 << " p.y() = " << p.y()/mm << " mm" << G4endl 00877 << " p.z() = " << p.z()/mm << " mm" << G4endl 00878 << "Direction:" << G4endl 00879 << " v.x() = " << v.x() << G4endl 00880 << " v.y() = " << v.y() << G4endl 00881 << " v.z() = " << v.z() << G4endl 00882 << "Proposed distance :" << G4endl 00883 << " distMin = " << distMin/mm << " mm"; 00884 message.precision(oldprc); 00885 G4Exception("G4EllipticalCone::DistanceToOut(p,v,..)", 00886 "GeomSolids1002", JustWarning, message); 00887 break; 00888 } 00889 } 00890 } 00891 00892 if (distMin<0.5*kCarTolerance) { distMin=0; } 00893 00894 return distMin; 00895 }
G4double G4EllipticalCone::GetCubicVolume | ( | ) | [inline, virtual] |
Reimplemented from G4VSolid.
Definition at line 92 of file G4EllipticalCone.icc.
References G4INCL::Math::pi, and sqr().
00093 { 00094 if(fCubicVolume != 0 ) {;} 00095 else 00096 { 00097 if (zTopCut > +zheight ) 00098 { 00099 fCubicVolume = (8./3.)*CLHEP::pi*xSemiAxis* 00100 ySemiAxis*zheight*zheight*zheight; 00101 } 00102 else 00103 { 00104 fCubicVolume = CLHEP::pi*xSemiAxis*ySemiAxis* 00105 (2./3.*std::pow(zTopCut,3.)+2.*sqr(zheight)*zTopCut); 00106 } 00107 } 00108 return fCubicVolume; 00109 }
G4GeometryType G4EllipticalCone::GetEntityType | ( | ) | const [virtual] |
Implements G4VSolid.
Definition at line 951 of file G4EllipticalCone.cc.
00952 { 00953 return G4String("G4EllipticalCone"); 00954 }
G4VisExtent G4EllipticalCone::GetExtent | ( | ) | const [virtual] |
Reimplemented from G4VSolid.
Definition at line 1057 of file G4EllipticalCone.cc.
01058 { 01059 // Define the sides of the box into which the solid instance would fit. 01060 // 01061 G4double maxDim; 01062 maxDim = xSemiAxis > ySemiAxis ? xSemiAxis : ySemiAxis; 01063 maxDim = maxDim > zTopCut ? maxDim : zTopCut; 01064 01065 return G4VisExtent (-maxDim, maxDim, 01066 -maxDim, maxDim, 01067 -maxDim, maxDim); 01068 }
G4ThreeVector G4EllipticalCone::GetPointOnSurface | ( | ) | const [virtual] |
Reimplemented from G4VSolid.
Definition at line 994 of file G4EllipticalCone.cc.
References G4INCL::Math::pi, and sqr().
00995 { 00996 00997 G4double phi, sinphi, cosphi, aOne, aTwo, aThree, 00998 chose, zRand, rRand1, rRand2; 00999 01000 G4double rOne = std::sqrt(sqr(xSemiAxis) 01001 + sqr(ySemiAxis))*(zheight - zTopCut); 01002 G4double rTwo = std::sqrt(sqr(xSemiAxis) 01003 + sqr(ySemiAxis))*(zheight + zTopCut); 01004 01005 aOne = pi*(rOne + rTwo)*std::sqrt(sqr(rOne - rTwo)+sqr(2.*zTopCut)); 01006 aTwo = pi*xSemiAxis*ySemiAxis*sqr(zheight+zTopCut); 01007 aThree = pi*xSemiAxis*ySemiAxis*sqr(zheight-zTopCut); 01008 01009 phi = RandFlat::shoot(0.,twopi); 01010 cosphi = std::cos(phi); 01011 sinphi = std::sin(phi); 01012 01013 if(zTopCut >= zheight) aThree = 0.; 01014 01015 chose = RandFlat::shoot(0.,aOne+aTwo+aThree); 01016 if((chose>=0.) && (chose<aOne)) 01017 { 01018 zRand = RandFlat::shoot(-zTopCut,zTopCut); 01019 return G4ThreeVector(xSemiAxis*(zheight-zRand)*cosphi, 01020 ySemiAxis*(zheight-zRand)*sinphi,zRand); 01021 } 01022 else if((chose>=aOne) && (chose<aOne+aTwo)) 01023 { 01024 do 01025 { 01026 rRand1 = RandFlat::shoot(0.,1.) ; 01027 rRand2 = RandFlat::shoot(0.,1.) ; 01028 } while ( rRand2 >= rRand1 ) ; 01029 01030 // rRand2 = RandFlat::shoot(0.,std::sqrt(1.-sqr(rRand1))); 01031 return G4ThreeVector(rRand1*xSemiAxis*(zheight+zTopCut)*cosphi, 01032 rRand1*ySemiAxis*(zheight+zTopCut)*sinphi, -zTopCut); 01033 01034 } 01035 // else 01036 // 01037 01038 do 01039 { 01040 rRand1 = RandFlat::shoot(0.,1.) ; 01041 rRand2 = RandFlat::shoot(0.,1.) ; 01042 } while ( rRand2 >= rRand1 ) ; 01043 01044 return G4ThreeVector(rRand1*xSemiAxis*(zheight-zTopCut)*cosphi, 01045 rRand1*ySemiAxis*(zheight-zTopCut)*sinphi, zTopCut); 01046 }
G4Polyhedron * G4EllipticalCone::GetPolyhedron | ( | ) | const [virtual] |
Reimplemented from G4VSolid.
Definition at line 1082 of file G4EllipticalCone.cc.
References CreatePolyhedron(), fpPolyhedron, and G4Polyhedron::GetNumberOfRotationStepsAtTimeOfCreation().
01083 { 01084 if ( (!fpPolyhedron) 01085 || (fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() != 01086 fpPolyhedron->GetNumberOfRotationSteps()) ) 01087 { 01088 delete fpPolyhedron; 01089 fpPolyhedron = CreatePolyhedron(); 01090 } 01091 return fpPolyhedron; 01092 }
G4double G4EllipticalCone::GetSemiAxisMax | ( | ) | const [inline] |
G4double G4EllipticalCone::GetSemiAxisX | ( | ) | const [inline] |
G4double G4EllipticalCone::GetSemiAxisY | ( | ) | const [inline] |
G4double G4EllipticalCone::GetSurfaceArea | ( | ) | [inline, virtual] |
Reimplemented from G4VSolid.
Definition at line 112 of file G4EllipticalCone.icc.
References G4VSolid::GetSurfaceArea().
00113 { 00114 if(fSurfaceArea != 0.) {;} 00115 else { fSurfaceArea = G4VSolid::GetSurfaceArea(); } 00116 return fSurfaceArea; 00117 }
G4double G4EllipticalCone::GetZMax | ( | ) | const [inline] |
G4double G4EllipticalCone::GetZTopCut | ( | ) | const [inline] |
EInside G4EllipticalCone::Inside | ( | const G4ThreeVector & | p | ) | const [virtual] |
Implements G4VSolid.
Definition at line 281 of file G4EllipticalCone.cc.
References G4VSolid::kCarTolerance, kInside, kOutside, kSurface, and sqr().
Referenced by DistanceToOut().
00282 { 00283 G4double rad2oo, // outside surface outer tolerance 00284 rad2oi; // outside surface inner tolerance 00285 00286 EInside in; 00287 00288 static const G4double halfRadTol = 0.5*kRadTolerance; 00289 static const G4double halfCarTol = 0.5*kCarTolerance; 00290 00291 // check this side of z cut first, because that's fast 00292 // 00293 00294 if ( (p.z() < -zTopCut - halfCarTol) 00295 || (p.z() > zTopCut + halfCarTol ) ) 00296 { 00297 return in = kOutside; 00298 } 00299 00300 rad2oo= sqr(p.x()/( xSemiAxis + halfRadTol )) 00301 + sqr(p.y()/( ySemiAxis + halfRadTol )); 00302 00303 if ( rad2oo > sqr( zheight-p.z() ) ) 00304 { 00305 return in = kOutside; 00306 } 00307 00308 // rad2oi= sqr( p.x()*(1.0 + 0.5*kRadTolerance/(xSemiAxis*xSemiAxis)) ) 00309 // + sqr( p.y()*(1.0 + 0.5*kRadTolerance/(ySemiAxis*ySemiAxis)) ); 00310 rad2oi = sqr(p.x()/( xSemiAxis - halfRadTol )) 00311 + sqr(p.y()/( ySemiAxis - halfRadTol )); 00312 00313 if (rad2oi < sqr( zheight-p.z() ) ) 00314 { 00315 in = ( ( p.z() < -zTopCut + halfRadTol ) 00316 || ( p.z() > zTopCut - halfRadTol ) ) ? kSurface : kInside; 00317 } 00318 else 00319 { 00320 in = kSurface; 00321 } 00322 00323 return in; 00324 }
G4EllipticalCone & G4EllipticalCone::operator= | ( | const G4EllipticalCone & | rhs | ) |
Definition at line 137 of file G4EllipticalCone.cc.
References fCubicVolume, fpPolyhedron, fSurfaceArea, kRadTolerance, G4VSolid::operator=(), semiAxisMax, xSemiAxis, ySemiAxis, zheight, and zTopCut.
00138 { 00139 // Check assignment to self 00140 // 00141 if (this == &rhs) { return *this; } 00142 00143 // Copy base class data 00144 // 00145 G4VSolid::operator=(rhs); 00146 00147 // Copy data 00148 // 00149 fpPolyhedron = 0; kRadTolerance = rhs.kRadTolerance; 00150 fCubicVolume = rhs.fCubicVolume; fSurfaceArea = rhs.fSurfaceArea; 00151 xSemiAxis = rhs.xSemiAxis; ySemiAxis = rhs.ySemiAxis; 00152 zheight = rhs.zheight; semiAxisMax = rhs.semiAxisMax; zTopCut = rhs.zTopCut; 00153 00154 return *this; 00155 }
Definition at line 71 of file G4EllipticalCone.icc.
Referenced by G4EllipticalCone().
00074 { 00075 xSemiAxis = newxSemiAxis; 00076 ySemiAxis = newySemiAxis; 00077 zheight = newzMax; 00078 semiAxisMax = xSemiAxis > ySemiAxis ? xSemiAxis : ySemiAxis; 00079 if (zTopCut > +zheight) { zTopCut = +zheight; } 00080 }
void G4EllipticalCone::SetZCut | ( | G4double | newzTopCut | ) | [inline] |
Definition at line 83 of file G4EllipticalCone.icc.
Referenced by G4EllipticalCone().
00084 { 00085 if (newzTopCut > +zheight) 00086 { zTopCut = +zheight; } 00087 else 00088 { zTopCut = newzTopCut; } 00089 }
std::ostream & G4EllipticalCone::StreamInfo | ( | std::ostream & | os | ) | const [virtual] |
Implements G4VSolid.
Definition at line 969 of file G4EllipticalCone.cc.
References G4VSolid::GetName().
00970 { 00971 G4int oldprc = os.precision(16); 00972 os << "-----------------------------------------------------------\n" 00973 << " *** Dump for solid - " << GetName() << " ***\n" 00974 << " ===================================================\n" 00975 << " Solid type: G4EllipticalCone\n" 00976 << " Parameters: \n" 00977 00978 << " semi-axis x: " << xSemiAxis/mm << " mm \n" 00979 << " semi-axis y: " << ySemiAxis/mm << " mm \n" 00980 << " height z: " << zheight/mm << " mm \n" 00981 << " half length in z: " << zTopCut/mm << " mm \n" 00982 << "-----------------------------------------------------------\n"; 00983 os.precision(oldprc); 00984 00985 return os; 00986 }
G4ThreeVector G4EllipticalCone::SurfaceNormal | ( | const G4ThreeVector & | p | ) | const [virtual] |
Implements G4VSolid.
Definition at line 330 of file G4EllipticalCone.cc.
References G4INCL::Math::pi, and sqr().
00331 { 00332 00333 G4double rx = sqr(p.x()/xSemiAxis), 00334 ry = sqr(p.y()/ySemiAxis); 00335 00336 G4double rds = std::sqrt(rx + ry); 00337 00338 G4ThreeVector norm; 00339 00340 if( (p.z() < -zTopCut) && ((rx+ry) < sqr(zTopCut + zheight)) ) 00341 { 00342 return G4ThreeVector( 0., 0., -1. ); 00343 } 00344 00345 if( (p.z() > (zheight > zTopCut ? zheight : zTopCut)) && 00346 ((rx+ry) < sqr(zheight-zTopCut)) ) 00347 { 00348 return G4ThreeVector( 0., 0., 1. ); 00349 } 00350 00351 if( p.z() > rds + 2.*zTopCut - zheight ) 00352 { 00353 if ( p.z() > zTopCut ) 00354 { 00355 if( p.x() == 0. ) 00356 { 00357 norm = G4ThreeVector( 0., p.y() < 0. ? -1. : 1., 1. ); 00358 return norm /= norm.mag(); 00359 } 00360 if( p.y() == 0. ) 00361 { 00362 norm = G4ThreeVector( p.x() < 0. ? -1. : 1., 0., 1. ); 00363 return norm /= norm.mag(); 00364 } 00365 00366 G4double k = std::fabs(p.x()/p.y()); 00367 G4double c2 = sqr(zheight-zTopCut)/(1./sqr(xSemiAxis)+sqr(k/ySemiAxis)); 00368 G4double x = std::sqrt(c2); 00369 G4double y = k*x; 00370 00371 x /= sqr(xSemiAxis); 00372 y /= sqr(ySemiAxis); 00373 00374 norm = G4ThreeVector( p.x() < 0. ? -x : x, 00375 p.y() < 0. ? -y : y, 00376 - ( zheight - zTopCut ) ); 00377 norm /= norm.mag(); 00378 norm += G4ThreeVector( 0., 0., 1. ); 00379 return norm /= norm.mag(); 00380 } 00381 00382 return G4ThreeVector( 0., 0., 1. ); 00383 } 00384 00385 if( p.z() < rds - 2.*zTopCut - zheight ) 00386 { 00387 if( p.x() == 0. ) 00388 { 00389 norm = G4ThreeVector( 0., p.y() < 0. ? -1. : 1., -1. ); 00390 return norm /= norm.mag(); 00391 } 00392 if( p.y() == 0. ) 00393 { 00394 norm = G4ThreeVector( p.x() < 0. ? -1. : 1., 0., -1. ); 00395 return norm /= norm.mag(); 00396 } 00397 00398 G4double k = std::fabs(p.x()/p.y()); 00399 G4double c2 = sqr(zheight+zTopCut)/(1./sqr(xSemiAxis)+sqr(k/ySemiAxis)); 00400 G4double x = std::sqrt(c2); 00401 G4double y = k*x; 00402 00403 x /= sqr(xSemiAxis); 00404 y /= sqr(ySemiAxis); 00405 00406 norm = G4ThreeVector( p.x() < 0. ? -x : x, 00407 p.y() < 0. ? -y : y, 00408 - ( zheight - zTopCut ) ); 00409 norm /= norm.mag(); 00410 norm += G4ThreeVector( 0., 0., -1. ); 00411 return norm /= norm.mag(); 00412 } 00413 00414 norm = G4ThreeVector(p.x()/sqr(xSemiAxis), p.y()/sqr(ySemiAxis), rds); 00415 00416 G4double k = std::tan(pi/8.); 00417 G4double c = -zTopCut - k*(zTopCut + zheight); 00418 00419 if( p.z() < -k*rds + c ) 00420 return G4ThreeVector (0.,0.,-1.); 00421 00422 return norm /= norm.mag(); 00423 }
G4Polyhedron* G4EllipticalCone::fpPolyhedron [mutable, protected] |