G4Ray.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$
00028 //
00029 // ----------------------------------------------------------------------
00030 // GEANT 4 class source file
00031 //
00032 // G4Ray.cc
00033 //
00034 // ----------------------------------------------------------------------
00035 
00036 #include "G4Ray.hh"
00037 #include "G4PointRat.hh"
00038 
00039 G4Ray::G4Ray()
00040   : r_min(0.), r_max(0.)
00041 {
00042 }
00043 
00044 G4Ray::G4Ray(const G4Point3D& start0, const G4Vector3D& dir0)
00045 {
00046   Init(start0, dir0);
00047 }
00048 
00049 G4Ray::~G4Ray()
00050 {
00051 }
00052 
00053 
00054 const G4Plane& G4Ray::GetPlane(G4int number_of_plane) const
00055 {
00056   if(number_of_plane==1)
00057     { return plane2; }
00058   else
00059     { return plane1; }
00060 }
00061 
00062 
00063 void G4Ray::CreatePlanes()
00064 {
00065   // Creates two orthogonal planes(plane1,plane2) the ray (rray) 
00066   // situated in the intersection of the planes. The planes are 
00067   // used to project the surface (nurb) in two dimensions.
00068   
00069   G4Vector3D RayDir    = dir;
00070   G4Point3D  RayOrigin = start;
00071 
00072   G4Point3D  p1, p2, p3, p4;
00073   G4Vector3D dir1, dir2;
00074   G4Vector3D invdir = G4Vector3D( PINFINITY );
00075   
00076   if(!NearZero(RayDir.x(), SQRT_SMALL_FASTF)) 
00077     { invdir.setX(1.0 / RayDir.x()); }
00078     
00079   if(!NearZero(RayDir.y(), SQRT_SMALL_FASTF)) 
00080     { invdir.setY(1.0 / RayDir.y()); }
00081     
00082   if(!NearZero(RayDir.z(), SQRT_SMALL_FASTF)) 
00083     { invdir.setZ(1.0 / RayDir.z()); }
00084 
00085   MatVecOrtho(dir1, RayDir);
00086   
00087   Vcross( dir2, RayDir, dir1);
00088   Vmove(p1, RayOrigin);
00089   Vadd2(p2, RayOrigin, RayDir);
00090   Vadd2(p3, RayOrigin, dir1);
00091   Vadd2(p4, RayOrigin, dir2);
00092 
00093   CalcPlane3Pts( plane1, p1, p3, p2);
00094   CalcPlane3Pts( plane2, p1, p2, p4);
00095 }
00096 
00097 
00098 void G4Ray::MatVecOrtho(register G4Vector3D &out,
00099                         register const G4Vector3D &in )
00100 {
00101   register G4double f;
00102   G4int             i_Which;
00103 
00104   if( NearZero(in.x(), 0.0001)
00105    && NearZero(in.y(), 0.0001)
00106    && NearZero(in.z(), 0.0001) )  
00107   {
00108     Vsetall( out, 0 );
00109     return;
00110   }
00111  
00112   //      Find component closest to zero 
00113   f = std::fabs(in.x());
00114   i_Which=0;
00115   
00116   if( std::fabs(in.y()) < f )
00117   {
00118     f = std::fabs(in.y());
00119     i_Which=1;
00120   }
00121   
00122   if( std::fabs(in.z()) < f )
00123   {
00124     i_Which=2;
00125   }
00126   
00127   if(!i_Which)
00128   {
00129     f = std::sqrt((in.y())*(in.y())+(in.z())*(in.z()));    // hypot(in.y(),in.z())
00130   }
00131   else
00132   {
00133     if(i_Which==1)
00134     {
00135       f = std::sqrt((in.z())*(in.z())+(in.x())*(in.x()));  // hypot(in.z(),in.x())
00136     }
00137     else
00138     {
00139       f = std::sqrt((in.x())*(in.x())+(in.y())*(in.y()));  // hypot(in.x(),in.y())
00140     }
00141   }
00142   if( NearZero( f, SMALL ) )
00143   {
00144     Vsetall( out, 0 );
00145     return;
00146   }
00147     
00148   f = 1.0/f;
00149     
00150   if(!i_Which)
00151   {
00152     out.setX(0.0);
00153     out.setY(-in.z()*f);
00154     out.setZ( in.y()*f);
00155   }
00156   else
00157   {
00158     if(i_Which==1)
00159     {
00160       out.setY(0.0);
00161       out.setZ(-in.x()*f);
00162       out.setX( in.y()*f);
00163     }
00164     else
00165     {
00166       out.setZ(0.0);
00167       out.setX(-in.z()*f);
00168       out.setY( in.y()*f);
00169     }
00170   }
00171 } 
00172 
00173 
00174 //                      CALC_PLANE_3PTS
00175 //
00176 //  Find the equation of a G4Plane that contains three points.
00177 //  Note that Normal vector created is expected to point out (see vmath.h),
00178 //  so the vector from A to C had better be counter-clockwise
00179 //  (about the point A) from the vector from A to B.
00180 //  This follows the outward-pointing Normal convention, and the
00181 //  right-hand rule for cross products.
00182 //
00183 /*
00184                         C
00185                         *
00186                         |\
00187                         | \
00188            ^     N      |  \
00189            |      \     |   \
00190            |       \    |    \
00191            |C-A     \   |     \
00192            |         \  |      \
00193            |          \ |       \
00194                        \|        \
00195                         *---------*
00196                         A         B
00197                            ----->
00198                             B-A
00199 */
00200 //  If the points are given in the order A B C (eg, *counter*-clockwise),
00201 //  then the outward pointing surface Normal N = (B-A) x (C-A).
00202 //
00203 //  Explicit Return -
00204 //       0      OK
00205 //      -1      Failure.  At least two of the points were not distinct,
00206 //              or all three were colinear.
00207 //
00208 //  Implicit Return -
00209 //      G4Plane   The G4Plane equation is stored here.
00210 
00211 
00212 G4int G4Ray::CalcPlane3Pts(G4Plane &plane1,
00213                            const G4Point3D& a,
00214                            const G4Point3D& b,
00215                            const G4Point3D& c )
00216 {
00217   // Creates the two orthogonal planes which are needed in projecting the
00218   // surface into 2D.
00219 
00220   G4Vector3D    B_A;
00221   G4Vector3D    C_A;
00222   G4Vector3D    C_B;
00223   
00224   register G4double mag;
00225 
00226   Vsub2( B_A, b, a );
00227   Vsub2( C_A, c, a );
00228   Vsub2( C_B, c, b );
00229 
00230   Vcross( plane1, B_A, C_A );
00231 
00232   //    Ensure unit length Normal 
00233   mag = Magnitude(plane1);
00234   if( mag  <= SQRT_SMALL_FASTF )
00235   {
00236     return(-1);//        FAIL 
00237   }
00238   
00239   mag = 1/mag;
00240   
00241   G4Plane pl2(plane1);
00242   Vscale( plane1, pl2, mag );
00243 
00244   //     Find distance from the origin to the G4Plane 
00245   plane1.d = Vdot( plane1, a );
00246   
00247   return(0);    //ok
00248 }
00249 
00250 
00251 void G4Ray::RayCheck()
00252 {
00253   // Check that the ray has a G4Vector3D...
00254   if (dir==G4Vector3D(0, 0, 0)) 
00255   {
00256     G4Exception("G4Ray::RayCheck()", "GeomSolids0002", FatalException,
00257                 "Invalid zero direction given !");
00258   }
00259 
00260   // Make sure that the vector is unit length
00261   dir= dir.unit();
00262   r_min = 0;
00263   r_max = 0;
00264 }
00265 
00266 
00267 void G4Ray::Vcross(G4Plane &a, const G4Vector3D &b, const G4Vector3D &c) 
00268 { 
00269   a.a = b.y()  * c.z()  - b.z()  * c.y() ;
00270   a.b = b.z()  * c.x()  - b.x()  * c.z() ;
00271   a.c = b.x()  * c.y()  - b.y()  * c.x() ;
00272 }
00273 
00274 
00275 void G4Ray::Vcross(G4Vector3D &a, const G4Vector3D &b, const G4Vector3D &c) 
00276 { 
00277   a.setX(b.y()  * c.z()  - b.z()  * c.y()) ;
00278   a.setY(b.z()  * c.x()  - b.x()  * c.z()) ;
00279   a.setZ(b.x()  * c.y()  - b.y()  * c.x()) ;
00280 }
00281 
00282 
00283 void G4Ray::Vmove(G4Point3D &a, const G4Point3D &b) 
00284 { 
00285   a.setX(b.x());
00286   a.setY(b.y());
00287   a.setZ(b.z());
00288 }
00289 
00290 
00291 void G4Ray::Vadd2(G4Point3D &a, const G4Point3D &b, const G4Vector3D &c) 
00292 {
00293   a.setX(b.x() + c.x()) ;
00294   a.setY(b.y() + c.y()) ;
00295   a.setZ(b.z() + c.z()) ;
00296 }       
00297   
00298 
00299 void G4Ray::Vsub2(G4Vector3D &a, const G4Point3D &b, const G4Point3D &c) 
00300 {
00301   a.setX(b.x() - c.x());
00302   a.setY(b.y() - c.y());
00303   a.setZ(b.z() - c.z());
00304 }
00305 
00306 
00307 void G4Ray::Vscale(G4Plane& a, const G4Plane& b, G4double c) 
00308 { 
00309   a.a = b.a * c;
00310   a.b = b.b * c;
00311   a.c = b.c * c;
00312 }
00313 
00314 
00315 G4double G4Ray::Vdot(const G4Plane &a, const G4Point3D &b) 
00316 {
00317   return (a.a * b.x() + 
00318           a.b * b.y() + 
00319           a.c * b.z());
00320 }
00321   
00322 
00323 G4double G4Ray::Magsq(const G4Plane &a) 
00324 {
00325   return ( a.a * a.a + a.b * a.b + a.c *a.c );
00326 }
00327   
00328 
00329 G4double G4Ray::Magnitude(const G4Plane &a) 
00330 {
00331   return (std::sqrt( Magsq( a )) );
00332 }

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