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;
void SetEpsilonStep(G4double newEps)
G4double GetCurveLength() const
virtual void TakeIntermediatePoint(G4ThreeVector newPoint)=0
void PrintStepLengthDiagnostic(G4double currentProposedStepLength, G4double decreaseFactor, G4double stepTrial, const G4FieldTrack &aFieldTrack)
void RefreshIntersectionLocator()
const G4ThreeVector & GetMomentumDir() const
G4double GetDeltaOneStep() 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)
G4double AdvanceChordLimited(G4FieldTrack &yCurrent, G4double stepInitial, G4double epsStep_Relative, const G4ThreeVector latestSafetyOrigin, G4double lasestSafetyRadius)
G4LogicalVolume * GetLogicalVolume() const
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)
G4VPhysicalVolume * GetWorldVolume() const
virtual void LocateGlobalPointWithinVolume(const G4ThreeVector &position)
void printStatus(const G4FieldTrack &startFT, const G4FieldTrack ¤tFT, G4double requestStep, G4double safety, G4int step, G4VPhysicalVolume *startVolume)
G4double GetMinimumEpsilonStep() const