45 for (
auto idepth=0; idepth<
max_depth+1; ++idepth )
60 for (
auto idepth=0; idepth<
max_depth+1; ++idepth )
117 G4bool& recalculatedEndPoint,
122 const char* MethodName=
"G4MultiLevelLocator::EstimateIntersectionPoint()";
124 G4bool found_approximate_intersection =
false;
125 G4bool there_is_no_intersection =
false;
127 G4FieldTrack CurrentA_PointVelocity = CurveStartPointVelocity;
128 G4FieldTrack CurrentB_PointVelocity = CurveEndPointVelocity;
130 G4bool validNormalAtE =
false;
133 G4FieldTrack ApproxIntersecPointV(CurveEndPointVelocity);
134 G4bool validIntersectP=
true;
136 G4bool last_AF_intersection =
false;
143 G4bool first_section =
true;
144 recalculatedEndPoint =
false;
145 G4bool restoredFullEndpoint =
false;
147 unsigned int substep_no = 0;
163 const G4int param_substeps = 5;
167 G4bool Second_half =
false;
175 unsigned int depth = 0;
179 unsigned int trigger_substepno_print = 0;
181 unsigned int biggest_depth = 0;
183 #if (G4DEBUG_FIELD>1)
185 if( (TrialPoint - StartPosition).mag2() <
sqr(tolerance))
200 for (
auto idepth=1; idepth<
max_depth+1; ++idepth )
208 for (
auto idepth=0; idepth<
max_depth; ++idepth )
210 fin_section_depth[idepth] =
true;
214 G4FieldTrack SubStart_PointVelocity = CurveStartPointVelocity;
218 unsigned int substep_no_p = 0;
219 G4bool sub_final_section =
false;
221 SubStart_PointVelocity = CurrentA_PointVelocity;
233 CurrentB_PointVelocity,
243 "Intermediate F point is past end B point" );
252 G4ThreeVector ChordEF_Vector = CurrentF_Point - CurrentE_Point;
255 G4double MomDir_dot_Norm = NewMomentumDir.
dot( NormalAtEntry );
263 MomDir_dot_ABchord = (1.0 / ABchord_length)
264 * NewMomentumDir.
dot( ChordAB );
266 ChordEF_Vector, NewMomentumDir, NormalAtEntry, validNormalAtE );
267 G4cout <<
" dot( MomentumDir, ABchord_unit ) = "
268 << MomDir_dot_ABchord <<
G4endl;
272 ( MomDir_dot_Norm >= 0.0 )
273 || (! validNormalAtE );
278 found_approximate_intersection =
true;
282 IntersectedOrRecalculatedFT = ApproxIntersecPointV;
283 IntersectedOrRecalculatedFT.
SetPosition( CurrentE_Point );
292 CurrentE_Point, CurrentF_Point, MomentumDir,
293 last_AF_intersection, IP, NewSafety,
294 previousSafety, previousSftOrigin );
295 if ( goodCorrection )
297 IntersectedOrRecalculatedFT = ApproxIntersecPointV;
320 NewSafety, previousSafety,
324 last_AF_intersection = Intersects_AF;
333 CurrentB_PointVelocity = ApproxIntersecPointV;
334 CurrentE_Point = PointG;
336 validIntersectP =
true;
340 validNormalAtE = validNormalLast;
345 fin_section_depth[depth] =
false;
350 G4cout <<
"G4PiF::LI> Investigating intermediate point"
352 <<
" on way to full s="
372 NewSafety, previousSafety,
390 CurrentA_PointVelocity = ApproxIntersecPointV;
391 CurrentE_Point = PointH;
393 validIntersectP =
true;
397 validNormalAtE = validNormalH;
401 if( fin_section_depth[depth] )
411 if( ((Second_half)&&(depth==0)) || (first_section) )
413 there_is_no_intersection =
true;
419 substep_no_p = param_substeps+2;
424 sub_final_section =
true;
429 CurrentA_PointVelocity = CurrentB_PointVelocity;
430 CurrentB_PointVelocity = (depth==0) ? CurveEndPointVelocity
432 SubStart_PointVelocity = CurrentA_PointVelocity;
433 restoredFullEndpoint =
true;
435 validIntersectP =
false;
440 G4int errorEndPt = 0;
442 G4bool recalculatedB=
false;
443 if( driverReIntegrates )
447 CurrentB_PointVelocity,
452 CurrentB_PointVelocity = RevisedB_FT;
459 if ( (fin_section_depth[depth])
460 &&( first_section || ((Second_half)&&(depth==0)) ) )
462 recalculatedEndPoint =
true;
463 IntersectedOrRecalculatedFT = RevisedB_FT;
481 std::ostringstream errmsg;
482 errmsg <<
"Location: " << MethodName
483 <<
"- After EndIf(Intersects_AF)" <<
G4endl;
485 CurveStartPointVelocity, CurveEndPointVelocity,
487 CurrentA_PointVelocity, CurrentB_PointVelocity,
488 SubStart_PointVelocity, CurrentE_Point,
489 ApproxIntersecPointV, substep_no, substep_no_p, depth);
492 if( restoredFullEndpoint )
494 fin_section_depth[depth] = restoredFullEndpoint;
495 restoredFullEndpoint =
false;
501 if( trigger_substepno_print == 0)
506 if( substep_no >= trigger_substepno_print )
508 G4cout <<
"Difficulty in converging in " << MethodName
510 <<
" Substep no = " << substep_no <<
G4endl;
511 if( substep_no == trigger_substepno_print )
513 printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
514 -1.0, NewSafety, 0 );
516 G4cout <<
" State of point A: ";
517 printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
518 -1.0, NewSafety, substep_no-1 );
519 G4cout <<
" State of point B: ";
520 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
521 -1.0, NewSafety, substep_no );
527 }
while ( ( ! found_approximate_intersection )
528 && ( ! there_is_no_intersection )
529 && ( substep_no_p <= param_substeps) );
532 if( (!found_approximate_intersection) && (!there_is_no_intersection) )
545 if( (did_len < fraction_done*all_len)
546 && (depth<
max_depth) && (!sub_final_section) )
550 biggest_depth =
std::max(depth, biggest_depth);
555 first_section =
false;
557 G4double Sub_len = (all_len-did_len)/(2.);
560 integrDriver->AccurateAdvance(midPoint, Sub_len,
fiEpsilonStep);
569 G4bool goodAdvance = (lenAchieved >= adequateFraction * Sub_len);
575 G4cout <<
"MLL> AdvanceChordLimited not full at depth=" << depth
576 <<
" total length achieved = " << lenAchieved <<
" of "
577 << Sub_len <<
" fraction= ";
578 if (Sub_len != 0.0 ) {
G4cout << lenAchieved / Sub_len; }
579 else {
G4cout <<
"DivByZero"; }
580 G4cout <<
" Good-enough fraction = " << adequateFraction;
582 <<
" iteration " << substep_no
583 <<
" inner = " << substep_no_p <<
G4endl;
585 G4cout <<
" State at start is = " << CurrentA_PointVelocity
587 <<
" at end (midpoint)= " << midPoint <<
G4endl;
593 G4cout <<
" Original Start = "
594 << CurveStartPointVelocity <<
G4endl;
595 G4cout <<
" Original End = "
596 << CurveEndPointVelocity <<
G4endl;
597 G4cout <<
" Original TrialPoint = "
599 G4cout <<
" (this is STRICT mode) "
600 <<
" num Splits= " << numSplits;
606 CurrentB_PointVelocity = midPoint;
610 SubStart_PointVelocity = CurrentA_PointVelocity;
619 NewSafety, previousSafety,
620 previousSftOrigin, distAB,
624 last_AF_intersection = Intersects_AB;
625 CurrentE_Point = PointGe;
626 fin_section_depth[depth] =
true;
628 validIntersectP =
true;
632 validNormalAtE = validNormalAB;
641 validIntersectP=
false;
645 unsigned int levelPops = 0;
647 G4bool unfinished = Second_half;
648 while ( unfinished && (depth>0) )
661 G4int errorEndPt = 0;
663 G4bool recalculatedB=
false;
664 if( driverReIntegrates )
669 G4FieldTrack RevisedEndPointFT = CurrentB_PointVelocity;
672 CurrentB_PointVelocity,
677 CurrentB_PointVelocity = RevisedEndPointFT;
681 recalculatedEndPoint =
true;
682 IntersectedOrRecalculatedFT = RevisedEndPointFT;
694 std::ostringstream errmsg;
695 errmsg <<
"Location: " << MethodName <<
"- Second-Half" <<
G4endl;
697 CurveStartPointVelocity, CurveEndPointVelocity,
699 CurrentA_PointVelocity, CurrentB_PointVelocity,
700 SubStart_PointVelocity, CurrentE_Point,
701 ApproxIntersecPointV, substep_no, substep_no_p, depth);
709 previousSftOrigin, distAB,
713 last_AF_intersection = Intersects_AB;
714 CurrentE_Point = PointGe;
716 validIntersectP =
true;
720 validNormalAtE = validNormalAB;
724 validIntersectP =
false;
727 there_is_no_intersection =
true;
731 fin_section_depth[depth] =
true;
732 unfinished = !validIntersectP;
735 if( ! ( validIntersectP || there_is_no_intersection ) )
738 G4cout <<
"MLL - WARNING Potential FAILURE: Conditions not met!"
740 <<
" Depth = " << depth <<
G4endl
741 <<
" Levels popped = " << levelPops
742 <<
" Num Substeps= " << substep_no <<
G4endl;
743 G4cout <<
" Found intersection= " << found_approximate_intersection
747 CurveStartPointVelocity, CurveEndPointVelocity,
748 substep_no, CurrentA_PointVelocity,
749 CurrentB_PointVelocity,
756 assert( validIntersectP || there_is_no_intersection
757 || found_approximate_intersection);
759 }
while ( ( ! found_approximate_intersection )
760 && ( ! there_is_no_intersection )
763 if( substep_no > max_no_seen )
765 max_no_seen = substep_no;
769 trigger_substepno_print = max_no_seen-20;
774 if( !there_is_no_intersection && !found_approximate_intersection )
780 recalculatedEndPoint =
true;
781 IntersectedOrRecalculatedFT = CurrentA_PointVelocity;
782 found_approximate_intersection =
false;
786 message <<
"Convergence is requiring too many substeps: "
787 << substep_no <<
" ( limit = "<<
fMaxSteps <<
")"
789 <<
" Abandoning effort to intersect. " << G4endl <<
G4endl;
790 message <<
" Number of calls to MLL: " <<
fNumCalls;
791 message <<
" iteration = " << substep_no <<G4endl <<
G4endl;
793 message.precision( 12 );
796 message <<
" Undertaken only length: " << done_len
797 <<
" out of " << full_len <<
" required." << G4endl
798 <<
" Remaining length = " << full_len - done_len;
800 message <<
" Start and end-point of requested Step:" <<
G4endl;
801 printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
802 -1.0, NewSafety, 0, message, -1 );
803 message <<
" Start and end-point of current Sub-Step:" <<
G4endl;
804 printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
805 -1.0, NewSafety, substep_no-1, message, -1 );
806 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
807 -1.0, NewSafety, substep_no, message, -1 );
814 message <<
"Many substeps while trying to locate intersection."
816 <<
" Undertaken length: "
818 <<
" - Needed: " << substep_no <<
" substeps." <<
G4endl
820 <<
" and maximum substeps = " <<
fMaxSteps;
825 return (!there_is_no_intersection) && found_approximate_intersection;
832 G4cout <<
" Number of split level ('advances'): "
834 G4cout <<
" Number of full advances: "
836 G4cout <<
" Number of good advances: "
844 enum { maxNumFieldComp = 24 };
847 G4double startPoint[4] = { position.
x(), position.
y(), position.
z(),
850 for (
auto i=0; i<maxNumFieldComp; ++i )
855 G4cout <<
" B-field value (" << nameLoc <<
")= "
856 << FieldVec[0] <<
" " << FieldVec[1] <<
" " << FieldVec[2];
859 FieldVec[5] ).
mag2();
862 G4cout <<
" Electric = " << FieldVec[3] <<
" "
863 << FieldVec[4] <<
" "