Context Navigation

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

Last change on this file since 1337 was 1337, checked in by garnier, 14 years ago
tag geant4.9.4 beta 1 + modifs locales
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-04-beta-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
38	G4BrentLocator::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
63	G4BrentLocator::~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	//
115	G4bool 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( curveDistcurveDist(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( curveDistcurveDist(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.

Download in other formats: