45 fNumAdvanceFull(0.), fNumAdvanceGood(0), fNumAdvanceTrials(0)
124 G4bool& recalculatedEndPoint,
129 const char* MethodName=
"G4MultiLevelLocator::EstimateIntersectionPoint()";
131 G4bool found_approximate_intersection =
false;
132 G4bool there_is_no_intersection =
false;
134 G4FieldTrack CurrentA_PointVelocity = CurveStartPointVelocity;
135 G4FieldTrack CurrentB_PointVelocity = CurveEndPointVelocity;
137 G4bool validNormalAtE =
false;
140 G4FieldTrack ApproxIntersecPointV(CurveEndPointVelocity);
142 G4bool validIntersectP=
true;
144 G4bool last_AF_intersection =
false;
148 G4bool first_section =
true;
150 recalculatedEndPoint =
false;
152 G4bool restoredFullEndpoint =
false;
154 unsigned int substep_no = 0;
170 const G4int param_substeps=5;
174 G4bool Second_half =
false;
182 unsigned int depth=0;
186 unsigned int trigger_substepno_print=0;
188 unsigned int biggest_depth= 0;
190 #if (G4DEBUG_FIELD>1)
192 if( (TrialPoint - StartPosition).mag2() < tolerance*tolerance)
217 fin_section_depth[idepth]=
true;
221 G4FieldTrack SubStart_PointVelocity = CurveStartPointVelocity;
225 unsigned int substep_no_p = 0;
226 G4bool sub_final_section =
false;
228 SubStart_PointVelocity = CurrentA_PointVelocity;
240 CurrentB_PointVelocity,
252 "Intermediate F point is past end B point" );
261 G4ThreeVector ChordEF_Vector = CurrentF_Point - CurrentE_Point;
264 G4double MomDir_dot_Norm= NewMomentumDir.
dot( NormalAtEntry ) ;
272 MomDir_dot_ABchord = (1.0 / ABchord_length)
273 * NewMomentumDir.
dot( ChordAB );
275 ChordEF_Vector, NewMomentumDir, NormalAtEntry, validNormalAtE );
276 G4cout <<
" dot( MomentumDir, ABchord_unit ) = "
277 << MomDir_dot_ABchord <<
G4endl;
281 ( MomDir_dot_Norm >= 0.0 )
282 || (! validNormalAtE );
287 found_approximate_intersection =
true;
291 IntersectedOrRecalculatedFT = ApproxIntersecPointV;
292 IntersectedOrRecalculatedFT.
SetPosition( CurrentE_Point );
301 CurrentE_Point, CurrentF_Point, MomentumDir,
302 last_AF_intersection, IP, NewSafety,
303 previousSafety, previousSftOrigin );
304 if ( goodCorrection )
306 IntersectedOrRecalculatedFT = ApproxIntersecPointV;
329 NewSafety, previousSafety,
333 last_AF_intersection = Intersects_AF;
342 CurrentB_PointVelocity = ApproxIntersecPointV;
343 CurrentE_Point = PointG;
345 validIntersectP=
true;
350 validNormalAtE = validNormalLast;
355 fin_section_depth[depth]=
false;
360 G4cout <<
"G4PiF::LI> Investigating intermediate point"
362 <<
" on way to full s="
382 NewSafety, previousSafety,
400 CurrentA_PointVelocity = ApproxIntersecPointV;
401 CurrentE_Point = PointH;
403 validIntersectP =
true;
408 validNormalAtE = validNormalH;
412 if(fin_section_depth[depth])
422 if( ((Second_half)&&(depth==0)) || (first_section) )
424 there_is_no_intersection =
true;
430 substep_no_p = param_substeps+2;
435 sub_final_section =
true;
440 CurrentA_PointVelocity = CurrentB_PointVelocity;
441 CurrentB_PointVelocity = (depth==0) ? CurveEndPointVelocity
443 SubStart_PointVelocity = CurrentA_PointVelocity;
444 restoredFullEndpoint =
true;
446 validIntersectP=
false;
455 CurrentB_PointVelocity,
460 CurrentB_PointVelocity= RevisedB_FT;
469 if ( (fin_section_depth[depth])
470 &&( first_section || ((Second_half)&&(depth==0)) ) )
472 recalculatedEndPoint =
true;
473 IntersectedOrRecalculatedFT = RevisedB_FT;
485 std::ostringstream errmsg;
486 errmsg <<
"Location: " << MethodName
487 <<
"- After EndIf(Intersects_AF)" <<
G4endl;
489 CurveStartPointVelocity, CurveEndPointVelocity,
491 CurrentA_PointVelocity, CurrentB_PointVelocity,
492 SubStart_PointVelocity, CurrentE_Point,
493 ApproxIntersecPointV, substep_no, substep_no_p, depth);
496 if( restoredFullEndpoint )
498 fin_section_depth[depth] = restoredFullEndpoint;
499 restoredFullEndpoint =
false;
505 if( trigger_substepno_print == 0)
510 if( substep_no >= trigger_substepno_print )
512 G4cout <<
"Difficulty in converging in " << MethodName
514 <<
" Substep no = " << substep_no <<
G4endl;
515 if( substep_no == trigger_substepno_print )
517 printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
520 G4cout <<
" State of point A: ";
521 printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
522 -1.0, NewSafety, substep_no-1);
523 G4cout <<
" State of point B: ";
524 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
525 -1.0, NewSafety, substep_no);
531 }
while ( ( ! found_approximate_intersection )
532 && ( ! there_is_no_intersection )
533 && ( substep_no_p <= param_substeps) );
536 if( (!found_approximate_intersection) && (!there_is_no_intersection) )
549 if( (did_len < fraction_done*all_len)
550 && (depth<
max_depth) && (!sub_final_section) )
554 biggest_depth=
std::max(depth, biggest_depth);
559 first_section =
false;
561 G4double Sub_len = (all_len-did_len)/(2.);
575 G4bool goodAdvance = (lenAchieved >= adequateFraction * Sub_len);
581 G4cout <<
"MLL> AdvanceChordLimited not full at depth=" << depth
582 <<
" total length achieved = " << lenAchieved <<
" of "
583 << Sub_len <<
" fraction= ";
584 if (Sub_len != 0.0 ) {
G4cout << lenAchieved / Sub_len; }
585 else {
G4cout <<
"DivByZero"; }
586 G4cout <<
" Good-enough fraction = " << adequateFraction;
588 <<
" iteration " << substep_no
589 <<
" inner = " << substep_no_p <<
G4endl;
591 G4cout <<
" State at start is = " << CurrentA_PointVelocity
593 <<
" at end (midpoint)= " << midPoint <<
G4endl;
599 G4cout <<
" Original Start = "
600 << CurveStartPointVelocity <<
G4endl;
601 G4cout <<
" Original End = "
602 << CurveEndPointVelocity <<
G4endl;
603 G4cout <<
" Original TrialPoint = "
605 G4cout <<
" (this is STRICT mode) "
606 <<
" num Splits= " << numSplits;
612 CurrentB_PointVelocity = midPoint;
616 SubStart_PointVelocity = CurrentA_PointVelocity;
625 NewSafety, previousSafety,
626 previousSftOrigin, distAB,
630 last_AF_intersection = Intersects_AB;
631 CurrentE_Point = PointGe;
632 fin_section_depth[depth]=
true;
634 validIntersectP=
true;
639 validNormalAtE = validNormalAB;
648 validIntersectP=
false;
653 unsigned int levelPops=0;
655 G4bool unfinished = Second_half;
656 while ( unfinished && (depth>0) )
676 CurrentB_PointVelocity,
681 CurrentB_PointVelocity= RevisedEndPointFT;
685 recalculatedEndPoint =
true;
686 IntersectedOrRecalculatedFT = RevisedEndPointFT;
692 std::ostringstream errmsg;
693 errmsg <<
"Location: " << MethodName <<
"- Second-Half" <<
G4endl;
695 CurveStartPointVelocity, CurveEndPointVelocity,
697 CurrentA_PointVelocity, CurrentB_PointVelocity,
698 SubStart_PointVelocity, CurrentE_Point,
699 ApproxIntersecPointV, substep_no, substep_no_p, depth);
707 previousSftOrigin, distAB,
711 last_AF_intersection = Intersects_AB;
712 CurrentE_Point = PointGe;
714 validIntersectP=
true;
719 validNormalAtE = validNormalAB;
723 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;
826 if( found_approximate_intersection )
828 assert( validApproxIntPV &&
829 "Approximate Intersection must not have been invalidated." );
833 return (!there_is_no_intersection) && found_approximate_intersection;
840 G4cout <<
" Number of split level ('advances'): "
842 G4cout <<
" Number of full advances: "
844 G4cout <<
" Number of good advances: "
852 enum { maxNumFieldComp= 24 };
855 G4double startPoint[4] = { position.
x(), position.
y(), position.
z(),
858 for (
unsigned int i=0; i<maxNumFieldComp; ++i )
863 G4cout <<
" B-field value (" << nameLoc <<
")= "
864 << FieldVec[0] <<
" " << FieldVec[1] <<
" " << FieldVec[2];
867 FieldVec[5] ).
mag2();
870 G4cout <<
" Electric = " << FieldVec[3] <<
" "
871 << FieldVec[4] <<
" "
G4VIntegrationDriver * GetIntegrationDriver()
void ReportTrialStep(G4int step_no, const G4ThreeVector &ChordAB_v, const G4ThreeVector &ChordEF_v, const G4ThreeVector &NewMomentumDir, const G4ThreeVector &NormalAtEntry, G4bool validNormal)
T max(const T t1, const T t2)
brief Return the largest of the two arguments
static constexpr double perThousand
unsigned long int fNumAdvanceGood
CLHEP::Hep3Vector G4ThreeVector
G4ChordFinder * GetChordFinderFor()
G4bool CheckAndReEstimateEndpoint(const G4FieldTrack &CurrentStartA, const G4FieldTrack &EstimatedEndB, G4FieldTrack &RevisedEndPoint, G4int &errorCode)
void printStatus(const G4FieldTrack &startFT, const G4FieldTrack ¤tFT, G4double requestStep, G4double safety, G4int stepNum)
static const G4int max_depth
virtual void LocateGlobalPointWithinVolume(const G4ThreeVector &position)
G4double GetRestMass() const
void message(RunManager *runmanager)
double dot(const Hep3Vector &) const
G4bool GetAdjustementOfFoundIntersection()
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)
G4FieldTrack ApproxCurvePointV(const G4FieldTrack &curveAPointVelocity, const G4FieldTrack &curveBPointVelocity, const G4ThreeVector ¤tEPoint, G4double epsStep)
void SetPosition(G4ThreeVector nPos)
G4bool IntersectChord(const G4ThreeVector &StartPointA, const G4ThreeVector &EndPointB, G4double &NewSafety, G4double &PreviousSafety, G4ThreeVector &PreviousSftOrigin, G4double &LinearStepLength, G4ThreeVector &IntersectionPoint, G4bool *calledNavigator=0)
G4ThreeVector GetPosition() const
unsigned long int fNumCalls
G4ThreeVector GetMomentumDirection() const
void ReportFieldValue(const G4FieldTrack &locationPV, const char *nameLoc, const G4EquationOfMotion *equation)
unsigned long int fNumAdvanceFull
G4Navigator * GetNavigatorFor()
G4bool EstimateIntersectionPoint(const G4FieldTrack &curveStartPointTangent, const G4FieldTrack &curveEndPointTangent, const G4ThreeVector &trialPoint, G4FieldTrack &intersectPointTangent, G4bool &recalculatedEndPoint, G4double &fPreviousSafety, G4ThreeVector &fPreviousSftOrigin)
void GetFieldValue(const G4double Point[4], G4double Field[]) const
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)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
G4MultiLevelLocator(G4Navigator *theNavigator)
void SetWarnSteps(unsigned int valWarn)
void ReportImmediateHit(const char *MethodName, const G4ThreeVector &StartPosition, const G4ThreeVector &TrialPoint, double tolerance, unsigned long int numCalls)
const G4ThreeVector & GetMomentumDir() const
static constexpr double mm
G4GLOB_DLL std::ostream G4cout
G4FieldTrack * ptrInterMedFT[max_depth+1]
unsigned long int fNumAdvanceTrials
void SetMaxSteps(unsigned int valMax)
G4double GetEpsilonStepFor()
G4double GetLabTimeOfFlight() const
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)
virtual G4bool AccurateAdvance(G4FieldTrack &track, G4double hstep, G4double eps, G4double hinitial=0)=0
G4double GetCurveLength() const
G4ThreeVector GetSurfaceNormal(const G4ThreeVector &CurrentInt_Point, G4bool &validNormal)
G4double fiDeltaIntersection