ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4SimpleLocator.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4SimpleLocator.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 G4SimpleLocator implementation
27 //
28 // 27.10.08 - Tatiana Nikitina, extracted from G4PropagatorInField class
29 // 04.10.11 - John Apostolakis, revised convergence to use Surface Normal
30 // ---------------------------------------------------------------------------
31 
32 #include <iomanip>
33 
34 #include "G4ios.hh"
35 #include "G4SimpleLocator.hh"
36 
37 G4SimpleLocator::G4SimpleLocator(G4Navigator *theNavigator)
38  : G4VIntersectionLocator(theNavigator)
39 {
40 }
41 
42 G4SimpleLocator::~G4SimpleLocator()
43 {
44 }
45 
46 // --------------------------------------------------------------------------
47 // G4bool G4PropagatorInField::LocateIntersectionPoint(
48 // const G4FieldTrack& CurveStartPointVelocity, // A
49 // const G4FieldTrack& CurveEndPointVelocity, // B
50 // const G4ThreeVector& TrialPoint, // E
51 // G4FieldTrack& IntersectedOrRecalculated // Output
52 // G4bool& recalculated ) // Out
53 // --------------------------------------------------------------------------
54 //
55 // Function that returns the intersection of the true path with the surface
56 // of the current volume (either the external one or the inner one with one
57 // of the daughters:
58 //
59 // A = Initial point
60 // B = another point
61 //
62 // Both A and B are assumed to be on the true path:
63 //
64 // E is the first point of intersection of the chord AB with
65 // a volume other than A (on the surface of A or of a daughter)
66 //
67 // Convention of Use :
68 // i) If it returns "true", then IntersectionPointVelocity is set
69 // to the approximate intersection point.
70 // ii) If it returns "false", no intersection was found.
71 // The validity of IntersectedOrRecalculated depends on 'recalculated'
72 // a) if latter is false, then IntersectedOrRecalculated is invalid.
73 // b) if latter is true, then IntersectedOrRecalculated is
74 // the new endpoint, due to a re-integration.
75 // --------------------------------------------------------------------------
76 // NOTE: implementation taken from G4PropagatorInField
77 //
78 G4bool G4SimpleLocator::EstimateIntersectionPoint(
79  const G4FieldTrack& CurveStartPointVelocity, // A
80  const G4FieldTrack& CurveEndPointVelocity, // B
81  const G4ThreeVector& TrialPoint, // E
82  G4FieldTrack& IntersectedOrRecalculatedFT, // Output
83  G4bool& recalculatedEndPoint, // Out
84  G4double &fPreviousSafety, //In/Out
85  G4ThreeVector &fPreviousSftOrigin) //In/Out
86 {
87  // Find Intersection Point ( A, B, E ) of true path AB - start at E.
88 
89  G4bool found_approximate_intersection = false;
90  G4bool there_is_no_intersection = false;
91 
92  G4FieldTrack CurrentA_PointVelocity = CurveStartPointVelocity;
93  G4FieldTrack CurrentB_PointVelocity = CurveEndPointVelocity;
94  G4ThreeVector CurrentE_Point = TrialPoint;
95  G4bool validNormalAtE = false;
96  G4ThreeVector NormalAtEntry;
97 
98  G4FieldTrack ApproxIntersecPointV(CurveEndPointVelocity); // FT-Def-Construct
99  G4double NewSafety = 0.0;
100  G4bool last_AF_intersection = false;
101  G4bool final_section = true; // Shows whether current section is last
102  // (i.e. B=full end)
103  recalculatedEndPoint = false;
104 
105  G4bool restoredFullEndpoint = false;
106 
107  G4int substep_no = 0;
108 
109  // Limits for substep number
110  //
111  const G4int max_substeps = 100000000; // Test 120 (old value 100 )
112  const G4int warn_substeps = 1000; // 100
113 
114  // Statistics for substeps
115  //
116  static G4ThreadLocal G4int max_no_seen= -1;
117 
118  NormalAtEntry = GetSurfaceNormal( CurrentE_Point, validNormalAtE);
119 
120 #ifdef G4DEBUG_FIELD
121  const G4double tolerance = 1.0e-8;
122  G4ThreeVector StartPosition= CurveStartPointVelocity.GetPosition();
123  if( (TrialPoint - StartPosition).mag() < tolerance * CLHEP::mm )
124  {
125  G4Exception("G4SimpleLocator::EstimateIntersectionPoint()",
126  "GeomNav1002", JustWarning,
127  "Intersection point F is exactly at start point A." );
128  }
129 #endif
130 
131  do
132  {
133  G4ThreeVector Point_A = CurrentA_PointVelocity.GetPosition();
134  G4ThreeVector Point_B = CurrentB_PointVelocity.GetPosition();
135 
136  // F = a point on true AB path close to point E
137  // (the closest if possible)
138  //
139  ApproxIntersecPointV = GetChordFinderFor()
140  ->ApproxCurvePointV( CurrentA_PointVelocity,
141  CurrentB_PointVelocity,
142  CurrentE_Point,
143  GetEpsilonStepFor());
144  // The above method is the key & most intuitive part ...
145 
146 #ifdef G4DEBUG_FIELD
147  if( ApproxIntersecPointV.GetCurveLength() >
148  CurrentB_PointVelocity.GetCurveLength() * (1.0 + tolerance) )
149  {
150  G4Exception("G4SimpleLocator::EstimateIntersectionPoint()",
151  "GeomNav0003", FatalException,
152  "Intermediate F point is past end B point!" );
153  }
154 #endif
155 
156  G4ThreeVector CurrentF_Point = ApproxIntersecPointV.GetPosition();
157 
158  // First check whether EF is small - then F is a good approx. point
159  // Calculate the length and direction of the chord AF
160  //
161  G4ThreeVector ChordEF_Vector = CurrentF_Point - CurrentE_Point;
162 
163  G4ThreeVector NewMomentumDir = ApproxIntersecPointV.GetMomentumDir();
164  G4double MomDir_dot_Norm = NewMomentumDir.dot( NormalAtEntry ) ;
165 
166  G4ThreeVector ChordAB = Point_B - Point_A;
167 
168 #ifdef G4DEBUG_FIELD
169  G4VIntersectionLocator::
170  ReportTrialStep( substep_no, ChordAB, ChordEF_Vector,
171  NewMomentumDir, NormalAtEntry, validNormalAtE );
172 #endif
173  // Check Sign is always exiting !! TODO
174  // Could ( > -epsilon) be used instead?
175  //
176  G4bool adequate_angle = ( MomDir_dot_Norm >= 0.0 )
177  || (! validNormalAtE ); // Invalid
178  G4double EF_dist2= ChordEF_Vector.mag2();
179  if ( ( EF_dist2 <= sqr(fiDeltaIntersection) && ( adequate_angle ) )
180  || ( EF_dist2 <= kCarTolerance*kCarTolerance ) )
181  {
182  found_approximate_intersection = true;
183 
184  // Create the "point" return value
185  //
186  IntersectedOrRecalculatedFT = ApproxIntersecPointV;
187  IntersectedOrRecalculatedFT.SetPosition( CurrentE_Point );
188 
189  if ( GetAdjustementOfFoundIntersection() )
190  {
191  // Try to Get Correction of IntersectionPoint using SurfaceNormal()
192  //
193  G4ThreeVector IP;
194  G4ThreeVector MomentumDir= ApproxIntersecPointV.GetMomentumDirection();
195  G4bool goodCorrection = AdjustmentOfFoundIntersection( Point_A,
196  CurrentE_Point, CurrentF_Point, MomentumDir,
197  last_AF_intersection, IP, NewSafety,
198  fPreviousSafety, fPreviousSftOrigin );
199 
200  if ( goodCorrection )
201  {
202  IntersectedOrRecalculatedFT = ApproxIntersecPointV;
203  IntersectedOrRecalculatedFT.SetPosition(IP);
204  }
205  }
206 
207  // Note: in order to return a point on the boundary,
208  // we must return E. But it is F on the curve.
209  // So we must "cheat": we are using the position at point E
210  // and the velocity at point F !!!
211  //
212  // This must limit the length we can allow for displacement!
213  }
214  else // E is NOT close enough to the curve (ie point F)
215  {
216  // Check whether any volumes are encountered by the chord AF
217  // ---------------------------------------------------------
218  // First relocate to restore any Voxel etc information
219  // in the Navigator before calling ComputeStep()
220  //
221  GetNavigatorFor()->LocateGlobalPointWithinVolume( Point_A );
222 
223  G4ThreeVector PointG; // Candidate intersection point
224  G4double stepLengthAF;
225  G4bool usedNavigatorAF = false;
226  G4bool Intersects_AF = IntersectChord( Point_A,
227  CurrentF_Point,
228  NewSafety,
229  fPreviousSafety,
230  fPreviousSftOrigin,
231  stepLengthAF,
232  PointG,
233  &usedNavigatorAF );
234  last_AF_intersection = Intersects_AF;
235  if( Intersects_AF )
236  {
237  // G is our new Candidate for the intersection point.
238  // It replaces "E" and we will repeat the test to see if
239  // it is a good enough approximate point for us.
240  // B <- F
241  // E <- G
242 
243  CurrentB_PointVelocity = ApproxIntersecPointV;
244  CurrentE_Point = PointG;
245 
246  // Need to recalculate the Exit Normal at the new PointG
247  // Relies on a call to Navigator::ComputeStep in IntersectChord above
248  // If the safety was adequate (for the step) this would NOT be called!
249  // But this must not occur, no intersection can be found in that case,
250  // so this branch, ie if( Intersects_AF ) would not be reached!
251  //
252  G4bool validNormalLast;
253  NormalAtEntry = GetSurfaceNormal( PointG, validNormalLast );
254  validNormalAtE = validNormalLast;
255 
256  // By moving point B, must take care if current
257  // AF has no intersection to try current FB!!
258  //
259  final_section = false;
260 
261 #ifdef G4VERBOSE
262  if( fVerboseLevel > 3 )
263  {
264  G4cout << "G4PiF::LI> Investigating intermediate point"
265  << " at s=" << ApproxIntersecPointV.GetCurveLength()
266  << " on way to full s="
267  << CurveEndPointVelocity.GetCurveLength() << G4endl;
268  }
269 #endif
270  }
271  else // not Intersects_AF
272  {
273  // In this case:
274  // There is NO intersection of AF with a volume boundary.
275  // We must continue the search in the segment FB!
276  //
277  GetNavigatorFor()->LocateGlobalPointWithinVolume( CurrentF_Point );
278 
279  G4double stepLengthFB;
280  G4ThreeVector PointH;
281  G4bool usedNavigatorFB = false;
282 
283  // Check whether any volumes are encountered by the chord FB
284  // ---------------------------------------------------------
285 
286  G4bool Intersects_FB = IntersectChord( CurrentF_Point, Point_B,
287  NewSafety,fPreviousSafety,
288  fPreviousSftOrigin,
289  stepLengthFB,
290  PointH, &usedNavigatorFB );
291  if( Intersects_FB )
292  {
293  // There is an intersection of FB with a volume boundary
294  // H <- First Intersection of Chord FB
295 
296  // H is our new Candidate for the intersection point.
297  // It replaces "E" and we will repeat the test to see if
298  // it is a good enough approximate point for us.
299 
300  // Note that F must be in volume volA (the same as A)
301  // (otherwise AF would meet a volume boundary!)
302  // A <- F
303  // E <- H
304  //
305  CurrentA_PointVelocity = ApproxIntersecPointV;
306  CurrentE_Point = PointH;
307 
308  // Need to recalculate the Exit Normal at the new PointG
309  // Relies on call to Navigator::ComputeStep in IntersectChord above
310  // If safety was adequate (for the step) this would NOT be called!
311  // But this must not occur, no intersection found in that case,
312  // so this branch, i.e. if( Intersects_AF ) would not be reached!
313  //
314  G4bool validNormalLast;
315  NormalAtEntry = GetSurfaceNormal( PointH, validNormalLast );
316  validNormalAtE = validNormalLast;
317  }
318  else // not Intersects_FB
319  {
320  // There is NO intersection of FB with a volume boundary
321 
322  if( final_section )
323  {
324  // If B is the original endpoint, this means that whatever
325  // volume(s) intersected the original chord, none touch the
326  // smaller chords we have used.
327  // The value of 'IntersectedOrRecalculatedFT' returned is
328  // likely not valid
329 
330  there_is_no_intersection = true; // real final_section
331  }
332  else
333  {
334  // We must restore the original endpoint
335 
336  CurrentA_PointVelocity = CurrentB_PointVelocity; // Got to B
337  CurrentB_PointVelocity = CurveEndPointVelocity;
338  restoredFullEndpoint = true;
339  }
340  } // Endif (Intersects_FB)
341  } // Endif (Intersects_AF)
342 
343  // Ensure that the new endpoints are not further apart in space
344  // than on the curve due to different errors in the integration
345  //
346  G4double linDistSq, curveDist;
347  linDistSq = ( CurrentB_PointVelocity.GetPosition()
348  - CurrentA_PointVelocity.GetPosition() ).mag2();
349  curveDist = CurrentB_PointVelocity.GetCurveLength()
350  - CurrentA_PointVelocity.GetCurveLength();
351 
352  // Change this condition for very strict parameters of propagation
353  //
354  if( curveDist*curveDist*(1+2* GetEpsilonStepFor()) < linDistSq )
355  {
356  // Re-integrate to obtain a new B
357  //
358  G4FieldTrack newEndPointFT =
359  ReEstimateEndpoint( CurrentA_PointVelocity,
360  CurrentB_PointVelocity,
361  linDistSq, // to avoid recalculation
362  curveDist );
363  G4FieldTrack oldPointVelB = CurrentB_PointVelocity;
364  CurrentB_PointVelocity = newEndPointFT;
365 
366  if( (final_section)) // real final section
367  {
368  recalculatedEndPoint = true;
369  IntersectedOrRecalculatedFT = newEndPointFT;
370  // So that we can return it, if it is the endpoint!
371  }
372  }
373  if( curveDist < 0.0 )
374  {
375  fVerboseLevel = 5; // Print out a maximum of information
376  printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
377  -1.0, NewSafety, substep_no );
378  std::ostringstream message;
379  message << "Error in advancing propagation." << G4endl
380  << " Point A (start) is " << CurrentA_PointVelocity
381  << G4endl
382  << " Point B (end) is " << CurrentB_PointVelocity
383  << G4endl
384  << " Curve distance is " << curveDist << G4endl
385  << G4endl
386  << "The final curve point is not further along"
387  << " than the original!" << G4endl;
388 
389  if( recalculatedEndPoint )
390  {
391  message << "Recalculation of EndPoint was called with fEpsStep= "
392  << GetEpsilonStepFor() << G4endl;
393  }
394  message.precision(20);
395  message << " Point A (Curve start) is " << CurveStartPointVelocity
396  << G4endl
397  << " Point B (Curve end) is " << CurveEndPointVelocity
398  << G4endl
399  << " Point A (Current start) is " << CurrentA_PointVelocity
400  << G4endl
401  << " Point B (Current end) is " << CurrentB_PointVelocity
402  << G4endl
403  << " Point E (Trial Point) is " << CurrentE_Point
404  << G4endl
405  << " Point F (Intersection) is " << ApproxIntersecPointV
406  << G4endl
407  << " LocateIntersection parameters are : Substep no= "
408  << substep_no;
409 
410  G4Exception("G4SimpleLocator::EstimateIntersectionPoint()",
411  "GeomNav0003", FatalException, message);
412  }
413 
414  if ( restoredFullEndpoint )
415  {
416  final_section = restoredFullEndpoint;
417  restoredFullEndpoint = false;
418  }
419  } // EndIf ( E is close enough to the curve, ie point F. )
420  // tests ChordAF_Vector.mag() <= maximum_lateral_displacement
421 
422 #ifdef G4DEBUG_LOCATE_INTERSECTION
423  G4int trigger_substepno_print= warn_substeps - 20;
424 
425  if( substep_no >= trigger_substepno_print )
426  {
427  G4cout << "Difficulty in converging in "
428  << "G4SimpleLocator::EstimateIntersectionPoint():"
429  << G4endl
430  << " Substep no = " << substep_no << G4endl;
431  if( substep_no == trigger_substepno_print )
432  {
433  printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
434  -1.0, NewSafety, 0);
435  }
436  G4cout << " State of point A: ";
437  printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
438  -1.0, NewSafety, substep_no-1, 0);
439  G4cout << " State of point B: ";
440  printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
441  -1.0, NewSafety, substep_no);
442  }
443 #endif
444  ++substep_no;
445 
446  } while ( ( ! found_approximate_intersection )
447  && ( ! there_is_no_intersection )
448  && ( substep_no <= max_substeps) ); // UNTIL found or failed
449 
450  if( substep_no > max_no_seen )
451  {
452  max_no_seen = substep_no;
453 #ifdef G4DEBUG_LOCATE_INTERSECTION
454  if( max_no_seen > warn_substeps )
455  {
456  trigger_substepno_print = max_no_seen-20; // Want to see last 20 steps
457  }
458 #endif
459  }
460 
461  if( ( substep_no >= max_substeps)
462  && !there_is_no_intersection
463  && !found_approximate_intersection )
464  {
465  G4cout << "ERROR - G4SimpleLocator::EstimateIntersectionPoint()" << G4endl
466  << " Start and Endpoint of Requested Step:" << G4endl;
467  printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
468  -1.0, NewSafety, 0);
469  G4cout << G4endl
470  << " Start and end-point of current Sub-Step:" << G4endl;
471  printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
472  -1.0, NewSafety, substep_no-1);
473  printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
474  -1.0, NewSafety, substep_no);
475 
476  std::ostringstream message;
477  message << "Convergence is requiring too many substeps: "
478  << substep_no << G4endl
479  << " Abandoning effort to intersect." << G4endl
480  << " Found intersection = "
481  << found_approximate_intersection << G4endl
482  << " Intersection exists = "
483  << !there_is_no_intersection << G4endl;
484  message.precision(10);
485  G4double done_len = CurrentA_PointVelocity.GetCurveLength();
486  G4double full_len = CurveEndPointVelocity.GetCurveLength();
487  message << " Undertaken only length: " << done_len
488  << " out of " << full_len << " required." << G4endl
489  << " Remaining length = " << full_len-done_len;
490 
491  G4Exception("G4SimpleLocator::EstimateIntersectionPoint()",
492  "GeomNav0003", FatalException, message);
493  }
494  else if( substep_no >= warn_substeps )
495  {
496  std::ostringstream message;
497  message.precision(10);
498 
499  message << "Many substeps while trying to locate intersection." << G4endl
500  << " Undertaken length: "
501  << CurrentB_PointVelocity.GetCurveLength()
502  << " - Needed: " << substep_no << " substeps." << G4endl
503  << " Warning level = " << warn_substeps
504  << " and maximum substeps = " << max_substeps;
505  G4Exception("G4SimpleLocator::EstimateIntersectionPoint()",
506  "GeomNav1002", JustWarning, message);
507  }
508  return !there_is_no_intersection; // Success or failure
509 }