49 fUseNormalCorrection(false),
51 fiNavigator(theNavigator),
54 fiDeltaIntersection(-1.0),
83 std::ostringstream os;
110 if( ((stepNo == 0) && (verboseLevel <3)) || (verboseLevel >= 3) )
112 oldprc = os.precision(4);
113 os << std::setw( 6) <<
" "
114 << std::setw( 25) <<
" Current Position and Direction" <<
" "
116 os << std::setw( 5) <<
"Step#"
117 << std::setw(10) <<
" s " <<
" "
118 << std::setw(10) <<
"X(mm)" <<
" "
119 << std::setw(10) <<
"Y(mm)" <<
" "
120 << std::setw(10) <<
"Z(mm)" <<
" "
121 << std::setw( 7) <<
" N_x " <<
" "
122 << std::setw( 7) <<
" N_y " <<
" "
123 << std::setw( 7) <<
" N_z " <<
" " ;
124 os << std::setw( 7) <<
" Delta|N|" <<
" "
125 << std::setw( 9) <<
"StepLen" <<
" "
126 << std::setw(12) <<
"StartSafety" <<
" "
127 << std::setw( 9) <<
"PhsStep" <<
" ";
129 os.precision(oldprc);
131 if((stepNo == 0) && (verboseLevel <=3))
135 printStatus( StartFT, StartFT, -1.0, safety, -1, os, verboseLevel);
137 if( verboseLevel <= 3 )
141 os << std::setw( 4) << stepNo <<
" ";
145 os << std::setw( 5) <<
"Start" ;
147 oldprc = os.precision(8);
149 os << std::setw(10) << CurrentPosition.
x() <<
" "
150 << std::setw(10) << CurrentPosition.
y() <<
" "
151 << std::setw(10) << CurrentPosition.
z() <<
" ";
153 os << std::setw( 7) << CurrentUnitVelocity.
x() <<
" "
154 << std::setw( 7) << CurrentUnitVelocity.
y() <<
" "
155 << std::setw( 7) << CurrentUnitVelocity.
z() <<
" ";
160 os << std::setw( 9) << step_len <<
" ";
161 os << std::setw(12) << safety <<
" ";
162 if( requestStep != -1.0 )
164 os << std::setw( 9) << requestStep <<
" ";
168 os << std::setw( 9) <<
"Init/NotKnown" <<
" ";
171 os.precision(oldprc);
177 os <<
"Step taken was " << step_len
178 <<
" out of PhysicalStep= " << requestStep <<
G4endl;
179 os <<
"Final safety is: " << safety <<
G4endl;
180 os <<
"Chord length = " << (CurrentPosition-StartPosition).mag()
203 const G4int no_trials=20;
211 G4double advanceLength= endCurveLen - currentCurveLen ;
218 goodAdvance= integrDriver->AccurateAdvance(newEndPoint, advanceLength,
222 while( !goodAdvance && (++itrial < no_trials) );
226 retEndPoint = newEndPoint;
230 retEndPoint = EstimatedEndStateB;
236 const G4String MethodName(
"G4VIntersectionLocator::ReEstimateEndpoint()");
239 G4int latest_good_trials = 0;
244 G4cout << MethodName <<
" called - goodAdv= " << goodAdvance
245 <<
" trials = " << itrial
246 <<
" previous good= " << latest_good_trials
249 latest_good_trials = 0;
253 latest_good_trials++;
264 G4cout << MethodName <<
"> AccurateAdvance failed " ;
265 G4cout <<
" in " << itrial <<
" integration trials/steps. " <<
G4endl;
266 G4cout <<
" It went only " << lengthDone <<
" instead of " << curveDist
267 <<
" -- a difference of " << curveDist - lengthDone <<
G4endl;
268 G4cout <<
" ReEstimateEndpoint> Reset endPoint to original value!"
274 static G4int noInaccuracyWarnings = 0;
275 G4int maxNoWarnings = 10;
276 if ( (noInaccuracyWarnings < maxNoWarnings )
282 message.precision(12);
283 message <<
" Integration inaccuracy requires"
284 <<
" an adjustment in the step's endpoint." <<
G4endl
285 <<
" Two mid-points are further apart than their"
286 <<
" curve length difference" <<
G4endl
287 <<
" Dist = " << linearDist
288 <<
" curve length = " << curveDist <<
G4endl;
289 message <<
" Correction applied is " << move.
mag() << G4endl
290 <<
" Old Estimated B position= "
292 <<
" Recalculated Position= "
294 <<
" Change ( new - old ) = " << move;
295 G4Exception(
"G4VIntersectionLocator::ReEstimateEndpoint()",
306 sumCorrectionsSq += (EstimatedEndStateB.
GetPosition() -
334 G4bool recalculated=
false;
342 && (curveDist*curveDist *(1.0+2.0*
fiEpsilonStep ) < linDistSq ) )
359 newEndPointFT= CurrentStartA;
363 G4Exception(
"G4MultiLevelLocator::EstimateIntersectionPoint()",
365 "A & B are at equal distance in 2nd half. A & B will coincide." );
371 if( curveDist < 0.0 )
410 if( (pLogical != 0) && ( (pSolid=pLogical->
GetSolid()) !=0 ) )
424 G4cerr <<
"PROBLEM in G4VIntersectionLocator::GetLocalSurfaceNormal."
427 G4cerr <<
" at trial local point " << CurrentE_Point <<
G4endl;
454 G4bool goodAdjust=
false, Intersects_FP=
false, validNormal=
false;
459 if(!validNormal) {
return false; }
469 G4Exception(
"G4VIntersectionLocator::AdjustmentOfFoundIntersection()",
471 "No intersection. Parallels lines!");
474 lambda =- Normal.
dot(CurrentF_Point-CurrentE_Point)/n_d_m;
478 NewPoint = CurrentF_Point+lambda*MomentumDir;
482 dist = std::abs(lambda);
494 fPreviousSafety, fPreviousSftOrigin,
495 stepLengthFP, Point_G );
502 Intersects_FP =
IntersectChord( CurrentF_Point, NewPoint, NewSafety,
503 fPreviousSafety, fPreviousSftOrigin,
504 stepLengthFP, Point_G );
509 IntersectionPoint = Point_G;
526 G4ThreeVector NormalAtEntryLast, NormalAtEntryGlobal, diffNormals;
543 message <<
"PROBLEM: Normal is not unit - magnitude = "
545 message <<
" at trial intersection point " << CurrentInt_Point <<
G4endl;
546 message <<
" Obtained from Get *Last* Surface Normal.";
547 G4Exception(
"G4VIntersectionLocator::GetSurfaceNormal()",
552 if( validNormalLast )
554 NormalAtEntry=NormalAtEntryLast;
556 validNormal = validNormalLast;
558 return NormalAtEntry;
577 if( validNormal && ( std::fabs(globalNormal.
mag2() - 1.0) >
perThousand ) )
580 message <<
"**************************************************************"
582 message <<
" Bad Normal in G4VIntersectionLocator::GetGlobalSurfaceNormal "
584 message <<
" * Constituents: " <<
G4endl;
585 message <<
" Local Normal= " << localNormal <<
G4endl;
586 message <<
" Transform: " << G4endl
589 <<
" Net Rotation = " << localToGlobal.
NetRotation()
591 message <<
" * Result: " <<
G4endl;
592 message <<
" Global Normal= " << localNormal <<
G4endl;
593 message <<
"**************************************************************";
594 G4Exception(
"G4VIntersectionLocator::GetGlobalSurfaceNormal()",
608 G4bool& normalIsValid)
const
613 normalIsValid= validNorm;
630 G4double MomDir_dot_Norm = NewMomentumDir.
dot( NormalAtEntry ) ;
632 MomDir_dot_ABchord= (1.0 / ABchord_length) * NewMomentumDir.
dot( ChordAB_v );
634 std::ostringstream outStream;
636 << std::setw(6) <<
" Step# "
637 << std::setw(17) <<
" |ChordEF|(mag)" <<
" "
638 << std::setw(18) <<
" uMomentum.Normal" <<
" "
639 << std::setw(18) <<
" uMomentum.ABdir " <<
" "
640 << std::setw(16) <<
" AB-dist " <<
" "
641 <<
" Chord Vector (EF) "
643 outStream.precision(7);
645 <<
" " << std::setw(5) << step_no
646 <<
" " << std::setw(18) << ChordEF_v.
mag()
647 <<
" " << std::setw(18) << MomDir_dot_Norm
648 <<
" " << std::setw(18) << MomDir_dot_ABchord
649 <<
" " << std::setw(12) << ABchord_length
653 <<
" MomentumDir= " <<
" " << NewMomentumDir
654 <<
" Normal at Entry E= " << NormalAtEntry
655 <<
" AB chord = " << ChordAB_v
657 G4cout << outStream.str();
662 message <<
"Normal is not unit - mag= " << NormalAtEntry.
mag() << G4endl
663 <<
" ValidNormalAtE = " << validNormal;
664 G4Exception(
"G4VIntersectionLocator::ReportTrialStep()",
687 MethodName(
"G4VIntersectionLocator::LocateGlobalPointWithinVolumeAndCheck()");
709 message <<
"Position located "
710 << ( inMother ==
kSurface ?
" on Surface " :
" outside " )
711 <<
"expected volume" <<
G4endl
712 <<
" Safety (from Outside) = "
714 G4Exception(
"G4VIntersectionLocator::LocateGlobalPointWithinVolumeAndCheck()",
723 if( (nextPhysical != motherPhys)
724 || (nextPhysical->
GetCopyNo() != motherCopyNo )
727 G4Exception(
"G4VIntersectionLocator::LocateGlobalPointWithinVolumeAndCheck()",
729 "Position located outside expected volume.");
756 message <<
"Failed point location." <<
G4endl
757 <<
" Code Location info: " << CodeLocationInfo;
758 G4Exception(
"G4VIntersectionLocator::LocateGlobalPointWithinVolumeCheckAndReport()",
791 G4int verboseLevel= 5;
794 -1.0, NewSafety, substep_no, msg, verboseLevel );
795 msg <<
"Error in advancing propagation." <<
G4endl
796 <<
" Point A (start) is " << A_PtVel <<
G4endl
797 <<
" Point B (end) is " << B_PtVel <<
G4endl
798 <<
" Curve distance is " << curveDist <<
G4endl
800 <<
"The final curve point is not further along"
801 <<
" than the original!" <<
G4endl;
802 msg <<
" Value of fEpsStep= " << epsStep <<
G4endl;
804 G4int oldprc = msg.precision(20);
805 msg <<
" Point A (Curve start) is " << StartPointVel << G4endl
806 <<
" Point B (Curve end) is " << EndPointVel << G4endl
807 <<
" Point A (Current start) is " << A_PtVel << G4endl
808 <<
" Point B (Current end) is " << B_PtVel << G4endl
809 <<
" Point S (Sub start) is " << SubStart_PtVel
810 <<
" Point E (Trial Point) is " << E_Point << G4endl
811 <<
" Point F (Intersection) is " << ApproxIntersecPointV
813 <<
" LocateIntersection parameters are : " << G4endl
814 <<
" Substep no (total) = " << substep_no << G4endl
815 <<
" Substep (depth= " << depth << substep_no_p;
816 msg.precision(oldprc);
833 oss <<
"ReportProgress: Current status of intersection search: " <<
G4endl;
834 if( depth > 0 ) oss <<
" Depth= " << depth;
835 oss <<
" Substep no = " << substep_no <<
G4endl;
836 G4int verboseLevel = 5;
841 printStatus( StartPointVel, EndPointVel, -1.0, -1.0, -1,
843 oss <<
" * Start and end-point of requested Step:" <<
G4endl;
844 oss <<
" ** State of point A: ";
845 printStatus( A_PtVel, A_PtVel, -1.0, safetyPrev, substep_no-1,
847 oss <<
" ** State of point B: ";
848 printStatus( A_PtVel, B_PtVel, -1.0, safetyLast, substep_no,
861 unsigned long int numCalls )
872 if( (TrialPoint - StartPosition).mag2() < tolerance*tolerance)
877 G4cout <<
"Intersection F == start A in " << MethodName;
880 G4cout <<
" Start-Trial: " << TrialPoint - StartPosition;
881 G4cout <<
" Start-last: " << StartPosition - lastStart;
883 if( (StartPosition - lastStart).mag() < tolerance )
888 G4cout <<
" { Unmoved: " <<
" still#= " << numStill
889 <<
" total # = " << numUnmoved <<
" } - ";
895 G4cout <<
" Occured: " << ++occurredOnTop;
896 G4cout <<
" out of total calls= " << numCalls;
898 lastStart = StartPosition;
G4VIntegrationDriver * GetIntegrationDriver()
void ReportTrialStep(G4int step_no, const G4ThreeVector &ChordAB_v, const G4ThreeVector &ChordEF_v, const G4ThreeVector &NewMomentumDir, const G4ThreeVector &NormalAtEntry, G4bool validNormal)
static constexpr double perThousand
CLHEP::Hep3Vector G4ThreeVector
G4ThreeVector GetLocalSurfaceNormal(const G4ThreeVector &CurrentE_Point, G4bool &validNormal)
void LocateGlobalPointWithinVolumeCheckAndReport(const G4ThreeVector &pos, const G4String &CodeLocationInfo, G4int CheckMode)
G4LogicalVolume * GetLogicalVolume() const
G4bool IsCheckModeActive() const
G4bool LocateGlobalPointWithinVolumeAndCheck(const G4ThreeVector &pos)
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)
virtual void LocateGlobalPointWithinVolume(const G4ThreeVector &position)
G4TouchableHistory * CreateTouchableHistory() const
void message(RunManager *runmanager)
double dot(const Hep3Vector &) const
virtual G4double DistanceToIn(const G4ThreeVector &p, const G4ThreeVector &v) const =0
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)
virtual G4ThreeVector SurfaceNormal(const G4ThreeVector &p) const =0
G4ThreeVector GetLastSurfaceNormal(const G4ThreeVector &intersectPoint, G4bool &validNormal) const
void CheckMode(G4bool mode)
virtual EInside Inside(const G4ThreeVector &p) const =0
const G4AffineTransform & GetGlobalToLocalTransform() const
const G4AffineTransform GetLocalToGlobalTransform() const
G4bool IntersectChord(const G4ThreeVector &StartPointA, const G4ThreeVector &EndPointB, G4double &NewSafety, G4double &PreviousSafety, G4ThreeVector &PreviousSftOrigin, G4double &LinearStepLength, G4ThreeVector &IntersectionPoint, G4bool *calledNavigator=0)
G4ThreeVector GetPosition() const
G4VSolid * GetSolid(G4int depth=0) const
G4ThreeVector GetGlobalSurfaceNormal(const G4ThreeVector &CurrentE_Point, G4bool &validNormal)
const G4NavigationHistory * GetHistory() const
#define fPreviousSftOrigin
G4VPhysicalVolume * GetVolume(G4int depth=0) const
G4Navigator * GetNavigatorFor()
void SetWorldVolume(G4VPhysicalVolume *pWorld)
virtual ~G4VIntersectionLocator()
const G4AffineTransform & GetTopTransform() const
G4ThreeVector GetMomentum() const
G4FieldTrack ReEstimateEndpoint(const G4FieldTrack &CurrentStateA, const G4FieldTrack &EstimtdEndStateB, G4double linearDistSq, G4double curveDist)
virtual G4int GetCopyNo() const =0
G4GLOB_DLL std::ostream G4cerr
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)
virtual G4TouchableHistoryHandle CreateTouchableHistoryHandle() const
G4Navigator * fHelpingNavigator
G4TouchableHistory * fpTouchable
G4Navigator * fiNavigator
void ReportImmediateHit(const char *MethodName, const G4ThreeVector &StartPosition, const G4ThreeVector &TrialPoint, double tolerance, unsigned long int numCalls)
const G4ThreeVector & GetMomentumDir() const
G4GLOB_DLL std::ostream G4cout
virtual G4double DistanceToOut(const G4ThreeVector &p, const G4ThreeVector &v, const G4bool calcNorm=false, G4bool *validNorm=0, G4ThreeVector *n=0) const =0
G4VSolid * GetSolid() const
virtual G4VPhysicalVolume * LocateGlobalPointAndSetup(const G4ThreeVector &point, const G4ThreeVector *direction=0, const G4bool pRelativeSearch=true, const G4bool ignoreDirection=true)
static G4GeometryTolerance * GetInstance()
G4double GetEpsilonStepFor()
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)
static constexpr double perThousand
G4double GetSurfaceTolerance() const
G4double GetCurveLength() const
G4VIntersectionLocator(G4Navigator *theNavigator)
G4ThreeVector GetSurfaceNormal(const G4ThreeVector &CurrentInt_Point, G4bool &validNormal)
void SetCheckMode(G4bool value)
virtual G4ThreeVector GetGlobalExitNormal(const G4ThreeVector &point, G4bool *valid)