72 G4cout <<
"G4BrentLocator::Location with Max Number of Steps="
74 G4cout <<
"G4BrentLocator::ReIntegrateEndPoint was called "
76 <<
" times and for depth algorithm "
120 G4bool& recalculatedEndPoint,
127 G4bool found_approximate_intersection =
false;
128 G4bool there_is_no_intersection =
false;
130 G4FieldTrack CurrentA_PointVelocity = CurveStartPointVelocity;
131 G4FieldTrack CurrentB_PointVelocity = CurveEndPointVelocity;
133 G4bool validNormalAtE =
false;
136 G4FieldTrack ApproxIntersecPointV(CurveEndPointVelocity);
138 G4bool last_AF_intersection =
false;
142 G4bool first_section =
true;
143 recalculatedEndPoint =
false;
145 G4bool restoredFullEndpoint =
false;
148 G4int substep_no = 0;
152 const G4int max_substeps= 10000;
153 const G4int warn_substeps= 1000;
173 const G4int param_substeps=50;
177 G4bool Second_half =
false;
192 if( (TrialPoint - StartPosition).mag() < tolerance *
CLHEP::mm )
194 G4Exception(
"G4BrentLocator::EstimateIntersectionPoint()",
196 "Intersection point F is exactly at start point A." );
214 fin_section_depth[idepth]=
true;
219 G4FieldTrack SubStart_PointVelocity = CurveStartPointVelocity;
223 G4int substep_no_p = 0;
224 G4bool sub_final_section =
false;
226 SubStart_PointVelocity = CurrentA_PointVelocity;
240 CurrentB_PointVelocity,
249 G4Exception(
"G4BrentLocator::EstimateIntersectionPoint()",
251 "Intermediate F point is past end B point" );
260 G4ThreeVector ChordEF_Vector = CurrentF_Point - CurrentE_Point;
262 G4double MomDir_dot_Norm= NewMomentumDir.
dot( NormalAtEntry ) ;
268 ChordEF_Vector, NewMomentumDir, NormalAtEntry, validNormalAtE );
272 adequate_angle = ( 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 fPreviousSafety, fPreviousSftOrigin );
295 if ( goodCorrection )
297 IntersectedOrRecalculatedFT = ApproxIntersecPointV;
320 G4bool usedNavigatorAF =
false;
322 NewSafety,fPreviousSafety,
327 last_AF_intersection = Intersects_AF;
338 CurrentA_PointVelocity, CurrentB_PointVelocity,
339 EndPoint,CurrentE_Point, CurrentF_Point,PointG,
341 CurrentB_PointVelocity = EndPoint;
342 CurrentE_Point = PointG;
350 validNormalAtE = validNormalLast;
355 fin_section_depth[depth] =
false;
359 G4cout <<
"G4PiF::LI> Investigating intermediate point"
361 <<
" on way to full s="
376 G4bool usedNavigatorFB =
false;
382 NewSafety,fPreviousSafety,
403 CurrentA_PointVelocity,CurrentB_PointVelocity,
404 InterMed,CurrentE_Point,CurrentF_Point,PointH,
406 CurrentA_PointVelocity = InterMed;
407 CurrentE_Point = PointH;
413 validNormalAtE= validNormalLast;
419 if( fin_section_depth[depth] )
429 if( ((Second_half)&&(depth==0)) || (first_section) )
431 there_is_no_intersection =
true;
438 substep_no_p = param_substeps+2;
444 sub_final_section =
true;
453 CurrentA_PointVelocity = CurrentB_PointVelocity;
454 CurrentB_PointVelocity = CurveEndPointVelocity;
455 SubStart_PointVelocity = CurrentA_PointVelocity;
458 CurrentB_PointVelocity,
462 restoredFullEndpoint =
true;
469 CurrentA_PointVelocity = CurrentB_PointVelocity;
471 SubStart_PointVelocity = CurrentA_PointVelocity;
474 CurrentB_PointVelocity,
477 restoredFullEndpoint =
true;
501 CurrentB_PointVelocity,
505 CurrentB_PointVelocity = newEndPointFT;
507 if ( (fin_section_depth[depth])
508 &&( first_section || ((Second_half)&&(depth==0)) ) )
510 recalculatedEndPoint =
true;
511 IntersectedOrRecalculatedFT = newEndPointFT;
515 if( curveDist < 0.0 )
518 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
519 -1.0, NewSafety, substep_no );
521 message <<
"Error in advancing propagation." <<
G4endl
522 <<
" Error in advancing propagation." <<
G4endl
523 <<
" Point A (start) is " << CurrentA_PointVelocity
525 <<
" Point B (end) is " << CurrentB_PointVelocity
527 <<
" Curve distance is " << curveDist <<
G4endl
529 <<
"The final curve point is not further along"
530 <<
" than the original!" <<
G4endl;
532 if( recalculatedEndPoint )
534 message <<
"Recalculation of EndPoint was called with fEpsStep= "
537 oldprc =
G4cerr.precision(20);
538 message <<
" Point A (Curve start) is " << CurveStartPointVelocity
540 <<
" Point B (Curve end) is " << CurveEndPointVelocity
542 <<
" Point A (Current start) is " << CurrentA_PointVelocity
544 <<
" Point B (Current end) is " << CurrentB_PointVelocity
546 <<
" Point S (Sub start) is " << SubStart_PointVelocity
548 <<
" Point E (Trial Point) is " << CurrentE_Point
550 <<
" Old Point F(Intersection) is " << CurrentF_Point
552 <<
" New Point F(Intersection) is " << ApproxIntersecPointV
554 <<
" LocateIntersection parameters are : Substep no= "
555 << substep_no << G4endl
556 <<
" Substep depth no= "<< substep_no_p <<
" Depth= "
558 <<
" Restarted no= "<< restartB <<
" Epsilon= "
561 G4cerr.precision( oldprc );
563 G4Exception(
"G4BrentLocator::EstimateIntersectionPoint()",
567 if(restoredFullEndpoint)
569 fin_section_depth[depth] = restoredFullEndpoint;
570 restoredFullEndpoint =
false;
575 #ifdef G4DEBUG_LOCATE_INTERSECTION
576 G4int trigger_substepno_print= warn_substeps - 20 ;
578 if( substep_no >= trigger_substepno_print )
580 G4cout <<
"Difficulty in converging in "
581 <<
"G4BrentLocator::EstimateIntersectionPoint()"
583 <<
" Substep no = " << substep_no <<
G4endl;
584 if( substep_no == trigger_substepno_print )
586 printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
589 G4cout <<
" State of point A: ";
590 printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
591 -1.0, NewSafety, substep_no-1, 0);
592 G4cout <<
" State of point B: ";
593 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
594 -1.0, NewSafety, substep_no);
600 }
while ( ( ! found_approximate_intersection )
601 && ( ! there_is_no_intersection )
602 && ( substep_no_p <= param_substeps) );
604 first_section =
false;
606 if( (!found_approximate_intersection) && (!there_is_no_intersection) )
619 if ( ( did_len < fraction_done*all_len )
620 && (depth<
max_depth) && (!sub_final_section) )
625 G4double Sub_len = (all_len-did_len)/(2.);
635 SubStart_PointVelocity = CurrentA_PointVelocity;
644 NewSafety, fPreviousSafety,
645 fPreviousSftOrigin,stepLengthAB,
649 last_AF_intersection = Intersects_AB;
650 CurrentE_Point = PointGe;
651 fin_section_depth[depth]=
true;
657 validNormalAtE= validNormalAB;
668 if( (Second_half)&&(depth!=0) )
694 CurrentB_PointVelocity,
698 CurrentB_PointVelocity = newEndPointFT;
701 recalculatedEndPoint =
true;
702 IntersectedOrRecalculatedFT = newEndPointFT;
713 fPreviousSftOrigin,stepLengthAB, PointGe);
716 last_AF_intersection = Intersects_AB;
717 CurrentE_Point = PointGe;
721 validNormalAtE = validNormalAB;
725 fin_section_depth[depth]=
true;
729 }
while ( ( ! found_approximate_intersection )
730 && ( ! there_is_no_intersection )
731 && ( substep_no <= max_substeps) );
733 if( substep_no > max_no_seen )
735 max_no_seen = substep_no;
736 #ifdef G4DEBUG_LOCATE_INTERSECTION
737 if( max_no_seen > warn_substeps )
739 trigger_substepno_print = max_no_seen-20;
744 if( ( substep_no >= max_substeps)
745 && !there_is_no_intersection
746 && !found_approximate_intersection )
748 G4cout <<
"ERROR - G4BrentLocator::EstimateIntersectionPoint()" <<
G4endl
749 <<
" Start and end-point of requested Step:" <<
G4endl;
750 printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
752 G4cout <<
" Start and end-point of current Sub-Step:" <<
G4endl;
753 printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
754 -1.0, NewSafety, substep_no-1);
755 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
756 -1.0, NewSafety, substep_no);
758 message <<
"Too many substeps!" << G4endl
759 <<
" Convergence is requiring too many substeps: "
760 << substep_no << G4endl
761 <<
" Abandoning effort to intersect. " << G4endl
762 <<
" Found intersection = "
763 << found_approximate_intersection << G4endl
764 <<
" Intersection exists = "
765 << !there_is_no_intersection <<
G4endl;
766 oldprc =
G4cout.precision( 10 );
769 message <<
" Undertaken only length: " << done_len
770 <<
" out of " << full_len <<
" required." << G4endl
771 <<
" Remaining length = " << full_len - done_len;
772 G4cout.precision( oldprc );
774 G4Exception(
"G4BrentLocator::EstimateIntersectionPoint()",
777 else if( substep_no >= warn_substeps )
779 oldprc=
G4cout.precision( 10 );
781 message <<
"Many substeps while trying to locate intersection."
783 <<
" Undertaken length: "
785 <<
" - Needed: " << substep_no <<
" substeps." <<
G4endl
786 <<
" Warning level = " << warn_substeps
787 <<
" and maximum substeps = " << max_substeps;
788 G4Exception(
"G4BrentLocator::EstimateIntersectionPoint()",
790 G4cout.precision( oldprc );
792 return !there_is_no_intersection;
G4VIntegrationDriver * GetIntegrationDriver()
void ReportTrialStep(G4int step_no, const G4ThreeVector &ChordAB_v, const G4ThreeVector &ChordEF_v, const G4ThreeVector &NewMomentumDir, const G4ThreeVector &NormalAtEntry, G4bool validNormal)
G4int maxNumberOfCallsToReIntegration_depth
G4ChordFinder * GetChordFinderFor()
void printStatus(const G4FieldTrack &startFT, const G4FieldTrack ¤tFT, G4double requestStep, G4double safety, G4int stepNum)
virtual void LocateGlobalPointWithinVolume(const G4ThreeVector &position)
void message(RunManager *runmanager)
G4int maxNumberOfStepsForIntersection
double dot(const Hep3Vector &) const
G4bool GetAdjustementOfFoundIntersection()
G4FieldTrack ApproxCurvePointS(const G4FieldTrack &curveAPointVelocity, const G4FieldTrack &curveBPointVelocity, const G4FieldTrack &ApproxCurveV, const G4ThreeVector ¤tEPoint, const G4ThreeVector ¤tFPoint, const G4ThreeVector &PointG, G4bool first, G4double epsStep)
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)
G4int maxNumberOfCallsToReIntegration
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)
G4FieldTrack * ptrInterMedFT[max_depth+1]
G4ThreeVector GetPosition() const
G4ThreeVector GetMomentumDirection() const
#define fPreviousSftOrigin
static const G4int max_depth
G4Navigator * GetNavigatorFor()
G4FieldTrack ReEstimateEndpoint(const G4FieldTrack &CurrentStateA, const G4FieldTrack &EstimtdEndStateB, G4double linearDistSq, G4double curveDist)
G4double GetDeltaIntersectionFor()
G4GLOB_DLL std::ostream G4cerr
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
G4BrentLocator(G4Navigator *theNavigator)
G4bool EstimateIntersectionPoint(const G4FieldTrack &curveStartPointTangent, const G4FieldTrack &curveEndPointTangent, const G4ThreeVector &trialPoint, G4FieldTrack &intersectPointTangent, G4bool &recalculatedEndPoint, G4double &fPreviousSafety, G4ThreeVector &fPreviousSftOrigin)
const G4ThreeVector & GetMomentumDir() const
static constexpr double mm
G4GLOB_DLL std::ostream G4cout
G4double GetEpsilonStepFor()
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