Context Navigation

G4MultiLevelLocator.cc @ 1347

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

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

Context Navigation

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

Download in other formats: