| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // Class G4SimpleLocator implementation
27 //
28 // 27.10.08 - Tatiana Nikitina, extracted from
29 // 04.10.11 - John Apostolakis, revised conver
30 // -------------------------------------------
31
32 #include <iomanip>
33
34 #include "G4ios.hh"
35 #include "G4SimpleLocator.hh"
36
37 G4SimpleLocator::G4SimpleLocator(G4Navigator *
38 : G4VIntersectionLocator(theNavigator)
39 {
40 }
41
42 G4SimpleLocator::~G4SimpleLocator() = default;
43
44 // -------------------------------------------
45 // G4bool G4PropagatorInField::LocateIntersect
46 // const G4FieldTrack& CurveStart
47 // const G4FieldTrack& CurveEndPo
48 // const G4ThreeVector& TrialPoint
49 // G4FieldTrack& Intersecte
50 // G4bool& recalculat
51 // -------------------------------------------
52 //
53 // Function that returns the intersection of t
54 // of the current volume (either the external
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
61 //
62 // E is the first point of intersection of
63 // a volume other than A (on the surface
64 //
65 // Convention of Use :
66 // i) If it returns "true", then Intersect
67 // to the approximate intersection point
68 // ii) If it returns "false", no intersecti
69 // The validity of IntersectedOrRecal
70 // a) if latter is false, then Intersec
71 // b) if latter is true, then Intersec
72 // the new endpoint, due to a re-i
73 // -------------------------------------------
74 // NOTE: implementation taken from G4Propagato
75 //
76 G4bool G4SimpleLocator::EstimateIntersectionPo
77 const G4FieldTrack& CurveStart
78 const G4FieldTrack& CurveEndPo
79 const G4ThreeVector& TrialPoint
80 G4FieldTrack& Intersecte
81 G4bool& recalculat
82 G4double &fPrevious
83 G4ThreeVector &fPrevious
84 {
85 // Find Intersection Point ( A, B, E ) of t
86
87 G4bool found_approximate_intersection = fals
88 G4bool there_is_no_intersection = fals
89
90 G4FieldTrack CurrentA_PointVelocity = Curve
91 G4FieldTrack CurrentB_PointVelocity = Curve
92 G4ThreeVector CurrentE_Point = TrialPoint;
93 G4bool validNormalAtE = false;
94 G4ThreeVector NormalAtEntry;
95
96 G4FieldTrack ApproxIntersecPointV(CurveEndP
97 G4double NewSafety = 0.0;
98 G4bool last_AF_intersection = false;
99 G4bool final_section = true; // Shows wheth
100 // (i.e. B=ful
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; // T
110 const G4int warn_substeps = 1000; //
111
112 // Statistics for substeps
113 //
114 static G4ThreadLocal G4int max_no_seen= -1;
115
116 NormalAtEntry = GetSurfaceNormal( CurrentE_P
117
118 #ifdef G4DEBUG_FIELD
119 const G4double tolerance = 1.0e-8;
120 G4ThreeVector StartPosition= CurveStartPoin
121 if( (TrialPoint - StartPosition).mag() < tol
122 {
123 G4Exception("G4SimpleLocator::EstimateInt
124 "GeomNav1002", JustWarning,
125 "Intersection point F is exac
126 }
127 #endif
128
129 do
130 {
131 G4ThreeVector Point_A = CurrentA_PointVel
132 G4ThreeVector Point_B = CurrentB_PointVel
133
134 // F = a point on true AB path close to p
135 // (the closest if possible)
136 //
137 ApproxIntersecPointV = GetChordFinderFor(
138 ->ApproxCurvePoint
139
140
141
142 // The above method is the key & most
143
144 #ifdef G4DEBUG_FIELD
145 if( ApproxIntersecPointV.GetCurveLength()
146 CurrentB_PointVelocity.GetCurveLength
147 {
148 G4Exception("G4SimpleLocator::EstimateI
149 "GeomNav0003", FatalExcepti
150 "Intermediate F point is pa
151 }
152 #endif
153
154 G4ThreeVector CurrentF_Point = ApproxInte
155
156 // First check whether EF is small - then
157 // Calculate the length and direction of
158 //
159 G4ThreeVector ChordEF_Vector = CurrentF_P
160
161 G4ThreeVector NewMomentumDir = ApproxInte
162 G4double MomDir_dot_Norm = NewMomentumDir
163
164 G4ThreeVector ChordAB = Point_B - Point_A
165
166 #ifdef G4DEBUG_FIELD
167 G4VIntersectionLocator::
168 ReportTrialStep( substep_no, ChordAB, C
169 NewMomentumDir, Normal
170 #endif
171 // Check Sign is always exiting !! TODO
172 // Could ( > -epsilon) be used instead?
173 //
174 G4bool adequate_angle = ( MomDir_dot_Norm
175 || (! validNormalAtE
176 G4double EF_dist2= ChordEF_Vector.mag2();
177 if ( ( EF_dist2 <= sqr(fiDeltaIntersecti
178 || ( EF_dist2 <= kCarTolerance*kCarTole
179 {
180 found_approximate_intersection = true;
181
182 // Create the "point" return value
183 //
184 IntersectedOrRecalculatedFT = ApproxInt
185 IntersectedOrRecalculatedFT.SetPosition
186
187 if ( GetAdjustementOfFoundIntersection(
188 {
189 // Try to Get Correction of Intersect
190 //
191 G4ThreeVector IP;
192 G4ThreeVector MomentumDir= ApproxInte
193 G4bool goodCorrection = AdjustmentOfF
194 CurrentE_Po
195 last_AF_int
196 fPreviousSa
197
198 if ( goodCorrection )
199 {
200 IntersectedOrRecalculatedFT = Appro
201 IntersectedOrRecalculatedFT.SetPosi
202 }
203 }
204
205 // Note: in order to return a point on
206 // we must return E. But it is F
207 // So we must "cheat": we are usi
208 // and the velocity at point F !!
209 //
210 // This must limit the length we can al
211 }
212 else // E is NOT close enough to the cur
213 {
214 // Check whether any volumes are encoun
215 // ------------------------------------
216 // First relocate to restore any Voxel
217 // in the Navigator before calling Comp
218 //
219 GetNavigatorFor()->LocateGlobalPointWit
220
221 G4ThreeVector PointG; // Candidate in
222 G4double stepLengthAF;
223 G4bool usedNavigatorAF = false;
224 G4bool Intersects_AF = IntersectChord(
225
226
227
228
229
230
231
232 last_AF_intersection = Intersects_AF;
233 if( Intersects_AF )
234 {
235 // G is our new Candidate for the int
236 // It replaces "E" and we will repea
237 // it is a good enough approximate po
238 // B <- F
239 // E <- G
240
241 CurrentB_PointVelocity = ApproxInters
242 CurrentE_Point = PointG;
243
244 // Need to recalculate the Exit Norma
245 // Relies on a call to Navigator::Com
246 // If the safety was adequate (for th
247 // But this must not occur, no inters
248 // so this branch, ie if( Intersects_
249 //
250 G4bool validNormalLast;
251 NormalAtEntry = GetSurfaceNormal( Po
252 validNormalAtE = validNormalLast;
253
254 // By moving point B, must take care
255 // AF has no intersection to try curr
256 //
257 final_section = false;
258
259 #ifdef G4VERBOSE
260 if( fVerboseLevel > 3 )
261 {
262 G4cout << "G4PiF::LI> Investigating
263 << " at s=" << ApproxInterse
264 << " on way to full s="
265 << CurveEndPointVelocity.Get
266 }
267 #endif
268 }
269 else // not Intersects_AF
270 {
271 // In this case:
272 // There is NO intersection of AF wit
273 // We must continue the search in the
274 //
275 GetNavigatorFor()->LocateGlobalPointW
276
277 G4double stepLengthFB;
278 G4ThreeVector PointH;
279 G4bool usedNavigatorFB = false;
280
281 // Check whether any volumes are enco
282 // ----------------------------------
283
284 G4bool Intersects_FB = IntersectChord
285
286
287
288
289 if( Intersects_FB )
290 {
291 // There is an intersection of FB w
292 // H <- First Intersection of Chord
293
294 // H is our new Candidate for the i
295 // It replaces "E" and we will rep
296 // it is a good enough approximate
297
298 // Note that F must be in volume vo
299 // (otherwise AF would meet a volum
300 // A <- F
301 // E <- H
302 //
303 CurrentA_PointVelocity = ApproxInte
304 CurrentE_Point = PointH;
305
306 // Need to recalculate the Exit Nor
307 // Relies on call to Navigator::Com
308 // If safety was adequate (for the
309 // But this must not occur, no inte
310 // so this branch, i.e. if( Interse
311 //
312 G4bool validNormalLast;
313 NormalAtEntry = GetSurfaceNormal(
314 validNormalAtE = validNormalLast;
315 }
316 else // not Intersects_FB
317 {
318 // There is NO intersection of FB w
319
320 if( final_section )
321 {
322 // If B is the original endpoint,
323 // volume(s) intersected the orig
324 // smaller chords we have used.
325 // The value of 'IntersectedOrRec
326 // likely not valid
327
328 there_is_no_intersection = true;
329 }
330 else
331 {
332 // We must restore the original e
333
334 CurrentA_PointVelocity = CurrentB
335 CurrentB_PointVelocity = CurveEnd
336 restoredFullEndpoint = true;
337 }
338 } // Endif (Intersects_FB)
339 } // Endif (Intersects_AF)
340
341 // Ensure that the new endpoints are no
342 // than on the curve due to different e
343 //
344 G4double linDistSq, curveDist;
345 linDistSq = ( CurrentB_PointVelocity.Ge
346 - CurrentA_PointVelocity.Ge
347 curveDist = CurrentB_PointVelocity.GetC
348 - CurrentA_PointVelocity.Ge
349
350 // Change this condition for very stric
351 //
352 if( curveDist*curveDist*(1+2* GetEpsilo
353 {
354 // Re-integrate to obtain a new B
355 //
356 G4FieldTrack newEndPointFT =
357 ReEstimateEndpoint( CurrentA_
358 CurrentB_
359 linDistSq
360 curveDist
361 G4FieldTrack oldPointVelB = CurrentB_
362 CurrentB_PointVelocity = newEndPointF
363
364 if( (final_section)) // real final se
365 {
366 recalculatedEndPoint = true;
367 IntersectedOrRecalculatedFT = newEn
368 // So that we can return it, if i
369 }
370 }
371 if( curveDist < 0.0 )
372 {
373 fVerboseLevel = 5; // Print out a max
374 printStatus( CurrentA_PointVelocity,
375 -1.0, NewSafety, subste
376 std::ostringstream message;
377 message << "Error in advancing propag
378 << " Point A (start) i
379 << G4endl
380 << " Point B (end) i
381 << G4endl
382 << " Curve distance is
383 << G4endl
384 << "The final curve point is
385 << " than the original!" << G
386
387 if( recalculatedEndPoint )
388 {
389 message << "Recalculation of EndPoi
390 << GetEpsilonStepFor() << G
391 }
392 message.precision(20);
393 message << " Point A (Curve start)
394 << G4endl
395 << " Point B (Curve end)
396 << G4endl
397 << " Point A (Current start)
398 << G4endl
399 << " Point B (Current end)
400 << G4endl
401 << " Point E (Trial Point)
402 << G4endl
403 << " Point F (Intersection)
404 << G4endl
405 << " LocateIntersectio
406 << substep_no;
407
408 G4Exception("G4SimpleLocator::Estimat
409 "GeomNav0003", FatalExcep
410 }
411
412 if ( restoredFullEndpoint )
413 {
414 final_section = restoredFullEndpoint;
415 restoredFullEndpoint = false;
416 }
417 } // EndIf ( E is close enough to the cur
418 // tests ChordAF_Vector.mag() <= maximu
419
420 #ifdef G4DEBUG_LOCATE_INTERSECTION
421 G4int trigger_substepno_print= warn_subst
422
423 if( substep_no >= trigger_substepno_print
424 {
425 G4cout << "Difficulty in converging in
426 << "G4SimpleLocator::EstimateInt
427 << G4endl
428 << " Substep no = " << subste
429 if( substep_no == trigger_substepno_pri
430 {
431 printStatus( CurveStartPointVelocity,
432 -1.0, NewSafety, 0);
433 }
434 G4cout << " State of point A: ";
435 printStatus( CurrentA_PointVelocity, Cu
436 -1.0, NewSafety, substep_n
437 G4cout << " State of point B: ";
438 printStatus( CurrentA_PointVelocity, Cu
439 -1.0, NewSafety, substep_n
440 }
441 #endif
442 ++substep_no;
443
444 } while ( ( ! found_approximate_intersection
445 && ( ! there_is_no_intersection )
446 && ( substep_no <= max_substeps) );
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
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::Estima
464 << " Start and Endpoint of R
465 printStatus( CurveStartPointVelocity, Curv
466 -1.0, NewSafety, 0);
467 G4cout << G4endl
468 << " Start and end-point of
469 printStatus( CurrentA_PointVelocity, Curre
470 -1.0, NewSafety, substep_no-1
471 printStatus( CurrentA_PointVelocity, Curre
472 -1.0, NewSafety, substep_no);
473
474 std::ostringstream message;
475 message << "Convergence is requiring too m
476 << substep_no << G4endl
477 << " Abandoning effort to
478 << " Found intersection =
479 << found_approximate_intersection
480 << " Intersection exists
481 << !there_is_no_intersection << G4
482 message.precision(10);
483 G4double done_len = CurrentA_PointVelocity
484 G4double full_len = CurveEndPointVelocity.
485 message << " Undertaken only leng
486 << " out of " << full_len << " req
487 << " Remaining length = "
488
489 G4Exception("G4SimpleLocator::EstimateInte
490 "GeomNav0003", FatalException,
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
498 << " Undertaken length: "
499 << CurrentB_PointVelocity.GetCurve
500 << " - Needed: " << substep_no <<
501 << " Warning level = " <<
502 << " and maximum substeps = " << m
503 G4Exception("G4SimpleLocator::EstimateInte
504 "GeomNav1002", JustWarning, me
505 }
506 return !there_is_no_intersection; // Succe
507 }
508