Context Navigation

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

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

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

Download in other formats: