47 for (
G4int idepth=0; idepth<max_depth+1; idepth++ )
49 ptrInterMedFT[ idepth ] =
new G4FieldTrack( zeroV, zeroV, 0., 0., 0., 0.);
55 for (
G4int idepth=0; idepth<max_depth+1; idepth++)
57 delete ptrInterMedFT[idepth];
98 G4bool& recalculatedEndPoint,
104 G4bool found_approximate_intersection =
false;
105 G4bool there_is_no_intersection =
false;
107 G4FieldTrack CurrentA_PointVelocity = CurveStartPointVelocity;
108 G4FieldTrack CurrentB_PointVelocity = CurveEndPointVelocity;
110 G4bool validNormalAtE =
false;
113 G4FieldTrack ApproxIntersecPointV(CurveEndPointVelocity);
115 G4bool last_AF_intersection =
false;
119 G4bool first_section =
true;
120 recalculatedEndPoint =
false;
122 G4bool restoredFullEndpoint =
false;
124 G4int substep_no = 0;
130 const G4int max_substeps= 10000;
131 const G4int warn_substeps= 1000;
147 const G4int param_substeps=5;
151 G4bool Second_half =
false;
162 static const G4double tolerance = 1.0e-8 *
mm;
164 if( (TrialPoint - StartPosition).mag() < tolerance)
166 G4Exception(
"G4MultiLevelLocator::EstimateIntersectionPoint()",
168 "Intersection point F is exactly at start point A." );
178 *ptrInterMedFT[0] = CurveEndPointVelocity;
179 for (
G4int idepth=1; idepth<max_depth+1; idepth++ )
181 *ptrInterMedFT[idepth]=CurveStartPointVelocity;
186 G4bool fin_section_depth[max_depth];
187 for (
G4int idepth=0; idepth<max_depth; idepth++ )
189 fin_section_depth[idepth]=
true;
193 G4FieldTrack SubStart_PointVelocity = CurveStartPointVelocity;
197 G4int substep_no_p = 0;
198 G4bool sub_final_section =
false;
200 SubStart_PointVelocity = CurrentA_PointVelocity;
211 CurrentB_PointVelocity,
220 G4Exception(
"G4multiLevelLocator::EstimateIntersectionPoint()",
222 "Intermediate F point is past end B point" );
231 G4ThreeVector ChordEF_Vector = CurrentF_Point - CurrentE_Point;
234 G4double MomDir_dot_Norm= NewMomentumDir.
dot( NormalAtEntry ) ;
237 if( VerboseLevel > 3 )
242 MomDir_dot_ABchord = (1.0 / ABchord_length)
243 * NewMomentumDir.
dot( ChordAB );
245 ChordEF_Vector, NewMomentumDir, NormalAtEntry, validNormalAtE );
246 G4cout <<
" dot( MomentumDir, ABchord_unit ) = "
247 << MomDir_dot_ABchord <<
G4endl;
251 ( MomDir_dot_Norm >= 0.0 )
252 || (! validNormalAtE );
257 found_approximate_intersection =
true;
261 IntersectedOrRecalculatedFT = ApproxIntersecPointV;
262 IntersectedOrRecalculatedFT.
SetPosition( CurrentE_Point );
271 CurrentE_Point, CurrentF_Point, MomentumDir,
272 last_AF_intersection, IP, NewSafety,
273 previousSafety, previousSftOrigin );
274 if ( goodCorrection )
276 IntersectedOrRecalculatedFT = ApproxIntersecPointV;
299 NewSafety, previousSafety,
303 last_AF_intersection = Intersects_AF;
312 CurrentB_PointVelocity = ApproxIntersecPointV;
313 CurrentE_Point = PointG;
317 validNormalAtE = validNormalLast;
322 fin_section_depth[depth]=
false;
328 G4cout <<
"G4PiF::LI> Investigating intermediate point"
330 <<
" on way to full s="
350 NewSafety, previousSafety,
368 CurrentA_PointVelocity = ApproxIntersecPointV;
369 CurrentE_Point = PointH;
373 validNormalAtE = validNormalH;
380 if(fin_section_depth[depth])
390 if( ((Second_half)&&(depth==0)) || (first_section) )
392 there_is_no_intersection =
true;
399 substep_no_p = param_substeps+2;
405 sub_final_section =
true;
415 CurrentA_PointVelocity = CurrentB_PointVelocity;
416 CurrentB_PointVelocity = CurveEndPointVelocity;
417 SubStart_PointVelocity = CurrentA_PointVelocity;
418 restoredFullEndpoint =
true;
424 CurrentA_PointVelocity = CurrentB_PointVelocity;
425 CurrentB_PointVelocity = *ptrInterMedFT[depth];
426 SubStart_PointVelocity = CurrentA_PointVelocity;
427 restoredFullEndpoint =
true;
450 CurrentB_PointVelocity,
454 CurrentB_PointVelocity = newEndPointFT;
456 if ( (fin_section_depth[depth])
457 &&( first_section || ((Second_half)&&(depth==0)) ) )
459 recalculatedEndPoint =
true;
460 IntersectedOrRecalculatedFT = newEndPointFT;
464 if( curveDist < 0.0 )
467 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
468 -1.0, NewSafety, substep_no );
469 std::ostringstream message;
470 message <<
"Error in advancing propagation." <<
G4endl
471 <<
" Point A (start) is " << CurrentA_PointVelocity
473 <<
" Point B (end) is " << CurrentB_PointVelocity
475 <<
" Curve distance is " << curveDist <<
G4endl
477 <<
"The final curve point is not further along"
478 <<
" than the original!" <<
G4endl;
480 if( recalculatedEndPoint )
482 message <<
"Recalculation of EndPoint was called with fEpsStep= "
485 oldprc =
G4cerr.precision(20);
486 message <<
" Point A (Curve start) is " << CurveStartPointVelocity
488 <<
" Point B (Curve end) is " << CurveEndPointVelocity
490 <<
" Point A (Current start) is " << CurrentA_PointVelocity
492 <<
" Point B (Current end) is " << CurrentB_PointVelocity
494 <<
" Point S (Sub start) is " << SubStart_PointVelocity
496 <<
" Point E (Trial Point) is " << CurrentE_Point
498 <<
" Point F (Intersection) is " << ApproxIntersecPointV
500 <<
" LocateIntersection parameters are : Substep no= "
501 << substep_no << G4endl
502 <<
" Substep depth no= "<< substep_no_p <<
" Depth= "
506 G4Exception(
"G4MultiLevelLocator::EstimateIntersectionPoint()",
509 if(restoredFullEndpoint)
511 fin_section_depth[depth] = restoredFullEndpoint;
512 restoredFullEndpoint =
false;
517 #ifdef G4DEBUG_LOCATE_INTERSECTION
518 static G4int trigger_substepno_print= warn_substeps - 20 ;
520 if( substep_no >= trigger_substepno_print )
522 G4cout <<
"Difficulty in converging in "
523 <<
"G4MultiLevelLocator::EstimateIntersectionPoint():"
525 <<
" Substep no = " << substep_no <<
G4endl;
526 if( substep_no == trigger_substepno_print )
528 printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
531 G4cout <<
" State of point A: ";
532 printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
533 -1.0, NewSafety, substep_no-1, 0);
534 G4cout <<
" State of point B: ";
535 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
536 -1.0, NewSafety, substep_no);
542 }
while ( ( ! found_approximate_intersection )
543 && ( ! there_is_no_intersection )
544 && ( substep_no_p <= param_substeps) );
546 first_section =
false;
548 if( (!found_approximate_intersection) && (!there_is_no_intersection) )
560 if( ( ( did_len )<fraction_done*all_len)
561 && (depth<max_depth) && (!sub_final_section) )
566 G4double Sub_len = (all_len-did_len)/(2.);
571 *ptrInterMedFT[depth] = start;
572 CurrentB_PointVelocity = *ptrInterMedFT[depth];
576 SubStart_PointVelocity = CurrentA_PointVelocity;
585 NewSafety, previousSafety,
586 previousSftOrigin, stepLengthAB,
590 last_AF_intersection = Intersects_AB;
591 CurrentE_Point = PointGe;
592 fin_section_depth[depth]=
true;
596 validNormalAtE = validNormalAB;
607 if( (Second_half)&&(depth!=0) )
615 SubStart_PointVelocity = *ptrInterMedFT[depth];
616 CurrentA_PointVelocity = *ptrInterMedFT[depth];
617 CurrentB_PointVelocity = *ptrInterMedFT[depth-1];
632 CurrentB_PointVelocity,
636 CurrentB_PointVelocity = newEndPointFT;
639 recalculatedEndPoint =
true;
640 IntersectedOrRecalculatedFT = newEndPointFT;
650 previousSftOrigin, stepLengthAB,
654 last_AF_intersection = Intersects_AB;
655 CurrentE_Point = PointGe;
659 validNormalAtE = validNormalAB;
662 fin_section_depth[depth]=
true;
666 }
while ( ( ! found_approximate_intersection )
667 && ( ! there_is_no_intersection )
668 && ( substep_no <= max_substeps) );
670 if( substep_no > max_no_seen )
672 max_no_seen = substep_no;
673 #ifdef G4DEBUG_LOCATE_INTERSECTION
674 if( max_no_seen > warn_substeps )
676 trigger_substepno_print = max_no_seen-20;
681 if( ( substep_no >= max_substeps)
682 && !there_is_no_intersection
683 && !found_approximate_intersection )
685 G4cout <<
"ERROR - G4MultiLevelLocator::EstimateIntersectionPoint()"
687 <<
" Start and end-point of requested Step:" <<
G4endl;
688 printStatus( CurveStartPointVelocity, CurveEndPointVelocity,
690 G4cout <<
" Start and end-point of current Sub-Step:" <<
G4endl;
691 printStatus( CurrentA_PointVelocity, CurrentA_PointVelocity,
692 -1.0, NewSafety, substep_no-1);
693 printStatus( CurrentA_PointVelocity, CurrentB_PointVelocity,
694 -1.0, NewSafety, substep_no);
695 std::ostringstream message;
696 message <<
"Too many substeps!" << G4endl
697 <<
" Convergence is requiring too many substeps: "
698 << substep_no << G4endl
699 <<
" Abandoning effort to intersect. " << G4endl
700 <<
" Found intersection = "
701 << found_approximate_intersection << G4endl
702 <<
" Intersection exists = "
703 << !there_is_no_intersection <<
G4endl;
705 #ifdef FUTURE_CORRECTION
708 if ( ! found_approximate_intersection )
710 recalculatedEndPoint =
true;
713 IntersectedOrRecalculatedFT = CurrentA_PointVelocity;
715 message <<
" Did not convergence after " << substep_no
716 <<
" substeps." << G4endl
717 <<
" The endpoint was adjused to pointA resulting"
719 <<
" from the last substep: " << CurrentA_PointVelocity
724 oldprc =
G4cout.precision( 10 );
727 message <<
" Undertaken only length: " << done_len
728 <<
" out of " << full_len <<
" required." << G4endl
729 <<
" Remaining length = " << full_len - done_len;
730 G4cout.precision( oldprc );
732 G4Exception(
"G4MultiLevelLocator::EstimateIntersectionPoint()",
735 else if( substep_no >= warn_substeps )
737 oldprc =
G4cout.precision( 10 );
738 std::ostringstream message;
739 message <<
"Many substeps while trying to locate intersection."
741 <<
" Undertaken length: "
743 <<
" - Needed: " << substep_no <<
" substeps." <<
G4endl
744 <<
" Warning level = " << warn_substeps
745 <<
" and maximum substeps = " << max_substeps;
746 G4Exception(
"G4MultiLevelLocator::EstimateIntersectionPoint()",
748 G4cout.precision( oldprc );
750 return !there_is_no_intersection;
G4FieldTrack ReEstimateEndpoint(const G4FieldTrack &CurrentStateA, const G4FieldTrack &EstimtdEndStateB, G4double linearDistSq, G4double curveDist)
void SetPosition(G4ThreeVector nPos)
G4double GetCurveLength() const
double dot(const Hep3Vector &) const
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)
G4double GetEpsilonStepFor()
G4ThreeVector GetPosition() const
G4GLOB_DLL std::ostream G4cout
G4MultiLevelLocator(G4Navigator *theNavigator)
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 printStatus(const G4FieldTrack &startFT, const G4FieldTrack ¤tFT, G4double requestStep, G4double safety, G4int step)
G4bool GetAdjustementOfFoundIntersection()
G4MagInt_Driver * GetIntegrationDriver()
G4double fiDeltaIntersection
G4ChordFinder * GetChordFinderFor()
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