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 static 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;
239 CurrentB_PointVelocity,
251 "Intermediate F point is past end B point" );
260 G4ThreeVector ChordEF_Vector = CurrentF_Point - CurrentE_Point;
263 G4double MomDir_dot_Norm= NewMomentumDir.dot( NormalAtEntry ) ;
269 G4double ABchord_length = ChordAB.mag();
271 MomDir_dot_ABchord = (1.0 / ABchord_length)
272 * NewMomentumDir.dot( ChordAB );
274 ChordEF_Vector, NewMomentumDir, NormalAtEntry, validNormalAtE );
275 G4cout <<
" dot( MomentumDir, ABchord_unit ) = "
276 << MomDir_dot_ABchord <<
G4endl;
280 ( MomDir_dot_Norm >= 0.0 )
281 || (! validNormalAtE );
282 G4double EF_dist2 = ChordEF_Vector.mag2();
286 found_approximate_intersection =
true;
290 IntersectedOrRecalculatedFT = ApproxIntersecPointV;
291 IntersectedOrRecalculatedFT.
SetPosition( CurrentE_Point );
300 CurrentE_Point, CurrentF_Point, MomentumDir,
301 last_AF_intersection, IP, NewSafety,
302 previousSafety, previousSftOrigin );
303 if ( goodCorrection )
305 IntersectedOrRecalculatedFT = ApproxIntersecPointV;
328 NewSafety, previousSafety,
332 last_AF_intersection = Intersects_AF;
341 CurrentB_PointVelocity = ApproxIntersecPointV;
342 CurrentE_Point = PointG;
344 validIntersectP=
true;
349 validNormalAtE = validNormalLast;
354 fin_section_depth[depth]=
false;
359 G4cout <<
"G4PiF::LI> Investigating intermediate point"
361 <<
" on way to full s="
381 NewSafety, previousSafety,
399 CurrentA_PointVelocity = ApproxIntersecPointV;
400 CurrentE_Point = PointH;
402 validIntersectP =
true;
407 validNormalAtE = validNormalH;
411 if(fin_section_depth[depth])
421 if( ((Second_half)&&(depth==0)) || (first_section) )
423 there_is_no_intersection =
true;
429 substep_no_p = param_substeps+2;
434 sub_final_section =
true;
439 CurrentA_PointVelocity = CurrentB_PointVelocity;
440 CurrentB_PointVelocity = (depth==0) ? CurveEndPointVelocity
442 SubStart_PointVelocity = CurrentA_PointVelocity;
443 restoredFullEndpoint =
true;
445 validIntersectP=
false;
454 CurrentB_PointVelocity,
459 CurrentB_PointVelocity= RevisedB_FT;
468 if ( (fin_section_depth[depth])
469 &&( first_section || ((Second_half)&&(depth==0)) ) )
471 recalculatedEndPoint =
true;
472 IntersectedOrRecalculatedFT = RevisedB_FT;
484 std::ostringstream errmsg;
485 errmsg <<
"Location: " << MethodName
486 <<
"- After EndIf(Intersects_AF)" <<
G4endl;
488 CurveStartPointVelocity, CurveEndPointVelocity,
490 CurrentA_PointVelocity, CurrentB_PointVelocity,
491 SubStart_PointVelocity, CurrentE_Point,
492 ApproxIntersecPointV, substep_no, substep_no_p, depth);
495 if( restoredFullEndpoint )
497 fin_section_depth[depth] = restoredFullEndpoint;
498 restoredFullEndpoint =
false;
504 if( trigger_substepno_print == 0)
509 if( substep_no >= trigger_substepno_print )
511 G4cout <<
"Difficulty in converging in " << MethodName
513 <<
" Substep no = " << substep_no <<
G4endl;
514 if( substep_no == trigger_substepno_print )
516 printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
519 G4cout <<
" State of point A: ";
520 printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
521 -1.0, NewSafety, substep_no-1);
522 G4cout <<
" State of point B: ";
523 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
524 -1.0, NewSafety, substep_no);
530 }
while ( ( ! found_approximate_intersection )
531 && ( ! there_is_no_intersection )
532 && ( substep_no_p <= param_substeps) );
535 if( (!found_approximate_intersection) && (!there_is_no_intersection) )
548 if( (did_len < fraction_done*all_len)
549 && (depth<
max_depth) && (!sub_final_section) )
553 biggest_depth=
std::max(depth, biggest_depth);
558 first_section =
false;
560 G4double Sub_len = (all_len-did_len)/(2.);
574 G4bool goodAdvance = (lenAchieved >= adequateFraction * Sub_len);
580 G4cout <<
"MLL> AdvanceChordLimited not full at depth=" << depth
581 <<
" total length achieved = " << lenAchieved <<
" of "
582 << Sub_len <<
" fraction= ";
583 if (Sub_len != 0.0 ) {
G4cout << lenAchieved / Sub_len; }
584 else {
G4cout <<
"DivByZero"; }
585 G4cout <<
" Good-enough fraction = " << adequateFraction;
587 <<
" iteration " << substep_no
588 <<
" inner = " << substep_no_p <<
G4endl;
590 G4cout <<
" State at start is = " << CurrentA_PointVelocity
592 <<
" 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, CurrentA_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;
784 std::ostringstream message;
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 );
813 std::ostringstream message;
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] <<
" "
G4FieldTrack * ptrInterMedFT[max_depth+1]
void SetPosition(G4ThreeVector nPos)
void ReportFieldValue(const G4FieldTrack &locationPV, const char *nameLoc, const G4EquationOfMotion *equation)
G4double GetCurveLength() const
CLHEP::Hep3Vector G4ThreeVector
const G4MagIntegratorStepper * GetStepper() const
static const G4float tolerance
G4ThreeVector GetSurfaceNormal(const G4ThreeVector &CurrentInt_Point, G4bool &validNormal)
G4bool IntersectChord(const G4ThreeVector &StartPointA, const G4ThreeVector &EndPointB, G4double &NewSafety, G4double &PreviousSafety, G4ThreeVector &PreviousSftOrigin, G4double &LinearStepLength, G4ThreeVector &IntersectionPoint, G4bool *calledNavigator=0)
G4FieldTrack ApproxCurvePointV(const G4FieldTrack &curveAPointVelocity, const G4FieldTrack &curveBPointVelocity, const G4ThreeVector ¤tEPoint, G4double epsStep)
const G4ThreeVector & GetMomentumDir() const
G4bool EstimateIntersectionPoint(const G4FieldTrack &curveStartPointTangent, const G4FieldTrack &curveEndPointTangent, const G4ThreeVector &trialPoint, G4FieldTrack &intersectPointTangent, G4bool &recalculatedEndPoint, G4double &fPreviousSafety, G4ThreeVector &fPreviousSftOrigin)
static const double perThousand
void SetWarnSteps(unsigned int valWarn)
void printStatus(const G4FieldTrack &startFT, const G4FieldTrack ¤tFT, G4double requestStep, G4double safety, G4int stepNum)
unsigned long int fNumAdvanceTrials
G4double GetEpsilonStepFor()
void ReportImmediateHit(const char *MethodName, const G4ThreeVector &StartPosition, const G4ThreeVector &TrialPoint, double tolerance, unsigned long int numCalls)
G4EquationOfMotion * GetEquationOfMotion()
G4ThreeVector GetPosition() const
G4GLOB_DLL std::ostream G4cout
static const G4int max_depth
G4MultiLevelLocator(G4Navigator *theNavigator)
G4double GetRestMass() 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 ReportTrialStep(G4int step_no, const G4ThreeVector &ChordAB_v, const G4ThreeVector &ChordEF_v, const G4ThreeVector &NewMomentumDir, const G4ThreeVector &NormalAtEntry, G4bool validNormal)
G4Navigator * GetNavigatorFor()
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
void GetFieldValue(const G4double Point[4], G4double Field[]) const
G4double GetLabTimeOfFlight() const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4bool CheckAndReEstimateEndpoint(const G4FieldTrack &CurrentStartA, const G4FieldTrack &EstimatedEndB, G4FieldTrack &RevisedEndPoint, G4int &errorCode)
G4bool GetAdjustementOfFoundIntersection()
void SetMaxSteps(unsigned int valMax)
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)
G4MagInt_Driver * GetIntegrationDriver()
unsigned long int fNumAdvanceGood
G4double fiDeltaIntersection
unsigned long int fNumAdvanceFull
G4ChordFinder * GetChordFinderFor()
unsigned long int fNumCalls
G4bool AccurateAdvance(G4FieldTrack &y_current, G4double hstep, G4double eps, G4double hinitial=0.0)
virtual void LocateGlobalPointWithinVolume(const G4ThreeVector &position)
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)
G4ThreeVector GetMomentumDirection() const