123{
124
125 const char* MethodName= "G4MultiLevelLocator::EstimateIntersectionPoint()";
126
127 G4bool found_approximate_intersection =
false;
128 G4bool there_is_no_intersection =
false;
129
130 G4FieldTrack CurrentA_PointVelocity = CurveStartPointVelocity;
131 G4FieldTrack CurrentB_PointVelocity = CurveEndPointVelocity;
133 G4bool validNormalAtE =
false;
135
136 G4FieldTrack ApproxIntersecPointV(CurveEndPointVelocity);
137 G4bool validIntersectP=
true;
139 G4bool last_AF_intersection =
false;
140
143
144 G4bool first_section =
true;
145 recalculatedEndPoint = false;
146 G4bool restoredFullEndpoint =
false;
147
148 unsigned int substep_no = 0;
149
150
152
153#ifdef G4DEBUG_FIELD
154 unsigned int trigger_substepno_print = 0;
156 unsigned int biggest_depth = 0;
157
158#endif
159
160
161
163 endChangeB("EndPointB"), recApproxPoint("ApproxPoint"),
164 pointH_logger("Trial points 'E': position, normal");
165 unsigned int eventCount = 0;
166
168 {
169
170 endChangeA.clear();
171 endChangeB.clear();
172 recApproxPoint.clear();
173 pointH_logger.clear();
174
175
176 ++eventCount;
178 eventCount, CurrentA_PointVelocity );
180 eventCount, CurrentB_PointVelocity );
181 }
182
183
184
185
186
187
188
189
190
191
192
193 const G4int param_substeps = 5;
194
196
197 G4bool Second_half =
false;
198
199
200
201
202
203
204
205 unsigned int depth = 0;
207
209
211 {
213 substep_no, eventCount,
216 #if (G4DEBUG_FIELD>1)
217 G4ThreeVector StartPosition = CurveStartPointVelocity.GetPosition();
218 if( (TrialPoint - StartPosition).mag2() <
sqr(tolerance))
219 {
222 }
223 #endif
224 }
225
226
227
228
229
231 for (
auto idepth=1; idepth<
max_depth+1; ++idepth )
232 {
234 }
235
236
237
239 for (
auto idepth=0; idepth<
max_depth; ++idepth )
240 {
241 fin_section_depth[idepth] = true;
242 }
243
244
245 G4FieldTrack SubStart_PointVelocity = CurveStartPointVelocity;
246
247 do
248 {
249 unsigned int substep_no_p = 0;
250 G4bool sub_final_section =
false;
251
252 SubStart_PointVelocity = CurrentA_PointVelocity;
253
254 do
255 {
256
257#ifdef G4DEBUG_FIELD
260 {
261 G4cerr <<
"ERROR> (Start) Point A coincides with or has gone past (end) point B"
262 <<
"MLL: iters = " << substep_no <<
G4endl;
263
264
266 }
267#endif
270
271
272
273
276 CurrentB_PointVelocity,
277 CurrentE_Point,
279
280
281#ifdef G4DEBUG_FIELD
283 substep_no, eventCount, ApproxIntersecPointV ) );
286 G4double checkVsEnd= lenB - lenIntsc;
287
288 if( lenIntsc > lenB )
289 {
290 std::ostringstream errmsg;
291 errmsg.precision(17);
293 G4double ratioTol = std::fabs(ratio) / tolerance;
294 errmsg <<
"Intermediate F point is past end B point" <<
G4endl
295 <<
" l( intersection ) = " << lenIntsc <<
G4endl
296 <<
" l( endpoint ) = " << lenB <<
G4endl;
297 errmsg.precision(8);
298 errmsg <<
" l_end - l_inters = " << checkVsEnd <<
G4endl
299 <<
" / l_end = " << ratio <<
G4endl
300 <<
" ratio / tolerance = " << ratioTol <<
G4endl;
301 if( ratioTol < 1.0 )
303 else
305 }
306#endif
307
308 G4ThreeVector CurrentF_Point= ApproxIntersecPointV.GetPosition();
309
310
311
312
313 G4ThreeVector ChordEF_Vector = CurrentF_Point - CurrentE_Point;
314
315 G4ThreeVector NewMomentumDir = ApproxIntersecPointV.GetMomentumDir();
316 G4double MomDir_dot_Norm = NewMomentumDir.
dot( NormalAtEntry );
317
318#ifdef G4DEBUG_FIELD
320 {
324 MomDir_dot_ABchord = (1.0 / ABchord_length)
325 * NewMomentumDir.
dot( ChordAB );
327 ChordEF_Vector, NewMomentumDir, NormalAtEntry, validNormalAtE );
328 G4cout <<
" dot( MomentumDir, ABchord_unit ) = "
329 << MomDir_dot_ABchord <<
G4endl;
330 }
331#endif
333 ( MomDir_dot_Norm >= 0.0 )
334 || (! validNormalAtE );
338 {
339 found_approximate_intersection = true;
340
341
342
343 IntersectedOrRecalculatedFT = ApproxIntersecPointV;
344 IntersectedOrRecalculatedFT.SetPosition( CurrentE_Point );
345
347 {
348
349
351 G4ThreeVector MomentumDir=ApproxIntersecPointV.GetMomentumDirection();
353 CurrentE_Point, CurrentF_Point, MomentumDir,
354 last_AF_intersection, IP, NewSafety,
355 previousSafety, previousSftOrigin );
356 if ( goodCorrection )
357 {
358 IntersectedOrRecalculatedFT = ApproxIntersecPointV;
359 IntersectedOrRecalculatedFT.SetPosition(IP);
360 }
361 }
362
363
364
365
366
367
368 }
369 else
370 {
371
372
373
374
375
377
380 G4bool Intersects_FB =
false;
382 NewSafety, previousSafety,
383 previousSftOrigin,
384 stepLengthAF,
385 PointG );
386 last_AF_intersection = Intersects_AF;
387 if( Intersects_AF )
388 {
389
390
391 CurrentB_PointVelocity = ApproxIntersecPointV;
392 CurrentE_Point = PointG;
393
394 validIntersectP = true;
395
398 validNormalAtE = validNormalLast;
399
400
401
402 fin_section_depth[depth] = false;
403
405 {
406 ++eventCount;
407 endChangeB.push_back(
409 substep_no, eventCount, CurrentB_PointVelocity) );
410 }
411#ifdef G4VERBOSE
413 {
414 G4cout <<
"G4PiF::LI> Investigating intermediate point"
415 << " at s=" << ApproxIntersecPointV.GetCurveLength()
416 << " on way to full s="
417 << CurveEndPointVelocity.GetCurveLength() <<
G4endl;
418 }
419#endif
420 }
421 else
422 {
423
424
425
426
428
431
432
433
434
436 NewSafety, previousSafety,
437 previousSftOrigin,
438 stepLengthFB,
439 PointH );
440 if( Intersects_FB )
441 {
442
443
444
445
446
447
448
449
450
451
452
453
454 CurrentA_PointVelocity = ApproxIntersecPointV;
455 CurrentE_Point = PointH;
456
457 validIntersectP = true;
458
461 validNormalAtE = validNormalH;
462
464 {
465 ++eventCount;
466 endChangeA.push_back(
468 substep_no, eventCount, CurrentA_PointVelocity) );
470
471 intersectH_pn.SetPosition( PointH );
472 intersectH_pn.SetMomentum( NormalAtEntry );
474 substep_no, eventCount, intersectH_pn );
475 }
476 }
477 else
478 {
479 if( fin_section_depth[depth] )
480 {
481
482
483
484
485
486
487
488
489 if( ((Second_half)&&(depth==0)) || (first_section) )
490 {
491 there_is_no_intersection = true;
492 }
493 else
494 {
495
496
497 substep_no_p = param_substeps+2;
498
499
500
501 Second_half = true;
502 sub_final_section = true;
503 }
504 }
505 else
506 {
507 CurrentA_PointVelocity = CurrentB_PointVelocity;
508 CurrentB_PointVelocity = (depth==0) ? CurveEndPointVelocity
510 SubStart_PointVelocity = CurrentA_PointVelocity;
511 restoredFullEndpoint = true;
512
513 validIntersectP = false;
514
516 {
517 ++eventCount;
518 endChangeA.push_back(
521 substep_no, eventCount, CurrentA_PointVelocity) );
522 endChangeB.push_back(
525 substep_no, eventCount, CurrentB_PointVelocity) );
526 }
527 }
528 }
529 }
530
531 G4int errorEndPt = 0;
532
533 G4bool recalculatedB=
false;
534 if( driverReIntegrates )
535 {
538 CurrentB_PointVelocity,
539 RevisedB_FT,
540 errorEndPt );
541 if( recalculatedB )
542 {
543 CurrentB_PointVelocity = RevisedB_FT;
544
545
546
547
548
549
550 if ( (fin_section_depth[depth])
551 &&( first_section || ((Second_half)&&(depth==0)) ) )
552 {
553 recalculatedEndPoint = true;
554 IntersectedOrRecalculatedFT = RevisedB_FT;
555
556 }
557
558
559
560
561
562 }
563
564
565
566
567
569 {
570 ++eventCount;
571 endChangeB.push_back(
573 substep_no, eventCount, RevisedB_FT ) );
574 }
575 }
576 else
577 {
579 errorEndPt = 2;
580 }
581
582 if( errorEndPt > 1 )
583 {
584 std::ostringstream errmsg;
586 CurveStartPointVelocity, CurveEndPointVelocity,
588 CurrentA_PointVelocity, CurrentB_PointVelocity,
589 SubStart_PointVelocity, CurrentE_Point,
590 ApproxIntersecPointV, substep_no, substep_no_p, depth);
591 errmsg <<
G4endl <<
" * Location: " << MethodName
592 <<
"- After EndIf(Intersects_AF)" <<
G4endl;
593 errmsg << " * Bool flags: Recalculated = " << recalculatedB
594 << " Intersects_AF = " << Intersects_AF
595 <<
" Intersects_FB = " << Intersects_FB <<
G4endl;
598 }
599 if( restoredFullEndpoint )
600 {
601 fin_section_depth[depth] = restoredFullEndpoint;
602 restoredFullEndpoint = false;
603 }
604 }
605
606
607#ifdef G4DEBUG_FIELD
608 if( trigger_substepno_print == 0)
609 {
611 }
612
613 if( substep_no >= trigger_substepno_print )
614 {
615 G4cout <<
"Difficulty in converging in " << MethodName
616 <<
" Substep no = " << substep_no <<
G4endl;
617 if( substep_no == trigger_substepno_print )
618 {
620 printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
621 -1.0, NewSafety, 0 );
622
624 G4cout <<
"endPoints A (start) and B (end): combined changes of AB intervals" <<
G4endl;
626 }
628 printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
629 -1.0, NewSafety, substep_no-1 );
631 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
632 -1.0, NewSafety, substep_no );
633 }
634#endif
635 ++substep_no;
636 ++substep_no_p;
637
638 } while ( ( ! found_approximate_intersection )
639 && ( ! there_is_no_intersection )
640 && ( substep_no_p <= param_substeps) );
641
642
643 if( (!found_approximate_intersection) && (!there_is_no_intersection) )
644 {
649
651
652
653
654
655 if( (did_len < fraction_done*all_len)
656 && (depth<
max_depth) && (!sub_final_section) )
657 {
658#ifdef G4DEBUG_FIELD
660 biggest_depth =
std::max(depth, biggest_depth);
661 ++numSplits;
662#endif
663 Second_half = false;
664 ++depth;
665 first_section = false;
666
667 G4double Sub_len = (all_len-did_len)/(2.);
670 integrDriver->AccurateAdvance(midPoint, Sub_len,
fiEpsilonStep);
671
674
677
679 G4bool goodAdvance = (lenAchieved >= adequateFraction * Sub_len);
681
682#ifdef G4DEBUG_FIELD
683 else
684 {
685 G4cout <<
"MLL> AdvanceChordLimited not full at depth=" << depth
686 << " total length achieved = " << lenAchieved << " of "
687 << Sub_len << " fraction= ";
688 if (Sub_len != 0.0 ) {
G4cout << lenAchieved / Sub_len; }
689 else {
G4cout <<
"DivByZero"; }
690 G4cout <<
" Good-enough fraction = " << adequateFraction;
692 << " iteration " << substep_no
693 <<
" inner = " << substep_no_p <<
G4endl;
695 G4cout <<
" State at start is = " << CurrentA_PointVelocity
697 <<
" at end (midpoint)= " << midPoint <<
G4endl;
699
703 G4cout <<
" Original Start = "
704 << CurveStartPointVelocity <<
G4endl;
705 G4cout <<
" Original End = "
706 << CurveEndPointVelocity <<
G4endl;
707 G4cout <<
" Original TrialPoint = "
709 G4cout <<
" (this is STRICT mode) "
710 << " num Splits= " << numSplits;
712 }
713#endif
714
716 CurrentB_PointVelocity = midPoint;
717
719 {
720 ++eventCount;
721 endChangeB.push_back(
723 substep_no, eventCount, midPoint) );
724 }
725
726
727
728 SubStart_PointVelocity = CurrentA_PointVelocity;
729
730
731
734
738 NewSafety, previousSafety,
739 previousSftOrigin, distAB,
740 PointGe);
741 if( Intersects_AB )
742 {
743 last_AF_intersection = Intersects_AB;
744 CurrentE_Point = PointGe;
745 fin_section_depth[depth] = true;
746
747 validIntersectP = true;
748
751 validNormalAtE = validNormalAB;
752 }
753 else
754 {
755
756
757
758 Second_half = true;
759
760 validIntersectP= false;
761 }
762 }
763
764 unsigned int levelPops = 0;
765
766 G4bool unfinished = Second_half;
767 while ( unfinished && (depth>0) )
768 {
769
770
771
772 ++levelPops;
773
774
775
779
781 {
782 ++eventCount;
784 substep_no, eventCount, CurrentA_PointVelocity);
785 endChangeA.push_back( chngPop_a );
787 substep_no, eventCount, CurrentB_PointVelocity);
788 endChangeB.push_back( chngPop_b );
789 }
790
791
792
793
794 G4int errorEndPt = -1;
795 G4bool recalculatedB=
false;
796 if( driverReIntegrates )
797 {
798
799
800
801 G4FieldTrack RevisedEndPointFT = CurrentB_PointVelocity;
802 recalculatedB =
804 CurrentB_PointVelocity,
805 RevisedEndPointFT,
806 errorEndPt );
807 if( recalculatedB )
808 {
809 CurrentB_PointVelocity = RevisedEndPointFT;
810
811 if ( depth == 1 )
812 {
813 recalculatedEndPoint = true;
814 IntersectedOrRecalculatedFT = RevisedEndPointFT;
815
816 }
817 }
818 else
819 {
821 errorEndPt = 2;
822 }
823
825 {
826 ++eventCount;
827 endChangeB.push_back(
829 substep_no, eventCount, RevisedEndPointFT));
830 }
831 }
832 if( errorEndPt > 1 )
833 {
834 std::ostringstream errmsg;
836 CurveStartPointVelocity, CurveEndPointVelocity,
838 CurrentA_PointVelocity, CurrentB_PointVelocity,
839 SubStart_PointVelocity, CurrentE_Point,
840 ApproxIntersecPointV, substep_no, substep_no_p, depth);
841 errmsg <<
" * Location: " << MethodName <<
"- Second-Half" <<
G4endl;
842 errmsg <<
" * Recalculated = " << recalculatedB <<
G4endl;
844 }
845
851 previousSafety,
852 previousSftOrigin, distAB,
853 PointGi);
854 if( Intersects_AB )
855 {
856 last_AF_intersection = Intersects_AB;
857 CurrentE_Point = PointGi;
858
859 validIntersectP = true;
861 }
862 else
863 {
864 validIntersectP = false;
865 if( depth == 1)
866 {
867 there_is_no_intersection = true;
868 }
869 }
870 depth--;
871 fin_section_depth[depth] = true;
872 unfinished = !validIntersectP;
873 }
874#ifdef G4DEBUG_FIELD
875 if( ! ( validIntersectP || there_is_no_intersection ) )
876 {
877
878 G4cout <<
"MLL - WARNING Potential FAILURE: Conditions not met!"
880 <<
" Depth = " << depth <<
G4endl
881 << " Levels popped = " << levelPops
882 <<
" Num Substeps= " << substep_no <<
G4endl;
883 G4cout <<
" Found intersection= " << found_approximate_intersection
887 CurveStartPointVelocity, CurveEndPointVelocity,
888 substep_no, CurrentA_PointVelocity,
889 CurrentB_PointVelocity,
890 NewSafety, depth );
892 }
893#endif
894 }
895
896 assert( validIntersectP || there_is_no_intersection
897 || found_approximate_intersection);
898
899 } while ( ( ! found_approximate_intersection )
900 && ( ! there_is_no_intersection )
902
903 if( substep_no > max_no_seen )
904 {
905 max_no_seen = substep_no;
906#ifdef G4DEBUG_FIELD
908 {
909 trigger_substepno_print = max_no_seen-20;
910 }
911#endif
912 }
913
914 if( !there_is_no_intersection && !found_approximate_intersection )
915 {
917 {
918
919
920 recalculatedEndPoint = true;
921 IntersectedOrRecalculatedFT = CurrentA_PointVelocity;
922 found_approximate_intersection = false;
923
924 std::ostringstream message;
926 message << "Convergence is requiring too many substeps: "
927 << substep_no <<
" ( limit = "<<
fMaxSteps <<
")"
929 <<
" Abandoning effort to intersect. " <<
G4endl <<
G4endl;
930 message <<
" Number of calls to MLL: " <<
fNumCalls;
931 message <<
" iteration = " << substep_no <<
G4endl <<
G4endl;
932
933 message.precision( 12 );
935 G4double full_len = CurveEndPointVelocity.GetCurveLength();
936 message << " Undertaken only length: " << done_len
937 <<
" out of " << full_len <<
" required." <<
G4endl
938 << " Remaining length = " << full_len - done_len;
939
940 message <<
" Start and end-point of requested Step:" <<
G4endl;
941 printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
942 -1.0, NewSafety, 0, message, -1 );
943 message <<
" Start and end-point of current Sub-Step:" <<
G4endl;
944 printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
945 -1.0, NewSafety, substep_no-1, message, -1 );
946 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
947 -1.0, NewSafety, substep_no, message, -1 );
948
950 }
952 {
953 std::ostringstream message;
954 message << "Many substeps while trying to locate intersection."
956 << " Undertaken length: "
958 <<
" - Needed: " << substep_no <<
" substeps." <<
G4endl
960 <<
" and maximum substeps = " <<
fMaxSteps;
962 }
963 }
964
965 return (!there_is_no_intersection) && found_approximate_intersection;
966
967}
CLHEP::Hep3Vector G4ThreeVector
G4GLOB_DLL std::ostream G4cerr
G4GLOB_DLL std::ostream G4cout
G4FieldTrack ApproxCurvePointV(const G4FieldTrack &curveAPointVelocity, const G4FieldTrack &curveBPointVelocity, const G4ThreeVector ¤tEPoint, G4double epsStep)
G4VIntegrationDriver * GetIntegrationDriver()
G4double GetRestMass() const
static std::ostream & ReportEndChanges(std::ostream &os, const G4LocatorChangeLogger &startA, const G4LocatorChangeLogger &endB)
static std::ostream & ReportEndChanges(std::ostream &os, const std::vector< G4LocatorChangeRecord > &startA, const std::vector< G4LocatorChangeRecord > &endB)
void ReportFieldValue(const G4FieldTrack &locationPV, const char *nameLoc, const G4EquationOfMotion *equation)
unsigned long int fNumCalls
unsigned long int fNumAdvanceTrials
unsigned long int fNumAdvanceGood
unsigned long int fNumAdvanceFull
virtual G4bool DoesReIntegrate() const =0
G4double fiDeltaIntersection
G4ChordFinder * GetChordFinderFor()
G4ThreeVector GetSurfaceNormal(const G4ThreeVector &CurrentInt_Point, G4bool &validNormal)
void ReportTrialStep(G4int step_no, const G4ThreeVector &ChordAB_v, const G4ThreeVector &ChordEF_v, const G4ThreeVector &NewMomentumDir, const G4ThreeVector &NormalAtEntry, G4bool validNormal)
G4bool CheckAndReEstimateEndpoint(const G4FieldTrack &CurrentStartA, const G4FieldTrack &EstimatedEndB, G4FieldTrack &RevisedEndPoint, G4int &errorCode)
void ReportProgress(std::ostream &oss, const G4FieldTrack &StartPointVel, const G4FieldTrack &EndPointVel, G4int substep_no, const G4FieldTrack &A_PtVel, const G4FieldTrack &B_PtVel, G4double safetyLast, G4int depth=-1)
G4double GetEpsilonStepFor()
void ReportImmediateHit(const char *MethodName, const G4ThreeVector &StartPosition, const G4ThreeVector &TrialPoint, G4double tolerance, unsigned long int numCalls)
G4bool GetAdjustementOfFoundIntersection()
void printStatus(const G4FieldTrack &startFT, const G4FieldTrack ¤tFT, G4double requestStep, G4double safety, G4int stepNum)
G4bool AdjustmentOfFoundIntersection(const G4ThreeVector &A, const G4ThreeVector &CurrentE_Point, const G4ThreeVector &CurrentF_Point, const G4ThreeVector &MomentumDir, const G4bool IntersectAF, G4ThreeVector &IntersectionPoint, G4double &NewSafety, G4double &fPrevSafety, G4ThreeVector &fPrevSftOrigin)
void ReportReversedPoints(std::ostringstream &ossMsg, const G4FieldTrack &StartPointVel, const G4FieldTrack &EndPointVel, G4double NewSafety, G4double epsStep, const G4FieldTrack &CurrentA_PointVelocity, const G4FieldTrack &CurrentB_PointVelocity, const G4FieldTrack &SubStart_PointVelocity, const G4ThreeVector &CurrentE_Point, const G4FieldTrack &ApproxIntersecPointV, G4int sbstp_no, G4int sbstp_no_p, G4int depth)
static constexpr double mm
static constexpr double perThousand
T max(const T t1, const T t2)
brief Return the largest of the two arguments