G4Trd.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: G4Trd.cc 69788 2013-05-15 12:06:57Z gcosmo $
00028 //
00029 //
00030 // Implementation for G4Trd class
00031 //
00032 // History:
00033 //
00034 // 28.04.05 V.Grichine: new SurfaceNormal according to J. Apostolakis proposal 
00035 // 26.04.05, V.Grichine, new SurfaceNoramal is default
00036 // 07.12.04, V.Grichine, SurfaceNoramal with edges/vertices.
00037 // 07.05.00, V.Grichine, in d = DistanceToIn(p,v), if d<0.5*kCarTolerance, d=0
00038 //    ~1996, V.Grichine, 1st implementation based on old code of P.Kent
00039 //
00041 
00042 #include "G4Trd.hh"
00043 
00044 #include "G4VPVParameterisation.hh"
00045 #include "G4VoxelLimits.hh"
00046 #include "G4AffineTransform.hh"
00047 #include "Randomize.hh"
00048 
00049 #include "G4VGraphicsScene.hh"
00050 #include "G4Polyhedron.hh"
00051 #include "G4NURBS.hh"
00052 #include "G4NURBSbox.hh"
00053 
00054 using namespace CLHEP;
00055 
00057 //
00058 // Constructor - check & set half widths
00059 
00060 G4Trd::G4Trd( const G4String& pName,
00061                     G4double pdx1,  G4double pdx2,
00062                     G4double pdy1,  G4double pdy2,
00063                     G4double pdz )
00064   : G4CSGSolid(pName)
00065 {
00066   CheckAndSetAllParameters (pdx1, pdx2, pdy1, pdy2, pdz);
00067 }
00068 
00070 //
00071 // Set and check (coplanarity) of trd parameters
00072 
00073 void G4Trd::CheckAndSetAllParameters ( G4double pdx1,  G4double pdx2,
00074                                        G4double pdy1,  G4double pdy2,
00075                                        G4double pdz ) 
00076 {
00077   if ( pdx1>0&&pdx2>0&&pdy1>0&&pdy2>0&&pdz>0 )
00078   {
00079     fDx1=pdx1; fDx2=pdx2;
00080     fDy1=pdy1; fDy2=pdy2;
00081     fDz=pdz;
00082   }
00083   else
00084   {
00085     if ( pdx1>=0 && pdx2>=0 && pdy1>=0 && pdy2>=0 && pdz>=0 )
00086     {
00087       // G4double  Minimum_length= (1+per_thousand) * kCarTolerance/2.;
00088       // FIX-ME : temporary solution for ZERO or very-small parameters
00089       //
00090       G4double  Minimum_length= kCarTolerance/2.;
00091       fDx1=std::max(pdx1,Minimum_length); 
00092       fDx2=std::max(pdx2,Minimum_length); 
00093       fDy1=std::max(pdy1,Minimum_length); 
00094       fDy2=std::max(pdy2,Minimum_length); 
00095       fDz=std::max(pdz,Minimum_length);
00096     }
00097     else
00098     {
00099       std::ostringstream message;
00100       message << "Invalid negative dimensions for Solid: " << GetName()
00101               << G4endl
00102               << "          X - " << pdx1 << ", " << pdx2 << G4endl
00103               << "          Y - " << pdy1 << ", " << pdy2 << G4endl
00104               << "          Z - " << pdz;
00105       G4Exception("G4Trd::CheckAndSetAllParameters()",
00106                   "GeomSolids0002", FatalException, message);
00107     }
00108   }
00109   fCubicVolume= 0.;
00110   fSurfaceArea= 0.;
00111   fpPolyhedron = 0;
00112 }
00113 
00115 //
00116 // Fake default constructor - sets only member data and allocates memory
00117 //                            for usage restricted to object persistency.
00118 //
00119 G4Trd::G4Trd( __void__& a )
00120   : G4CSGSolid(a), fDx1(0.), fDx2(0.), fDy1(0.), fDy2(0.), fDz(0.)
00121 {
00122 }
00123 
00125 //
00126 // Destructor
00127 
00128 G4Trd::~G4Trd()
00129 {
00130 }
00131 
00133 //
00134 // Copy constructor
00135 
00136 G4Trd::G4Trd(const G4Trd& rhs)
00137   : G4CSGSolid(rhs), fDx1(rhs.fDx1), fDx2(rhs.fDx2),
00138     fDy1(rhs.fDy1), fDy2(rhs.fDy2), fDz(rhs.fDz)
00139 {
00140 }
00141 
00143 //
00144 // Assignment operator
00145 
00146 G4Trd& G4Trd::operator = (const G4Trd& rhs) 
00147 {
00148    // Check assignment to self
00149    //
00150    if (this == &rhs)  { return *this; }
00151 
00152    // Copy base class data
00153    //
00154    G4CSGSolid::operator=(rhs);
00155 
00156    // Copy data
00157    //
00158    fDx1 = rhs.fDx1; fDx2 = rhs.fDx2;
00159    fDy1 = rhs.fDy1; fDy2 = rhs.fDy2;
00160    fDz = rhs.fDz;
00161 
00162    return *this;
00163 }
00164 
00166 //
00167 //
00168 
00169 void G4Trd::SetAllParameters ( G4double pdx1, G4double pdx2, G4double pdy1, 
00170                                G4double pdy2, G4double pdz ) 
00171 {
00172   CheckAndSetAllParameters (pdx1, pdx2, pdy1, pdy2, pdz);
00173 }
00174 
00175 
00177 //
00178 // Dispatch to parameterisation for replication mechanism dimension
00179 // computation & modification.
00180 
00181 void G4Trd::ComputeDimensions(       G4VPVParameterisation* p,
00182                                const G4int n,
00183                                const G4VPhysicalVolume* pRep )
00184 {
00185   p->ComputeDimensions(*this,n,pRep);
00186 }
00187 
00188 
00190 //
00191 // Calculate extent under transform and specified limit
00192 
00193 G4bool G4Trd::CalculateExtent( const EAxis pAxis,
00194                                const G4VoxelLimits& pVoxelLimit,
00195                                const G4AffineTransform& pTransform,
00196                                      G4double& pMin, G4double& pMax ) const
00197 {
00198   if (!pTransform.IsRotated())
00199   {
00200     // Special case handling for unrotated solids
00201     // Compute x/y/z mins and maxs respecting limits, with early returns
00202     // if outside limits. Then switch() on pAxis
00203 
00204     G4double xoffset,xMin,xMax;
00205     G4double yoffset,yMin,yMax;
00206     G4double zoffset,zMin,zMax;
00207 
00208     zoffset=pTransform.NetTranslation().z();
00209     zMin=zoffset-fDz;
00210     zMax=zoffset+fDz;
00211     if (pVoxelLimit.IsZLimited())
00212     {
00213       if ( (zMin>pVoxelLimit.GetMaxZExtent()+kCarTolerance)
00214         || (zMax<pVoxelLimit.GetMinZExtent()-kCarTolerance) )
00215       {
00216         return false;
00217       }
00218         else
00219       {
00220         if (zMin<pVoxelLimit.GetMinZExtent())
00221         {
00222           zMin=pVoxelLimit.GetMinZExtent();
00223         }
00224         if (zMax>pVoxelLimit.GetMaxZExtent())
00225         {
00226           zMax=pVoxelLimit.GetMaxZExtent();
00227         }
00228       }
00229     }
00230     xoffset=pTransform.NetTranslation().x();
00231     if (fDx2 >= fDx1)
00232     { 
00233       xMax =  xoffset+(fDx1+fDx2)/2+(zMax-zoffset)*(fDx2-fDx1)/(2*fDz) ;
00234       xMin = 2*xoffset - xMax ;
00235     }
00236     else
00237     {
00238       xMax =  xoffset+(fDx1+fDx2)/2+(zMin-zoffset)*(fDx2-fDx1)/(2*fDz) ;
00239       xMin =  2*xoffset - xMax ;
00240     }   
00241     if (pVoxelLimit.IsXLimited())
00242     {
00243       if ( (xMin>pVoxelLimit.GetMaxXExtent()+kCarTolerance)
00244         || (xMax<pVoxelLimit.GetMinXExtent()-kCarTolerance) )
00245       {
00246         return false;
00247       }
00248       else
00249       {
00250         if (xMin<pVoxelLimit.GetMinXExtent())
00251         {
00252           xMin=pVoxelLimit.GetMinXExtent();
00253         }
00254         if (xMax>pVoxelLimit.GetMaxXExtent())
00255         {
00256           xMax=pVoxelLimit.GetMaxXExtent();
00257         }
00258       }
00259     }
00260     yoffset= pTransform.NetTranslation().y() ;
00261     if(fDy2 >= fDy1)
00262     {
00263       yMax = yoffset+(fDy2+fDy1)/2+(zMax-zoffset)*(fDy2-fDy1)/(2*fDz) ;
00264       yMin = 2*yoffset - yMax ;
00265     }
00266     else
00267     {
00268       yMax = yoffset+(fDy2+fDy1)/2+(zMin-zoffset)*(fDy2-fDy1)/(2*fDz) ;
00269       yMin = 2*yoffset - yMax ;  
00270     }  
00271     if (pVoxelLimit.IsYLimited())
00272     {
00273       if ( (yMin>pVoxelLimit.GetMaxYExtent()+kCarTolerance)
00274         || (yMax<pVoxelLimit.GetMinYExtent()-kCarTolerance) )
00275       {
00276         return false;
00277       }
00278       else
00279       {
00280         if (yMin<pVoxelLimit.GetMinYExtent())
00281         {
00282           yMin=pVoxelLimit.GetMinYExtent();
00283         }
00284         if (yMax>pVoxelLimit.GetMaxYExtent())
00285         {
00286           yMax=pVoxelLimit.GetMaxYExtent();
00287         }
00288       }
00289     }
00290 
00291     switch (pAxis)
00292     {
00293       case kXAxis:
00294         pMin=xMin;
00295         pMax=xMax;
00296         break;
00297       case kYAxis:
00298         pMin=yMin;
00299         pMax=yMax;
00300         break;
00301       case kZAxis:
00302         pMin=zMin;
00303         pMax=zMax;
00304         break;
00305       default:
00306         break;
00307     }
00308 
00309     // Add 2*Tolerance to avoid precision troubles ?
00310     //
00311     pMin-=kCarTolerance;
00312     pMax+=kCarTolerance;
00313 
00314     return true;
00315   }
00316   else
00317   {
00318     // General rotated case - create and clip mesh to boundaries
00319 
00320     G4bool existsAfterClip=false;
00321     G4ThreeVectorList *vertices;
00322 
00323     pMin=+kInfinity;
00324     pMax=-kInfinity;
00325 
00326     // Calculate rotated vertex coordinates
00327     //
00328     vertices=CreateRotatedVertices(pTransform);
00329     ClipCrossSection(vertices,0,pVoxelLimit,pAxis,pMin,pMax);
00330     ClipCrossSection(vertices,4,pVoxelLimit,pAxis,pMin,pMax);
00331     ClipBetweenSections(vertices,0,pVoxelLimit,pAxis,pMin,pMax);
00332       
00333     if (pMin!=kInfinity||pMax!=-kInfinity)
00334     {
00335       existsAfterClip=true;
00336 
00337       // Add 2*tolerance to avoid precision troubles
00338       //
00339       pMin-=kCarTolerance;
00340       pMax+=kCarTolerance;
00341         
00342     }
00343     else
00344     {
00345       // Check for case where completely enveloping clipping volume
00346       // If point inside then we are confident that the solid completely
00347       // envelopes the clipping volume. Hence set min/max extents according
00348       // to clipping volume extents along the specified axis.
00349 
00350       G4ThreeVector clipCentre(
00351          (pVoxelLimit.GetMinXExtent()+pVoxelLimit.GetMaxXExtent())*0.5,
00352          (pVoxelLimit.GetMinYExtent()+pVoxelLimit.GetMaxYExtent())*0.5,
00353          (pVoxelLimit.GetMinZExtent()+pVoxelLimit.GetMaxZExtent())*0.5);
00354         
00355       if (Inside(pTransform.Inverse().TransformPoint(clipCentre))!=kOutside)
00356       {
00357         existsAfterClip=true;
00358         pMin=pVoxelLimit.GetMinExtent(pAxis);
00359         pMax=pVoxelLimit.GetMaxExtent(pAxis);
00360       }
00361     }
00362     delete vertices;
00363     return existsAfterClip;
00364   }
00365 }
00366 
00368 //
00369 // Return whether point inside/outside/on surface, using tolerance
00370 
00371 EInside G4Trd::Inside( const G4ThreeVector& p ) const
00372 {  
00373   EInside in=kOutside;
00374   G4double x,y,zbase1,zbase2;
00375     
00376   if (std::fabs(p.z())<=fDz-kCarTolerance/2)
00377   {
00378     zbase1=p.z()+fDz;  // Dist from -ve z plane
00379     zbase2=fDz-p.z();  // Dist from +ve z plane
00380 
00381     // Check whether inside x tolerance
00382     //
00383     x=0.5*(fDx2*zbase1+fDx1*zbase2)/fDz - kCarTolerance/2;
00384     if (std::fabs(p.x())<=x)
00385     {
00386       y=0.5*((fDy2*zbase1+fDy1*zbase2))/fDz - kCarTolerance/2;
00387       if (std::fabs(p.y())<=y)
00388       {
00389         in=kInside;
00390       }
00391       else if (std::fabs(p.y())<=y+kCarTolerance)
00392       {
00393         in=kSurface;
00394       }
00395     }
00396     else if (std::fabs(p.x())<=x+kCarTolerance)
00397     {
00398       // y = y half width of shape at z of point + tolerant boundary
00399       //
00400       y=0.5*((fDy2*zbase1+fDy1*zbase2))/fDz + kCarTolerance/2;
00401       if (std::fabs(p.y())<=y)
00402       {
00403         in=kSurface;
00404       }
00405     }
00406   }
00407   else if (std::fabs(p.z())<=fDz+kCarTolerance/2)
00408   {
00409     // Only need to check outer tolerant boundaries
00410     //
00411     zbase1=p.z()+fDz;  // Dist from -ve z plane
00412     zbase2=fDz-p.z();   // Dist from +ve z plane
00413 
00414     // x = x half width of shape at z of point plus tolerance
00415     //
00416     x=0.5*(fDx2*zbase1+fDx1*zbase2)/fDz + kCarTolerance/2;
00417     if (std::fabs(p.x())<=x)
00418     {
00419       // y = y half width of shape at z of point
00420       //
00421       y=0.5*((fDy2*zbase1+fDy1*zbase2))/fDz + kCarTolerance/2;
00422       if (std::fabs(p.y())<=y) in=kSurface;
00423     }
00424   }  
00425   return in;
00426 }
00427 
00429 //
00430 // Calculate side nearest to p, and return normal
00431 // If two sides are equidistant, normal of first side (x/y/z) 
00432 // encountered returned
00433 
00434 G4ThreeVector G4Trd::SurfaceNormal( const G4ThreeVector& p ) const
00435 {
00436   G4ThreeVector norm, sumnorm(0.,0.,0.);
00437   G4int noSurfaces = 0; 
00438   G4double z = 2.0*fDz, tanx, secx, newpx, widx;
00439   G4double tany, secy, newpy, widy;
00440   G4double distx, disty, distz, fcos;
00441   G4double delta = 0.5*kCarTolerance;
00442 
00443   tanx  = (fDx2 - fDx1)/z;
00444   secx  = std::sqrt(1.0+tanx*tanx);
00445   newpx = std::fabs(p.x())-p.z()*tanx;
00446   widx  = fDx2 - fDz*tanx;
00447 
00448   tany  = (fDy2 - fDy1)/z;
00449   secy  = std::sqrt(1.0+tany*tany);
00450   newpy = std::fabs(p.y())-p.z()*tany;
00451   widy  = fDy2 - fDz*tany;
00452 
00453   distx = std::fabs(newpx-widx)/secx;       // perp. distance to x side
00454   disty = std::fabs(newpy-widy)/secy;       //                to y side
00455   distz = std::fabs(std::fabs(p.z())-fDz);  //                to z side
00456 
00457   fcos              = 1.0/secx;
00458   G4ThreeVector nX  = G4ThreeVector( fcos,0,-tanx*fcos);
00459   G4ThreeVector nmX = G4ThreeVector(-fcos,0,-tanx*fcos);
00460 
00461   fcos              = 1.0/secy;
00462   G4ThreeVector nY  = G4ThreeVector(0, fcos,-tany*fcos);
00463   G4ThreeVector nmY = G4ThreeVector(0,-fcos,-tany*fcos);
00464   G4ThreeVector nZ  = G4ThreeVector( 0, 0,  1.0);
00465  
00466   if (distx <= delta)      
00467   {
00468     noSurfaces ++;
00469     if ( p.x() >= 0.) sumnorm += nX;
00470     else              sumnorm += nmX;   
00471   }
00472   if (disty <= delta)
00473   {
00474     noSurfaces ++;
00475     if ( p.y() >= 0.) sumnorm += nY;
00476     else              sumnorm += nmY;   
00477   }
00478   if (distz <= delta)  
00479   {
00480     noSurfaces ++;
00481     if ( p.z() >= 0.) sumnorm += nZ;
00482     else              sumnorm -= nZ; 
00483   }
00484   if ( noSurfaces == 0 )
00485   {
00486 #ifdef G4CSGDEBUG
00487     G4Exception("G4Trd::SurfaceNormal(p)", "GeomSolids1002", JustWarning, 
00488                 "Point p is not on surface !?" );
00489 #endif 
00490      norm = ApproxSurfaceNormal(p);
00491   }
00492   else if ( noSurfaces == 1 ) norm = sumnorm;
00493   else                        norm = sumnorm.unit();
00494   return norm;   
00495 }
00496 
00497 
00499 //
00500 // Algorithm for SurfaceNormal() following the original specification
00501 // for points not on the surface
00502 
00503 G4ThreeVector G4Trd::ApproxSurfaceNormal( const G4ThreeVector& p ) const
00504 {
00505   G4ThreeVector norm;
00506   G4double z,tanx,secx,newpx,widx;
00507   G4double tany,secy,newpy,widy;
00508   G4double distx,disty,distz,fcos;
00509 
00510   z=2.0*fDz;
00511 
00512   tanx=(fDx2-fDx1)/z;
00513   secx=std::sqrt(1.0+tanx*tanx);
00514   newpx=std::fabs(p.x())-p.z()*tanx;
00515   widx=fDx2-fDz*tanx;
00516 
00517   tany=(fDy2-fDy1)/z;
00518   secy=std::sqrt(1.0+tany*tany);
00519   newpy=std::fabs(p.y())-p.z()*tany;
00520   widy=fDy2-fDz*tany;
00521 
00522   distx=std::fabs(newpx-widx)/secx;  // perpendicular distance to x side
00523   disty=std::fabs(newpy-widy)/secy;  //                        to y side
00524   distz=std::fabs(std::fabs(p.z())-fDz);  //                        to z side
00525 
00526   // find closest side
00527   //
00528   if (distx<=disty)
00529   { 
00530     if (distx<=distz) 
00531     {
00532       // Closest to X
00533       //
00534       fcos=1.0/secx;
00535       // normal=(+/-std::cos(ang),0,-std::sin(ang))
00536       if (p.x()>=0)
00537         norm=G4ThreeVector(fcos,0,-tanx*fcos);
00538       else
00539         norm=G4ThreeVector(-fcos,0,-tanx*fcos);
00540     }
00541     else
00542     {
00543       // Closest to Z
00544       //
00545       if (p.z()>=0)
00546         norm=G4ThreeVector(0,0,1);
00547       else
00548         norm=G4ThreeVector(0,0,-1);
00549     }
00550   }
00551   else
00552   {  
00553     if (disty<=distz)
00554     {
00555       // Closest to Y
00556       //
00557       fcos=1.0/secy;
00558       if (p.y()>=0)
00559         norm=G4ThreeVector(0,fcos,-tany*fcos);
00560       else
00561         norm=G4ThreeVector(0,-fcos,-tany*fcos);
00562     }
00563     else 
00564     {
00565       // Closest to Z
00566       //
00567       if (p.z()>=0)
00568         norm=G4ThreeVector(0,0,1);
00569       else
00570         norm=G4ThreeVector(0,0,-1);
00571     }
00572   }
00573   return norm;   
00574 }
00575 
00577 //
00578 // Calculate distance to shape from outside
00579 // - return kInfinity if no intersection
00580 //
00581 // ALGORITHM:
00582 // For each component, calculate pair of minimum and maximum intersection
00583 // values for which the particle is in the extent of the shape
00584 // - The smallest (MAX minimum) allowed distance of the pairs is intersect
00585 // - Z plane intersectin uses tolerance
00586 // - XZ YZ planes use logic & *SLIGHTLY INCORRECT* tolerance
00587 //   (this saves at least 1 sqrt, 1 multiply and 1 divide... in applicable
00588 //    cases)
00589 // - Note: XZ and YZ planes each divide space into four regions,
00590 //   characterised by ss1 ss2
00591 // NOTE:
00592 //
00593 // `Inside' safe - meaningful answers given if point is inside the exact
00594 // shape.
00595 
00596 G4double G4Trd::DistanceToIn( const G4ThreeVector& p,
00597                               const G4ThreeVector& v ) const
00598 {  
00599   G4double snxt = kInfinity ;    // snxt = default return value
00600   G4double smin,smax;
00601   G4double s1,s2,tanxz,tanyz,ds1,ds2;
00602   G4double ss1,ss2,sn1=0.,sn2=0.,Dist;
00603 
00604   if ( v.z() )  // Calculate valid z intersect range
00605   {
00606     if ( v.z() > 0 )   // Calculate smax: must be +ve or no intersection.
00607     {
00608       Dist = fDz - p.z() ;  // to plane at +dz
00609 
00610       if (Dist >= 0.5*kCarTolerance)
00611       {
00612         smax = Dist/v.z() ;
00613         smin = -(fDz + p.z())/v.z() ;
00614       }
00615       else  return snxt ;
00616     }
00617     else // v.z <0
00618     {
00619       Dist=fDz+p.z();  // plane at -dz
00620 
00621       if ( Dist >= 0.5*kCarTolerance )
00622       {
00623         smax = -Dist/v.z() ;
00624         smin = (fDz - p.z())/v.z() ;
00625       }
00626       else return snxt ; 
00627     }
00628     if (smin < 0 ) smin = 0 ;
00629   }
00630   else // v.z=0
00631   {
00632     if (std::fabs(p.z()) >= fDz ) return snxt ;     // Outside & no intersect
00633     else
00634     {
00635       smin = 0 ;    // Always inside z range
00636       smax = kInfinity;
00637     }
00638   }
00639 
00640   // Calculate x intersection range
00641   //
00642   // Calc half width at p.z, and components towards planes
00643 
00644   tanxz = (fDx2 - fDx1)*0.5/fDz ;
00645   s1    = 0.5*(fDx1+fDx2) + tanxz*p.z() ;  // x half width at p.z
00646   ds1   = v.x() - tanxz*v.z() ;       // Components of v towards faces at +-x
00647   ds2   = v.x() + tanxz*v.z() ;
00648   ss1   = s1 - p.x() ;         // -delta x to +ve plane
00649                                // -ve when outside
00650   ss2   = -s1 - p.x() ;        // -delta x to -ve plane
00651                                // +ve when outside
00652     
00653   if (ss1 < 0 && ss2 <= 0 )
00654   {
00655     if (ds1 < 0)   // In +ve coord Area
00656     {
00657       sn1 = ss1/ds1 ;
00658 
00659       if ( ds2 < 0 ) sn2 = ss2/ds2 ;           
00660       else           sn2 = kInfinity ;
00661     }
00662     else return snxt ;
00663   }
00664   else if ( ss1 >= 0 && ss2 > 0 )
00665   {
00666     if ( ds2 > 0 )  // In -ve coord Area
00667     {
00668       sn1 = ss2/ds2 ;
00669 
00670       if (ds1 > 0)  sn2 = ss1/ds1 ;      
00671       else          sn2 = kInfinity;      
00672         
00673     }
00674     else   return snxt ;
00675   }
00676   else if (ss1 >= 0 && ss2 <= 0 )
00677   {
00678     // Inside Area - calculate leaving distance
00679     // *Don't* use exact distance to side for tolerance
00680     //                                             = ss1*std::cos(ang xz)
00681     //                                             = ss1/std::sqrt(1.0+tanxz*tanxz)
00682     sn1 = 0 ;
00683 
00684     if ( ds1 > 0 )
00685     {
00686       if (ss1 > 0.5*kCarTolerance) sn2 = ss1/ds1 ; // Leave +ve side extent
00687       else                         return snxt ;   // Leave immediately by +ve 
00688     }
00689     else  sn2 = kInfinity ;
00690       
00691     if ( ds2 < 0 )
00692     {
00693       if ( ss2 < -0.5*kCarTolerance )
00694       {
00695         Dist = ss2/ds2 ;            // Leave -ve side extent
00696         if ( Dist < sn2 ) sn2 = Dist ;
00697       }
00698       else  return snxt ;
00699     }    
00700   }
00701   else if (ss1 < 0 && ss2 > 0 )
00702   {
00703     // Within +/- plane cross-over areas (not on boundaries ss1||ss2==0)
00704 
00705     if ( ds1 >= 0 || ds2 <= 0 )
00706     {   
00707       return snxt ;
00708     }
00709     else  // Will intersect & stay inside
00710     {
00711       sn1  = ss1/ds1 ;
00712       Dist = ss2/ds2 ;
00713       if (Dist > sn1 ) sn1 = Dist ;
00714       sn2 = kInfinity ;
00715     }
00716   }
00717 
00718   // Reduce allowed range of distances as appropriate
00719 
00720   if ( sn1 > smin ) smin = sn1 ;
00721   if ( sn2 < smax ) smax = sn2 ;
00722 
00723   // Check for incompatible ranges (eg z intersects between 50 ->100 and x
00724   // only 10-40 -> no intersection)
00725 
00726   if ( smax < smin ) return snxt ;
00727 
00728   // Calculate valid y intersection range 
00729   // (repeat of x intersection code)
00730 
00731   tanyz = (fDy2-fDy1)*0.5/fDz ;
00732   s2    = 0.5*(fDy1+fDy2) + tanyz*p.z() ;  // y half width at p.z
00733   ds1   = v.y() - tanyz*v.z() ;       // Components of v towards faces at +-y
00734   ds2   = v.y() + tanyz*v.z() ;
00735   ss1   = s2 - p.y() ;         // -delta y to +ve plane
00736   ss2   = -s2 - p.y() ;        // -delta y to -ve plane
00737     
00738   if ( ss1 < 0 && ss2 <= 0 )
00739   {
00740     if (ds1 < 0 ) // In +ve coord Area
00741     {
00742       sn1 = ss1/ds1 ;
00743       if ( ds2 < 0 )  sn2 = ss2/ds2 ;
00744       else            sn2 = kInfinity ;
00745     }
00746     else   return snxt ;
00747   }
00748   else if ( ss1 >= 0 && ss2 > 0 )
00749   {
00750     if ( ds2 > 0 )  // In -ve coord Area
00751     {
00752       sn1 = ss2/ds2 ;
00753       if ( ds1 > 0 )  sn2 = ss1/ds1 ;
00754       else            sn2 = kInfinity ;      
00755     }
00756     else   return snxt ;
00757   }
00758   else if (ss1 >= 0 && ss2 <= 0 )
00759   {
00760     // Inside Area - calculate leaving distance
00761     // *Don't* use exact distance to side for tolerance
00762     //                                          = ss1*std::cos(ang yz)
00763     //                                          = ss1/std::sqrt(1.0+tanyz*tanyz)
00764     sn1 = 0 ;
00765 
00766     if ( ds1 > 0 )
00767     {
00768       if (ss1 > 0.5*kCarTolerance) sn2 = ss1/ds1 ; // Leave +ve side extent
00769       else                         return snxt ;   // Leave immediately by +ve
00770     }
00771     else  sn2 = kInfinity ;
00772       
00773     if ( ds2 < 0 )
00774     {
00775       if ( ss2 < -0.5*kCarTolerance )
00776       {
00777         Dist = ss2/ds2 ; // Leave -ve side extent
00778         if (Dist < sn2) sn2=Dist;
00779       }
00780       else  return snxt ;
00781     }    
00782   }
00783   else if (ss1 < 0 && ss2 > 0 )
00784   {
00785     // Within +/- plane cross-over areas (not on boundaries ss1||ss2==0)
00786 
00787     if (ds1 >= 0 || ds2 <= 0 )  
00788     {
00789       return snxt ;
00790     }
00791     else  // Will intersect & stay inside
00792     {
00793       sn1 = ss1/ds1 ;
00794       Dist = ss2/ds2 ;
00795       if (Dist > sn1 ) sn1 = Dist ;
00796       sn2 = kInfinity ;
00797     }
00798   }
00799   
00800   // Reduce allowed range of distances as appropriate
00801 
00802   if ( sn1 > smin) smin = sn1 ;
00803   if ( sn2 < smax) smax = sn2 ;
00804 
00805   // Check for incompatible ranges (eg x intersects between 50 ->100 and y
00806   // only 10-40 -> no intersection). Set snxt if ok
00807 
00808   if ( smax > smin ) snxt = smin ;
00809   if (snxt < 0.5*kCarTolerance ) snxt = 0.0 ;
00810 
00811   return snxt ;
00812 }
00813 
00815 //
00816 // Approximate distance to shape
00817 // Calculate perpendicular distances to z/x/y surfaces, return largest
00818 // which is the most fast estimation of shortest distance to Trd
00819 //  - Safe underestimate
00820 //  - If point within exact shape, return 0 
00821 
00822 G4double G4Trd::DistanceToIn( const G4ThreeVector& p ) const
00823 {
00824   G4double safe=0.0;
00825   G4double tanxz,distx,safx;
00826   G4double tanyz,disty,safy;
00827   G4double zbase;
00828 
00829   safe=std::fabs(p.z())-fDz;
00830   if (safe<0) safe=0;      // Also used to ensure x/y distances
00831                            // POSITIVE 
00832 
00833   zbase=fDz+p.z();
00834 
00835   // Find distance along x direction to closest x plane
00836   //
00837   tanxz=(fDx2-fDx1)*0.5/fDz;
00838   //    widx=fDx1+tanxz*(fDz+p.z()); // x width at p.z
00839   //    distx=std::fabs(p.x())-widx;      // distance to plane
00840   distx=std::fabs(p.x())-(fDx1+tanxz*zbase);
00841   if (distx>safe)
00842   {
00843     safx=distx/std::sqrt(1.0+tanxz*tanxz); // vector Dist=Dist*std::cos(ang)
00844     if (safx>safe) safe=safx;
00845   }
00846 
00847   // Find distance along y direction to slanted wall
00848   tanyz=(fDy2-fDy1)*0.5/fDz;
00849   //    widy=fDy1+tanyz*(fDz+p.z()); // y width at p.z
00850   //    disty=std::fabs(p.y())-widy;      // distance to plane
00851   disty=std::fabs(p.y())-(fDy1+tanyz*zbase);
00852   if (disty>safe)    
00853   {
00854     safy=disty/std::sqrt(1.0+tanyz*tanyz); // distance along vector
00855     if (safy>safe) safe=safy;
00856   }
00857   return safe;
00858 }
00859 
00861 //
00862 // Calcluate distance to surface of shape from inside
00863 // Calculate distance to x/y/z planes - smallest is exiting distance
00864 // - z planes have std. check for tolerance
00865 // - xz yz planes have check based on distance || to x or y axis
00866 //   (not corrected for slope of planes)
00867 // ?BUG? If v.z==0 are there cases when snside not set????
00868 
00869 G4double G4Trd::DistanceToOut( const G4ThreeVector& p,
00870                                const G4ThreeVector& v,
00871                                const G4bool calcNorm,
00872                                      G4bool *validNorm,
00873                                      G4ThreeVector *n ) const
00874 {
00875   ESide side = kUndefined, snside = kUndefined;
00876   G4double snxt,pdist;
00877   G4double central,ss1,ss2,ds1,ds2,sn=0.,sn2=0.;
00878   G4double tanxz=0.,cosxz=0.,tanyz=0.,cosyz=0.;
00879 
00880   if (calcNorm) *validNorm=true; // All normals are valid
00881 
00882   // Calculate z plane intersection
00883   if (v.z()>0)
00884   {
00885     pdist=fDz-p.z();
00886     if (pdist>kCarTolerance/2)
00887     {
00888       snxt=pdist/v.z();
00889       side=kPZ;
00890     }
00891     else
00892     {
00893       if (calcNorm)
00894       {
00895         *n=G4ThreeVector(0,0,1);
00896       }
00897       return snxt=0;
00898     }
00899   }
00900   else if (v.z()<0) 
00901   {
00902     pdist=fDz+p.z();
00903     if (pdist>kCarTolerance/2)
00904     {
00905       snxt=-pdist/v.z();
00906       side=kMZ;
00907     }
00908     else
00909     {
00910       if (calcNorm)
00911       {
00912         *n=G4ThreeVector(0,0,-1);
00913       }
00914       return snxt=0;
00915     }
00916   }
00917   else
00918   {
00919     snxt=kInfinity;
00920   }
00921 
00922   //
00923   // Calculate x intersection
00924   //
00925   tanxz=(fDx2-fDx1)*0.5/fDz;
00926   central=0.5*(fDx1+fDx2);
00927 
00928   // +ve plane (1)
00929   //
00930   ss1=central+tanxz*p.z()-p.x();  // distance || x axis to plane
00931                                   // (+ve if point inside)
00932   ds1=v.x()-tanxz*v.z();    // component towards plane at +x
00933                             // (-ve if +ve -> -ve direction)
00934   // -ve plane (2)
00935   //
00936   ss2=-tanxz*p.z()-p.x()-central;  //distance || x axis to plane
00937                                    // (-ve if point inside)
00938   ds2=tanxz*v.z()+v.x();    // component towards plane at -x
00939 
00940   if (ss1>0&&ss2<0)
00941   {
00942     // Normal case - entirely inside region
00943     if (ds1<=0&&ds2<0)
00944     {   
00945       if (ss2<-kCarTolerance/2)
00946       {
00947         sn=ss2/ds2;  // Leave by -ve side
00948         snside=kMX;
00949       }
00950       else
00951       {
00952         sn=0; // Leave immediately by -ve side
00953         snside=kMX;
00954       }
00955     }
00956     else if (ds1>0&&ds2>=0)
00957     {
00958       if (ss1>kCarTolerance/2)
00959       {
00960         sn=ss1/ds1;  // Leave by +ve side
00961         snside=kPX;
00962       }
00963       else
00964       {
00965         sn=0; // Leave immediately by +ve side
00966         snside=kPX;
00967       }
00968     }
00969     else if (ds1>0&&ds2<0)
00970     {
00971       if (ss1>kCarTolerance/2)
00972       {
00973         // sn=ss1/ds1;  // Leave by +ve side
00974         if (ss2<-kCarTolerance/2)
00975         {
00976           sn=ss1/ds1;  // Leave by +ve side
00977           sn2=ss2/ds2;
00978           if (sn2<sn)
00979           {
00980             sn=sn2;
00981             snside=kMX;
00982           }
00983           else
00984           {
00985             snside=kPX;
00986           }
00987         }
00988         else
00989         {
00990           sn=0; // Leave immediately by -ve
00991           snside=kMX;
00992         }      
00993       }
00994       else
00995       {
00996         sn=0; // Leave immediately by +ve side
00997         snside=kPX;
00998       }
00999     }
01000     else
01001     {
01002       // Must be || to both
01003       //
01004       sn=kInfinity;    // Don't leave by either side
01005     }
01006   }
01007   else if (ss1<=0&&ss2<0)
01008   {
01009     // Outside, in +ve Area
01010     
01011     if (ds1>0)
01012     {
01013       sn=0;       // Away from shape
01014                   // Left by +ve side
01015       snside=kPX;
01016     }
01017     else
01018     {
01019       if (ds2<0)
01020       {
01021         // Ignore +ve plane and use -ve plane intersect
01022         //
01023         sn=ss2/ds2; // Leave by -ve side
01024         snside=kMX;
01025       }
01026       else
01027       {
01028         // Must be || to both -> exit determined by other axes
01029         //
01030         sn=kInfinity; // Don't leave by either side
01031       }
01032     }
01033   }
01034   else if (ss1>0&&ss2>=0)
01035   {
01036     // Outside, in -ve Area
01037 
01038     if (ds2<0)
01039     {
01040       sn=0;       // away from shape
01041                   // Left by -ve side
01042       snside=kMX;
01043     }
01044     else
01045     {
01046       if (ds1>0)
01047       {
01048         // Ignore +ve plane and use -ve plane intersect
01049         //
01050         sn=ss1/ds1; // Leave by +ve side
01051         snside=kPX;
01052       }
01053       else
01054       {
01055         // Must be || to both -> exit determined by other axes
01056         //
01057         sn=kInfinity; // Don't leave by either side
01058       }
01059     }
01060   }
01061 
01062   // Update minimum exit distance
01063 
01064   if (sn<snxt)
01065   {
01066     snxt=sn;
01067     side=snside;
01068   }
01069   if (snxt>0)
01070   {
01071     // Calculate y intersection
01072 
01073     tanyz=(fDy2-fDy1)*0.5/fDz;
01074     central=0.5*(fDy1+fDy2);
01075 
01076     // +ve plane (1)
01077     //
01078     ss1=central+tanyz*p.z()-p.y(); // distance || y axis to plane
01079                                    // (+ve if point inside)
01080     ds1=v.y()-tanyz*v.z();  // component towards +ve plane
01081                             // (-ve if +ve -> -ve direction)
01082     // -ve plane (2)
01083     //
01084     ss2=-tanyz*p.z()-p.y()-central; // distance || y axis to plane
01085                                     // (-ve if point inside)
01086     ds2=tanyz*v.z()+v.y();  // component towards -ve plane
01087 
01088     if (ss1>0&&ss2<0)
01089     {
01090       // Normal case - entirely inside region
01091 
01092       if (ds1<=0&&ds2<0)
01093       {   
01094         if (ss2<-kCarTolerance/2)
01095         {
01096           sn=ss2/ds2;  // Leave by -ve side
01097           snside=kMY;
01098         }
01099         else
01100         {
01101           sn=0; // Leave immediately by -ve side
01102           snside=kMY;
01103         }
01104       }
01105       else if (ds1>0&&ds2>=0)
01106       {
01107         if (ss1>kCarTolerance/2)
01108         {
01109           sn=ss1/ds1;  // Leave by +ve side
01110           snside=kPY;
01111         }
01112         else
01113         {
01114           sn=0; // Leave immediately by +ve side
01115           snside=kPY;
01116         }
01117       }
01118       else if (ds1>0&&ds2<0)
01119       {
01120         if (ss1>kCarTolerance/2)
01121         {
01122           // sn=ss1/ds1;  // Leave by +ve side
01123           if (ss2<-kCarTolerance/2)
01124           {
01125             sn=ss1/ds1;  // Leave by +ve side
01126             sn2=ss2/ds2;
01127             if (sn2<sn)
01128             {
01129               sn=sn2;
01130               snside=kMY;
01131             }
01132             else
01133             {
01134               snside=kPY;
01135             }
01136           }
01137           else
01138           {
01139             sn=0; // Leave immediately by -ve
01140             snside=kMY;
01141           }
01142         }
01143         else
01144         {
01145           sn=0; // Leave immediately by +ve side
01146           snside=kPY;
01147         }
01148       }
01149       else
01150       {
01151         // Must be || to both
01152         //
01153         sn=kInfinity;    // Don't leave by either side
01154       }
01155     }
01156     else if (ss1<=0&&ss2<0)
01157     {
01158       // Outside, in +ve Area
01159 
01160       if (ds1>0)
01161       {
01162         sn=0;       // Away from shape
01163                     // Left by +ve side
01164         snside=kPY;
01165       }
01166       else
01167       {
01168         if (ds2<0)
01169         {
01170           // Ignore +ve plane and use -ve plane intersect
01171           //
01172           sn=ss2/ds2; // Leave by -ve side
01173           snside=kMY;
01174         }
01175         else
01176         {
01177           // Must be || to both -> exit determined by other axes
01178           //
01179           sn=kInfinity; // Don't leave by either side
01180         }
01181       }
01182     }
01183     else if (ss1>0&&ss2>=0)
01184     {
01185       // Outside, in -ve Area
01186       if (ds2<0)
01187       {
01188         sn=0;       // away from shape
01189                     // Left by -ve side
01190         snside=kMY;
01191       }
01192       else
01193       {
01194         if (ds1>0)
01195         {
01196           // Ignore +ve plane and use -ve plane intersect
01197           //
01198           sn=ss1/ds1; // Leave by +ve side
01199           snside=kPY;
01200         }
01201         else
01202         {
01203           // Must be || to both -> exit determined by other axes
01204           //
01205           sn=kInfinity; // Don't leave by either side
01206         }
01207       }
01208     }
01209 
01210     // Update minimum exit distance
01211 
01212     if (sn<snxt)
01213     {
01214       snxt=sn;
01215       side=snside;
01216     }
01217   }
01218 
01219   if (calcNorm)
01220   {
01221     switch (side)
01222     {
01223       case kPX:
01224         cosxz=1.0/std::sqrt(1.0+tanxz*tanxz);
01225         *n=G4ThreeVector(cosxz,0,-tanxz*cosxz);
01226         break;
01227       case kMX:
01228         cosxz=-1.0/std::sqrt(1.0+tanxz*tanxz);
01229         *n=G4ThreeVector(cosxz,0,tanxz*cosxz);
01230         break;
01231       case kPY:
01232         cosyz=1.0/std::sqrt(1.0+tanyz*tanyz);
01233         *n=G4ThreeVector(0,cosyz,-tanyz*cosyz);
01234         break;
01235       case kMY:
01236         cosyz=-1.0/std::sqrt(1.0+tanyz*tanyz);
01237         *n=G4ThreeVector(0,cosyz,tanyz*cosyz);
01238         break;
01239       case kPZ:
01240         *n=G4ThreeVector(0,0,1);
01241         break;
01242       case kMZ:
01243         *n=G4ThreeVector(0,0,-1);
01244         break;
01245       default:
01246         DumpInfo();
01247         G4Exception("G4Trd::DistanceToOut(p,v,..)",
01248                     "GeomSolids1002", JustWarning, 
01249                     "Undefined side for valid surface normal to solid.");
01250         break;
01251     }
01252   }
01253   return snxt; 
01254 }
01255 
01257 //
01258 // Calculate exact shortest distance to any boundary from inside
01259 // - Returns 0 is point outside
01260 
01261 G4double G4Trd::DistanceToOut( const G4ThreeVector& p ) const
01262 {
01263   G4double safe=0.0;
01264   G4double tanxz,xdist,saf1;
01265   G4double tanyz,ydist,saf2;
01266   G4double zbase;
01267 
01268 #ifdef G4CSGDEBUG
01269   if( Inside(p) == kOutside )
01270   {
01271      G4int oldprc = G4cout.precision(16) ;
01272      G4cout << G4endl ;
01273      DumpInfo();
01274      G4cout << "Position:"  << G4endl << G4endl ;
01275      G4cout << "p.x() = "   << p.x()/mm << " mm" << G4endl ;
01276      G4cout << "p.y() = "   << p.y()/mm << " mm" << G4endl ;
01277      G4cout << "p.z() = "   << p.z()/mm << " mm" << G4endl << G4endl ;
01278      G4cout.precision(oldprc) ;
01279      G4Exception("G4Trd::DistanceToOut(p)", "GeomSolids1002", JustWarning, 
01280                  "Point p is outside !?" );
01281   }
01282 #endif
01283 
01284   safe=fDz-std::fabs(p.z());  // z perpendicular Dist
01285 
01286   zbase=fDz+p.z();
01287 
01288   // xdist = distance perpendicular to z axis to closest x plane from p
01289   //       = (x half width of shape at p.z) - std::fabs(p.x)
01290   //
01291   tanxz=(fDx2-fDx1)*0.5/fDz;
01292   xdist=fDx1+tanxz*zbase-std::fabs(p.x());
01293   saf1=xdist/std::sqrt(1.0+tanxz*tanxz); // x*std::cos(ang_xz) =
01294                                     // shortest (perpendicular)
01295                                     // distance to plane
01296   tanyz=(fDy2-fDy1)*0.5/fDz;
01297   ydist=fDy1+tanyz*zbase-std::fabs(p.y());
01298   saf2=ydist/std::sqrt(1.0+tanyz*tanyz);
01299 
01300   // Return minimum x/y/z distance
01301   //
01302   if (safe>saf1) safe=saf1;
01303   if (safe>saf2) safe=saf2;
01304 
01305   if (safe<0) safe=0;
01306   return safe;     
01307 }
01308 
01310 //
01311 // Create a List containing the transformed vertices
01312 // Ordering [0-3] -fDz cross section
01313 //          [4-7] +fDz cross section such that [0] is below [4],
01314 //                                             [1] below [5] etc.
01315 // Note:
01316 //  Caller has deletion resposibility
01317 
01318 G4ThreeVectorList*
01319 G4Trd::CreateRotatedVertices( const G4AffineTransform& pTransform ) const
01320 {
01321   G4ThreeVectorList *vertices;
01322   vertices=new G4ThreeVectorList();
01323   if (vertices)
01324   {
01325     vertices->reserve(8);
01326     G4ThreeVector vertex0(-fDx1,-fDy1,-fDz);
01327     G4ThreeVector vertex1(fDx1,-fDy1,-fDz);
01328     G4ThreeVector vertex2(fDx1,fDy1,-fDz);
01329     G4ThreeVector vertex3(-fDx1,fDy1,-fDz);
01330     G4ThreeVector vertex4(-fDx2,-fDy2,fDz);
01331     G4ThreeVector vertex5(fDx2,-fDy2,fDz);
01332     G4ThreeVector vertex6(fDx2,fDy2,fDz);
01333     G4ThreeVector vertex7(-fDx2,fDy2,fDz);
01334 
01335     vertices->push_back(pTransform.TransformPoint(vertex0));
01336     vertices->push_back(pTransform.TransformPoint(vertex1));
01337     vertices->push_back(pTransform.TransformPoint(vertex2));
01338     vertices->push_back(pTransform.TransformPoint(vertex3));
01339     vertices->push_back(pTransform.TransformPoint(vertex4));
01340     vertices->push_back(pTransform.TransformPoint(vertex5));
01341     vertices->push_back(pTransform.TransformPoint(vertex6));
01342     vertices->push_back(pTransform.TransformPoint(vertex7));
01343   }
01344   else
01345   {
01346     DumpInfo();
01347     G4Exception("G4Trd::CreateRotatedVertices()",
01348                 "GeomSolids0003", FatalException,
01349                 "Error in allocation of vertices. Out of memory !");
01350   }
01351   return vertices;
01352 }
01353 
01355 //
01356 // GetEntityType
01357 
01358 G4GeometryType G4Trd::GetEntityType() const
01359 {
01360   return G4String("G4Trd");
01361 }
01362 
01364 //
01365 // Make a clone of the object
01366 //
01367 G4VSolid* G4Trd::Clone() const
01368 {
01369   return new G4Trd(*this);
01370 }
01371 
01373 //
01374 // Stream object contents to an output stream
01375 
01376 std::ostream& G4Trd::StreamInfo( std::ostream& os ) const
01377 {
01378   G4int oldprc = os.precision(16);
01379   os << "-----------------------------------------------------------\n"
01380      << "    *** Dump for solid - " << GetName() << " ***\n"
01381      << "    ===================================================\n"
01382      << " Solid type: G4Trd\n"
01383      << " Parameters: \n"
01384      << "    half length X, surface -dZ: " << fDx1/mm << " mm \n"
01385      << "    half length X, surface +dZ: " << fDx2/mm << " mm \n"
01386      << "    half length Y, surface -dZ: " << fDy1/mm << " mm \n"
01387      << "    half length Y, surface +dZ: " << fDy2/mm << " mm \n"
01388      << "    half length Z             : " << fDz/mm << " mm \n"
01389      << "-----------------------------------------------------------\n";
01390   os.precision(oldprc);
01391 
01392   return os;
01393 }
01394 
01395 
01397 //
01398 // GetPointOnSurface
01399 //
01400 // Return a point (G4ThreeVector) randomly and uniformly
01401 // selected on the solid surface
01402 
01403 G4ThreeVector G4Trd::GetPointOnSurface() const
01404 {
01405   G4double px, py, pz, tgX, tgY, secX, secY, select, sumS, tmp;
01406   G4double Sxy1, Sxy2, Sxy, Sxz, Syz;
01407 
01408   tgX  = 0.5*(fDx2-fDx1)/fDz;
01409   secX = std::sqrt(1+tgX*tgX);
01410   tgY  = 0.5*(fDy2-fDy1)/fDz;
01411   secY = std::sqrt(1+tgY*tgY);
01412 
01413   // calculate 0.25 of side surfaces, sumS is 0.25 of total surface
01414 
01415   Sxy1 = fDx1*fDy1; 
01416   Sxy2 = fDx2*fDy2;
01417   Sxy  = Sxy1 + Sxy2; 
01418   Sxz  = (fDx1 + fDx2)*fDz*secY; 
01419   Syz  = (fDy1 + fDy2)*fDz*secX;
01420   sumS = Sxy + Sxz + Syz;
01421 
01422   select = sumS*G4UniformRand();
01423  
01424   if( select < Sxy )                  // Sxy1 or Sxy2
01425   {
01426     if( select < Sxy1 ) 
01427     {
01428       pz = -fDz;
01429       px = -fDx1 + 2*fDx1*G4UniformRand();
01430       py = -fDy1 + 2*fDy1*G4UniformRand();
01431     }
01432     else      
01433     {
01434       pz =  fDz;
01435       px = -fDx2 + 2*fDx2*G4UniformRand();
01436       py = -fDy2 + 2*fDy2*G4UniformRand();
01437     }
01438   }
01439   else if ( ( select - Sxy ) < Sxz )    // Sxz
01440   {
01441     pz  = -fDz  + 2*fDz*G4UniformRand();
01442     tmp =  fDx1 + (pz + fDz)*tgX;
01443     px  = -tmp  + 2*tmp*G4UniformRand();
01444     tmp =  fDy1 + (pz + fDz)*tgY;
01445 
01446     if(G4UniformRand() > 0.5) { py =  tmp; }
01447     else                      { py = -tmp; }
01448   }
01449   else                                   // Syz
01450   {
01451     pz  = -fDz  + 2*fDz*G4UniformRand();
01452     tmp =  fDy1 + (pz + fDz)*tgY;
01453     py  = -tmp  + 2*tmp*G4UniformRand();
01454     tmp =  fDx1 + (pz + fDz)*tgX;
01455 
01456     if(G4UniformRand() > 0.5) { px =  tmp; }
01457     else                      { px = -tmp; }
01458   } 
01459   return G4ThreeVector(px,py,pz);
01460 }
01461 
01463 //
01464 // Methods for visualisation
01465 
01466 void G4Trd::DescribeYourselfTo ( G4VGraphicsScene& scene ) const
01467 {
01468   scene.AddSolid (*this);
01469 }
01470 
01471 G4Polyhedron* G4Trd::CreatePolyhedron () const
01472 {
01473   return new G4PolyhedronTrd2 (fDx1, fDx2, fDy1, fDy2, fDz);
01474 }
01475 
01476 G4NURBS* G4Trd::CreateNURBS () const
01477 {
01478   //  return new G4NURBSbox (fDx, fDy, fDz);
01479   return 0;
01480 }

Generated on Mon May 27 17:50:03 2013 for Geant4 by  doxygen 1.4.7