Context Navigation

G4SimpleLocator.cc @ 1347

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

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

Context Navigation

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

Download in other formats: