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