ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4MultiLevelLocator.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4MultiLevelLocator.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // Class G4MultiLevelLocator implementation
27 //
28 // 27.10.08 - Tatiana Nikitina.
29 // 04.10.11 - John Apostolakis, revised convergence to use Surface Normal
30 // ---------------------------------------------------------------------------
31 
32 #include <iomanip>
33 
34 #include "G4ios.hh"
35 #include "G4MultiLevelLocator.hh"
36 
38  : G4VIntersectionLocator(theNavigator)
39 {
40  // In case of too slow progress in finding Intersection Point
41  // intermediates Points on the Track must be stored.
42  // Initialise the array of Pointers [max_depth+1] to do this
43 
44  G4ThreeVector zeroV(0.0,0.0,0.0);
45  for ( auto idepth=0; idepth<max_depth+1; ++idepth )
46  {
47  ptrInterMedFT[ idepth ] = new G4FieldTrack( zeroV, zeroV, 0., 0., 0., 0.);
48  }
49 
50 #ifdef G4DEBUG_FIELD
51  // Trial values Loose Tight
52  // To happen: Infrequent Often
53  SetMaxSteps(50); // 300 25
54  SetWarnSteps(40); // 250 15
55 #endif
56 }
57 
59 {
60  for ( auto idepth=0; idepth<max_depth+1; ++idepth )
61  {
62  delete ptrInterMedFT[idepth];
63  }
64 #ifdef G4DEBUG_FIELD
66 #endif
67 }
68 
69 
70 // --------------------------------------------------------------------------
71 // G4bool G4PropagatorInField::LocateIntersectionPoint(
72 // const G4FieldTrack& CurveStartPointVelocity, // A
73 // const G4FieldTrack& CurveEndPointVelocity, // B
74 // const G4ThreeVector& TrialPoint, // E
75 // G4FieldTrack& IntersectedOrRecalculated // Output
76 // G4bool& recalculated ) // Out
77 // --------------------------------------------------------------------------
78 //
79 // Function that returns the intersection of the true path with the surface
80 // of the current volume (either the external one or the inner one with one
81 // of the daughters:
82 //
83 // A = Initial point
84 // B = another point
85 //
86 // Both A and B are assumed to be on the true path:
87 //
88 // E is the first point of intersection of the chord AB with
89 // a volume other than A (on the surface of A or of a daughter)
90 //
91 // Convention of Use :
92 // i) If it returns "true", then IntersectionPointVelocity is set
93 // to the approximate intersection point.
94 // ii) If it returns "false", no intersection was found.
95 // Potential reasons:
96 // a) no segment found an intersection
97 // b) too many steps were required - after that it abandoned the effort
98 // and is returning how far it could go. (New - 29 Oct 2015)
99 // (If so, it must set 'recalculated' to true.)
100 // TODO/idea: add a new flag: 'unfinished' to identify these cases.
101 //
102 // IntersectedOrRecalculated means different things:
103 // a) if it is the same curve lenght along, it is a revision of the
104 // original enpdoint due to the need for re-integration.
105 // b) if it is at a shorter curve length, it is 'end of what it could do'
106 // i.e. as far as it could go, because it took too many steps!
107 // Note: IntersectedOrRecalculated is valid only if 'recalculated' is
108 // 'true'.
109 // --------------------------------------------------------------------------
110 // NOTE: implementation taken from G4PropagatorInField
111 //
113  const G4FieldTrack& CurveStartPointVelocity, // A
114  const G4FieldTrack& CurveEndPointVelocity, // B
115  const G4ThreeVector& TrialPoint, // E
116  G4FieldTrack& IntersectedOrRecalculatedFT, // Output
117  G4bool& recalculatedEndPoint, // Out
118  G4double& previousSafety, // In/Out
119  G4ThreeVector& previousSftOrigin) // In/Out
120 {
121  // Find Intersection Point ( A, B, E ) of true path AB - start at E.
122  const char* MethodName= "G4MultiLevelLocator::EstimateIntersectionPoint()";
123 
124  G4bool found_approximate_intersection = false;
125  G4bool there_is_no_intersection = false;
126 
127  G4FieldTrack CurrentA_PointVelocity = CurveStartPointVelocity;
128  G4FieldTrack CurrentB_PointVelocity = CurveEndPointVelocity;
129  G4ThreeVector CurrentE_Point = TrialPoint;
130  G4bool validNormalAtE = false;
131  G4ThreeVector NormalAtEntry;
132 
133  G4FieldTrack ApproxIntersecPointV(CurveEndPointVelocity); // FT-Def-Construct
134  G4bool validIntersectP= true; // Is it current ?
135  G4double NewSafety = 0.0;
136  G4bool last_AF_intersection = false;
137 
138  auto integrDriver = GetChordFinderFor()->GetIntegrationDriver();
139  G4bool driverReIntegrates = integrDriver->DoesReIntegrate();
140 
141  // G4bool final_section= true; // Shows whether current section is last
142  // (i.e. B=full end)
143  G4bool first_section = true;
144  recalculatedEndPoint = false;
145  G4bool restoredFullEndpoint = false;
146 
147  unsigned int substep_no = 0;
148 
149  // Statistics for substeps
150  //
151  static G4ThreadLocal unsigned int max_no_seen= 0;
152 
153  //--------------------------------------------------------------------------
154  // Algorithm for the case if progress in founding intersection is too slow.
155  // Process is defined too slow if after N=param_substeps advances on the
156  // path, it will be only 'fraction_done' of the total length.
157  // In this case the remaining length is divided in two half and
158  // the loop is restarted for each half.
159  // If progress is still too slow, the division in two halfs continue
160  // until 'max_depth'.
161  //--------------------------------------------------------------------------
162 
163  const G4int param_substeps = 5; // Test value for the maximum number
164  // of substeps
165  const G4double fraction_done = 0.3;
166 
167  G4bool Second_half = false; // First half or second half of divided step
168 
169  // We need to know this for the 'final_section':
170  // real 'final_section' or first half 'final_section'
171  // In algorithm it is considered that the 'Second_half' is true
172  // and it becomes false only if we are in the first-half of level
173  // depthness or if we are in the first section
174 
175  unsigned int depth = 0; // Depth counts subdivisions of initial step made
176  ++fNumCalls;
177 
178 #ifdef G4DEBUG_FIELD
179  unsigned int trigger_substepno_print = 0;
180  const G4double tolerance = 1.0e-8 * CLHEP::mm;
181  unsigned int biggest_depth = 0;
182 
183 #if (G4DEBUG_FIELD>1)
184  G4ThreeVector StartPosition = CurveStartPointVelocity.GetPosition();
185  if( (TrialPoint - StartPosition).mag2() < sqr(tolerance))
186  {
187  ReportImmediateHit( MethodName, StartPosition, TrialPoint,
188  tolerance, fNumCalls);
189  }
190 #endif
191 #endif
192 
193  NormalAtEntry = GetSurfaceNormal(CurrentE_Point, validNormalAtE);
194 
195  // Intermediates Points on the Track = Subdivided Points must be stored.
196  // Give the initial values to 'InterMedFt'
197  // Important is 'ptrInterMedFT[0]', it saves the 'EndCurvePoint'
198  //
199  *ptrInterMedFT[0] = CurveEndPointVelocity;
200  for ( auto idepth=1; idepth<max_depth+1; ++idepth )
201  {
202  *ptrInterMedFT[idepth] = CurveStartPointVelocity;
203  }
204 
205  // Final_section boolean store
206  //
207  G4bool fin_section_depth[max_depth];
208  for ( auto idepth=0; idepth<max_depth; ++idepth )
209  {
210  fin_section_depth[idepth] = true;
211  }
212  // 'SubStartPoint' is needed to calculate the length of the divided step
213  //
214  G4FieldTrack SubStart_PointVelocity = CurveStartPointVelocity;
215 
216  do // Loop checking, 07.10.2016, J.Apostolakis
217  {
218  unsigned int substep_no_p = 0;
219  G4bool sub_final_section = false; // the same as final_section,
220  // but for 'sub_section'
221  SubStart_PointVelocity = CurrentA_PointVelocity;
222 
223  do // Loop checking, 07.10.2016, J.Apostolakis
224  { // REPEAT param
225  G4ThreeVector Point_A = CurrentA_PointVelocity.GetPosition();
226  G4ThreeVector Point_B = CurrentB_PointVelocity.GetPosition();
227 
228  // F = a point on true AB path close to point E
229  // (the closest if possible)
230  //
231  ApproxIntersecPointV = GetChordFinderFor()
232  ->ApproxCurvePointV( CurrentA_PointVelocity,
233  CurrentB_PointVelocity,
234  CurrentE_Point,
236  // The above method is the key & most intuitive part ...
237 
238 #ifdef G4DEBUG_FIELD
239  if( ApproxIntersecPointV.GetCurveLength() >
240  CurrentB_PointVelocity.GetCurveLength() * (1.0 + tolerance) )
241  {
242  G4Exception(MethodName, "GeomNav0003", FatalException,
243  "Intermediate F point is past end B point" );
244  }
245 #endif
246 
247  G4ThreeVector CurrentF_Point= ApproxIntersecPointV.GetPosition();
248 
249  // First check whether EF is small - then F is a good approx. point
250  // Calculate the length and direction of the chord AF
251  //
252  G4ThreeVector ChordEF_Vector = CurrentF_Point - CurrentE_Point;
253 
254  G4ThreeVector NewMomentumDir = ApproxIntersecPointV.GetMomentumDir();
255  G4double MomDir_dot_Norm = NewMomentumDir.dot( NormalAtEntry );
256 
257 #ifdef G4DEBUG_FIELD
258  if( fVerboseLevel > 3 )
259  {
260  G4ThreeVector ChordAB = Point_B - Point_A;
261  G4double ABchord_length = ChordAB.mag();
262  G4double MomDir_dot_ABchord;
263  MomDir_dot_ABchord = (1.0 / ABchord_length)
264  * NewMomentumDir.dot( ChordAB );
265  G4VIntersectionLocator::ReportTrialStep( substep_no, ChordAB,
266  ChordEF_Vector, NewMomentumDir, NormalAtEntry, validNormalAtE );
267  G4cout << " dot( MomentumDir, ABchord_unit ) = "
268  << MomDir_dot_ABchord << G4endl;
269  }
270 #endif
271  G4bool adequate_angle =
272  ( MomDir_dot_Norm >= 0.0 ) // Can use ( > -epsilon) instead
273  || (! validNormalAtE ); // Invalid, cannot use this criterion
274  G4double EF_dist2 = ChordEF_Vector.mag2();
275  if ( ( EF_dist2 <= sqr(fiDeltaIntersection) && ( adequate_angle ) )
276  || ( EF_dist2 <= kCarTolerance*kCarTolerance ) )
277  {
278  found_approximate_intersection = true;
279 
280  // Create the "point" return value
281  //
282  IntersectedOrRecalculatedFT = ApproxIntersecPointV;
283  IntersectedOrRecalculatedFT.SetPosition( CurrentE_Point );
284 
286  {
287  // Try to Get Correction of IntersectionPoint using SurfaceNormal()
288  //
289  G4ThreeVector IP;
290  G4ThreeVector MomentumDir=ApproxIntersecPointV.GetMomentumDirection();
291  G4bool goodCorrection = AdjustmentOfFoundIntersection(Point_A,
292  CurrentE_Point, CurrentF_Point, MomentumDir,
293  last_AF_intersection, IP, NewSafety,
294  previousSafety, previousSftOrigin );
295  if ( goodCorrection )
296  {
297  IntersectedOrRecalculatedFT = ApproxIntersecPointV;
298  IntersectedOrRecalculatedFT.SetPosition(IP);
299  }
300  }
301  // Note: in order to return a point on the boundary,
302  // we must return E. But it is F on the curve.
303  // So we must "cheat": we are using the position at point E
304  // and the velocity at point F !!!
305  //
306  // This must limit the length we can allow for displacement!
307  }
308  else // E is NOT close enough to the curve (ie point F)
309  {
310  // Check whether any volumes are encountered by the chord AF
311  // ---------------------------------------------------------
312  // First relocate to restore any Voxel etc information
313  // in the Navigator before calling ComputeStep()
314  //
316 
317  G4ThreeVector PointG; // Candidate intersection point
318  G4double stepLengthAF;
319  G4bool Intersects_AF = IntersectChord( Point_A, CurrentF_Point,
320  NewSafety, previousSafety,
321  previousSftOrigin,
322  stepLengthAF,
323  PointG );
324  last_AF_intersection = Intersects_AF;
325  if( Intersects_AF )
326  {
327  // G is our new Candidate for the intersection point.
328  // It replaces "E" and we will repeat the test to see if
329  // it is a good enough approximate point for us.
330  // B <- F
331  // E <- G
332  //
333  CurrentB_PointVelocity = ApproxIntersecPointV;
334  CurrentE_Point = PointG;
335 
336  validIntersectP = true; // 'E' has been updated.
337 
338  G4bool validNormalLast;
339  NormalAtEntry = GetSurfaceNormal( PointG, validNormalLast );
340  validNormalAtE = validNormalLast;
341 
342  // By moving point B, must take care if current
343  // AF has no intersection to try current FB!!
344  //
345  fin_section_depth[depth] = false;
346 
347 #ifdef G4VERBOSE
348  if( fVerboseLevel > 3 )
349  {
350  G4cout << "G4PiF::LI> Investigating intermediate point"
351  << " at s=" << ApproxIntersecPointV.GetCurveLength()
352  << " on way to full s="
353  << CurveEndPointVelocity.GetCurveLength() << G4endl;
354  }
355 #endif
356  }
357  else // not Intersects_AF
358  {
359  // In this case:
360  // There is NO intersection of AF with a volume boundary.
361  // We must continue the search in the segment FB!
362  //
363  GetNavigatorFor()->LocateGlobalPointWithinVolume( CurrentF_Point );
364 
365  G4double stepLengthFB;
366  G4ThreeVector PointH;
367 
368  // Check whether any volumes are encountered by the chord FB
369  // ---------------------------------------------------------
370 
371  G4bool Intersects_FB = IntersectChord( CurrentF_Point, Point_B,
372  NewSafety, previousSafety,
373  previousSftOrigin,
374  stepLengthFB,
375  PointH );
376  if( Intersects_FB )
377  {
378  // There is an intersection of FB with a volume boundary
379  // H <- First Intersection of Chord FB
380 
381  // H is our new Candidate for the intersection point.
382  // It replaces "E" and we will repeat the test to see if
383  // it is a good enough approximate point for us.
384 
385  // Note that F must be in volume volA (the same as A)
386  // (otherwise AF would meet a volume boundary!)
387  // A <- F
388  // E <- H
389  //
390  CurrentA_PointVelocity = ApproxIntersecPointV;
391  CurrentE_Point = PointH;
392 
393  validIntersectP = true; // 'E' has been updated.
394 
395  G4bool validNormalH;
396  NormalAtEntry = GetSurfaceNormal( PointH, validNormalH );
397  validNormalAtE = validNormalH;
398  }
399  else // not Intersects_FB
400  {
401  if( fin_section_depth[depth] )
402  {
403  // If B is the original endpoint, this means that whatever
404  // volume(s) intersected the original chord, none touch the
405  // smaller chords we have used.
406  // The value of 'IntersectedOrRecalculatedFT' returned is
407  // likely not valid
408 
409  // Check on real final_section or SubEndSection
410  //
411  if( ((Second_half)&&(depth==0)) || (first_section) )
412  {
413  there_is_no_intersection = true; // real final_section
414  }
415  else
416  {
417  // end of subsection, not real final section
418  // exit from the and go to the depth-1 level
419  substep_no_p = param_substeps+2; // exit from the loop
420 
421  // but 'Second_half' is still true because we need to find
422  // the 'CurrentE_point' for the next loop
423  Second_half = true;
424  sub_final_section = true;
425  }
426  }
427  else
428  {
429  CurrentA_PointVelocity = CurrentB_PointVelocity; // Got to B
430  CurrentB_PointVelocity = (depth==0) ? CurveEndPointVelocity
431  : *ptrInterMedFT[depth] ;
432  SubStart_PointVelocity = CurrentA_PointVelocity;
433  restoredFullEndpoint = true;
434 
435  validIntersectP = false; // 'E' has NOT been updated.
436  }
437  } // Endif (Intersects_FB)
438  } // Endif (Intersects_AF)
439 
440  G4int errorEndPt = 0; // Default: no error (if not calling CheckAnd...
441 
442  G4bool recalculatedB= false;
443  if( driverReIntegrates )
444  {
445  G4FieldTrack RevisedB_FT = CurrentB_PointVelocity;
446  recalculatedB= CheckAndReEstimateEndpoint(CurrentA_PointVelocity,
447  CurrentB_PointVelocity,
448  RevisedB_FT,
449  errorEndPt );
450  if( recalculatedB )
451  {
452  CurrentB_PointVelocity = RevisedB_FT; // Use it !
453  // Do not invalidate intersection F -- it is still roughly OK.
454  //
455  // The best course would be to invalidate (reset validIntersectP)
456  // BUT if we invalidate it, we must re-estimate it somewhere! E.g.
457  // validIntersectP = false; // 'E' has NOT been updated.
458 
459  if ( (fin_section_depth[depth]) // real final section
460  &&( first_section || ((Second_half)&&(depth==0)) ) )
461  {
462  recalculatedEndPoint = true;
463  IntersectedOrRecalculatedFT = RevisedB_FT;
464  // So that we can return it, if it is the endpoint!
465  }
466  // else
467  // Move forward the other points
468  // - or better flag it, so that they are re-computed when next used
469  // [ Implementation: a counter for # of recomputations
470  // => avoids extra work]
471  }
472  }
473  else
474  {
475  if( CurrentB_PointVelocity.GetCurveLength() < CurrentA_PointVelocity.GetCurveLength() )
476  errorEndPt = 2;
477  }
478 
479  if( errorEndPt > 1 ) // errorEndPt = 1 is milder, just: len(B)=len(A)
480  {
481  std::ostringstream errmsg;
482  errmsg << "Location: " << MethodName
483  << "- After EndIf(Intersects_AF)" << G4endl;
484  ReportReversedPoints(errmsg,
485  CurveStartPointVelocity, CurveEndPointVelocity,
486  NewSafety, fiEpsilonStep,
487  CurrentA_PointVelocity, CurrentB_PointVelocity,
488  SubStart_PointVelocity, CurrentE_Point,
489  ApproxIntersecPointV, substep_no, substep_no_p, depth);
490  G4Exception(MethodName, "GeomNav0003", FatalException, errmsg);
491  }
492  if( restoredFullEndpoint )
493  {
494  fin_section_depth[depth] = restoredFullEndpoint;
495  restoredFullEndpoint = false;
496  }
497  } // EndIf ( E is close enough to the curve, ie point F. )
498  // tests ChordAF_Vector.mag() <= maximum_lateral_displacement
499 
500 #ifdef G4DEBUG_FIELD
501  if( trigger_substepno_print == 0)
502  {
503  trigger_substepno_print= fWarnSteps - 20;
504  }
505 
506  if( substep_no >= trigger_substepno_print )
507  {
508  G4cout << "Difficulty in converging in " << MethodName
509  << G4endl
510  << " Substep no = " << substep_no << G4endl;
511  if( substep_no == trigger_substepno_print )
512  {
513  printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
514  -1.0, NewSafety, 0 );
515  }
516  G4cout << " State of point A: ";
517  printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
518  -1.0, NewSafety, substep_no-1 );
519  G4cout << " State of point B: ";
520  printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
521  -1.0, NewSafety, substep_no );
522  }
523 #endif
524  ++substep_no;
525  ++substep_no_p;
526 
527  } while ( ( ! found_approximate_intersection )
528  && ( ! there_is_no_intersection )
529  && ( substep_no_p <= param_substeps) ); // UNTIL found or
530  // failed param substep
531 
532  if( (!found_approximate_intersection) && (!there_is_no_intersection) )
533  {
534  G4double did_len = std::abs( CurrentA_PointVelocity.GetCurveLength()
535  - SubStart_PointVelocity.GetCurveLength());
536  G4double all_len = std::abs( CurrentB_PointVelocity.GetCurveLength()
537  - SubStart_PointVelocity.GetCurveLength());
538 
539  G4double distAB = -1;
540  G4ThreeVector PointGe;
541  //
542  // Is progress is too slow, and is it possible to go deeper?
543  // If so, then *halve the step*
544  // ==============
545  if( (did_len < fraction_done*all_len)
546  && (depth<max_depth) && (!sub_final_section) )
547  {
548 #ifdef G4DEBUG_FIELD
549  static G4ThreadLocal unsigned int numSplits=0; // For debugging only
550  biggest_depth = std::max(depth, biggest_depth);
551  ++numSplits;
552 #endif
553  Second_half = false;
554  ++depth;
555  first_section = false;
556 
557  G4double Sub_len = (all_len-did_len)/(2.);
558  G4FieldTrack midPoint = CurrentA_PointVelocity;
559  G4bool fullAdvance=
560  integrDriver->AccurateAdvance(midPoint, Sub_len, fiEpsilonStep);
561 
563  if( fullAdvance ) { ++fNumAdvanceFull; }
564 
565  G4double lenAchieved=
566  midPoint.GetCurveLength()-CurrentA_PointVelocity.GetCurveLength();
567 
568  const G4double adequateFraction = (1.0-CLHEP::perThousand);
569  G4bool goodAdvance = (lenAchieved >= adequateFraction * Sub_len);
570  if ( goodAdvance ) { ++fNumAdvanceGood; }
571 
572 #ifdef G4DEBUG_FIELD
573  else // !goodAdvance
574  {
575  G4cout << "MLL> AdvanceChordLimited not full at depth=" << depth
576  << " total length achieved = " << lenAchieved << " of "
577  << Sub_len << " fraction= ";
578  if (Sub_len != 0.0 ) { G4cout << lenAchieved / Sub_len; }
579  else { G4cout << "DivByZero"; }
580  G4cout << " Good-enough fraction = " << adequateFraction;
581  G4cout << " Number of call to mll = " << fNumCalls
582  << " iteration " << substep_no
583  << " inner = " << substep_no_p << G4endl;
584  G4cout << " Epsilon of integration = " << fiEpsilonStep << G4endl;
585  G4cout << " State at start is = " << CurrentA_PointVelocity
586  << G4endl
587  << " at end (midpoint)= " << midPoint << G4endl;
588  G4cout << " Particle mass = " << midPoint.GetRestMass() << G4endl;
589 
590  G4EquationOfMotion *equation = integrDriver->GetEquationOfMotion();
591  ReportFieldValue( CurrentA_PointVelocity, "start", equation );
592  ReportFieldValue( midPoint, "midPoint", equation );
593  G4cout << " Original Start = "
594  << CurveStartPointVelocity << G4endl;
595  G4cout << " Original End = "
596  << CurveEndPointVelocity << G4endl;
597  G4cout << " Original TrialPoint = "
598  << TrialPoint << G4endl;
599  G4cout << " (this is STRICT mode) "
600  << " num Splits= " << numSplits;
601  G4cout << G4endl;
602  }
603 #endif
604 
605  *ptrInterMedFT[depth] = midPoint;
606  CurrentB_PointVelocity = midPoint;
607 
608  // Adjust 'SubStartPoint' to calculate the 'did_length' in next loop
609  //
610  SubStart_PointVelocity = CurrentA_PointVelocity;
611 
612  // Find new trial intersection point needed at start of the loop
613  //
614  G4ThreeVector Point_A = CurrentA_PointVelocity.GetPosition();
615  G4ThreeVector Point_B = CurrentB_PointVelocity.GetPosition();
616 
618  G4bool Intersects_AB = IntersectChord(Point_A, Point_B,
619  NewSafety, previousSafety,
620  previousSftOrigin, distAB,
621  PointGe);
622  if( Intersects_AB )
623  {
624  last_AF_intersection = Intersects_AB;
625  CurrentE_Point = PointGe;
626  fin_section_depth[depth] = true;
627 
628  validIntersectP = true; // 'E' has been updated.
629 
630  G4bool validNormalAB;
631  NormalAtEntry = GetSurfaceNormal( PointGe, validNormalAB );
632  validNormalAtE = validNormalAB;
633  }
634  else
635  {
636  // No intersection found for first part of curve
637  // (CurrentA,InterMedPoint[depth]). Go to the second part
638  //
639  Second_half = true;
640 
641  validIntersectP= false; // No new 'E' chord intersection found
642  }
643  } // if did_len
644 
645  unsigned int levelPops = 0;
646 
647  G4bool unfinished = Second_half;
648  while ( unfinished && (depth>0) ) // Loop checking, 07.10.2016, JA
649  {
650  // Second part of curve (InterMed[depth],Intermed[depth-1]))
651  // On the depth-1 level normally we are on the 'second_half'
652 
653  ++levelPops;
654 
655  // Find new trial intersection point needed at start of the loop
656  //
657  SubStart_PointVelocity = *ptrInterMedFT[depth];
658  CurrentA_PointVelocity = *ptrInterMedFT[depth];
659  CurrentB_PointVelocity = *ptrInterMedFT[depth-1];
660 
661  G4int errorEndPt = 0; // Default: no error (if not calling CheckAnd...
662 
663  G4bool recalculatedB= false;
664  if( driverReIntegrates )
665  {
666  // Ensure that the new endpoints are not further apart in space
667  // than on the curve due to different errors in the integration
668  //
669  G4FieldTrack RevisedEndPointFT = CurrentB_PointVelocity;
670  recalculatedB =
671  CheckAndReEstimateEndpoint( CurrentA_PointVelocity,
672  CurrentB_PointVelocity,
673  RevisedEndPointFT,
674  errorEndPt );
675  if( recalculatedB )
676  {
677  CurrentB_PointVelocity = RevisedEndPointFT; // Use it !
678 
679  if ( depth == 1 )
680  {
681  recalculatedEndPoint = true;
682  IntersectedOrRecalculatedFT = RevisedEndPointFT;
683  // So that we can return it, if it is the endpoint!
684  }
685  }
686  else
687  {
688  if( CurrentB_PointVelocity.GetCurveLength() < CurrentA_PointVelocity.GetCurveLength() )
689  errorEndPt = 2;
690  }
691  }
692  if( errorEndPt > 1 ) // errorEndPt = 1 is milder, just: len(B)=len(A)
693  {
694  std::ostringstream errmsg;
695  errmsg << "Location: " << MethodName << "- Second-Half" << G4endl;
696  ReportReversedPoints(errmsg,
697  CurveStartPointVelocity, CurveEndPointVelocity,
698  NewSafety, fiEpsilonStep,
699  CurrentA_PointVelocity, CurrentB_PointVelocity,
700  SubStart_PointVelocity, CurrentE_Point,
701  ApproxIntersecPointV, substep_no, substep_no_p, depth);
702  G4Exception(MethodName, "GeomNav0003", FatalException, errmsg);
703  }
704  G4ThreeVector Point_A = CurrentA_PointVelocity.GetPosition();
705  G4ThreeVector Point_B = CurrentB_PointVelocity.GetPosition();
707  G4bool Intersects_AB = IntersectChord(Point_A, Point_B, NewSafety,
708  previousSafety,
709  previousSftOrigin, distAB,
710  PointGe);
711  if( Intersects_AB )
712  {
713  last_AF_intersection = Intersects_AB;
714  CurrentE_Point = PointGe;
715 
716  validIntersectP = true; // 'E' has been updated.
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  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  return (!there_is_no_intersection) && found_approximate_intersection;
826  // Success or failure
827 }
828 
830 {
831  G4cout << " Number of calls = " << fNumCalls << G4endl;
832  G4cout << " Number of split level ('advances'): "
834  G4cout << " Number of full advances: "
835  << fNumAdvanceGood << G4endl;
836  G4cout << " Number of good advances: "
837  << fNumAdvanceFull << G4endl;
838 }
839 
841  const char* nameLoc,
842  const G4EquationOfMotion* equation )
843 {
844  enum { maxNumFieldComp = 24 };
845 
846  G4ThreeVector position = locationPV.GetPosition();
847  G4double startPoint[4] = { position.x(), position.y(), position.z(),
848  locationPV.GetLabTimeOfFlight() };
849  G4double FieldVec[maxNumFieldComp]; // 24 ;
850  for (auto i=0; i<maxNumFieldComp; ++i )
851  {
852  FieldVec[i] = 0.0;
853  }
854  equation->GetFieldValue( startPoint, FieldVec);
855  G4cout << " B-field value (" << nameLoc << ")= "
856  << FieldVec[0] << " " << FieldVec[1] << " " << FieldVec[2];
857  G4double Emag2= G4ThreeVector( FieldVec[3],
858  FieldVec[4],
859  FieldVec[5] ).mag2();
860  if( Emag2 > 0.0 )
861  {
862  G4cout << " Electric = " << FieldVec[3] << " "
863  << FieldVec[4] << " "
864  << FieldVec[5]<< G4endl;
865  }
866  return;
867 }