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

Last change on this file since 1355 was 1347, checked in by garnier, 15 years ago

geant4 tag 9.4
File size: 20.2 KB
RevLine 
[985]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//
[1347]26// $Id: G4SimpleLocator.cc,v 1.6 2010/07/13 15:59:42 gcosmo Exp $
27// GEANT4 tag $Name: geant4-09-04-ref-00 $
[985]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
39G4SimpleLocator::G4SimpleLocator(G4Navigator *theNavigator)
40 : G4VIntersectionLocator(theNavigator)
41{
42}
43
44G4SimpleLocator::~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//
80G4bool 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
[1347]108 G4int oldprc; // cout/cerr precision settings
109
[985]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 {
[1058]125 G4cerr << "WARNING - G4SimpleLocator::EstimateIntersectionPoint()"
[985]126 << G4endl
127 << " Intermediate F point is on top of starting point A."
128 << G4endl;
[1058]129 G4Exception("G4SimpleLocator::EstimateIntersectionPoint()",
[985]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 {
[1058]154 G4cerr << "ERROR - G4SimpleLocator::EstimateIntersectionPoint()"
[985]155 << G4endl
156 << " Intermediate F point is more advanced than"
157 << " endpoint B." << G4endl;
[1058]158 G4Exception("G4SimpleLocator::EstimateIntersectionPoint()",
[985]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( curveDist*curveDist*(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 {
[1058]342 G4cerr << "ERROR - G4SimpleLocator::EstimateIntersectionPoint()"
[985]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 }
[1347]362 oldprc = G4cerr.precision(20);
[985]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;
[1347]377 G4cerr.precision(oldprc);
[985]378
[1058]379 G4Exception("G4SimpleLocator::EstimateIntersectionPoint()",
[985]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 "
[1058]396 << "G4SimpleLocator::EstimateIntersectionPoint():"
[985]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 {
[1058]431 G4cerr << "WARNING - G4SimpleLocator::EstimateIntersectionPoint()"
[985]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;
[1058]438 G4cout << "WARNING - G4SimpleLocator::EstimateIntersectionPoint()"
[985]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;
[1347]457 oldprc = G4cout.precision(10);
[985]458 G4double done_len = CurrentA_PointVelocity.GetCurveLength();
459 G4double full_len = CurveEndPointVelocity.GetCurveLength();
[1058]460 G4cout << "ERROR - G4SimpleLocator::EstimateIntersectionPoint()"
[985]461 << G4endl
462 << " Undertaken only length: " << done_len
463 << " out of " << full_len << " required." << G4endl;
464 G4cout << " Remaining length = " << full_len - done_len << G4endl;
[1347]465 G4cout.precision(oldprc);
[985]466
[1058]467 G4Exception("G4SimpleLocator::EstimateIntersectionPoint()",
[985]468 "UnableToLocateIntersection", FatalException,
469 "Too many substeps while trying to locate intersection.");
470 }
471 else if( substep_no >= warn_substeps )
472 {
[1347]473 oldprc= G4cout.precision(10);
[1058]474 G4cout << "WARNING - G4SimpleLocator::EstimateIntersectionPoint()"
[985]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;
[1347]481 G4cout.precision(oldprc);
[1058]482 G4Exception("G4SimpleLocator::EstimateIntersectionPoint()",
[985]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.