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

Last change on this file since 1242 was 1228, checked in by garnier, 16 years ago

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