Context Navigation

G4BrentLocator.cc @ 1358

Last change on this file since 1358 was 1347, checked in by garnier, 14 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
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 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( curveDistcurveDist(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( curveDistcurveDist(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.

Context Navigation

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

Download in other formats: