Geant4  10.00.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 66872 2013-01-15 01:25:57Z japost $
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 
40  : G4VIntersectionLocator(theNavigator)
41 {
42  // In case of too slow progress in finding Intersection Point
43  // intermediates Points on the Track must be stored.
44  // Initialise the array of Pointers [max_depth+1] to do this
45 
46  G4ThreeVector zeroV(0.0,0.0,0.0);
47  for (G4int idepth=0; idepth<max_depth+1; idepth++ )
48  {
49  ptrInterMedFT[ idepth ] = new G4FieldTrack( zeroV, zeroV, 0., 0., 0., 0.);
50  }
51 }
52 
54 {
55  for ( G4int idepth=0; idepth<max_depth+1; idepth++)
56  {
57  delete ptrInterMedFT[idepth];
58  }
59 }
60 
61 // --------------------------------------------------------------------------
62 // G4bool G4PropagatorInField::LocateIntersectionPoint(
63 // const G4FieldTrack& CurveStartPointVelocity, // A
64 // const G4FieldTrack& CurveEndPointVelocity, // B
65 // const G4ThreeVector& TrialPoint, // E
66 // G4FieldTrack& IntersectedOrRecalculated // Output
67 // G4bool& recalculated ) // Out
68 // --------------------------------------------------------------------------
69 //
70 // Function that returns the intersection of the true path with the surface
71 // of the current volume (either the external one or the inner one with one
72 // of the daughters:
73 //
74 // A = Initial point
75 // B = another point
76 //
77 // Both A and B are assumed to be on the true path:
78 //
79 // E is the first point of intersection of the chord AB with
80 // a volume other than A (on the surface of A or of a daughter)
81 //
82 // Convention of Use :
83 // i) If it returns "true", then IntersectionPointVelocity is set
84 // to the approximate intersection point.
85 // ii) If it returns "false", no intersection was found.
86 // The validity of IntersectedOrRecalculated depends on 'recalculated'
87 // a) if latter is false, then IntersectedOrRecalculated is invalid.
88 // b) if latter is true, then IntersectedOrRecalculated is
89 // the new endpoint, due to a re-integration.
90 // --------------------------------------------------------------------------
91 // NOTE: implementation taken from G4PropagatorInField
92 //
94  const G4FieldTrack& CurveStartPointVelocity, // A
95  const G4FieldTrack& CurveEndPointVelocity, // B
96  const G4ThreeVector& TrialPoint, // E
97  G4FieldTrack& IntersectedOrRecalculatedFT, // Output
98  G4bool& recalculatedEndPoint, // Out
99  G4double& previousSafety, // In/Out
100  G4ThreeVector& previousSftOrigin) // In/Out
101 {
102  // Find Intersection Point ( A, B, E ) of true path AB - start at E.
103 
104  G4bool found_approximate_intersection = false;
105  G4bool there_is_no_intersection = false;
106 
107  G4FieldTrack CurrentA_PointVelocity = CurveStartPointVelocity;
108  G4FieldTrack CurrentB_PointVelocity = CurveEndPointVelocity;
109  G4ThreeVector CurrentE_Point = TrialPoint;
110  G4bool validNormalAtE = false;
111  G4ThreeVector NormalAtEntry;
112 
113  G4FieldTrack ApproxIntersecPointV(CurveEndPointVelocity); // FT-Def-Construct
114  G4double NewSafety = 0.0;
115  G4bool last_AF_intersection = false;
116 
117  // G4bool final_section= true; // Shows whether current section is last
118  // (i.e. B=full end)
119  G4bool first_section = true;
120  recalculatedEndPoint = false;
121 
122  G4bool restoredFullEndpoint = false;
123 
124  G4int substep_no = 0;
125 
126  G4int oldprc; // cout/cerr precision settings
127 
128  // Limits for substep number
129  //
130  const G4int max_substeps= 10000; // Test 120 (old value 100 )
131  const G4int warn_substeps= 1000; // 100
132 
133  // Statistics for substeps
134  //
135  static G4ThreadLocal G4int max_no_seen= -1;
136 
137  //--------------------------------------------------------------------------
138  // Algorithm for the case if progress in founding intersection is too slow.
139  // Process is defined too slow if after N=param_substeps advances on the
140  // path, it will be only 'fraction_done' of the total length.
141  // In this case the remaining length is divided in two half and
142  // the loop is restarted for each half.
143  // If progress is still too slow, the division in two halfs continue
144  // until 'max_depth'.
145  //--------------------------------------------------------------------------
146 
147  const G4int param_substeps=5; // Test value for the maximum number
148  // of substeps
149  const G4double fraction_done=0.3;
150 
151  G4bool Second_half = false; // First half or second half of divided step
152 
153  // We need to know this for the 'final_section':
154  // real 'final_section' or first half 'final_section'
155  // In algorithm it is considered that the 'Second_half' is true
156  // and it becomes false only if we are in the first-half of level
157  // depthness or if we are in the first section
158 
159  G4int depth=0; // Depth counts how many subdivisions of initial step made
160 
161 #ifdef G4DEBUG_FIELD
162  static const G4double tolerance = 1.0e-8 * mm;
163  G4ThreeVector StartPosition= CurveStartPointVelocity.GetPosition();
164  if( (TrialPoint - StartPosition).mag() < tolerance)
165  {
166  G4Exception("G4MultiLevelLocator::EstimateIntersectionPoint()",
167  "GeomNav1002", JustWarning,
168  "Intersection point F is exactly at start point A." );
169  }
170 #endif
171 
172  NormalAtEntry = GetSurfaceNormal(CurrentE_Point, validNormalAtE);
173 
174  // Intermediates Points on the Track = Subdivided Points must be stored.
175  // Give the initial values to 'InterMedFt'
176  // Important is 'ptrInterMedFT[0]', it saves the 'EndCurvePoint'
177  //
178  *ptrInterMedFT[0] = CurveEndPointVelocity;
179  for (G4int idepth=1; idepth<max_depth+1; idepth++ )
180  {
181  *ptrInterMedFT[idepth]=CurveStartPointVelocity;
182  }
183 
184  // Final_section boolean store
185  //
186  G4bool fin_section_depth[max_depth];
187  for (G4int idepth=0; idepth<max_depth; idepth++ )
188  {
189  fin_section_depth[idepth]=true;
190  }
191  // 'SubStartPoint' is needed to calculate the length of the divided step
192  //
193  G4FieldTrack SubStart_PointVelocity = CurveStartPointVelocity;
194 
195  do
196  {
197  G4int substep_no_p = 0;
198  G4bool sub_final_section = false; // the same as final_section,
199  // but for 'sub_section'
200  SubStart_PointVelocity = CurrentA_PointVelocity;
201  do // REPEAT param
202  {
203  G4ThreeVector Point_A = CurrentA_PointVelocity.GetPosition();
204  G4ThreeVector Point_B = CurrentB_PointVelocity.GetPosition();
205 
206  // F = a point on true AB path close to point E
207  // (the closest if possible)
208  //
209  ApproxIntersecPointV = GetChordFinderFor()
210  ->ApproxCurvePointV( CurrentA_PointVelocity,
211  CurrentB_PointVelocity,
212  CurrentE_Point,
214  // The above method is the key & most intuitive part ...
215 
216 #ifdef G4DEBUG_FIELD
217  if( ApproxIntersecPointV.GetCurveLength() >
218  CurrentB_PointVelocity.GetCurveLength() * (1.0 + tolerance) )
219  {
220  G4Exception("G4multiLevelLocator::EstimateIntersectionPoint()",
221  "GeomNav0003", FatalException,
222  "Intermediate F point is past end B point" );
223  }
224 #endif
225 
226  G4ThreeVector CurrentF_Point= ApproxIntersecPointV.GetPosition();
227 
228  // First check whether EF is small - then F is a good approx. point
229  // Calculate the length and direction of the chord AF
230  //
231  G4ThreeVector ChordEF_Vector = CurrentF_Point - CurrentE_Point;
232 
233  G4ThreeVector NewMomentumDir= ApproxIntersecPointV.GetMomentumDir();
234  G4double MomDir_dot_Norm= NewMomentumDir.dot( NormalAtEntry ) ;
235 
236 #ifdef G4DEBUG_FIELD
237  if( VerboseLevel > 3 )
238  {
239  G4ThreeVector ChordAB = Point_B - Point_A;
240  G4double ABchord_length = ChordAB.mag();
241  G4double MomDir_dot_ABchord;
242  MomDir_dot_ABchord = (1.0 / ABchord_length)
243  * NewMomentumDir.dot( ChordAB );
244  G4VIntersectionLocator::ReportTrialStep( substep_no, ChordAB,
245  ChordEF_Vector, NewMomentumDir, NormalAtEntry, validNormalAtE );
246  G4cout << " dot( MomentumDir, ABchord_unit ) = "
247  << MomDir_dot_ABchord << G4endl;
248  }
249 #endif
250  G4bool adequate_angle =
251  ( MomDir_dot_Norm >= 0.0 ) // Can use ( > -epsilon) instead
252  || (! validNormalAtE ); // Invalid, cannot use this criterion
253  G4double EF_dist2 = ChordEF_Vector.mag2();
254  if ( ( EF_dist2 <= sqr(fiDeltaIntersection) && ( adequate_angle ) )
255  || ( EF_dist2 <= kCarTolerance*kCarTolerance ) )
256  {
257  found_approximate_intersection = true;
258 
259  // Create the "point" return value
260  //
261  IntersectedOrRecalculatedFT = ApproxIntersecPointV;
262  IntersectedOrRecalculatedFT.SetPosition( CurrentE_Point );
263 
265  {
266  // Try to Get Correction of IntersectionPoint using SurfaceNormal()
267  //
268  G4ThreeVector IP;
269  G4ThreeVector MomentumDir=ApproxIntersecPointV.GetMomentumDirection();
270  G4bool goodCorrection = AdjustmentOfFoundIntersection(Point_A,
271  CurrentE_Point, CurrentF_Point, MomentumDir,
272  last_AF_intersection, IP, NewSafety,
273  previousSafety, previousSftOrigin );
274  if ( goodCorrection )
275  {
276  IntersectedOrRecalculatedFT = ApproxIntersecPointV;
277  IntersectedOrRecalculatedFT.SetPosition(IP);
278  }
279  }
280  // Note: in order to return a point on the boundary,
281  // we must return E. But it is F on the curve.
282  // So we must "cheat": we are using the position at point E
283  // and the velocity at point F !!!
284  //
285  // This must limit the length we can allow for displacement!
286  }
287  else // E is NOT close enough to the curve (ie point F)
288  {
289  // Check whether any volumes are encountered by the chord AF
290  // ---------------------------------------------------------
291  // First relocate to restore any Voxel etc information
292  // in the Navigator before calling ComputeStep()
293  //
295 
296  G4ThreeVector PointG; // Candidate intersection point
297  G4double stepLengthAF;
298  G4bool Intersects_AF = IntersectChord( Point_A, CurrentF_Point,
299  NewSafety, previousSafety,
300  previousSftOrigin,
301  stepLengthAF,
302  PointG );
303  last_AF_intersection = Intersects_AF;
304  if( Intersects_AF )
305  {
306  // G is our new Candidate for the intersection point.
307  // It replaces "E" and we will repeat the test to see if
308  // it is a good enough approximate point for us.
309  // B <- F
310  // E <- G
311  //
312  CurrentB_PointVelocity = ApproxIntersecPointV;
313  CurrentE_Point = PointG;
314 
315  G4bool validNormalLast;
316  NormalAtEntry = GetSurfaceNormal( PointG, validNormalLast );
317  validNormalAtE = validNormalLast;
318 
319  // By moving point B, must take care if current
320  // AF has no intersection to try current FB!!
321  //
322  fin_section_depth[depth]=false;
323 
324 
325 #ifdef G4VERBOSE
326  if( fVerboseLevel > 3 )
327  {
328  G4cout << "G4PiF::LI> Investigating intermediate point"
329  << " at s=" << ApproxIntersecPointV.GetCurveLength()
330  << " on way to full s="
331  << CurveEndPointVelocity.GetCurveLength() << G4endl;
332  }
333 #endif
334  }
335  else // not Intersects_AF
336  {
337  // In this case:
338  // There is NO intersection of AF with a volume boundary.
339  // We must continue the search in the segment FB!
340  //
341  GetNavigatorFor()->LocateGlobalPointWithinVolume( CurrentF_Point );
342 
343  G4double stepLengthFB;
344  G4ThreeVector PointH;
345 
346  // Check whether any volumes are encountered by the chord FB
347  // ---------------------------------------------------------
348 
349  G4bool Intersects_FB = IntersectChord( CurrentF_Point, Point_B,
350  NewSafety, previousSafety,
351  previousSftOrigin,
352  stepLengthFB,
353  PointH );
354  if( Intersects_FB )
355  {
356  // There is an intersection of FB with a volume boundary
357  // H <- First Intersection of Chord FB
358 
359  // H is our new Candidate for the intersection point.
360  // It replaces "E" and we will repeat the test to see if
361  // it is a good enough approximate point for us.
362 
363  // Note that F must be in volume volA (the same as A)
364  // (otherwise AF would meet a volume boundary!)
365  // A <- F
366  // E <- H
367  //
368  CurrentA_PointVelocity = ApproxIntersecPointV;
369  CurrentE_Point = PointH;
370 
371  G4bool validNormalH;
372  NormalAtEntry = GetSurfaceNormal( PointH, validNormalH );
373  validNormalAtE = validNormalH;
374 
375  }
376  else // not Intersects_FB
377  {
378  // There is NO intersection of FB with a volume boundary
379 
380  if(fin_section_depth[depth])
381  {
382  // If B is the original endpoint, this means that whatever
383  // volume(s) intersected the original chord, none touch the
384  // smaller chords we have used.
385  // The value of 'IntersectedOrRecalculatedFT' returned is
386  // likely not valid
387 
388  // Check on real final_section or SubEndSection
389  //
390  if( ((Second_half)&&(depth==0)) || (first_section) )
391  {
392  there_is_no_intersection = true; // real final_section
393  }
394  else
395  {
396  // end of subsection, not real final section
397  // exit from the and go to the depth-1 level
398 
399  substep_no_p = param_substeps+2; // exit from the loop
400 
401  // but 'Second_half' is still true because we need to find
402  // the 'CurrentE_point' for the next loop
403  //
404  Second_half = true;
405  sub_final_section = true;
406 
407  }
408  }
409  else
410  {
411  if(depth==0)
412  {
413  // We must restore the original endpoint
414  //
415  CurrentA_PointVelocity = CurrentB_PointVelocity; // Got to B
416  CurrentB_PointVelocity = CurveEndPointVelocity;
417  SubStart_PointVelocity = CurrentA_PointVelocity;
418  restoredFullEndpoint = true;
419  }
420  else
421  {
422  // We must restore the depth endpoint
423  //
424  CurrentA_PointVelocity = CurrentB_PointVelocity; // Got to B
425  CurrentB_PointVelocity = *ptrInterMedFT[depth];
426  SubStart_PointVelocity = CurrentA_PointVelocity;
427  restoredFullEndpoint = true;
428  }
429  }
430  } // Endif (Intersects_FB)
431  } // Endif (Intersects_AF)
432 
433  // Ensure that the new endpoints are not further apart in space
434  // than on the curve due to different errors in the integration
435  //
436  G4double linDistSq, curveDist;
437  linDistSq = ( CurrentB_PointVelocity.GetPosition()
438  - CurrentA_PointVelocity.GetPosition() ).mag2();
439  curveDist = CurrentB_PointVelocity.GetCurveLength()
440  - CurrentA_PointVelocity.GetCurveLength();
441 
442  // Change this condition for very strict parameters of propagation
443  //
444  if( curveDist*curveDist*(1+2* GetEpsilonStepFor()) < linDistSq )
445  {
446  // Re-integrate to obtain a new B
447  //
448  G4FieldTrack newEndPointFT=
449  ReEstimateEndpoint( CurrentA_PointVelocity,
450  CurrentB_PointVelocity,
451  linDistSq, // to avoid recalculation
452  curveDist );
453  G4FieldTrack oldPointVelB = CurrentB_PointVelocity;
454  CurrentB_PointVelocity = newEndPointFT;
455 
456  if ( (fin_section_depth[depth]) // real final section
457  &&( first_section || ((Second_half)&&(depth==0)) ) )
458  {
459  recalculatedEndPoint = true;
460  IntersectedOrRecalculatedFT = newEndPointFT;
461  // So that we can return it, if it is the endpoint!
462  }
463  }
464  if( curveDist < 0.0 )
465  {
466  fVerboseLevel = 5; // Print out a maximum of information
467  printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
468  -1.0, NewSafety, substep_no );
469  std::ostringstream message;
470  message << "Error in advancing propagation." << G4endl
471  << " Point A (start) is " << CurrentA_PointVelocity
472  << G4endl
473  << " Point B (end) is " << CurrentB_PointVelocity
474  << G4endl
475  << " Curve distance is " << curveDist << G4endl
476  << G4endl
477  << "The final curve point is not further along"
478  << " than the original!" << G4endl;
479 
480  if( recalculatedEndPoint )
481  {
482  message << "Recalculation of EndPoint was called with fEpsStep= "
483  << GetEpsilonStepFor() << G4endl;
484  }
485  oldprc = G4cerr.precision(20);
486  message << " Point A (Curve start) is " << CurveStartPointVelocity
487  << G4endl
488  << " Point B (Curve end) is " << CurveEndPointVelocity
489  << G4endl
490  << " Point A (Current start) is " << CurrentA_PointVelocity
491  << G4endl
492  << " Point B (Current end) is " << CurrentB_PointVelocity
493  << G4endl
494  << " Point S (Sub start) is " << SubStart_PointVelocity
495  << G4endl
496  << " Point E (Trial Point) is " << CurrentE_Point
497  << G4endl
498  << " Point F (Intersection) is " << ApproxIntersecPointV
499  << G4endl
500  << " LocateIntersection parameters are : Substep no= "
501  << substep_no << G4endl
502  << " Substep depth no= "<< substep_no_p << " Depth= "
503  << depth;
504  G4cerr.precision(oldprc);
505 
506  G4Exception("G4MultiLevelLocator::EstimateIntersectionPoint()",
507  "GeomNav0003", FatalException, message);
508  }
509  if(restoredFullEndpoint)
510  {
511  fin_section_depth[depth] = restoredFullEndpoint;
512  restoredFullEndpoint = false;
513  }
514  } // EndIf ( E is close enough to the curve, ie point F. )
515  // tests ChordAF_Vector.mag() <= maximum_lateral_displacement
516 
517 #ifdef G4DEBUG_LOCATE_INTERSECTION
518  static G4int trigger_substepno_print= warn_substeps - 20 ;
519 
520  if( substep_no >= trigger_substepno_print )
521  {
522  G4cout << "Difficulty in converging in "
523  << "G4MultiLevelLocator::EstimateIntersectionPoint():"
524  << G4endl
525  << " Substep no = " << substep_no << G4endl;
526  if( substep_no == trigger_substepno_print )
527  {
528  printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
529  -1.0, NewSafety, 0);
530  }
531  G4cout << " State of point A: ";
532  printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
533  -1.0, NewSafety, substep_no-1, 0);
534  G4cout << " State of point B: ";
535  printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
536  -1.0, NewSafety, substep_no);
537  }
538 #endif
539  substep_no++;
540  substep_no_p++;
541 
542  } while ( ( ! found_approximate_intersection )
543  && ( ! there_is_no_intersection )
544  && ( substep_no_p <= param_substeps) ); // UNTIL found or
545  // failed param substep
546  first_section = false;
547 
548  if( (!found_approximate_intersection) && (!there_is_no_intersection) )
549  {
550  G4double did_len = std::abs( CurrentA_PointVelocity.GetCurveLength()
551  - SubStart_PointVelocity.GetCurveLength());
552  G4double all_len = std::abs( CurrentB_PointVelocity.GetCurveLength()
553  - SubStart_PointVelocity.GetCurveLength());
554 
555  G4double stepLengthAB;
556  G4ThreeVector PointGe;
557  // Check if progress is too slow and if it possible to go deeper,
558  // then halve the step if so
559  //
560  if( ( ( did_len )<fraction_done*all_len)
561  && (depth<max_depth) && (!sub_final_section) )
562  {
563  Second_half=false;
564  depth++;
565 
566  G4double Sub_len = (all_len-did_len)/(2.);
567  G4FieldTrack start = CurrentA_PointVelocity;
568  G4MagInt_Driver* integrDriver
570  integrDriver->AccurateAdvance(start, Sub_len, GetEpsilonStepFor());
571  *ptrInterMedFT[depth] = start;
572  CurrentB_PointVelocity = *ptrInterMedFT[depth];
573 
574  // Adjust 'SubStartPoint' to calculate the 'did_length' in next loop
575  //
576  SubStart_PointVelocity = CurrentA_PointVelocity;
577 
578  // Find new trial intersection point needed at start of the loop
579  //
580  G4ThreeVector Point_A = CurrentA_PointVelocity.GetPosition();
581  G4ThreeVector SubE_point = CurrentB_PointVelocity.GetPosition();
582 
584  G4bool Intersects_AB = IntersectChord(Point_A, SubE_point,
585  NewSafety, previousSafety,
586  previousSftOrigin, stepLengthAB,
587  PointGe);
588  if( Intersects_AB )
589  {
590  last_AF_intersection = Intersects_AB;
591  CurrentE_Point = PointGe;
592  fin_section_depth[depth]=true;
593 
594  G4bool validNormalAB;
595  NormalAtEntry = GetSurfaceNormal( PointGe, validNormalAB );
596  validNormalAtE = validNormalAB;
597  }
598  else
599  {
600  // No intersection found for first part of curve
601  // (CurrentA,InterMedPoint[depth]). Go to the second part
602  //
603  Second_half = true;
604  }
605  } // if did_len
606 
607  if( (Second_half)&&(depth!=0) )
608  {
609  // Second part of curve (InterMed[depth],Intermed[depth-1]))
610  // On the depth-1 level normally we are on the 'second_half'
611 
612  Second_half = true;
613  // Find new trial intersection point needed at start of the loop
614  //
615  SubStart_PointVelocity = *ptrInterMedFT[depth];
616  CurrentA_PointVelocity = *ptrInterMedFT[depth];
617  CurrentB_PointVelocity = *ptrInterMedFT[depth-1];
618  // Ensure that the new endpoints are not further apart in space
619  // than on the curve due to different errors in the integration
620  //
621  G4double linDistSq, curveDist;
622  linDistSq = ( CurrentB_PointVelocity.GetPosition()
623  - CurrentA_PointVelocity.GetPosition() ).mag2();
624  curveDist = CurrentB_PointVelocity.GetCurveLength()
625  - CurrentA_PointVelocity.GetCurveLength();
626  if( curveDist*curveDist*(1+2*GetEpsilonStepFor() ) < linDistSq )
627  {
628  // Re-integrate to obtain a new B
629  //
630  G4FieldTrack newEndPointFT=
631  ReEstimateEndpoint( CurrentA_PointVelocity,
632  CurrentB_PointVelocity,
633  linDistSq, // to avoid recalculation
634  curveDist );
635  G4FieldTrack oldPointVelB = CurrentB_PointVelocity;
636  CurrentB_PointVelocity = newEndPointFT;
637  if (depth==1)
638  {
639  recalculatedEndPoint = true;
640  IntersectedOrRecalculatedFT = newEndPointFT;
641  // So that we can return it, if it is the endpoint!
642  }
643  }
644 
645  G4ThreeVector Point_A = CurrentA_PointVelocity.GetPosition();
646  G4ThreeVector SubE_point = CurrentB_PointVelocity.GetPosition();
648  G4bool Intersects_AB = IntersectChord(Point_A, SubE_point, NewSafety,
649  previousSafety,
650  previousSftOrigin, stepLengthAB,
651  PointGe);
652  if( Intersects_AB )
653  {
654  last_AF_intersection = Intersects_AB;
655  CurrentE_Point = PointGe;
656 
657  G4bool validNormalAB;
658  NormalAtEntry = GetSurfaceNormal( PointGe, validNormalAB );
659  validNormalAtE = validNormalAB;
660  }
661  depth--;
662  fin_section_depth[depth]=true;
663  }
664  } // if(!found_aproximate_intersection)
665 
666  } while ( ( ! found_approximate_intersection )
667  && ( ! there_is_no_intersection )
668  && ( substep_no <= max_substeps) ); // UNTIL found or failed
669 
670  if( substep_no > max_no_seen )
671  {
672  max_no_seen = substep_no;
673 #ifdef G4DEBUG_LOCATE_INTERSECTION
674  if( max_no_seen > warn_substeps )
675  {
676  trigger_substepno_print = max_no_seen-20; // Want to see last 20 steps
677  }
678 #endif
679  }
680 
681  if( ( substep_no >= max_substeps)
682  && !there_is_no_intersection
683  && !found_approximate_intersection )
684  {
685  G4cout << "ERROR - G4MultiLevelLocator::EstimateIntersectionPoint()"
686  << G4endl
687  << " Start and end-point of requested Step:" << G4endl;
688  printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
689  -1.0, NewSafety, 0);
690  G4cout << " Start and end-point of current Sub-Step:" << G4endl;
691  printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
692  -1.0, NewSafety, substep_no-1);
693  printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
694  -1.0, NewSafety, substep_no);
695  std::ostringstream message;
696  message << "Too many substeps!" << G4endl
697  << " Convergence is requiring too many substeps: "
698  << substep_no << G4endl
699  << " Abandoning effort to intersect. " << G4endl
700  << " Found intersection = "
701  << found_approximate_intersection << G4endl
702  << " Intersection exists = "
703  << !there_is_no_intersection << G4endl;
704 
705 #ifdef FUTURE_CORRECTION
706  // Attempt to correct the results of the method // FIX - TODO
707 
708  if ( ! found_approximate_intersection )
709  {
710  recalculatedEndPoint = true;
711  // Return the further valid intersection point -- potentially A ??
712  // JA/19 Jan 2006
713  IntersectedOrRecalculatedFT = CurrentA_PointVelocity;
714 
715  message << " Did not convergence after " << substep_no
716  << " substeps." << G4endl
717  << " The endpoint was adjused to pointA resulting"
718  << G4endl
719  << " from the last substep: " << CurrentA_PointVelocity
720  << G4endl;
721  }
722 #endif
723 
724  oldprc = G4cout.precision( 10 );
725  G4double done_len = CurrentA_PointVelocity.GetCurveLength();
726  G4double full_len = CurveEndPointVelocity.GetCurveLength();
727  message << " Undertaken only length: " << done_len
728  << " out of " << full_len << " required." << G4endl
729  << " Remaining length = " << full_len - done_len;
730  G4cout.precision( oldprc );
731 
732  G4Exception("G4MultiLevelLocator::EstimateIntersectionPoint()",
733  "GeomNav0003", FatalException, message);
734  }
735  else if( substep_no >= warn_substeps )
736  {
737  oldprc = G4cout.precision( 10 );
738  std::ostringstream message;
739  message << "Many substeps while trying to locate intersection."
740  << G4endl
741  << " Undertaken length: "
742  << CurrentB_PointVelocity.GetCurveLength()
743  << " - Needed: " << substep_no << " substeps." << G4endl
744  << " Warning level = " << warn_substeps
745  << " and maximum substeps = " << max_substeps;
746  G4Exception("G4MultiLevelLocator::EstimateIntersectionPoint()",
747  "GeomNav1002", JustWarning, message);
748  G4cout.precision( oldprc );
749  }
750  return !there_is_no_intersection; // Success or failure
751 }
G4FieldTrack * ptrInterMedFT[max_depth+1]
G4FieldTrack ReEstimateEndpoint(const G4FieldTrack &CurrentStateA, const G4FieldTrack &EstimtdEndStateB, G4double linearDistSq, G4double curveDist)
void SetPosition(G4ThreeVector nPos)
G4double GetCurveLength() const
CLHEP::Hep3Vector G4ThreeVector
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)
#define G4ThreadLocal
Definition: tls.hh:52
int G4int
Definition: G4Types.hh:78
G4double GetEpsilonStepFor()
G4ThreeVector GetPosition() const
G4GLOB_DLL std::ostream G4cout
static const G4int max_depth
G4MultiLevelLocator(G4Navigator *theNavigator)
bool G4bool
Definition: G4Types.hh:79
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 printStatus(const G4FieldTrack &startFT, const G4FieldTrack &currentFT, G4double requestStep, G4double safety, G4int step)
G4bool GetAdjustementOfFoundIntersection()
#define G4endl
Definition: G4ios.hh:61
T sqr(const T &x)
Definition: templates.hh:145
double G4double
Definition: G4Types.hh:76
G4MagInt_Driver * GetIntegrationDriver()
static const double mm
Definition: G4SIunits.hh:102
G4ChordFinder * GetChordFinderFor()
G4bool AccurateAdvance(G4FieldTrack &y_current, G4double hstep, G4double eps, G4double hinitial=0.0)
virtual void LocateGlobalPointWithinVolume(const G4ThreeVector &position)
Definition: G4Navigator.cc:550
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)
G4GLOB_DLL std::ostream G4cerr
G4ThreeVector GetMomentumDirection() const