G4Para.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: G4Para.cc 69788 2013-05-15 12:06:57Z gcosmo $
00028 //
00029 // class G4Para
00030 //
00031 // Implementation for G4Para class
00032 //
00033 // History:
00034 //
00035 // 23.10.05 V.Grichine: bug fixed in DistanceToOut(p,v,...) for the v.x()<0 case 
00036 // 28.04.05 V.Grichine: new SurfaceNormal according to J. Apostolakis proposal 
00037 // 30.11.04 V.Grichine: modifications in SurfaceNormal for edges/vertices and
00038 //                      in constructor with vertices
00039 // 14.02.02 V.Grichine: bug fixed in Inside according to proposal of D.Wright
00040 // 18.11.99 V.Grichine: kUndef was added to ESide
00041 // 31.10.96 V.Grichine: Modifications according G4Box/Tubs before to commit
00042 // 21.03.95 P.Kent: Modified for `tolerant' geom
00043 //
00045 
00046 #include "G4Para.hh"
00047 
00048 #include "G4VoxelLimits.hh"
00049 #include "G4AffineTransform.hh"
00050 #include "Randomize.hh"
00051 
00052 #include "G4VPVParameterisation.hh"
00053 
00054 #include "G4VGraphicsScene.hh"
00055 #include "G4Polyhedron.hh"
00056 #include "G4NURBS.hh"
00057 #include "G4NURBSbox.hh"
00058 
00059 using namespace CLHEP;
00060 
00061 // Private enum: Not for external use 
00062     
00063 enum ESide {kUndef,kPX,kMX,kPY,kMY,kPZ,kMZ};
00064 
00065 // used internally for normal routine
00066 
00067 enum ENSide {kNZ,kNX,kNY};
00068 
00070 //
00071 // Constructor - check and set half-widths
00072 
00073 void G4Para::SetAllParameters( G4double pDx, G4double pDy, G4double pDz, 
00074                                G4double pAlpha, G4double pTheta, G4double pPhi )
00075 {
00076   if ( pDx > 0 && pDy > 0 && pDz > 0 )
00077   {
00078     fDx         = pDx;
00079     fDy         = pDy;
00080     fDz         = pDz;
00081     fTalpha     = std::tan(pAlpha);
00082     fTthetaCphi = std::tan(pTheta)*std::cos(pPhi);
00083     fTthetaSphi = std::tan(pTheta)*std::sin(pPhi);
00084   }
00085   else
00086   {
00087     std::ostringstream message;
00088     message << "Invalid Length Parameters for Solid: " << GetName() << G4endl
00089             << "        pDx, pDy, pDz = "
00090             << pDx << ", " << pDy << ", " << pDz;
00091     G4Exception("G4Para::SetAllParameters()", "GeomSolids0002",
00092                 FatalException, message);
00093   }
00094   fCubicVolume = 0.;
00095   fSurfaceArea = 0.;
00096   fpPolyhedron = 0;
00097 }
00098 
00100 //
00101 
00102 G4Para::G4Para(const G4String& pName,
00103                      G4double pDx, G4double pDy, G4double pDz,
00104                      G4double pAlpha, G4double pTheta, G4double pPhi)
00105   : G4CSGSolid(pName)
00106 {
00107   if ((pDx<=0) || (pDy<=0) || (pDz<=0))
00108   {
00109     std::ostringstream message;
00110     message << "Invalid Length Parameters for Solid: " << GetName() << G4endl
00111             << "        pDx, pDy, pDz = "
00112             << pDx << ", " << pDy << ", " << pDz;
00113     G4Exception("G4Para::G4Para()", "GeomSolids0002",
00114                 FatalException, message);
00115   }
00116   SetAllParameters( pDx, pDy, pDz, pAlpha, pTheta, pPhi);
00117 }
00118 
00120 //
00121 // Constructor - Design of trapezoid based on 8 G4ThreeVector parameters, 
00122 // which are its vertices. Checking of planarity with preparation of 
00123 // fPlanes[] and than calculation of other members
00124 
00125 G4Para::G4Para( const G4String& pName,
00126                 const G4ThreeVector pt[8] )
00127   : G4CSGSolid(pName)
00128 {
00129   if (!( pt[0].z()<0 && pt[0].z()==pt[1].z() && pt[0].z()==pt[2].z() &&
00130          pt[0].z()==pt[3].z() && pt[4].z()>0 && pt[4].z()==pt[5].z() &&
00131          pt[4].z()==pt[6].z() && pt[4].z()==pt[7].z()           &&
00132         (pt[0].z()+pt[4].z())==0                                &&
00133          pt[0].y()==pt[1].y() && pt[2].y()==pt[3].y()           &&
00134          pt[4].y()==pt[5].y() && pt[6].y()==pt[7].y()           &&
00135        ( pt[0].y() + pt[2].y() + pt[4].y() + pt[6].y() ) == 0   && 
00136        ( pt[0].x() + pt[1].x() + pt[4].x() + pt[5].x() ) == 0) )
00137   {
00138     std::ostringstream message;
00139     message << "Invalid vertice coordinates for Solid: " << GetName();
00140     G4Exception("G4Para::G4Para()", "GeomSolids0002",
00141                 FatalException, message);
00142   }    
00143   fDx = ((pt[3]).x()-(pt[2]).x())*0.5;
00144   fDy = ((pt[2]).y()-(pt[1]).y())*0.5;
00145   fDz = (pt[7]).z();
00146   fTalpha = ((pt[2]).x()+(pt[3]).x()-(pt[1]).x()-(pt[0]).x())*0.25/fDy ;
00147   fTthetaCphi = ((pt[4]).x()+fDy*fTalpha+fDx)/fDz ;
00148   fTthetaSphi = ((pt[4]).y()+fDy)/fDz ;
00149   fCubicVolume = 0.;
00150   fSurfaceArea = 0.;
00151   fpPolyhedron = 0;
00152 }
00153 
00155 //
00156 // Fake default constructor - sets only member data and allocates memory
00157 //                            for usage restricted to object persistency.
00158 //
00159 G4Para::G4Para( __void__& a )
00160   : G4CSGSolid(a), fDx(0.), fDy(0.), fDz(0.),
00161     fTalpha(0.), fTthetaCphi(0.), fTthetaSphi(0.)
00162 {
00163 }
00164 
00166 //
00167 
00168 G4Para::~G4Para()
00169 {
00170 }
00171 
00173 //
00174 // Copy constructor
00175 
00176 G4Para::G4Para(const G4Para& rhs)
00177   : G4CSGSolid(rhs), fDx(rhs.fDx), fDy(rhs.fDy), fDz(rhs.fDz),
00178     fTalpha(rhs.fTalpha), fTthetaCphi(rhs.fTthetaCphi),
00179     fTthetaSphi(rhs.fTthetaSphi)
00180 {
00181 }
00182 
00184 //
00185 // Assignment operator
00186 
00187 G4Para& G4Para::operator = (const G4Para& rhs) 
00188 {
00189    // Check assignment to self
00190    //
00191    if (this == &rhs)  { return *this; }
00192 
00193    // Copy base class data
00194    //
00195    G4CSGSolid::operator=(rhs);
00196 
00197    // Copy data
00198    //
00199    fDx = rhs.fDx; fDy = rhs.fDy; fDz = rhs.fDz;
00200    fTalpha = rhs.fTalpha; fTthetaCphi = rhs.fTthetaCphi;
00201    fTthetaSphi = rhs.fTthetaSphi;
00202 
00203    return *this;
00204 }
00205 
00207 //
00208 // Dispatch to parameterisation for replication mechanism dimension
00209 // computation & modification.
00210 
00211 void G4Para::ComputeDimensions(      G4VPVParameterisation* p,
00212                                 const G4int n,
00213                                 const G4VPhysicalVolume* pRep )
00214 {
00215   p->ComputeDimensions(*this,n,pRep);
00216 }
00217 
00218 
00220 //
00221 // Calculate extent under transform and specified limit
00222 
00223 G4bool G4Para::CalculateExtent( const EAxis pAxis,
00224                                 const G4VoxelLimits& pVoxelLimit,
00225                                 const G4AffineTransform& pTransform,
00226                                      G4double& pMin, G4double& pMax ) const
00227 {
00228   G4bool flag;
00229 
00230   if (!pTransform.IsRotated())
00231   {  
00232     // Special case handling for unrotated trapezoids
00233     // Compute z/x/y/ mins and maxs respecting limits, with early returns
00234     // if outside limits. Then switch() on pAxis
00235 
00236     G4int i ; 
00237     G4double xoffset,xMin,xMax;
00238     G4double yoffset,yMin,yMax;
00239     G4double zoffset,zMin,zMax;
00240     G4double temp[8] ;       // some points for intersection with zMin/zMax
00241     
00242     xoffset=pTransform.NetTranslation().x();      
00243     yoffset=pTransform.NetTranslation().y();
00244     zoffset=pTransform.NetTranslation().z();
00245  
00246     G4ThreeVector pt[8];   // vertices after translation
00247     pt[0]=G4ThreeVector(xoffset-fDz*fTthetaCphi-fDy*fTalpha-fDx,
00248                         yoffset-fDz*fTthetaSphi-fDy,zoffset-fDz);
00249     pt[1]=G4ThreeVector(xoffset-fDz*fTthetaCphi-fDy*fTalpha+fDx,
00250                         yoffset-fDz*fTthetaSphi-fDy,zoffset-fDz);
00251     pt[2]=G4ThreeVector(xoffset-fDz*fTthetaCphi+fDy*fTalpha-fDx,
00252                         yoffset-fDz*fTthetaSphi+fDy,zoffset-fDz);
00253     pt[3]=G4ThreeVector(xoffset-fDz*fTthetaCphi+fDy*fTalpha+fDx,
00254                         yoffset-fDz*fTthetaSphi+fDy,zoffset-fDz);
00255     pt[4]=G4ThreeVector(xoffset+fDz*fTthetaCphi-fDy*fTalpha-fDx,
00256                         yoffset+fDz*fTthetaSphi-fDy,zoffset+fDz);
00257     pt[5]=G4ThreeVector(xoffset+fDz*fTthetaCphi-fDy*fTalpha+fDx,
00258                         yoffset+fDz*fTthetaSphi-fDy,zoffset+fDz);
00259     pt[6]=G4ThreeVector(xoffset+fDz*fTthetaCphi+fDy*fTalpha-fDx,
00260                         yoffset+fDz*fTthetaSphi+fDy,zoffset+fDz);
00261     pt[7]=G4ThreeVector(xoffset+fDz*fTthetaCphi+fDy*fTalpha+fDx,
00262                         yoffset+fDz*fTthetaSphi+fDy,zoffset+fDz);
00263     zMin=zoffset-fDz;
00264     zMax=zoffset+fDz;
00265     if ( pVoxelLimit.IsZLimited() )
00266     {
00267       if   ( (zMin>pVoxelLimit.GetMaxZExtent()+kCarTolerance)
00268           || (zMax<pVoxelLimit.GetMinZExtent()-kCarTolerance) )
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 
00285     temp[0] = pt[0].y()+(pt[4].y()-pt[0].y())
00286                        *(zMin-pt[0].z())/(pt[4].z()-pt[0].z()) ;
00287     temp[1] = pt[0].y()+(pt[4].y()-pt[0].y())
00288                        *(zMax-pt[0].z())/(pt[4].z()-pt[0].z()) ;
00289     temp[2] = pt[2].y()+(pt[6].y()-pt[2].y())
00290                        *(zMin-pt[2].z())/(pt[6].z()-pt[2].z()) ;
00291     temp[3] = pt[2].y()+(pt[6].y()-pt[2].y())
00292                        *(zMax-pt[2].z())/(pt[6].z()-pt[2].z()) ;        
00293     yMax = yoffset - std::fabs(fDz*fTthetaSphi) - fDy - fDy ;
00294     yMin = -yMax ;
00295     for(i=0;i<4;i++)
00296     {
00297       if(temp[i] > yMax) yMax = temp[i] ;
00298       if(temp[i] < yMin) yMin = temp[i] ;
00299     }
00300       
00301     if (pVoxelLimit.IsYLimited())
00302     {
00303       if ( (yMin>pVoxelLimit.GetMaxYExtent()+kCarTolerance)
00304         || (yMax<pVoxelLimit.GetMinYExtent()-kCarTolerance) )
00305       {
00306         return false;
00307       }
00308       else
00309       {
00310         if (yMin<pVoxelLimit.GetMinYExtent())
00311         {
00312           yMin=pVoxelLimit.GetMinYExtent();
00313         }
00314         if (yMax>pVoxelLimit.GetMaxYExtent())
00315         {
00316           yMax=pVoxelLimit.GetMaxYExtent();
00317         }
00318       }
00319     }
00320 
00321     temp[0] = pt[0].x()+(pt[4].x()-pt[0].x())
00322                        *(zMin-pt[0].z())/(pt[4].z()-pt[0].z()) ;
00323     temp[1] = pt[0].x()+(pt[4].x()-pt[0].x())
00324                        *(zMax-pt[0].z())/(pt[4].z()-pt[0].z()) ;
00325     temp[2] = pt[2].x()+(pt[6].x()-pt[2].x())
00326                        *(zMin-pt[2].z())/(pt[6].z()-pt[2].z()) ;
00327     temp[3] = pt[2].x()+(pt[6].x()-pt[2].x())
00328                        *(zMax-pt[2].z())/(pt[6].z()-pt[2].z()) ;
00329     temp[4] = pt[3].x()+(pt[7].x()-pt[3].x())
00330                        *(zMin-pt[3].z())/(pt[7].z()-pt[3].z()) ;
00331     temp[5] = pt[3].x()+(pt[7].x()-pt[3].x())
00332                        *(zMax-pt[3].z())/(pt[7].z()-pt[3].z()) ;
00333     temp[6] = pt[1].x()+(pt[5].x()-pt[1].x())
00334                        *(zMin-pt[1].z())/(pt[5].z()-pt[1].z()) ;
00335     temp[7] = pt[1].x()+(pt[5].x()-pt[1].x())
00336                        *(zMax-pt[1].z())/(pt[5].z()-pt[1].z()) ;
00337 
00338     xMax = xoffset - std::fabs(fDz*fTthetaCphi) - fDx - fDx -fDx - fDx;
00339     xMin = -xMax ;
00340     for(i=0;i<8;i++)
00341     {
00342       if(temp[i] > xMax) xMax = temp[i] ;
00343       if(temp[i] < xMin) xMin = temp[i] ;
00344     }
00345       // xMax/Min = f(yMax/Min) ?
00346     if (pVoxelLimit.IsXLimited())
00347     {
00348       if ( (xMin>pVoxelLimit.GetMaxXExtent()+kCarTolerance)
00349         || (xMax<pVoxelLimit.GetMinXExtent()-kCarTolerance) )
00350       {
00351         return false;
00352       }
00353       else
00354       {
00355         if (xMin<pVoxelLimit.GetMinXExtent())
00356         {
00357           xMin=pVoxelLimit.GetMinXExtent();
00358         }
00359         if (xMax>pVoxelLimit.GetMaxXExtent())
00360         {
00361           xMax=pVoxelLimit.GetMaxXExtent();
00362         }
00363       }
00364     }
00365 
00366     switch (pAxis)
00367     {
00368       case kXAxis:
00369         pMin=xMin;
00370         pMax=xMax;
00371         break;
00372       case kYAxis:
00373         pMin=yMin;
00374         pMax=yMax;
00375         break;
00376       case kZAxis:
00377         pMin=zMin;
00378         pMax=zMax;
00379         break;
00380       default:
00381         break;
00382     }
00383 
00384     pMin-=kCarTolerance;
00385     pMax+=kCarTolerance;
00386     flag = true;
00387   }
00388   else
00389   {
00390     // General rotated case - create and clip mesh to boundaries
00391 
00392     G4bool existsAfterClip=false;
00393     G4ThreeVectorList *vertices;
00394 
00395     pMin=+kInfinity;
00396     pMax=-kInfinity;
00397 
00398     // Calculate rotated vertex coordinates
00399 
00400     vertices=CreateRotatedVertices(pTransform);
00401     ClipCrossSection(vertices,0,pVoxelLimit,pAxis,pMin,pMax);
00402     ClipCrossSection(vertices,4,pVoxelLimit,pAxis,pMin,pMax);
00403     ClipBetweenSections(vertices,0,pVoxelLimit,pAxis,pMin,pMax);
00404       
00405     if (pMin!=kInfinity||pMax!=-kInfinity)
00406     {
00407       existsAfterClip=true;
00408         
00409       // Add 2*tolerance to avoid precision troubles
00410       //
00411       pMin-=kCarTolerance;
00412       pMax+=kCarTolerance;
00413     }
00414     else
00415     {
00416       // Check for case where completely enveloping clipping volume
00417       // If point inside then we are confident that the solid completely
00418       // envelopes the clipping volume. Hence set min/max extents according
00419       // to clipping volume extents along the specified axis.
00420        
00421       G4ThreeVector clipCentre(
00422         (pVoxelLimit.GetMinXExtent()+pVoxelLimit.GetMaxXExtent())*0.5,
00423         (pVoxelLimit.GetMinYExtent()+pVoxelLimit.GetMaxYExtent())*0.5,
00424         (pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5);
00425         
00426       if (Inside(pTransform.Inverse().TransformPoint(clipCentre))!=kOutside)
00427       {
00428         existsAfterClip=true;
00429         pMin=pVoxelLimit.GetMinExtent(pAxis);
00430         pMax=pVoxelLimit.GetMaxExtent(pAxis);
00431       }
00432     }
00433     delete vertices ;          //  'new' in the function called
00434     flag = existsAfterClip ;
00435   }
00436   return flag;
00437 }
00438 
00440 //
00441 // Check in p is inside/on surface/outside solid
00442 
00443 EInside G4Para::Inside( const G4ThreeVector& p ) const
00444 {
00445   G4double xt, yt, yt1;
00446   EInside  in = kOutside;
00447 
00448   yt1 = p.y() - fTthetaSphi*p.z();
00449   yt  = std::fabs(yt1) ;
00450 
00451   // xt = std::fabs( p.x() - fTthetaCphi*p.z() - fTalpha*yt );
00452 
00453   xt = std::fabs( p.x() - fTthetaCphi*p.z() - fTalpha*yt1 );
00454 
00455   if ( std::fabs( p.z() ) <= fDz - kCarTolerance*0.5)
00456   {
00457     if (yt <= fDy - kCarTolerance*0.5)
00458     {
00459       if      ( xt <= fDx - kCarTolerance*0.5 ) in = kInside;
00460       else if ( xt <= fDx + kCarTolerance*0.5 ) in = kSurface;
00461     }
00462     else if ( yt <= fDy + kCarTolerance*0.5)
00463     {
00464       if ( xt <= fDx + kCarTolerance*0.5 ) in = kSurface;  
00465     }
00466   }
00467   else  if ( std::fabs(p.z()) <= fDz + kCarTolerance*0.5 )
00468   {
00469     if ( yt <= fDy + kCarTolerance*0.5)
00470     {
00471       if ( xt <= fDx + kCarTolerance*0.5 ) in = kSurface;  
00472     }
00473   }
00474   return in;
00475 }
00476 
00478 //
00479 // Calculate side nearest to p, and return normal
00480 // If 2+ sides equidistant, first side's normal returned (arbitrarily)
00481 
00482 G4ThreeVector G4Para::SurfaceNormal( const G4ThreeVector& p ) const
00483 {
00484   G4ThreeVector norm, sumnorm(0.,0.,0.);
00485   G4int noSurfaces = 0; 
00486   G4double distx,disty,distz;
00487   G4double newpx,newpy,xshift;
00488   G4double calpha,salpha;      // Sin/Cos(alpha) - needed to recalc G4Parameter
00489   G4double tntheta,cosntheta;  // tan and cos of normal's theta component
00490   G4double ycomp;
00491   G4double delta = 0.5*kCarTolerance;
00492 
00493   newpx  = p.x()-fTthetaCphi*p.z();
00494   newpy  = p.y()-fTthetaSphi*p.z();
00495 
00496   calpha = 1/std::sqrt(1+fTalpha*fTalpha);
00497   if (fTalpha) {salpha = -calpha*fTalpha;} // NOTE: using MINUS std::sin(alpha)
00498   else         {salpha = 0.;}
00499   
00500   //  xshift = newpx*calpha+newpy*salpha;
00501   xshift = newpx - newpy*fTalpha;
00502 
00503   //  distx  = std::fabs(std::fabs(xshift)-fDx*calpha);
00504   distx  = std::fabs(std::fabs(xshift)-fDx);
00505   disty  = std::fabs(std::fabs(newpy)-fDy);
00506   distz  = std::fabs(std::fabs(p.z())-fDz);
00507 
00508   tntheta   = fTthetaCphi*calpha + fTthetaSphi*salpha;
00509   cosntheta = 1/std::sqrt(1+tntheta*tntheta);
00510   ycomp     = 1/std::sqrt(1+fTthetaSphi*fTthetaSphi);
00511 
00512   G4ThreeVector nX  = G4ThreeVector( calpha*cosntheta,
00513                                      salpha*cosntheta,
00514                                     -tntheta*cosntheta);
00515   G4ThreeVector nY  = G4ThreeVector( 0, ycomp,-fTthetaSphi*ycomp);
00516   G4ThreeVector nZ  = G4ThreeVector( 0, 0,  1.0);
00517 
00518   if (distx <= delta)      
00519   {
00520     noSurfaces ++;
00521     if ( xshift >= 0.) {sumnorm += nX;}
00522     else               {sumnorm -= nX;}
00523   }
00524   if (disty <= delta)
00525   {
00526     noSurfaces ++;
00527     if ( newpy >= 0.)  {sumnorm += nY;}
00528     else               {sumnorm -= nY;}
00529   }
00530   if (distz <= delta)  
00531   {
00532     noSurfaces ++;
00533     if ( p.z() >= 0.)  {sumnorm += nZ;}
00534     else               {sumnorm -= nZ;}
00535   }
00536   if ( noSurfaces == 0 )
00537   {
00538 #ifdef G4CSGDEBUG
00539     G4Exception("G4Para::SurfaceNormal(p)", "GeomSolids1002",
00540                 JustWarning, "Point p is not on surface !?" );
00541 #endif 
00542      norm = ApproxSurfaceNormal(p);
00543   }
00544   else if ( noSurfaces == 1 ) {norm = sumnorm;}
00545   else                        {norm = sumnorm.unit();}
00546 
00547   return norm;
00548 }
00549 
00550 
00552 //
00553 // Algorithm for SurfaceNormal() following the original specification
00554 // for points not on the surface
00555 
00556 G4ThreeVector G4Para::ApproxSurfaceNormal( const G4ThreeVector& p ) const
00557 {
00558   ENSide  side;
00559   G4ThreeVector norm;
00560   G4double distx,disty,distz;
00561   G4double newpx,newpy,xshift;
00562   G4double calpha,salpha;  // Sin/Cos(alpha) - needed to recalc G4Parameter 
00563   G4double tntheta,cosntheta;  // tan and cos of normal's theta component
00564   G4double ycomp;
00565 
00566   newpx=p.x()-fTthetaCphi*p.z();
00567   newpy=p.y()-fTthetaSphi*p.z();
00568 
00569   calpha=1/std::sqrt(1+fTalpha*fTalpha);
00570   if (fTalpha)
00571   {
00572     salpha=-calpha/fTalpha;  // NOTE: actually use MINUS std::sin(alpha)
00573   }
00574   else
00575   {
00576     salpha=0;
00577   }
00578 
00579   xshift=newpx*calpha+newpy*salpha;
00580 
00581   distx=std::fabs(std::fabs(xshift)-fDx*calpha);
00582   disty=std::fabs(std::fabs(newpy)-fDy);
00583   distz=std::fabs(std::fabs(p.z())-fDz);
00584     
00585   if (distx<disty)
00586   {
00587     if (distx<distz) {side=kNX;}
00588     else {side=kNZ;}
00589   }
00590   else
00591   {
00592     if (disty<distz) {side=kNY;}
00593     else {side=kNZ;}
00594   }
00595 
00596   switch (side)
00597   {
00598     case kNX:
00599       tntheta=fTthetaCphi*calpha+fTthetaSphi*salpha;
00600       if (xshift<0)
00601       {
00602         cosntheta=-1/std::sqrt(1+tntheta*tntheta);
00603       }
00604       else
00605       {
00606         cosntheta=1/std::sqrt(1+tntheta*tntheta);
00607       }
00608       norm=G4ThreeVector(calpha*cosntheta,salpha*cosntheta,-tntheta*cosntheta);
00609       break;
00610     case kNY:
00611       if (newpy<0)
00612       {
00613         ycomp=-1/std::sqrt(1+fTthetaSphi*fTthetaSphi);
00614       }
00615       else
00616       {
00617         ycomp=1/std::sqrt(1+fTthetaSphi*fTthetaSphi);
00618       }
00619       norm=G4ThreeVector(0,ycomp,-fTthetaSphi*ycomp);
00620       break;
00621     case kNZ:           // Closest to Z
00622       if (p.z()>=0)
00623       {
00624         norm=G4ThreeVector(0,0,1);
00625       }
00626       else
00627       {
00628         norm=G4ThreeVector(0,0,-1);
00629       }
00630       break;
00631   }
00632   return norm;
00633 }
00634 
00636 //
00637 // Calculate distance to shape from outside
00638 // - return kInfinity if no intersection
00639 //
00640 // ALGORITHM:
00641 // For each component, calculate pair of minimum and maximum intersection
00642 // values for which the particle is in the extent of the shape
00643 // - The smallest (MAX minimum) allowed distance of the pairs is intersect
00644 // - Z plane intersectin uses tolerance
00645 // - XZ YZ planes use logic & *SLIGHTLY INCORRECT* tolerance
00646 //   (this saves at least 1 sqrt, 1 multiply and 1 divide... in applicable
00647 //    cases)
00648 // - Note: XZ and YZ planes each divide space into four regions,
00649 //   characterised by ss1 ss2
00650 
00651 G4double G4Para::DistanceToIn( const G4ThreeVector& p,
00652                                const G4ThreeVector& v ) const
00653 {
00654   G4double snxt;    // snxt = default return value
00655   G4double smin,smax;
00656   G4double tmin,tmax;
00657   G4double yt,vy,xt,vx;
00658   G4double max;
00659   //
00660   // Z Intersection range
00661   //
00662   if (v.z()>0)
00663   {
00664     max=fDz-p.z();
00665     if (max>kCarTolerance*0.5)
00666     {
00667       smax=max/v.z();
00668       smin=(-fDz-p.z())/v.z();
00669     }
00670     else
00671     {
00672       return snxt=kInfinity;
00673     }
00674   }
00675   else if (v.z()<0)
00676   {
00677     max=-fDz-p.z();
00678     if (max<-kCarTolerance*0.5)
00679     {
00680       smax=max/v.z();
00681       smin=(fDz-p.z())/v.z();
00682     }
00683     else
00684     {
00685       return snxt=kInfinity;
00686     }
00687   }
00688   else
00689   {
00690     if (std::fabs(p.z())<=fDz) // Inside
00691     {
00692       smin=0;
00693       smax=kInfinity;
00694     }
00695     else
00696     {
00697       return snxt=kInfinity;
00698     }
00699   }
00700     
00701   //
00702   // Y G4Parallel planes intersection
00703   //
00704 
00705   yt=p.y()-fTthetaSphi*p.z();
00706   vy=v.y()-fTthetaSphi*v.z();
00707 
00708   if (vy>0)
00709   {
00710     max=fDy-yt;
00711     if (max>kCarTolerance*0.5)
00712     {
00713       tmax=max/vy;
00714       tmin=(-fDy-yt)/vy;
00715     }
00716     else
00717     {
00718       return snxt=kInfinity;
00719     }
00720   }
00721   else if (vy<0)
00722   {
00723     max=-fDy-yt;
00724     if (max<-kCarTolerance*0.5)
00725     {
00726       tmax=max/vy;
00727       tmin=(fDy-yt)/vy;
00728     }
00729     else
00730     {
00731       return snxt=kInfinity;
00732     }
00733   }
00734   else
00735   {
00736     if (std::fabs(yt)<=fDy)
00737     {
00738       tmin=0;
00739       tmax=kInfinity;
00740     }
00741     else
00742     {
00743       return snxt=kInfinity;
00744     }
00745   }
00746 
00747   // Re-Calc valid intersection range
00748   //
00749   if (tmin>smin) smin=tmin;
00750   if (tmax<smax) smax=tmax;
00751   if (smax<=smin)
00752   {
00753     return snxt=kInfinity;
00754   }
00755   else
00756   {
00757     //
00758     // X G4Parallel planes intersection
00759     //
00760     xt=p.x()-fTthetaCphi*p.z()-fTalpha*yt;
00761     vx=v.x()-fTthetaCphi*v.z()-fTalpha*vy;
00762     if (vx>0)
00763     {
00764       max=fDx-xt;
00765       if (max>kCarTolerance*0.5)
00766       {
00767         tmax=max/vx;
00768         tmin=(-fDx-xt)/vx;
00769       }
00770       else
00771       {
00772         return snxt=kInfinity;
00773       }
00774     }
00775     else if (vx<0)
00776     {
00777       max=-fDx-xt;
00778       if (max<-kCarTolerance*0.5)
00779       {
00780         tmax=max/vx;
00781         tmin=(fDx-xt)/vx;
00782       }
00783       else
00784       {
00785         return snxt=kInfinity;
00786       }
00787     }
00788     else
00789     {
00790       if (std::fabs(xt)<=fDx)
00791       {
00792         tmin=0;
00793         tmax=kInfinity;
00794       }
00795       else
00796       {
00797         return snxt=kInfinity;
00798       }
00799     }
00800     if (tmin>smin) smin=tmin;
00801     if (tmax<smax) smax=tmax;
00802   }
00803 
00804   if (smax>0&&smin<smax)
00805   {
00806     if (smin>0)
00807     {
00808       snxt=smin;
00809     }
00810     else
00811     {
00812       snxt=0;
00813     }
00814   }
00815   else
00816   {
00817     snxt=kInfinity;
00818   }
00819   return snxt;
00820 }
00821 
00823 //
00824 // Calculate exact shortest distance to any boundary from outside
00825 // - Returns 0 is point inside
00826 
00827 G4double G4Para::DistanceToIn( const G4ThreeVector& p ) const
00828 {
00829   G4double safe=0.0;
00830   G4double distz1,distz2,disty1,disty2,distx1,distx2;
00831   G4double trany,cosy,tranx,cosx;
00832 
00833   // Z planes
00834   //
00835   distz1=p.z()-fDz;
00836   distz2=-fDz-p.z();
00837   if (distz1>distz2)
00838   {
00839     safe=distz1;
00840   }
00841   else
00842   {
00843     safe=distz2;
00844   }
00845 
00846   trany=p.y()-fTthetaSphi*p.z(); // Transformed y into `box' system
00847 
00848   // Transformed x into `box' system
00849   //
00850   cosy=1.0/std::sqrt(1.0+fTthetaSphi*fTthetaSphi);
00851   disty1=(trany-fDy)*cosy;
00852   disty2=(-fDy-trany)*cosy;
00853     
00854   if (disty1>safe) safe=disty1;
00855   if (disty2>safe) safe=disty2;
00856 
00857   tranx=p.x()-fTthetaCphi*p.z()-fTalpha*trany;
00858   cosx=1.0/std::sqrt(1.0+fTalpha*fTalpha+fTthetaCphi*fTthetaCphi);
00859   distx1=(tranx-fDx)*cosx;
00860   distx2=(-fDx-tranx)*cosx;
00861     
00862   if (distx1>safe) safe=distx1;
00863   if (distx2>safe) safe=distx2;
00864     
00865   if (safe<0) safe=0;
00866   return safe;  
00867 }
00868 
00870 //
00871 // Calculate distance to surface of shape from inside
00872 // Calculate distance to x/y/z planes - smallest is exiting distance
00873 
00874 G4double G4Para::DistanceToOut(const G4ThreeVector& p, const G4ThreeVector& v,
00875                                const G4bool calcNorm,
00876                                G4bool *validNorm, G4ThreeVector *n) const
00877 {
00878   ESide side = kUndef;
00879   G4double snxt;    // snxt = return value
00880   G4double max,tmax;
00881   G4double yt,vy,xt,vx;
00882 
00883   G4double ycomp,calpha,salpha,tntheta,cosntheta;
00884 
00885   //
00886   // Z Intersections
00887   //
00888 
00889   if (v.z()>0)
00890   {
00891     max=fDz-p.z();
00892     if (max>kCarTolerance*0.5)
00893     {
00894       snxt=max/v.z();
00895       side=kPZ;
00896     }
00897     else
00898     {
00899       if (calcNorm)
00900       {
00901         *validNorm=true;
00902         *n=G4ThreeVector(0,0,1);
00903       }
00904       return snxt=0;
00905     }
00906   }
00907   else if (v.z()<0)
00908   {
00909     max=-fDz-p.z();
00910     if (max<-kCarTolerance*0.5)
00911     {
00912       snxt=max/v.z();
00913       side=kMZ;
00914     }
00915     else
00916     {
00917       if (calcNorm)
00918       {
00919         *validNorm=true;
00920         *n=G4ThreeVector(0,0,-1);
00921       }
00922       return snxt=0;
00923     }
00924   }
00925   else
00926   {
00927     snxt=kInfinity;
00928   }
00929     
00930   //
00931   // Y plane intersection
00932   //
00933 
00934   yt=p.y()-fTthetaSphi*p.z();
00935   vy=v.y()-fTthetaSphi*v.z();
00936 
00937   if (vy>0)
00938   {
00939     max=fDy-yt;
00940     if (max>kCarTolerance*0.5)
00941     {
00942       tmax=max/vy;
00943       if (tmax<snxt)
00944       {
00945         snxt=tmax;
00946         side=kPY;
00947       }
00948     }
00949     else
00950     {
00951       if (calcNorm)
00952       {      
00953         *validNorm=true; // Leaving via plus Y
00954         ycomp=1/std::sqrt(1+fTthetaSphi*fTthetaSphi);
00955         *n=G4ThreeVector(0,ycomp,-fTthetaSphi*ycomp);
00956       }
00957       return snxt=0;
00958     }
00959   }
00960   else if (vy<0)
00961   {
00962     max=-fDy-yt;
00963     if (max<-kCarTolerance*0.5)
00964     {
00965       tmax=max/vy;
00966       if (tmax<snxt)
00967       {
00968         snxt=tmax;
00969         side=kMY;
00970       }
00971     }
00972     else
00973     {
00974       if (calcNorm)
00975       {
00976         *validNorm=true; // Leaving via minus Y
00977         ycomp=-1/std::sqrt(1+fTthetaSphi*fTthetaSphi);
00978         *n=G4ThreeVector(0,ycomp,-fTthetaSphi*ycomp);
00979       }
00980       return snxt=0;
00981     }
00982   }
00983 
00984   //
00985   // X plane intersection
00986   //
00987 
00988   xt=p.x()-fTthetaCphi*p.z()-fTalpha*yt;
00989   vx=v.x()-fTthetaCphi*v.z()-fTalpha*vy;
00990   if (vx>0)
00991   {
00992     max=fDx-xt;
00993     if (max>kCarTolerance*0.5)
00994     {
00995       tmax=max/vx;
00996       if (tmax<snxt)
00997       {
00998         snxt=tmax;
00999         side=kPX;
01000       }
01001     }
01002     else
01003     {
01004       if (calcNorm)
01005       {
01006         *validNorm=true; // Leaving via plus X
01007         calpha=1/std::sqrt(1+fTalpha*fTalpha);
01008         if (fTalpha)
01009         {
01010           salpha=-calpha/fTalpha;  // NOTE: actually use MINUS std::sin(alpha)
01011         }
01012         else
01013         {
01014           salpha=0;
01015         }
01016         tntheta=fTthetaCphi*calpha+fTthetaSphi*salpha;
01017         cosntheta=1/std::sqrt(1+tntheta*tntheta);
01018         *n=G4ThreeVector(calpha*cosntheta,salpha*cosntheta,-tntheta*cosntheta);
01019       }
01020       return snxt=0;
01021     }
01022   }
01023   else if (vx<0)
01024   {
01025     max=-fDx-xt;
01026     if (max<-kCarTolerance*0.5)
01027     {
01028       tmax=max/vx;
01029       if (tmax<snxt)
01030       {
01031         snxt=tmax;
01032         side=kMX;
01033       }
01034     }
01035     else
01036     {
01037       if (calcNorm)
01038       {
01039         *validNorm=true; // Leaving via minus X
01040         calpha=1/std::sqrt(1+fTalpha*fTalpha);
01041         if (fTalpha)
01042         {
01043           salpha=-calpha/fTalpha;  // NOTE: actually use MINUS std::sin(alpha)
01044         }
01045         else
01046         {
01047           salpha=0;
01048         }
01049         tntheta=fTthetaCphi*calpha+fTthetaSphi*salpha;
01050         cosntheta=-1/std::sqrt(1+tntheta*tntheta);
01051         *n=G4ThreeVector(calpha*cosntheta,salpha*cosntheta,-tntheta*cosntheta);
01052       }
01053       return snxt=0;
01054     }
01055   }
01056 
01057   if (calcNorm)
01058   {
01059     *validNorm=true;
01060     switch (side)
01061     {
01062       case kMZ:
01063         *n=G4ThreeVector(0,0,-1);
01064         break;
01065       case kPZ:
01066         *n=G4ThreeVector(0,0,1);
01067         break;
01068       case kMY:
01069         ycomp=-1/std::sqrt(1+fTthetaSphi*fTthetaSphi);
01070         *n=G4ThreeVector(0,ycomp,-fTthetaSphi*ycomp);
01071         break;        
01072       case kPY:
01073         ycomp=1/std::sqrt(1+fTthetaSphi*fTthetaSphi);
01074         *n=G4ThreeVector(0,ycomp,-fTthetaSphi*ycomp);
01075         break;        
01076       case kMX:
01077         calpha=1/std::sqrt(1+fTalpha*fTalpha);
01078         if (fTalpha)
01079         {
01080           salpha=-calpha/fTalpha;  // NOTE: actually use MINUS std::sin(alpha)
01081         }
01082         else
01083         {
01084           salpha=0;
01085         }
01086         tntheta=fTthetaCphi*calpha+fTthetaSphi*salpha;
01087         cosntheta=-1/std::sqrt(1+tntheta*tntheta);
01088         *n=G4ThreeVector(calpha*cosntheta,salpha*cosntheta,-tntheta*cosntheta);
01089         break;
01090       case kPX:
01091         calpha=1/std::sqrt(1+fTalpha*fTalpha);
01092         if (fTalpha)
01093         {
01094           salpha=-calpha/fTalpha;  // NOTE: actually use MINUS std::sin(alpha)
01095         }
01096         else
01097         {
01098           salpha=0;
01099         }
01100         tntheta=fTthetaCphi*calpha+fTthetaSphi*salpha;
01101         cosntheta=1/std::sqrt(1+tntheta*tntheta);
01102         *n=G4ThreeVector(calpha*cosntheta,salpha*cosntheta,-tntheta*cosntheta);
01103         break;
01104       default:
01105         DumpInfo();
01106         G4Exception("G4Para::DistanceToOut(p,v,..)",
01107                     "GeomSolids1002",JustWarning,
01108                     "Undefined side for valid surface normal to solid.");
01109         break;
01110     }
01111   }
01112   return snxt;
01113 }
01114 
01116 //
01117 // Calculate exact shortest distance to any boundary from inside
01118 // - Returns 0 is point outside
01119 
01120 G4double G4Para::DistanceToOut( const G4ThreeVector& p ) const
01121 {
01122   G4double safe=0.0;
01123   G4double distz1,distz2,disty1,disty2,distx1,distx2;
01124   G4double trany,cosy,tranx,cosx;
01125 
01126 #ifdef G4CSGDEBUG
01127   if( Inside(p) == kOutside )
01128   {
01129      G4int oldprc = G4cout.precision(16) ;
01130      G4cout << G4endl ;
01131      DumpInfo();
01132      G4cout << "Position:"  << G4endl << G4endl ;
01133      G4cout << "p.x() = "   << p.x()/mm << " mm" << G4endl ;
01134      G4cout << "p.y() = "   << p.y()/mm << " mm" << G4endl ;
01135      G4cout << "p.z() = "   << p.z()/mm << " mm" << G4endl << G4endl ;
01136      G4cout.precision(oldprc) ;
01137      G4Exception("G4Para::DistanceToOut(p)", "GeomSolids1002",
01138                  JustWarning, "Point p is outside !?" );
01139   }
01140 #endif
01141 
01142   // Z planes
01143   //
01144   distz1=fDz-p.z();
01145   distz2=fDz+p.z();
01146   if (distz1<distz2)
01147   {
01148     safe=distz1;
01149   }
01150   else
01151   {
01152     safe=distz2;
01153   }
01154 
01155   trany=p.y()-fTthetaSphi*p.z(); // Transformed y into `box' system
01156 
01157   // Transformed x into `box' system
01158   //
01159   cosy=1.0/std::sqrt(1.0+fTthetaSphi*fTthetaSphi);
01160   disty1=(fDy-trany)*cosy;
01161   disty2=(fDy+trany)*cosy;
01162     
01163   if (disty1<safe) safe=disty1;
01164   if (disty2<safe) safe=disty2;
01165 
01166   tranx=p.x()-fTthetaCphi*p.z()-fTalpha*trany;
01167   cosx=1.0/std::sqrt(1.0+fTalpha*fTalpha+fTthetaCphi*fTthetaCphi);
01168   distx1=(fDx-tranx)*cosx;
01169   distx2=(fDx+tranx)*cosx;
01170     
01171   if (distx1<safe) safe=distx1;
01172   if (distx2<safe) safe=distx2;
01173     
01174   if (safe<0) safe=0;
01175   return safe;  
01176 }
01177 
01179 //
01180 // Create a List containing the transformed vertices
01181 // Ordering [0-3] -fDz cross section
01182 //          [4-7] +fDz cross section such that [0] is below [4],
01183 //                                             [1] below [5] etc.
01184 // Note:
01185 //  Caller has deletion resposibility
01186 
01187 G4ThreeVectorList*
01188 G4Para::CreateRotatedVertices( const G4AffineTransform& pTransform ) const
01189 {
01190   G4ThreeVectorList *vertices;
01191   vertices=new G4ThreeVectorList();
01192   if (vertices)
01193   {
01194     vertices->reserve(8);
01195     G4ThreeVector vertex0(-fDz*fTthetaCphi-fDy*fTalpha-fDx,
01196                           -fDz*fTthetaSphi-fDy, -fDz);
01197     G4ThreeVector vertex1(-fDz*fTthetaCphi-fDy*fTalpha+fDx,
01198                           -fDz*fTthetaSphi-fDy, -fDz);
01199     G4ThreeVector vertex2(-fDz*fTthetaCphi+fDy*fTalpha-fDx,
01200                           -fDz*fTthetaSphi+fDy, -fDz);
01201     G4ThreeVector vertex3(-fDz*fTthetaCphi+fDy*fTalpha+fDx,
01202                           -fDz*fTthetaSphi+fDy, -fDz);
01203     G4ThreeVector vertex4(+fDz*fTthetaCphi-fDy*fTalpha-fDx,
01204                           +fDz*fTthetaSphi-fDy, +fDz);
01205     G4ThreeVector vertex5(+fDz*fTthetaCphi-fDy*fTalpha+fDx,
01206                           +fDz*fTthetaSphi-fDy, +fDz);
01207     G4ThreeVector vertex6(+fDz*fTthetaCphi+fDy*fTalpha-fDx,
01208                           +fDz*fTthetaSphi+fDy, +fDz);
01209     G4ThreeVector vertex7(+fDz*fTthetaCphi+fDy*fTalpha+fDx,
01210                           +fDz*fTthetaSphi+fDy, +fDz);
01211 
01212     vertices->push_back(pTransform.TransformPoint(vertex0));
01213     vertices->push_back(pTransform.TransformPoint(vertex1));
01214     vertices->push_back(pTransform.TransformPoint(vertex2));
01215     vertices->push_back(pTransform.TransformPoint(vertex3));
01216     vertices->push_back(pTransform.TransformPoint(vertex4));
01217     vertices->push_back(pTransform.TransformPoint(vertex5));
01218     vertices->push_back(pTransform.TransformPoint(vertex6));
01219     vertices->push_back(pTransform.TransformPoint(vertex7));
01220   }
01221   else
01222   {
01223     DumpInfo();
01224     G4Exception("G4Para::CreateRotatedVertices()",
01225                 "GeomSolids0003", FatalException,
01226                 "Error in allocation of vertices. Out of memory !");
01227   }
01228   return vertices;
01229 }
01230 
01232 //
01233 // GetEntityType
01234 
01235 G4GeometryType G4Para::GetEntityType() const
01236 {
01237   return G4String("G4Para");
01238 }
01239 
01241 //
01242 // Make a clone of the object
01243 //
01244 G4VSolid* G4Para::Clone() const
01245 {
01246   return new G4Para(*this);
01247 }
01248 
01250 //
01251 // Stream object contents to an output stream
01252 
01253 std::ostream& G4Para::StreamInfo( std::ostream& os ) const
01254 {
01255   G4int oldprc = os.precision(16);
01256   os << "-----------------------------------------------------------\n"
01257      << "    *** Dump for solid - " << GetName() << " ***\n"
01258      << "    ===================================================\n"
01259      << " Solid type: G4Para\n"
01260      << " Parameters: \n"
01261      << "    half length X: " << fDx/mm << " mm \n"
01262      << "    half length Y: " << fDy/mm << " mm \n"
01263      << "    half length Z: " << fDz/mm << " mm \n"
01264      << "    std::tan(alpha)         : " << fTalpha/degree << " degrees \n"
01265      << "    std::tan(theta)*std::cos(phi): " << fTthetaCphi/degree
01266      << " degrees \n"
01267      << "    std::tan(theta)*std::sin(phi): " << fTthetaSphi/degree
01268      << " degrees \n"
01269      << "-----------------------------------------------------------\n";
01270   os.precision(oldprc);
01271 
01272   return os;
01273 }
01274 
01276 //
01277 // GetPointOnPlane
01278 // Auxiliary method for Get Point on Surface
01279 //
01280 
01281 G4ThreeVector G4Para::GetPointOnPlane(G4ThreeVector p0, G4ThreeVector p1, 
01282                                       G4ThreeVector p2, G4ThreeVector p3, 
01283                                       G4double& area) const
01284 {
01285   G4double lambda1, lambda2, chose, aOne, aTwo;
01286   G4ThreeVector t, u, v, w, Area, normal;
01287   
01288   t = p1 - p0;
01289   u = p2 - p1;
01290   v = p3 - p2;
01291   w = p0 - p3;
01292 
01293   Area = G4ThreeVector(w.y()*v.z() - w.z()*v.y(),
01294                        w.z()*v.x() - w.x()*v.z(),
01295                        w.x()*v.y() - w.y()*v.x());
01296   
01297   aOne = 0.5*Area.mag();
01298   
01299   Area = G4ThreeVector(t.y()*u.z() - t.z()*u.y(),
01300                        t.z()*u.x() - t.x()*u.z(),
01301                        t.x()*u.y() - t.y()*u.x());
01302   
01303   aTwo = 0.5*Area.mag();
01304   
01305   area = aOne + aTwo;
01306   
01307   chose = RandFlat::shoot(0.,aOne+aTwo);
01308 
01309   if( (chose>=0.) && (chose < aOne) )
01310   {
01311     lambda1 = RandFlat::shoot(0.,1.);
01312     lambda2 = RandFlat::shoot(0.,lambda1);
01313     return (p2+lambda1*v+lambda2*w);    
01314   }
01315 
01316   // else
01317 
01318   lambda1 = RandFlat::shoot(0.,1.);
01319   lambda2 = RandFlat::shoot(0.,lambda1);
01320   return (p0+lambda1*t+lambda2*u);    
01321 }
01322 
01324 //
01325 // GetPointOnSurface
01326 //
01327 // Return a point (G4ThreeVector) randomly and uniformly
01328 // selected on the solid surface
01329 
01330 G4ThreeVector G4Para::GetPointOnSurface() const
01331 {
01332   G4ThreeVector One, Two, Three, Four, Five, Six;
01333   G4ThreeVector pt[8] ;
01334   G4double chose, aOne, aTwo, aThree, aFour, aFive, aSix;
01335 
01336   pt[0] = G4ThreeVector(-fDz*fTthetaCphi-fDy*fTalpha-fDx,
01337                         -fDz*fTthetaSphi-fDy, -fDz);
01338   pt[1] = G4ThreeVector(-fDz*fTthetaCphi-fDy*fTalpha+fDx,
01339                         -fDz*fTthetaSphi-fDy, -fDz);
01340   pt[2] = G4ThreeVector(-fDz*fTthetaCphi+fDy*fTalpha-fDx,
01341                         -fDz*fTthetaSphi+fDy, -fDz);
01342   pt[3] = G4ThreeVector(-fDz*fTthetaCphi+fDy*fTalpha+fDx,
01343                         -fDz*fTthetaSphi+fDy, -fDz);
01344   pt[4] = G4ThreeVector(+fDz*fTthetaCphi-fDy*fTalpha-fDx,
01345                         +fDz*fTthetaSphi-fDy, +fDz);
01346   pt[5] = G4ThreeVector(+fDz*fTthetaCphi-fDy*fTalpha+fDx,
01347                         +fDz*fTthetaSphi-fDy, +fDz);
01348   pt[6] = G4ThreeVector(+fDz*fTthetaCphi+fDy*fTalpha-fDx,
01349                         +fDz*fTthetaSphi+fDy, +fDz);
01350   pt[7] = G4ThreeVector(+fDz*fTthetaCphi+fDy*fTalpha+fDx,
01351                         +fDz*fTthetaSphi+fDy, +fDz);
01352 
01353   // make sure we provide the points in a clockwise fashion
01354 
01355   One   = GetPointOnPlane(pt[0],pt[1],pt[3],pt[2], aOne);
01356   Two   = GetPointOnPlane(pt[4],pt[5],pt[7],pt[6], aTwo);
01357   Three = GetPointOnPlane(pt[6],pt[7],pt[3],pt[2], aThree);
01358   Four  = GetPointOnPlane(pt[4],pt[5],pt[1],pt[0], aFour); 
01359   Five  = GetPointOnPlane(pt[0],pt[2],pt[6],pt[4], aFive);
01360   Six   = GetPointOnPlane(pt[1],pt[3],pt[7],pt[5], aSix);
01361 
01362   chose = RandFlat::shoot(0.,aOne+aTwo+aThree+aFour+aFive+aSix);
01363   
01364   if( (chose>=0.) && (chose<aOne) )                    
01365     { return One; }
01366   else if(chose>=aOne && chose<aOne+aTwo)  
01367     { return Two; }
01368   else if(chose>=aOne+aTwo && chose<aOne+aTwo+aThree)
01369     { return Three; }
01370   else if(chose>=aOne+aTwo+aThree && chose<aOne+aTwo+aThree+aFour)
01371     { return Four; }
01372   else if(chose>=aOne+aTwo+aThree+aFour && chose<aOne+aTwo+aThree+aFour+aFive)
01373     { return Five; }
01374   return Six;
01375 }
01376 
01378 //
01379 // Methods for visualisation
01380 
01381 void G4Para::DescribeYourselfTo ( G4VGraphicsScene& scene ) const
01382 {
01383   scene.AddSolid (*this);
01384 }
01385 
01386 G4Polyhedron* G4Para::CreatePolyhedron () const
01387 {
01388   G4double phi = std::atan2(fTthetaSphi, fTthetaCphi);
01389   G4double alpha = std::atan(fTalpha);
01390   G4double theta = std::atan(std::sqrt(fTthetaCphi*fTthetaCphi
01391                             +fTthetaSphi*fTthetaSphi));
01392     
01393   return new G4PolyhedronPara(fDx, fDy, fDz, alpha, theta, phi);
01394 }
01395 
01396 G4NURBS* G4Para::CreateNURBS () const
01397 {
01398   // return new G4NURBSbox (fDx, fDy, fDz);
01399   return 0 ;
01400 }

Generated on Mon May 27 17:49:14 2013 for Geant4 by  doxygen 1.4.7