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