Context Navigation

source: trunk/source/geometry/navigation/src/G4MultiLevelLocator.cc @ 1202

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: