Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4MultiLevelLocator Class Reference

#include <G4MultiLevelLocator.hh>

+ Inheritance diagram for G4MultiLevelLocator:

Public Member Functions

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

Additional Inherited Members

- Protected Member Functions inherited from G4VIntersectionLocator
G4FieldTrack ReEstimateEndpoint (const G4FieldTrack &CurrentStateA, const G4FieldTrack &EstimtdEndStateB, G4double linearDistSq, G4double curveDist)
 
G4ThreeVector GetSurfaceNormal (const G4ThreeVector &CurrentInt_Point, G4bool &validNormal)
 
G4ThreeVector GetGlobalSurfaceNormal (const G4ThreeVector &CurrentE_Point, G4bool &validNormal)
 
G4bool AdjustmentOfFoundIntersection (const G4ThreeVector &A, const G4ThreeVector &CurrentE_Point, const G4ThreeVector &CurrentF_Point, const G4ThreeVector &MomentumDir, const G4bool IntersectAF, G4ThreeVector &IntersectionPoint, G4double &NewSafety, G4double &fPrevSafety, G4ThreeVector &fPrevSftOrigin)
 
void ReportTrialStep (G4int step_no, const G4ThreeVector &ChordAB_v, const G4ThreeVector &ChordEF_v, const G4ThreeVector &NewMomentumDir, const G4ThreeVector &NormalAtEntry, G4bool validNormal)
 
- Protected Attributes inherited from G4VIntersectionLocator
G4double kCarTolerance
 
G4int fVerboseLevel
 
G4bool fUseNormalCorrection
 
G4NavigatorfiNavigator
 
G4ChordFinderfiChordFinder
 
G4double fiEpsilonStep
 
G4double fiDeltaIntersection
 
G4bool fiUseSafety
 
G4NavigatorfHelpingNavigator
 
G4TouchableHistoryfpTouchable
 

Detailed Description

Definition at line 51 of file G4MultiLevelLocator.hh.

Constructor & Destructor Documentation

◆ G4MultiLevelLocator()

G4MultiLevelLocator::G4MultiLevelLocator ( G4Navigator theNavigator)

Definition at line 39 of file G4MultiLevelLocator.cc.

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}
int G4int
Definition: G4Types.hh:66

◆ ~G4MultiLevelLocator()

G4MultiLevelLocator::~G4MultiLevelLocator ( )

Definition at line 53 of file G4MultiLevelLocator.cc.

54{
55 for ( G4int idepth=0; idepth<max_depth+1; idepth++)
56 {
57 delete ptrInterMedFT[idepth];
58 }
59}

Member Function Documentation

◆ EstimateIntersectionPoint()

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

Implements G4VIntersectionLocator.

Definition at line 93 of file G4MultiLevelLocator.cc.

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 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 //
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= "
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}
@ JustWarning
@ FatalException
double G4double
Definition: G4Types.hh:64
bool G4bool
Definition: G4Types.hh:67
#define G4endl
Definition: G4ios.hh:52
G4DLLIMPORT std::ostream G4cerr
G4DLLIMPORT std::ostream G4cout
double mag2() const
double dot(const Hep3Vector &) const
double mag() const
G4FieldTrack ApproxCurvePointV(const G4FieldTrack &curveAPointVelocity, const G4FieldTrack &curveBPointVelocity, const G4ThreeVector &currentEPoint, G4double epsStep)
G4MagInt_Driver * GetIntegrationDriver()
G4double GetCurveLength() const
G4ThreeVector GetPosition() const
G4bool AccurateAdvance(G4FieldTrack &y_current, G4double hstep, G4double eps, G4double hinitial=0.0)
virtual void LocateGlobalPointWithinVolume(const G4ThreeVector &position)
Definition: G4Navigator.cc:548
G4Navigator * GetNavigatorFor()
G4ChordFinder * GetChordFinderFor()
G4ThreeVector GetSurfaceNormal(const G4ThreeVector &CurrentInt_Point, G4bool &validNormal)
void ReportTrialStep(G4int step_no, const G4ThreeVector &ChordAB_v, const G4ThreeVector &ChordEF_v, const G4ThreeVector &NewMomentumDir, const G4ThreeVector &NormalAtEntry, G4bool validNormal)
G4bool IntersectChord(const G4ThreeVector &StartPointA, const G4ThreeVector &EndPointB, G4double &NewSafety, G4double &PreviousSafety, G4ThreeVector &PreviousSftOrigin, G4double &LinearStepLength, G4ThreeVector &IntersectionPoint, G4bool *calledNavigator=0)
G4double GetEpsilonStepFor()
G4FieldTrack ReEstimateEndpoint(const G4FieldTrack &CurrentStateA, const G4FieldTrack &EstimtdEndStateB, G4double linearDistSq, G4double curveDist)
G4bool GetAdjustementOfFoundIntersection()
G4bool AdjustmentOfFoundIntersection(const G4ThreeVector &A, const G4ThreeVector &CurrentE_Point, const G4ThreeVector &CurrentF_Point, const G4ThreeVector &MomentumDir, const G4bool IntersectAF, G4ThreeVector &IntersectionPoint, G4double &NewSafety, G4double &fPrevSafety, G4ThreeVector &fPrevSftOrigin)
void printStatus(const G4FieldTrack &startFT, const G4FieldTrack &currentFT, G4double requestStep, G4double safety, G4int step)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
T sqr(const T &x)
Definition: templates.hh:145

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