Geant4  10.02.p02
G4MultiLevelLocator.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: G4MultiLevelLocator.cc 97507 2016-06-03 12:48:42Z gcosmo $
27 //
28 // Class G4MultiLevelLocator implementation
29 //
30 // 27.10.08 - Tatiana Nikitina.
31 // 04.10.11 - John Apostolakis, revised convergence to use Surface Normal
32 // ---------------------------------------------------------------------------
33 
34 #include <iomanip>
35 
36 #include "G4ios.hh"
37 #include "G4MultiLevelLocator.hh"
38 
39 
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 }
64 
66 {
67  for ( G4int idepth=0; idepth<max_depth+1; idepth++)
68  {
69  delete ptrInterMedFT[idepth];
70  }
71 #ifdef G4DEBUG_FIELD
73 #endif
74 }
75 
76 
77 // --------------------------------------------------------------------------
78 // G4bool G4PropagatorInField::LocateIntersectionPoint(
79 // const G4FieldTrack& CurveStartPointVelocity, // A
80 // const G4FieldTrack& CurveEndPointVelocity, // B
81 // const G4ThreeVector& TrialPoint, // E
82 // G4FieldTrack& IntersectedOrRecalculated // Output
83 // G4bool& recalculated ) // Out
84 // --------------------------------------------------------------------------
85 //
86 // Function that returns the intersection of the true path with the surface
87 // of the current volume (either the external one or the inner one with one
88 // of the daughters:
89 //
90 // A = Initial point
91 // B = another point
92 //
93 // Both A and B are assumed to be on the true path:
94 //
95 // E is the first point of intersection of the chord AB with
96 // a volume other than A (on the surface of A or of a daughter)
97 //
98 // Convention of Use :
99 // i) If it returns "true", then IntersectionPointVelocity is set
100 // to the approximate intersection point.
101 // ii) If it returns "false", no intersection was found.
102 // Potential reasons:
103 // a) no segment found an intersection
104 // b) too many steps were required - after that it abandoned the effort
105 // and is returning how far it could go. (New - 29 Oct 2015)
106 // (If so, it must set 'recalculated' to true.)
107 // TODO/idea: add a new flag: 'unfinished' to identify these cases.
108 //
109 // IntersectedOrRecalculated means different things:
110 // a) if it is the same curve lenght along, it is a revision of the
111 // original enpdoint due to the need for re-integration.
112 // b) if it is at a shorter curve length, it is 'end of what it could do'
113 // i.e. as far as it could go, because it took too many steps!
114 // Note: IntersectedOrRecalculated is valid only if 'recalculated' is
115 // 'true'.
116 // --------------------------------------------------------------------------
117 // NOTE: implementation taken from G4PropagatorInField
118 //
120  const G4FieldTrack& CurveStartPointVelocity, // A
121  const G4FieldTrack& CurveEndPointVelocity, // B
122  const G4ThreeVector& TrialPoint, // E
123  G4FieldTrack& IntersectedOrRecalculatedFT, // Output
124  G4bool& recalculatedEndPoint, // Out
125  G4double& previousSafety, // In/Out
126  G4ThreeVector& previousSftOrigin) // In/Out
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
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  do // REPEAT param
230  {
231  G4ThreeVector Point_A = CurrentA_PointVelocity.GetPosition();
232  G4ThreeVector Point_B = CurrentB_PointVelocity.GetPosition();
233 
234  // F = a point on true AB path close to point E
235  // (the closest if possible)
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  // validApproxIntPV = true;
245 
246 #ifdef G4DEBUG_FIELD
247  if( ApproxIntersecPointV.GetCurveLength() >
248  CurrentB_PointVelocity.GetCurveLength() * (1.0 + tolerance) )
249  {
250  G4Exception(MethodName, "GeomNav0003", FatalException,
251  "Intermediate F point is past end B point" );
252  }
253 #endif
254 
255  G4ThreeVector CurrentF_Point= ApproxIntersecPointV.GetPosition();
256 
257  // First check whether EF is small - then F is a good approx. point
258  // Calculate the length and direction of the chord AF
259  //
260  G4ThreeVector ChordEF_Vector = CurrentF_Point - CurrentE_Point;
261 
262  G4ThreeVector NewMomentumDir= ApproxIntersecPointV.GetMomentumDir();
263  G4double MomDir_dot_Norm= NewMomentumDir.dot( NormalAtEntry ) ;
264 
265 #ifdef G4DEBUG_FIELD
266  if( fVerboseLevel > 3 )
267  {
268  G4ThreeVector ChordAB = Point_B - Point_A;
269  G4double ABchord_length = ChordAB.mag();
270  G4double MomDir_dot_ABchord;
271  MomDir_dot_ABchord = (1.0 / ABchord_length)
272  * NewMomentumDir.dot( ChordAB );
273  G4VIntersectionLocator::ReportTrialStep( substep_no, ChordAB,
274  ChordEF_Vector, NewMomentumDir, NormalAtEntry, validNormalAtE );
275  G4cout << " dot( MomentumDir, ABchord_unit ) = "
276  << MomDir_dot_ABchord << G4endl;
277  }
278 #endif
279  G4bool adequate_angle =
280  ( MomDir_dot_Norm >= 0.0 ) // Can use ( > -epsilon) instead
281  || (! validNormalAtE ); // Invalid, cannot use this criterion
282  G4double EF_dist2 = ChordEF_Vector.mag2();
283  if ( ( EF_dist2 <= sqr(fiDeltaIntersection) && ( adequate_angle ) )
284  || ( EF_dist2 <= kCarTolerance*kCarTolerance ) )
285  {
286  found_approximate_intersection = true;
287 
288  // Create the "point" return value
289  //
290  IntersectedOrRecalculatedFT = ApproxIntersecPointV;
291  IntersectedOrRecalculatedFT.SetPosition( CurrentE_Point );
292 
294  {
295  // Try to Get Correction of IntersectionPoint using SurfaceNormal()
296  //
297  G4ThreeVector IP;
298  G4ThreeVector MomentumDir=ApproxIntersecPointV.GetMomentumDirection();
299  G4bool goodCorrection = AdjustmentOfFoundIntersection(Point_A,
300  CurrentE_Point, CurrentF_Point, MomentumDir,
301  last_AF_intersection, IP, NewSafety,
302  previousSafety, previousSftOrigin );
303  if ( goodCorrection )
304  {
305  IntersectedOrRecalculatedFT = ApproxIntersecPointV;
306  IntersectedOrRecalculatedFT.SetPosition(IP);
307  }
308  }
309  // Note: in order to return a point on the boundary,
310  // we must return E. But it is F on the curve.
311  // So we must "cheat": we are using the position at point E
312  // and the velocity at point F !!!
313  //
314  // This must limit the length we can allow for displacement!
315  }
316  else // E is NOT close enough to the curve (ie point F)
317  {
318  // Check whether any volumes are encountered by the chord AF
319  // ---------------------------------------------------------
320  // First relocate to restore any Voxel etc information
321  // in the Navigator before calling ComputeStep()
322  //
324 
325  G4ThreeVector PointG; // Candidate intersection point
326  G4double stepLengthAF;
327  G4bool Intersects_AF = IntersectChord( Point_A, CurrentF_Point,
328  NewSafety, previousSafety,
329  previousSftOrigin,
330  stepLengthAF,
331  PointG );
332  last_AF_intersection = Intersects_AF;
333  if( Intersects_AF )
334  {
335  // G is our new Candidate for the intersection point.
336  // It replaces "E" and we will repeat the test to see if
337  // it is a good enough approximate point for us.
338  // B <- F
339  // E <- G
340  //
341  CurrentB_PointVelocity = ApproxIntersecPointV;
342  CurrentE_Point = PointG;
343 
344  validIntersectP= true; // 'E' has been updated.
345  // validApproxIntPV= false; // 'F' is no longer valid, as B changed
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 
356 #ifdef G4VERBOSE
357  if( fVerboseLevel > 3 )
358  {
359  G4cout << "G4PiF::LI> Investigating intermediate point"
360  << " at s=" << ApproxIntersecPointV.GetCurveLength()
361  << " on way to full s="
362  << CurveEndPointVelocity.GetCurveLength() << G4endl;
363  }
364 #endif
365  }
366  else // not Intersects_AF
367  {
368  // In this case:
369  // There is NO intersection of AF with a volume boundary.
370  // We must continue the search in the segment FB!
371  //
372  GetNavigatorFor()->LocateGlobalPointWithinVolume( CurrentF_Point );
373 
374  G4double stepLengthFB;
375  G4ThreeVector PointH;
376 
377  // Check whether any volumes are encountered by the chord FB
378  // ---------------------------------------------------------
379 
380  G4bool Intersects_FB = IntersectChord( CurrentF_Point, Point_B,
381  NewSafety, previousSafety,
382  previousSftOrigin,
383  stepLengthFB,
384  PointH );
385  if( Intersects_FB )
386  {
387  // There is an intersection of FB with a volume boundary
388  // H <- First Intersection of Chord FB
389 
390  // H is our new Candidate for the intersection point.
391  // It replaces "E" and we will repeat the test to see if
392  // it is a good enough approximate point for us.
393 
394  // Note that F must be in volume volA (the same as A)
395  // (otherwise AF would meet a volume boundary!)
396  // A <- F
397  // E <- H
398  //
399  CurrentA_PointVelocity = ApproxIntersecPointV;
400  CurrentE_Point = PointH;
401 
402  validIntersectP = true; // 'E' has been updated.
403  // validApproxIntPV = false; // 'F' is no longer valid, as A changed
404 
405  G4bool validNormalH;
406  NormalAtEntry = GetSurfaceNormal( PointH, validNormalH );
407  validNormalAtE = validNormalH;
408  }
409  else // not Intersects_FB
410  {
411  if(fin_section_depth[depth])
412  {
413  // If B is the original endpoint, this means that whatever
414  // volume(s) intersected the original chord, none touch the
415  // smaller chords we have used.
416  // The value of 'IntersectedOrRecalculatedFT' returned is
417  // likely not valid
418 
419  // Check on real final_section or SubEndSection
420  //
421  if( ((Second_half)&&(depth==0)) || (first_section) )
422  {
423  there_is_no_intersection = true; // real final_section
424  }
425  else
426  {
427  // end of subsection, not real final section
428  // exit from the and go to the depth-1 level
429  substep_no_p = param_substeps+2; // exit from the loop
430 
431  // but 'Second_half' is still true because we need to find
432  // the 'CurrentE_point' for the next loop
433  Second_half = true;
434  sub_final_section = true;
435  }
436  }
437  else
438  {
439  CurrentA_PointVelocity = CurrentB_PointVelocity; // Got to B
440  CurrentB_PointVelocity = (depth==0) ? CurveEndPointVelocity
441  : *ptrInterMedFT[depth] ;
442  SubStart_PointVelocity = CurrentA_PointVelocity;
443  restoredFullEndpoint = true;
444 
445  validIntersectP= false; // 'E' has NOT been updated.
446  // validApproxIntPV= false; // 'F' is no longer valid, A changed
447  }
448  } // Endif (Intersects_FB)
449  } // Endif (Intersects_AF)
450 
451  G4FieldTrack RevisedB_FT= CurrentB_PointVelocity;
452  G4int errorEndPt;
453  G4bool recalculatedB= CheckAndReEstimateEndpoint(CurrentA_PointVelocity,
454  CurrentB_PointVelocity,
455  RevisedB_FT,
456  errorEndPt );
457  if( recalculatedB )
458  {
459  CurrentB_PointVelocity= RevisedB_FT; // Use it !
460  // Do not invalidate intersection F -- it is still roughly OK.
461  //
462  // The best course would be to invalidate (reset validIntersectP)
463  // BUT if we invalidate it, we must re-estimate it somewhere!
464 
465  // validApproxIntPV= false; // 'F' is no longer valid, as B changed
466  // validIntersectP= false; // 'E' has NOT been updated.
467 
468  if ( (fin_section_depth[depth]) // real final section
469  &&( first_section || ((Second_half)&&(depth==0)) ) )
470  {
471  recalculatedEndPoint = true;
472  IntersectedOrRecalculatedFT = RevisedB_FT;
473  // So that we can return it, if it is the endpoint!
474  }
475  // else
476  // Move forward the other points
477  // - or better flag it, so that they are re-computed when next used
478  // [ Implementation: a counter for # of recomputations
479  // => avoids extra work]
480 
481  }
482  if( errorEndPt > 1 ) // errorEndPt = 1 is milder, just: len(B)=len(A)
483  {
484  std::ostringstream errmsg;
485  errmsg << "Location: " << MethodName
486  << "- After EndIf(Intersects_AF)" << G4endl;
487  ReportReversedPoints(errmsg,
488  CurveStartPointVelocity, CurveEndPointVelocity,
489  NewSafety, fiEpsilonStep,
490  CurrentA_PointVelocity, CurrentB_PointVelocity,
491  SubStart_PointVelocity, CurrentE_Point,
492  ApproxIntersecPointV, substep_no, substep_no_p, depth);
493  G4Exception(MethodName, "GeomNav0003", FatalException, errmsg);
494  }
495  if( restoredFullEndpoint )
496  {
497  fin_section_depth[depth] = restoredFullEndpoint;
498  restoredFullEndpoint = false;
499  }
500  } // EndIf ( E is close enough to the curve, ie point F. )
501  // tests ChordAF_Vector.mag() <= maximum_lateral_displacement
502 
503 #ifdef G4DEBUG_FIELD
504  if( trigger_substepno_print == 0)
505  {
506  trigger_substepno_print= fWarnSteps - 20;
507  }
508 
509  if( substep_no >= trigger_substepno_print )
510  {
511  G4cout << "Difficulty in converging in " << MethodName
512  << G4endl
513  << " Substep no = " << substep_no << G4endl;
514  if( substep_no == trigger_substepno_print )
515  {
516  printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
517  -1.0, NewSafety, 0);
518  }
519  G4cout << " State of point A: ";
520  printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
521  -1.0, NewSafety, substep_no-1);
522  G4cout << " State of point B: ";
523  printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
524  -1.0, NewSafety, substep_no);
525  }
526 #endif
527  substep_no++;
528  substep_no_p++;
529 
530  } while ( ( ! found_approximate_intersection )
531  && ( ! there_is_no_intersection )
532  && ( substep_no_p <= param_substeps) ); // UNTIL found or
533  // failed param substep
534 
535  if( (!found_approximate_intersection) && (!there_is_no_intersection) )
536  {
537  G4double did_len = std::abs( CurrentA_PointVelocity.GetCurveLength()
538  - SubStart_PointVelocity.GetCurveLength());
539  G4double all_len = std::abs( CurrentB_PointVelocity.GetCurveLength()
540  - SubStart_PointVelocity.GetCurveLength());
541 
542  G4double distAB= -1;
543  G4ThreeVector PointGe;
544  //
545  // Is progress is too slow, and is it possible to go deeper?
546  // If so, then *halve the step*
547  // ==============
548  if( (did_len < fraction_done*all_len)
549  && (depth<max_depth) && (!sub_final_section) )
550  {
551 #ifdef G4DEBUG_FIELD
552  static G4ThreadLocal unsigned int numSplits=0; // For debugging only
553  biggest_depth= std::max(depth, biggest_depth);
554  numSplits++;
555 #endif
556  Second_half=false;
557  depth++;
558  first_section = false;
559 
560  G4double Sub_len = (all_len-did_len)/(2.);
561  G4FieldTrack midPoint = CurrentA_PointVelocity;
562  G4MagInt_Driver* integrDriver
564  G4bool fullAdvance=
565  integrDriver->AccurateAdvance(midPoint, Sub_len, fiEpsilonStep);
566 
568  if( fullAdvance ) { fNumAdvanceFull++; }
569 
570  G4double lenAchieved=
571  midPoint.GetCurveLength()-CurrentA_PointVelocity.GetCurveLength();
572 
573  const G4double adequateFraction = (1.0-CLHEP::perThousand);
574  G4bool goodAdvance = (lenAchieved >= adequateFraction * Sub_len);
575  if ( goodAdvance ) { fNumAdvanceGood++; }
576 
577 #ifdef G4DEBUG_FIELD
578  else // !goodAdvance
579  {
580  G4cout << "MLL> AdvanceChordLimited not full at depth=" << depth
581  << " total length achieved = " << lenAchieved << " of "
582  << Sub_len << " fraction= ";
583  if (Sub_len != 0.0 ) { G4cout << lenAchieved / Sub_len; }
584  else { G4cout << "DivByZero"; }
585  G4cout << " Good-enough fraction = " << adequateFraction;
586  G4cout << " Number of call to mll = " << fNumCalls
587  << " iteration " << substep_no
588  << " inner = " << substep_no_p << G4endl;
589  G4cout << " Epsilon of integration = " << fiEpsilonStep << G4endl;
590  G4cout << " State at start is = " << CurrentA_PointVelocity
591  << G4endl
592  << " at end (midpoint)= " << midPoint << G4endl;
593  G4cout << " Particle mass = " << midPoint.GetRestMass() << G4endl;
594 
595  G4EquationOfMotion *equation
596  = integrDriver->GetStepper()->GetEquationOfMotion();
597  ReportFieldValue( CurrentA_PointVelocity, "start", equation );
598  ReportFieldValue( midPoint, "midPoint", equation );
599  G4cout << " Original Start = "
600  << CurveStartPointVelocity << G4endl;
601  G4cout << " Original End = "
602  << CurveEndPointVelocity << G4endl;
603  G4cout << " Original TrialPoint = "
604  << TrialPoint << G4endl;
605  G4cout << " (this is STRICT mode) "
606  << " num Splits= " << numSplits;
607  G4cout << G4endl;
608  }
609 #endif
610 
611  *ptrInterMedFT[depth] = midPoint;
612  CurrentB_PointVelocity = midPoint;
613 
614  // Adjust 'SubStartPoint' to calculate the 'did_length' in next loop
615  //
616  SubStart_PointVelocity = CurrentA_PointVelocity;
617 
618  // Find new trial intersection point needed at start of the loop
619  //
620  G4ThreeVector Point_A = CurrentA_PointVelocity.GetPosition();
621  G4ThreeVector Point_B = CurrentB_PointVelocity.GetPosition();
622 
624  G4bool Intersects_AB = IntersectChord(Point_A, Point_B,
625  NewSafety, previousSafety,
626  previousSftOrigin, distAB,
627  PointGe);
628  if( Intersects_AB )
629  {
630  last_AF_intersection = Intersects_AB;
631  CurrentE_Point = PointGe;
632  fin_section_depth[depth]=true;
633 
634  validIntersectP= true; // 'E' has been updated.
635  // validApproxIntPV= false; // 'F' is no longer valid, as E changed
636 
637  G4bool validNormalAB;
638  NormalAtEntry = GetSurfaceNormal( PointGe, validNormalAB );
639  validNormalAtE = validNormalAB;
640  }
641  else
642  {
643  // No intersection found for first part of curve
644  // (CurrentA,InterMedPoint[depth]). Go to the second part
645  //
646  Second_half = true;
647 
648  validIntersectP= false; // No new 'E' chord intersection found
649  // validApproxIntPV= false; // So also 'F' is invalid
650  }
651  } // if did_len
652 
653  unsigned int levelPops=0;
654 
655  G4bool unfinished = Second_half;
656  while ( unfinished && (depth>0) )
657  {
658  // Second part of curve (InterMed[depth],Intermed[depth-1]))
659  // On the depth-1 level normally we are on the 'second_half'
660 
661  levelPops++;
662 
663  // Find new trial intersection point needed at start of the loop
664  //
665  SubStart_PointVelocity = *ptrInterMedFT[depth];
666  CurrentA_PointVelocity = *ptrInterMedFT[depth];
667  CurrentB_PointVelocity = *ptrInterMedFT[depth-1];
668 
669  // Ensure that the new endpoints are not further apart in space
670  // than on the curve due to different errors in the integration
671  //
672  G4FieldTrack RevisedEndPointFT= CurrentB_PointVelocity;
673  G4int errorEndPt;
674  G4bool recalculatedB=
675  CheckAndReEstimateEndpoint( CurrentA_PointVelocity,
676  CurrentB_PointVelocity,
677  RevisedEndPointFT,
678  errorEndPt );
679  if( recalculatedB )
680  {
681  CurrentB_PointVelocity= RevisedEndPointFT; // Use it !
682 
683  if (depth==1)
684  {
685  recalculatedEndPoint = true;
686  IntersectedOrRecalculatedFT = RevisedEndPointFT;
687  // So that we can return it, if it is the endpoint!
688  }
689  }
690  if( errorEndPt > 1 ) // errorEndPt = 1 is milder, just: len(B)=len(A)
691  {
692  std::ostringstream errmsg;
693  errmsg << "Location: " << MethodName << "- Second-Half" << G4endl;
694  ReportReversedPoints(errmsg,
695  CurveStartPointVelocity, CurveEndPointVelocity,
696  NewSafety, fiEpsilonStep,
697  CurrentA_PointVelocity, CurrentA_PointVelocity,
698  SubStart_PointVelocity, CurrentE_Point,
699  ApproxIntersecPointV, substep_no, substep_no_p, depth);
700  G4Exception(MethodName, "GeomNav0003", FatalException, errmsg);
701  }
702  G4ThreeVector Point_A = CurrentA_PointVelocity.GetPosition();
703  G4ThreeVector Point_B = CurrentB_PointVelocity.GetPosition();
705  G4bool Intersects_AB = IntersectChord(Point_A, Point_B, NewSafety,
706  previousSafety,
707  previousSftOrigin, distAB,
708  PointGe);
709  if( Intersects_AB )
710  {
711  last_AF_intersection = Intersects_AB;
712  CurrentE_Point = PointGe;
713 
714  validIntersectP= true; // 'E' has been updated.
715  // validApproxIntPV= false; // 'F' is no longer valid, as E changed
716 
717  G4bool validNormalAB;
718  NormalAtEntry = GetSurfaceNormal( PointGe, validNormalAB );
719  validNormalAtE = validNormalAB;
720  }
721  else
722  {
723  validIntersectP= false; // No new 'E' chord intersection found
724  // validApproxIntPV= false; // So also 'F' is invalid
725  if( depth == 1)
726  {
727  there_is_no_intersection = true;
728  }
729  }
730  depth--;
731  fin_section_depth[depth]=true;
732  unfinished = !validIntersectP;
733  }
734 #ifdef G4DEBUG_FIELD
735  if( ! ( validIntersectP || there_is_no_intersection ) )
736  {
737  // What happened ??
738  G4cout << "MLL - WARNING Potential FAILURE: Conditions not met!"
739  << G4endl
740  << " Depth = " << depth << G4endl
741  << " Levels popped = " << levelPops
742  << " Num Substeps= " << substep_no << G4endl;
743  G4cout << " Found intersection= " << found_approximate_intersection
744  << G4endl;
745  G4cout << " Progress report: -- " << G4endl;
747  CurveStartPointVelocity, CurveEndPointVelocity,
748  substep_no, CurrentA_PointVelocity,
749  CurrentB_PointVelocity,
750  NewSafety, depth );
751  G4cout << G4endl;
752  }
753 #endif
754  } // if(!found_aproximate_intersection)
755 
756  assert( validIntersectP || there_is_no_intersection
757  || found_approximate_intersection);
758 
759  } while ( ( ! found_approximate_intersection )
760  && ( ! there_is_no_intersection )
761  && ( substep_no <= fMaxSteps) ); // UNTIL found or failed
762 
763  if( substep_no > max_no_seen )
764  {
765  max_no_seen = substep_no;
766 #ifdef G4DEBUG_FIELD
767  if( max_no_seen > fWarnSteps )
768  {
769  trigger_substepno_print = max_no_seen-20; // Want to see last 20 steps
770  }
771 #endif
772  }
773 
774  if( !there_is_no_intersection && !found_approximate_intersection )
775  {
776  if( substep_no >= fMaxSteps)
777  {
778  // Since we cannot go further (yet), we return as far as we have gone
779 
780  recalculatedEndPoint = true;
781  IntersectedOrRecalculatedFT = CurrentA_PointVelocity;
782  found_approximate_intersection = false;
783 
784  std::ostringstream message;
785  message << G4endl;
786  message << "Convergence is requiring too many substeps: "
787  << substep_no << " ( limit = "<< fMaxSteps << ")"
788  << G4endl
789  << " Abandoning effort to intersect. " << G4endl << G4endl;
790  message << " Number of calls to MLL: " << fNumCalls;
791  message << " iteration = " << substep_no <<G4endl << G4endl;
792 
793  message.precision( 12 );
794  G4double done_len = CurrentA_PointVelocity.GetCurveLength();
795  G4double full_len = CurveEndPointVelocity.GetCurveLength();
796  message << " Undertaken only length: " << done_len
797  << " out of " << full_len << " required." << G4endl
798  << " Remaining length = " << full_len - done_len;
799 
800  message << " Start and end-point of requested Step:" << G4endl;
801  printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
802  -1.0, NewSafety, 0, message, -1 );
803  message << " Start and end-point of current Sub-Step:" << G4endl;
804  printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
805  -1.0, NewSafety, substep_no-1, message, -1 );
806  printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
807  -1.0, NewSafety, substep_no, message, -1 );
808 
809  G4Exception(MethodName, "GeomNav0003", JustWarning, message);
810  }
811  else if( substep_no >= fWarnSteps)
812  {
813  std::ostringstream message;
814  message << "Many substeps while trying to locate intersection."
815  << G4endl
816  << " Undertaken length: "
817  << CurrentB_PointVelocity.GetCurveLength()
818  << " - Needed: " << substep_no << " substeps." << G4endl
819  << " Warning number = " << fWarnSteps
820  << " and maximum substeps = " << fMaxSteps;
821  G4Exception(MethodName, "GeomNav1002", JustWarning, message);
822  }
823  }
824 
825 #ifdef G4DEBUG_FIELD
826  if( found_approximate_intersection )
827  {
828  assert( validApproxIntPV &&
829  "Approximate Intersection must not have been invalidated." );
830  }
831 #endif
832 
833  return (!there_is_no_intersection) && found_approximate_intersection;
834  // Success or failure
835 }
836 
838 {
839  G4cout << " Number of calls = " << fNumCalls << G4endl;
840  G4cout << " Number of split level ('advances'): "
842  G4cout << " Number of full advances: "
843  << fNumAdvanceGood << G4endl;
844  G4cout << " Number of good advances: "
845  << fNumAdvanceFull << G4endl;
846 }
847 
849  const char* nameLoc,
850  const G4EquationOfMotion* equation )
851 {
852  enum { maxNumFieldComp= 24 };
853 
854  G4ThreeVector position = locationPV.GetPosition();
855  G4double startPoint[4] = { position.x(), position.y(), position.z(),
856  locationPV.GetLabTimeOfFlight() };
857  G4double FieldVec[maxNumFieldComp]; // 24 ;
858  for (unsigned int i=0; i<maxNumFieldComp; ++i )
859  {
860  FieldVec[i]= 0.0;
861  }
862  equation->GetFieldValue( startPoint, FieldVec);
863  G4cout << " B-field value (" << nameLoc << ")= "
864  << FieldVec[0] << " " << FieldVec[1] << " " << FieldVec[2];
865  G4double Emag2= G4ThreeVector( FieldVec[3],
866  FieldVec[4],
867  FieldVec[5] ).mag2();
868  if( Emag2 > 0.0 )
869  {
870  G4cout << " Electric = " << FieldVec[3] << " "
871  << FieldVec[4] << " "
872  << FieldVec[5]<< G4endl;
873  }
874  return;
875 }
G4FieldTrack * ptrInterMedFT[max_depth+1]
void SetPosition(G4ThreeVector nPos)
void ReportFieldValue(const G4FieldTrack &locationPV, const char *nameLoc, const G4EquationOfMotion *equation)
G4double GetCurveLength() const
CLHEP::Hep3Vector G4ThreeVector
const G4MagIntegratorStepper * GetStepper() const
static const G4float tolerance
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
G4bool EstimateIntersectionPoint(const G4FieldTrack &curveStartPointTangent, const G4FieldTrack &curveEndPointTangent, const G4ThreeVector &trialPoint, G4FieldTrack &intersectPointTangent, G4bool &recalculatedEndPoint, G4double &fPreviousSafety, G4ThreeVector &fPreviousSftOrigin)
static const double perThousand
Definition: G4SIunits.hh:330
#define G4ThreadLocal
Definition: tls.hh:89
void SetWarnSteps(unsigned int valWarn)
int G4int
Definition: G4Types.hh:78
void printStatus(const G4FieldTrack &startFT, const G4FieldTrack &currentFT, G4double requestStep, G4double safety, G4int stepNum)
unsigned long int fNumAdvanceTrials
G4double GetEpsilonStepFor()
void ReportImmediateHit(const char *MethodName, const G4ThreeVector &StartPosition, const G4ThreeVector &TrialPoint, double tolerance, unsigned long int numCalls)
G4EquationOfMotion * GetEquationOfMotion()
#define position
Definition: xmlparse.cc:622
G4ThreeVector GetPosition() const
G4GLOB_DLL std::ostream G4cout
static const G4int max_depth
G4MultiLevelLocator(G4Navigator *theNavigator)
bool G4bool
Definition: G4Types.hh:79
G4double GetRestMass() const
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 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
void GetFieldValue(const G4double Point[4], G4double Field[]) const
G4double GetLabTimeOfFlight() const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4bool CheckAndReEstimateEndpoint(const G4FieldTrack &CurrentStartA, const G4FieldTrack &EstimatedEndB, G4FieldTrack &RevisedEndPoint, G4int &errorCode)
G4bool GetAdjustementOfFoundIntersection()
void SetMaxSteps(unsigned int valMax)
#define G4endl
Definition: G4ios.hh:61
T sqr(const T &x)
Definition: templates.hh:145
double G4double
Definition: G4Types.hh:76
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)
G4MagInt_Driver * GetIntegrationDriver()
unsigned long int fNumAdvanceGood
static const double mm
Definition: G4SIunits.hh:114
unsigned long int fNumAdvanceFull
G4ChordFinder * GetChordFinderFor()
unsigned long int fNumCalls
G4bool AccurateAdvance(G4FieldTrack &y_current, G4double hstep, G4double eps, G4double hinitial=0.0)
virtual void LocateGlobalPointWithinVolume(const G4ThreeVector &position)
Definition: G4Navigator.cc:582
G4bool AdjustmentOfFoundIntersection(const G4ThreeVector &A, const G4ThreeVector &CurrentE_Point, const G4ThreeVector &CurrentF_Point, const G4ThreeVector &MomentumDir, const G4bool IntersectAF, G4ThreeVector &IntersectionPoint, G4double &NewSafety, G4double &fPrevSafety, G4ThreeVector &fPrevSftOrigin)
G4ThreeVector GetMomentumDirection() const