source: trunk/source/geometry/navigation/src/G4BrentLocator.cc@ 1036

Last change on this file since 1036 was 985, checked in by garnier, 17 years ago

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