G4FPlane.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 // G4FPlane.cc
00033 //
00034 // ----------------------------------------------------------------------
00035 // Corrections by S.Giani:
00036 // - The constructor using iVec now properly stores both the internal and
00037 //   external boundaries in the bounds vector.
00038 // - Proper initialization of sameSense in both the constructors. 
00039 // - Addition of third argument (sense) in the second constructor to ensure
00040 //   consistent setting of the normal in all the client code.
00041 // - Proper use of the tolerance in the Intersect function.
00042 // ----------------------------------------------------------------------
00043 
00044 #include "G4FPlane.hh"
00045 #include "G4SystemOfUnits.hh"
00046 #include "G4CompositeCurve.hh"
00047 
00048 
00049 G4FPlane::G4FPlane( const G4Vector3D& direction,
00050                     const G4Vector3D& axis     , 
00051                     const G4Point3D&  Pt0, G4int sense )
00052   : pplace(direction, axis, Pt0), Convex(0), projectedBoundary(0)
00053 {
00054   G4Point3D Pt1 = G4Point3D( Pt0 + direction );
00055 
00056   // The plane include direction and axis is the normal,
00057   // so axis^direction is included in the plane
00058   G4Point3D Pt2 = G4Point3D( Pt0 + axis.cross(direction) );
00059 
00060   G4Ray::CalcPlane3Pts( Pl, Pt0, Pt1, Pt2 );
00061 
00062   active   = 1;
00063   sameSense = sense;
00064   CalcNormal();
00065   distance = kInfinity;
00066   Type     = 1;
00067 }
00068 
00069 
00070 G4FPlane::G4FPlane(const G4Point3DVector* pVec,
00071                    const G4Point3DVector* iVec,
00072                    G4int sense)
00073   : pplace( (*pVec)[0]-(*pVec)[1],                    // direction
00074             ((*pVec)[pVec->size()-1]-(*pVec)[0])
00075             .cross((*pVec)[0]-(*pVec)[1]),            // axis
00076             (*pVec)[0]                             )  // location
00077 
00078 {
00079   G4Ray::CalcPlane3Pts( Pl, (*pVec)[0], (*pVec)[1], (*pVec)[2] );
00080  
00081   G4CurveVector bounds;
00082   G4CompositeCurve* polygon;
00083 
00084   projectedBoundary = new G4SurfaceBoundary;
00085 
00086   sameSense = sense;
00087 
00088   // Outer boundary
00089 
00090   polygon= new G4CompositeCurve(*pVec);
00091  
00092   for (size_t i=0; i< polygon->GetSegments().size(); i++) 
00093     polygon->GetSegments()[i]->SetSameSense(sameSense);
00094 
00095   bounds.push_back(polygon);
00096   
00097   // Eventual inner boundary
00098   
00099   if (iVec) 
00100   {
00101     polygon= new G4CompositeCurve(*iVec);
00102 
00103     for (size_t i=0; i< polygon->GetSegments().size(); i++) 
00104     polygon->GetSegments()[i]->SetSameSense(sameSense);
00105 
00106     bounds.push_back(polygon);
00107   }
00108   
00109   // Set sense for boundaries  
00110   
00111   for (size_t j=0; j< bounds.size(); j++) 
00112     bounds[j]->SetSameSense(sameSense);
00113   
00114 
00115   SetBoundaries(&bounds);
00116       
00117   CalcNormal();
00118   IsConvex();
00119   distance = kInfinity;
00120   Type=1;
00121 }
00122 
00123 
00124 G4FPlane::~G4FPlane()
00125 {
00126   delete NormalX;
00127 }
00128 
00129 
00130 void G4FPlane::CalcBBox()
00131 {
00132   // This is needed since the bounds are used for the Solid
00133   // bbox calculation. The bbox test is NOT performed for
00134   // planar surfaces.
00135 
00136   // Finds the bounds of the G4Plane surface iow
00137   // calculates the bounds for a bounding box
00138   // to the surface. The bounding box is used
00139   // for a preliminary check of intersection.
00140   
00141   bbox= new G4BoundingBox3D(surfaceBoundary.BBox().GetBoxMin(), 
00142                             surfaceBoundary.BBox().GetBoxMax());
00143 
00144 }
00145 
00146 
00147 void G4FPlane::CalcNormal()
00148 {
00149 /* 
00150   // Calc Normal for surface which is used for the projection
00151   // Make planes
00152   G4Vector3D norm;
00153 
00154   G4Vector3D RefDirection = pplace.GetRefDirection();
00155   G4Vector3D Axis = pplace.GetAxis();
00156 
00157   // L. Broglia : before in G4Placement
00158   if( RefDirection == Axis )
00159     norm = RefDirection;
00160   else
00161   {
00162     // L. Broglia : error on setY, and it`s better to use cross function
00163     // norm.setX( RefDirection.y() * Axis.z() -  RefDirection.z() * Axis.y() );
00164     // norm.setY( RefDirection.x() * Axis.z() -  RefDirection.z() * Axis.x() );
00165     // norm.setZ( RefDirection.x() * Axis.y() -  RefDirection.y() * Axis.x() );
00166        
00167     norm = RefDirection.cross(Axis);
00168   }
00169   
00170   //  const G4Point3D& tmp = pplace.GetSrfPoint();
00171   const G4Point3D tmp = pplace.GetLocation();
00172 */  
00173 
00174   // L. Broglia
00175   // The direction of the normal is the axis of his location
00176   // Its sense depend on the orientation of the bounded curve
00177   const G4Point3D tmp = pplace.GetLocation();
00178   G4Vector3D norm;
00179   G4int sense = GetSameSense();
00180   
00181   if (sense)
00182     norm = pplace.GetAxis();
00183   else
00184     norm = - pplace.GetAxis();
00185 
00186   NormalX =  new G4Ray(tmp, norm);
00187   NormalX->RayCheck();
00188   NormalX->CreatePlanes();
00189 }
00190 
00191 
00192 void G4FPlane::Project()
00193 {
00194     // Project
00195     // const G4Plane& Plane1 = NormalX->GetPlane(1);
00196     // const G4Plane& Plane2 = NormalX->GetPlane(2);
00197 
00198     // probably not necessary
00199     // projections of the boundary should be handled by the intersection
00200     //    OuterBoundary->ProjectBoundaryTo2D(Plane1, Plane2, 0);
00201 }
00202 
00203 
00204 G4int G4FPlane::IsConvex() const
00205 {
00206   return -1;  
00207 }
00208 
00209 
00210 G4int G4FPlane::Intersect(const G4Ray& rayref)
00211 {
00212   // This function count the number of intersections of a 
00213   // bounded surface by a ray.
00214   
00215 
00216   // Find the intersection with the infinite plane
00217   Intersected =1;
00218 
00219   // s is solution, line is p + tq, n is G4Plane Normal, r is point on G4Plane 
00220   // all parameters are pointers to arrays of three elements
00221  
00222   hitpoint = PINFINITY;
00223   register G4double a, b, t;
00224 
00225   register const G4Vector3D& RayDir   = rayref.GetDir();
00226   register const G4Point3D&  RayStart = rayref.GetStart();
00227 
00228   G4double dirx =  RayDir.x();
00229   G4double diry =  RayDir.y();
00230   G4double dirz =  RayDir.z();
00231 
00232   G4Vector3D norm = (*NormalX).GetDir();
00233   G4Point3D  srf_point = pplace.GetLocation();
00234 
00235   b = norm.x() * dirx + norm.y() * diry + norm.z() * dirz;
00236 
00237   if ( std::fabs(b) < perMillion )    
00238   {
00239     // G4cout << "\nLine is parallel to G4Plane.No Hit.";
00240   }  
00241   else
00242   {
00243     G4double startx =  RayStart.x();
00244     G4double starty =  RayStart.y();
00245     G4double startz =  RayStart.z();    
00246     
00247     a = norm.x() * (srf_point.x() - startx) + 
00248         norm.y() * (srf_point.y() - starty) + 
00249         norm.z() * (srf_point.z() - startz)   ;
00250     
00251     t = a/b;
00252     
00253     // substitute t into line equation
00254     // to calculate final solution     
00255     G4double solx,soly,solz;
00256     solx = startx + t * dirx;
00257     soly = starty + t * diry;
00258     solz = startz + t * dirz;
00259 
00260     // solve tolerance problem
00261     if( (t*dirx >= -kCarTolerance/2) && (t*dirx <= kCarTolerance/2) )
00262       solx = startx;
00263 
00264     if( (t*diry >= -kCarTolerance/2) && (t*diry <= kCarTolerance/2) )
00265       soly = starty;
00266 
00267     if( (t*dirz >= -kCarTolerance/2) && (t*dirz <= kCarTolerance/2) )
00268       solz = startz;
00269     
00270     G4bool xhit = (dirx < 0 && solx <= startx) || (dirx >= 0 && solx >= startx);
00271     G4bool yhit = (diry < 0 && soly <= starty) || (diry >= 0 && soly >= starty);
00272     G4bool zhit = (dirz < 0 && solz <= startz) || (dirz >= 0 && solz >= startz);
00273     
00274     if( xhit && yhit && zhit ) {
00275       hitpoint= G4Point3D(solx, soly, solz);
00276     }
00277   }
00278    
00279   // closest_hit is a public Point3D in G4Surface
00280   closest_hit = hitpoint;
00281   
00282   if(closest_hit.x() == kInfinity)
00283   {
00284     // no hit
00285     active=0;
00286     SetDistance(kInfinity);
00287     return 0;
00288   }
00289   else
00290   {
00291     // calculate the squared distance from the point to the intersection
00292     // and set it in the distance data member (all clients know they have
00293     // to take the sqrt)
00294     SetDistance( RayStart.distance2(closest_hit) );   
00295 
00296     // now, we have to verify that the hit point founded
00297     // is included into the G4FPlane boundaries
00298 
00299     // project the hit to the xy plane,
00300     // with the same projection that took the boundary
00301     // into projectedBoundary
00302     G4Point3D projectedHit= pplace.GetToPlacementCoordinates() * closest_hit;
00303     
00304     // test ray from the hit on the xy plane
00305     G4Ray testRay( projectedHit, G4Vector3D(1, 0.01, 0) );
00306 
00307     // check if it intersects the boundary
00308     G4int nbinter = projectedBoundary->IntersectRay2D(testRay);
00309 
00310     // If this number is par, it`s signify that the projected point  
00311     // is outside the projected surface, so the hit point is outside
00312     // the bounded surface
00313     if(nbinter&1)
00314     {
00315       // the intersection point is into the boundaries
00316       // check if the intersection point is on the surface
00317       if(distance <= kCarTolerance*0.5*kCarTolerance*0.5)
00318       {
00319         // the point is on the surface, set the distance to 0            
00320         SetDistance(0);         
00321       }
00322       else
00323       {
00324         // the point is outside the surface
00325       }
00326       
00327       return 1 ;      
00328     }
00329     else
00330     {
00331       // the intersection point is out the boundaries
00332       // it is not a real intersection
00333       active=0;
00334       SetDistance(kInfinity);
00335       return 0;
00336     }
00337   }
00338 }
00339 
00340 
00341 G4double G4FPlane::ClosestDistanceToPoint(const G4Point3D& Pt)
00342 {
00343   // Calculates signed distance of point Pt to G4Plane Pl
00344   // Be careful, the equation of the plane is :
00345   // ax + by + cz = d
00346   G4double dist = Pt.x()*Pl.a + Pt.y()*Pl.b + Pt.z()*Pl.c - Pl.d;
00347 
00348   return dist;
00349 }
00350 
00351 
00352 void G4FPlane::InitBounded()
00353 {
00354   // L. Broglia
00355 
00356   projectedBoundary =
00357     surfaceBoundary.Project( pplace.GetToPlacementCoordinates() );
00358 }
00359 
00360 G4double G4FPlane::HowNear( const G4Vector3D& Pt ) const
00361 {
00362   G4double hownear = Pt.x()*Pl.a + Pt.y()*Pl.b + Pt.z()*Pl.c - Pl.d;
00363 
00364   return hownear;
00365 }

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