Geant4  10.01.p03
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 90836 2015-06-10 09:31:06Z 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  //
56 
57  // Counter for Number Of Calls to ReIntegrationEndPoint Method
58  //
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="
74  G4cout << "G4BrentLocator::ReIntegrateEndPoint was called "
76  << " times and for depth algorithm "
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  static 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
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  do // REPEAT param
228  {
229  G4ThreeVector Point_A = CurrentA_PointVelocity.GetPosition();
230  G4ThreeVector Point_B = CurrentB_PointVelocity.GetPosition();
231 
232  // F = a point on true AB path close to point E
233  // (the closest if possible)
234  //
235  if(substep_no_p==0)
236  {
237  ApproxIntersecPointV = GetChordFinderFor()
238  ->ApproxCurvePointV( CurrentA_PointVelocity,
239  CurrentB_PointVelocity,
240  CurrentE_Point,
242  // The above method is the key & most intuitive part ...
243  }
244 #ifdef G4DEBUG_FIELD
245  if( ApproxIntersecPointV.GetCurveLength() >
246  CurrentB_PointVelocity.GetCurveLength() * (1.0 + tolerance) )
247  {
248  G4Exception("G4BrentLocator::EstimateIntersectionPoint()",
249  "GeomNav0003", FatalException,
250  "Intermediate F point is past end B point" );
251  }
252 #endif
253 
254  G4ThreeVector CurrentF_Point= ApproxIntersecPointV.GetPosition();
255 
256  // First check whether EF is small - then F is a good approx. point
257  // Calculate the length and direction of the chord AF
258  //
259  G4ThreeVector ChordEF_Vector = CurrentF_Point - CurrentE_Point;
260  G4ThreeVector NewMomentumDir= ApproxIntersecPointV.GetMomentumDir();
261  G4double MomDir_dot_Norm= NewMomentumDir.dot( NormalAtEntry ) ;
262 
263 #ifdef G4DEBUG_FIELD
264  G4ThreeVector ChordAB = Point_B - Point_A;
265 
266  G4VIntersectionLocator::ReportTrialStep( substep_no, ChordAB,
267  ChordEF_Vector, NewMomentumDir, NormalAtEntry, validNormalAtE );
268 #endif
269 
270  G4bool adequate_angle;
271  adequate_angle = ( MomDir_dot_Norm >= 0.0 ) // Can use -epsilon instead.
272  || (! validNormalAtE ); // Makes criterion invalid
273  G4double EF_dist2 = ChordEF_Vector.mag2();
274  if ( ( EF_dist2 <= sqr(fiDeltaIntersection) && ( adequate_angle ) )
275  || ( EF_dist2 <= kCarTolerance*kCarTolerance ) )
276  {
277  found_approximate_intersection = true;
278 
279  // Create the "point" return value
280  //
281  IntersectedOrRecalculatedFT = ApproxIntersecPointV;
282  IntersectedOrRecalculatedFT.SetPosition( CurrentE_Point );
283 
285  {
286  // Try to Get Correction of IntersectionPoint using SurfaceNormal()
287  //
288  G4ThreeVector IP;
289  G4ThreeVector MomentumDir=ApproxIntersecPointV.GetMomentumDirection();
290  G4bool goodCorrection = AdjustmentOfFoundIntersection( Point_A,
291  CurrentE_Point, CurrentF_Point, MomentumDir,
292  last_AF_intersection, IP, NewSafety,
293  fPreviousSafety, fPreviousSftOrigin );
294  if ( goodCorrection )
295  {
296  IntersectedOrRecalculatedFT = ApproxIntersecPointV;
297  IntersectedOrRecalculatedFT.SetPosition(IP);
298  }
299  }
300 
301  // Note: in order to return a point on the boundary,
302  // we must return E. But it is F on the curve.
303  // So we must "cheat": we are using the position at point E
304  // and the velocity at point F !!!
305  //
306  // This must limit the length we can allow for displacement!
307  }
308  else // E is NOT close enough to the curve (ie point F)
309  {
310  // Check whether any volumes are encountered by the chord AF
311  // ---------------------------------------------------------
312  // First relocate to restore any Voxel etc information
313  // in the Navigator before calling ComputeStep()
314  //
316 
317  G4ThreeVector PointG; // Candidate intersection point
318  G4double stepLengthAF;
319  G4bool usedNavigatorAF = false;
320  G4bool Intersects_AF = IntersectChord( Point_A, CurrentF_Point,
321  NewSafety,fPreviousSafety,
322  fPreviousSftOrigin,
323  stepLengthAF,
324  PointG,
325  &usedNavigatorAF);
326  last_AF_intersection = Intersects_AF;
327  if( Intersects_AF )
328  {
329  // G is our new Candidate for the intersection point.
330  // It replaces "E" and we will repeat the test to see if
331  // it is a good enough approximate point for us.
332  // B <- F
333  // E <- G
334  //
335  G4FieldTrack EndPoint = ApproxIntersecPointV;
336  ApproxIntersecPointV = GetChordFinderFor()->ApproxCurvePointS(
337  CurrentA_PointVelocity, CurrentB_PointVelocity,
338  EndPoint,CurrentE_Point, CurrentF_Point,PointG,
339  true, GetEpsilonStepFor() );
340  CurrentB_PointVelocity = EndPoint;
341  CurrentE_Point = PointG;
342 
343  // Need to recalculate the Exit Normal at the new PointG
344  // Know that a call was made to Navigator::ComputeStep in
345  // IntersectChord above.
346  //
347  G4bool validNormalLast;
348  NormalAtEntry = GetSurfaceNormal( PointG, validNormalLast );
349  validNormalAtE = validNormalLast;
350 
351  // By moving point B, must take care if current
352  // AF has no intersection to try current FB!!
353  //
354  fin_section_depth[depth] = false;
355 #ifdef G4VERBOSE
356  if( fVerboseLevel > 3 )
357  {
358  G4cout << "G4PiF::LI> Investigating intermediate point"
359  << " at s=" << ApproxIntersecPointV.GetCurveLength()
360  << " on way to full s="
361  << CurveEndPointVelocity.GetCurveLength() << G4endl;
362  }
363 #endif
364  }
365  else // not Intersects_AF
366  {
367  // In this case:
368  // There is NO intersection of AF with a volume boundary.
369  // We must continue the search in the segment FB!
370  //
371  GetNavigatorFor()->LocateGlobalPointWithinVolume( CurrentF_Point );
372 
373  G4double stepLengthFB;
374  G4ThreeVector PointH;
375  G4bool usedNavigatorFB = false;
376 
377  // Check whether any volumes are encountered by the chord FB
378  // ---------------------------------------------------------
379 
380  G4bool Intersects_FB = IntersectChord( CurrentF_Point, Point_B,
381  NewSafety,fPreviousSafety,
382  fPreviousSftOrigin,
383  stepLengthFB,
384  PointH,
385  &usedNavigatorFB);
386  if( Intersects_FB )
387  {
388  // There is an intersection of FB with a volume boundary
389  // H <- First Intersection of Chord FB
390 
391  // H is our new Candidate for the intersection point.
392  // It replaces "E" and we will repeat the test to see if
393  // it is a good enough approximate point for us.
394 
395  // Note that F must be in volume volA (the same as A)
396  // (otherwise AF would meet a volume boundary!)
397  // A <- F
398  // E <- H
399  //
400  G4FieldTrack InterMed=ApproxIntersecPointV;
401  ApproxIntersecPointV = GetChordFinderFor()->ApproxCurvePointS(
402  CurrentA_PointVelocity,CurrentB_PointVelocity,
403  InterMed,CurrentE_Point,CurrentF_Point,PointH,
404  false,GetEpsilonStepFor());
405  CurrentA_PointVelocity = InterMed;
406  CurrentE_Point = PointH;
407 
408  // Need to recalculate the Exit Normal at the new PointG
409  //
410  G4bool validNormalLast;
411  NormalAtEntry = GetSurfaceNormal( PointH, validNormalLast );
412  validNormalAtE= validNormalLast;
413  }
414  else // not Intersects_FB
415  {
416  // There is NO intersection of FB with a volume boundary
417 
418  if( fin_section_depth[depth] )
419  {
420  // If B is the original endpoint, this means that whatever
421  // volume(s) intersected the original chord, none touch the
422  // smaller chords we have used.
423  // The value of 'IntersectedOrRecalculatedFT' returned is
424  // likely not valid
425 
426  // Check on real final_section or SubEndSection
427  //
428  if( ((Second_half)&&(depth==0)) || (first_section) )
429  {
430  there_is_no_intersection = true; // real final_section
431  }
432  else
433  {
434  // end of subsection, not real final section
435  // exit from the and go to the depth-1 level
436 
437  substep_no_p = param_substeps+2; // exit from the loop
438 
439  // but 'Second_half' is still true because we need to find
440  // the 'CurrentE_point' for the next loop
441  //
442  Second_half = true;
443  sub_final_section = true;
444  }
445  }
446  else
447  {
448  if(depth==0)
449  {
450  // We must restore the original endpoint
451  //
452  CurrentA_PointVelocity = CurrentB_PointVelocity; // Got to B
453  CurrentB_PointVelocity = CurveEndPointVelocity;
454  SubStart_PointVelocity = CurrentA_PointVelocity;
455  ApproxIntersecPointV = GetChordFinderFor()
456  ->ApproxCurvePointV( CurrentA_PointVelocity,
457  CurrentB_PointVelocity,
458  CurrentE_Point,
460 
461  restoredFullEndpoint = true;
462  restartB++; // counter
463  }
464  else
465  {
466  // We must restore the depth endpoint
467  //
468  CurrentA_PointVelocity = CurrentB_PointVelocity; // Got to B
469  CurrentB_PointVelocity = *ptrInterMedFT[depth];
470  SubStart_PointVelocity = CurrentA_PointVelocity;
471  ApproxIntersecPointV = GetChordFinderFor()
472  ->ApproxCurvePointV( CurrentA_PointVelocity,
473  CurrentB_PointVelocity,
474  CurrentE_Point,
476  restoredFullEndpoint = true;
477  restartB++; // counter
478  }
479  }
480  } // Endif (Intersects_FB)
481  } // Endif (Intersects_AF)
482 
483  // Ensure that the new endpoints are not further apart in space
484  // than on the curve due to different errors in the integration
485  //
486  G4double linDistSq, curveDist;
487  linDistSq = ( CurrentB_PointVelocity.GetPosition()
488  - CurrentA_PointVelocity.GetPosition() ).mag2();
489  curveDist = CurrentB_PointVelocity.GetCurveLength()
490  - CurrentA_PointVelocity.GetCurveLength();
491 
492  // Change this condition for very strict parameters of propagation
493  //
494  if( curveDist*curveDist*(1+2* GetEpsilonStepFor()) < linDistSq )
495  {
496  // Re-integrate to obtain a new B
497  //
498  G4FieldTrack newEndPointFT=
499  ReEstimateEndpoint( CurrentA_PointVelocity,
500  CurrentB_PointVelocity,
501  linDistSq, // to avoid recalculation
502  curveDist );
503  G4FieldTrack oldPointVelB = CurrentB_PointVelocity;
504  CurrentB_PointVelocity = newEndPointFT;
505 
506  if ( (fin_section_depth[depth]) // real final section
507  &&( first_section || ((Second_half)&&(depth==0)) ) )
508  {
509  recalculatedEndPoint = true;
510  IntersectedOrRecalculatedFT = newEndPointFT;
511  // So that we can return it, if it is the endpoint!
512  }
513  }
514  if( curveDist < 0.0 )
515  {
516  fVerboseLevel = 5; // Print out a maximum of information
517  printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
518  -1.0, NewSafety, substep_no );
519  std::ostringstream message;
520  message << "Error in advancing propagation." << G4endl
521  << " Error in advancing propagation." << G4endl
522  << " Point A (start) is " << CurrentA_PointVelocity
523  << G4endl
524  << " Point B (end) is " << CurrentB_PointVelocity
525  << G4endl
526  << " Curve distance is " << curveDist << G4endl
527  << G4endl
528  << "The final curve point is not further along"
529  << " than the original!" << G4endl;
530 
531  if( recalculatedEndPoint )
532  {
533  message << "Recalculation of EndPoint was called with fEpsStep= "
534  << GetEpsilonStepFor() << G4endl;
535  }
536  oldprc = G4cerr.precision(20);
537  message << " Point A (Curve start) is " << CurveStartPointVelocity
538  << G4endl
539  << " Point B (Curve end) is " << CurveEndPointVelocity
540  << G4endl
541  << " Point A (Current start) is " << CurrentA_PointVelocity
542  << G4endl
543  << " Point B (Current end) is " << CurrentB_PointVelocity
544  << G4endl
545  << " Point S (Sub start) is " << SubStart_PointVelocity
546  << G4endl
547  << " Point E (Trial Point) is " << CurrentE_Point
548  << G4endl
549  << " Old Point F(Intersection) is " << CurrentF_Point
550  << G4endl
551  << " New Point F(Intersection) is " << ApproxIntersecPointV
552  << G4endl
553  << " LocateIntersection parameters are : Substep no= "
554  << substep_no << G4endl
555  << " Substep depth no= "<< substep_no_p << " Depth= "
556  << depth << G4endl
557  << " Restarted no= "<< restartB << " Epsilon= "
558  << GetEpsilonStepFor() <<" DeltaInters= "
560  G4cerr.precision( oldprc );
561 
562  G4Exception("G4BrentLocator::EstimateIntersectionPoint()",
563  "GeomNav0003", FatalException, message);
564  }
565 
566  if(restoredFullEndpoint)
567  {
568  fin_section_depth[depth] = restoredFullEndpoint;
569  restoredFullEndpoint = false;
570  }
571  } // EndIf ( E is close enough to the curve, ie point F. )
572  // tests ChordAF_Vector.mag() <= maximum_lateral_displacement
573 
574 #ifdef G4DEBUG_LOCATE_INTERSECTION
575  static G4int trigger_substepno_print= warn_substeps - 20 ;
576 
577  if( substep_no >= trigger_substepno_print )
578  {
579  G4cout << "Difficulty in converging in "
580  << "G4BrentLocator::EstimateIntersectionPoint()"
581  << G4endl
582  << " Substep no = " << substep_no << G4endl;
583  if( substep_no == trigger_substepno_print )
584  {
585  printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
586  -1.0, NewSafety, 0);
587  }
588  G4cout << " State of point A: ";
589  printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
590  -1.0, NewSafety, substep_no-1, 0);
591  G4cout << " State of point B: ";
592  printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
593  -1.0, NewSafety, substep_no);
594  }
595 #endif
596  substep_no++;
597  substep_no_p++;
598 
599  } while ( ( ! found_approximate_intersection )
600  && ( ! there_is_no_intersection )
601  && ( substep_no_p <= param_substeps) ); // UNTIL found or
602  // failed param substep
603  first_section = false;
604 
605  if( (!found_approximate_intersection) && (!there_is_no_intersection) )
606  {
607  G4double did_len = std::abs( CurrentA_PointVelocity.GetCurveLength()
608  - SubStart_PointVelocity.GetCurveLength());
609  G4double all_len = std::abs( CurrentB_PointVelocity.GetCurveLength()
610  - SubStart_PointVelocity.GetCurveLength());
611 
612  G4double stepLengthAB;
613  G4ThreeVector PointGe;
614 
615  // Check if progress is too slow and if it possible to go deeper,
616  // then halve the step if so
617  //
618  if ( ( did_len < fraction_done*all_len )
619  && (depth<max_depth) && (!sub_final_section) )
620  {
621  Second_half=false;
622  depth++;
623 
624  G4double Sub_len = (all_len-did_len)/(2.);
625  G4FieldTrack start = CurrentA_PointVelocity;
626  G4MagInt_Driver* integrDriver =
628  integrDriver->AccurateAdvance(start, Sub_len, GetEpsilonStepFor());
629  *ptrInterMedFT[depth] = start;
630  CurrentB_PointVelocity = *ptrInterMedFT[depth];
631 
632  // Adjust 'SubStartPoint' to calculate the 'did_length' in next loop
633  //
634  SubStart_PointVelocity = CurrentA_PointVelocity;
635 
636  // Find new trial intersection point needed at start of the loop
637  //
638  G4ThreeVector Point_A = CurrentA_PointVelocity.GetPosition();
639  G4ThreeVector SubE_point = CurrentB_PointVelocity.GetPosition();
640 
642  G4bool Intersects_AB = IntersectChord(Point_A, SubE_point,
643  NewSafety, fPreviousSafety,
644  fPreviousSftOrigin,stepLengthAB,
645  PointGe);
646  if( Intersects_AB )
647  {
648  last_AF_intersection = Intersects_AB;
649  CurrentE_Point = PointGe;
650  fin_section_depth[depth]=true;
651 
652  // Need to recalculate the Exit Normal at the new PointG
653  //
654  G4bool validNormalAB;
655  NormalAtEntry = GetSurfaceNormal( PointGe, validNormalAB );
656  validNormalAtE= validNormalAB;
657  }
658  else
659  {
660  // No intersection found for first part of curve
661  // (CurrentA,InterMedPoint[depth]). Go to the second part
662  //
663  Second_half = true;
664  }
665  } // if did_len
666 
667  if( (Second_half)&&(depth!=0) )
668  {
669  // Second part of curve (InterMed[depth],Intermed[depth-1]) )
670  // On the depth-1 level normally we are on the 'second_half'
671 
672  Second_half = true;
673 
674  // Find new trial intersection point needed at start of the loop
675  //
676  SubStart_PointVelocity = *ptrInterMedFT[depth];
677  CurrentA_PointVelocity = *ptrInterMedFT[depth];
678  CurrentB_PointVelocity = *ptrInterMedFT[depth-1];
679  // Ensure that the new endpoints are not further apart in space
680  // than on the curve due to different errors in the integration
681  //
682  G4double linDistSq, curveDist;
683  linDistSq = ( CurrentB_PointVelocity.GetPosition()
684  - CurrentA_PointVelocity.GetPosition() ).mag2();
685  curveDist = CurrentB_PointVelocity.GetCurveLength()
686  - CurrentA_PointVelocity.GetCurveLength();
687  if( curveDist*curveDist*(1+2*GetEpsilonStepFor() ) < linDistSq )
688  {
689  // Re-integrate to obtain a new B
690  //
691  G4FieldTrack newEndPointFT=
692  ReEstimateEndpoint( CurrentA_PointVelocity,
693  CurrentB_PointVelocity,
694  linDistSq, // to avoid recalculation
695  curveDist );
696  G4FieldTrack oldPointVelB = CurrentB_PointVelocity;
697  CurrentB_PointVelocity = newEndPointFT;
698  if (depth==1)
699  {
700  recalculatedEndPoint = true;
701  IntersectedOrRecalculatedFT = newEndPointFT;
702  // So that we can return it, if it is the endpoint!
703  }
704  }
705 
706 
707  G4ThreeVector Point_A = CurrentA_PointVelocity.GetPosition();
708  G4ThreeVector SubE_point = CurrentB_PointVelocity.GetPosition();
710  G4bool Intersects_AB = IntersectChord(Point_A, SubE_point, NewSafety,
711  fPreviousSafety,
712  fPreviousSftOrigin,stepLengthAB, PointGe);
713  if( Intersects_AB )
714  {
715  last_AF_intersection = Intersects_AB;
716  CurrentE_Point = PointGe;
717 
718  G4bool validNormalAB;
719  NormalAtEntry = GetSurfaceNormal( PointGe, validNormalAB );
720  validNormalAtE = validNormalAB;
721  }
722 
723  depth--;
724  fin_section_depth[depth]=true;
725  }
726  } // if(!found_aproximate_intersection)
727 
728  } while ( ( ! found_approximate_intersection )
729  && ( ! there_is_no_intersection )
730  && ( substep_no <= max_substeps) ); // UNTIL found or failed
731 
732  if( substep_no > max_no_seen )
733  {
734  max_no_seen = substep_no;
735 #ifdef G4DEBUG_LOCATE_INTERSECTION
736  if( max_no_seen > warn_substeps )
737  {
738  trigger_substepno_print = max_no_seen-20; // Want to see last 20 steps
739  }
740 #endif
741  }
742 
743  if( ( substep_no >= max_substeps)
744  && !there_is_no_intersection
745  && !found_approximate_intersection )
746  {
747  G4cout << "ERROR - G4BrentLocator::EstimateIntersectionPoint()" << G4endl
748  << " Start and end-point of requested Step:" << G4endl;
749  printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
750  -1.0, NewSafety, 0);
751  G4cout << " Start and end-point of current Sub-Step:" << G4endl;
752  printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
753  -1.0, NewSafety, substep_no-1);
754  printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
755  -1.0, NewSafety, substep_no);
756  std::ostringstream message;
757  message << "Too many substeps!" << G4endl
758  << " Convergence is requiring too many substeps: "
759  << substep_no << G4endl
760  << " Abandoning effort to intersect. " << G4endl
761  << " Found intersection = "
762  << found_approximate_intersection << G4endl
763  << " Intersection exists = "
764  << !there_is_no_intersection << G4endl;
765  oldprc = G4cout.precision( 10 );
766  G4double done_len = CurrentA_PointVelocity.GetCurveLength();
767  G4double full_len = CurveEndPointVelocity.GetCurveLength();
768  message << " Undertaken only length: " << done_len
769  << " out of " << full_len << " required." << G4endl
770  << " Remaining length = " << full_len - done_len;
771  G4cout.precision( oldprc );
772 
773  G4Exception("G4BrentLocator::EstimateIntersectionPoint()",
774  "GeomNav0003", FatalException, message);
775  }
776  else if( substep_no >= warn_substeps )
777  {
778  oldprc= G4cout.precision( 10 );
779  std::ostringstream message;
780  message << "Many substeps while trying to locate intersection."
781  << G4endl
782  << " Undertaken length: "
783  << CurrentB_PointVelocity.GetCurveLength()
784  << " - Needed: " << substep_no << " substeps." << G4endl
785  << " Warning level = " << warn_substeps
786  << " and maximum substeps = " << max_substeps;
787  G4Exception("G4BrentLocator::EstimateIntersectionPoint()",
788  "GeomNav1002", JustWarning, message);
789  G4cout.precision( oldprc );
790  }
791  return !there_is_no_intersection; // Success or failure
792 }
G4FieldTrack ReEstimateEndpoint(const G4FieldTrack &CurrentStateA, const G4FieldTrack &EstimtdEndStateB, G4double linearDistSq, G4double curveDist)
void SetPosition(G4ThreeVector nPos)
G4double GetCurveLength() const
CLHEP::Hep3Vector G4ThreeVector
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)
G4FieldTrack * ptrInterMedFT[max_depth+1]
const G4ThreeVector & GetMomentumDir() const
static const G4double tolerance
#define G4ThreadLocal
Definition: tls.hh:89
int G4int
Definition: G4Types.hh:78
G4int maxNumberOfStepsForIntersection
void printStatus(const G4FieldTrack &startFT, const G4FieldTrack &currentFT, G4double requestStep, G4double safety, G4int stepNum)
G4int maxNumberOfCallsToReIntegration_depth
G4double GetEpsilonStepFor()
G4ThreeVector GetPosition() const
G4GLOB_DLL std::ostream G4cout
bool G4bool
Definition: G4Types.hh:79
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
G4bool GetAdjustementOfFoundIntersection()
#define G4endl
Definition: G4ios.hh:61
static const G4int max_depth
T sqr(const T &x)
Definition: templates.hh:145
double G4double
Definition: G4Types.hh:76
G4int maxNumberOfCallsToReIntegration
G4MagInt_Driver * GetIntegrationDriver()
static const double mm
Definition: G4SIunits.hh:102
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:579
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)