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