Context Navigation

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

Last change on this file since 1337 was 1337, checked in by garnier, 14 years ago
tag geant4.9.4 beta 1 + modifs locales
File size: 20.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: G4SimpleLocator.cc,v 1.5 2008/12/11 10:27:58 tnikitin Exp $
27	// GEANT4 tag $Name: geant4-09-04-beta-01 $
28	//
29	// Class G4SimpleLocator implementation
30	//
31	// 27.10.08 - Tatiana Nikitina.
32	// ---------------------------------------------------------------------------
33
34	#include <iomanip>
35
36	#include "G4ios.hh"
37	#include "G4SimpleLocator.hh"
38
39	G4SimpleLocator::G4SimpleLocator(G4Navigator *theNavigator)
40	: G4VIntersectionLocator(theNavigator)
41	{
42	}
43
44	G4SimpleLocator::~G4SimpleLocator()
45	{
46	}
47
48	// --------------------------------------------------------------------------
49	// G4bool G4PropagatorInField::LocateIntersectionPoint(
50	// const G4FieldTrack& CurveStartPointVelocity, // A
51	// const G4FieldTrack& CurveEndPointVelocity, // B
52	// const G4ThreeVector& TrialPoint, // E
53	// G4FieldTrack& IntersectedOrRecalculated // Output
54	// G4bool& recalculated ) // Out
55	// --------------------------------------------------------------------------
56	//
57	// Function that returns the intersection of the true path with the surface
58	// of the current volume (either the external one or the inner one with one
59	// of the daughters:
60	//
61	// A = Initial point
62	// B = another point
63	//
64	// Both A and B are assumed to be on the true path:
65	//
66	// E is the first point of intersection of the chord AB with
67	// a volume other than A (on the surface of A or of a daughter)
68	//
69	// Convention of Use :
70	// i) If it returns "true", then IntersectionPointVelocity is set
71	// to the approximate intersection point.
72	// ii) If it returns "false", no intersection was found.
73	// The validity of IntersectedOrRecalculated depends on 'recalculated'
74	// a) if latter is false, then IntersectedOrRecalculated is invalid.
75	// b) if latter is true, then IntersectedOrRecalculated is
76	// the new endpoint, due to a re-integration.
77	// --------------------------------------------------------------------------
78	// NOTE: implementation taken from G4PropagatorInField
79	//
80	G4bool G4SimpleLocator::EstimateIntersectionPoint(
81	const G4FieldTrack& CurveStartPointVelocity, // A
82	const G4FieldTrack& CurveEndPointVelocity, // B
83	const G4ThreeVector& TrialPoint, // E
84	G4FieldTrack& IntersectedOrRecalculatedFT, // Output
85	G4bool& recalculatedEndPoint, // Out
86	G4double &fPreviousSafety, //In/Out
87	G4ThreeVector &fPreviousSftOrigin) //In/Out
88	{
89	// Find Intersection Point ( A, B, E ) of true path AB - start at E.
90
91	G4bool found_approximate_intersection = false;
92	G4bool there_is_no_intersection = false;
93
94	G4FieldTrack CurrentA_PointVelocity = CurveStartPointVelocity;
95	G4FieldTrack CurrentB_PointVelocity = CurveEndPointVelocity;
96	G4ThreeVector CurrentE_Point = TrialPoint;
97	G4FieldTrack ApproxIntersecPointV(CurveEndPointVelocity); // FT-Def-Construct
98	G4double NewSafety = 0.0;
99	G4bool last_AF_intersection = false;
100	G4bool final_section = true; // Shows whether current section is last
101	// (i.e. B=full end)
102	recalculatedEndPoint = false;
103
104	G4bool restoredFullEndpoint = false;
105
106	G4int substep_no = 0;
107
108	// Limits for substep number
109	//
110	const G4int max_substeps = 100000000; // Test 120 (old value 100 )
111	const G4int warn_substeps = 1000; // 100
112
113	// Statistics for substeps
114	//
115	static G4int max_no_seen= -1;
116	static G4int trigger_substepno_print= warn_substeps - 20;
117
118	#ifdef G4DEBUG_FIELD
119	static G4double tolerance= 1.0e-8;
120	G4ThreeVector StartPosition= CurveStartPointVelocity.GetPosition();
121	if( (TrialPoint - StartPosition).mag() < tolerance * mm )
122	{
123	G4cerr << "WARNING - G4SimpleLocator::EstimateIntersectionPoint()"
124	<< G4endl
125	<< " Intermediate F point is on top of starting point A."
126	<< G4endl;
127	G4Exception("G4SimpleLocator::EstimateIntersectionPoint()",
128	"IntersectionPointIsAtStart", JustWarning,
129	"Intersection point F is exactly at start point A." );
130	}
131	#endif
132
133	do
134	{
135	G4ThreeVector Point_A = CurrentA_PointVelocity.GetPosition();
136	G4ThreeVector Point_B = CurrentB_PointVelocity.GetPosition();
137
138	// F = a point on true AB path close to point E
139	// (the closest if possible)
140	//
141	ApproxIntersecPointV = GetChordFinderFor()
142	->ApproxCurvePointV( CurrentA_PointVelocity,
143	CurrentB_PointVelocity,
144	CurrentE_Point,
145	GetEpsilonStepFor());
146	// The above method is the key & most intuitive part ...
147
148	#ifdef G4DEBUG_FIELD
149	if( ApproxIntersecPointV.GetCurveLength() >
150	CurrentB_PointVelocity.GetCurveLength() * (1.0 + tolerance) )
151	{
152	G4cerr << "ERROR - G4SimpleLocator::EstimateIntersectionPoint()"
153	<< G4endl
154	<< " Intermediate F point is more advanced than"
155	<< " endpoint B." << G4endl;
156	G4Exception("G4SimpleLocator::EstimateIntersectionPoint()",
157	"IntermediatePointConfusion", FatalException,
158	"Intermediate F point is past end B point" );
159	}
160	#endif
161
162	G4ThreeVector CurrentF_Point= ApproxIntersecPointV.GetPosition();
163
164	// First check whether EF is small - then F is a good approx. point
165	// Calculate the length and direction of the chord AF
166	//
167	G4ThreeVector ChordEF_Vector = CurrentF_Point - CurrentE_Point;
168
169	if ( ChordEF_Vector.mag2() <= sqr(GetDeltaIntersectionFor()) )
170	{
171	found_approximate_intersection = true;
172
173	// Create the "point" return value
174	//
175
176	IntersectedOrRecalculatedFT = ApproxIntersecPointV;
177	IntersectedOrRecalculatedFT.SetPosition( CurrentE_Point );
178
179	if ( GetAdjustementOfFoundIntersection() )
180	{
181	// Try to Get Correction of IntersectionPoint using SurfaceNormal()
182	//
183	G4ThreeVector IP;
184	G4ThreeVector MomentumDir= ApproxIntersecPointV.GetMomentumDirection();
185	G4bool goodCorrection = AdjustmentOfFoundIntersection( Point_A,
186	CurrentE_Point, CurrentF_Point, MomentumDir,
187	last_AF_intersection, IP, NewSafety,
188	fPreviousSafety, fPreviousSftOrigin );
189
190	if(goodCorrection)
191	{
192	IntersectedOrRecalculatedFT = ApproxIntersecPointV;
193	IntersectedOrRecalculatedFT.SetPosition(IP);
194	}
195	}
196
197	// Note: in order to return a point on the boundary,
198	// we must return E. But it is F on the curve.
199	// So we must "cheat": we are using the position at point E
200	// and the velocity at point F !!!
201	//
202	// This must limit the length we can allow for displacement!
203	}
204	else // E is NOT close enough to the curve (ie point F)
205	{
206	// Check whether any volumes are encountered by the chord AF
207	// ---------------------------------------------------------
208	// First relocate to restore any Voxel etc information
209	// in the Navigator before calling ComputeStep()
210	//
211	GetNavigatorFor()->LocateGlobalPointWithinVolume( Point_A );
212
213	G4ThreeVector PointG; // Candidate intersection point
214	G4double stepLengthAF;
215	G4bool Intersects_AF = IntersectChord( Point_A, CurrentF_Point,
216	NewSafety,fPreviousSafety,
217	fPreviousSftOrigin,
218	stepLengthAF,
219	PointG );
220	last_AF_intersection = Intersects_AF;
221	if( Intersects_AF )
222	{
223	// G is our new Candidate for the intersection point.
224	// It replaces "E" and we will repeat the test to see if
225	// it is a good enough approximate point for us.
226	// B <- F
227	// E <- G
228
229	CurrentB_PointVelocity = ApproxIntersecPointV;
230	CurrentE_Point = PointG;
231
232	// By moving point B, must take care if current
233	// AF has no intersection to try current FB!!
234	//
235	final_section= false;
236
237	#ifdef G4VERBOSE
238	if( fVerboseLevel > 3 )
239	{
240	G4cout << "G4PiF::LI> Investigating intermediate point"
241	<< " at s=" << ApproxIntersecPointV.GetCurveLength()
242	<< " on way to full s="
243	<< CurveEndPointVelocity.GetCurveLength() << G4endl;
244	}
245	#endif
246	}
247	else // not Intersects_AF
248	{
249	// In this case:
250	// There is NO intersection of AF with a volume boundary.
251	// We must continue the search in the segment FB!
252	//
253	GetNavigatorFor()->LocateGlobalPointWithinVolume( CurrentF_Point );
254
255	G4double stepLengthFB;
256	G4ThreeVector PointH;
257
258	// Check whether any volumes are encountered by the chord FB
259	// ---------------------------------------------------------
260
261	G4bool Intersects_FB = IntersectChord( CurrentF_Point, Point_B,
262	NewSafety,fPreviousSafety,
263	fPreviousSftOrigin,
264	stepLengthFB,
265	PointH );
266	if( Intersects_FB )
267	{
268	// There is an intersection of FB with a volume boundary
269	// H <- First Intersection of Chord FB
270
271	// H is our new Candidate for the intersection point.
272	// It replaces "E" and we will repeat the test to see if
273	// it is a good enough approximate point for us.
274
275	// Note that F must be in volume volA (the same as A)
276	// (otherwise AF would meet a volume boundary!)
277	// A <- F
278	// E <- H
279	//
280	CurrentA_PointVelocity = ApproxIntersecPointV;
281	CurrentE_Point = PointH;
282	}
283	else // not Intersects_FB
284	{
285	// There is NO intersection of FB with a volume boundary
286
287	if( final_section )
288	{
289	// If B is the original endpoint, this means that whatever
290	// volume(s) intersected the original chord, none touch the
291	// smaller chords we have used.
292	// The value of 'IntersectedOrRecalculatedFT' returned is
293	// likely not valid
294
295	there_is_no_intersection = true; // real final_section
296	}
297	else
298	{
299	// We must restore the original endpoint
300
301	CurrentA_PointVelocity = CurrentB_PointVelocity; // Got to B
302	CurrentB_PointVelocity = CurveEndPointVelocity;
303	restoredFullEndpoint = true;
304	}
305	} // Endif (Intersects_FB)
306	} // Endif (Intersects_AF)
307
308	// Ensure that the new endpoints are not further apart in space
309	// than on the curve due to different errors in the integration
310	//
311	G4double linDistSq, curveDist;
312	linDistSq = ( CurrentB_PointVelocity.GetPosition()
313	- CurrentA_PointVelocity.GetPosition() ).mag2();
314	curveDist = CurrentB_PointVelocity.GetCurveLength()
315	- CurrentA_PointVelocity.GetCurveLength();
316
317	// Change this condition for very strict parameters of propagation
318	//
319	if( curveDistcurveDist(1+2* GetEpsilonStepFor()) < linDistSq )
320	{
321	// Re-integrate to obtain a new B
322	//
323	G4FieldTrack newEndPointFT =
324	ReEstimateEndpoint( CurrentA_PointVelocity,
325	CurrentB_PointVelocity,
326	linDistSq, // to avoid recalculation
327	curveDist );
328	G4FieldTrack oldPointVelB = CurrentB_PointVelocity;
329	CurrentB_PointVelocity = newEndPointFT;
330
331	if( (final_section)) // real final section
332	{
333	recalculatedEndPoint = true;
334	IntersectedOrRecalculatedFT = newEndPointFT;
335	// So that we can return it, if it is the endpoint!
336	}
337	}
338	if( curveDist < 0.0 )
339	{
340	G4cerr << "ERROR - G4SimpleLocator::EstimateIntersectionPoint()"
341	<< G4endl
342	<< " Error in advancing propagation." << G4endl;
343	fVerboseLevel = 5; // Print out a maximum of information
344	printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
345	-1.0, NewSafety, substep_no );
346	G4cerr << " Point A (start) is " << CurrentA_PointVelocity
347	<< G4endl;
348	G4cerr << " Point B (end) is " << CurrentB_PointVelocity
349	<< G4endl;
350	G4cerr << " Curve distance is " << curveDist << G4endl;
351	G4cerr << G4endl
352	<< "The final curve point is not further along"
353	<< " than the original!" << G4endl;
354
355	if( recalculatedEndPoint )
356	{
357	G4cerr << "Recalculation of EndPoint was called with fEpsStep= "
358	<< GetEpsilonStepFor() << G4endl;
359	}
360	G4cerr.precision(20);
361	G4cerr << " Point A (Curve start) is " << CurveStartPointVelocity
362	<< G4endl;
363	G4cerr << " Point B (Curve end) is " << CurveEndPointVelocity
364	<< G4endl;
365	G4cerr << " Point A (Current start) is " << CurrentA_PointVelocity
366	<< G4endl;
367	G4cerr << " Point B (Current end) is " << CurrentB_PointVelocity
368	<< G4endl;
369	G4cerr << " Point E (Trial Point) is " << CurrentE_Point
370	<< G4endl;
371	G4cerr << " Point F (Intersection) is " << ApproxIntersecPointV
372	<< G4endl;
373	G4cerr << " LocateIntersection parameters are : Substep no= "
374	<< substep_no << G4endl;
375
376	G4Exception("G4SimpleLocator::EstimateIntersectionPoint()",
377	"FatalError", FatalException,
378	"Error in advancing propagation.");
379	}
380
381	if(restoredFullEndpoint)
382	{
383	final_section = restoredFullEndpoint;
384	restoredFullEndpoint = false;
385	}
386	} // EndIf ( E is close enough to the curve, ie point F. )
387	// tests ChordAF_Vector.mag() <= maximum_lateral_displacement
388
389	#ifdef G4DEBUG_LOCATE_INTERSECTION
390	if( substep_no >= trigger_substepno_print )
391	{
392	G4cout << "Difficulty in converging in "
393	<< "G4SimpleLocator::EstimateIntersectionPoint():"
394	<< G4endl
395	<< " Substep no = " << substep_no << G4endl;
396	if( substep_no == trigger_substepno_print )
397	{
398	printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
399	-1.0, NewSafety, 0);
400	}
401	G4cout << " State of point A: ";
402	printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
403	-1.0, NewSafety, substep_no-1, 0);
404	G4cout << " State of point B: ";
405	printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
406	-1.0, NewSafety, substep_no);
407	}
408	#endif
409	substep_no++;
410
411	} while ( ( ! found_approximate_intersection )
412	&& ( ! there_is_no_intersection )
413	&& ( substep_no <= max_substeps) ); // UNTIL found or failed
414
415	if( substep_no > max_no_seen )
416	{
417	max_no_seen = substep_no;
418	if( max_no_seen > warn_substeps )
419	{
420	trigger_substepno_print = max_no_seen-20; // Want to see last 20 steps
421	}
422	}
423
424	if( ( substep_no >= max_substeps)
425	&& !there_is_no_intersection
426	&& !found_approximate_intersection )
427	{
428	G4cerr << "WARNING - G4SimpleLocator::EstimateIntersectionPoint()"
429	<< G4endl
430	<< " Convergence is requiring too many substeps: "
431	<< substep_no << G4endl;
432	G4cerr << " Abandoning effort to intersect. " << G4endl;
433	G4cerr << " Information on start & current step follows in cout."
434	<< G4endl;
435	G4cout << "WARNING - G4SimpleLocator::EstimateIntersectionPoint()"
436	<< G4endl
437	<< " Convergence is requiring too many substeps: "
438	<< substep_no << G4endl;
439	G4cout << " Found intersection = "
440	<< found_approximate_intersection << G4endl
441	<< " Intersection exists = "
442	<< !there_is_no_intersection << G4endl;
443	G4cout << " Start and Endpoint of Requested Step:" << G4endl;
444	printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
445	-1.0, NewSafety, 0);
446	G4cout << G4endl;
447	G4cout << " 'Bracketing' starting and endpoint of current Sub-Step"
448	<< G4endl;
449	printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
450	-1.0, NewSafety, substep_no-1);
451	printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
452	-1.0, NewSafety, substep_no);
453	G4cout << G4endl;
454	G4cout.precision( 10 );
455	G4double done_len = CurrentA_PointVelocity.GetCurveLength();
456	G4double full_len = CurveEndPointVelocity.GetCurveLength();
457	G4cout << "ERROR - G4SimpleLocator::EstimateIntersectionPoint()"
458	<< G4endl
459	<< " Undertaken only length: " << done_len
460	<< " out of " << full_len << " required." << G4endl;
461	G4cout << " Remaining length = " << full_len - done_len << G4endl;
462
463	G4Exception("G4SimpleLocator::EstimateIntersectionPoint()",
464	"UnableToLocateIntersection", FatalException,
465	"Too many substeps while trying to locate intersection.");
466	}
467	else if( substep_no >= warn_substeps )
468	{
469	G4int oldprc= G4cout.precision( 10 );
470	G4cout << "WARNING - G4SimpleLocator::EstimateIntersectionPoint()"
471	<< G4endl
472	<< " Undertaken length: "
473	<< CurrentB_PointVelocity.GetCurveLength();
474	G4cout << " - Needed: " << substep_no << " substeps." << G4endl
475	<< " Warning level = " << warn_substeps
476	<< " and maximum substeps = " << max_substeps << G4endl;
477	G4Exception("G4SimpleLocator::EstimateIntersectionPoint()",
478	"DifficultyToLocateIntersection", JustWarning,
479	"Many substeps while trying to locate intersection.");
480	G4cout.precision( oldprc );
481	}
482	return !there_is_no_intersection; // Success or failure
483	}

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

Download in other formats: