Geant4  10.03.p01
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Public Member Functions | List of all members
G4MultiLevelLocator Class Reference

#include <G4MultiLevelLocator.hh>

Inheritance diagram for G4MultiLevelLocator:
Inheritance graph
[legend]
Collaboration diagram for G4MultiLevelLocator:
Collaboration graph
[legend]

Public Member Functions

 G4MultiLevelLocator (G4Navigator *theNavigator)
 
 ~G4MultiLevelLocator ()
 
G4bool EstimateIntersectionPoint (const G4FieldTrack &curveStartPointTangent, const G4FieldTrack &curveEndPointTangent, const G4ThreeVector &trialPoint, G4FieldTrack &intersectPointTangent, G4bool &recalculatedEndPoint, G4double &fPreviousSafety, G4ThreeVector &fPreviousSftOrigin)
 
void ReportStatistics ()
 
void SetMaxSteps (unsigned int valMax)
 
void SetWarnSteps (unsigned int valWarn)
 
- Public Member Functions inherited from G4VIntersectionLocator
 G4VIntersectionLocator (G4Navigator *theNavigator)
 
virtual ~G4VIntersectionLocator ()
 
void printStatus (const G4FieldTrack &startFT, const G4FieldTrack &currentFT, G4double requestStep, G4double safety, G4int stepNum)
 
G4bool IntersectChord (const G4ThreeVector &StartPointA, const G4ThreeVector &EndPointB, G4double &NewSafety, G4double &PreviousSafety, G4ThreeVector &PreviousSftOrigin, G4double &LinearStepLength, G4ThreeVector &IntersectionPoint, G4bool *calledNavigator=0)
 
void SetEpsilonStepFor (G4double EpsilonStep)
 
void SetDeltaIntersectionFor (G4double deltaIntersection)
 
void SetNavigatorFor (G4Navigator *fNavigator)
 
void SetChordFinderFor (G4ChordFinder *fCFinder)
 
void SetVerboseFor (G4int fVerbose)
 
G4int GetVerboseFor ()
 
G4double GetDeltaIntersectionFor ()
 
G4double GetEpsilonStepFor ()
 
G4Navigator * GetNavigatorFor ()
 
G4ChordFinder * GetChordFinderFor ()
 
void SetSafetyParametersFor (G4bool UseSafety)
 
void AddAdjustementOfFoundIntersection (G4bool UseCorrection)
 
G4bool GetAdjustementOfFoundIntersection ()
 
void AdjustIntersections (G4bool UseCorrection)
 
G4bool AreIntersectionsAdjusted ()
 

Additional Inherited Members

- Static Public Member Functions inherited from G4VIntersectionLocator
static void printStatus (const G4FieldTrack &startFT, const G4FieldTrack &currentFT, G4double requestStep, G4double safety, G4int stepNum, std::ostream &oss, G4int verboseLevel)
 
- Protected Member Functions inherited from G4VIntersectionLocator
G4FieldTrack ReEstimateEndpoint (const G4FieldTrack &CurrentStateA, const G4FieldTrack &EstimtdEndStateB, G4double linearDistSq, G4double curveDist)
 
G4bool CheckAndReEstimateEndpoint (const G4FieldTrack &CurrentStartA, const G4FieldTrack &EstimatedEndB, G4FieldTrack &RevisedEndPoint, G4int &errorCode)
 
G4ThreeVector GetSurfaceNormal (const G4ThreeVector &CurrentInt_Point, G4bool &validNormal)
 
G4ThreeVector GetGlobalSurfaceNormal (const G4ThreeVector &CurrentE_Point, G4bool &validNormal)
 
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)
 
void ReportTrialStep (G4int step_no, const G4ThreeVector &ChordAB_v, const G4ThreeVector &ChordEF_v, const G4ThreeVector &NewMomentumDir, const G4ThreeVector &NormalAtEntry, G4bool validNormal)
 
G4bool LocateGlobalPointWithinVolumeAndCheck (const G4ThreeVector &pos)
 
void LocateGlobalPointWithinVolumeCheckAndReport (const G4ThreeVector &pos, const G4String &CodeLocationInfo, G4int CheckMode)
 
void SetCheckMode (G4bool value)
 
G4bool GetCheckMode ()
 
void ReportReversedPoints (std::ostringstream &ossMsg, const G4FieldTrack &StartPointVel, const G4FieldTrack &EndPointVel, G4double NewSafety, G4double epsStep, const G4FieldTrack &CurrentA_PointVelocity, const G4FieldTrack &CurrentB_PointVelocity, const G4FieldTrack &SubStart_PointVelocity, const G4ThreeVector &CurrentE_Point, const G4FieldTrack &ApproxIntersecPointV, G4int sbstp_no, G4int sbstp_no_p, G4int depth)
 
void ReportProgress (std::ostream &oss, const G4FieldTrack &StartPointVel, const G4FieldTrack &EndPointVel, G4int substep_no, const G4FieldTrack &A_PtVel, const G4FieldTrack &B_PtVel, G4double safetyLast, G4int depth=-1)
 
void ReportImmediateHit (const char *MethodName, const G4ThreeVector &StartPosition, const G4ThreeVector &TrialPoint, double tolerance, unsigned long int numCalls)
 
- Protected Attributes inherited from G4VIntersectionLocator
G4double kCarTolerance
 
G4int fVerboseLevel
 
G4bool fUseNormalCorrection
 
G4bool fCheckMode
 
G4Navigator * fiNavigator
 
G4ChordFinder * fiChordFinder
 
G4double fiEpsilonStep
 
G4double fiDeltaIntersection
 
G4bool fiUseSafety
 
G4Navigator * fHelpingNavigator
 
G4TouchableHistory * fpTouchable
 

Detailed Description

Definition at line 51 of file G4MultiLevelLocator.hh.

Constructor & Destructor Documentation

G4MultiLevelLocator::G4MultiLevelLocator ( G4Navigator *  theNavigator)

Definition at line 40 of file G4MultiLevelLocator.cc.

41  : G4VIntersectionLocator(theNavigator),
42  fMaxSteps(10000), // Very loose - allows many steps (looping will be rare)
43  fWarnSteps(1000), //
44  fNumCalls(0),
45  fNumAdvanceFull(0.), fNumAdvanceGood(0), fNumAdvanceTrials(0)
46 {
47  // In case of too slow progress in finding Intersection Point
48  // intermediates Points on the Track must be stored.
49  // Initialise the array of Pointers [max_depth+1] to do this
50 
51  G4ThreeVector zeroV(0.0,0.0,0.0);
52  for (G4int idepth=0; idepth<max_depth+1; idepth++ )
53  {
54  ptrInterMedFT[ idepth ] = new G4FieldTrack( zeroV, zeroV, 0., 0., 0., 0.);
55  }
56 
57 #ifdef G4DEBUG_FIELD
58  // Trial values Loose Tight
59  // To happen: Infrequent Often
60  SetMaxSteps(50); // 300 25
61  SetWarnSteps(40); // 250 15
62 #endif
63 }
G4MultiLevelLocator::SetWarnSteps
void SetWarnSteps(unsigned int valWarn)
Definition: G4MultiLevelLocator.hh:76
G4int
int G4int
Definition: G4Types.hh:78
G4VIntersectionLocator::G4VIntersectionLocator
G4VIntersectionLocator(G4Navigator *theNavigator)
Definition: G4VIntersectionLocator.cc:47
G4MultiLevelLocator::SetMaxSteps
void SetMaxSteps(unsigned int valMax)
Definition: G4MultiLevelLocator.hh:75

Here is the call graph for this function:

G4MultiLevelLocator::~G4MultiLevelLocator ( )

Definition at line 65 of file G4MultiLevelLocator.cc.

66 {
67  for ( G4int idepth=0; idepth<max_depth+1; idepth++)
68  {
69  delete ptrInterMedFT[idepth];
70  }
71 #ifdef G4DEBUG_FIELD
72  ReportStatistics();
73 #endif
74 }
G4int
int G4int
Definition: G4Types.hh:78

Here is the call graph for this function:

Member Function Documentation

G4bool G4MultiLevelLocator::EstimateIntersectionPoint ( const G4FieldTrack &  curveStartPointTangent,
const G4FieldTrack &  curveEndPointTangent,
const G4ThreeVector &  trialPoint,
G4FieldTrack &  intersectPointTangent,
G4bool &  recalculatedEndPoint,
G4double &  fPreviousSafety,
G4ThreeVector &  fPreviousSftOrigin 
)
virtual

Implements G4VIntersectionLocator.

Definition at line 119 of file G4MultiLevelLocator.cc.

