G4GeomTestSegment.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 // G4GeomTestSegment
00033 //
00034 // Author: D.C.Williams, UCSC (davidw@scipp.ucsc.edu)
00035 // --------------------------------------------------------------------
00036 
00037 #include "G4GeomTestSegment.hh"
00038 
00039 #include "G4VSolid.hh"
00040 #include "G4GeomTestLogger.hh"
00041 #include "G4GeometryTolerance.hh"
00042 
00043 //
00044 // Constructor
00045 //
00046 G4GeomTestSegment::G4GeomTestSegment( const G4VSolid *theSolid,
00047                                       const G4ThreeVector &theP,
00048                                       const G4ThreeVector &theV,
00049                                             G4GeomTestLogger *logger )
00050   : solid(theSolid),
00051     p0(theP),
00052     v(theV)
00053 {
00054   kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
00055   FindPoints(logger);
00056 }
00057 
00058 
00059 //
00060 // Simple accessors
00061 //
00062 const G4VSolid * G4GeomTestSegment::GetSolid() const { return solid; }
00063 
00064 const G4ThreeVector &G4GeomTestSegment::GetP() const { return p0; }
00065 
00066 const G4ThreeVector &G4GeomTestSegment::GetV() const { return v; }
00067 
00068 
00069 //
00070 // Return number points
00071 //
00072 G4int G4GeomTestSegment::GetNumberPoints() const
00073 {
00074   return points.size();
00075 }
00076 
00077 
00078 //
00079 // Return ith point
00080 //
00081 const G4GeomTestPoint &G4GeomTestSegment::GetPoint( G4int i ) const
00082 {
00083   return points[i];
00084 }
00085 
00086 
00087 //
00088 // FindPoints
00089 //
00090 // Do the dirty work
00091 //
00092 void G4GeomTestSegment::FindPoints( G4GeomTestLogger *logger )
00093 {
00094   FindSomePoints( logger, true );
00095   FindSomePoints( logger, false );
00096   
00097   PatchInconsistencies( logger );
00098 }
00099 
00100 
00101 //
00102 // PatchInconsistencies
00103 //
00104 // Search for inconsistancies in the hit list and patch
00105 // them up, if possible
00106 //
00107 void G4GeomTestSegment::PatchInconsistencies(  G4GeomTestLogger *logger )
00108 {
00109   if (points.size() == 0) return;
00110   
00111   //
00112   // Sort
00113   //
00114   std::sort( points.begin(), points.end() );
00115   
00116   //
00117   // Loop over entering/leaving pairs
00118   //
00119   std::vector<G4GeomTestPoint>::iterator curr = points.begin();
00120   do {
00121     
00122     std::vector<G4GeomTestPoint>::iterator next = curr + 1;
00123   
00124     //
00125     // Is the next point close by?
00126     //
00127     while (next != points.end() && 
00128            next->GetDistance()-curr->GetDistance() < kCarTolerance) {
00129       //
00130       // Yup. If we have an entering/leaving pair, all is okay
00131       //
00132       if (next->Entering() != curr->Entering()) {
00133         curr = next + 1;
00134         next = curr + 1;
00135         
00136         break;
00137       }
00138       
00139       //
00140       // Otherwise, they are duplicate, and just delete one
00141       //
00142       next = points.erase(next);
00143       curr = next - 1;
00144      
00145     }
00146     
00147     if (curr == points.end()) break;
00148   
00149     //
00150     // The next point should be entering...
00151     //
00152     
00153     if (!curr->Entering()) {
00154       //
00155       // Point is out of sequence:
00156       // Test solid just before this point
00157       //
00158       G4double ds = curr->GetDistance();
00159       G4ThreeVector p = p0 + ds*v;
00160       
00161       G4ThreeVector p1 = p - 10*kCarTolerance*v;
00162       
00163       if (solid->Inside(p1) == kOutside||(solid->Inside(p1)== kSurface)) {
00164         //
00165         // We are missing an entrance point near the current
00166         // point. Add one.
00167         //
00168         curr = points.insert(curr, G4GeomTestPoint( p, ds, true ) );
00169         ++curr;
00170         break;
00171       }
00172       
00173       //
00174       // Test solid just after last point
00175       //
00176       if (curr != points.begin()) {
00177         std::vector<G4GeomTestPoint>::iterator prev = curr - 1;
00178 
00179         ds = prev->GetDistance();
00180         p = p0 + ds*v;
00181         
00182         p1 = p + 10*kCarTolerance*v;
00183         if ((solid->Inside(p1) == kOutside)||(solid->Inside(p1)== kSurface)) {
00184           
00185           //
00186           // We are missing an entrance point near the previous
00187           // point. Add one.
00188           //
00189           curr = points.insert(curr, G4GeomTestPoint( p, ds, true ) );
00190           ++curr;
00191           break;
00192         }
00193       }
00194       
00195       //
00196       // Oh oh. No solution. Delete the current point and complain.
00197       //
00198       logger->SolidProblem( solid, "Spurious exiting intersection point", p );
00199       curr = points.erase(curr);
00200       break;
00201     }  
00202 
00203     //
00204     // The following point should be leaving
00205     //
00206      
00207     if (next == points.end() || next->Entering() ) {
00208       //
00209       // But is missing. Check immediately after this point.
00210       //
00211       G4double ds = curr->GetDistance();
00212       G4ThreeVector p = p0 + ds*v;
00213       G4ThreeVector p1 = p + 10*kCarTolerance*v;
00214       
00215       if (solid->Inside(p1) == kOutside) {
00216         //
00217         // We are missing an exit point near the current
00218         // point. Add one.
00219         //
00220         curr = points.insert(next, G4GeomTestPoint( p, ds, false ) );
00221         break;
00222       }
00223       
00224       if (next != points.end()) {
00225         //
00226         // Check just before next point
00227         //
00228         ds = next->GetDistance();
00229         p = p0 + ds*v;
00230         p1 = p - 10*kCarTolerance*v;
00231         if (solid->Inside(p1) == kOutside) {
00232           //
00233           // We are missing an exit point before the next
00234           // point. Add one.
00235           //
00236           curr = points.insert(next, G4GeomTestPoint( p, ds, false ) );
00237           break;
00238         }
00239       }        
00240       
00241       //
00242       // Oh oh. Hanging entering point. Delete and complain.
00243       //
00244       logger->SolidProblem( solid, "Spurious entering intersection point", p );
00245       curr = points.erase(curr);
00246     }
00247     
00248     if(curr!=points.end()){curr = next + 1;}
00249     
00250   } while( curr != points.end() );
00251   
00252   //
00253   // Double check
00254   //
00255   if (points.size()&1) 
00256     logger->SolidProblem( solid,
00257                           "Solid has odd number of intersection points", p0 );
00258   
00259   return;
00260 }  
00261 
00262 
00263 //
00264 // Find intersections either in front or behind the point
00265 //
00266 void G4GeomTestSegment::FindSomePoints( G4GeomTestLogger *logger,
00267                                         G4bool forward )
00268 { 
00269   G4double sign = forward ? +1 : -1;
00270 
00271   G4ThreeVector p(p0);
00272   G4ThreeVector vSearch(sign*v);
00273   G4double ds(0);
00274   G4bool entering;
00275   G4double vSurfN;
00276  
00277   // Look for nearest intersection point in the specified
00278   // direction and return if there isn't one
00279   //
00280   G4double dist;
00281   switch(solid->Inside(p)) {
00282     case kInside:
00283       dist = solid->DistanceToOut(p,vSearch);
00284       if (dist >= kInfinity) {
00285         logger->SolidProblem( solid,
00286                 "DistanceToOut(p,v) = kInfinity for point inside", p );
00287         return;
00288       }
00289       ds += sign*dist;
00290       entering = false;
00291       break;
00292     case kOutside:
00293       dist = solid->DistanceToIn(p,vSearch);
00294       if (dist >= kInfinity) return;
00295       ds += sign*dist;
00296       entering = true;
00297       break;
00298     case kSurface:
00299       vSurfN=v.dot(solid->SurfaceNormal(p));
00300       if(std::fabs(vSurfN)<kCarTolerance)vSurfN=0;
00301       entering = (vSurfN < 0);
00302       break;
00303     default:
00304       logger->SolidProblem( solid,
00305               "Inside returns illegal enumerated value", p );
00306       return;
00307   }
00308   
00309   //
00310   // Loop
00311   //
00312   // nzero = the number of consecutive zeros
00313   //
00314   G4int nzero=0;
00315   
00316   for(;;) {
00317     //
00318     // Locate point
00319     //
00320     p = p0 + ds*v;
00321     
00322     if (nzero > 2) {
00323       //
00324       // Oops. In a loop. Probably along a spherical or cylindrical surface.
00325       // Let's give the tool a little help with a push
00326       //
00327       G4double push = 1E-6;
00328       ds += sign*push;
00329       for(;;) {
00330         p = p0 + ds*v;
00331         EInside inside = solid->Inside(p);
00332         if (inside == kInside) {
00333           entering = true;
00334           break;
00335         }
00336         else if (inside == kOutside) {
00337           entering = false;
00338           break;
00339         }
00340 
00341         push = 2*push;
00342         if (push > 0.1) {
00343           logger->SolidProblem( solid,
00344                   "Push fails to fix geometry inconsistency", p );
00345           return;
00346         }
00347         ds += sign*push;
00348       }
00349     }
00350     else {
00351 
00352       //
00353       // Record point
00354       //
00355       points.push_back( G4GeomTestPoint( p, ds, entering==forward ) );
00356       
00357     }
00358     
00359     //
00360     // Find next intersection
00361     //
00362     if (entering) {
00363       dist = solid->DistanceToOut(p,vSearch);
00364       if (dist >= kInfinity) {
00365         logger->SolidProblem( solid,
00366                 "DistanceToOut(p,v) = kInfinity for point inside", p );
00367         return;
00368       }
00369       
00370       if ( (dist > kCarTolerance)
00371         && (solid->Inside(p + (dist*0.999)*vSearch) == kOutside) ) {
00372         //
00373         // This shouldn't have happened
00374         //
00375         if (solid->Inside(p) == kOutside)
00376           logger->SolidProblem(solid,
00377                   "Entering point is outside (possible roundoff error)",p);
00378         else
00379           logger->SolidProblem(solid,
00380                   "DistanceToOut(p,v) brings trajectory well outside solid",p);
00381         return;
00382       }
00383             
00384       entering = false;
00385     }
00386     else {
00387       dist = solid->DistanceToIn(p,vSearch);
00388       if (dist >= kInfinity) return;
00389       
00390       if ( (dist > kCarTolerance)
00391         && (solid->Inside(p + (dist*0.999)*vSearch) == kInside) ) {
00392         //
00393         // This shouldn't have happened
00394         //
00395         if (solid->Inside(p) == kInside)
00396           logger->SolidProblem(solid,
00397                   "Exiting point is inside (possible roundoff error)", p);
00398         else
00399           logger->SolidProblem(solid,
00400                   "DistanceToIn(p,v) brings trajectory well inside solid", p);
00401         return;
00402       }
00403       
00404       entering = true;
00405     }
00406     
00407     //
00408     // Update ds
00409     //
00410     if (dist <= 0) {
00411       nzero++; 
00412     }
00413     else {
00414       nzero=0;
00415       ds += sign*dist;
00416     }
00417   }
00418 }

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