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

Last change on this file since 1355 was 1347, checked in by garnier, 15 years ago

geant4 tag 9.4
File size: 31.9 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.9 2010/07/13 15:59:42 gcosmo Exp $
27// GEANT4 tag $Name: geant4-09-04-ref-00 $
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 oldprc; // cout, cerr precision
145 G4int substep_no = 0;
146
147 // Limits for substep number
148 //
149 const G4int max_substeps= 10000; // Test 120 (old value 100 )
150 const G4int warn_substeps= 1000; // 100
151
152 // Statistics for substeps
153 //
154 static G4int max_no_seen= -1;
155 static G4int trigger_substepno_print= warn_substeps - 20 ;
156
157 // Counter for restarting Bintermed
158 //
159 G4int restartB = 0;
160
161 //--------------------------------------------------------------------------
162 // Algorithm for the case if progress in founding intersection is too slow.
163 // Process is defined too slow if after N=param_substeps advances on the
164 // path, it will be only 'fraction_done' of the total length.
165 // In this case the remaining length is divided in two half and
166 // the loop is restarted for each half.
167 // If progress is still too slow, the division in two halfs continue
168 // until 'max_depth'.
169 //--------------------------------------------------------------------------
170
171 const G4int param_substeps=50; // Test value for the maximum number
172 // of substeps
173 const G4double fraction_done=0.3;
174
175 G4bool Second_half = false; // First half or second half of divided step
176
177 // We need to know this for the 'final_section':
178 // real 'final_section' or first half 'final_section'
179 // In algorithm it is considered that the 'Second_half' is true
180 // and it becomes false only if we are in the first-half of level
181 // depthness or if we are in the first section
182
183 G4int depth=0; // Depth counts how many subdivisions of initial step made
184
185#ifdef G4DEBUG_FIELD
186 static G4double tolerance= 1.0e-8;
187 G4ThreeVector StartPosition= CurveStartPointVelocity.GetPosition();
188 if( (TrialPoint - StartPosition).mag() < tolerance * mm )
189 {
190 G4cerr << "WARNING - G4BrentLocator::EstimateIntersectionPoint()"
191 << G4endl
192 << " Intermediate F point is on top of starting point A."
193 << G4endl;
194 G4Exception("G4BrentLocator::EstimateIntersectionPoint()",
195 "IntersectionPointIsAtStart", JustWarning,
196 "Intersection point F is exactly at start point A." );
197 }
198#endif
199
200 // Intermediates Points on the Track = Subdivided Points must be stored.
201 // Give the initial values to 'InterMedFt'
202 // Important is 'ptrInterMedFT[0]', it saves the 'EndCurvePoint'
203 //
204 *ptrInterMedFT[0] = CurveEndPointVelocity;
205 for (G4int idepth=1; idepth<max_depth+1; idepth++ )
206 {
207 *ptrInterMedFT[idepth]=CurveStartPointVelocity;
208 }
209
210 //Final_section boolean store
211 G4bool fin_section_depth[max_depth];
212 for (G4int idepth=0; idepth<max_depth; idepth++ )
213 {
214 fin_section_depth[idepth]=true;
215 }
216
217 // 'SubStartPoint' is needed to calculate the length of the divided step
218 //
219 G4FieldTrack SubStart_PointVelocity = CurveStartPointVelocity;
220
221 do
222 {
223 G4int substep_no_p = 0;
224 G4bool sub_final_section = false; // the same as final_section,
225 // but for 'sub_section'
226 SubStart_PointVelocity = CurrentA_PointVelocity;
227 do // REPEAT param
228 {
229 G4ThreeVector Point_A = CurrentA_PointVelocity.GetPosition();
230 G4ThreeVector Point_B = CurrentB_PointVelocity.GetPosition();
231
232 // F = a point on true AB path close to point E
233 // (the closest if possible)
234 //
235 if(substep_no_p==0)
236 {
237 ApproxIntersecPointV = GetChordFinderFor()
238 ->ApproxCurvePointV( CurrentA_PointVelocity,
239 CurrentB_PointVelocity,
240 CurrentE_Point,
241 GetEpsilonStepFor());
242 // The above method is the key & most intuitive part ...
243 }
244#ifdef G4DEBUG_FIELD
245 if( ApproxIntersecPointV.GetCurveLength() >
246 CurrentB_PointVelocity.GetCurveLength() * (1.0 + tolerance) )
247 {
248 G4cerr << "ERROR - G4BrentLocator::EstimateIntersectionPoint()"
249 << G4endl
250 << " Intermediate F point is more advanced than"
251 << " endpoint B." << G4endl;
252 G4Exception("G4BrentLocator::EstimateIntersectionPoint()",
253 "IntermediatePointConfusion", FatalException,
254 "Intermediate F point is past end B point" );
255 }
256#endif
257
258 G4ThreeVector CurrentF_Point= ApproxIntersecPointV.GetPosition();
259
260 // First check whether EF is small - then F is a good approx. point
261 // Calculate the length and direction of the chord AF
262 //
263 G4ThreeVector ChordEF_Vector = CurrentF_Point - CurrentE_Point;
264
265 if ( ChordEF_Vector.mag2() <= sqr(GetDeltaIntersectionFor()) )
266 {
267 found_approximate_intersection = true;
268
269 // Create the "point" return value
270 //
271 IntersectedOrRecalculatedFT = ApproxIntersecPointV;
272 IntersectedOrRecalculatedFT.SetPosition( CurrentE_Point );
273
274 if ( GetAdjustementOfFoundIntersection() )
275 {
276 // Try to Get Correction of IntersectionPoint using SurfaceNormal()
277 //
278 G4ThreeVector IP;
279 G4ThreeVector MomentumDir=ApproxIntersecPointV.GetMomentumDirection();
280 G4bool goodCorrection = AdjustmentOfFoundIntersection( Point_A,
281 CurrentE_Point, CurrentF_Point, MomentumDir,
282 last_AF_intersection, IP, NewSafety,
283 fPreviousSafety, fPreviousSftOrigin );
284 if ( goodCorrection )
285 {
286 IntersectedOrRecalculatedFT = ApproxIntersecPointV;
287 IntersectedOrRecalculatedFT.SetPosition(IP);
288 }
289 }
290
291 // Note: in order to return a point on the boundary,
292 // we must return E. But it is F on the curve.
293 // So we must "cheat": we are using the position at point E
294 // and the velocity at point F !!!
295 //
296 // This must limit the length we can allow for displacement!
297 }
298 else // E is NOT close enough to the curve (ie point F)
299 {
300 // Check whether any volumes are encountered by the chord AF
301 // ---------------------------------------------------------
302 // First relocate to restore any Voxel etc information
303 // in the Navigator before calling ComputeStep()
304 //
305 GetNavigatorFor()->LocateGlobalPointWithinVolume( Point_A );
306
307 G4ThreeVector PointG; // Candidate intersection point
308 G4double stepLengthAF;
309 G4bool Intersects_AF = IntersectChord( Point_A, CurrentF_Point,
310 NewSafety,fPreviousSafety,
311 fPreviousSftOrigin,
312 stepLengthAF,
313 PointG );
314 last_AF_intersection = Intersects_AF;
315 if( Intersects_AF )
316 {
317 // G is our new Candidate for the intersection point.
318 // It replaces "E" and we will repeat the test to see if
319 // it is a good enough approximate point for us.
320 // B <- F
321 // E <- G
322 //
323 G4FieldTrack EndPoint = ApproxIntersecPointV;
324 ApproxIntersecPointV = GetChordFinderFor()->ApproxCurvePointS(
325 CurrentA_PointVelocity, CurrentB_PointVelocity,
326 EndPoint,CurrentE_Point, CurrentF_Point,PointG,
327 true, GetEpsilonStepFor() );
328 CurrentB_PointVelocity = EndPoint;
329 CurrentE_Point = PointG;
330
331 // By moving point B, must take care if current
332 // AF has no intersection to try current FB!!
333 //
334 fin_section_depth[depth] = false;
335#ifdef G4VERBOSE
336 if( fVerboseLevel > 3 )
337 {
338 G4cout << "G4PiF::LI> Investigating intermediate point"
339 << " at s=" << ApproxIntersecPointV.GetCurveLength()
340 << " on way to full s="
341 << CurveEndPointVelocity.GetCurveLength() << G4endl;
342 }
343#endif
344 }
345 else // not Intersects_AF
346 {
347 // In this case:
348 // There is NO intersection of AF with a volume boundary.
349 // We must continue the search in the segment FB!
350 //
351 GetNavigatorFor()->LocateGlobalPointWithinVolume( CurrentF_Point );
352
353 G4double stepLengthFB;
354 G4ThreeVector PointH;
355
356 // Check whether any volumes are encountered by the chord FB
357 // ---------------------------------------------------------
358
359 G4bool Intersects_FB = IntersectChord( CurrentF_Point, Point_B,
360 NewSafety,fPreviousSafety,
361 fPreviousSftOrigin,
362 stepLengthFB,
363 PointH );
364 if( Intersects_FB )
365 {
366 // There is an intersection of FB with a volume boundary
367 // H <- First Intersection of Chord FB
368
369 // H is our new Candidate for the intersection point.
370 // It replaces "E" and we will repeat the test to see if
371 // it is a good enough approximate point for us.
372
373 // Note that F must be in volume volA (the same as A)
374 // (otherwise AF would meet a volume boundary!)
375 // A <- F
376 // E <- H
377 //
378 G4FieldTrack InterMed=ApproxIntersecPointV;
379 ApproxIntersecPointV = GetChordFinderFor()->ApproxCurvePointS(
380 CurrentA_PointVelocity,CurrentB_PointVelocity,
381 InterMed,CurrentE_Point,CurrentF_Point,PointH,
382 false,GetEpsilonStepFor());
383 CurrentA_PointVelocity = InterMed;
384 CurrentE_Point = PointH;
385 }
386 else // not Intersects_FB
387 {
388 // There is NO intersection of FB with a volume boundary
389
390 if( fin_section_depth[depth] )
391 {
392 // If B is the original endpoint, this means that whatever
393 // volume(s) intersected the original chord, none touch the
394 // smaller chords we have used.
395 // The value of 'IntersectedOrRecalculatedFT' returned is
396 // likely not valid
397
398 // Check on real final_section or SubEndSection
399 //
400 if( ((Second_half)&&(depth==0)) || (first_section) )
401 {
402 there_is_no_intersection = true; // real final_section
403 }
404 else
405 {
406 // end of subsection, not real final section
407 // exit from the and go to the depth-1 level
408
409 substep_no_p = param_substeps+2; // exit from the loop
410
411 // but 'Second_half' is still true because we need to find
412 // the 'CurrentE_point' for the next loop
413 //
414 Second_half = true;
415 sub_final_section = true;
416 }
417 }
418 else
419 {
420 if(depth==0)
421 {
422 // We must restore the original endpoint
423 //
424 CurrentA_PointVelocity = CurrentB_PointVelocity; // Got to B
425 CurrentB_PointVelocity = CurveEndPointVelocity;
426 SubStart_PointVelocity = CurrentA_PointVelocity;
427 ApproxIntersecPointV = GetChordFinderFor()
428 ->ApproxCurvePointV( CurrentA_PointVelocity,
429 CurrentB_PointVelocity,
430 CurrentE_Point,
431 GetEpsilonStepFor());
432
433 restoredFullEndpoint = true;
434 restartB++; // counter
435 }
436 else
437 {
438 // We must restore the depth endpoint
439 //
440 CurrentA_PointVelocity = CurrentB_PointVelocity; // Got to B
441 CurrentB_PointVelocity = *ptrInterMedFT[depth];
442 SubStart_PointVelocity = CurrentA_PointVelocity;
443 ApproxIntersecPointV = GetChordFinderFor()
444 ->ApproxCurvePointV( CurrentA_PointVelocity,
445 CurrentB_PointVelocity,
446 CurrentE_Point,
447 GetEpsilonStepFor());
448 restoredFullEndpoint = true;
449 restartB++; // counter
450 }
451 }
452 } // Endif (Intersects_FB)
453 } // Endif (Intersects_AF)
454
455 // Ensure that the new endpoints are not further apart in space
456 // than on the curve due to different errors in the integration
457 //
458 G4double linDistSq, curveDist;
459 linDistSq = ( CurrentB_PointVelocity.GetPosition()
460 - CurrentA_PointVelocity.GetPosition() ).mag2();
461 curveDist = CurrentB_PointVelocity.GetCurveLength()
462 - CurrentA_PointVelocity.GetCurveLength();
463
464 // Change this condition for very strict parameters of propagation
465 //
466 if( curveDist*curveDist*(1+2* GetEpsilonStepFor()) < linDistSq )
467 {
468 // Re-integrate to obtain a new B
469 //
470 G4FieldTrack newEndPointFT=
471 ReEstimateEndpoint( CurrentA_PointVelocity,
472 CurrentB_PointVelocity,
473 linDistSq, // to avoid recalculation
474 curveDist );
475 G4FieldTrack oldPointVelB = CurrentB_PointVelocity;
476 CurrentB_PointVelocity = newEndPointFT;
477
478 if ( (fin_section_depth[depth]) // real final section
479 &&( first_section || ((Second_half)&&(depth==0)) ) )
480 {
481 recalculatedEndPoint = true;
482 IntersectedOrRecalculatedFT = newEndPointFT;
483 // So that we can return it, if it is the endpoint!
484 }
485 }
486 if( curveDist < 0.0 )
487 {
488 G4cerr << "ERROR - G4BrentLocator::EstimateIntersectionPoint()"
489 << G4endl
490 << " Error in advancing propagation." << G4endl;
491 fVerboseLevel = 5; // Print out a maximum of information
492 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
493 -1.0, NewSafety, substep_no );
494 G4cerr << " Point A (start) is " << CurrentA_PointVelocity
495 << G4endl;
496 G4cerr << " Point B (end) is " << CurrentB_PointVelocity
497 << G4endl;
498 G4cerr << " Curve distance is " << curveDist << G4endl;
499 G4cerr << G4endl
500 << "The final curve point is not further along"
501 << " than the original!" << G4endl;
502
503 if( recalculatedEndPoint )
504 {
505 G4cerr << "Recalculation of EndPoint was called with fEpsStep= "
506 << GetEpsilonStepFor() << G4endl;
507 }
508 oldprc = G4cerr.precision(20);
509 G4cerr << " Point A (Curve start) is " << CurveStartPointVelocity
510 << G4endl;
511 G4cerr << " Point B (Curve end) is " << CurveEndPointVelocity
512 << G4endl;
513 G4cerr << " Point A (Current start) is " << CurrentA_PointVelocity
514 << G4endl;
515 G4cerr << " Point B (Current end) is " << CurrentB_PointVelocity
516 << G4endl;
517 G4cerr << " Point S (Sub start) is " << SubStart_PointVelocity
518 << G4endl;
519 G4cerr << " Point E (Trial Point) is " << CurrentE_Point
520 << G4endl;
521 G4cerr << " Old Point F(Intersection) is " << CurrentF_Point
522 << G4endl;
523 G4cerr << " New Point F(Intersection) is " << ApproxIntersecPointV
524 << G4endl;
525 G4cerr << " LocateIntersection parameters are : Substep no= "
526 << substep_no << G4endl;
527 G4cerr << " Substep depth no= "<< substep_no_p << " Depth= "
528 << depth << G4endl;
529 G4cerr << " Restarted no= "<< restartB << " Epsilon= "
530 << GetEpsilonStepFor() <<" DeltaInters= "
531 << GetDeltaIntersectionFor() << G4endl;
532 G4cerr.precision( oldprc );
533
534 G4Exception("G4BrentLocator::EstimateIntersectionPoint()",
535 "FatalError", FatalException,
536 "Error in advancing propagation.");
537 }
538
539 if(restoredFullEndpoint)
540 {
541 fin_section_depth[depth] = restoredFullEndpoint;
542 restoredFullEndpoint = false;
543 }
544 } // EndIf ( E is close enough to the curve, ie point F. )
545 // tests ChordAF_Vector.mag() <= maximum_lateral_displacement
546
547#ifdef G4DEBUG_LOCATE_INTERSECTION
548 if( substep_no >= trigger_substepno_print )
549 {
550 G4cout << "Difficulty in converging in "
551 << "G4BrentLocator::EstimateIntersectionPoint()"
552 << G4endl
553 << " Substep no = " << substep_no << G4endl;
554 if( substep_no == trigger_substepno_print )
555 {
556 printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
557 -1.0, NewSafety, 0);
558 }
559 G4cout << " State of point A: ";
560 printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
561 -1.0, NewSafety, substep_no-1, 0);
562 G4cout << " State of point B: ";
563 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
564 -1.0, NewSafety, substep_no);
565 }
566#endif
567 substep_no++;
568 substep_no_p++;
569
570 } while ( ( ! found_approximate_intersection )
571 && ( ! there_is_no_intersection )
572 && ( substep_no_p <= param_substeps) ); // UNTIL found or
573 // failed param substep
574 first_section = false;
575
576 if( (!found_approximate_intersection) && (!there_is_no_intersection) )
577 {
578 G4double did_len = std::abs( CurrentA_PointVelocity.GetCurveLength()
579 - SubStart_PointVelocity.GetCurveLength());
580 G4double all_len = std::abs( CurrentB_PointVelocity.GetCurveLength()
581 - SubStart_PointVelocity.GetCurveLength());
582
583 G4double stepLengthAB;
584 G4ThreeVector PointGe;
585
586 // Check if progress is too slow and if it possible to go deeper,
587 // then halve the step if so
588 //
589 if ( ( did_len < fraction_done*all_len )
590 && (depth<max_depth) && (!sub_final_section) )
591 {
592 Second_half=false;
593 depth++;
594
595 G4double Sub_len = (all_len-did_len)/(2.);
596 G4FieldTrack start = CurrentA_PointVelocity;
597 G4MagInt_Driver* integrDriver =
598 GetChordFinderFor()->GetIntegrationDriver();
599 integrDriver->AccurateAdvance(start, Sub_len, GetEpsilonStepFor());
600 *ptrInterMedFT[depth] = start;
601 CurrentB_PointVelocity = *ptrInterMedFT[depth];
602
603 // Adjust 'SubStartPoint' to calculate the 'did_length' in next loop
604 //
605 SubStart_PointVelocity = CurrentA_PointVelocity;
606
607 // Find new trial intersection point needed at start of the loop
608 //
609 G4ThreeVector Point_A = CurrentA_PointVelocity.GetPosition();
610 G4ThreeVector SubE_point = CurrentB_PointVelocity.GetPosition();
611
612 GetNavigatorFor()->LocateGlobalPointWithinVolume(Point_A);
613 G4bool Intersects_AB = IntersectChord(Point_A, SubE_point,
614 NewSafety, fPreviousSafety,
615 fPreviousSftOrigin,stepLengthAB,
616 PointGe);
617 if( Intersects_AB )
618 {
619 last_AF_intersection = Intersects_AB;
620 CurrentE_Point = PointGe;
621 fin_section_depth[depth]=true;
622 }
623 else
624 {
625 // No intersection found for first part of curve
626 // (CurrentA,InterMedPoint[depth]). Go to the second part
627 //
628 Second_half = true;
629 }
630 } // if did_len
631
632 if( (Second_half)&&(depth!=0) )
633 {
634 // Second part of curve (InterMed[depth],Intermed[depth-1]) )
635 // On the depth-1 level normally we are on the 'second_half'
636
637 Second_half = true;
638
639 // Find new trial intersection point needed at start of the loop
640 //
641 SubStart_PointVelocity = *ptrInterMedFT[depth];
642 CurrentA_PointVelocity = *ptrInterMedFT[depth];
643 CurrentB_PointVelocity = *ptrInterMedFT[depth-1];
644 // Ensure that the new endpoints are not further apart in space
645 // than on the curve due to different errors in the integration
646 //
647 G4double linDistSq, curveDist;
648 linDistSq = ( CurrentB_PointVelocity.GetPosition()
649 - CurrentA_PointVelocity.GetPosition() ).mag2();
650 curveDist = CurrentB_PointVelocity.GetCurveLength()
651 - CurrentA_PointVelocity.GetCurveLength();
652 if( curveDist*curveDist*(1+2*GetEpsilonStepFor() ) < linDistSq )
653 {
654 // Re-integrate to obtain a new B
655 //
656 G4FieldTrack newEndPointFT=
657 ReEstimateEndpoint( CurrentA_PointVelocity,
658 CurrentB_PointVelocity,
659 linDistSq, // to avoid recalculation
660 curveDist );
661 G4FieldTrack oldPointVelB = CurrentB_PointVelocity;
662 CurrentB_PointVelocity = newEndPointFT;
663 if (depth==1)
664 {
665 recalculatedEndPoint = true;
666 IntersectedOrRecalculatedFT = newEndPointFT;
667 // So that we can return it, if it is the endpoint!
668 }
669 }
670
671
672 G4ThreeVector Point_A = CurrentA_PointVelocity.GetPosition();
673 G4ThreeVector SubE_point = CurrentB_PointVelocity.GetPosition();
674 GetNavigatorFor()->LocateGlobalPointWithinVolume(Point_A);
675 G4bool Intersects_AB = IntersectChord(Point_A, SubE_point, NewSafety,
676 fPreviousSafety,
677 fPreviousSftOrigin,stepLengthAB, PointGe);
678 if( Intersects_AB )
679 {
680 last_AF_intersection = Intersects_AB;
681 CurrentE_Point = PointGe;
682 }
683
684 depth--;
685 fin_section_depth[depth]=true;
686 }
687 } // if(!found_aproximate_intersection)
688
689 } while ( ( ! found_approximate_intersection )
690 && ( ! there_is_no_intersection )
691 && ( substep_no <= max_substeps) ); // UNTIL found or failed
692
693 if( substep_no > max_no_seen )
694 {
695 max_no_seen = substep_no;
696 if( max_no_seen > warn_substeps )
697 {
698 trigger_substepno_print = max_no_seen-20; // Want to see last 20 steps
699 }
700 }
701
702 if( ( substep_no >= max_substeps)
703 && !there_is_no_intersection
704 && !found_approximate_intersection )
705 {
706 G4cerr << "WARNING - G4BrentLocator::EstimateIntersectionPoint()"
707 << G4endl
708 << " Convergence is requiring too many substeps: "
709 << substep_no << G4endl;
710 G4cerr << " Abandoning effort to intersect. " << G4endl;
711 G4cerr << " Information on start & current step follows in cout."
712 << G4endl;
713 G4cout << "WARNING - G4BrentLocator::EstimateIntersectionPoint()"
714 << G4endl
715 << " Convergence is requiring too many substeps: "
716 << substep_no << G4endl;
717 G4cout << " Found intersection = "
718 << found_approximate_intersection << G4endl
719 << " Intersection exists = "
720 << !there_is_no_intersection << G4endl;
721 G4cout << " Start and Endpoint of Requested Step:" << G4endl;
722 printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
723 -1.0, NewSafety, 0);
724 G4cout << G4endl;
725 G4cout << " 'Bracketing' starting and endpoint of current Sub-Step"
726 << G4endl;
727 printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
728 -1.0, NewSafety, substep_no-1);
729 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
730 -1.0, NewSafety, substep_no);
731 G4cout << G4endl;
732 oldprc = G4cout.precision( 10 );
733 G4double done_len = CurrentA_PointVelocity.GetCurveLength();
734 G4double full_len = CurveEndPointVelocity.GetCurveLength();
735 G4cout << "ERROR - G4BrentLocator::EstimateIntersectionPoint()"
736 << G4endl
737 << " Undertaken only length: " << done_len
738 << " out of " << full_len << " required." << G4endl;
739 G4cout << " Remaining length = " << full_len - done_len << G4endl;
740 G4cout.precision( oldprc );
741
742 G4Exception("G4BrentLocator::EstimateIntersectionPoint()",
743 "UnableToLocateIntersection", FatalException,
744 "Too many substeps while trying to locate intersection.");
745 }
746 else if( substep_no >= warn_substeps )
747 {
748 oldprc= G4cout.precision( 10 );
749 G4cout << "WARNING - G4BrentLocator::EstimateIntersectionPoint()"
750 << G4endl
751 << " Undertaken length: "
752 << CurrentB_PointVelocity.GetCurveLength();
753 G4cout << " - Needed: " << substep_no << " substeps." << G4endl
754 << " Warning level = " << warn_substeps
755 << " and maximum substeps = " << max_substeps << G4endl;
756 G4Exception("G4BrentLocator::EstimateIntersectionPoint()",
757 "DifficultyToLocateIntersection", JustWarning,
758 "Many substeps while trying to locate intersection.");
759 G4cout.precision( oldprc );
760 }
761 return !there_is_no_intersection; // Success or failure
762}
Note: See TracBrowser for help on using the repository browser.