127 {
128  // Find Intersection Point ( A, B, E ) of true path AB - start at E.
129  const char* MethodName= "G4MultiLevelLocator::EstimateIntersectionPoint()";
130 
131  G4bool found_approximate_intersection = false;
132  G4bool there_is_no_intersection = false;
133 
134  G4FieldTrack CurrentA_PointVelocity = CurveStartPointVelocity;
135  G4FieldTrack CurrentB_PointVelocity = CurveEndPointVelocity;
136  G4ThreeVector CurrentE_Point = TrialPoint;
137  G4bool validNormalAtE = false;
138  G4ThreeVector NormalAtEntry;
139 
140  G4FieldTrack ApproxIntersecPointV(CurveEndPointVelocity); // FT-Def-Construct
141  // G4bool validApproxIntPV= false; // Is it current: valid and up-to-date?
142  G4bool validIntersectP= true; // Is it current ?
143  G4double NewSafety = 0.0;
144  G4bool last_AF_intersection = false;
145 
146  // G4bool final_section= true; // Shows whether current section is last
147  // (i.e. B=full end)
148  G4bool first_section = true;
149 
150  recalculatedEndPoint = false;
151 
152  G4bool restoredFullEndpoint = false;
153 
154  unsigned int substep_no = 0;
155 
156  // Statistics for substeps
157  //
158  static G4ThreadLocal unsigned int max_no_seen= 0;
159 
160  //--------------------------------------------------------------------------
161  // Algorithm for the case if progress in founding intersection is too slow.
162  // Process is defined too slow if after N=param_substeps advances on the
163  // path, it will be only 'fraction_done' of the total length.
164  // In this case the remaining length is divided in two half and
165  // the loop is restarted for each half.
166  // If progress is still too slow, the division in two halfs continue
167  // until 'max_depth'.
168  //--------------------------------------------------------------------------
169 
170  const G4int param_substeps=5; // Test value for the maximum number
171  // of substeps
172  const G4double fraction_done=0.3;
173 
174  G4bool Second_half = false; // First half or second half of divided step
175 
176  // We need to know this for the 'final_section':
177  // real 'final_section' or first half 'final_section'
178  // In algorithm it is considered that the 'Second_half' is true
179  // and it becomes false only if we are in the first-half of level
180  // depthness or if we are in the first section
181 
182  unsigned int depth=0; // Depth counts subdivisions of initial step made
183  fNumCalls++;
184 
185 #ifdef G4DEBUG_FIELD
186  unsigned int trigger_substepno_print=0;
187  const G4double tolerance = 1.0e-8 * CLHEP::mm;
188  unsigned int biggest_depth= 0;
189 
190 #if (G4DEBUG_FIELD>1)
191  G4ThreeVector StartPosition= CurveStartPointVelocity.GetPosition();
192  if( (TrialPoint - StartPosition).mag2() < tolerance*tolerance)
193  {
194  ReportImmediateHit( MethodName, StartPosition, TrialPoint,
195  tolerance, fNumCalls);
196  }
197 #endif
198 #endif
199 
200  NormalAtEntry = GetSurfaceNormal(CurrentE_Point, validNormalAtE);
201 
202  // Intermediates Points on the Track = Subdivided Points must be stored.
203  // Give the initial values to 'InterMedFt'
204  // Important is 'ptrInterMedFT[0]', it saves the 'EndCurvePoint'
205  //
206  *ptrInterMedFT[0] = CurveEndPointVelocity;
207  for (G4int idepth=1; idepth<max_depth+1; idepth++ )
208  {
209  *ptrInterMedFT[idepth]=CurveStartPointVelocity;
210  }
211 
212  // Final_section boolean store
213  //
214  G4bool fin_section_depth[max_depth];
215  for (G4int idepth=0; idepth<max_depth; idepth++ )
216  {
217  fin_section_depth[idepth]=true;
218  }
219  // 'SubStartPoint' is needed to calculate the length of the divided step
220  //
221  G4FieldTrack SubStart_PointVelocity = CurveStartPointVelocity;
222 
223  do // Loop checking, 07.10.2016, J.Apostolakis
224  {
225  unsigned int substep_no_p = 0;
226  G4bool sub_final_section = false; // the same as final_section,
227  // but for 'sub_section'
228  SubStart_PointVelocity = CurrentA_PointVelocity;
229 
230  do // Loop checking, 07.10.2016, J.Apostolakis
231  { // REPEAT param
232  G4ThreeVector Point_A = CurrentA_PointVelocity.GetPosition();
233  G4ThreeVector Point_B = CurrentB_PointVelocity.GetPosition();
234 
235  // F = a point on true AB path close to point E
236  // (the closest if possible)
237  //
238  ApproxIntersecPointV = GetChordFinderFor()
239  ->ApproxCurvePointV( CurrentA_PointVelocity,
240  CurrentB_PointVelocity,
241  CurrentE_Point,
242  GetEpsilonStepFor());
243  // The above method is the key & most intuitive part ...
244 
245  // validApproxIntPV = true;
246 
247 #ifdef G4DEBUG_FIELD
248  if( ApproxIntersecPointV.GetCurveLength() >
249  CurrentB_PointVelocity.GetCurveLength() * (1.0 + tolerance) )
250  {
251  G4Exception(MethodName, "GeomNav0003", FatalException,
252  "Intermediate F point is past end B point" );
253  }
254 #endif
255 
256  G4ThreeVector CurrentF_Point= ApproxIntersecPointV.GetPosition();
257 
258  // First check whether EF is small - then F is a good approx. point
259  // Calculate the length and direction of the chord AF
260  //
261  G4ThreeVector ChordEF_Vector = CurrentF_Point - CurrentE_Point;
262 
263  G4ThreeVector NewMomentumDir= ApproxIntersecPointV.GetMomentumDir();
264  G4double MomDir_dot_Norm= NewMomentumDir.dot( NormalAtEntry ) ;
265 
266 #ifdef G4DEBUG_FIELD
267  if( fVerboseLevel > 3 )
268  {
269  G4ThreeVector ChordAB = Point_B - Point_A;
270  G4double ABchord_length = ChordAB.mag();
271  G4double MomDir_dot_ABchord;
272  MomDir_dot_ABchord = (1.0 / ABchord_length)
273  * NewMomentumDir.dot( ChordAB );
274  G4VIntersectionLocator::ReportTrialStep( substep_no, ChordAB,
275  ChordEF_Vector, NewMomentumDir, NormalAtEntry, validNormalAtE );
276  G4cout << " dot( MomentumDir, ABchord_unit ) = "
277  << MomDir_dot_ABchord << G4endl;
278  }
279 #endif
280  G4bool adequate_angle =
281  ( MomDir_dot_Norm >= 0.0 ) // Can use ( > -epsilon) instead
282  || (! validNormalAtE ); // Invalid, cannot use this criterion
283  G4double EF_dist2 = ChordEF_Vector.mag2();
284  if ( ( EF_dist2 <= sqr(fiDeltaIntersection) && ( adequate_angle ) )
285  || ( EF_dist2 <= kCarTolerance*kCarTolerance ) )
286  {
287  found_approximate_intersection = true;
288 
289  // Create the "point" return value
290  //
291  IntersectedOrRecalculatedFT = ApproxIntersecPointV;
292  IntersectedOrRecalculatedFT.SetPosition( CurrentE_Point );
293 
294  if ( GetAdjustementOfFoundIntersection() )
295  {
296  // Try to Get Correction of IntersectionPoint using SurfaceNormal()
297  //
298  G4ThreeVector IP;
299  G4ThreeVector MomentumDir=ApproxIntersecPointV.GetMomentumDirection();
300  G4bool goodCorrection = AdjustmentOfFoundIntersection(Point_A,
301  CurrentE_Point, CurrentF_Point, MomentumDir,
302  last_AF_intersection, IP, NewSafety,
303  previousSafety, previousSftOrigin );
304  if ( goodCorrection )
305  {
306  IntersectedOrRecalculatedFT = ApproxIntersecPointV;
307  IntersectedOrRecalculatedFT.SetPosition(IP);
308  }
309  }
310  // Note: in order to return a point on the boundary,
311  // we must return E. But it is F on the curve.
312  // So we must "cheat": we are using the position at point E
313  // and the velocity at point F !!!
314  //
315  // This must limit the length we can allow for displacement!
316  }
317  else // E is NOT close enough to the curve (ie point F)
318  {
319  // Check whether any volumes are encountered by the chord AF
320  // ---------------------------------------------------------
321  // First relocate to restore any Voxel etc information
322  // in the Navigator before calling ComputeStep()
323  //
324  GetNavigatorFor()->LocateGlobalPointWithinVolume( Point_A );
325 
326  G4ThreeVector PointG; // Candidate intersection point
327  G4double stepLengthAF;
328  G4bool Intersects_AF = IntersectChord( Point_A, CurrentF_Point,
329  NewSafety, previousSafety,
330  previousSftOrigin,
331  stepLengthAF,
332  PointG );
333  last_AF_intersection = Intersects_AF;
334  if( Intersects_AF )
335  {
336  // G is our new Candidate for the intersection point.
337  // It replaces "E" and we will repeat the test to see if
338  // it is a good enough approximate point for us.
339  // B <- F
340  // E <- G
341  //
342  CurrentB_PointVelocity = ApproxIntersecPointV;
343  CurrentE_Point = PointG;
344 
345  validIntersectP= true; // 'E' has been updated.
346  // validApproxIntPV= false; // 'F' is no longer valid, as B changed
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 
357 #ifdef G4VERBOSE
358  if( fVerboseLevel > 3 )
359  {
360  G4cout << "G4PiF::LI> Investigating intermediate point"
361  << " at s=" << ApproxIntersecPointV.GetCurveLength()
362  << " on way to full s="
363  << CurveEndPointVelocity.GetCurveLength() << G4endl;
364  }
365 #endif
366  }
367  else // not Intersects_AF
368  {
369  // In this case:
370  // There is NO intersection of AF with a volume boundary.
371  // We must continue the search in the segment FB!
372  //
373  GetNavigatorFor()->LocateGlobalPointWithinVolume( CurrentF_Point );
374 
375  G4double stepLengthFB;
376  G4ThreeVector PointH;
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, previousSafety,
383  previousSftOrigin,
384  stepLengthFB,
385  PointH );
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  CurrentA_PointVelocity = ApproxIntersecPointV;
401  CurrentE_Point = PointH;
402 
403  validIntersectP = true; // 'E' has been updated.
404  // validApproxIntPV = false; // 'F' is no longer valid, as A changed
405 
406  G4bool validNormalH;
407  NormalAtEntry = GetSurfaceNormal( PointH, validNormalH );
408  validNormalAtE = validNormalH;
409  }
410  else // not Intersects_FB
411  {
412  if(fin_section_depth[depth])
413  {
414  // If B is the original endpoint, this means that whatever
415  // volume(s) intersected the original chord, none touch the
416  // smaller chords we have used.
417  // The value of 'IntersectedOrRecalculatedFT' returned is
418  // likely not valid
419 
420  // Check on real final_section or SubEndSection
421  //
422  if( ((Second_half)&&(depth==0)) || (first_section) )
423  {
424  there_is_no_intersection = true; // real final_section
425  }
426  else
427  {
428  // end of subsection, not real final section
429  // exit from the and go to the depth-1 level
430  substep_no_p = param_substeps+2; // exit from the loop
431 
432  // but 'Second_half' is still true because we need to find
433  // the 'CurrentE_point' for the next loop
434  Second_half = true;
435  sub_final_section = true;
436  }
437  }
438  else
439  {
440  CurrentA_PointVelocity = CurrentB_PointVelocity; // Got to B
441  CurrentB_PointVelocity = (depth==0) ? CurveEndPointVelocity
442  : *ptrInterMedFT[depth] ;
443  SubStart_PointVelocity = CurrentA_PointVelocity;
444  restoredFullEndpoint = true;
445 
446  validIntersectP= false; // 'E' has NOT been updated.
447  // validApproxIntPV= false; // 'F' is no longer valid, A changed
448  }
449  } // Endif (Intersects_FB)
450  } // Endif (Intersects_AF)
451 
452  G4FieldTrack RevisedB_FT= CurrentB_PointVelocity;
453  G4int errorEndPt;
454  G4bool recalculatedB= CheckAndReEstimateEndpoint(CurrentA_PointVelocity,
455  CurrentB_PointVelocity,
456  RevisedB_FT,
457  errorEndPt );
458  if( recalculatedB )
459  {
460  CurrentB_PointVelocity= RevisedB_FT; // Use it !
461  // Do not invalidate intersection F -- it is still roughly OK.
462  //
463  // The best course would be to invalidate (reset validIntersectP)
464  // BUT if we invalidate it, we must re-estimate it somewhere!
465 
466  // validApproxIntPV= false; // 'F' is no longer valid, as B changed
467  // validIntersectP= false; // 'E' has NOT been updated.
468 
469  if ( (fin_section_depth[depth]) // real final section
470  &&( first_section || ((Second_half)&&(depth==0)) ) )
471  {
472  recalculatedEndPoint = true;
473  IntersectedOrRecalculatedFT = RevisedB_FT;
474  // So that we can return it, if it is the endpoint!
475  }
476  // else
477  // Move forward the other points
478  // - or better flag it, so that they are re-computed when next used
479  // [ Implementation: a counter for # of recomputations
480  // => avoids extra work]
481 
482  }
483  if( errorEndPt > 1 ) // errorEndPt = 1 is milder, just: len(B)=len(A)
484  {
485  std::ostringstream errmsg;
486  errmsg << "Location: " << MethodName
487  << "- After EndIf(Intersects_AF)" << G4endl;
488  ReportReversedPoints(errmsg,
489  CurveStartPointVelocity, CurveEndPointVelocity,
490  NewSafety, fiEpsilonStep,
491  CurrentA_PointVelocity, CurrentB_PointVelocity,
492  SubStart_PointVelocity, CurrentE_Point,
493  ApproxIntersecPointV, substep_no, substep_no_p, depth);
494  G4Exception(MethodName, "GeomNav0003", FatalException, errmsg);
495  }
496  if( restoredFullEndpoint )
497  {
498  fin_section_depth[depth] = restoredFullEndpoint;
499  restoredFullEndpoint = false;
500  }
501  } // EndIf ( E is close enough to the curve, ie point F. )
502  // tests ChordAF_Vector.mag() <= maximum_lateral_displacement
503 
504 #ifdef G4DEBUG_FIELD
505  if( trigger_substepno_print == 0)
506  {
507  trigger_substepno_print= fWarnSteps - 20;
508  }
509 
510  if( substep_no >= trigger_substepno_print )
511  {
512  G4cout << "Difficulty in converging in " << MethodName
513  << G4endl
514  << " Substep no = " << substep_no << G4endl;
515  if( substep_no == trigger_substepno_print )
516  {
517  printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
518  -1.0, NewSafety, 0);
519  }
520  G4cout << " State of point A: ";
521  printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
522  -1.0, NewSafety, substep_no-1);
523  G4cout << " State of point B: ";
524  printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
525  -1.0, NewSafety, substep_no);
526  }
527 #endif
528  substep_no++;
529  substep_no_p++;
530 
531  } while ( ( ! found_approximate_intersection )
532  && ( ! there_is_no_intersection )
533  && ( substep_no_p <= param_substeps) ); // UNTIL found or
534  // failed param substep
535 
536  if( (!found_approximate_intersection) && (!there_is_no_intersection) )
537  {
538  G4double did_len = std::abs( CurrentA_PointVelocity.GetCurveLength()
539  - SubStart_PointVelocity.GetCurveLength());
540  G4double all_len = std::abs( CurrentB_PointVelocity.GetCurveLength()
541  - SubStart_PointVelocity.GetCurveLength());
542 
543  G4double distAB= -1;
544  G4ThreeVector PointGe;
545  //
546  // Is progress is too slow, and is it possible to go deeper?
547  // If so, then *halve the step*
548  // ==============
549  if( (did_len < fraction_done*all_len)
550  && (depth<max_depth) && (!sub_final_section) )
551  {
552 #ifdef G4DEBUG_FIELD
553  static G4ThreadLocal unsigned int numSplits=0; // For debugging only
554  biggest_depth= std::max(depth, biggest_depth);
555  numSplits++;
556 #endif
557  Second_half=false;
558  depth++;
559  first_section = false;
560 
561  G4double Sub_len = (all_len-did_len)/(2.);
562  G4FieldTrack midPoint = CurrentA_PointVelocity;
563  G4MagInt_Driver* integrDriver
564  = GetChordFinderFor()->GetIntegrationDriver();
565  G4bool fullAdvance=
566  integrDriver->AccurateAdvance(midPoint, Sub_len, fiEpsilonStep);
567 
568  fNumAdvanceTrials++;
569  if( fullAdvance ) { fNumAdvanceFull++; }
570 
571  G4double lenAchieved=
572  midPoint.GetCurveLength()-CurrentA_PointVelocity.GetCurveLength();
573 
574  const G4double adequateFraction = (1.0-CLHEP::perThousand);
575  G4bool goodAdvance = (lenAchieved >= adequateFraction * Sub_len);
576  if ( goodAdvance ) { fNumAdvanceGood++; }
577 
578 #ifdef G4DEBUG_FIELD
579  else // !goodAdvance
580  {
581  G4cout << "MLL> AdvanceChordLimited not full at depth=" << depth
582  << " total length achieved = " << lenAchieved << " of "
583  << Sub_len << " fraction= ";
584  if (Sub_len != 0.0 ) { G4cout << lenAchieved / Sub_len; }
585  else { G4cout << "DivByZero"; }
586  G4cout << " Good-enough fraction = " << adequateFraction;
587  G4cout << " Number of call to mll = " << fNumCalls
588  << " iteration " << substep_no
589  << " inner = " << substep_no_p << G4endl;
590  G4cout << " Epsilon of integration = " << fiEpsilonStep << G4endl;
591  G4cout << " State at start is = " << CurrentA_PointVelocity
592  << G4endl
593  << " at end (midpoint)= " << midPoint << G4endl;
594  G4cout << " Particle mass = " << midPoint.GetRestMass() << G4endl;
595 
596  G4EquationOfMotion *equation
597  = integrDriver->GetStepper()->GetEquationOfMotion();
598  ReportFieldValue( CurrentA_PointVelocity, "start", equation );
599  ReportFieldValue( midPoint, "midPoint", equation );
600  G4cout << " Original Start = "
601  << CurveStartPointVelocity << G4endl;
602  G4cout << " Original End = "
603  << CurveEndPointVelocity << G4endl;
604  G4cout << " Original TrialPoint = "
605  << TrialPoint << G4endl;
606  G4cout << " (this is STRICT mode) "
607  << " num Splits= " << numSplits;
608  G4cout << G4endl;
609  }
610 #endif
611 
612  *ptrInterMedFT[depth] = midPoint;
613  CurrentB_PointVelocity = midPoint;
614 
615  // Adjust 'SubStartPoint' to calculate the 'did_length' in next loop
616  //
617  SubStart_PointVelocity = CurrentA_PointVelocity;
618 
619  // Find new trial intersection point needed at start of the loop
620  //
621  G4ThreeVector Point_A = CurrentA_PointVelocity.GetPosition();
622  G4ThreeVector Point_B = CurrentB_PointVelocity.GetPosition();
623 
624  GetNavigatorFor()->LocateGlobalPointWithinVolume(Point_A);
625  G4bool Intersects_AB = IntersectChord(Point_A, Point_B,
626  NewSafety, previousSafety,
627  previousSftOrigin, distAB,
628  PointGe);
629  if( Intersects_AB )
630  {
631  last_AF_intersection = Intersects_AB;
632  CurrentE_Point = PointGe;
633  fin_section_depth[depth]=true;
634 
635  validIntersectP= true; // 'E' has been updated.
636  // validApproxIntPV= false; // 'F' is no longer valid, as E changed
637 
638  G4bool validNormalAB;
639  NormalAtEntry = GetSurfaceNormal( PointGe, validNormalAB );
640  validNormalAtE = validNormalAB;
641  }
642  else
643  {
644  // No intersection found for first part of curve
645  // (CurrentA,InterMedPoint[depth]). Go to the second part
646  //
647  Second_half = true;
648 
649  validIntersectP= false; // No new 'E' chord intersection found
650  // validApproxIntPV= false; // So also 'F' is invalid
651  }
652  } // if did_len
653 
654  unsigned int levelPops=0;
655 
656  G4bool unfinished = Second_half;
657  while ( unfinished && (depth>0) ) // Loop checking, 07.10.2016, J. Apostolakis
658  {
659  // Second part of curve (InterMed[depth],Intermed[depth-1]))
660  // On the depth-1 level normally we are on the 'second_half'
661 
662  levelPops++;
663 
664  // Find new trial intersection point needed at start of the loop
665  //
666  SubStart_PointVelocity = *ptrInterMedFT[depth];
667  CurrentA_PointVelocity = *ptrInterMedFT[depth];
668  CurrentB_PointVelocity = *ptrInterMedFT[depth-1];
669 
670  // Ensure that the new endpoints are not further apart in space
671  // than on the curve due to different errors in the integration
672  //
673  G4FieldTrack RevisedEndPointFT= CurrentB_PointVelocity;
674  G4int errorEndPt;
675  G4bool recalculatedB=
676  CheckAndReEstimateEndpoint( CurrentA_PointVelocity,
677  CurrentB_PointVelocity,
678  RevisedEndPointFT,
679  errorEndPt );
680  if( recalculatedB )
681  {
682  CurrentB_PointVelocity= RevisedEndPointFT; // Use it !
683 
684  if (depth==1)
685  {
686  recalculatedEndPoint = true;
687  IntersectedOrRecalculatedFT = RevisedEndPointFT;
688  // So that we can return it, if it is the endpoint!
689  }
690  }
691  if( errorEndPt > 1 ) // errorEndPt = 1 is milder, just: len(B)=len(A)
692  {
693  std::ostringstream errmsg;
694  errmsg << "Location: " << MethodName << "- Second-Half" << G4endl;
695  ReportReversedPoints(errmsg,
696  CurveStartPointVelocity, CurveEndPointVelocity,
697  NewSafety, fiEpsilonStep,
698  CurrentA_PointVelocity, CurrentA_PointVelocity,
699  SubStart_PointVelocity, CurrentE_Point,
700  ApproxIntersecPointV, substep_no, substep_no_p, depth);
701  G4Exception(MethodName, "GeomNav0003", FatalException, errmsg);
702  }
703  G4ThreeVector Point_A = CurrentA_PointVelocity.GetPosition();
704  G4ThreeVector Point_B = CurrentB_PointVelocity.GetPosition();
705  GetNavigatorFor()->LocateGlobalPointWithinVolume(Point_A);
706  G4bool Intersects_AB = IntersectChord(Point_A, Point_B, NewSafety,
707  previousSafety,
708  previousSftOrigin, distAB,
709  PointGe);
710  if( Intersects_AB )
711  {
712  last_AF_intersection = Intersects_AB;
713  CurrentE_Point = PointGe;
714 
715  validIntersectP= true; // 'E' has been updated.
716  // validApproxIntPV= false; // 'F' is no longer valid, as E changed
717 
718  G4bool validNormalAB;
719  NormalAtEntry = GetSurfaceNormal( PointGe, validNormalAB );
720  validNormalAtE = validNormalAB;
721  }
722  else
723  {
724  validIntersectP= false; // No new 'E' chord intersection found
725  // validApproxIntPV= false; // So also 'F' is invalid
726  if( depth == 1)
727  {
728  there_is_no_intersection = true;
729  }
730  }
731  depth--;
732  fin_section_depth[depth]=true;
733  unfinished = !validIntersectP;
734  }
735 #ifdef G4DEBUG_FIELD
736  if( ! ( validIntersectP || there_is_no_intersection ) )
737  {
738  // What happened ??
739  G4cout << "MLL - WARNING Potential FAILURE: Conditions not met!"
740  << G4endl
741  << " Depth = " << depth << G4endl
742  << " Levels popped = " << levelPops
743  << " Num Substeps= " << substep_no << G4endl;
744  G4cout << " Found intersection= " << found_approximate_intersection
745  << G4endl;
746  G4cout << " Progress report: -- " << G4endl;
747  ReportProgress(G4cout,
748  CurveStartPointVelocity, CurveEndPointVelocity,
749  substep_no, CurrentA_PointVelocity,
750  CurrentB_PointVelocity,
751  NewSafety, depth );
752  G4cout << G4endl;
753  }
754 #endif
755  } // if(!found_aproximate_intersection)
756 
757  assert( validIntersectP || there_is_no_intersection
758  || found_approximate_intersection);
759 
760  } while ( ( ! found_approximate_intersection )
761  && ( ! there_is_no_intersection )
762  && ( substep_no <= fMaxSteps) ); // UNTIL found or failed
763 
764  if( substep_no > max_no_seen )
765  {
766  max_no_seen = substep_no;
767 #ifdef G4DEBUG_FIELD
768  if( max_no_seen > fWarnSteps )
769  {
770  trigger_substepno_print = max_no_seen-20; // Want to see last 20 steps
771  }
772 #endif
773  }
774 
775  if( !there_is_no_intersection && !found_approximate_intersection )
776  {
777  if( substep_no >= fMaxSteps)
778  {
779  // Since we cannot go further (yet), we return as far as we have gone
780 
781  recalculatedEndPoint = true;
782  IntersectedOrRecalculatedFT = CurrentA_PointVelocity;
783  found_approximate_intersection = false;
784 
785  std::ostringstream message;
786  message << G4endl;
787  message << "Convergence is requiring too many substeps: "
788  << substep_no << " ( limit = "<< fMaxSteps << ")"
789  << G4endl
790  << " Abandoning effort to intersect. " << G4endl << G4endl;
791  message << " Number of calls to MLL: " << fNumCalls;
792  message << " iteration = " << substep_no <<G4endl << G4endl;
793 
794  message.precision( 12 );
795  G4double done_len = CurrentA_PointVelocity.GetCurveLength();
796  G4double full_len = CurveEndPointVelocity.GetCurveLength();
797  message << " Undertaken only length: " << done_len
798  << " out of " << full_len << " required." << G4endl
799  << " Remaining length = " << full_len - done_len;
800 
801  message << " Start and end-point of requested Step:" << G4endl;
802  printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
803  -1.0, NewSafety, 0, message, -1 );
804  message << " Start and end-point of current Sub-Step:" << G4endl;
805  printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
806  -1.0, NewSafety, substep_no-1, message, -1 );
807  printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
808  -1.0, NewSafety, substep_no, message, -1 );
809 
810  G4Exception(MethodName, "GeomNav0003", JustWarning, message);
811  }
812  else if( substep_no >= fWarnSteps)
813  {
814  std::ostringstream message;
815  message << "Many substeps while trying to locate intersection."
816  << G4endl
817  << " Undertaken length: "
818  << CurrentB_PointVelocity.GetCurveLength()
819  << " - Needed: " << substep_no << " substeps." << G4endl
820  << " Warning number = " << fWarnSteps
821  << " and maximum substeps = " << fMaxSteps;
822  G4Exception(MethodName, "GeomNav1002", JustWarning, message);
823  }
824  }
825 
826 #ifdef G4DEBUG_FIELD
827  if( found_approximate_intersection )
828  {
829  assert( validApproxIntPV &&
830  "Approximate Intersection must not have been invalidated." );
831  }
832 #endif
833 
834  return (!there_is_no_intersection) && found_approximate_intersection;
835  // Success or failure
836 }
G4FieldTrack::GetCurveLength
G4double GetCurveLength() const
CLHEP::Hep3Vector::dot
double dot(const Hep3Vector &) const
G4MagInt_Driver::GetStepper
const G4MagIntegratorStepper * GetStepper() const
G4VIntersectionLocator::GetSurfaceNormal
G4ThreeVector GetSurfaceNormal(const G4ThreeVector &CurrentInt_Point, G4bool &validNormal)
Definition: G4VIntersectionLocator.cc:522
G4VIntersectionLocator::IntersectChord
G4bool IntersectChord(const G4ThreeVector &StartPointA, const G4ThreeVector &EndPointB, G4double &NewSafety, G4double &PreviousSafety, G4ThreeVector &PreviousSftOrigin, G4double &LinearStepLength, G4ThreeVector &IntersectionPoint, G4bool *calledNavigator=0)
G4ChordFinder::ApproxCurvePointV
G4FieldTrack ApproxCurvePointV(const G4FieldTrack &curveAPointVelocity, const G4FieldTrack &curveBPointVelocity, const G4ThreeVector &currentEPoint, G4double epsStep)
Definition: G4ChordFinder.cc:549
G4ThreadLocal
#define G4ThreadLocal
Definition: tls.hh:89
JustWarning
Definition: G4ExceptionSeverity.hh:64
G4int
int G4int
Definition: G4Types.hh:78
G4VIntersectionLocator::printStatus
void printStatus(const G4FieldTrack &startFT, const G4FieldTrack &currentFT, G4double requestStep, G4double safety, G4int stepNum)
Definition: G4VIntersectionLocator.cc:77
CLHEP::perThousand
static constexpr double perThousand
Definition: SystemOfUnits.h:313
G4VIntersectionLocator::GetEpsilonStepFor
G4double GetEpsilonStepFor()
G4VIntersectionLocator::ReportImmediateHit
void ReportImmediateHit(const char *MethodName, const G4ThreeVector &StartPosition, const G4ThreeVector &TrialPoint, double tolerance, unsigned long int numCalls)
Definition: G4VIntersectionLocator.cc:849
G4MagIntegratorStepper::GetEquationOfMotion
G4EquationOfMotion * GetEquationOfMotion()
G4FieldTrack::GetPosition
G4ThreeVector GetPosition() const
G4cout
G4GLOB_DLL std::ostream G4cout
G4bool
bool G4bool
Definition: G4Types.hh:79
CLHEP::mm
static constexpr double mm
Definition: SystemOfUnits.h:95
G4MagInt_Driver
Definition: G4MagIntegratorDriver.hh:48
G4FieldTrack::GetRestMass
G4double GetRestMass() const
Definition: G4FieldTrack.hh:126
G4VIntersectionLocator::ReportReversedPoints
void ReportReversedPoints(std::ostringstream &ossMsg, const G4FieldTrack &StartPointVel, const G4FieldTrack &EndPointVel, G4double NewSafety, G4double epsStep, const G4FieldTrack &CurrentA_PointVelocity, const G4FieldTrack &CurrentB_PointVelocity, const G4FieldTrack &SubStart_PointVelocity, const G4ThreeVector &CurrentE_Point, const G4FieldTrack &ApproxIntersecPointV, G4int sbstp_no, G4int sbstp_no_p, G4int depth)
Definition: G4VIntersectionLocator.cc:768
G4VIntersectionLocator::ReportTrialStep
void ReportTrialStep(G4int step_no, const G4ThreeVector &ChordAB_v, const G4ThreeVector &ChordEF_v, const G4ThreeVector &NewMomentumDir, const G4ThreeVector &NormalAtEntry, G4bool validNormal)
Definition: G4VIntersectionLocator.cc:624
G4VIntersectionLocator::GetNavigatorFor
G4Navigator * GetNavigatorFor()
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4INCL::Math::max
T max(const T t1, const T t2)
brief Return the largest of the two arguments
Definition: G4INCLGlobals.hh:112
FatalException
Definition: G4ExceptionSeverity.hh:60
CLHEP::Hep3Vector::mag2
double mag2() const
G4VIntersectionLocator::CheckAndReEstimateEndpoint
G4bool CheckAndReEstimateEndpoint(const G4FieldTrack &CurrentStartA, const G4FieldTrack &EstimatedEndB, G4FieldTrack &RevisedEndPoint, G4int &errorCode)
Definition: G4VIntersectionLocator.cc:327
G4VIntersectionLocator::GetAdjustementOfFoundIntersection
G4bool GetAdjustementOfFoundIntersection()
G4endl
#define G4endl
Definition: G4ios.hh:61
sqr
T sqr(const T &x)
Definition: templates.hh:145
G4double
double G4double
Definition: G4Types.hh:76
G4VIntersectionLocator::ReportProgress
void ReportProgress(std::ostream &oss, const G4FieldTrack &StartPointVel, const G4FieldTrack &EndPointVel, G4int substep_no, const G4FieldTrack &A_PtVel, const G4FieldTrack &B_PtVel, G4double safetyLast, G4int depth=-1)
Definition: G4VIntersectionLocator.cc:815
G4ChordFinder::GetIntegrationDriver
G4MagInt_Driver * GetIntegrationDriver()
CLHEP::Hep3Vector::mag
double mag() const
G4VIntersectionLocator::GetChordFinderFor
G4ChordFinder * GetChordFinderFor()
G4MagInt_Driver::AccurateAdvance
G4bool AccurateAdvance(G4FieldTrack &y_current, G4double hstep, G4double eps, G4double hinitial=0.0)
Definition: G4MagIntegratorDriver.cc:127
G4Navigator::LocateGlobalPointWithinVolume
virtual void LocateGlobalPointWithinVolume(const G4ThreeVector &position)
Definition: G4Navigator.cc:582
G4VIntersectionLocator::AdjustmentOfFoundIntersection
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)
Definition: G4VIntersectionLocator.cc:442

Here is the call graph for this function:

void G4MultiLevelLocator::ReportStatistics ( )

Definition at line 838 of file G4MultiLevelLocator.cc.

839 {
840  G4cout << " Number of calls = " << fNumCalls << G4endl;
841  G4cout << " Number of split level ('advances'): "
842  << fNumAdvanceTrials << G4endl;
843  G4cout << " Number of full advances: "
844  << fNumAdvanceGood << G4endl;
845  G4cout << " Number of good advances: "
846  << fNumAdvanceFull << G4endl;
847 }
G4cout
G4GLOB_DLL std::ostream G4cout
G4endl
#define G4endl
Definition: G4ios.hh:61

Here is the caller graph for this function:

void G4MultiLevelLocator::SetMaxSteps ( unsigned int  valMax)
inline

Definition at line 75 of file G4MultiLevelLocator.hh.

75 { fMaxSteps= valMax; }

Here is the caller graph for this function:

void G4MultiLevelLocator::SetWarnSteps ( unsigned int  valWarn)
inline

Definition at line 76 of file G4MultiLevelLocator.hh.

76 { fWarnSteps= valWarn; }

Here is the caller graph for this function:

The documentation for this class was generated from the following files: