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

Generated on Mon May 27 17:47:59 2013 for Geant4 by  doxygen 1.4.7