G4EllipticalCone.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 // $Id: G4EllipticalCone.cc 67011 2013-01-29 16:17:41Z gcosmo $
00027 //
00028 // Implementation of G4EllipticalCone class
00029 //
00030 // This code implements an Elliptical Cone given explicitly by the
00031 // equation:
00032 //   x^2/a^2 + y^2/b^2 = (z-h)^2
00033 // and specified by the parameters (a,b,h) and a cut parallel to the
00034 // xy plane above z = 0.
00035 //
00036 // Author: Dionysios Anninos
00037 //
00038 // --------------------------------------------------------------------
00039 
00040 #include "globals.hh"
00041 
00042 #include "G4EllipticalCone.hh"
00043 
00044 #include "G4ClippablePolygon.hh"
00045 #include "G4SolidExtentList.hh"
00046 #include "G4VoxelLimits.hh"
00047 #include "G4AffineTransform.hh"
00048 #include "G4GeometryTolerance.hh"
00049 
00050 #include "meshdefs.hh"
00051 
00052 #include "Randomize.hh"
00053 
00054 #include "G4VGraphicsScene.hh"
00055 #include "G4Polyhedron.hh"
00056 #include "G4NURBS.hh"
00057 #include "G4NURBSbox.hh"
00058 #include "G4VisExtent.hh"
00059 
00060 //#define G4SPECSDEBUG 1    
00061 
00062 using namespace CLHEP;
00063 
00065 //
00066 // Constructor - check parameters
00067 //
00068 G4EllipticalCone::G4EllipticalCone(const G4String& pName,
00069                                          G4double  pxSemiAxis,
00070                                          G4double  pySemiAxis,
00071                                          G4double  pzMax,
00072                                          G4double  pzTopCut)
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 }
00099 
00101 //
00102 // Fake default constructor - sets only member data and allocates memory
00103 //                            for usage restricted to object persistency.
00104 //
00105 G4EllipticalCone::G4EllipticalCone( __void__& a )
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 }
00111 
00113 //
00114 // Destructor
00115 //
00116 G4EllipticalCone::~G4EllipticalCone()
00117 {
00118 }
00119 
00121 //
00122 // Copy constructor
00123 //
00124 G4EllipticalCone::G4EllipticalCone(const G4EllipticalCone& rhs)
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 }
00132 
00134 //
00135 // Assignment operator
00136 //
00137 G4EllipticalCone& G4EllipticalCone::operator = (const G4EllipticalCone& rhs) 
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 }
00156 
00158 //
00159 // Calculate extent under transform and specified limit
00160 //
00161 G4bool
00162 G4EllipticalCone::CalculateExtent( const EAxis axis,
00163                                    const G4VoxelLimits &voxelLimit,
00164                                    const G4AffineTransform &transform,
00165                                          G4double &min, G4double &max ) const
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 }
00274 
00276 //
00277 // Return whether point inside/outside/on surface
00278 // Split into radius, phi, theta checks
00279 // Each check modifies `in', or returns as approprate
00280 //
00281 EInside G4EllipticalCone::Inside(const G4ThreeVector& p) const
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 }
00325 
00327 //
00328 // Return unit normal of surface closest to p not protected against p=0
00329 //
00330 G4ThreeVector G4EllipticalCone::SurfaceNormal( const G4ThreeVector& p) const
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 }
00424 
00426 //
00427 // Calculate distance to shape from outside, along normalised vector
00428 // return kInfinity if no intersection, or intersection distance <= tolerance
00429 //
00430 G4double G4EllipticalCone::DistanceToIn( const G4ThreeVector& p,
00431                                          const G4ThreeVector& v  ) const
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 }
00675 
00677 //
00678 // Calculate distance (<= actual) to closest surface of shape from outside
00679 // Return 0 if point inside
00680 //
00681 G4double G4EllipticalCone::DistanceToIn(const G4ThreeVector& p) const
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 }
00733 
00735 //
00736 // Calculate distance to surface of shape from `inside',
00737 // allowing for tolerance
00738 //
00739 G4double G4EllipticalCone::DistanceToOut(const G4ThreeVector& p,
00740                                          const G4ThreeVector& v,
00741                                          const G4bool calcNorm,
00742                                                G4bool *validNorm,
00743                                                G4ThreeVector *n  ) const
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 }
00896 
00898 //
00899 // Calculate distance (<=actual) to closest surface of shape from inside
00900 //
00901 G4double G4EllipticalCone::DistanceToOut(const G4ThreeVector& p) const
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 }
00946 
00948 //
00949 // GetEntityType
00950 //
00951 G4GeometryType G4EllipticalCone::GetEntityType() const
00952 {
00953   return G4String("G4EllipticalCone");
00954 }
00955 
00957 //
00958 // Make a clone of the object
00959 //
00960 G4VSolid* G4EllipticalCone::Clone() const
00961 {
00962   return new G4EllipticalCone(*this);
00963 }
00964 
00966 //
00967 // Stream object contents to an output stream
00968 //
00969 std::ostream& G4EllipticalCone::StreamInfo( std::ostream& os ) const
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 }
00987 
00989 //
00990 // GetPointOnSurface
00991 //
00992 // returns quasi-uniformly distributed point on surface of elliptical cone
00993 //
00994 G4ThreeVector G4EllipticalCone::GetPointOnSurface() const
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 }
01047 
01048 //
01049 // Methods for visualisation
01050 //
01051 
01052 void G4EllipticalCone::DescribeYourselfTo (G4VGraphicsScene& scene) const
01053 {
01054   scene.AddSolid(*this);
01055 }
01056 
01057 G4VisExtent G4EllipticalCone::GetExtent() const
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 }
01069 
01070 G4NURBS* G4EllipticalCone::CreateNURBS () const
01071 {
01072   // Box for now!!!
01073   //
01074   return new G4NURBSbox(xSemiAxis, ySemiAxis,zheight);
01075 }
01076 
01077 G4Polyhedron* G4EllipticalCone::CreatePolyhedron () const
01078 {
01079   return new G4PolyhedronEllipticalCone(xSemiAxis, ySemiAxis, zheight, zTopCut);
01080 }
01081 
01082 G4Polyhedron* G4EllipticalCone::GetPolyhedron () const
01083 {
01084   if ( (!fpPolyhedron)
01085     || (fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() !=
01086         fpPolyhedron->GetNumberOfRotationSteps()) )
01087     {
01088       delete fpPolyhedron;
01089       fpPolyhedron = CreatePolyhedron();
01090     }
01091   return fpPolyhedron;
01092 }

Generated on Mon May 27 17:48:07 2013 for Geant4 by  doxygen 1.4.7