G4Hype.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 //
00027 // $Id: G4Hype.cc 67011 2013-01-29 16:17:41Z gcosmo $
00028 // $Original: G4Hype.cc,v 1.0 1998/06/09 16:57:50 safai Exp $
00029 //
00030 // 
00031 // --------------------------------------------------------------------
00032 // GEANT 4 class source file
00033 //
00034 //
00035 // G4Hype.cc
00036 //
00037 // --------------------------------------------------------------------
00038 //
00039 // Authors: 
00040 //      Ernesto Lamanna (Ernesto.Lamanna@roma1.infn.it) &
00041 //      Francesco Safai Tehrani (Francesco.SafaiTehrani@roma1.infn.it)
00042 //      Rome, INFN & University of Rome "La Sapienza",  9 June 1998.
00043 //
00044 // --------------------------------------------------------------------
00045 
00046 #include "G4Hype.hh"
00047 
00048 #include "G4VoxelLimits.hh"
00049 #include "G4AffineTransform.hh"
00050 #include "G4SolidExtentList.hh"
00051 #include "G4ClippablePolygon.hh"
00052 
00053 #include "G4VPVParameterisation.hh"
00054 
00055 #include "meshdefs.hh"
00056 
00057 #include <cmath>
00058 
00059 #include "Randomize.hh"
00060 
00061 #include "G4VGraphicsScene.hh"
00062 #include "G4Polyhedron.hh"
00063 #include "G4VisExtent.hh"
00064 #include "G4NURBS.hh"
00065 #include "G4NURBStube.hh"
00066 #include "G4NURBScylinder.hh"
00067 #include "G4NURBStubesector.hh"
00068 
00069 using namespace CLHEP;
00070 
00071 // Constructor - check parameters, and fills protected data members
00072 G4Hype::G4Hype(const G4String& pName,
00073                      G4double newInnerRadius,
00074                      G4double newOuterRadius,
00075                      G4double newInnerStereo,
00076                      G4double newOuterStereo,
00077                      G4double newHalfLenZ)
00078   : G4VSolid(pName), fCubicVolume(0.), fSurfaceArea(0.), fpPolyhedron(0)
00079 {
00080   // Check z-len
00081   //
00082   if (newHalfLenZ<=0)
00083   {
00084     std::ostringstream message;
00085     message << "Invalid Z half-length - " << GetName() << G4endl
00086             << "        Invalid Z half-length: "
00087             << newHalfLenZ/mm << " mm";
00088     G4Exception("G4Hype::G4Hype()", "GeomSolids0002",
00089                 FatalErrorInArgument, message);
00090   }
00091   halfLenZ=newHalfLenZ;   
00092 
00093   // Check radii
00094   //
00095   if (newInnerRadius<0 || newOuterRadius<0)
00096   {
00097     std::ostringstream message;
00098     message << "Invalid radii - " << GetName() << G4endl
00099             << "        Invalid radii !  Inner radius: "
00100             << newInnerRadius/mm << " mm" << G4endl
00101             << "                         Outer radius: "
00102             << newOuterRadius/mm << " mm";
00103     G4Exception("G4Hype::G4Hype()", "GeomSolids0002",
00104                 FatalErrorInArgument, message);
00105   }
00106   if (newInnerRadius >= newOuterRadius)
00107   {
00108     std::ostringstream message;
00109     message << "Outer > inner radius - " << GetName() << G4endl
00110             << "        Invalid radii !  Inner radius: "
00111             << newInnerRadius/mm << " mm" << G4endl
00112             << "                         Outer radius: "
00113             << newOuterRadius/mm << " mm";
00114     G4Exception("G4Hype::G4Hype()", "GeomSolids0002",
00115                 FatalErrorInArgument, message);
00116   }
00117 
00118   innerRadius=newInnerRadius;
00119   outerRadius=newOuterRadius;
00120 
00121   innerRadius2=innerRadius*innerRadius;
00122   outerRadius2=outerRadius*outerRadius;
00123     
00124   SetInnerStereo( newInnerStereo );
00125   SetOuterStereo( newOuterStereo );
00126 }
00127 
00128 
00129 //
00130 // Fake default constructor - sets only member data and allocates memory
00131 //                            for usage restricted to object persistency.
00132 //
00133 G4Hype::G4Hype( __void__& a  )
00134   : G4VSolid(a), innerRadius(0.), outerRadius(0.), halfLenZ(0.), innerStereo(0.),
00135     outerStereo(0.), tanInnerStereo(0.), tanOuterStereo(0.), tanInnerStereo2(0.),
00136     tanOuterStereo2(0.), innerRadius2(0.), outerRadius2(0.), endInnerRadius2(0.),
00137     endOuterRadius2(0.), endInnerRadius(0.), endOuterRadius(0.),
00138     fCubicVolume(0.), fSurfaceArea(0.), fpPolyhedron(0)
00139 {
00140 }
00141 
00142 
00143 //
00144 // Destructor
00145 //
00146 G4Hype::~G4Hype()
00147 {
00148   delete fpPolyhedron;
00149 }
00150 
00151 
00152 //
00153 // Copy constructor
00154 //
00155 G4Hype::G4Hype(const G4Hype& rhs)
00156   : G4VSolid(rhs), innerRadius(rhs.innerRadius),
00157     outerRadius(rhs.outerRadius), halfLenZ(rhs.halfLenZ),
00158     innerStereo(rhs.innerStereo), outerStereo(rhs.outerStereo),
00159     tanInnerStereo(rhs.tanInnerStereo), tanOuterStereo(rhs.tanOuterStereo),
00160     tanInnerStereo2(rhs.tanInnerStereo2), tanOuterStereo2(rhs.tanOuterStereo2),
00161     innerRadius2(rhs.innerRadius2), outerRadius2(rhs.outerRadius2),
00162     endInnerRadius2(rhs.endInnerRadius2), endOuterRadius2(rhs.endOuterRadius2),
00163     endInnerRadius(rhs.endInnerRadius), endOuterRadius(rhs.endOuterRadius),
00164     fCubicVolume(rhs.fCubicVolume), fSurfaceArea(rhs.fSurfaceArea),
00165     fpPolyhedron(0)
00166 {
00167 }
00168 
00169 
00170 //
00171 // Assignment operator
00172 //
00173 G4Hype& G4Hype::operator = (const G4Hype& rhs) 
00174 {
00175    // Check assignment to self
00176    //
00177    if (this == &rhs)  { return *this; }
00178 
00179    // Copy base class data
00180    //
00181    G4VSolid::operator=(rhs);
00182 
00183    // Copy data
00184    //
00185    innerRadius = rhs.innerRadius; outerRadius = rhs.outerRadius;
00186    halfLenZ = rhs.halfLenZ;
00187    innerStereo = rhs.innerStereo; outerStereo = rhs.outerStereo;
00188    tanInnerStereo = rhs.tanInnerStereo; tanOuterStereo = rhs.tanOuterStereo;
00189    tanInnerStereo2 = rhs.tanInnerStereo2; tanOuterStereo2 = rhs.tanOuterStereo2;
00190    innerRadius2 = rhs.innerRadius2; outerRadius2 = rhs.outerRadius2;
00191    endInnerRadius2 = rhs.endInnerRadius2; endOuterRadius2 = rhs.endOuterRadius2;
00192    endInnerRadius = rhs.endInnerRadius; endOuterRadius = rhs.endOuterRadius;
00193    fCubicVolume = rhs.fCubicVolume; fSurfaceArea = rhs.fSurfaceArea;
00194    fpPolyhedron = 0; 
00195 
00196    return *this;
00197 }
00198 
00199 
00200 //
00201 // Dispatch to parameterisation for replication mechanism dimension
00202 // computation & modification.
00203 //
00204 void G4Hype::ComputeDimensions(G4VPVParameterisation* p,
00205                               const G4int n,
00206                               const G4VPhysicalVolume* pRep)
00207 {
00208   p->ComputeDimensions(*this,n,pRep);
00209 }
00210 
00211 
00212 //
00213 // CalculateExtent
00214 //
00215 G4bool G4Hype::CalculateExtent( const EAxis axis,
00216                                 const G4VoxelLimits &voxelLimit,
00217                                 const G4AffineTransform &transform,
00218                                 G4double &min, G4double &max ) const
00219 {
00220   G4SolidExtentList  extentList( axis, voxelLimit );
00221   
00222   //
00223   // Choose phi size of our segment(s) based on constants as
00224   // defined in meshdefs.hh
00225   //
00226   G4int numPhi = kMaxMeshSections;
00227   G4double sigPhi = twopi/numPhi;
00228   G4double rFudge = 1.0/std::cos(0.5*sigPhi);
00229   
00230   //
00231   // We work around in phi building polygons along the way.
00232   // As a reasonable compromise between accuracy and
00233   // complexity (=cpu time), the following facets are chosen:
00234   //
00235   //   1. If outerRadius/endOuterRadius > 0.95, approximate
00236   //      the outer surface as a cylinder, and use one
00237   //      rectangular polygon (0-1) to build its mesh.
00238   //
00239   //      Otherwise, use two trapazoidal polygons that 
00240   //      meet at z = 0 (0-4-1)
00241   //
00242   //   2. If there is no inner surface, then use one
00243   //      polygon for each entire endcap.  (0) and (1)
00244   //
00245   //      Otherwise, use a trapazoidal polygon for each
00246   //      phi segment of each endcap.    (0-2) and (1-3)
00247   //
00248   //   3. For the inner surface, if innerRadius/endInnerRadius > 0.95,
00249   //      approximate the inner surface as a cylinder of
00250   //      radius innerRadius and use one rectangular polygon
00251   //      to build each phi segment of its mesh.   (2-3)
00252   //
00253   //      Otherwise, use one rectangular polygon centered
00254   //      at z = 0 (5-6) and two connecting trapazoidal polygons
00255   //      for each phi segment (2-5) and (3-6).
00256   //
00257   
00258   G4bool splitOuter = (outerRadius/endOuterRadius < 0.95);
00259   G4bool splitInner = 0;
00260   if (InnerSurfaceExists())
00261   {
00262     splitInner = (innerRadius/endInnerRadius < 0.95);
00263   }
00264   
00265   //
00266   // Vertex assignments (v and w arrays)
00267   // [0] and [1] are mandatory
00268   // the rest are optional
00269   //
00270   //     +                     -
00271   //      [0]------[4]------[1]      <--- outer radius
00272   //       |                 |       
00273   //       |                 |       
00274   //      [2]---[5]---[6]---[3]      <--- inner radius
00275   //
00276 
00277 
00278   G4ClippablePolygon endPoly1, endPoly2;
00279   
00280   G4double phi = 0, 
00281      cosPhi = std::cos(phi),
00282      sinPhi = std::sin(phi);
00283   G4ThreeVector v0( rFudge*endOuterRadius*cosPhi,
00284                     rFudge*endOuterRadius*sinPhi,
00285                     +halfLenZ ),
00286                 v1( rFudge*endOuterRadius*cosPhi,
00287                     rFudge*endOuterRadius*sinPhi,
00288                     -halfLenZ ),
00289                 v2, v3, v4, v5, v6,
00290                 w0, w1, w2, w3, w4, w5, w6;
00291   transform.ApplyPointTransform( v0 );
00292   transform.ApplyPointTransform( v1 );
00293   
00294   G4double zInnerSplit=0.;
00295   if (InnerSurfaceExists())
00296   {
00297     if (splitInner)
00298     {
00299       v2 = transform.TransformPoint( 
00300         G4ThreeVector( endInnerRadius*cosPhi,
00301                        endInnerRadius*sinPhi, +halfLenZ ) );
00302       v3 = transform.TransformPoint( 
00303         G4ThreeVector( endInnerRadius*cosPhi,
00304                        endInnerRadius*sinPhi, -halfLenZ ) );
00305       //
00306       // Find intersection of line normal to inner
00307       // surface at z = halfLenZ and line r=innerRadius
00308       //
00309       G4double rn = halfLenZ*tanInnerStereo2;
00310       G4double zn = endInnerRadius;
00311 
00312       zInnerSplit = halfLenZ + (innerRadius - endInnerRadius)*zn/rn;
00313 
00314       //
00315       // Build associated vertices
00316       //      
00317       v5 = transform.TransformPoint( 
00318         G4ThreeVector( innerRadius*cosPhi,
00319                        innerRadius*sinPhi, +zInnerSplit ) );
00320       v6 = transform.TransformPoint( 
00321         G4ThreeVector( innerRadius*cosPhi,
00322                        innerRadius*sinPhi, -zInnerSplit ) );
00323     }
00324     else
00325     {
00326       v2 = transform.TransformPoint( 
00327         G4ThreeVector( innerRadius*cosPhi,
00328                        innerRadius*sinPhi, +halfLenZ ) );
00329       v3 = transform.TransformPoint( 
00330         G4ThreeVector( innerRadius*cosPhi,
00331                        innerRadius*sinPhi, -halfLenZ ) );
00332     }
00333   }
00334   
00335   if (splitOuter)
00336   {
00337     v4 = transform.TransformPoint( 
00338       G4ThreeVector( rFudge*outerRadius*cosPhi,
00339                      rFudge*outerRadius*sinPhi, 0 ) );
00340   }
00341   
00342   //
00343   // Loop over phi segments
00344   //
00345   do
00346   {
00347     phi += sigPhi;
00348     if (numPhi == 1) phi = 0;  // Try to avoid roundoff
00349           cosPhi = std::cos(phi), 
00350     sinPhi = std::sin(phi);
00351     
00352     G4double r(rFudge*endOuterRadius);
00353     w0 = G4ThreeVector( r*cosPhi, r*sinPhi, +halfLenZ );
00354     w1 = G4ThreeVector( r*cosPhi, r*sinPhi, -halfLenZ );
00355     transform.ApplyPointTransform( w0 );
00356     transform.ApplyPointTransform( w1 );
00357   
00358     //
00359     // Outer hyperbolic surface
00360     //
00361     if (splitOuter)
00362     {
00363       r = rFudge*outerRadius;
00364       w4 = G4ThreeVector( r*cosPhi, r*sinPhi, 0 );
00365       transform.ApplyPointTransform( w4 );
00366       
00367       AddPolyToExtent( v0, v4, w4, w0, voxelLimit, axis, extentList );
00368       AddPolyToExtent( v4, v1, w1, w4, voxelLimit, axis, extentList );
00369     }
00370     else
00371     {
00372       AddPolyToExtent( v0, v1, w1, w0, voxelLimit, axis, extentList );
00373     }
00374   
00375     if (InnerSurfaceExists())
00376     {
00377       //
00378       // Inner hyperbolic surface
00379       //
00380       if (splitInner)
00381       {
00382         w2 = G4ThreeVector( endInnerRadius*cosPhi,
00383                             endInnerRadius*sinPhi, +halfLenZ );
00384         w3 = G4ThreeVector( endInnerRadius*cosPhi,
00385                             endInnerRadius*sinPhi, -halfLenZ );
00386         transform.ApplyPointTransform( w2 );
00387         transform.ApplyPointTransform( w3 );
00388 
00389         w5 = G4ThreeVector( innerRadius*cosPhi,
00390                             innerRadius*sinPhi, +zInnerSplit );
00391         w6 = G4ThreeVector( innerRadius*cosPhi,
00392                             innerRadius*sinPhi, -zInnerSplit );
00393         transform.ApplyPointTransform( w5 );
00394         transform.ApplyPointTransform( w6 );
00395         AddPolyToExtent( v3, v6, w6, w3, voxelLimit, axis, extentList );
00396         AddPolyToExtent( v6, v5, w5, w6, voxelLimit, axis, extentList );
00397         AddPolyToExtent( v5, v2, w2, w5, voxelLimit, axis, extentList );
00398       }
00399       else
00400       {
00401         w2 = G4ThreeVector( innerRadius*cosPhi,
00402                             innerRadius*sinPhi, +halfLenZ );
00403         w3 = G4ThreeVector( innerRadius*cosPhi,
00404                             innerRadius*sinPhi, -halfLenZ );
00405         transform.ApplyPointTransform( w2 );
00406         transform.ApplyPointTransform( w3 );
00407 
00408         AddPolyToExtent( v3, v2, w2, w3, voxelLimit, axis, extentList );
00409       }
00410 
00411       //
00412       // Endplate segments
00413       //
00414       AddPolyToExtent( v1, v3, w3, w1, voxelLimit, axis, extentList );
00415       AddPolyToExtent( v2, v0, w0, w2, voxelLimit, axis, extentList );
00416     }
00417     else
00418     {
00419       //
00420       // Continue building endplate polygons
00421       //
00422       endPoly1.AddVertexInOrder( v0 );
00423       endPoly2.AddVertexInOrder( v1 );
00424     }
00425 
00426     //
00427     // Next phi segments
00428     //    
00429     v0 = w0;
00430     v1 = w1;
00431     if (InnerSurfaceExists())
00432     {
00433       v2 = w2;
00434       v3 = w3;
00435       if (splitInner)
00436       {
00437         v5 = w5;
00438         v6 = w6;
00439       }
00440     }
00441     if (splitOuter) v4 = w4;
00442     
00443   } while( --numPhi > 0 );
00444   
00445   
00446   //
00447   // Don't forget about the endplate polygons, if
00448   // we use them
00449   //
00450   if (!InnerSurfaceExists())
00451   {
00452     if (endPoly1.PartialClip( voxelLimit, axis ))
00453     {
00454       static const G4ThreeVector normal(0,0,+1);
00455       endPoly1.SetNormal( transform.TransformAxis(normal) );
00456       extentList.AddSurface( endPoly1 );
00457     }
00458 
00459     if (endPoly2.PartialClip( voxelLimit, axis ))
00460     {
00461       static const G4ThreeVector normal(0,0,-1);
00462       endPoly2.SetNormal( transform.TransformAxis(normal) );
00463       extentList.AddSurface( endPoly2 );
00464     }
00465   }
00466   
00467   //
00468   // Return min/max value
00469   //
00470   return extentList.GetExtent( min, max );
00471 }
00472 
00473 
00474 //
00475 // AddPolyToExtent (static)
00476 //
00477 // Utility function for CalculateExtent
00478 //
00479 void G4Hype::AddPolyToExtent( const G4ThreeVector &v0,
00480                               const G4ThreeVector &v1,
00481                               const G4ThreeVector &w1,
00482                               const G4ThreeVector &w0,
00483                               const G4VoxelLimits &voxelLimit,
00484                               const EAxis axis,
00485                               G4SolidExtentList &extentList ) 
00486 {
00487   G4ClippablePolygon phiPoly;
00488 
00489   phiPoly.AddVertexInOrder( v0 );
00490   phiPoly.AddVertexInOrder( v1 );
00491   phiPoly.AddVertexInOrder( w1 );
00492   phiPoly.AddVertexInOrder( w0 );
00493 
00494   if (phiPoly.PartialClip( voxelLimit, axis ))
00495   {
00496     phiPoly.SetNormal( (v1-v0).cross(w0-v0).unit() );
00497     extentList.AddSurface( phiPoly );
00498   }
00499 }
00500 
00501 
00502 //
00503 // Decides whether point is inside,outside or on the surface
00504 //
00505 EInside G4Hype::Inside(const G4ThreeVector& p) const
00506 {
00507   static const G4double halfTol = 0.5*kCarTolerance;
00508   
00509   //
00510   // Check z extents: are we outside?
00511   //
00512   const G4double absZ(std::fabs(p.z()));
00513   if (absZ > halfLenZ + halfTol) return kOutside;
00514   
00515   //
00516   // Check outer radius
00517   //
00518   const G4double oRad2(HypeOuterRadius2(absZ));
00519   const G4double xR2( p.x()*p.x()+p.y()*p.y() );
00520   
00521   if (xR2 > oRad2 + kCarTolerance*endOuterRadius) return kOutside;
00522   
00523   if (xR2 > oRad2 - kCarTolerance*endOuterRadius) return kSurface;
00524   
00525   if (InnerSurfaceExists())
00526   {
00527     //
00528     // Check inner radius
00529     //
00530     const G4double iRad2(HypeInnerRadius2(absZ));
00531     
00532     if (xR2 < iRad2 - kCarTolerance*endInnerRadius) return kOutside;
00533     
00534     if (xR2 < iRad2 + kCarTolerance*endInnerRadius) return kSurface;
00535   }
00536   
00537   //
00538   // We are inside in radius, now check endplate surface
00539   //
00540   if (absZ > halfLenZ - halfTol) return kSurface;
00541   
00542   return kInside;
00543 }
00544 
00545 
00546 
00547 //
00548 // return the normal unit vector to the Hyperbolical Surface at a point 
00549 // p on (or nearly on) the surface
00550 //
00551 G4ThreeVector G4Hype::SurfaceNormal( const G4ThreeVector& p ) const
00552 {
00553   //
00554   // Which of the three or four surfaces are we closest to?
00555   //
00556   const G4double absZ(std::fabs(p.z()));
00557   const G4double distZ(absZ - halfLenZ);
00558   const G4double dist2Z(distZ*distZ);
00559   
00560   const G4double xR2( p.x()*p.x()+p.y()*p.y() );
00561   const G4double dist2Outer( std::fabs(xR2 - HypeOuterRadius2(absZ)) );
00562   
00563   if (InnerSurfaceExists())
00564   {
00565     //
00566     // Has inner surface: is this closest?
00567     //
00568     const G4double dist2Inner( std::fabs(xR2 - HypeInnerRadius2(absZ)) );
00569     if (dist2Inner < dist2Z && dist2Inner < dist2Outer)
00570       return G4ThreeVector( -p.x(), -p.y(), p.z()*tanInnerStereo2 ).unit();
00571   }
00572 
00573   //
00574   // Do the "endcaps" win?
00575   //
00576   if (dist2Z < dist2Outer) 
00577     return G4ThreeVector( 0.0, 0.0, p.z() < 0 ? -1.0 : 1.0 );
00578     
00579     
00580   //
00581   // Outer surface wins
00582   //
00583   return G4ThreeVector( p.x(), p.y(), -p.z()*tanOuterStereo2 ).unit();
00584 }
00585 
00586 
00587 //
00588 // Calculate distance to shape from outside, along normalised vector
00589 // - return kInfinity if no intersection,
00590 //   or intersection distance <= tolerance
00591 //
00592 // Calculating the intersection of a line with the surfaces
00593 // is fairly straight forward. The difficult problem is dealing
00594 // with the intersections of the surfaces in a consistent manner, 
00595 // and this accounts for the complicated logic.
00596 //
00597 G4double G4Hype::DistanceToIn( const G4ThreeVector& p,
00598                                const G4ThreeVector& v ) const
00599 {
00600   static const G4double halfTol = 0.5*kCarTolerance;
00601   
00602   //
00603   // Quick test. Beware! This assumes v is a unit vector!
00604   //
00605   if (std::fabs(p.x()*v.y() - p.y()*v.x()) > endOuterRadius+kCarTolerance)
00606     return kInfinity;
00607   
00608   //
00609   // Take advantage of z symmetry, and reflect throught the
00610   // z=0 plane so that pz is always positive
00611   //
00612   G4double pz(p.z()), vz(v.z());
00613   if (pz < 0)
00614   {
00615     pz = -pz;
00616     vz = -vz;
00617   }
00618 
00619   //
00620   // We must be very careful if we don't want to
00621   // create subtle leaks at the edges where the
00622   // hyperbolic surfaces connect to the endplate.
00623   // The only reliable way to do so is to make sure
00624   // that the decision as to when a track passes
00625   // over the edge of one surface is exactly the
00626   // same decision as to when a track passes into the
00627   // other surface. By "exact", we don't mean algebraicly
00628   // exact, but we mean the same machine instructions
00629   // should be used.
00630   //
00631   G4bool couldMissOuter(true),
00632          couldMissInner(true),
00633          cantMissInnerCylinder(false);
00634   
00635   //
00636   // Check endplate intersection
00637   //
00638   G4double sigz = pz-halfLenZ;
00639   
00640   if (sigz > -halfTol)    // equivalent to: if (pz > halfLenZ - halfTol)
00641   {
00642     //
00643     // We start in front of the endplate (within roundoff)
00644     // Correct direction to intersect endplate?
00645     //
00646     if (vz >= 0)
00647     {
00648       //
00649       // Nope. As long as we are far enough away, we
00650       // can't intersect anything
00651       //
00652       if (sigz > 0) return kInfinity;
00653       
00654       //
00655       // Otherwise, we may still hit a hyperbolic surface
00656       // if the point is on the hyperbolic surface (within tolerance)
00657       //
00658       G4double pr2 = p.x()*p.x() + p.y()*p.y();
00659       if (pr2 > endOuterRadius2 + kCarTolerance*endOuterRadius)
00660         return kInfinity;
00661       
00662       if (InnerSurfaceExists())
00663       {
00664         if (pr2 < endInnerRadius2 - kCarTolerance*endInnerRadius)
00665           return kInfinity;
00666         if ( (pr2 < endOuterRadius2 - kCarTolerance*endOuterRadius)
00667           && (pr2 > endInnerRadius2 + kCarTolerance*endInnerRadius) )
00668           return kInfinity;
00669       }
00670       else
00671       {
00672         if (pr2 < endOuterRadius2 - kCarTolerance*endOuterRadius)
00673           return kInfinity;
00674       }
00675     }
00676     else
00677     {
00678       //
00679       // Where do we intersect at z = halfLenZ?
00680       //
00681       G4double q = -sigz/vz;
00682       G4double xi = p.x() + q*v.x(),
00683                yi = p.y() + q*v.y();
00684          
00685       //
00686       // Is this on the endplate? If so, return s, unless
00687       // we are on the tolerant surface, in which case return 0
00688       //
00689       G4double pr2 = xi*xi + yi*yi;
00690       if (pr2 <= endOuterRadius2)
00691       {
00692         if (InnerSurfaceExists())
00693         {
00694           if (pr2 >= endInnerRadius2) return (sigz < halfTol) ? 0 : q;
00695           //
00696           // This test is sufficient to ensure that the
00697           // trajectory cannot miss the inner hyperbolic surface
00698           // for z > 0, if the normal is correct.
00699           //
00700           G4double dot1 = (xi*v.x() + yi*v.y())*endInnerRadius/std::sqrt(pr2);
00701           couldMissInner = (dot1 - halfLenZ*tanInnerStereo2*vz <= 0);
00702           
00703           if (pr2 > endInnerRadius2*(1 - 2*DBL_EPSILON) )
00704           {
00705             //
00706             // There is a potential leak if the inner
00707             // surface is a cylinder
00708             //
00709             if ( (innerStereo < DBL_MIN)
00710               && ((std::fabs(v.x()) > DBL_MIN) || (std::fabs(v.y()) > DBL_MIN)) )
00711               cantMissInnerCylinder = true;
00712           }
00713         }
00714         else
00715         {
00716           return (sigz < halfTol) ? 0 : q;
00717         }
00718       }
00719       else
00720       {
00721         G4double dotR( xi*v.x() + yi*v.y() );
00722         if (dotR >= 0)
00723         {
00724           //
00725           // Otherwise, if we are traveling outwards, we know
00726           // we must miss the hyperbolic surfaces also, so
00727           // we need not bother checking
00728           //
00729           return kInfinity;
00730         }
00731         else
00732         {
00733           //
00734           // This test is sufficient to ensure that the
00735           // trajectory cannot miss the outer hyperbolic surface
00736           // for z > 0, if the normal is correct.
00737           //
00738           G4double dot1 = dotR*endOuterRadius/std::sqrt(pr2);
00739           couldMissOuter = (dot1 - halfLenZ*tanOuterStereo2*vz>= 0);
00740         }
00741       }
00742     }
00743   }
00744     
00745   //
00746   // Check intersection with outer hyperbolic surface, save
00747   // distance to valid intersection into "best".
00748   //    
00749   G4double best = kInfinity;
00750   
00751   G4double q[2];
00752   G4int n = IntersectHype( p, v, outerRadius2, tanOuterStereo2, q );
00753   
00754   if (n > 0)
00755   {
00756     //
00757     // Potential intersection: is p on this surface?
00758     //
00759     if (pz < halfLenZ+halfTol)
00760     {
00761       G4double dr2 = p.x()*p.x() + p.y()*p.y() - HypeOuterRadius2(pz);
00762       if (std::fabs(dr2) < kCarTolerance*endOuterRadius)
00763       {
00764         //
00765         // Sure, but make sure we're traveling inwards at
00766         // this point
00767         //
00768         if (p.x()*v.x() + p.y()*v.y() - pz*tanOuterStereo2*vz < 0)
00769           return 0;
00770       }
00771     }
00772     
00773     //
00774     // We are now certain that p is not on the tolerant surface.
00775     // Accept only position distance q
00776     //
00777     G4int i;
00778     for( i=0; i<n; i++ )
00779     {
00780       if (q[i] >= 0)
00781       {
00782         //
00783         // Check to make sure this intersection point is
00784         // on the surface, but only do so if we haven't
00785         // checked the endplate intersection already
00786         //
00787         G4double zi = pz + q[i]*vz;
00788         
00789         if (zi < -halfLenZ) continue;
00790         if (zi > +halfLenZ && couldMissOuter) continue;
00791         
00792         //
00793         // Check normal
00794         //
00795         G4double xi = p.x() + q[i]*v.x(),
00796            yi = p.y() + q[i]*v.y();
00797            
00798         if (xi*v.x() + yi*v.y() - zi*tanOuterStereo2*vz > 0) continue;
00799 
00800         best = q[i];
00801         break;
00802       }
00803     }
00804   }
00805   
00806   if (!InnerSurfaceExists()) return best;    
00807   
00808   //
00809   // Check intersection with inner hyperbolic surface
00810   //
00811   n = IntersectHype( p, v, innerRadius2, tanInnerStereo2, q );  
00812   if (n == 0)
00813   {
00814     if (cantMissInnerCylinder) return (sigz < halfTol) ? 0 : -sigz/vz;
00815         
00816     return best;
00817   }
00818   
00819   //
00820   // P on this surface?
00821   //
00822   if (pz < halfLenZ+halfTol)
00823   {
00824     G4double dr2 = p.x()*p.x() + p.y()*p.y() - HypeInnerRadius2(pz);
00825     if (std::fabs(dr2) < kCarTolerance*endInnerRadius)
00826     {
00827       //
00828       // Sure, but make sure we're traveling outwards at
00829       // this point
00830       //
00831       if (p.x()*v.x() + p.y()*v.y() - pz*tanInnerStereo2*vz > 0) return 0;
00832     }
00833   }
00834   
00835   //
00836   // No, so only positive q is valid. Search for a valid intersection
00837   // that is closer than the outer intersection (if it exists)
00838   //
00839   G4int i;
00840   for( i=0; i<n; i++ )
00841   {
00842     if (q[i] > best) break;
00843     if (q[i] >= 0)
00844     {
00845       //
00846       // Check to make sure this intersection point is
00847       // on the surface, but only do so if we haven't
00848       // checked the endplate intersection already
00849       //
00850       G4double zi = pz + q[i]*vz;
00851 
00852       if (zi < -halfLenZ) continue;
00853       if (zi > +halfLenZ && couldMissInner) continue;
00854 
00855       //
00856       // Check normal
00857       //
00858       G4double xi = p.x() + q[i]*v.x(),
00859          yi = p.y() + q[i]*v.y();
00860 
00861       if (xi*v.x() + yi*v.y() - zi*tanOuterStereo2*vz < 0) continue;
00862 
00863       best = q[i];
00864       break;
00865     }
00866   }
00867     
00868   //
00869   // Done
00870   //
00871   return best;
00872 }
00873  
00874 
00875 //
00876 // Calculate distance to shape from outside, along perpendicular direction 
00877 // (if one exists). May be an underestimate.
00878 //
00879 // There are five (r,z) regions:
00880 //    1. a point that is beyond the endcap but within the
00881 //       endcap radii
00882 //    2. a point with r > outer endcap radius and with
00883 //       a z position that is beyond the cone formed by the
00884 //       normal of the outer hyperbolic surface at the 
00885 //       edge at which it meets the endcap. 
00886 //    3. a point that is outside the outer surface and not in (1 or 2)
00887 //    4. a point that is inside the inner surface and not in (5)
00888 //    5. a point with radius < inner endcap radius and
00889 //       with a z position beyond the cone formed by the
00890 //       normal of the inner hyperbolic surface at the
00891 //       edge at which it meets the endcap.
00892 // (regions 4 and 5 only exist if there is an inner surface)
00893 //
00894 G4double G4Hype::DistanceToIn(const G4ThreeVector& p) const
00895 {
00896   static const G4double halfTol(0.5*kCarTolerance);
00897   
00898   G4double absZ(std::fabs(p.z()));
00899   
00900   //
00901   // Check region
00902   //
00903   G4double r2 = p.x()*p.x() + p.y()*p.y();
00904   G4double r = std::sqrt(r2);
00905   
00906   G4double sigz = absZ - halfLenZ;
00907   
00908   if (r < endOuterRadius)
00909   {
00910     if (sigz > -halfTol)
00911     {
00912       if (InnerSurfaceExists())
00913       {
00914         if (r > endInnerRadius) 
00915           return sigz < halfTol ? 0 : sigz;  // Region 1
00916         
00917         G4double dr = endInnerRadius - r;
00918         if (sigz > dr*tanInnerStereo2)
00919         {
00920           //
00921           // In region 5
00922           //
00923           G4double answer = std::sqrt( dr*dr + sigz*sigz );
00924           return answer < halfTol ? 0 : answer;
00925         }
00926       }
00927       else
00928       {
00929         //
00930         // In region 1 (no inner surface)
00931         //
00932         return sigz < halfTol ? 0 : sigz;
00933       }
00934     }
00935   }
00936   else
00937   {
00938     G4double dr = r - endOuterRadius;
00939     if (sigz > -dr*tanOuterStereo2)
00940     {
00941       //
00942       // In region 2
00943       //
00944       G4double answer = std::sqrt( dr*dr + sigz*sigz );
00945       return answer < halfTol ? 0 : answer;
00946     }
00947   }
00948   
00949   if (InnerSurfaceExists())
00950   {
00951     if (r2 < HypeInnerRadius2(absZ)+kCarTolerance*endInnerRadius)
00952     {
00953        //
00954       // In region 4
00955       //
00956       G4double answer = ApproxDistInside( r,absZ,innerRadius,tanInnerStereo2 );
00957       return answer < halfTol ? 0 : answer;
00958     }
00959   }
00960   
00961   //
00962   // We are left by elimination with region 3
00963   //
00964   G4double answer = ApproxDistOutside( r, absZ, outerRadius, tanOuterStereo );
00965   return answer < halfTol ? 0 : answer;
00966 }
00967 
00968 
00969 //
00970 // Calculate distance to surface of shape from `inside', allowing for tolerance
00971 //
00972 // The situation here is much simplier than DistanceToIn(p,v). For
00973 // example, there is no need to even check whether an intersection
00974 // point is inside the boundary of a surface, as long as all surfaces 
00975 // are checked and the smallest distance is used.
00976 //
00977 G4double G4Hype::DistanceToOut( const G4ThreeVector& p, const G4ThreeVector& v,
00978                                 const G4bool calcNorm,
00979                                 G4bool *validNorm, G4ThreeVector *norm ) const
00980 {
00981   static const G4double halfTol = 0.5*kCarTolerance;
00982   
00983   
00984   static const G4ThreeVector normEnd1(0.0,0.0,+1.0);
00985   static const G4ThreeVector normEnd2(0.0,0.0,-1.0);
00986   
00987   //
00988   // Keep track of closest surface
00989   //
00990   G4double sBest;        // distance to
00991   const G4ThreeVector *nBest;    // normal vector
00992   G4bool vBest;        // whether "valid"
00993 
00994   //
00995   // Check endplate, taking advantage of symmetry.
00996   // Note that the endcap is the only surface which
00997   // has a "valid" normal, i.e. is a surface of which
00998   // the entire solid is behind.
00999   //
01000   G4double pz(p.z()), vz(v.z());
01001   if (vz < 0)
01002   {
01003     pz = -pz;
01004     vz = -vz;
01005     nBest = &normEnd2;
01006   }
01007   else
01008     nBest = &normEnd1;
01009 
01010   //
01011   // Possible intercept. Are we on the surface?
01012   //
01013   if (pz > halfLenZ-halfTol)
01014   {
01015     if (calcNorm) { *norm = *nBest; *validNorm = true; }
01016     return 0;
01017   }
01018 
01019   //
01020   // Nope. Get distance. Beware of zero vz.
01021   //
01022   sBest = (vz > DBL_MIN) ? (halfLenZ - pz)/vz : kInfinity;
01023   vBest = true;
01024   
01025   //
01026   // Check outer surface
01027   //
01028   G4double r2 = p.x()*p.x() + p.y()*p.y();
01029   
01030   G4double q[2];
01031   G4int n = IntersectHype( p, v, outerRadius2, tanOuterStereo2, q );
01032   
01033   G4ThreeVector norm1, norm2;
01034 
01035   if (n > 0)
01036   {
01037     //
01038     // We hit somewhere. Are we on the surface?
01039     //  
01040     G4double dr2 = r2 - HypeOuterRadius2(pz);
01041     if (std::fabs(dr2) < endOuterRadius*kCarTolerance)
01042     {
01043       G4ThreeVector normHere( p.x(), p.y(), -p.z()*tanOuterStereo2 );
01044       //
01045       // Sure. But are we going the right way?
01046       //
01047       if (normHere.dot(v) > 0)
01048       {
01049         if (calcNorm) { *norm = normHere.unit(); *validNorm = false; }
01050         return 0;
01051       }
01052     }
01053     
01054     //
01055     // Nope. Check closest positive intercept.
01056     //
01057     G4int i;
01058     for( i=0; i<n; i++ )
01059     {
01060       if (q[i] > sBest) break;
01061       if (q[i] > 0)
01062       {
01063         //
01064         // Make sure normal is correct (that this
01065         // solution is an outgoing solution)
01066         //
01067         G4ThreeVector pk(p+q[i]*v);
01068         norm1 = G4ThreeVector( pk.x(), pk.y(), -pk.z()*tanOuterStereo2 );
01069         if (norm1.dot(v) > 0)
01070         {
01071           sBest = q[i];
01072           nBest = &norm1;
01073           vBest = false;
01074           break;
01075         }
01076       }
01077     }
01078   }
01079   
01080   if (InnerSurfaceExists())
01081   {
01082     //
01083     // Check inner surface
01084     //
01085     n = IntersectHype( p, v, innerRadius2, tanInnerStereo2, q );
01086     if (n > 0)
01087     {
01088       //
01089       // On surface?
01090       //
01091       G4double dr2 = r2 - HypeInnerRadius2(pz);
01092       if (std::fabs(dr2) < endInnerRadius*kCarTolerance)
01093       {
01094         G4ThreeVector normHere( -p.x(), -p.y(), p.z()*tanInnerStereo2 );
01095         if (normHere.dot(v) > 0)
01096         {
01097           if (calcNorm)
01098           {
01099             *norm = normHere.unit();
01100             *validNorm = false;
01101           }
01102           return 0;
01103         }
01104       }
01105       
01106       //
01107       // Check closest positive
01108       //
01109       G4int i;
01110       for( i=0; i<n; i++ )
01111       {
01112         if (q[i] > sBest) break;
01113         if (q[i] > 0)
01114         {
01115           G4ThreeVector pk(p+q[i]*v);
01116           norm2 = G4ThreeVector( -pk.x(), -pk.y(), pk.z()*tanInnerStereo2 );
01117           if (norm2.dot(v) > 0)
01118           {
01119             sBest = q[i];
01120             nBest = &norm2;
01121             vBest = false;
01122             break;
01123           }
01124         }
01125       }
01126     }
01127   }
01128   
01129   //
01130   // Done!
01131   //
01132   if (calcNorm)
01133   {
01134     *validNorm = vBest;
01135     
01136     if (nBest == &norm1 || nBest == &norm2) 
01137       *norm = nBest->unit();
01138     else
01139       *norm = *nBest;
01140   }
01141   
01142   return sBest;
01143 }
01144 
01145 
01146 //
01147 // Calculate distance (<=actual) to closest surface of shape from inside
01148 //
01149 // May be an underestimate
01150 //
01151 G4double G4Hype::DistanceToOut(const G4ThreeVector& p) const
01152 {
01153   //
01154   // Try each surface and remember the closest
01155   //
01156   G4double absZ(std::fabs(p.z()));
01157   G4double r(p.perp());
01158   
01159   G4double sBest = halfLenZ - absZ;
01160   
01161   G4double tryOuter = ApproxDistInside( r, absZ, outerRadius, tanOuterStereo2 );
01162   if (tryOuter < sBest)
01163     sBest = tryOuter;
01164   
01165   if (InnerSurfaceExists())
01166   {
01167     G4double tryInner = ApproxDistOutside( r,absZ,innerRadius,tanInnerStereo );
01168     if (tryInner < sBest) sBest = tryInner;
01169   }
01170   
01171   return sBest < 0.5*kCarTolerance ? 0 : sBest;
01172 }
01173 
01174 
01175 //
01176 // IntersectHype (static)
01177 //
01178 // Decide if and where a line intersects with a hyperbolic
01179 // surface (of infinite extent)
01180 //
01181 // Arguments:
01182 //     p       - (in) Point on trajectory
01183 //     v       - (in) Vector along trajectory
01184 //     r2      - (in) Square of radius at z = 0
01185 //     tan2phi - (in) std::tan(phi)**2
01186 //     q       - (out) Up to two points of intersection, where the
01187 //                     intersection point is p + q*v, and if there are
01188 //                     two intersections, q[0] < q[1]. May be negative.
01189 // Returns:
01190 //     The number of intersections. If 0, the trajectory misses. 
01191 //
01192 //
01193 // Equation of a line:
01194 //
01195 //       x = x0 + q*tx      y = y0 + q*ty      z = z0 + q*tz
01196 //
01197 // Equation of a hyperbolic surface:
01198 //
01199 //       x**2 + y**2 = r**2 + (z*tanPhi)**2
01200 //
01201 // Solution is quadratic:
01202 //
01203 //  a*q**2 + b*q + c = 0
01204 //
01205 // where:
01206 //
01207 //  a = tx**2 + ty**2 - (tz*tanPhi)**2
01208 //
01209 //  b = 2*( x0*tx + y0*ty - z0*tz*tanPhi**2 )
01210 //
01211 //  c = x0**2 + y0**2 - r**2 - (z0*tanPhi)**2
01212 //
01213 // 
01214 G4int G4Hype::IntersectHype( const G4ThreeVector &p, const G4ThreeVector &v, 
01215                              G4double r2, G4double tan2Phi, G4double ss[2] )
01216 {
01217   G4double x0 = p.x(), y0 = p.y(), z0 = p.z();
01218   G4double tx = v.x(), ty = v.y(), tz = v.z();
01219 
01220   G4double a = tx*tx + ty*ty - tz*tz*tan2Phi;
01221   G4double b = 2*( x0*tx + y0*ty - z0*tz*tan2Phi );
01222   G4double c = x0*x0 + y0*y0 - r2 - z0*z0*tan2Phi;
01223   
01224   if (std::fabs(a) < DBL_MIN)
01225   {
01226     //
01227     // The trajectory is parallel to the asympotic limit of
01228     // the surface: single solution
01229     //
01230     if (std::fabs(b) < DBL_MIN) return 0;  // Unless we travel through exact center
01231     
01232     ss[0] = c/b;
01233     return 1;
01234   }
01235     
01236   
01237   G4double radical = b*b - 4*a*c;
01238   
01239   if (radical < -DBL_MIN) return 0;    // No solution
01240   
01241   if (radical < DBL_MIN)
01242   {
01243     //
01244     // Grazes surface
01245     //
01246     ss[0] = -b/a/2.0;
01247     return 1;
01248   }
01249   
01250   radical = std::sqrt(radical);
01251   
01252   G4double q = -0.5*( b + (b < 0 ? -radical : +radical) );
01253   G4double sa = q/a;
01254   G4double sb = c/q;    
01255   if (sa < sb) { ss[0] = sa; ss[1] = sb; } else { ss[0] = sb; ss[1] = sa; }
01256   return 2;
01257 }
01258   
01259   
01260 //
01261 // ApproxDistOutside (static)
01262 //
01263 // Find the approximate distance of a point outside
01264 // (greater radius) of a hyperbolic surface. The distance
01265 // must be an underestimate. It will also be nice (although
01266 // not necesary) that the estimate is always finite no
01267 // matter how close the point is.
01268 //
01269 // Our hyperbola approaches the asymptotic limit at z = +/- infinity
01270 // to the lines r = z*tanPhi. We call these lines the 
01271 // asymptotic limit line.
01272 //
01273 // We need the distance of the 2d point p(r,z) to the
01274 // hyperbola r**2 = r0**2 + (z*tanPhi)**2. Find two
01275 // points that bracket the true normal and use the 
01276 // distance to the line that connects these two points.
01277 // The first such point is z=p.z. The second point is
01278 // the z position on the asymptotic limit line that
01279 // contains the normal on the line through the point p.
01280 //
01281 G4double G4Hype::ApproxDistOutside( G4double pr, G4double pz,
01282                                     G4double r0, G4double tanPhi )
01283 {
01284   if (tanPhi < DBL_MIN) return pr-r0;
01285 
01286   G4double tan2Phi = tanPhi*tanPhi;
01287 
01288   //
01289   // First point
01290   //
01291   G4double z1 = pz;
01292   G4double r1 = std::sqrt( r0*r0 + z1*z1*tan2Phi );
01293   
01294   //
01295   // Second point
01296   //
01297   G4double z2 = (pr*tanPhi + pz)/(1 + tan2Phi);
01298   G4double r2 = std::sqrt( r0*r0 + z2*z2*tan2Phi );
01299   
01300   //
01301   // Line between them
01302   //
01303   G4double dr = r2-r1;
01304   G4double dz = z2-z1;
01305   
01306   G4double len = std::sqrt(dr*dr + dz*dz);
01307   if (len < DBL_MIN)
01308   {
01309     //
01310     // The two points are the same?? I guess we
01311     // must have really bracketed the normal
01312     //
01313     dr = pr-r1;
01314     dz = pz-z1;
01315     return std::sqrt( dr*dr + dz*dz );
01316   }
01317   
01318   //
01319   // Distance
01320   //
01321   return std::fabs((pr-r1)*dz - (pz-z1)*dr)/len;
01322 }
01323 
01324 //
01325 // ApproxDistInside (static)
01326 //
01327 // Find the approximate distance of a point inside
01328 // of a hyperbolic surface. The distance
01329 // must be an underestimate. It will also be nice (although
01330 // not necesary) that the estimate is always finite no
01331 // matter how close the point is.
01332 //
01333 // This estimate uses the distance to a line tangent to
01334 // the hyperbolic function. The point of tangent is chosen
01335 // by the z position point
01336 //
01337 // Assumes pr and pz are positive
01338 //
01339 G4double G4Hype::ApproxDistInside( G4double pr, G4double pz,
01340                                    G4double r0, G4double tan2Phi )
01341 {
01342   if (tan2Phi < DBL_MIN) return r0 - pr;
01343 
01344   //
01345   // Corresponding position and normal on hyperbolic
01346   //
01347   G4double rh = std::sqrt( r0*r0 + pz*pz*tan2Phi );
01348   
01349   G4double dr = -rh;
01350   G4double dz = pz*tan2Phi;
01351   G4double len = std::sqrt(dr*dr + dz*dz);
01352   
01353   //
01354   // Answer
01355   //
01356   return std::fabs((pr-rh)*dr)/len;
01357 }
01358 
01359 
01360 //
01361 // GetEntityType
01362 //
01363 G4GeometryType G4Hype::GetEntityType() const
01364 {
01365   return G4String("G4Hype");
01366 }
01367 
01368 
01369 //
01370 // Clone
01371 //
01372 G4VSolid* G4Hype::Clone() const
01373 {
01374   return new G4Hype(*this);
01375 }
01376 
01377 
01378 //
01379 // GetCubicVolume
01380 //
01381 G4double G4Hype::GetCubicVolume()
01382 {
01383   if(fCubicVolume != 0.) {;}
01384     else { fCubicVolume = G4VSolid::GetCubicVolume(); }
01385   return fCubicVolume;
01386 }
01387 
01388 
01389 //
01390 // GetSurfaceArea
01391 //
01392 G4double G4Hype::GetSurfaceArea()
01393 {
01394   if(fSurfaceArea != 0.) {;}
01395   else   { fSurfaceArea = G4VSolid::GetSurfaceArea(); }
01396   return fSurfaceArea;
01397 }
01398 
01399 
01400 //
01401 // Stream object contents to an output stream
01402 //
01403 std::ostream& G4Hype::StreamInfo(std::ostream& os) const
01404 {
01405   G4int oldprc = os.precision(16);
01406   os << "-----------------------------------------------------------\n"
01407      << "    *** Dump for solid - " << GetName() << " ***\n"
01408      << "    ===================================================\n"
01409      << " Solid type: G4Hype\n"
01410      << " Parameters: \n"
01411      << "    half length Z: " << halfLenZ/mm << " mm \n"
01412      << "    inner radius : " << innerRadius/mm << " mm \n"
01413      << "    outer radius : " << outerRadius/mm << " mm \n"
01414      << "    inner stereo angle : " << innerStereo/degree << " degrees \n"
01415      << "    outer stereo angle : " << outerStereo/degree << " degrees \n"
01416      << "-----------------------------------------------------------\n";
01417   os.precision(oldprc);
01418 
01419   return os;
01420 }
01421 
01422 
01423 
01424 //
01425 // GetPointOnSurface
01426 //
01427 G4ThreeVector G4Hype::GetPointOnSurface() const
01428 {
01429   G4double xRand, yRand, zRand, r2 , aOne, aTwo, aThree, chose, sinhu;
01430   G4double phi, cosphi, sinphi, rBar2Out, rBar2In, alpha, t, rOut, rIn2, rOut2;
01431 
01432   // we use the formula of the area of a surface of revolution to compute 
01433   // the areas, using the equation of the hyperbola:
01434   // x^2 + y^2 = (z*tanphi)^2 + r^2
01435 
01436   rBar2Out = outerRadius2;
01437   alpha = 2.*pi*rBar2Out*std::cos(outerStereo)/tanOuterStereo;
01438   t     = halfLenZ*tanOuterStereo/(outerRadius*std::cos(outerStereo));
01439   t     = std::log(t+std::sqrt(sqr(t)+1));
01440   aOne  = std::fabs(2.*alpha*(std::sinh(2.*t)/4.+t/2.));
01441 
01442 
01443   rBar2In = innerRadius2;
01444   alpha = 2.*pi*rBar2In*std::cos(innerStereo)/tanInnerStereo;
01445   t     = halfLenZ*tanInnerStereo/(innerRadius*std::cos(innerStereo));
01446   t     = std::log(t+std::sqrt(sqr(t)+1));
01447   aTwo  = std::fabs(2.*alpha*(std::sinh(2.*t)/4.+t/2.));
01448 
01449   aThree = pi*((outerRadius2+sqr(halfLenZ*tanOuterStereo)
01450               -(innerRadius2+sqr(halfLenZ*tanInnerStereo))));
01451   
01452   if(outerStereo == 0.) {aOne = std::fabs(2.*pi*outerRadius*2.*halfLenZ);}
01453   if(innerStereo == 0.) {aTwo = std::fabs(2.*pi*innerRadius*2.*halfLenZ);}
01454   
01455   phi = RandFlat::shoot(0.,2.*pi);
01456   cosphi = std::cos(phi);
01457   sinphi = std::sin(phi);
01458   sinhu = RandFlat::shoot(-1.*halfLenZ*tanOuterStereo/outerRadius,
01459                           halfLenZ*tanOuterStereo/outerRadius);
01460 
01461   chose = RandFlat::shoot(0.,aOne+aTwo+2.*aThree);
01462   if(chose>=0. && chose < aOne)
01463   {
01464     if(outerStereo != 0.)
01465     {
01466       zRand = outerRadius*sinhu/tanOuterStereo;
01467       xRand = std::sqrt(sqr(sinhu)+1)*outerRadius*cosphi;
01468       yRand = std::sqrt(sqr(sinhu)+1)*outerRadius*sinphi;
01469       return G4ThreeVector (xRand, yRand, zRand);
01470     }
01471     else
01472     {
01473       return G4ThreeVector(outerRadius*cosphi,outerRadius*sinphi,
01474                            RandFlat::shoot(-halfLenZ,halfLenZ));
01475     }
01476   }
01477   else if(chose>=aOne && chose<aOne+aTwo)
01478   {
01479     if(innerStereo != 0.)
01480     {
01481       sinhu = RandFlat::shoot(-1.*halfLenZ*tanInnerStereo/innerRadius,
01482                               halfLenZ*tanInnerStereo/innerRadius);
01483       zRand = innerRadius*sinhu/tanInnerStereo;
01484       xRand = std::sqrt(sqr(sinhu)+1)*innerRadius*cosphi;
01485       yRand = std::sqrt(sqr(sinhu)+1)*innerRadius*sinphi;
01486       return G4ThreeVector (xRand, yRand, zRand);
01487     }
01488     else 
01489     {
01490       return G4ThreeVector(innerRadius*cosphi,innerRadius*sinphi,
01491                            RandFlat::shoot(-1.*halfLenZ,halfLenZ));
01492     }
01493   }
01494   else if(chose>=aOne+aTwo && chose<aOne+aTwo+aThree)
01495   {
01496     rIn2  = innerRadius2+tanInnerStereo2*halfLenZ*halfLenZ;
01497     rOut2 = outerRadius2+tanOuterStereo2*halfLenZ*halfLenZ;
01498     rOut  = std::sqrt(rOut2) ;
01499  
01500     do {
01501       xRand = RandFlat::shoot(-rOut,rOut) ;
01502       yRand = RandFlat::shoot(-rOut,rOut) ;
01503       r2 = xRand*xRand + yRand*yRand ;
01504     } while ( ! ( r2 >= rIn2 && r2 <= rOut2 ) ) ;
01505 
01506     zRand = halfLenZ;
01507     return G4ThreeVector (xRand, yRand, zRand);
01508   }
01509   else
01510   {
01511     rIn2  = innerRadius2+tanInnerStereo2*halfLenZ*halfLenZ;
01512     rOut2 = outerRadius2+tanOuterStereo2*halfLenZ*halfLenZ;
01513     rOut  = std::sqrt(rOut2) ;
01514  
01515     do {
01516       xRand = RandFlat::shoot(-rOut,rOut) ;
01517       yRand = RandFlat::shoot(-rOut,rOut) ;
01518       r2 = xRand*xRand + yRand*yRand ;
01519     } while ( ! ( r2 >= rIn2 && r2 <= rOut2 ) ) ;
01520 
01521     zRand = -1.*halfLenZ;
01522     return G4ThreeVector (xRand, yRand, zRand);
01523   }
01524 }
01525 
01526 
01527 //
01528 // DescribeYourselfTo
01529 //
01530 void G4Hype::DescribeYourselfTo (G4VGraphicsScene& scene) const 
01531 {
01532   scene.AddSolid (*this);
01533 }
01534 
01535 
01536 //
01537 // GetExtent
01538 //
01539 G4VisExtent G4Hype::GetExtent() const 
01540 {
01541   // Define the sides of the box into which the G4Tubs instance would fit.
01542   //
01543   return G4VisExtent( -endOuterRadius, endOuterRadius, 
01544                       -endOuterRadius, endOuterRadius, 
01545                       -halfLenZ, halfLenZ );
01546 }
01547 
01548 
01549 //
01550 // CreatePolyhedron
01551 //
01552 G4Polyhedron* G4Hype::CreatePolyhedron() const 
01553 {
01554    return new G4PolyhedronHype(innerRadius, outerRadius,
01555                                tanInnerStereo2, tanOuterStereo2, halfLenZ);
01556 }
01557 
01558 
01559 //
01560 // GetPolyhedron
01561 //
01562 G4Polyhedron* G4Hype::GetPolyhedron () const
01563 {
01564   if (!fpPolyhedron ||
01565       fpPolyhedron->GetNumberOfRotationStepsAtTimeOfCreation() !=
01566       fpPolyhedron->GetNumberOfRotationSteps())
01567     {
01568       delete fpPolyhedron;
01569       fpPolyhedron = CreatePolyhedron();
01570     }
01571   return fpPolyhedron;
01572 }
01573 
01574 
01575 //
01576 // CreateNURBS
01577 //
01578 G4NURBS* G4Hype::CreateNURBS() const 
01579 {
01580   // Tube for now!!!
01581   //
01582   return new G4NURBStube(endInnerRadius, endOuterRadius, halfLenZ);
01583 }
01584 
01585 
01586 //
01587 //  asinh
01588 //
01589 G4double G4Hype::asinh(G4double arg)
01590 {
01591   return std::log(arg+std::sqrt(sqr(arg)+1));
01592 }

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