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.2.p3)


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