G4SimpleLocator.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 // $Id$
00027 //
00028 // Class G4SimpleLocator implementation
00029 //
00030 // 27.10.08 - Tatiana Nikitina, extracted from G4PropagatorInField class
00031 // 04.10.11 - John Apostolakis, revised convergence to use Surface Normal
00032 // ---------------------------------------------------------------------------
00033 
00034 #include <iomanip>
00035 
00036 #include "G4ios.hh"
00037 #include "G4SimpleLocator.hh"
00038 
00039 G4SimpleLocator::G4SimpleLocator(G4Navigator *theNavigator)
00040   : G4VIntersectionLocator(theNavigator)
00041 {
00042 }      
00043 
00044 G4SimpleLocator::~G4SimpleLocator()
00045 {
00046 }
00047 
00048 // --------------------------------------------------------------------------
00049 // G4bool G4PropagatorInField::LocateIntersectionPoint( 
00050 //        const G4FieldTrack&       CurveStartPointVelocity,   // A
00051 //        const G4FieldTrack&       CurveEndPointVelocity,     // B
00052 //        const G4ThreeVector&      TrialPoint,                // E
00053 //              G4FieldTrack&       IntersectedOrRecalculated  // Output
00054 //              G4bool&             recalculated )             // Out
00055 // --------------------------------------------------------------------------
00056 //
00057 // Function that returns the intersection of the true path with the surface
00058 // of the current volume (either the external one or the inner one with one
00059 // of the daughters:
00060 //
00061 //     A = Initial point
00062 //     B = another point 
00063 //
00064 // Both A and B are assumed to be on the true path:
00065 //
00066 //     E is the first point of intersection of the chord AB with 
00067 //     a volume other than A  (on the surface of A or of a daughter)
00068 //
00069 // Convention of Use :
00070 //     i) If it returns "true", then IntersectionPointVelocity is set
00071 //       to the approximate intersection point.
00072 //    ii) If it returns "false", no intersection was found.
00073 //          The validity of IntersectedOrRecalculated depends on 'recalculated'
00074 //        a) if latter is false, then IntersectedOrRecalculated is invalid. 
00075 //        b) if latter is true,  then IntersectedOrRecalculated is
00076 //             the new endpoint, due to a re-integration.
00077 // --------------------------------------------------------------------------
00078 // NOTE: implementation taken from G4PropagatorInField
00079 //
00080 G4bool G4SimpleLocator::EstimateIntersectionPoint( 
00081          const  G4FieldTrack&       CurveStartPointVelocity,     // A
00082          const  G4FieldTrack&       CurveEndPointVelocity,       // B
00083          const  G4ThreeVector&      TrialPoint,                  // E
00084                 G4FieldTrack&       IntersectedOrRecalculatedFT, // Output
00085                 G4bool&             recalculatedEndPoint,        // Out
00086                 G4double            &fPreviousSafety,            //In/Out
00087                 G4ThreeVector       &fPreviousSftOrigin)         //In/Out
00088 {
00089   // Find Intersection Point ( A, B, E )  of true path AB - start at E.
00090 
00091   G4bool found_approximate_intersection = false;
00092   G4bool there_is_no_intersection       = false;
00093   
00094   G4FieldTrack  CurrentA_PointVelocity = CurveStartPointVelocity; 
00095   G4FieldTrack  CurrentB_PointVelocity = CurveEndPointVelocity;
00096   G4ThreeVector CurrentE_Point = TrialPoint;
00097   G4bool        validNormalAtE = false;
00098   G4ThreeVector NormalAtEntry;
00099 
00100   G4FieldTrack  ApproxIntersecPointV(CurveEndPointVelocity); // FT-Def-Construct
00101   G4double      NewSafety = 0.0;
00102   G4bool last_AF_intersection = false;
00103   G4bool final_section = true;  // Shows whether current section is last
00104                                 // (i.e. B=full end)
00105   recalculatedEndPoint = false; 
00106 
00107   G4bool restoredFullEndpoint = false;
00108 
00109   G4int substep_no = 0;
00110 
00111   // Limits for substep number
00112   //
00113   const G4int max_substeps  = 100000000;  // Test 120  (old value 100 )
00114   const G4int warn_substeps = 1000;       //      100  
00115 
00116   // Statistics for substeps
00117   //
00118   static G4int max_no_seen= -1; 
00119 
00120   NormalAtEntry = GetSurfaceNormal( CurrentE_Point, validNormalAtE); 
00121 
00122 #ifdef G4DEBUG_FIELD
00123   static G4double tolerance = 1.0e-8; 
00124   G4ThreeVector  StartPosition= CurveStartPointVelocity.GetPosition(); 
00125   if( (TrialPoint - StartPosition).mag() < tolerance * mm ) 
00126   {
00127      G4Exception("G4SimpleLocator::EstimateIntersectionPoint()", 
00128                  "GeomNav1002", JustWarning,
00129                  "Intersection point F is exactly at start point A." ); 
00130   }
00131 #endif
00132 
00133   do 
00134   {
00135      G4ThreeVector Point_A = CurrentA_PointVelocity.GetPosition();  
00136      G4ThreeVector Point_B = CurrentB_PointVelocity.GetPosition();
00137        
00138      // F = a point on true AB path close to point E 
00139      // (the closest if possible)
00140      //
00141      ApproxIntersecPointV = GetChordFinderFor()
00142                             ->ApproxCurvePointV( CurrentA_PointVelocity, 
00143                                                  CurrentB_PointVelocity, 
00144                                                  CurrentE_Point,
00145                                                  GetEpsilonStepFor());
00146        //  The above method is the key & most intuitive part ...
00147       
00148 #ifdef G4DEBUG_FIELD
00149      if( ApproxIntersecPointV.GetCurveLength() > 
00150          CurrentB_PointVelocity.GetCurveLength() * (1.0 + tolerance) )
00151      {
00152        G4Exception("G4SimpleLocator::EstimateIntersectionPoint()", 
00153                    "GeomNav0003", FatalException,
00154                    "Intermediate F point is past end B point!" ); 
00155      }
00156 #endif
00157 
00158      G4ThreeVector CurrentF_Point= ApproxIntersecPointV.GetPosition();
00159 
00160      // First check whether EF is small - then F is a good approx. point 
00161      // Calculate the length and direction of the chord AF
00162      //
00163      G4ThreeVector  ChordEF_Vector = CurrentF_Point - CurrentE_Point;
00164 
00165      G4ThreeVector  NewMomentumDir= ApproxIntersecPointV.GetMomentumDir(); 
00166      G4double       MomDir_dot_Norm= NewMomentumDir.dot( NormalAtEntry ) ;
00167 
00168      G4ThreeVector  ChordAB           = Point_B - Point_A;
00169 
00170 #ifdef G4DEBUG_FIELD
00171      G4VIntersectionLocator::
00172        ReportTrialStep( substep_no, ChordAB, ChordEF_Vector, 
00173                       NewMomentumDir, NormalAtEntry, validNormalAtE ); 
00174 #endif
00175      // Check Sign is always exiting !! TODO
00176      // Could ( > -epsilon) be used instead?
00177      //
00178      G4bool adequate_angle = ( MomDir_dot_Norm >= 0.0 ) 
00179                           || (! validNormalAtE );  // Invalid
00180      G4double EF_dist2= ChordEF_Vector.mag2();
00181      if ( ( EF_dist2  <= sqr(fiDeltaIntersection) && ( adequate_angle ) )
00182        || ( EF_dist2 <= kCarTolerance*kCarTolerance ) )
00183      {
00184        found_approximate_intersection = true;
00185         
00186        // Create the "point" return value
00187        //
00188        IntersectedOrRecalculatedFT = ApproxIntersecPointV;
00189        IntersectedOrRecalculatedFT.SetPosition( CurrentE_Point );
00190 
00191        if ( GetAdjustementOfFoundIntersection() )
00192        {
00193          // Try to Get Correction of IntersectionPoint using SurfaceNormal()
00194          //  
00195          G4ThreeVector IP;
00196          G4ThreeVector MomentumDir= ApproxIntersecPointV.GetMomentumDirection();
00197          G4bool goodCorrection = AdjustmentOfFoundIntersection( Point_A,
00198                                    CurrentE_Point, CurrentF_Point, MomentumDir,
00199                                    last_AF_intersection, IP, NewSafety,
00200                                    fPreviousSafety, fPreviousSftOrigin );
00201 
00202          if(goodCorrection)
00203          {
00204            IntersectedOrRecalculatedFT = ApproxIntersecPointV;
00205            IntersectedOrRecalculatedFT.SetPosition(IP);
00206          }
00207        }
00208 
00209        // Note: in order to return a point on the boundary, 
00210        //       we must return E. But it is F on the curve.
00211        //       So we must "cheat": we are using the position at point E
00212        //       and the velocity at point F !!!
00213        //
00214        // This must limit the length we can allow for displacement!
00215      }
00216      else  // E is NOT close enough to the curve (ie point F)
00217      {
00218        // Check whether any volumes are encountered by the chord AF
00219        // ---------------------------------------------------------
00220        // First relocate to restore any Voxel etc information
00221        // in the Navigator before calling ComputeStep()
00222        //
00223        GetNavigatorFor()->LocateGlobalPointWithinVolume( Point_A );
00224 
00225        G4ThreeVector PointG;   // Candidate intersection point
00226        G4double stepLengthAF; 
00227        G4bool usedNavigatorAF = false; 
00228        G4bool Intersects_AF = IntersectChord( Point_A,   
00229                                               CurrentF_Point,
00230                                               NewSafety,
00231                                               fPreviousSafety,
00232                                               fPreviousSftOrigin,
00233                                               stepLengthAF,
00234                                               PointG,
00235                                               &usedNavigatorAF );
00236        last_AF_intersection = Intersects_AF;
00237        if( Intersects_AF )
00238        {
00239          // G is our new Candidate for the intersection point.
00240          // It replaces  "E" and we will repeat the test to see if
00241          // it is a good enough approximate point for us.
00242          //       B    <- F
00243          //       E    <- G
00244 
00245          CurrentB_PointVelocity = ApproxIntersecPointV;
00246          CurrentE_Point = PointG;
00247 
00248          // Need to recalculate the Exit Normal at the new PointG 
00249          // Relies on a call to Navigator::ComputeStep in IntersectChord above
00250          // If the safety was adequate (for the step) this would NOT be called!
00251          // But this must not occur, no intersection can be found in that case,
00252          // so this branch, ie if( Intersects_AF ) would not be reached!
00253          //
00254          G4bool validNormalLast; 
00255          NormalAtEntry  = GetSurfaceNormal( PointG, validNormalLast ); 
00256          validNormalAtE = validNormalLast; 
00257 
00258          // By moving point B, must take care if current
00259          // AF has no intersection to try current FB!!
00260          //
00261          final_section= false; 
00262           
00263 #ifdef G4VERBOSE
00264          if( fVerboseLevel > 3 )
00265          {
00266            G4cout << "G4PiF::LI> Investigating intermediate point"
00267                   << " at s=" << ApproxIntersecPointV.GetCurveLength()
00268                   << " on way to full s="
00269                   << CurveEndPointVelocity.GetCurveLength() << G4endl;
00270          }
00271 #endif
00272        }
00273        else  // not Intersects_AF
00274        {  
00275          // In this case:
00276          // There is NO intersection of AF with a volume boundary.
00277          // We must continue the search in the segment FB!
00278          //
00279          GetNavigatorFor()->LocateGlobalPointWithinVolume( CurrentF_Point );
00280 
00281          G4double stepLengthFB;
00282          G4ThreeVector PointH;
00283          G4bool usedNavigatorFB=false; 
00284 
00285          // Check whether any volumes are encountered by the chord FB
00286          // ---------------------------------------------------------
00287 
00288          G4bool Intersects_FB = IntersectChord( CurrentF_Point, Point_B, 
00289                                                 NewSafety,fPreviousSafety,
00290                                                 fPreviousSftOrigin,
00291                                                 stepLengthFB,
00292                                                 PointH, &usedNavigatorFB );
00293          if( Intersects_FB )
00294          { 
00295            // There is an intersection of FB with a volume boundary
00296            // H <- First Intersection of Chord FB 
00297 
00298            // H is our new Candidate for the intersection point.
00299            // It replaces  "E" and we will repeat the test to see if
00300            // it is a good enough approximate point for us.
00301 
00302            // Note that F must be in volume volA  (the same as A)
00303            // (otherwise AF would meet a volume boundary!)
00304            //   A    <- F 
00305            //   E    <- H
00306            //
00307            CurrentA_PointVelocity = ApproxIntersecPointV;
00308            CurrentE_Point = PointH;
00309 
00310            // Need to recalculate the Exit Normal at the new PointG
00311            // Relies on call to Navigator::ComputeStep in IntersectChord above
00312            // If safety was adequate (for the step) this would NOT be called!
00313            // But this must not occur, no intersection found in that case,
00314            // so this branch, i.e. if( Intersects_AF ) would not be reached!
00315            //
00316            G4bool validNormalLast; 
00317            NormalAtEntry  = GetSurfaceNormal( PointH, validNormalLast ); 
00318            validNormalAtE = validNormalLast;
00319          }
00320          else  // not Intersects_FB
00321          {
00322            // There is NO intersection of FB with a volume boundary
00323 
00324            if( final_section  )
00325            {
00326              // If B is the original endpoint, this means that whatever
00327              // volume(s) intersected the original chord, none touch the
00328              // smaller chords we have used.
00329              // The value of 'IntersectedOrRecalculatedFT' returned is
00330              // likely not valid 
00331               
00332              there_is_no_intersection = true;   // real final_section
00333            }
00334            else
00335            {
00336              // We must restore the original endpoint
00337 
00338              CurrentA_PointVelocity = CurrentB_PointVelocity;  // Got to B
00339              CurrentB_PointVelocity = CurveEndPointVelocity;
00340              restoredFullEndpoint = true;
00341            }
00342          } // Endif (Intersects_FB)
00343        } // Endif (Intersects_AF)
00344 
00345        // Ensure that the new endpoints are not further apart in space
00346        // than on the curve due to different errors in the integration
00347        //
00348        G4double linDistSq, curveDist; 
00349        linDistSq = ( CurrentB_PointVelocity.GetPosition() 
00350                    - CurrentA_PointVelocity.GetPosition() ).mag2(); 
00351        curveDist = CurrentB_PointVelocity.GetCurveLength()
00352                    - CurrentA_PointVelocity.GetCurveLength();
00353 
00354        // Change this condition for very strict parameters of propagation 
00355        //
00356        if( curveDist*curveDist*(1+2* GetEpsilonStepFor()) < linDistSq )
00357        {
00358          // Re-integrate to obtain a new B
00359          //
00360          G4FieldTrack newEndPointFT =
00361                  ReEstimateEndpoint( CurrentA_PointVelocity,
00362                                      CurrentB_PointVelocity,
00363                                      linDistSq,    // to avoid recalculation
00364                                      curveDist );
00365          G4FieldTrack oldPointVelB = CurrentB_PointVelocity; 
00366          CurrentB_PointVelocity = newEndPointFT;
00367 
00368          if( (final_section)) // real final section
00369          {
00370            recalculatedEndPoint = true;
00371            IntersectedOrRecalculatedFT = newEndPointFT;
00372              // So that we can return it, if it is the endpoint!
00373          }
00374        }
00375        if( curveDist < 0.0 )
00376        {
00377          fVerboseLevel = 5; // Print out a maximum of information
00378          printStatus( CurrentA_PointVelocity,  CurrentB_PointVelocity,
00379                       -1.0, NewSafety,  substep_no );
00380          std::ostringstream message;
00381          message << "Error in advancing propagation." << G4endl
00382                  << "        Point A (start) is " << CurrentA_PointVelocity
00383                  << G4endl
00384                  << "        Point B (end)   is " << CurrentB_PointVelocity
00385                  << G4endl
00386                  << "        Curve distance is " << curveDist << G4endl
00387                  << G4endl
00388                  << "The final curve point is not further along"
00389                  << " than the original!" << G4endl;
00390 
00391          if( recalculatedEndPoint )
00392          {
00393            message << "Recalculation of EndPoint was called with fEpsStep= "
00394                    << GetEpsilonStepFor() << G4endl;
00395          }
00396          message.precision(20);
00397          message << " Point A (Curve start)   is " << CurveStartPointVelocity
00398                  << G4endl
00399                  << " Point B (Curve   end)   is " << CurveEndPointVelocity
00400                  << G4endl
00401                  << " Point A (Current start) is " << CurrentA_PointVelocity
00402                  << G4endl
00403                  << " Point B (Current end)   is " << CurrentB_PointVelocity
00404                  << G4endl
00405                  << " Point E (Trial Point)   is " << CurrentE_Point
00406                  << G4endl
00407                  << " Point F (Intersection)  is " << ApproxIntersecPointV
00408                  << G4endl
00409                  << "        LocateIntersection parameters are : Substep no= "
00410                  << substep_no;
00411 
00412          G4Exception("G4SimpleLocator::EstimateIntersectionPoint()",
00413                      "GeomNav0003", FatalException, message);
00414        }
00415 
00416        if(restoredFullEndpoint)
00417        {
00418          final_section = restoredFullEndpoint;
00419          restoredFullEndpoint = false;
00420        }
00421      } // EndIf ( E is close enough to the curve, ie point F. )
00422        // tests ChordAF_Vector.mag() <= maximum_lateral_displacement 
00423 
00424 #ifdef G4DEBUG_LOCATE_INTERSECTION  
00425      static G4int trigger_substepno_print= warn_substeps - 20;
00426 
00427      if( substep_no >= trigger_substepno_print )
00428      {
00429        G4cout << "Difficulty in converging in "
00430               << "G4SimpleLocator::EstimateIntersectionPoint():"
00431               << G4endl
00432               << "    Substep no = " << substep_no << G4endl;
00433        if( substep_no == trigger_substepno_print )
00434        {
00435          printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
00436                       -1.0, NewSafety, 0);
00437        }
00438        G4cout << " State of point A: "; 
00439        printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
00440                     -1.0, NewSafety, substep_no-1, 0);
00441        G4cout << " State of point B: "; 
00442        printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
00443                     -1.0, NewSafety, substep_no);
00444      }
00445 #endif
00446      substep_no++; 
00447 
00448   } while ( ( ! found_approximate_intersection )
00449            && ( ! there_is_no_intersection )     
00450            && ( substep_no <= max_substeps) ); // UNTIL found or failed
00451 
00452   if( substep_no > max_no_seen )
00453   {
00454     max_no_seen = substep_no; 
00455 #ifdef G4DEBUG_LOCATE_INTERSECTION  
00456     if( max_no_seen > warn_substeps )
00457     {
00458       trigger_substepno_print = max_no_seen-20; // Want to see last 20 steps 
00459     } 
00460 #endif
00461   }
00462 
00463   if(  ( substep_no >= max_substeps)
00464       && !there_is_no_intersection
00465       && !found_approximate_intersection )
00466   {
00467     G4cout << "ERROR - G4SimpleLocator::EstimateIntersectionPoint()" << G4endl
00468            << "        Start and Endpoint of Requested Step:" << G4endl;
00469     printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
00470                  -1.0, NewSafety, 0);
00471     G4cout << G4endl
00472            << "        Start and end-point of current Sub-Step:" << G4endl;
00473     printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
00474                  -1.0, NewSafety, substep_no-1);
00475     printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
00476                  -1.0, NewSafety, substep_no);
00477 
00478     std::ostringstream message;
00479     message << "Convergence is requiring too many substeps: "
00480             << substep_no << G4endl
00481             << "          Abandoning effort to intersect." << G4endl
00482             << "          Found intersection = "
00483             << found_approximate_intersection << G4endl
00484             << "          Intersection exists = "
00485             << !there_is_no_intersection << G4endl;
00486     message.precision(10); 
00487     G4double done_len = CurrentA_PointVelocity.GetCurveLength(); 
00488     G4double full_len = CurveEndPointVelocity.GetCurveLength();
00489     message << "          Undertaken only length: " << done_len
00490             << " out of " << full_len << " required." << G4endl
00491             << "          Remaining length = " << full_len-done_len; 
00492 
00493     G4Exception("G4SimpleLocator::EstimateIntersectionPoint()",
00494                 "GeomNav0003", FatalException, message);
00495   }
00496   else if( substep_no >= warn_substeps )
00497   {  
00498     std::ostringstream message;
00499     message.precision(10); 
00500 
00501     message << "Many substeps while trying to locate intersection." << G4endl
00502             << "          Undertaken length: "  
00503             << CurrentB_PointVelocity.GetCurveLength() 
00504             << " - Needed: "  << substep_no << " substeps." << G4endl
00505             << "          Warning level = " << warn_substeps
00506             << " and maximum substeps = " << max_substeps;
00507     G4Exception("G4SimpleLocator::EstimateIntersectionPoint()",
00508                 "GeomNav1002", JustWarning, message);
00509   }
00510   return  !there_is_no_intersection; //  Success or failure
00511 }

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