Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/geometry/navigation/src/G4BrentLocator.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /geometry/navigation/src/G4BrentLocator.cc (Version 11.3.0) and /geometry/navigation/src/G4BrentLocator.cc (Version 11.2.2)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
 26 // class G4BrentLocator implementation             26 // class G4BrentLocator implementation
 27 //                                                 27 //
 28 // 27.10.08 - Tatiana Nikitina.                    28 // 27.10.08 - Tatiana Nikitina.
 29 // 04.10.11 - John Apostolakis, revised conver     29 // 04.10.11 - John Apostolakis, revised convergence to use Surface Normal
 30 // -------------------------------------------     30 // ---------------------------------------------------------------------------
 31                                                    31 
 32 #include <iomanip>                                 32 #include <iomanip>
 33                                                    33 
 34 #include "G4BrentLocator.hh"                       34 #include "G4BrentLocator.hh"
 35 #include "G4ios.hh"                                35 #include "G4ios.hh"
 36                                                    36 
 37 G4BrentLocator::G4BrentLocator(G4Navigator *th     37 G4BrentLocator::G4BrentLocator(G4Navigator *theNavigator)
 38   : G4VIntersectionLocator(theNavigator)           38   : G4VIntersectionLocator(theNavigator)
 39 {                                                  39 {
 40   // In case of too slow progress in finding I     40   // In case of too slow progress in finding Intersection Point
 41   // intermediates Points on the Track must be     41   // intermediates Points on the Track must be stored.
 42   // Initialise the array of Pointers [max_dep     42   // Initialise the array of Pointers [max_depth+1] to do this  
 43                                                    43   
 44   G4ThreeVector zeroV(0.0,0.0,0.0);                44   G4ThreeVector zeroV(0.0,0.0,0.0);
 45   for (auto & idepth : ptrInterMedFT)              45   for (auto & idepth : ptrInterMedFT)
 46   {                                                46   {
 47     idepth = new G4FieldTrack( zeroV, zeroV, 0     47     idepth = new G4FieldTrack( zeroV, zeroV, 0., 0., 0., 0.);
 48   }                                                48   }
 49 }                                                  49 }
 50                                                    50 
 51 G4BrentLocator::~G4BrentLocator()                  51 G4BrentLocator::~G4BrentLocator()
 52 {                                                  52 {
 53   for (auto & idepth : ptrInterMedFT)              53   for (auto & idepth : ptrInterMedFT)
 54   {                                                54   {
 55     delete idepth;                                 55     delete idepth;
 56   }                                                56   }
 57 }                                                  57 }
 58                                                    58 
 59 // -------------------------------------------     59 // --------------------------------------------------------------------------
 60 // G4bool G4PropagatorInField::LocateIntersect     60 // G4bool G4PropagatorInField::LocateIntersectionPoint( 
 61 //        const G4FieldTrack&       CurveStart     61 //        const G4FieldTrack&       CurveStartPointVelocity,   //  A
 62 //        const G4FieldTrack&       CurveEndPo     62 //        const G4FieldTrack&       CurveEndPointVelocity,     //  B
 63 //        const G4ThreeVector&      TrialPoint     63 //        const G4ThreeVector&      TrialPoint,                //  E
 64 //              G4FieldTrack&       Intersecte     64 //              G4FieldTrack&       IntersectedOrRecalculated  // Output
 65 //              G4bool&             recalculat     65 //              G4bool&             recalculated)              // Out
 66 // -------------------------------------------     66 // --------------------------------------------------------------------------
 67 //                                                 67 //
 68 // Function that returns the intersection of t     68 // Function that returns the intersection of the true path with the surface
 69 // of the current volume (either the external      69 // of the current volume (either the external one or the inner one with one
 70 // of the daughters:                               70 // of the daughters:
 71 //                                                 71 //
 72 //     A = Initial point                           72 //     A = Initial point
 73 //     B = another point                           73 //     B = another point 
 74 //                                                 74 //
 75 // Both A and B are assumed to be on the true      75 // Both A and B are assumed to be on the true path:
 76 //                                                 76 //
 77 //     E is the first point of intersection of     77 //     E is the first point of intersection of the chord AB with 
 78 //     a volume other than A  (on the surface      78 //     a volume other than A  (on the surface of A or of a daughter)
 79 //                                                 79 //
 80 // Convention of Use :                             80 // Convention of Use :
 81 //     i) If it returns "true", then Intersect     81 //     i) If it returns "true", then IntersectionPointVelocity is set
 82 //        to the approximate intersection poin     82 //        to the approximate intersection point.
 83 //    ii) If it returns "false", no intersecti     83 //    ii) If it returns "false", no intersection was found.
 84 //        The validity of IntersectedOrRecalcu     84 //        The validity of IntersectedOrRecalculated depends on 'recalculated'
 85 //        a) if latter is false, then Intersec     85 //        a) if latter is false, then IntersectedOrRecalculated is invalid. 
 86 //        b) if latter is true,  then Intersec     86 //        b) if latter is true,  then IntersectedOrRecalculated is
 87 //           the new endpoint, due to a re-int     87 //           the new endpoint, due to a re-integration.
 88 // -------------------------------------------     88 // --------------------------------------------------------------------------
 89 // NOTE: implementation taken from G4Propagato     89 // NOTE: implementation taken from G4PropagatorInField
 90 //       New second order locator is added         90 //       New second order locator is added 
 91 //                                                 91 //
 92 G4bool G4BrentLocator::EstimateIntersectionPoi     92 G4bool G4BrentLocator::EstimateIntersectionPoint( 
 93          const  G4FieldTrack&       CurveStart     93          const  G4FieldTrack&       CurveStartPointVelocity,       // A
 94          const  G4FieldTrack&       CurveEndPo     94          const  G4FieldTrack&       CurveEndPointVelocity,         // B
 95          const  G4ThreeVector&      TrialPoint     95          const  G4ThreeVector&      TrialPoint,                    // E
 96                 G4FieldTrack&       Intersecte     96                 G4FieldTrack&       IntersectedOrRecalculatedFT,   // Output
 97                 G4bool&             recalculat     97                 G4bool&             recalculatedEndPoint,          // Out
 98                 G4double&           fPreviousS     98                 G4double&           fPreviousSafety,               // In/Out
 99                 G4ThreeVector&      fPreviousS     99                 G4ThreeVector&      fPreviousSftOrigin)            // In/Out 
100                                                   100               
101 {                                                 101 {   
102   // Find Intersection Point ( A, B, E )  of t    102   // Find Intersection Point ( A, B, E )  of true path AB - start at E.
103                                                   103 
104   G4bool found_approximate_intersection = fals    104   G4bool found_approximate_intersection = false;
105   G4bool there_is_no_intersection       = fals    105   G4bool there_is_no_intersection       = false;
106                                                   106   
107   G4FieldTrack  CurrentA_PointVelocity = Curve    107   G4FieldTrack  CurrentA_PointVelocity = CurveStartPointVelocity; 
108   G4FieldTrack  CurrentB_PointVelocity = Curve    108   G4FieldTrack  CurrentB_PointVelocity = CurveEndPointVelocity;
109   G4ThreeVector CurrentE_Point = TrialPoint;      109   G4ThreeVector CurrentE_Point = TrialPoint;
110   G4bool        validNormalAtE = false;           110   G4bool        validNormalAtE = false;
111   G4ThreeVector NormalAtEntry;                    111   G4ThreeVector NormalAtEntry;
112                                                   112 
113   G4FieldTrack  ApproxIntersecPointV(CurveEndP    113   G4FieldTrack  ApproxIntersecPointV(CurveEndPointVelocity); // FT-Def-Construct
114   G4double      NewSafety = 0.0;                  114   G4double      NewSafety = 0.0;
115   G4bool        last_AF_intersection = false;     115   G4bool        last_AF_intersection = false; 
116                                                   116 
117   // G4bool final_section= true;  // Shows whe    117   // G4bool final_section= true;  // Shows whether current section is last
118                                   // (i.e. B=f    118                                   // (i.e. B=full end)
119   G4bool first_section = true;                    119   G4bool first_section = true;
120   recalculatedEndPoint = false;                   120   recalculatedEndPoint = false; 
121                                                   121 
122   G4bool restoredFullEndpoint = false;            122   G4bool restoredFullEndpoint = false;
123                                                   123 
124   G4long oldprc;  // cout, cerr precision         124   G4long oldprc;  // cout, cerr precision
125   G4int substep_no = 0;                           125   G4int substep_no = 0;
126                                                   126    
127   // Limits for substep number                    127   // Limits for substep number
128   //                                              128   //
129   const G4int max_substeps=   10000;  // Test     129   const G4int max_substeps=   10000;  // Test 120  (old value 100 )
130   const G4int warn_substeps=   1000;  //          130   const G4int warn_substeps=   1000;  //      100  
131                                                   131 
132   // Statistics for substeps                      132   // Statistics for substeps
133   //                                              133   //
134   static G4ThreadLocal G4int max_no_seen= -1;     134   static G4ThreadLocal G4int max_no_seen= -1; 
135                                                   135 
136   // Counter for restarting Bintermed             136   // Counter for restarting Bintermed
137   //                                              137   //
138   G4int restartB = 0;                             138   G4int restartB = 0;
139                                                   139 
140   //------------------------------------------    140   //--------------------------------------------------------------------------  
141   //  Algorithm for the case if progress in fo    141   //  Algorithm for the case if progress in founding intersection is too slow.
142   //  Process is defined too slow if after N=p    142   //  Process is defined too slow if after N=param_substeps advances on the
143   //  path, it will be only 'fraction_done' of    143   //  path, it will be only 'fraction_done' of the total length.
144   //  In this case the remaining length is div    144   //  In this case the remaining length is divided in two half and 
145   //  the loop is restarted for each half.        145   //  the loop is restarted for each half.  
146   //  If progress is still too slow, the divis    146   //  If progress is still too slow, the division in two halfs continue
147   //  until 'max_depth'.                          147   //  until 'max_depth'.
148   //------------------------------------------    148   //--------------------------------------------------------------------------
149                                                   149 
150   const G4int param_substeps = 50; // Test val    150   const G4int param_substeps = 50; // Test value for the maximum number
151                                    // of subst    151                                    // of substeps
152   const G4double fraction_done = 0.3;             152   const G4double fraction_done = 0.3;
153                                                   153 
154   G4bool Second_half = false;     // First hal    154   G4bool Second_half = false;     // First half or second half of divided step
155                                                   155 
156   NormalAtEntry = GetSurfaceNormal(CurrentE_Po    156   NormalAtEntry = GetSurfaceNormal(CurrentE_Point, validNormalAtE); 
157                                                   157 
158   // We need to know this for the 'final_secti    158   // We need to know this for the 'final_section':
159   // real 'final_section' or first half 'final    159   // real 'final_section' or first half 'final_section'
160   // In algorithm it is considered that the 'S    160   // In algorithm it is considered that the 'Second_half' is true
161   // and it becomes false only if we are in th    161   // and it becomes false only if we are in the first-half of level
162   // depthness or if we are in the first secti    162   // depthness or if we are in the first section
163                                                   163 
164   G4int depth = 0; // Depth counts how many su    164   G4int depth = 0; // Depth counts how many subdivisions of initial step made
165                                                   165 
166 #ifdef G4DEBUG_FIELD                              166 #ifdef G4DEBUG_FIELD
167   const G4double tolerance = 1.0e-8;              167   const G4double tolerance = 1.0e-8; 
168   G4ThreeVector  StartPosition = CurveStartPoi    168   G4ThreeVector  StartPosition = CurveStartPointVelocity.GetPosition(); 
169   if( (TrialPoint - StartPosition).mag() < tol    169   if( (TrialPoint - StartPosition).mag() < tolerance * CLHEP::mm ) 
170   {                                               170   {
171      G4Exception("G4BrentLocator::EstimateInte    171      G4Exception("G4BrentLocator::EstimateIntersectionPoint()", 
172                  "GeomNav1002", JustWarning,      172                  "GeomNav1002", JustWarning,
173                  "Intersection point F is exac    173                  "Intersection point F is exactly at start point A." ); 
174   }                                               174   }
175 #endif                                            175 #endif
176                                                   176 
177   // Intermediates Points on the Track = Subdi    177   // Intermediates Points on the Track = Subdivided Points must be stored.
178   // Give the initial values to 'InterMedFt'      178   // Give the initial values to 'InterMedFt'
179   // Important is 'ptrInterMedFT[0]', it saves    179   // Important is 'ptrInterMedFT[0]', it saves the 'EndCurvePoint'
180   //                                              180   //
181   *ptrInterMedFT[0] = CurveEndPointVelocity;      181   *ptrInterMedFT[0] = CurveEndPointVelocity;
182   for (auto idepth=1; idepth<max_depth+1; ++id    182   for (auto idepth=1; idepth<max_depth+1; ++idepth )
183   {                                               183   {
184     *ptrInterMedFT[idepth] = CurveStartPointVe    184     *ptrInterMedFT[idepth] = CurveStartPointVelocity;
185   }                                               185   }
186                                                   186 
187   //Final_section boolean store                   187   //Final_section boolean store
188   G4bool fin_section_depth[max_depth];            188   G4bool fin_section_depth[max_depth];
189   for (bool & idepth : fin_section_depth)         189   for (bool & idepth : fin_section_depth)
190   {                                               190   {
191     idepth = true;                                191     idepth = true;
192   }                                               192   }
193                                                   193 
194   // 'SubStartPoint' is needed to calculate th    194   // 'SubStartPoint' is needed to calculate the length of the divided step
195   //                                              195   //
196   G4FieldTrack SubStart_PointVelocity = CurveS    196   G4FieldTrack SubStart_PointVelocity = CurveStartPointVelocity;
197                                                   197    
198   do   // Loop checking, 07.10.2016, J.Apostol    198   do   // Loop checking, 07.10.2016, J.Apostolakis
199   {                                               199   {
200     G4int substep_no_p = 0;                       200     G4int substep_no_p = 0;
201     G4bool sub_final_section = false; // the s    201     G4bool sub_final_section = false; // the same as final_section,
202                                       // but f    202                                       // but for 'sub_section'
203     SubStart_PointVelocity = CurrentA_PointVel    203     SubStart_PointVelocity = CurrentA_PointVelocity;
204                                                   204 
205     do   // Loop checking, 07.10.2016, J.Apost    205     do   // Loop checking, 07.10.2016, J.Apostolakis
206     { // REPEAT param                             206     { // REPEAT param 
207       G4ThreeVector Point_A = CurrentA_PointVe    207       G4ThreeVector Point_A = CurrentA_PointVelocity.GetPosition();  
208       G4ThreeVector Point_B = CurrentB_PointVe    208       G4ThreeVector Point_B = CurrentB_PointVelocity.GetPosition();
209                                                   209        
210       // F = a point on true AB path close to     210       // F = a point on true AB path close to point E 
211       // (the closest if possible)                211       // (the closest if possible)
212       //                                          212       //
213       if(substep_no_p==0)                         213       if(substep_no_p==0)
214       {                                           214       {
215         ApproxIntersecPointV = GetChordFinderF    215         ApproxIntersecPointV = GetChordFinderFor()
216                                ->ApproxCurvePo    216                                ->ApproxCurvePointV( CurrentA_PointVelocity, 
217                                                   217                                                     CurrentB_PointVelocity, 
218                                                   218                                                     CurrentE_Point,
219                                                   219                                                     GetEpsilonStepFor());
220           //  The above method is the key & mo    220           //  The above method is the key & most intuitive part ...
221       }                                           221       }
222 #ifdef G4DEBUG_FIELD                              222 #ifdef G4DEBUG_FIELD
223       if( ApproxIntersecPointV.GetCurveLength(    223       if( ApproxIntersecPointV.GetCurveLength() > 
224           CurrentB_PointVelocity.GetCurveLengt    224           CurrentB_PointVelocity.GetCurveLength() * (1.0 + tolerance) )
225       {                                           225       {
226         G4Exception("G4BrentLocator::EstimateI    226         G4Exception("G4BrentLocator::EstimateIntersectionPoint()", 
227                     "GeomNav0003", FatalExcept    227                     "GeomNav0003", FatalException,
228                     "Intermediate F point is p    228                     "Intermediate F point is past end B point" ); 
229       }                                           229       }
230 #endif                                            230 #endif
231                                                   231 
232       G4ThreeVector CurrentF_Point = ApproxInt    232       G4ThreeVector CurrentF_Point = ApproxIntersecPointV.GetPosition();
233                                                   233 
234       // First check whether EF is small - the    234       // First check whether EF is small - then F is a good approx. point 
235       // Calculate the length and direction of    235       // Calculate the length and direction of the chord AF
236       //                                          236       //
237       G4ThreeVector  ChordEF_Vector = CurrentF    237       G4ThreeVector  ChordEF_Vector = CurrentF_Point - CurrentE_Point;
238       G4ThreeVector  NewMomentumDir = ApproxIn    238       G4ThreeVector  NewMomentumDir = ApproxIntersecPointV.GetMomentumDir(); 
239       G4double       MomDir_dot_Norm = NewMome    239       G4double       MomDir_dot_Norm = NewMomentumDir.dot( NormalAtEntry ) ;
240                                                   240      
241 #ifdef G4DEBUG_FIELD                              241 #ifdef G4DEBUG_FIELD
242       G4ThreeVector  ChordAB = Point_B - Point    242       G4ThreeVector  ChordAB = Point_B - Point_A;
243                                                   243 
244       G4VIntersectionLocator::ReportTrialStep(    244       G4VIntersectionLocator::ReportTrialStep( substep_no, ChordAB,
245                ChordEF_Vector, NewMomentumDir,    245                ChordEF_Vector, NewMomentumDir, NormalAtEntry, validNormalAtE ); 
246 #endif                                            246 #endif
247                                                   247 
248       G4bool adequate_angle;                      248       G4bool adequate_angle;
249       adequate_angle =  ( MomDir_dot_Norm >= 0    249       adequate_angle =  ( MomDir_dot_Norm >= 0.0 ) // Can use -epsilon instead.
250                     || (! validNormalAtE );       250                     || (! validNormalAtE );        //  Makes criterion invalid
251       G4double EF_dist2 = ChordEF_Vector.mag2(    251       G4double EF_dist2 = ChordEF_Vector.mag2();
252       if ( ( EF_dist2 <= sqr(fiDeltaIntersecti    252       if ( ( EF_dist2 <= sqr(fiDeltaIntersection) && ( adequate_angle ) )
253         || ( EF_dist2 <= kCarTolerance*kCarTol    253         || ( EF_dist2 <= kCarTolerance*kCarTolerance ) )
254       {                                           254       {
255         found_approximate_intersection = true;    255         found_approximate_intersection = true;
256                                                   256     
257         // Create the "point" return value        257         // Create the "point" return value
258         //                                        258         //
259         IntersectedOrRecalculatedFT = ApproxIn    259         IntersectedOrRecalculatedFT = ApproxIntersecPointV;
260         IntersectedOrRecalculatedFT.SetPositio    260         IntersectedOrRecalculatedFT.SetPosition( CurrentE_Point );
261                                                   261         
262         if ( GetAdjustementOfFoundIntersection    262         if ( GetAdjustementOfFoundIntersection() )
263         {                                         263         {
264           // Try to Get Correction of Intersec    264           // Try to Get Correction of IntersectionPoint using SurfaceNormal()
265           //                                      265           //  
266           G4ThreeVector IP;                       266           G4ThreeVector IP;
267           G4ThreeVector MomentumDir=ApproxInte    267           G4ThreeVector MomentumDir=ApproxIntersecPointV.GetMomentumDirection();
268           G4bool goodCorrection = AdjustmentOf    268           G4bool goodCorrection = AdjustmentOfFoundIntersection( Point_A,
269                                     CurrentE_P    269                                     CurrentE_Point, CurrentF_Point, MomentumDir,
270                                     last_AF_in    270                                     last_AF_intersection, IP, NewSafety,
271                                     fPreviousS    271                                     fPreviousSafety, fPreviousSftOrigin );
272           if ( goodCorrection )                   272           if ( goodCorrection )
273           {                                       273           {
274             IntersectedOrRecalculatedFT = Appr    274             IntersectedOrRecalculatedFT = ApproxIntersecPointV;
275             IntersectedOrRecalculatedFT.SetPos    275             IntersectedOrRecalculatedFT.SetPosition(IP);
276           }                                       276           }
277         }                                         277         }
278                                                   278        
279         // Note: in order to return a point on    279         // Note: in order to return a point on the boundary, 
280         //       we must return E. But it is F    280         //       we must return E. But it is F on the curve.
281         //       So we must "cheat": we are us    281         //       So we must "cheat": we are using the position at point E
282         //       and the velocity at point F !    282         //       and the velocity at point F !!!
283         //                                        283         //
284         // This must limit the length we can a    284         // This must limit the length we can allow for displacement!
285       }                                           285       }
286       else  // E is NOT close enough to the cu    286       else  // E is NOT close enough to the curve (ie point F)
287       {                                           287       {
288         // Check whether any volumes are encou    288         // Check whether any volumes are encountered by the chord AF
289         // -----------------------------------    289         // ---------------------------------------------------------
290         // First relocate to restore any Voxel    290         // First relocate to restore any Voxel etc information
291         // in the Navigator before calling Com    291         // in the Navigator before calling ComputeStep()
292         //                                        292         //
293         GetNavigatorFor()->LocateGlobalPointWi    293         GetNavigatorFor()->LocateGlobalPointWithinVolume( Point_A );
294                                                   294 
295         G4ThreeVector PointG;   // Candidate i    295         G4ThreeVector PointG;   // Candidate intersection point
296         G4double stepLengthAF;                    296         G4double stepLengthAF; 
297         G4bool usedNavigatorAF = false;           297         G4bool usedNavigatorAF = false; 
298         G4bool Intersects_AF = IntersectChord(    298         G4bool Intersects_AF = IntersectChord( Point_A,   CurrentF_Point,
299                                                   299                                                NewSafety,fPreviousSafety,
300                                                   300                                                fPreviousSftOrigin,
301                                                   301                                                stepLengthAF,
302                                                   302                                                PointG,
303                                                   303                                                &usedNavigatorAF);
304         last_AF_intersection = Intersects_AF;     304         last_AF_intersection = Intersects_AF;
305         if( Intersects_AF )                       305         if( Intersects_AF )
306         {                                         306         {
307           // G is our new Candidate for the in    307           // G is our new Candidate for the intersection point.
308           // It replaces  "E" and we will repe    308           // It replaces  "E" and we will repeat the test to see if
309           // it is a good enough approximate p    309           // it is a good enough approximate point for us.
310           //       B    <- F                      310           //       B    <- F
311           //       E    <- G                      311           //       E    <- G
312           //                                      312           //
313           G4FieldTrack EndPoint = ApproxInters    313           G4FieldTrack EndPoint = ApproxIntersecPointV;
314           ApproxIntersecPointV = GetChordFinde    314           ApproxIntersecPointV = GetChordFinderFor()->ApproxCurvePointS(
315                                  CurrentA_Poin    315                                  CurrentA_PointVelocity, CurrentB_PointVelocity,
316                                  EndPoint,Curr    316                                  EndPoint,CurrentE_Point, CurrentF_Point,PointG,
317                                  true, GetEpsi    317                                  true, GetEpsilonStepFor() );
318           CurrentB_PointVelocity = EndPoint;      318           CurrentB_PointVelocity = EndPoint;
319           CurrentE_Point = PointG;                319           CurrentE_Point = PointG;
320                                                   320 
321           // Need to recalculate the Exit Norm    321           // Need to recalculate the Exit Normal at the new PointG 
322           // Know that a call was made to Navi    322           // Know that a call was made to Navigator::ComputeStep in
323           // IntersectChord above.                323           // IntersectChord above.
324           //                                      324           //
325           G4bool validNormalLast;                 325           G4bool validNormalLast; 
326           NormalAtEntry  = GetSurfaceNormal( P    326           NormalAtEntry  = GetSurfaceNormal( PointG, validNormalLast ); 
327           validNormalAtE = validNormalLast;       327           validNormalAtE = validNormalLast; 
328                                                   328             
329           // By moving point B, must take care    329           // By moving point B, must take care if current
330           // AF has no intersection to try cur    330           // AF has no intersection to try current FB!!
331           //                                      331           //
332           fin_section_depth[depth] = false;       332           fin_section_depth[depth] = false;
333 #ifdef G4VERBOSE                                  333 #ifdef G4VERBOSE
334           if( fVerboseLevel > 3 )                 334           if( fVerboseLevel > 3 )
335           {                                       335           {
336             G4cout << "G4PiF::LI> Investigatin    336             G4cout << "G4PiF::LI> Investigating intermediate point"
337                    << " at s=" << ApproxInters    337                    << " at s=" << ApproxIntersecPointV.GetCurveLength()
338                    << " on way to full s="        338                    << " on way to full s="
339                    << CurveEndPointVelocity.Ge    339                    << CurveEndPointVelocity.GetCurveLength() << G4endl;
340           }                                       340           }
341 #endif                                            341 #endif
342         }                                         342         }
343         else  // not Intersects_AF                343         else  // not Intersects_AF
344         {                                         344         {  
345           // In this case:                        345           // In this case:
346           // There is NO intersection of AF wi    346           // There is NO intersection of AF with a volume boundary.
347           // We must continue the search in th    347           // We must continue the search in the segment FB!
348           //                                      348           //
349           GetNavigatorFor()->LocateGlobalPoint    349           GetNavigatorFor()->LocateGlobalPointWithinVolume( CurrentF_Point );
350                                                   350 
351           G4double stepLengthFB;                  351           G4double stepLengthFB;
352           G4ThreeVector PointH;                   352           G4ThreeVector PointH;
353           G4bool usedNavigatorFB = false;         353           G4bool usedNavigatorFB = false; 
354                                                   354 
355           // Check whether any volumes are enc    355           // Check whether any volumes are encountered by the chord FB
356           // ---------------------------------    356           // ---------------------------------------------------------
357                                                   357 
358           G4bool Intersects_FB = IntersectChor    358           G4bool Intersects_FB = IntersectChord( CurrentF_Point, Point_B, 
359                                                   359                                                  NewSafety,fPreviousSafety,
360                                                   360                                                  fPreviousSftOrigin,
361                                                   361                                                  stepLengthFB,
362                                                   362                                                  PointH,
363                                                   363                                                  &usedNavigatorFB);
364           if( Intersects_FB )                     364           if( Intersects_FB )
365           {                                       365           { 
366             // There is an intersection of FB     366             // There is an intersection of FB with a volume boundary
367             // H <- First Intersection of Chor    367             // H <- First Intersection of Chord FB 
368                                                   368 
369             // H is our new Candidate for the     369             // H is our new Candidate for the intersection point.
370             // It replaces  "E" and we will re    370             // It replaces  "E" and we will repeat the test to see if
371             // it is a good enough approximate    371             // it is a good enough approximate point for us.
372                                                   372 
373             // Note that F must be in volume v    373             // Note that F must be in volume volA  (the same as A)
374             // (otherwise AF would meet a volu    374             // (otherwise AF would meet a volume boundary!)
375             //   A    <- F                        375             //   A    <- F 
376             //   E    <- H                        376             //   E    <- H
377             //                                    377             //
378             G4FieldTrack InterMed = ApproxInte    378             G4FieldTrack InterMed = ApproxIntersecPointV;
379             ApproxIntersecPointV = GetChordFin    379             ApproxIntersecPointV = GetChordFinderFor()->ApproxCurvePointS(
380                           CurrentA_PointVeloci    380                           CurrentA_PointVelocity,CurrentB_PointVelocity,
381                           InterMed,CurrentE_Po    381                           InterMed,CurrentE_Point,CurrentF_Point,PointH,
382                           false,GetEpsilonStep    382                           false,GetEpsilonStepFor());
383             CurrentA_PointVelocity = InterMed;    383             CurrentA_PointVelocity = InterMed;
384             CurrentE_Point = PointH;              384             CurrentE_Point = PointH;
385                                                   385 
386             // Need to recalculate the Exit No    386             // Need to recalculate the Exit Normal at the new PointG
387             //                                    387             //
388             G4bool validNormalLast;               388             G4bool validNormalLast; 
389             NormalAtEntry = GetSurfaceNormal(     389             NormalAtEntry = GetSurfaceNormal( PointH, validNormalLast ); 
390             validNormalAtE = validNormalLast;     390             validNormalAtE = validNormalLast;
391           }                                       391           }
392           else  // not Intersects_FB              392           else  // not Intersects_FB
393           {                                       393           {
394             // There is NO intersection of FB     394             // There is NO intersection of FB with a volume boundary
395                                                   395 
396             if( fin_section_depth[depth]  )       396             if( fin_section_depth[depth]  )
397             {                                     397             { 
398               // If B is the original endpoint    398               // If B is the original endpoint, this means that whatever
399               // volume(s) intersected the ori    399               // volume(s) intersected the original chord, none touch the
400               // smaller chords we have used.     400               // smaller chords we have used.
401               // The value of 'IntersectedOrRe    401               // The value of 'IntersectedOrRecalculatedFT' returned is
402               // likely not valid                 402               // likely not valid 
403                                                   403 
404               // Check on real final_section o    404               // Check on real final_section or SubEndSection
405               //                                  405               //
406               if( ((Second_half)&&(depth==0))     406               if( ((Second_half)&&(depth==0)) || (first_section) )
407               {                                   407               {
408                 there_is_no_intersection = tru    408                 there_is_no_intersection = true;   // real final_section
409               }                                   409               }
410               else                                410               else
411               {                                   411               {
412                 // end of subsection, not real    412                 // end of subsection, not real final section 
413                 // exit from the and go to the    413                 // exit from the and go to the depth-1 level 
414                                                   414 
415                 substep_no_p = param_substeps+    415                 substep_no_p = param_substeps+2;  // exit from the loop
416                                                   416 
417                 // but 'Second_half' is still     417                 // but 'Second_half' is still true because we need to find
418                 // the 'CurrentE_point' for th    418                 // the 'CurrentE_point' for the next loop
419                 //                                419                 //
420                 Second_half = true;               420                 Second_half = true; 
421                 sub_final_section = true;         421                 sub_final_section = true;
422               }                                   422               }
423             }                                     423             }
424             else                                  424             else
425             {                                     425             {
426               if( depth==0 )                      426               if( depth==0 )
427               {                                   427               {
428                 // We must restore the origina    428                 // We must restore the original endpoint
429                 //                                429                 //
430                 CurrentA_PointVelocity = Curre    430                 CurrentA_PointVelocity = CurrentB_PointVelocity;  // Got to B
431                 CurrentB_PointVelocity = Curve    431                 CurrentB_PointVelocity = CurveEndPointVelocity;
432                 SubStart_PointVelocity = Curre    432                 SubStart_PointVelocity = CurrentA_PointVelocity;
433                 ApproxIntersecPointV = GetChor    433                 ApproxIntersecPointV = GetChordFinderFor()
434                                ->ApproxCurvePo    434                                ->ApproxCurvePointV( CurrentA_PointVelocity, 
435                                                   435                                                     CurrentB_PointVelocity, 
436                                                   436                                                     CurrentE_Point,
437                                                   437                                                     GetEpsilonStepFor());
438                                                   438 
439                 restoredFullEndpoint = true;      439                 restoredFullEndpoint = true;
440                 ++restartB; // counter            440                 ++restartB; // counter
441               }                                   441               }
442               else                                442               else
443               {                                   443               {
444                 // We must restore the depth e    444                 // We must restore the depth endpoint
445                 //                                445                 //
446                 CurrentA_PointVelocity = Curre    446                 CurrentA_PointVelocity = CurrentB_PointVelocity;  // Got to B
447                 CurrentB_PointVelocity =  *ptr    447                 CurrentB_PointVelocity =  *ptrInterMedFT[depth];
448                 SubStart_PointVelocity = Curre    448                 SubStart_PointVelocity = CurrentA_PointVelocity;
449                 ApproxIntersecPointV = GetChor    449                 ApproxIntersecPointV = GetChordFinderFor()
450                                ->ApproxCurvePo    450                                ->ApproxCurvePointV( CurrentA_PointVelocity, 
451                                                   451                                                     CurrentB_PointVelocity, 
452                                                   452                                                     CurrentE_Point,
453                                                   453                                                     GetEpsilonStepFor());
454                 restoredFullEndpoint = true;      454                 restoredFullEndpoint = true;
455                 ++restartB; // counter            455                 ++restartB; // counter
456               }                                   456               }
457             }                                     457             }
458           } // Endif (Intersects_FB)              458           } // Endif (Intersects_FB)
459         } // Endif (Intersects_AF)                459         } // Endif (Intersects_AF)
460                                                   460 
461         // Ensure that the new endpoints are n    461         // Ensure that the new endpoints are not further apart in space
462         // than on the curve due to different     462         // than on the curve due to different errors in the integration
463         //                                        463         //
464         G4double linDistSq, curveDist;            464         G4double linDistSq, curveDist; 
465         linDistSq = ( CurrentB_PointVelocity.G    465         linDistSq = ( CurrentB_PointVelocity.GetPosition() 
466                     - CurrentA_PointVelocity.G    466                     - CurrentA_PointVelocity.GetPosition() ).mag2(); 
467         curveDist = CurrentB_PointVelocity.Get    467         curveDist = CurrentB_PointVelocity.GetCurveLength()
468                     - CurrentA_PointVelocity.G    468                     - CurrentA_PointVelocity.GetCurveLength();
469                                                   469 
470         // Change this condition for very stri    470         // Change this condition for very strict parameters of propagation 
471         //                                        471         //
472         if( curveDist*curveDist*(1+2* GetEpsil    472         if( curveDist*curveDist*(1+2* GetEpsilonStepFor()) < linDistSq )
473         {                                         473         {
474           // Re-integrate to obtain a new B       474           // Re-integrate to obtain a new B
475           //                                      475           //
476           G4FieldTrack newEndPointFT=             476           G4FieldTrack newEndPointFT=
477                   ReEstimateEndpoint( CurrentA    477                   ReEstimateEndpoint( CurrentA_PointVelocity,
478                                       CurrentB    478                                       CurrentB_PointVelocity,
479                                       linDistS    479                                       linDistSq,    // to avoid recalculation
480                                       curveDis    480                                       curveDist );
481           G4FieldTrack oldPointVelB = CurrentB    481           G4FieldTrack oldPointVelB = CurrentB_PointVelocity; 
482           CurrentB_PointVelocity = newEndPoint    482           CurrentB_PointVelocity = newEndPointFT;
483                                                   483          
484           if ( (fin_section_depth[depth])         484           if ( (fin_section_depth[depth])           // real final section
485              &&( first_section  || ((Second_ha    485              &&( first_section  || ((Second_half)&&(depth==0)) ) )
486           {                                       486           {
487             recalculatedEndPoint = true;          487             recalculatedEndPoint = true;
488             IntersectedOrRecalculatedFT = newE    488             IntersectedOrRecalculatedFT = newEndPointFT;
489               // So that we can return it, if     489               // So that we can return it, if it is the endpoint!
490           }                                       490           }
491         }                                         491         }
492         if( curveDist < 0.0 )                     492         if( curveDist < 0.0 )
493         {                                         493         {
494           fVerboseLevel = 5; // Print out a ma    494           fVerboseLevel = 5; // Print out a maximum of information
495           printStatus( CurrentA_PointVelocity,    495           printStatus( CurrentA_PointVelocity,  CurrentB_PointVelocity,
496                        -1.0, NewSafety,  subst    496                        -1.0, NewSafety,  substep_no );
497           std::ostringstream message;             497           std::ostringstream message;
498           message << "Error in advancing propa    498           message << "Error in advancing propagation." << G4endl
499                   << "        Error in advanci    499                   << "        Error in advancing propagation." << G4endl
500                   << "        Point A (start)     500                   << "        Point A (start) is " << CurrentA_PointVelocity
501                   << G4endl                       501                   << G4endl
502                   << "        Point B (end)       502                   << "        Point B (end)   is " << CurrentB_PointVelocity
503                   << G4endl                       503                   << G4endl
504                   << "        Curve distance i    504                   << "        Curve distance is " << curveDist << G4endl
505                   << G4endl                       505                   << G4endl
506                   << "The final curve point is    506                   << "The final curve point is not further along"
507                   << " than the original!" <<     507                   << " than the original!" << G4endl;
508                                                   508 
509           if( recalculatedEndPoint )              509           if( recalculatedEndPoint )
510           {                                       510           {
511             message << "Recalculation of EndPo    511             message << "Recalculation of EndPoint was called with fEpsStep= "
512                     << GetEpsilonStepFor() <<     512                     << GetEpsilonStepFor() << G4endl;
513           }                                       513           }
514           oldprc = G4cerr.precision(20);          514           oldprc = G4cerr.precision(20);
515           message << " Point A (Curve start)      515           message << " Point A (Curve start)     is " << CurveStartPointVelocity
516                   << G4endl                       516                   << G4endl
517                   << " Point B (Curve   end)      517                   << " Point B (Curve   end)     is " << CurveEndPointVelocity
518                   << G4endl                       518                   << G4endl
519                   << " Point A (Current start)    519                   << " Point A (Current start)   is " << CurrentA_PointVelocity
520                   << G4endl                       520                   << G4endl
521                   << " Point B (Current end)      521                   << " Point B (Current end)     is " << CurrentB_PointVelocity
522                   << G4endl                       522                   << G4endl
523                   << " Point S (Sub start)        523                   << " Point S (Sub start)       is " << SubStart_PointVelocity
524                   << G4endl                       524                   << G4endl
525                   << " Point E (Trial Point)      525                   << " Point E (Trial Point)     is " << CurrentE_Point
526                   << G4endl                       526                   << G4endl
527                   << " Old Point F(Intersectio    527                   << " Old Point F(Intersection) is " << CurrentF_Point
528                   << G4endl                       528                   << G4endl
529                   << " New Point F(Intersectio    529                   << " New Point F(Intersection) is " << ApproxIntersecPointV
530                   << G4endl                       530                   << G4endl
531                   << "        LocateIntersecti    531                   << "        LocateIntersection parameters are : Substep no= "
532                   << substep_no << G4endl         532                   << substep_no << G4endl
533                   << "        Substep depth no    533                   << "        Substep depth no= "<< substep_no_p  << " Depth= "
534                   << depth << G4endl              534                   << depth << G4endl
535                   << "        Restarted no= "<    535                   << "        Restarted no= "<< restartB  << " Epsilon= "
536                   << GetEpsilonStepFor() <<" D    536                   << GetEpsilonStepFor() <<" DeltaInters= "
537                   << GetDeltaIntersectionFor()    537                   << GetDeltaIntersectionFor();
538           G4cerr.precision( oldprc );             538           G4cerr.precision( oldprc ); 
539                                                   539 
540           G4Exception("G4BrentLocator::Estimat    540           G4Exception("G4BrentLocator::EstimateIntersectionPoint()",
541                       "GeomNav0003", FatalExce    541                       "GeomNav0003", FatalException, message);
542         }                                         542         }
543                                                   543 
544         if( restoredFullEndpoint )                544         if( restoredFullEndpoint )
545         {                                         545         {
546           fin_section_depth[depth] = restoredF    546           fin_section_depth[depth] = restoredFullEndpoint;
547           restoredFullEndpoint = false;           547           restoredFullEndpoint = false;
548         }                                         548         }
549       } // EndIf ( E is close enough to the cu    549       } // EndIf ( E is close enough to the curve, ie point F. )
550         // tests ChordAF_Vector.mag() <= maxim    550         // tests ChordAF_Vector.mag() <= maximum_lateral_displacement 
551                                                   551 
552 #ifdef G4DEBUG_LOCATE_INTERSECTION                552 #ifdef G4DEBUG_LOCATE_INTERSECTION  
553       G4int trigger_substepno_print= warn_subs    553       G4int trigger_substepno_print= warn_substeps - 20 ;
554                                                   554 
555       if( substep_no >= trigger_substepno_prin    555       if( substep_no >= trigger_substepno_print )
556       {                                           556       {
557         G4cout << "Difficulty in converging in    557         G4cout << "Difficulty in converging in "
558                << "G4BrentLocator::EstimateInt    558                << "G4BrentLocator::EstimateIntersectionPoint()"
559                << G4endl                          559                << G4endl
560                << "    Substep no = " << subst    560                << "    Substep no = " << substep_no << G4endl;
561         if( substep_no == trigger_substepno_pr    561         if( substep_no == trigger_substepno_print )
562         {                                         562         {
563           printStatus( CurveStartPointVelocity    563           printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
564                        -1.0, NewSafety, 0);       564                        -1.0, NewSafety, 0);
565         }                                         565         }
566         G4cout << " State of point A: ";          566         G4cout << " State of point A: "; 
567         printStatus( CurrentA_PointVelocity, C    567         printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
568                      -1.0, NewSafety, substep_    568                      -1.0, NewSafety, substep_no-1, 0);
569         G4cout << " State of point B: ";          569         G4cout << " State of point B: "; 
570         printStatus( CurrentA_PointVelocity, C    570         printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
571                      -1.0, NewSafety, substep_    571                      -1.0, NewSafety, substep_no);
572       }                                           572       }
573 #endif                                            573 #endif
574       ++substep_no;                               574       ++substep_no; 
575       ++substep_no_p;                             575       ++substep_no_p;
576                                                   576 
577     } while (  ( ! found_approximate_intersect    577     } while (  ( ! found_approximate_intersection )
578             && ( ! there_is_no_intersection )     578             && ( ! there_is_no_intersection )     
579             && ( substep_no_p <= param_substep    579             && ( substep_no_p <= param_substeps) );  // UNTIL found or
580                                                   580                                                      // failed param substep
581     first_section = false;                        581     first_section = false;
582                                                   582 
583     if( (!found_approximate_intersection) && (    583     if( (!found_approximate_intersection) && (!there_is_no_intersection) )
584     {                                             584     {
585       G4double did_len = std::abs( CurrentA_Po    585       G4double did_len = std::abs( CurrentA_PointVelocity.GetCurveLength()
586                        - SubStart_PointVelocit    586                        - SubStart_PointVelocity.GetCurveLength()); 
587       G4double all_len = std::abs( CurrentB_Po    587       G4double all_len = std::abs( CurrentB_PointVelocity.GetCurveLength()
588                        - SubStart_PointVelocit    588                        - SubStart_PointVelocity.GetCurveLength());
589                                                   589    
590       G4double stepLengthAB;                      590       G4double stepLengthAB;
591       G4ThreeVector PointGe;                      591       G4ThreeVector PointGe;
592                                                   592 
593       // Check if progress is too slow and if     593       // Check if progress is too slow and if it possible to go deeper,
594       // then halve the step if so                594       // then halve the step if so
595       //                                          595       //
596       if ( ( did_len < fraction_done*all_len )    596       if ( ( did_len < fraction_done*all_len )
597         && (depth < max_depth) && (!sub_final_    597         && (depth < max_depth) && (!sub_final_section) )
598       {                                           598       {
599         Second_half=false;                        599         Second_half=false;
600         ++depth;                                  600         ++depth;
601                                                   601 
602         G4double Sub_len = (all_len-did_len)/(    602         G4double Sub_len = (all_len-did_len)/(2.);
603         G4FieldTrack start = CurrentA_PointVel    603         G4FieldTrack start = CurrentA_PointVelocity;
604         auto integrDriver =                       604         auto integrDriver =
605                          GetChordFinderFor()->    605                          GetChordFinderFor()->GetIntegrationDriver();
606         integrDriver->AccurateAdvance(start, S    606         integrDriver->AccurateAdvance(start, Sub_len, GetEpsilonStepFor());
607         *ptrInterMedFT[depth] = start;            607         *ptrInterMedFT[depth] = start;
608         CurrentB_PointVelocity = *ptrInterMedF    608         CurrentB_PointVelocity = *ptrInterMedFT[depth];
609                                                   609  
610         // Adjust 'SubStartPoint' to calculate    610         // Adjust 'SubStartPoint' to calculate the 'did_length' in next loop
611         //                                        611         //
612         SubStart_PointVelocity = CurrentA_Poin    612         SubStart_PointVelocity = CurrentA_PointVelocity;
613                                                   613 
614         // Find new trial intersection point n    614         // Find new trial intersection point needed at start of the loop
615         //                                        615         //
616         G4ThreeVector Point_A = CurrentA_Point    616         G4ThreeVector Point_A = CurrentA_PointVelocity.GetPosition();
617         G4ThreeVector SubE_point = CurrentB_Po    617         G4ThreeVector SubE_point = CurrentB_PointVelocity.GetPosition();   
618                                                   618      
619         GetNavigatorFor()->LocateGlobalPointWi    619         GetNavigatorFor()->LocateGlobalPointWithinVolume(Point_A);
620         G4bool Intersects_AB = IntersectChord(    620         G4bool Intersects_AB = IntersectChord(Point_A, SubE_point,
621                                                   621                                               NewSafety, fPreviousSafety,
622                                                   622                                               fPreviousSftOrigin,stepLengthAB,
623                                                   623                                               PointGe);
624         if( Intersects_AB )                       624         if( Intersects_AB )
625         {                                         625         {
626           last_AF_intersection = Intersects_AB    626           last_AF_intersection = Intersects_AB;
627           CurrentE_Point = PointGe;               627           CurrentE_Point = PointGe;
628           fin_section_depth[depth] = true;        628           fin_section_depth[depth] = true;
629                                                   629 
630           // Need to recalculate the Exit Norm    630           // Need to recalculate the Exit Normal at the new PointG
631           //                                      631           //
632           G4bool validNormalAB;                   632           G4bool validNormalAB; 
633           NormalAtEntry = GetSurfaceNormal( Po    633           NormalAtEntry = GetSurfaceNormal( PointGe, validNormalAB ); 
634           validNormalAtE = validNormalAB;         634           validNormalAtE = validNormalAB;  
635         }                                         635         }
636         else                                      636         else
637         {                                         637         {
638           // No intersection found for first p    638           // No intersection found for first part of curve
639           // (CurrentA,InterMedPoint[depth]).     639           // (CurrentA,InterMedPoint[depth]). Go to the second part
640           //                                      640           //
641           Second_half = true;                     641           Second_half = true;
642         }                                         642         }
643       } // if did_len                             643       } // if did_len
644                                                   644 
645       if( (Second_half)&&(depth!=0) )             645       if( (Second_half)&&(depth!=0) )
646       {                                           646       {
647         // Second part of curve (InterMed[dept    647         // Second part of curve (InterMed[depth],Intermed[depth-1])                       ) 
648         // On the depth-1 level normally we ar    648         // On the depth-1 level normally we are on the 'second_half'
649                                                   649 
650         Second_half = true;                       650         Second_half = true;
651                                                   651 
652         //  Find new trial intersection point     652         //  Find new trial intersection point needed at start of the loop
653         //                                        653         //
654         SubStart_PointVelocity = *ptrInterMedF    654         SubStart_PointVelocity = *ptrInterMedFT[depth];
655         CurrentA_PointVelocity = *ptrInterMedF    655         CurrentA_PointVelocity = *ptrInterMedFT[depth];
656         CurrentB_PointVelocity = *ptrInterMedF    656         CurrentB_PointVelocity = *ptrInterMedFT[depth-1];
657          // Ensure that the new endpoints are     657          // Ensure that the new endpoints are not further apart in space
658         // than on the curve due to different     658         // than on the curve due to different errors in the integration
659         //                                        659         //
660         G4double linDistSq, curveDist;            660         G4double linDistSq, curveDist; 
661         linDistSq = ( CurrentB_PointVelocity.G    661         linDistSq = ( CurrentB_PointVelocity.GetPosition() 
662                     - CurrentA_PointVelocity.G    662                     - CurrentA_PointVelocity.GetPosition() ).mag2(); 
663         curveDist = CurrentB_PointVelocity.Get    663         curveDist = CurrentB_PointVelocity.GetCurveLength()
664                     - CurrentA_PointVelocity.G    664                     - CurrentA_PointVelocity.GetCurveLength();
665         if( curveDist*curveDist*(1+2*GetEpsilo    665         if( curveDist*curveDist*(1+2*GetEpsilonStepFor() ) < linDistSq )
666         {                                         666         {
667           // Re-integrate to obtain a new B       667           // Re-integrate to obtain a new B
668           //                                      668           //
669           G4FieldTrack newEndPointFT =            669           G4FieldTrack newEndPointFT =
670                   ReEstimateEndpoint( CurrentA    670                   ReEstimateEndpoint( CurrentA_PointVelocity,
671                                       CurrentB    671                                       CurrentB_PointVelocity,
672                                       linDistS    672                                       linDistSq,    // to avoid recalculation
673                                       curveDis    673                                       curveDist );
674           G4FieldTrack oldPointVelB = CurrentB    674           G4FieldTrack oldPointVelB = CurrentB_PointVelocity; 
675           CurrentB_PointVelocity = newEndPoint    675           CurrentB_PointVelocity = newEndPointFT;
676           if ( depth==1 )                         676           if ( depth==1 )
677           {                                       677           {
678             recalculatedEndPoint = true;          678             recalculatedEndPoint = true;
679             IntersectedOrRecalculatedFT = newE    679             IntersectedOrRecalculatedFT = newEndPointFT;
680             // So that we can return it, if it    680             // So that we can return it, if it is the endpoint!
681           }                                       681           }
682         }                                         682         }
683                                                   683 
684                                                   684 
685         G4ThreeVector Point_A    = CurrentA_Po    685         G4ThreeVector Point_A    = CurrentA_PointVelocity.GetPosition();
686         G4ThreeVector SubE_point = CurrentB_Po    686         G4ThreeVector SubE_point = CurrentB_PointVelocity.GetPosition();   
687         GetNavigatorFor()->LocateGlobalPointWi    687         GetNavigatorFor()->LocateGlobalPointWithinVolume(Point_A);
688         G4bool Intersects_AB = IntersectChord(    688         G4bool Intersects_AB = IntersectChord(Point_A, SubE_point, NewSafety,
689                                                   689                                               fPreviousSafety,
690                                                   690                                                fPreviousSftOrigin,stepLengthAB, PointGe);
691         if( Intersects_AB )                       691         if( Intersects_AB )
692         {                                         692         {
693           last_AF_intersection = Intersects_AB    693           last_AF_intersection = Intersects_AB;
694           CurrentE_Point = PointGe;               694           CurrentE_Point = PointGe;
695                                                   695 
696           G4bool validNormalAB;                   696           G4bool validNormalAB; 
697           NormalAtEntry  = GetSurfaceNormal( P    697           NormalAtEntry  = GetSurfaceNormal( PointGe, validNormalAB ); 
698           validNormalAtE = validNormalAB;         698           validNormalAtE = validNormalAB;
699         }                                         699         }
700                                                   700        
701         depth--;                                  701         depth--;
702         fin_section_depth[depth]=true;            702         fin_section_depth[depth]=true;
703       }                                           703       }
704     }  // if(!found_aproximate_intersection)      704     }  // if(!found_aproximate_intersection)
705                                                   705 
706   } while ( ( ! found_approximate_intersection    706   } while ( ( ! found_approximate_intersection )
707             && ( ! there_is_no_intersection )     707             && ( ! there_is_no_intersection )     
708             && ( substep_no <= max_substeps) )    708             && ( substep_no <= max_substeps) ); // UNTIL found or failed
709                                                   709 
710   if( substep_no > max_no_seen )                  710   if( substep_no > max_no_seen )
711   {                                               711   {
712     max_no_seen = substep_no;                     712     max_no_seen = substep_no; 
713 #ifdef G4DEBUG_LOCATE_INTERSECTION                713 #ifdef G4DEBUG_LOCATE_INTERSECTION
714     if( max_no_seen > warn_substeps )             714     if( max_no_seen > warn_substeps )
715     {                                             715     {
716       trigger_substepno_print = max_no_seen-20    716       trigger_substepno_print = max_no_seen-20; // Want to see last 20 steps 
717     }                                             717     }
718 #endif                                            718 #endif
719   }                                               719   }
720                                                   720 
721   if(  ( substep_no >= max_substeps)              721   if(  ( substep_no >= max_substeps)
722       && !there_is_no_intersection                722       && !there_is_no_intersection
723       && !found_approximate_intersection )        723       && !found_approximate_intersection )
724   {                                               724   {
725     G4cout << "ERROR - G4BrentLocator::Estimat    725     G4cout << "ERROR - G4BrentLocator::EstimateIntersectionPoint()" << G4endl
726            << "        Start and end-point of     726            << "        Start and end-point of requested Step:" << G4endl;
727     printStatus( CurveStartPointVelocity, Curv    727     printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
728                  -1.0, NewSafety, 0);             728                  -1.0, NewSafety, 0);
729     G4cout << "        Start and end-point of     729     G4cout << "        Start and end-point of current Sub-Step:" << G4endl;
730     printStatus( CurrentA_PointVelocity, Curre    730     printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
731                  -1.0, NewSafety, substep_no-1    731                  -1.0, NewSafety, substep_no-1);
732     printStatus( CurrentA_PointVelocity, Curre    732     printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
733                  -1.0, NewSafety, substep_no);    733                  -1.0, NewSafety, substep_no);
734     std::ostringstream message;                   734     std::ostringstream message;
735     message << "Too many substeps!" << G4endl     735     message << "Too many substeps!" << G4endl
736             << "          Convergence is requi    736             << "          Convergence is requiring too many substeps: "
737             << substep_no << G4endl               737             << substep_no << G4endl
738             << "          Abandoning effort to    738             << "          Abandoning effort to intersect. " << G4endl
739             << "          Found intersection =    739             << "          Found intersection = "
740             << found_approximate_intersection     740             << found_approximate_intersection << G4endl
741             << "          Intersection exists     741             << "          Intersection exists = "
742             << !there_is_no_intersection << G4    742             << !there_is_no_intersection << G4endl;
743     oldprc = G4cout.precision( 10 );              743     oldprc = G4cout.precision( 10 ); 
744     G4double done_len = CurrentA_PointVelocity    744     G4double done_len = CurrentA_PointVelocity.GetCurveLength(); 
745     G4double full_len = CurveEndPointVelocity.    745     G4double full_len = CurveEndPointVelocity.GetCurveLength();
746     message << "        Undertaken only length    746     message << "        Undertaken only length: " << done_len
747             << " out of " << full_len << " req    747             << " out of " << full_len << " required." << G4endl
748             << "        Remaining length = " <    748             << "        Remaining length = " << full_len - done_len; 
749     G4cout.precision( oldprc );                   749     G4cout.precision( oldprc ); 
750                                                   750 
751     G4Exception("G4BrentLocator::EstimateInter    751     G4Exception("G4BrentLocator::EstimateIntersectionPoint()",
752                 "GeomNav0003", FatalException,    752                 "GeomNav0003", FatalException, message);
753   }                                               753   }
754   else if( substep_no >= warn_substeps )          754   else if( substep_no >= warn_substeps )
755   {                                               755   {  
756     oldprc = G4cout.precision( 10 );              756     oldprc = G4cout.precision( 10 ); 
757     std::ostringstream message;                   757     std::ostringstream message;
758     message << "Many substeps while trying to     758     message << "Many substeps while trying to locate intersection."
759             << G4endl                             759             << G4endl
760             << "          Undertaken length: "    760             << "          Undertaken length: "  
761             << CurrentB_PointVelocity.GetCurve    761             << CurrentB_PointVelocity.GetCurveLength()
762             << " - Needed: "  << substep_no <<    762             << " - Needed: "  << substep_no << " substeps." << G4endl
763             << "          Warning level = " <<    763             << "          Warning level = " << warn_substeps
764             << " and maximum substeps = " << m    764             << " and maximum substeps = " << max_substeps;
765     G4Exception("G4BrentLocator::EstimateInter    765     G4Exception("G4BrentLocator::EstimateIntersectionPoint()",
766                 "GeomNav1002", JustWarning, me    766                 "GeomNav1002", JustWarning, message);
767     G4cout.precision( oldprc );                   767     G4cout.precision( oldprc ); 
768   }                                               768   }
769   return  !there_is_no_intersection; //  Succe    769   return  !there_is_no_intersection; //  Success or failure
770 }                                                 770 }
771                                                   771