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 9.4)


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