154 if(CurrentProposedStepLength<kCarTolerance)
161 if (fpTrajectoryFilter)
166 fFirstStepInVolume = fNewTrack ?
true : fLastStepInVolume;
167 fLastStepInVolume=
false;
170 if( fVerboseLevel > 2 )
172 G4cout <<
"G4PropagatorInField::ComputeStep() called" <<
G4endl;
173 G4cout <<
" Starting FT: " << pFieldTrack;
174 G4cout <<
" Requested length = " << CurrentProposedStepLength <<
G4endl;
186 G4double TruePathLength = CurrentProposedStepLength;
190 G4bool first_substep =
true;
193 fParticleIsLooping =
false;
209 if( CurrentProposedStepLength >= fLargestAcceptableStep )
212 StartPointA = pFieldTrack.GetPosition();
213 VelocityUnit = pFieldTrack.GetMomentumDir();
215 G4double trialProposedStep = 1.e2 * ( 10.0 *
cm +
217 GetSolid()->DistanceToOut(StartPointA, VelocityUnit) );
218 CurrentProposedStepLength=
std::min( trialProposedStep,
219 fLargestAcceptableStep );
221 epsilon = fCurrentFieldMgr->
GetDeltaOneStep() / CurrentProposedStepLength;
225 if( epsilon < epsilonMin ) epsilon = epsilonMin;
226 if( epsilon > epsilonMax ) epsilon = epsilonMax;
239 if( fNoZeroStep > fActionThreshold_NoZeroSteps )
243 stepTrial= fFull_CurveLen_of_LastAttempt;
244 if( (stepTrial <= 0.0) && (fLast_ProposedStepLength > 0.0) )
245 stepTrial= fLast_ProposedStepLength;
248 if( (fNoZeroStep < fSevereActionThreshold_NoZeroSteps)
249 && (stepTrial > 100.0*fZeroStepThreshold) )
253 decreaseFactor= 0.25;
261 if( stepTrial > 100.0*fZeroStepThreshold )
262 decreaseFactor = 0.35;
263 else if( stepTrial > 30.0*fZeroStepThreshold )
265 else if( stepTrial > 10.0*fZeroStepThreshold )
266 decreaseFactor= 0.75;
270 stepTrial *= decreaseFactor;
273 G4cerr <<
" G4PropagatorInField::ComputeStep(): " << G4endl
274 <<
" Decreasing step - in volume " << pPhysVol;
276 G4cerr <<
" with name " << pPhysVol->GetName();
279 stepTrial, pFieldTrack);
281 if( stepTrial == 0.0 )
283 std::ostringstream message;
284 message <<
"Particle abandoned due to lack of progress in field."
286 <<
" Properties : " << pFieldTrack << G4endl
287 <<
" Attempting a zero step = " << stepTrial << G4endl
288 <<
" while attempting to progress after " << fNoZeroStep
289 <<
" trial steps. Will abandon step.";
290 G4Exception(
"G4PropagatorInField::ComputeStep()",
"GeomNav1002",
292 fParticleIsLooping=
true;
295 if( stepTrial < CurrentProposedStepLength )
296 CurrentProposedStepLength = stepTrial;
298 fLast_ProposedStepLength = CurrentProposedStepLength;
300 G4int do_loop_count = 0;
308 if( fVerboseLevel > 4 )
310 G4cout <<
" PiF: Calling Nav/Locate Global Point within-Volume "
318 h_TrialStepSize = CurrentProposedStepLength - StepTaken;
331 fFull_CurveLen_of_LastAttempt = s_length_taken;
339 NewSafety, LinearStepLength,
340 InterSectionPointE );
346 currentSafety = NewSafety;
355 G4bool recalculatedEndPt=
false;
357 G4bool found_intersection = fIntersectionLocator->
358 EstimateIntersectionPoint( SubStepStartState, CurrentState,
359 InterSectionPointE, IntersectPointVelct_G,
360 recalculatedEndPt, fPreviousSafety,
362 intersects = found_intersection;
363 if( found_intersection )
365 End_PointAndTangent= IntersectPointVelct_G;
366 StepTaken = TruePathLength = IntersectPointVelct_G.
GetCurveLength()
374 if( recalculatedEndPt )
376 G4double endAchieved = IntersectPointVelct_G.GetCurveLength();
377 G4double endExpected = CurrentState.GetCurveLength();
380 G4bool shortEnd = endAchieved
389 CurrentState= IntersectPointVelct_G;
390 s_length_taken = stepAchieved;
393 fParticleIsLooping =
true;
400 StepTaken += s_length_taken;
402 if (fpTrajectoryFilter) {
406 first_substep =
false;
409 if( fNoZeroStep > fActionThreshold_NoZeroSteps )
412 CurrentState, CurrentProposedStepLength,
413 NewSafety, do_loop_count, pPhysVol );
415 if( (fVerboseLevel > 1) && (do_loop_count > fMax_loop_count-10 ))
417 if( do_loop_count == fMax_loop_count-9 )
419 G4cout <<
" G4PropagatorInField::ComputeStep(): " << G4endl
420 <<
" Difficult track - taking many sub steps." <<
G4endl;
422 printStatus( SubStepStartState, CurrentState, CurrentProposedStepLength,
423 NewSafety, do_loop_count, pPhysVol );
429 }
while( (!intersects )
430 && (!fParticleIsLooping)
431 && (StepTaken + kCarTolerance < CurrentProposedStepLength)
432 && ( do_loop_count < fMax_loop_count ) );
434 if( do_loop_count >= fMax_loop_count )
436 fParticleIsLooping =
true;
438 if ( fParticleIsLooping && (fVerboseLevel > 0) )
448 End_PointAndTangent = CurrentState;
449 TruePathLength = StepTaken;
454 fLastStepInVolume = intersects;
458 pFieldTrack = End_PointAndTangent;
464 - End_PointAndTangent.
GetCurveLength()) > 3.e-4 * TruePathLength )
466 std::ostringstream message;
467 message <<
"Curve length mis-match between original state "
468 <<
"and proposed endpoint of propagation." << G4endl
469 <<
" The curve length of the endpoint should be: "
471 <<
" and it is instead: "
473 <<
" A difference of: "
476 <<
" Original state = " << OriginalState << G4endl
477 <<
" Proposed state = " << End_PointAndTangent;
483 if( TruePathLength+kCarTolerance >= CurrentProposedStepLength )
492 if( TruePathLength <
std::max( fZeroStepThreshold, 0.5*kCarTolerance ) )
500 if( fNoZeroStep > fAbandonThreshold_NoZeroSteps )
502 fParticleIsLooping =
true;
503 ReportStuckParticle( fNoZeroStep, CurrentProposedStepLength, fFull_CurveLen_of_LastAttempt,
508 return TruePathLength;
void SetEpsilonStep(G4double newEps)
static const G4double kInfinity
G4double GetCurveLength() const
virtual void TakeIntermediatePoint(G4ThreeVector newPoint)=0
void PrintStepLengthDiagnostic(G4double currentProposedStepLength, G4double decreaseFactor, G4double stepTrial, const G4FieldTrack &aFieldTrack)
void RefreshIntersectionLocator()
G4double GetDeltaOneStep() const
G4GLOB_DLL std::ostream G4cout
const G4String & GetName() const
void ReportLoopingParticle(G4int count, double StepTaken, G4VPhysicalVolume *pPhysVol)
G4double GetMaximumEpsilonStep() const
static constexpr double cm
void CreateNewTrajectorySegment()
static constexpr double perMillion
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
T max(const T t1, const T t2)
brief Return the largest of the two arguments
void ReportStuckParticle(G4int noZeroSteps, G4double proposedStep, G4double lastTriedStep, G4VPhysicalVolume *physVol)
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
double epsilon(double density, double temperature)
virtual void LocateGlobalPointWithinVolume(const G4ThreeVector &position)
G4GLOB_DLL std::ostream G4cerr
void printStatus(const G4FieldTrack &startFT, const G4FieldTrack ¤tFT, G4double requestStep, G4double safety, G4int step, G4VPhysicalVolume *startVolume)
G4double GetMinimumEpsilonStep() const