63 fMax_loop_count(1000),
64 fUseSafetyForOptimisation(true),
65 fZeroStepThreshold( 0.0 ),
66 fDetectorFieldMgr(detectorFieldMgr),
67 fpTrajectoryFilter( 0 ),
68 fNavigator(theNavigator),
69 fCurrentFieldMgr(detectorFieldMgr),
73 fParticleIsLooping(false),
78 else { fEpsilonStep= 1.0e-5; }
79 fActionThreshold_NoZeroSteps = 2;
80 fSevereActionThreshold_NoZeroSteps = 10;
81 fAbandonThreshold_NoZeroSteps = 50;
82 fFull_CurveLen_of_LastAttempt = -1;
83 fLast_ProposedStepLength = -1;
84 fLargestAcceptableStep = 1000.0 *
meter;
92 G4cout <<
" PiF: Zero Step Threshold set to "
95 G4cout <<
" PiF: Value of kCarTolerance = "
103 fAllocatedLocator=
true;
105 fIntersectionLocator=vLocator;
106 fAllocatedLocator=
false;
113 if(fAllocatedLocator)
delete fIntersectionLocator;
139 if(CurrentProposedStepLength<kCarTolerance)
146 if (fpTrajectoryFilter)
154 G4double TruePathLength = CurrentProposedStepLength;
158 G4bool first_substep =
true;
161 fParticleIsLooping =
false;
181 if( CurrentProposedStepLength >= fLargestAcceptableStep )
187 G4double trialProposedStep = 1.e2 * ( 10.0 *
cm +
189 GetSolid()->DistanceToOut(StartPointA, VelocityUnit) );
190 CurrentProposedStepLength=
std::min( trialProposedStep,
191 fLargestAcceptableStep );
193 epsilon = fCurrentFieldMgr->
GetDeltaOneStep() / CurrentProposedStepLength;
197 if( epsilon < epsilonMin ) epsilon = epsilonMin;
198 if( epsilon > epsilonMax ) epsilon = epsilonMax;
206 if( fNoZeroStep > fActionThreshold_NoZeroSteps )
210 stepTrial= fFull_CurveLen_of_LastAttempt;
211 if( (stepTrial <= 0.0) && (fLast_ProposedStepLength > 0.0) )
212 stepTrial= fLast_ProposedStepLength;
215 if( (fNoZeroStep < fSevereActionThreshold_NoZeroSteps)
216 && (stepTrial > 100.0*fZeroStepThreshold) )
220 decreaseFactor= 0.25;
228 if( stepTrial > 100.0*fZeroStepThreshold )
229 decreaseFactor = 0.35;
230 else if( stepTrial > 30.0*fZeroStepThreshold )
232 else if( stepTrial > 10.0*fZeroStepThreshold )
233 decreaseFactor= 0.75;
237 stepTrial *= decreaseFactor;
240 G4cout <<
" G4PropagatorInField::ComputeStep(): " <<
G4endl
241 <<
" Decreasing step - "
242 <<
" decreaseFactor= " << std::setw(8) << decreaseFactor
243 <<
" stepTrial = " << std::setw(18) << stepTrial <<
" "
244 <<
" fZeroStepThreshold = " << fZeroStepThreshold <<
G4endl;
246 stepTrial, pFieldTrack);
248 if( stepTrial == 0.0 )
250 std::ostringstream message;
251 message <<
"Particle abandoned due to lack of progress in field."
253 <<
" Properties : " << pFieldTrack << G4endl
254 <<
" Attempting a zero step = " << stepTrial << G4endl
255 <<
" while attempting to progress after " << fNoZeroStep
256 <<
" trial steps. Will abandon step.";
257 G4Exception(
"G4PropagatorInField::ComputeStep()",
"GeomNav1002",
259 fParticleIsLooping=
true;
262 if( stepTrial < CurrentProposedStepLength )
263 CurrentProposedStepLength = stepTrial;
265 fLast_ProposedStepLength = CurrentProposedStepLength;
267 G4int do_loop_count = 0;
273 if( !first_substep) {
279 h_TrialStepSize = CurrentProposedStepLength - StepTaken;
292 fFull_CurveLen_of_LastAttempt = s_length_taken;
300 NewSafety, LinearStepLength,
301 InterSectionPointE );
305 if( first_substep ) {
306 currentSafety = NewSafety;
315 G4bool recalculatedEndPt=
false;
317 G4bool found_intersection = fIntersectionLocator->
318 EstimateIntersectionPoint( SubStepStartState, CurrentState,
319 InterSectionPointE, IntersectPointVelct_G,
320 recalculatedEndPt,fPreviousSafety,fPreviousSftOrigin);
321 intersects = intersects && found_intersection;
322 if( found_intersection ) {
323 End_PointAndTangent= IntersectPointVelct_G;
324 StepTaken = TruePathLength = IntersectPointVelct_G.
GetCurveLength()
328 if( recalculatedEndPt ){
329 CurrentState= IntersectPointVelct_G;
335 StepTaken += s_length_taken;
337 if (fpTrajectoryFilter) {
341 first_substep =
false;
344 if( fNoZeroStep > fActionThreshold_NoZeroSteps ) {
346 CurrentState, CurrentProposedStepLength,
347 NewSafety, do_loop_count, pPhysVol );
349 if( (fVerboseLevel > 1) && (do_loop_count > fMax_loop_count-10 )) {
350 if( do_loop_count == fMax_loop_count-9 ){
351 G4cout <<
" G4PropagatorInField::ComputeStep(): " <<
G4endl
352 <<
" Difficult track - taking many sub steps." <<
G4endl;
354 printStatus( SubStepStartState, CurrentState, CurrentProposedStepLength,
355 NewSafety, do_loop_count, pPhysVol );
361 }
while( (!intersects )
362 && (StepTaken + kCarTolerance < CurrentProposedStepLength)
363 && ( do_loop_count < fMax_loop_count ) );
365 if( do_loop_count >= fMax_loop_count )
367 fParticleIsLooping =
true;
369 if ( fVerboseLevel > 0 )
371 G4cout <<
" G4PropagateInField::ComputeStep(): " <<
G4endl
372 <<
" Killing looping particle "
374 <<
" after " << do_loop_count <<
" field substeps "
375 <<
" totaling " << StepTaken /
mm <<
" mm " ;
379 G4cout <<
" in unknown or null volume. " ;
389 End_PointAndTangent = CurrentState;
390 TruePathLength = StepTaken;
395 pFieldTrack = End_PointAndTangent;
401 - End_PointAndTangent.
GetCurveLength()) > 3.e-4 * TruePathLength )
403 std::ostringstream message;
404 message <<
"Curve length mis-match between original state "
405 <<
"and proposed endpoint of propagation." <<
G4endl
406 <<
" The curve length of the endpoint should be: "
408 <<
" and it is instead: "
410 <<
" A difference of: "
413 <<
" Original state = " << OriginalState <<
G4endl
414 <<
" Proposed state = " << End_PointAndTangent;
424 if( ( (TruePathLength < fZeroStepThreshold)
425 && ( TruePathLength+kCarTolerance < CurrentProposedStepLength )
427 || ( TruePathLength < 0.5*kCarTolerance )
436 if( fNoZeroStep > fAbandonThreshold_NoZeroSteps )
438 fParticleIsLooping =
true;
439 std::ostringstream message;
440 message <<
"Particle is stuck; it will be killed." <<
G4endl
441 <<
" Zero progress for " << fNoZeroStep <<
" attempted steps."
443 <<
" Proposed Step is " << CurrentProposedStepLength
444 <<
" but Step Taken is "<< fFull_CurveLen_of_LastAttempt <<
G4endl;
446 message <<
" in volume " << pPhysVol->
GetName() ;
448 message <<
" in unknown or null volume. " ;
454 return TruePathLength;
469 const G4int verboseLevel=fVerboseLevel;
479 if( ((stepNo == 0) && (verboseLevel <3)) || (verboseLevel >= 3) )
481 oldprec =
G4cout.precision(4);
482 G4cout << std::setw( 6) <<
" "
483 << std::setw( 25) <<
" Current Position and Direction" <<
" "
485 G4cout << std::setw( 5) <<
"Step#"
486 << std::setw(10) <<
" s " <<
" "
487 << std::setw(10) <<
"X(mm)" <<
" "
488 << std::setw(10) <<
"Y(mm)" <<
" "
489 << std::setw(10) <<
"Z(mm)" <<
" "
490 << std::setw( 7) <<
" N_x " <<
" "
491 << std::setw( 7) <<
" N_y " <<
" "
492 << std::setw( 7) <<
" N_z " <<
" " ;
493 G4cout << std::setw( 7) <<
" Delta|N|" <<
" "
494 << std::setw( 9) <<
"StepLen" <<
" "
495 << std::setw(12) <<
"StartSafety" <<
" "
496 << std::setw( 9) <<
"PhsStep" <<
" ";
498 {
G4cout << std::setw(18) <<
"NextVolume" <<
" "; }
499 G4cout.precision(oldprec);
502 if((stepNo == 0) && (verboseLevel <=3))
506 printStatus( StartFT, StartFT, -1.0, safety, -1, startVolume);
508 if( verboseLevel <= 3 )
511 {
G4cout << std::setw( 4) << stepNo <<
" "; }
513 {
G4cout << std::setw( 5) <<
"Start" ; }
514 oldprec =
G4cout.precision(8);
517 G4cout << std::setw(10) << CurrentPosition.
x() <<
" "
518 << std::setw(10) << CurrentPosition.
y() <<
" "
519 << std::setw(10) << CurrentPosition.
z() <<
" ";
521 G4cout << std::setw( 7) << CurrentUnitVelocity.
x() <<
" "
522 << std::setw( 7) << CurrentUnitVelocity.
y() <<
" "
523 << std::setw( 7) << CurrentUnitVelocity.
z() <<
" ";
527 G4cout << std::setw( 9) << step_len <<
" ";
528 G4cout << std::setw(12) << safety <<
" ";
529 if( requestStep != -1.0 )
530 {
G4cout << std::setw( 9) << requestStep <<
" "; }
532 {
G4cout << std::setw( 9) <<
"Init/NotKnown" <<
" "; }
533 if( startVolume != 0)
534 {
G4cout << std::setw(12) << startVolume->
GetName() <<
" "; }
535 G4cout.precision(oldprec);
542 G4cout <<
"Step taken was " << step_len
543 <<
" out of PhysicalStep = " << requestStep <<
G4endl;
545 G4cout <<
"Chord length = " << (CurrentPosition-StartPosition).mag()
563 G4cout <<
" " << std::setw(12) <<
" PiF: NoZeroStep "
564 <<
" " << std::setw(20) <<
" CurrentProposed len "
565 <<
" " << std::setw(18) <<
" Full_curvelen_last"
566 <<
" " << std::setw(18) <<
" last proposed len "
567 <<
" " << std::setw(18) <<
" decrease factor "
568 <<
" " << std::setw(15) <<
" step trial "
571 G4cout <<
" " << std::setw(10) << fNoZeroStep <<
" "
572 <<
" " << std::setw(20) << CurrentProposedStepLength
573 <<
" " << std::setw(18) << fFull_CurveLen_of_LastAttempt
574 <<
" " << std::setw(18) << fLast_ProposedStepLength
575 <<
" " << std::setw(18) << decreaseFactor
576 <<
" " << std::setw(15) << stepTrial
578 G4cout.precision( iprec );
591 std::vector<G4ThreeVector>*
597 if (fpTrajectoryFilter)
610 fpTrajectoryFilter =
filter;
617 fParticleIsLooping=
false;
622 0.0,0.0,0.0,0.0,0.0);
623 fFull_CurveLen_of_LastAttempt = -1;
624 fLast_ProposedStepLength = -1;
627 fPreviousSafety= 0.0;
635 currentFieldMgr = fDetectorFieldMgr;
636 if( pCurrentPhysicalVolume)
646 if( pRegionFieldMgr )
647 currentFieldMgr= pRegionFieldMgr;
653 currentFieldMgr = localFieldMgr;
656 fCurrentFieldMgr= currentFieldMgr;
661 return currentFieldMgr;
666 G4int oldval= fVerboseLevel;
667 fVerboseLevel= level;
674 G4cout <<
"Set Driver verbosity to " << fVerboseLevel - 2 <<
G4endl;
void SetEpsilonStep(G4double newEps)
G4double GetCurveLength() const
std::vector< G4ThreeVector > * GimmeTrajectoryVectorAndForgetIt() const
CLHEP::Hep3Vector G4ThreeVector
void SetVerboseLevel(G4int newLevel)
virtual void TakeIntermediatePoint(G4ThreeVector newPoint)=0
void PrintStepLengthDiagnostic(G4double currentProposedStepLength, G4double decreaseFactor, G4double stepTrial, const G4FieldTrack &aFieldTrack)
void RefreshIntersectionLocator()
G4double GetSurfaceTolerance() const
const G4ThreeVector & GetMomentumDir() const
G4double GetDeltaOneStep() const
std::vector< G4ThreeVector > * GimmeThePointsAndForgetThem()
G4PropagatorInField(G4Navigator *theNavigator, G4FieldManager *detectorFieldMgr, G4VIntersectionLocator *vLocator=0)
void SetChordFinderFor(G4ChordFinder *fCFinder)
G4Region * GetRegion() const
G4int SetVerboseLevel(G4int verbose)
G4FieldManager * GetFieldManager() const
G4ThreeVector GetPosition() const
G4GLOB_DLL std::ostream G4cout
const G4String & GetName() const
G4double GetMaximumEpsilonStep() const
void CreateNewTrajectorySegment()
G4FieldManager * FindAndSetFieldManager(G4VPhysicalVolume *pCurrentPhysVol)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
G4FieldManager * GetFieldManager() const
G4double AdvanceChordLimited(G4FieldTrack &yCurrent, G4double stepInitial, G4double epsStep_Relative, const G4ThreeVector latestSafetyOrigin, G4double lasestSafetyRadius)
G4double ComputeStep(G4FieldTrack &pFieldTrack, G4double pCurrentProposedStepLength, G4double &pNewSafety, G4VPhysicalVolume *pPhysVol=0)
G4double GetDeltaIntersection() const
G4LogicalVolume * GetLogicalVolume() const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4ChordFinder * GetChordFinder()
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
G4bool IntersectChord(const G4ThreeVector &StartPointA, const G4ThreeVector &EndPointB, G4double &NewSafety, G4double &LinearStepLength, G4ThreeVector &IntersectionPoint)
void SetTrajectoryFilter(G4VCurvedTrajectoryFilter *filter)
void ClearPropagatorState()
void SetSafetyParametersFor(G4bool UseSafety)
G4ThreeVector GetMomentum() const
G4MagInt_Driver * GetIntegrationDriver()
G4VPhysicalVolume * GetWorldVolume() const
void SetEpsilonStepFor(G4double EpsilonStep)
static G4GeometryTolerance * GetInstance()
virtual void LocateGlobalPointWithinVolume(const G4ThreeVector &position)
void SetDeltaIntersectionFor(G4double deltaIntersection)
void printStatus(const G4FieldTrack &startFT, const G4FieldTrack ¤tFT, G4double requestStep, G4double safety, G4int step, G4VPhysicalVolume *startVolume)
G4double GetMinimumEpsilonStep() const