| 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 G4MultiLevelLocator implementation
27 //
28 // 27.10.08 - Tatiana Nikitina.
29 // 04.10.11 - John Apostolakis, revised conver
30 // -------------------------------------------
31
32 #include <iomanip>
33
34 #include "G4ios.hh"
35 #include "G4MultiLevelLocator.hh"
36 #include "G4LocatorChangeRecord.hh"
37 #include "G4LocatorChangeLogger.hh"
38
39 G4MultiLevelLocator::G4MultiLevelLocator(G4Nav
40 : G4VIntersectionLocator(theNavigator)
41 {
42 // In case of too slow progress in finding I
43 // intermediates Points on the Track must be
44 // Initialise the array of Pointers [max_dep
45
46 G4ThreeVector zeroV(0.0,0.0,0.0);
47 for (auto & idepth : ptrInterMedFT)
48 {
49 idepth = new G4FieldTrack( zeroV, zeroV, 0
50 }
51
52 if (fCheckMode)
53 {
54 // Trial values Loose Mediu
55 // To happen: Infrequent Occas
56 SetMaxSteps(150); // 300 85
57 SetWarnSteps(80); // 250 60
58 }
59 }
60
61 G4MultiLevelLocator::~G4MultiLevelLocator()
62 {
63 for (auto & idepth : ptrInterMedFT)
64 {
65 delete idepth;
66 }
67 #ifdef G4DEBUG_FIELD
68 ReportStatistics();
69 #endif
70 }
71
72
73 // -------------------------------------------
74 // G4bool G4PropagatorInField::LocateIntersect
75 // const G4FieldTrack& CurveStart
76 // const G4FieldTrack& CurveEndPo
77 // const G4ThreeVector& TrialPoint
78 // G4FieldTrack& Intersecte
79 // G4bool& recalculat
80 // -------------------------------------------
81 //
82 // Function that returns the intersection of t
83 // of the current volume (either the external
84 // of the daughters:
85 //
86 // A = Initial point
87 // B = another point
88 //
89 // Both A and B are assumed to be on the true
90 //
91 // E is the first point of intersection of
92 // a volume other than A (on the surface
93 //
94 // Convention of Use :
95 // i) If it returns "true", then Intersect
96 // to the approximate intersection point
97 // ii) If it returns "false", no intersecti
98 // Potential reasons:
99 // a) no segment found an intersection
100 // b) too many steps were required - af
101 // and is returning how far it could
102 // (If so, it must set 'recalculated
103 // TODO/idea: add a new flag: 'unfinish
104 //
105 // IntersectedOrRecalculated means diff
106 // a) if it is the same curve lenght al
107 // original enpdoint due to the nee
108 // b) if it is at a shorter curve lengt
109 // i.e. as far as it could go, becau
110 // Note: IntersectedOrRecalculated is v
111 // 'true'.
112 // -------------------------------------------
113 // NOTE: implementation taken from G4Propagato
114 //
115 G4bool G4MultiLevelLocator::EstimateIntersecti
116 const G4FieldTrack& CurveStart
117 const G4FieldTrack& CurveEndPo
118 const G4ThreeVector& TrialPoint
119 G4FieldTrack& Intersecte
120 G4bool& recalculat
121 G4double& previousSa
122 G4ThreeVector& previousSf
123 {
124 // Find Intersection Point ( A, B, E ) of t
125 const char* MethodName= "G4MultiLevelLocator
126
127 G4bool found_approximate_intersection = fals
128 G4bool there_is_no_intersection = fals
129
130 G4FieldTrack CurrentA_PointVelocity = Curve
131 G4FieldTrack CurrentB_PointVelocity = Curve
132 G4ThreeVector CurrentE_Point = TrialPoint;
133 G4bool validNormalAtE = false;
134 G4ThreeVector NormalAtEntry;
135
136 G4FieldTrack ApproxIntersecPointV(CurveEndP
137 G4bool validIntersectP= true; // Is
138 G4double NewSafety = 0.0;
139 G4bool last_AF_intersection = false;
140
141 auto integrDriver = GetChordFinderFor()->
142 G4bool driverReIntegrates = integrDriver->D
143
144 G4bool first_section = true;
145 recalculatedEndPoint = false;
146 G4bool restoredFullEndpoint = false;
147
148 unsigned int substep_no = 0;
149
150 // Statistics for substeps
151 static G4ThreadLocal unsigned int max_no_see
152
153 #ifdef G4DEBUG_FIELD
154 unsigned int trigger_substepno_print = 0;
155 const G4double tolerance = 1.0e-8 * CLHEP::m
156 unsigned int biggest_depth = 0;
157 // using kInitialisingCL = G4LocatorChange
158 #endif
159
160 // Log the location, iteration of changes in
161 //------------------------------------------
162 static G4ThreadLocal G4LocatorChangeLogger e
163 endChangeB("EndPointB"), recApproxPoint("A
164 pointH_logger("Trial points 'E': position,
165 unsigned int eventCount = 0;
166
167 if (fCheckMode)
168 {
169 // Clear previous call's data
170 endChangeA.clear();
171 endChangeB.clear();
172 recApproxPoint.clear();
173 pointH_logger.clear();
174
175 // Record the initialisation
176 ++eventCount;
177 endChangeA.AddRecord( G4LocatorChangeRecor
178 eventCount, CurrentA
179 endChangeB.AddRecord( G4LocatorChangeRecor
180 eventCount, CurrentB
181 }
182
183 //------------------------------------------
184 // Algorithm for the case if progress in fo
185 // Process is defined too slow if after N=p
186 // path, it will be only 'fraction_done' of
187 // In this case the remaining length is div
188 // the loop is restarted for each half.
189 // If progress is still too slow, the divis
190 // until 'max_depth'.
191 //------------------------------------------
192
193 const G4int param_substeps = 5; // Test val
194 // of subst
195 const G4double fraction_done = 0.3;
196
197 G4bool Second_half = false; // First half
198
199 // We need to know this for the 'final_secti
200 // real 'final_section' or first half 'final
201 // In algorithm it is considered that the 'S
202 // and it becomes false only if we are in th
203 // depthness or if we are in the first secti
204
205 unsigned int depth = 0; // Depth counts subd
206 ++fNumCalls;
207
208 NormalAtEntry = GetSurfaceNormal(CurrentE_Po
209
210 if (fCheckMode)
211 {
212 pointH_logger.AddRecord( G4LocatorChangeRe
213 substep_no, event
214 G4FieldTrack( Cur
215 0.,
216 #if (G4DEBUG_FIELD>1)
217 G4ThreeVector StartPosition = CurveStartP
218 if( (TrialPoint - StartPosition).mag2() <
219 {
220 ReportImmediateHit( MethodName, StartPo
221 tolerance, fNumCall
222 }
223 #endif
224 }
225
226 // Intermediates Points on the Track = Subdi
227 // Give the initial values to 'InterMedFt'
228 // Important is 'ptrInterMedFT[0]', it saves
229 //
230 *ptrInterMedFT[0] = CurveEndPointVelocity;
231 for ( auto idepth=1; idepth<max_depth+1; ++i
232 {
233 *ptrInterMedFT[idepth] = CurveStartPointVe
234 }
235
236 // Final_section boolean store
237 //
238 G4bool fin_section_depth[max_depth];
239 for (bool & idepth : fin_section_depth)
240 {
241 idepth = true;
242 }
243 // 'SubStartPoint' is needed to calculate th
244 //
245 G4FieldTrack SubStart_PointVelocity = CurveS
246
247 do // Loop checking, 07.10.2016, J.Apostola
248 {
249 unsigned int substep_no_p = 0;
250 G4bool sub_final_section = false; // the s
251 // but f
252 SubStart_PointVelocity = CurrentA_PointVel
253
254 do // Loop checking, 07.10.2016, J.Apostol
255 { // REPEAT param
256 G4ThreeVector Point_A = CurrentA_PointVe
257 G4ThreeVector Point_B = CurrentB_PointVe
258
259 #ifdef G4DEBUG_FIELD
260 const G4double lenA = CurrentA_PointVelo
261 const G4double lenB = CurrentB_PointVelo
262 G4double curv_lenAB = lenB - lenA;
263 G4double distAB = (Point_B - Point_A
264 if( curv_lenAB < distAB * ( 1. - 10.*fiE
265 {
266 G4cerr << "ERROR> (Start) Point A coin
267 << "MLL: iters = " << substep_n
268 G4long op=G4cerr.precision(6);
269 G4cerr << " Difference = " << di
270 << " exceeds limit of relative
271 << " i.e. limit = " << 10 * fi
272 G4cerr.precision(9);
273 G4cerr << " Len A, B = " << len
274 << " Position A: " << Po
275 << " Position B: " << Po
276 G4cerr.precision(op);
277 // G4LocatorChangeRecord::ReportVector
278 // G4cerr<<"EndPoints A(start) and B(e
279 if (fCheckMode)
280 {
281 G4LocatorChangeLogger::ReportEndChan
282 }
283 }
284
285 if( !validIntersectP )
286 {
287 G4ExceptionDescription errmsg;
288 errmsg << "Assertion FAILURE - invalid
289 << substep_no << " call: " << f
290 if (fCheckMode)
291 {
292 G4LocatorChangeRecord::ReportEndChan
293 }
294 G4Exception("G4MultiLevelLocator::Esti
295 JustWarning, errmsg);
296 }
297 #endif
298
299 // F = a point on true AB path close to
300 // (the closest if possible)
301 //
302 ApproxIntersecPointV = GetChordFinderFor
303 ->ApproxCurvePoin
304
305
306
307 // The above method is the key & most
308
309 #ifdef G4DEBUG_FIELD
310 recApproxPoint.push_back(G4LocatorChange
311
312 G4double lenIntsc= ApproxIntersecPointV.
313 G4double checkVsEnd= lenB - lenIntsc;
314
315 if( lenIntsc > lenB )
316 {
317 std::ostringstream errmsg;
318 errmsg.precision(17);
319 G4double ratio = checkVsEnd / lenB;
320 G4double ratioTol = std::fabs(ratio) /
321 errmsg << "Intermediate F point is pas
322 << " l( intersection ) = " << lenInt
323 << " l( endpoint ) = " << lenB
324 errmsg.precision(8);
325 errmsg << " l_end - l_inters = " <<
326 << " / l_end = " << rati
327 << " ratio / tolerance = " << rati
328 if( ratioTol < 1.0 )
329 G4Exception(MethodName, "GeomNav0003
330 else
331 G4Exception(MethodName, "GeomNav0003
332 }
333 #endif
334
335 G4ThreeVector CurrentF_Point= ApproxInte
336
337 // First check whether EF is small - the
338 // Calculate the length and direction of
339 //
340 G4ThreeVector ChordEF_Vector = CurrentF
341
342 G4ThreeVector NewMomentumDir = ApproxIn
343 G4double MomDir_dot_Norm = NewMome
344
345 #ifdef G4DEBUG_FIELD
346 if( fVerboseLevel > 3 )
347 {
348 G4ThreeVector ChordAB = Po
349 G4double ABchord_length = Ch
350 G4double MomDir_dot_ABchord;
351 MomDir_dot_ABchord = (1.0 / ABchord_l
352 * NewMomentumDir.d
353 G4VIntersectionLocator::ReportTrialSt
354 ChordEF_Vector, NewMomentumDir,
355 G4cout << " dot( MomentumDir, ABchord
356 << MomDir_dot_ABchord << G4end
357 }
358 #endif
359 G4bool adequate_angle =
360 ( MomDir_dot_Norm >= 0.0 ) // Can
361 || (! validNormalAtE ); // Inv
362 G4double EF_dist2 = ChordEF_Vector.mag2(
363 if ( ( EF_dist2 <= sqr(fiDeltaIntersecti
364 || ( EF_dist2 <= kCarTolerance*kCarTol
365 {
366 found_approximate_intersection = true;
367
368 // Create the "point" return value
369 //
370 IntersectedOrRecalculatedFT = ApproxIn
371 IntersectedOrRecalculatedFT.SetPositio
372
373 if ( GetAdjustementOfFoundIntersection
374 {
375 // Try to Get Correction of Intersec
376 //
377 G4ThreeVector IP;
378 G4ThreeVector MomentumDir=ApproxInte
379 G4bool goodCorrection = AdjustmentOf
380 CurrentE_P
381 last_AF_in
382 previousSa
383 if ( goodCorrection )
384 {
385 IntersectedOrRecalculatedFT = Appr
386 IntersectedOrRecalculatedFT.SetPos
387 }
388 }
389 // Note: in order to return a point on
390 // we must return E. But it is F
391 // So we must "cheat": we are us
392 // and the velocity at point F !
393 //
394 // This must limit the length we can a
395 }
396 else // E is NOT close enough to the cu
397 {
398 // Check whether any volumes are encou
399 // -----------------------------------
400 // First relocate to restore any Voxel
401 // in the Navigator before calling Com
402 //
403 GetNavigatorFor()->LocateGlobalPointWi
404
405 G4ThreeVector PointG; // Candidate i
406 G4double stepLengthAF;
407 G4bool Intersects_FB = false;
408 G4bool Intersects_AF = IntersectChord(
409
410
411
412
413 last_AF_intersection = Intersects_AF;
414 if( Intersects_AF )
415 {
416 // G is our new Candidate for the in
417 // It replaces "E" and we will see
418 CurrentB_PointVelocity = ApproxInter
419 CurrentE_Point = PointG;
420
421 validIntersectP = true; // 'E' ha
422
423 G4bool validNormalLast;
424 NormalAtEntry = GetSurfaceNormal( P
425 validNormalAtE = validNormalLast;
426
427 // As we move point B, must take car
428 // AF has no intersection to try cur
429 fin_section_depth[depth] = false;
430
431 if (fCheckMode)
432 {
433 ++eventCount;
434 endChangeB.push_back(
435 G4LocatorChangeRecord(G4LocatorC
436 substep_no, event
437 }
438 #ifdef G4VERBOSE
439 if( fVerboseLevel > 3 )
440 {
441 G4cout << "G4PiF::LI> Investigatin
442 << " at s=" << ApproxInters
443 << " on way to full s="
444 << CurveEndPointVelocity.Ge
445 }
446 #endif
447 }
448 else // not Intersects_AF
449 {
450 // In this case:
451 // There is NO intersection of AF wi
452 // We must continue the search in th
453 //
454 GetNavigatorFor()->LocateGlobalPoint
455
456 G4double stepLengthFB;
457 G4ThreeVector PointH;
458
459 // Check whether any volumes are enc
460 // ---------------------------------
461
462 Intersects_FB = IntersectChord( Curr
463 NewS
464 prev
465 step
466 Poin
467 if( Intersects_FB )
468 {
469 // There is an intersection of FB
470 // H <- First Intersection of Chor
471
472 // H is our new Candidate for the
473 // It replaces "E" and we will re
474 // it is a good enough approximate
475
476 // Note that F must be in volume v
477 // (otherwise AF would meet a volu
478 // A <- F
479 // E <- H
480 //
481 CurrentA_PointVelocity = ApproxInt
482 CurrentE_Point = PointH;
483
484 validIntersectP = true; // 'E'
485
486 G4bool validNormalH;
487 NormalAtEntry = GetSurfaceNormal(
488 validNormalAtE = validNormalH;
489
490 if (fCheckMode)
491 {
492 ++eventCount;
493 endChangeA.push_back(
494 G4LocatorChangeRecord(G4Locat
495 substep_no, event
496 G4FieldTrack intersectH_pn('0');
497
498 intersectH_pn.SetPosition( Point
499 intersectH_pn.SetMomentum( Norma
500 pointH_logger.AddRecord(G4Locato
501 substep_no
502 }
503 }
504 else // not Intersects_FB
505 {
506 validIntersectP = false; // Int
507 if( fin_section_depth[depth] )
508 {
509 // If B is the original endpoint
510 // volume(s) intersected the ori
511 // smaller chords we have used.
512 // The value of 'IntersectedOrRe
513 // likely not valid
514
515 // Check on real final_section o
516 //
517 if( ((Second_half)&&(depth==0))
518 {
519 there_is_no_intersection = tru
520 }
521 else
522 {
523 // end of subsection, not real
524 // exit from the and go to the
525 substep_no_p = param_substeps+
526
527 // but 'Second_half' is still
528 // the 'CurrentE_point' for th
529 Second_half = true;
530 sub_final_section = true;
531 }
532 }
533 else
534 {
535 CurrentA_PointVelocity = Current
536 CurrentB_PointVelocity = (depth=
537
538 SubStart_PointVelocity = Current
539 restoredFullEndpoint = true;
540
541 validIntersectP = false; // 'E'
542
543 if (fCheckMode)
544 {
545 ++eventCount;
546 endChangeA.push_back(
547 G4LocatorChangeRecord(
548 G4LocatorChangeRecord::kNo
549 substep_no, event
550 endChangeB.push_back(
551 G4LocatorChangeRecord (
552 G4LocatorChangeRecord::kNo
553 substep_no, event
554 }
555 }
556 } // Endif (Intersects_FB)
557 } // Endif (Intersects_AF)
558
559 G4int errorEndPt = 0; // Default: no e
560
561 G4bool recalculatedB= false;
562 if( driverReIntegrates )
563 {
564 G4FieldTrack RevisedB_FT = CurrentB_
565 recalculatedB= CheckAndReEstimateEnd
566
567
568
569 if( recalculatedB )
570 {
571 CurrentB_PointVelocity = RevisedB_
572 // Do not invalidate intersection
573 //
574 // The best course would be to inv
575 // BUT if we invalidate it, we mus
576 // validIntersectP = false; //
577
578 if ( (fin_section_depth[depth])
579 &&( first_section || ((Second_h
580 {
581 recalculatedEndPoint = true;
582 IntersectedOrRecalculatedFT = Re
583 // So that we can return it, if
584 }
585 // else
586 // Move forward the other points
587 // - or better flag it, so that
588 // [ Implementation: a counter
589 // => avoids extra work]
590 }
591 if (fCheckMode)
592 {
593 ++eventCount;
594 endChangeB.push_back(
595 G4LocatorChangeRecord( G4Locator
596 substep_n
597 }
598 }
599 else
600 {
601 if( CurrentB_PointVelocity.GetCurveL
602 < CurrentA_PointVelocity.GetCurveL
603 {
604 errorEndPt = 2;
605 }
606 }
607
608 if( errorEndPt > 1 ) // errorEndPt =
609 {
610 std::ostringstream errmsg;
611
612 ReportReversedPoints(errmsg,
613 CurveStartPoint
614 NewSafety, fiEp
615 CurrentA_PointV
616 SubStart_PointV
617 ApproxIntersecP
618
619 if (fCheckMode)
620 {
621 G4LocatorChangeRecord::ReportEndCh
622 }
623
624 errmsg << G4endl << " * Location: "
625 << "- After EndIf(Intersects_
626 errmsg << " * Bool flags: Recalcula
627 << " Intersects_AF = " << I
628 << " Intersects_FB = " << I
629 errmsg << " * Number of calls to MLL
630 G4Exception(MethodName, "GeomNav0003
631 }
632 if( restoredFullEndpoint )
633 {
634 fin_section_depth[depth] = restoredF
635 restoredFullEndpoint = false;
636 }
637 } // EndIf ( E is close enough to the cu
638 // tests ChordAF_Vector.mag() <= maxim
639
640 #ifdef G4DEBUG_FIELD
641 if( trigger_substepno_print == 0)
642 {
643 trigger_substepno_print= fWarnSteps -
644 }
645
646 if( substep_no >= trigger_substepno_prin
647 {
648 G4cout << "Difficulty in converging in
649 << " Substep no = " << subst
650 if( substep_no == trigger_substepno_pr
651 {
652 G4cout << " Start: ";
653 printStatus( CurveStartPointVelocity
654 -1.0, NewSafety, 0 );
655 if( fCheckMode ) {
656 G4LocatorChangeRecord::ReportEndCh
657 } else {
658 G4cout << " ** For more informatio
659 << "-- (it saves and can ou
660 }
661 }
662 G4cout << " Point A: ";
663 printStatus( CurrentA_PointVelocity, C
664 -1.0, NewSafety, substep_
665 G4cout << " Point B: ";
666 printStatus( CurrentA_PointVelocity, C
667 -1.0, NewSafety, substep_
668 }
669 #endif
670 ++substep_no;
671 ++substep_no_p;
672
673 } while ( ( ! found_approximate_intersect
674 && ( ! there_is_no_intersection )
675 && validIntersectP // New c
676 && ( substep_no_p <= param_substep
677
678
679 if( (!found_approximate_intersection) && (
680 {
681 G4double did_len = std::abs( CurrentA_Po
682 - SubStart_PointVelocit
683 G4double all_len = std::abs( CurrentB_Po
684 - SubStart_PointVelocit
685
686 G4double distAB = -1;
687 //
688 // Is progress is too slow, and is it po
689 // If so, then *halve the step*
690 // ==============
691 if( (did_len < fraction_done*all_len)
692 && (depth<max_depth) && (!sub_final_s
693 {
694 #ifdef G4DEBUG_FIELD
695 static G4ThreadLocal unsigned int numS
696 biggest_depth = std::max(depth, bigges
697 ++numSplits;
698 #endif
699 Second_half = false;
700 ++depth;
701 first_section = false;
702
703 G4double Sub_len = (all_len-did_len)/(
704 G4FieldTrack midPoint = CurrentA_Point
705 G4bool fullAdvance=
706 integrDriver->AccurateAdvance(midPo
707
708 ++fNumAdvanceTrials;
709 if( fullAdvance ) { ++fNumAdvanceFull
710
711 G4double lenAchieved=
712 midPoint.GetCurveLength()-CurrentA_
713
714 const G4double adequateFraction = (1.0
715 G4bool goodAdvance = (lenAchieved >= a
716 if ( goodAdvance ) { ++fNumAdvanceGoo
717
718 #ifdef G4DEBUG_FIELD
719 else // !goodAdvance
720 {
721 G4cout << "MLL> AdvanceChordLimited
722 << " total length achieved
723 << Sub_len << " fraction= "
724 if (Sub_len != 0.0 ) { G4cout << le
725 else { G4cout << "D
726 G4cout << " Good-enough fraction =
727 G4cout << " Number of call to mll
728 << " iteration " << subste
729 << " inner = " << substep_
730 G4cout << " Epsilon of integration
731 G4cout << " State at start is
732 << G4endl
733 << " at end (midpoint
734 G4cout << " Particle mass = " << m
735
736 G4EquationOfMotion *equation = inte
737 ReportFieldValue( CurrentA_PointVel
738 ReportFieldValue( midPoint, "midPoi
739 G4cout << " Original Start = "
740 << CurveStartPointVelocity <
741 G4cout << " Original End = "
742 << CurveEndPointVelocity <<
743 G4cout << " Original TrialPoint =
744 << TrialPoint << G4endl;
745 G4cout << " (this is STRICT mode)
746 << " num Splits= " << numSp
747 G4cout << G4endl;
748 }
749 #endif
750
751 *ptrInterMedFT[depth] = midPoint;
752 CurrentB_PointVelocity = midPoint;
753
754 if (fCheckMode)
755 {
756 ++eventCount;
757 endChangeB.push_back(
758 G4LocatorChangeRecord( G4LocatorCh
759 substep_no,
760 }
761
762 // Adjust 'SubStartPoint' to calculate
763 //
764 SubStart_PointVelocity = CurrentA_Poin
765
766 // Find new trial intersection point n
767 //
768 G4ThreeVector Point_A = CurrentA_Point
769 G4ThreeVector Point_B = CurrentB_Point
770
771 G4ThreeVector PointGe;
772 GetNavigatorFor()->LocateGlobalPointWi
773 G4bool Intersects_AB = IntersectChord(
774
775
776
777 if( Intersects_AB )
778 {
779 last_AF_intersection = Intersects_AB
780 CurrentE_Point = PointGe;
781 fin_section_depth[depth] = true;
782
783 validIntersectP = true; // 'E' has
784
785 G4bool validNormalAB;
786 NormalAtEntry = GetSurfaceNormal( P
787 validNormalAtE = validNormalAB;
788 }
789 else
790 {
791 // No intersection found for first p
792 // (CurrentA,InterMedPoint[depth]).
793 //
794 Second_half = true;
795
796 validIntersectP= false; // No new
797 }
798 } // if did_len
799
800 G4bool unfinished = Second_half;
801 while ( unfinished && (depth>0) ) // Lo
802 {
803 // Second part of curve (InterMed[dept
804 // On the depth-1 level normally we ar
805
806 // Find new trial intersection point
807 //
808 SubStart_PointVelocity = *ptrInterMedF
809 CurrentA_PointVelocity = *ptrInterMedF
810 CurrentB_PointVelocity = *ptrInterMedF
811
812 if (fCheckMode)
813 {
814 ++eventCount;
815 G4LocatorChangeRecord chngPop_a( G4L
816 substep_no, even
817 endChangeA.push_back( chngPop_a );
818 G4LocatorChangeRecord chngPop_b( G4L
819 substep_no, even
820 endChangeB.push_back( chngPop_b );
821 }
822
823 // Ensure that the new endpoints are n
824 // than on the curve due to different
825 //
826 G4int errorEndPt = -1;
827 G4bool recalculatedB= false;
828 if( driverReIntegrates )
829 {
830 // Ensure that the new endpoints are
831 // than on the curve due to differen
832 //
833 G4FieldTrack RevisedEndPointFT = Cur
834 recalculatedB =
835 CheckAndReEstimateEndpoint( Curre
836 Curre
837 Revis
838 error
839 if( recalculatedB )
840 {
841 CurrentB_PointVelocity = RevisedEn
842
843 if ( depth == 1 )
844 {
845 recalculatedEndPoint = true;
846 IntersectedOrRecalculatedFT = R
847 // So that we can return it, if
848 }
849 }
850 else
851 {
852 if( CurrentB_PointVelocity.GetCurv
853 < CurrentA_PointVelocity.GetCurv
854 {
855 errorEndPt = 2;
856 }
857 }
858
859 if (fCheckMode)
860 {
861 ++eventCount;
862 endChangeB.push_back(
863 G4LocatorChangeRecord(G4LocatorC
864 substep_no
865 }
866 }
867 if( errorEndPt > 1 ) // errorEndPt =
868 {
869 std::ostringstream errmsg;
870 ReportReversedPoints(errmsg,
871 CurveStartPointVelocity, C
872 NewSafety, fiEpsilonStep,
873 CurrentA_PointVelocity, Cu
874 SubStart_PointVelocity, Cu
875 ApproxIntersecPointV, subs
876 errmsg << " * Location: " << Method
877 errmsg << " * Recalculated = " << re
878 G4Exception(MethodName, "GeomNav0003
879 }
880
881 G4ThreeVector Point_A = CurrentA_Point
882 G4ThreeVector Point_B = CurrentB_Point
883 G4ThreeVector PointGi;
884 GetNavigatorFor()->LocateGlobalPointWi
885 G4bool Intersects_AB = IntersectChord(
886
887
888
889 if( Intersects_AB )
890 {
891 last_AF_intersection = Intersects_AB
892 CurrentE_Point = PointGi;
893
894 validIntersectP = true; // 'E' has
895 NormalAtEntry = GetSurfaceNormal( P
896 }
897 else
898 {
899 validIntersectP = false; // No new
900 if( depth == 1)
901 {
902 there_is_no_intersection = true;
903 }
904 }
905 depth--;
906 fin_section_depth[depth] = true;
907 unfinished = !validIntersectP;
908 }
909 #ifdef G4DEBUG_FIELD
910 if( ! ( validIntersectP || there_is_no_i
911 {
912 // What happened ??
913 G4cout << "MLL - WARNING Potential FA
914 << G4endl
915 << " Depth = " << depth << G4e
916 << " Num Substeps= " << subste
917 G4cout << " Found intersection= " <<
918 << G4endl;
919 G4cout << " Progress report: -- " <<
920 ReportProgress(G4cout,
921 CurveStartPointVelocit
922 substep_no, CurrentA_P
923 CurrentB_PointVelocity
924 NewSafety, depth );
925 G4cout << G4endl;
926 }
927 #endif
928 } // if(!found_aproximate_intersection)
929
930 assert( validIntersectP || there_is_no_int
931 || found_approxima
932
933 } while ( ( ! found_approximate_intersection
934 && ( ! there_is_no_intersection )
935 && ( substep_no <= fMaxSteps) ); /
936
937 if( substep_no > max_no_seen )
938 {
939 max_no_seen = substep_no;
940 #ifdef G4DEBUG_FIELD
941 if( max_no_seen > fWarnSteps )
942 {
943 trigger_substepno_print = max_no_seen-20
944 }
945 #endif
946 }
947
948 if( !there_is_no_intersection && !found_appr
949 {
950 if( substep_no >= fMaxSteps)
951 {
952 // Since we cannot go further (yet), w
953
954 recalculatedEndPoint = true;
955 IntersectedOrRecalculatedFT = CurrentA
956 found_approximate_intersection = false
957
958 std::ostringstream message;
959 message << G4endl;
960 message << "Convergence is requiring t
961 << substep_no << " ( limit =
962 << G4endl
963 << " Abandoning effort to int
964 message << " Number of calls to MLL
965 message << " iteration = " << substep
966
967 message.precision( 12 );
968 G4double done_len = CurrentA_PointVelo
969 G4double full_len = CurveEndPointVeloc
970 message << " Undertaken only le
971 << " out of " << full_len << "
972 << " Remaining length =
973
974 message << " Start and end-point o
975 printStatus( CurveStartPointVelocity,
976 -1.0, NewSafety, 0, m
977 message << " Start and end-poin
978 printStatus( CurrentA_PointVelocity, C
979 -1.0, NewSafety, substep_
980 printStatus( CurrentA_PointVelocity, C
981 -1.0, NewSafety, substep_
982
983 G4Exception(MethodName, "GeomNav0003",
984 }
985 else if( substep_no >= fWarnSteps)
986 {
987 std::ostringstream message;
988 message << "Many substeps while trying
989 << G4endl
990 << " Undertaken lengt
991 << CurrentB_PointVelocity.GetC
992 << " - Needed: " << substep_n
993 << " Warning number =
994 << " and maximum substeps = "
995 G4Exception(MethodName, "GeomNav1002",
996 }
997 }
998
999 return (!there_is_no_intersection) && found
1000 // Success or failure
1001 }
1002
1003 void G4MultiLevelLocator::ReportStatistics()
1004 {
1005 G4cout << " Number of calls = " << fNumCal
1006 G4cout << " Number of split level ('advanc
1007 << fNumAdvanceTrials << G4endl;
1008 G4cout << " Number of full advances:
1009 << fNumAdvanceGood << G4endl;
1010 G4cout << " Number of good advances:
1011 << fNumAdvanceFull << G4endl;
1012 }
1013
1014 void G4MultiLevelLocator::ReportFieldValue( c
1015 c
1016 c
1017 {
1018 enum { maxNumFieldComp = 24 };
1019
1020 G4ThreeVector position = locationPV.GetPos
1021 G4double startPoint[4] = { position.x(), p
1022 locationPV.GetL
1023 G4double FieldVec[maxNumFieldComp]; // 24
1024 for (double & i : FieldVec)
1025 {
1026 i = 0.0;
1027 }
1028 equation->GetFieldValue( startPoint, Field
1029 G4cout << " B-field value (" << nameLoc <
1030 << FieldVec[0] << " " << FieldVec[1
1031 G4double Emag2= G4ThreeVector( FieldVec[3]
1032 FieldVec[4]
1033 FieldVec[5]
1034 if( Emag2 > 0.0 )
1035 {
1036 G4cout << " Electric = " << FieldVec[3]
1037 << FieldVec[4]
1038 << FieldVec[5]
1039 }
1040 return;
1041 }
1042