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;
239 CurrentB_PointVelocity,
248 G4Exception(
"G4BrentLocator::EstimateIntersectionPoint()",
250 "Intermediate F point is past end B point" );
259 G4ThreeVector ChordEF_Vector = CurrentF_Point - CurrentE_Point;
261 G4double MomDir_dot_Norm= NewMomentumDir.dot( NormalAtEntry ) ;
267 ChordEF_Vector, NewMomentumDir, NormalAtEntry, validNormalAtE );
271 adequate_angle = ( MomDir_dot_Norm >= 0.0 )
272 || (! validNormalAtE );
273 G4double EF_dist2 = ChordEF_Vector.mag2();
277 found_approximate_intersection =
true;
281 IntersectedOrRecalculatedFT = ApproxIntersecPointV;
282 IntersectedOrRecalculatedFT.
SetPosition( CurrentE_Point );
291 CurrentE_Point, CurrentF_Point, MomentumDir,
292 last_AF_intersection, IP, NewSafety,
293 fPreviousSafety, fPreviousSftOrigin );
294 if ( goodCorrection )
296 IntersectedOrRecalculatedFT = ApproxIntersecPointV;
319 G4bool usedNavigatorAF =
false;
321 NewSafety,fPreviousSafety,
326 last_AF_intersection = Intersects_AF;
337 CurrentA_PointVelocity, CurrentB_PointVelocity,
338 EndPoint,CurrentE_Point, CurrentF_Point,PointG,
340 CurrentB_PointVelocity = EndPoint;
341 CurrentE_Point = PointG;
349 validNormalAtE = validNormalLast;
354 fin_section_depth[depth] =
false;
358 G4cout <<
"G4PiF::LI> Investigating intermediate point"
360 <<
" on way to full s="
375 G4bool usedNavigatorFB =
false;
381 NewSafety,fPreviousSafety,
402 CurrentA_PointVelocity,CurrentB_PointVelocity,
403 InterMed,CurrentE_Point,CurrentF_Point,PointH,
405 CurrentA_PointVelocity = InterMed;
406 CurrentE_Point = PointH;
412 validNormalAtE= validNormalLast;
418 if( fin_section_depth[depth] )
428 if( ((Second_half)&&(depth==0)) || (first_section) )
430 there_is_no_intersection =
true;
437 substep_no_p = param_substeps+2;
443 sub_final_section =
true;
452 CurrentA_PointVelocity = CurrentB_PointVelocity;
453 CurrentB_PointVelocity = CurveEndPointVelocity;
454 SubStart_PointVelocity = CurrentA_PointVelocity;
457 CurrentB_PointVelocity,
461 restoredFullEndpoint =
true;
468 CurrentA_PointVelocity = CurrentB_PointVelocity;
470 SubStart_PointVelocity = CurrentA_PointVelocity;
473 CurrentB_PointVelocity,
476 restoredFullEndpoint =
true;
500 CurrentB_PointVelocity,
504 CurrentB_PointVelocity = newEndPointFT;
506 if ( (fin_section_depth[depth])
507 &&( first_section || ((Second_half)&&(depth==0)) ) )
509 recalculatedEndPoint =
true;
510 IntersectedOrRecalculatedFT = newEndPointFT;
514 if( curveDist < 0.0 )
517 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
518 -1.0, NewSafety, substep_no );
519 std::ostringstream message;
520 message <<
"Error in advancing propagation." <<
G4endl
521 <<
" Error in advancing propagation." <<
G4endl
522 <<
" Point A (start) is " << CurrentA_PointVelocity
524 <<
" Point B (end) is " << CurrentB_PointVelocity
526 <<
" Curve distance is " << curveDist <<
G4endl
528 <<
"The final curve point is not further along"
529 <<
" than the original!" <<
G4endl;
531 if( recalculatedEndPoint )
533 message <<
"Recalculation of EndPoint was called with fEpsStep= "
536 oldprc =
G4cerr.precision(20);
537 message <<
" Point A (Curve start) is " << CurveStartPointVelocity
539 <<
" Point B (Curve end) is " << CurveEndPointVelocity
541 <<
" Point A (Current start) is " << CurrentA_PointVelocity
543 <<
" Point B (Current end) is " << CurrentB_PointVelocity
545 <<
" Point S (Sub start) is " << SubStart_PointVelocity
547 <<
" Point E (Trial Point) is " << CurrentE_Point
549 <<
" Old Point F(Intersection) is " << CurrentF_Point
551 <<
" New Point F(Intersection) is " << ApproxIntersecPointV
553 <<
" LocateIntersection parameters are : Substep no= "
554 << substep_no << G4endl
555 <<
" Substep depth no= "<< substep_no_p <<
" Depth= "
557 <<
" Restarted no= "<< restartB <<
" Epsilon= "
560 G4cerr.precision( oldprc );
562 G4Exception(
"G4BrentLocator::EstimateIntersectionPoint()",
566 if(restoredFullEndpoint)
568 fin_section_depth[depth] = restoredFullEndpoint;
569 restoredFullEndpoint =
false;
574 #ifdef G4DEBUG_LOCATE_INTERSECTION
575 static G4int trigger_substepno_print= warn_substeps - 20 ;
577 if( substep_no >= trigger_substepno_print )
579 G4cout <<
"Difficulty in converging in "
580 <<
"G4BrentLocator::EstimateIntersectionPoint()"
582 <<
" Substep no = " << substep_no <<
G4endl;
583 if( substep_no == trigger_substepno_print )
585 printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
588 G4cout <<
" State of point A: ";
589 printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
590 -1.0, NewSafety, substep_no-1, 0);
591 G4cout <<
" State of point B: ";
592 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
593 -1.0, NewSafety, substep_no);
599 }
while ( ( ! found_approximate_intersection )
600 && ( ! there_is_no_intersection )
601 && ( substep_no_p <= param_substeps) );
603 first_section =
false;
605 if( (!found_approximate_intersection) && (!there_is_no_intersection) )
618 if ( ( did_len < fraction_done*all_len )
619 && (depth<
max_depth) && (!sub_final_section) )
624 G4double Sub_len = (all_len-did_len)/(2.);
634 SubStart_PointVelocity = CurrentA_PointVelocity;
643 NewSafety, fPreviousSafety,
644 fPreviousSftOrigin,stepLengthAB,
648 last_AF_intersection = Intersects_AB;
649 CurrentE_Point = PointGe;
650 fin_section_depth[depth]=
true;
656 validNormalAtE= validNormalAB;
667 if( (Second_half)&&(depth!=0) )
693 CurrentB_PointVelocity,
697 CurrentB_PointVelocity = newEndPointFT;
700 recalculatedEndPoint =
true;
701 IntersectedOrRecalculatedFT = newEndPointFT;
712 fPreviousSftOrigin,stepLengthAB, PointGe);
715 last_AF_intersection = Intersects_AB;
716 CurrentE_Point = PointGe;
720 validNormalAtE = validNormalAB;
724 fin_section_depth[depth]=
true;
728 }
while ( ( ! found_approximate_intersection )
729 && ( ! there_is_no_intersection )
730 && ( substep_no <= max_substeps) );
732 if( substep_no > max_no_seen )
734 max_no_seen = substep_no;
735 #ifdef G4DEBUG_LOCATE_INTERSECTION
736 if( max_no_seen > warn_substeps )
738 trigger_substepno_print = max_no_seen-20;
743 if( ( substep_no >= max_substeps)
744 && !there_is_no_intersection
745 && !found_approximate_intersection )
747 G4cout <<
"ERROR - G4BrentLocator::EstimateIntersectionPoint()" <<
G4endl
748 <<
" Start and end-point of requested Step:" <<
G4endl;
749 printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
751 G4cout <<
" Start and end-point of current Sub-Step:" <<
G4endl;
752 printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
753 -1.0, NewSafety, substep_no-1);
754 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
755 -1.0, NewSafety, substep_no);
756 std::ostringstream message;
757 message <<
"Too many substeps!" << G4endl
758 <<
" Convergence is requiring too many substeps: "
759 << substep_no << G4endl
760 <<
" Abandoning effort to intersect. " << G4endl
761 <<
" Found intersection = "
762 << found_approximate_intersection << G4endl
763 <<
" Intersection exists = "
764 << !there_is_no_intersection <<
G4endl;
765 oldprc =
G4cout.precision( 10 );
768 message <<
" Undertaken only length: " << done_len
769 <<
" out of " << full_len <<
" required." << G4endl
770 <<
" Remaining length = " << full_len - done_len;
771 G4cout.precision( oldprc );
773 G4Exception(
"G4BrentLocator::EstimateIntersectionPoint()",
776 else if( substep_no >= warn_substeps )
778 oldprc=
G4cout.precision( 10 );
779 std::ostringstream message;
780 message <<
"Many substeps while trying to locate intersection."
782 <<
" Undertaken length: "
784 <<
" - Needed: " << substep_no <<
" substeps." <<
G4endl
785 <<
" Warning level = " << warn_substeps
786 <<
" and maximum substeps = " << max_substeps;
787 G4Exception(
"G4BrentLocator::EstimateIntersectionPoint()",
789 G4cout.precision( oldprc );
791 return !there_is_no_intersection;
G4FieldTrack ReEstimateEndpoint(const G4FieldTrack &CurrentStateA, const G4FieldTrack &EstimtdEndStateB, G4double linearDistSq, G4double curveDist)
void SetPosition(G4ThreeVector nPos)
G4double GetCurveLength() const
CLHEP::Hep3Vector G4ThreeVector
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)
G4FieldTrack * ptrInterMedFT[max_depth+1]
const G4ThreeVector & GetMomentumDir() const
G4int maxNumberOfStepsForIntersection
void printStatus(const G4FieldTrack &startFT, const G4FieldTrack ¤tFT, G4double requestStep, G4double safety, G4int stepNum)
G4int maxNumberOfCallsToReIntegration_depth
G4double GetEpsilonStepFor()
G4ThreeVector GetPosition() const
G4GLOB_DLL std::ostream G4cout
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)
G4bool EstimateIntersectionPoint(const G4FieldTrack &curveStartPointTangent, const G4FieldTrack &curveEndPointTangent, const G4ThreeVector &trialPoint, G4FieldTrack &intersectPointTangent, G4bool &recalculatedEndPoint, G4double &fPreviousSafety, G4ThreeVector &fPreviousSftOrigin)
#define fPreviousSftOrigin
G4double GetDeltaIntersectionFor()
G4bool GetAdjustementOfFoundIntersection()
static const G4int max_depth
G4int maxNumberOfCallsToReIntegration
G4MagInt_Driver * GetIntegrationDriver()
G4double fiDeltaIntersection
G4ChordFinder * GetChordFinderFor()
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 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)
G4GLOB_DLL std::ostream G4cerr
G4ThreeVector GetMomentumDirection() const
G4BrentLocator(G4Navigator *theNavigator)