64 fMax_loop_count(1000),
65 fUseSafetyForOptimisation(true),
66 fZeroStepThreshold( 0.0 ),
67 fDetectorFieldMgr(detectorFieldMgr),
68 fpTrajectoryFilter( 0 ),
69 fNavigator(theNavigator),
70 fCurrentFieldMgr(detectorFieldMgr),
74 fParticleIsLooping(false),
78 fFirstStepInVolume(true),
79 fLastStepInVolume(true),
97 G4cout <<
" PiF: Zero Step Threshold set to "
100 G4cout <<
" PiF: Value of kCarTolerance = "
154 const char* methodName=
"G4PropagatorInField::ComputeStep";
174 G4cout <<
" Starting FT: " << pFieldTrack;
175 G4cout <<
" Requested length = " << CurrentProposedStepLength <<
G4endl;
187 G4double TruePathLength = CurrentProposedStepLength;
191 G4bool first_substep =
true;
216 G4double trialProposedStep = 1.e2 * ( 10.0 *
cm +
218 GetSolid()->DistanceToOut(StartPointA, VelocityUnit) );
219 CurrentProposedStepLength=
std::min( trialProposedStep,
226 if( epsilon < epsilonMin ) epsilon = epsilonMin;
227 if( epsilon > epsilonMax ) epsilon = epsilonMax;
254 decreaseFactor= 0.25;
263 decreaseFactor = 0.35;
267 decreaseFactor= 0.75;
271 stepTrial *= decreaseFactor;
278 G4cerr <<
" " << methodName
279 <<
" Decreasing step after " <<
fNoZeroStep <<
" zero steps "
280 <<
" - in volume " << pPhysVol;
282 G4cerr <<
" with name " << pPhysVol->GetName();
284 G4cerr <<
" i.e. *unknown* volume.";
287 stepTrial, pFieldTrack);
290 if( stepTrial == 0.0 )
293 message <<
"Particle abandoned due to lack of progress in field."
295 <<
" Properties : " << pFieldTrack <<
G4endl
296 <<
" Attempting a zero step = " << stepTrial <<
G4endl
297 <<
" while attempting to progress after " <<
fNoZeroStep
298 <<
" trial steps. Will abandon step.";
303 if( stepTrial < CurrentProposedStepLength )
304 CurrentProposedStepLength = stepTrial;
308 G4int do_loop_count = 0;
318 G4cout <<
" PiF: Calling Nav/Locate Global Point within-Volume "
326 h_TrialStepSize = CurrentProposedStepLength - StepTaken;
347 NewSafety, LinearStepLength,
348 InterSectionPointE );
354 currentSafety = NewSafety;
363 G4bool recalculatedEndPt=
false;
366 EstimateIntersectionPoint( SubStepStartState, CurrentState,
367 InterSectionPointE, IntersectPointVelct_G,
370 intersects = found_intersection;
371 if( found_intersection )
374 StepTaken = TruePathLength = IntersectPointVelct_G.
GetCurveLength()
382 if( recalculatedEndPt )
388 G4bool shortEnd = endAchieved
397 CurrentState= IntersectPointVelct_G;
398 s_length_taken = stepAchieved;
408 StepTaken += s_length_taken;
414 first_substep =
false;
420 G4cout <<
" Above 'Severe Action' threshold -- for Zero steps. ";
422 G4cout <<
" Above 'action' threshold -- for Zero steps. ";
425 CurrentState, CurrentProposedStepLength,
426 NewSafety, do_loop_count, pPhysVol );
432 G4cout <<
" G4PropagatorInField::ComputeStep(): " <<
G4endl
433 <<
" Difficult track - taking many sub steps." <<
G4endl;
434 printStatus( SubStepStartState, SubStepStartState, CurrentProposedStepLength,
435 NewSafety, 0, pPhysVol );
437 printStatus( SubStepStartState, CurrentState, CurrentProposedStepLength,
438 NewSafety, do_loop_count, pPhysVol );
444 }
while( (!intersects )
458 CurrentProposedStepLength, methodName,
468 TruePathLength = StepTaken;
486 message <<
"Curve length mis-match between original state "
487 <<
"and proposed endpoint of propagation." <<
G4endl
488 <<
" The curve length of the endpoint should be: "
490 <<
" and it is instead: "
492 <<
" A difference of: "
495 <<
" Original state = " << OriginalState <<
G4endl
501 if( TruePathLength+
kCarTolerance >= CurrentProposedStepLength )
526 return TruePathLength;
551 if( ((stepNo == 0) && (verboseLevel <3)) || (verboseLevel >= 3) )
553 oldprec =
G4cout.precision(4);
557 G4cout << std::setw( 5) <<
"Step#"
558 << std::setw(10) <<
" s " <<
" "
559 << std::setw(10) <<
"X(mm)" <<
" "
560 << std::setw(10) <<
"Y(mm)" <<
" "
561 << std::setw(10) <<
"Z(mm)" <<
" "
562 << std::setw( 7) <<
" N_x " <<
" "
563 << std::setw( 7) <<
" N_y " <<
" "
564 << std::setw( 7) <<
" N_z " <<
" " ;
565 G4cout << std::setw( 7) <<
" Delta|N|" <<
" "
566 << std::setw( 9) <<
"StepLen" <<
" "
567 << std::setw(12) <<
"StartSafety" <<
" "
568 << std::setw( 9) <<
"PhsStep" <<
" ";
570 {
G4cout << std::setw(18) <<
"NextVolume" <<
" "; }
571 G4cout.precision(oldprec);
574 if((stepNo == 0) && (verboseLevel <=3))
578 printStatus( StartFT, StartFT, -1.0, safety, -1, startVolume);
580 if( verboseLevel <= 3 )
583 {
G4cout << std::setw( 4) << stepNo <<
" "; }
585 {
G4cout << std::setw( 5) <<
"Start" ; }
586 oldprec =
G4cout.precision(8);
589 G4cout << std::setw(10) << CurrentPosition.
x() <<
" "
590 << std::setw(10) << CurrentPosition.
y() <<
" "
591 << std::setw(10) << CurrentPosition.
z() <<
" ";
593 G4cout << std::setw( 7) << CurrentUnitVelocity.
x() <<
" "
594 << std::setw( 7) << CurrentUnitVelocity.
y() <<
" "
595 << std::setw( 7) << CurrentUnitVelocity.
z() <<
" ";
599 G4cout << std::setw( 9) << step_len <<
" ";
600 G4cout << std::setw(12) << safety <<
" ";
601 if( requestStep != -1.0 )
602 {
G4cout << std::setw( 9) << requestStep <<
" "; }
604 {
G4cout << std::setw( 9) <<
"Init/NotKnown" <<
" "; }
605 if( startVolume != 0)
606 {
G4cout << std::setw(12) << startVolume->
GetName() <<
" "; }
607 G4cout.precision(oldprec);
614 G4cout <<
"Step taken was " << step_len
615 <<
" out of PhysicalStep = " << requestStep <<
G4endl;
617 G4cout <<
"Chord length = " << (CurrentPosition-StartPosition).mag()
635 G4cout <<
" " << std::setw(12) <<
" PiF: NoZeroStep "
636 <<
" " << std::setw(20) <<
" CurrentProposed len "
637 <<
" " << std::setw(18) <<
" Full_curvelen_last"
638 <<
" " << std::setw(18) <<
" last proposed len "
639 <<
" " << std::setw(18) <<
" decrease factor "
640 <<
" " << std::setw(15) <<
" step trial "
644 <<
" " << std::setw(20) << CurrentProposedStepLength
647 <<
" " << std::setw(18) << decreaseFactor
648 <<
" " << std::setw(15) << stepTrial
650 G4cout.precision( iprec );
663 std::vector<G4ThreeVector>*
696 0.0,0.0,0.0,0.0,0.0);
710 if( pCurrentPhysicalVolume)
720 if( pRegionFieldMgr )
721 currentFieldMgr= pRegionFieldMgr;
727 currentFieldMgr = localFieldMgr;
736 return currentFieldMgr;
759 const char* methodName,
764 G4double fraction = StepTaken / StepRequested;
765 message <<
" Unfinished integration of track (likely looping particle) "
766 <<
" of momentum " << momentumVec <<
" ( magnitude = " << momentumVec.
mag() <<
" ) "
768 <<
" after " << count <<
" field substeps "
769 <<
" totaling " << std::setprecision(12) << StepTaken /
mm <<
" mm "
770 <<
" out of requested step " << std::setprecision(12) << StepRequested /
mm <<
" mm ";
771 message <<
" a fraction of ";
773 if( fraction > 0.99 )
776 if (fraction > 0.97 )
778 message << std::setprecision(prec)
779 << 100. * StepTaken / StepRequested <<
" % " <<
G4endl ;
782 message <<
" in volume " << pPhysVol->
GetName() ;
785 message <<
" with material " << material->
GetName()
792 message <<
" in unknown (null) volume. " ;
803 message <<
"Particle is stuck; it will be killed." <<
G4endl
804 <<
" Zero progress for " << noZeroSteps <<
" attempted steps."
806 <<
" Proposed Step is " << proposedStep
807 <<
" but Step Taken is "<< lastTriedStep <<
G4endl;
809 message <<
" in volume " << physVol->
GetName() ;
811 message <<
" in unknown or null volume. " ;
G4ThreeVector fPreviousSftOrigin
G4VIntegrationDriver * GetIntegrationDriver()
G4double AdvanceChordLimited(G4FieldTrack &yCurrent, G4double stepInitial, G4double epsStep_Relative, const G4ThreeVector &latestSafetyOrigin, G4double lasestSafetyRadius)
static constexpr double micrometer
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4int fAbandonThreshold_NoZeroSteps
void PrintStepLengthDiagnostic(G4double currentProposedStepLength, G4double decreaseFactor, G4double stepTrial, const G4FieldTrack &aFieldTrack)
CLHEP::Hep3Vector G4ThreeVector
G4double fLast_ProposedStepLength
G4LogicalVolume * GetLogicalVolume() const
static const G4double kInfinity
G4double GetDeltaOneStep() const
G4double fFull_CurveLen_of_LastAttempt
G4VIntersectionLocator * fIntersectionLocator
virtual void SetVerboseLevel(G4int level)=0
std::vector< G4ThreeVector > * GimmeThePointsAndForgetThem()
static constexpr double mm
virtual void LocateGlobalPointWithinVolume(const G4ThreeVector &position)
void SetEpsilonStepFor(G4double EpsilonStep)
void ReportLoopingParticle(G4int count, G4double StepTaken, G4double stepRequest, const char *methodName, G4ThreeVector momentumVec, G4VPhysicalVolume *physVol)
G4double fLargestAcceptableStep
G4double GetDeltaIntersection() const
void message(RunManager *runmanager)
G4FieldManager * fDetectorFieldMgr
G4bool fParticleIsLooping
void CreateNewTrajectorySegment()
void SetTrajectoryFilter(G4VCurvedTrajectoryFilter *filter)
G4double ComputeStep(G4FieldTrack &pFieldTrack, G4double pCurrentProposedStepLength, G4double &pNewSafety, G4VPhysicalVolume *pPhysVol=0)
G4FieldTrack End_PointAndTangent
G4bool IntersectChord(const G4ThreeVector &StartPointA, const G4ThreeVector &EndPointB, G4double &NewSafety, G4double &LinearStepLength, G4ThreeVector &IntersectionPoint)
G4VCurvedTrajectoryFilter * fpTrajectoryFilter
G4Material * GetMaterial() const
G4FieldManager * GetFieldManager() const
void SetSafetyParametersFor(G4bool UseSafety)
const G4String & GetName() const
G4ThreeVector GetPosition() const
G4FieldManager * FindAndSetFieldManager(G4VPhysicalVolume *pCurrentPhysVol)
G4bool fFirstStepInVolume
G4double GetMaximumEpsilonStep() const
G4int fSevereActionThreshold_NoZeroSteps
G4ThreeVector GetMomentum() const
static constexpr double perMillion
std::vector< G4ThreeVector > * GimmeTrajectoryVectorAndForgetIt() const
void RefreshIntersectionLocator()
G4FieldManager * fCurrentFieldMgr
G4GLOB_DLL std::ostream G4cerr
G4PropagatorInField(G4Navigator *theNavigator, G4FieldManager *detectorFieldMgr, G4VIntersectionLocator *vLocator=0)
G4bool fUseSafetyForOptimisation
double epsilon(double density, double temperature)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
static constexpr double meter
void SetEpsilonStep(G4double newEps)
const G4ThreeVector & GetMomentumDir() const
virtual void TakeIntermediatePoint(G4ThreeVector newPoint)=0
G4double fZeroStepThreshold
G4double GetMinimumEpsilonStep() const
static constexpr double centimeter
static constexpr double cm
G4GLOB_DLL std::ostream G4cout
G4FieldManager * GetFieldManager() const
G4int fActionThreshold_NoZeroSteps
void SetDeltaIntersectionFor(G4double deltaIntersection)
void ReportStuckParticle(G4int noZeroSteps, G4double proposedStep, G4double lastTriedStep, G4VPhysicalVolume *physVol)
static G4GeometryTolerance * GetInstance()
G4ChordFinder * GetChordFinder()
void SetChordFinderFor(G4ChordFinder *fCFinder)
void printStatus(const G4FieldTrack &startFT, const G4FieldTrack ¤tFT, G4double requestStep, G4double safety, G4int step, G4VPhysicalVolume *startVolume)
G4int SetVerboseLevel(G4int verbose)
G4VPhysicalVolume * GetWorldVolume() const
void ClearPropagatorState()
static constexpr double millimeter
G4Region * GetRegion() const
G4double GetSurfaceTolerance() const
G4double GetCurveLength() const
const G4String & GetName() const
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
static constexpr double gram