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