128{
131
132
133
134
135 const char* methodName = "G4PropagatorInField::ComputeStep";
136 if (CurrentProposedStepLength<kCarTolerance)
137 {
138 return kInfinity;
139 }
140
141
142
143 if (fpTrajectoryFilter != nullptr)
144 {
146 }
147
148 fFirstStepInVolume = fNewTrack ? true : fLastStepInVolume;
149 fLastStepInVolume = false;
150 fNewTrack = false;
151
152 if( fVerboseLevel > 2 )
153 {
155 G4cout <<
" Starting FT: " << pFieldTrack;
156 G4cout <<
" Requested length = " << CurrentProposedStepLength <<
G4endl;
158 if( pPhysVol != nullptr )
159 {
161 }
162 else
163 {
165 }
167 }
168
169
170
172 G4double TruePathLength = CurrentProposedStepLength;
176 G4bool first_substep =
true;
177
179 fParticleIsLooping = false;
180
181
182
183
184
185 if( !fSetFieldMgr )
186 {
188 }
189 fSetFieldMgr = false;
190
193
194
195
196
197 if( CurrentProposedStepLength >= fLargestAcceptableStep )
198 {
200 StartPointA = pFieldTrack.GetPosition();
201 VelocityUnit = pFieldTrack.GetMomentumDir();
202
203 G4double trialProposedStep = fMaxStepSizeMultiplier * ( fMinBigDistance +
205 GetSolid()->DistanceToOut(StartPointA, VelocityUnit) );
206 CurrentProposedStepLength = std::min( trialProposedStep,
207 fLargestAcceptableStep );
208 }
215
216
217
218
220
221
222
223 if( fNoZeroStep > fActionThreshold_NoZeroSteps )
224 {
226
227 stepTrial = fFull_CurveLen_of_LastAttempt;
228 if( (stepTrial <= 0.0) && (fLast_ProposedStepLength > 0.0) )
229 {
230 stepTrial = fLast_ProposedStepLength;
231 }
232
234 if( (fNoZeroStep < fSevereActionThreshold_NoZeroSteps)
235 && (stepTrial > 100.0*fZeroStepThreshold) )
236 {
237
238
239 decreaseFactor= 0.25;
240 }
241 else
242 {
243
244
245
246
247 if( stepTrial > 100.0*fZeroStepThreshold ) {
248 decreaseFactor = 0.35;
249 } else if( stepTrial > 30.0*fZeroStepThreshold ) {
250 decreaseFactor= 0.5;
251 } else if( stepTrial > 10.0*fZeroStepThreshold ) {
252 decreaseFactor= 0.75;
253 } else {
254 decreaseFactor= 0.9;
255 }
256 }
257 stepTrial *= decreaseFactor;
258
259#ifdef G4DEBUG_FIELD
260 if( fVerboseLevel > 2
261 || (fNoZeroStep >= fSevereActionThreshold_NoZeroSteps) )
262 {
263 G4cerr <<
" " << methodName
264 << " Decreasing step after " << fNoZeroStep << " zero steps "
265 << " - in volume " << pPhysVol;
266 if( pPhysVol )
268 else
269 G4cerr <<
" i.e. *unknown* volume.";
272 stepTrial, pFieldTrack);
273 }
274#endif
275 if( stepTrial == 0.0 )
276 {
277 std::ostringstream message;
278 message << "Particle abandoned due to lack of progress in field."
280 <<
" Properties : " << pFieldTrack <<
G4endl
281 <<
" Attempting a zero step = " << stepTrial <<
G4endl
282 << " while attempting to progress after " << fNoZeroStep
283 << " trial steps. Will abandon step.";
285 fParticleIsLooping = true;
286 return 0;
287 }
288 if( stepTrial < CurrentProposedStepLength )
289 {
290 CurrentProposedStepLength = stepTrial;
291 }
292 }
293 fLast_ProposedStepLength = CurrentProposedStepLength;
294
295 G4int do_loop_count = 0;
296 do
297 {
300
301 if(!first_substep)
302 {
303 if( fVerboseLevel > 4 )
304 {
305 G4cout <<
" PiF: Calling Nav/Locate Global Point within-Volume "
307 }
309 }
310
311
312
313 h_TrialStepSize = CurrentProposedStepLength - StepTaken;
314
315 if (canRelaxDeltaChord &&
316 fIncreaseChordDistanceThreshold > 0 &&
317 do_loop_count > fIncreaseChordDistanceThreshold &&
318 do_loop_count % fIncreaseChordDistanceThreshold == 0)
319 {
322 );
323 }
324
325
326
328 CurrentState,
329 h_TrialStepSize,
330 fEpsilonStep,
331 fPreviousSftOrigin,
332 fPreviousSafety );
333
334
335 fFull_CurveLen_of_LastAttempt = s_length_taken;
336
340
341
342
344 NewSafety, LinearStepLength,
345 InterSectionPointE );
346
347
348
349 if( first_substep )
350 {
351 currentSafety = NewSafety;
352 }
353
354 if( intersects )
355 {
357
358
359
360 G4bool recalculatedEndPt =
false;
361
362 G4bool found_intersection = fIntersectionLocator->
363 EstimateIntersectionPoint( SubStepStartState, CurrentState,
364 InterSectionPointE, IntersectPointVelct_G,
365 recalculatedEndPt, fPreviousSafety,
366 fPreviousSftOrigin);
367 intersects = found_intersection;
368 if( found_intersection )
369 {
370 End_PointAndTangent= IntersectPointVelct_G;
371 StepTaken = TruePathLength = IntersectPointVelct_G.
GetCurveLength()
373 }
374 else
375 {
376
377
378
379 if( recalculatedEndPt )
380 {
381 G4double endAchieved = IntersectPointVelct_G.GetCurveLength();
382 G4double endExpected = CurrentState.GetCurveLength();
383
384
385 G4bool shortEnd = endAchieved
386 < (endExpected*(1.0-CLHEP::perMillion));
387
390
391
392
393
394 CurrentState = IntersectPointVelct_G;
395 s_length_taken = stepAchieved;
396 if( shortEnd )
397 {
398 fParticleIsLooping = true;
399 }
400 }
401 }
402 }
403 if( !intersects )
404 {
405 StepTaken += s_length_taken;
406
407 if (fpTrajectoryFilter != nullptr)
408 {
410 }
411 }
412 first_substep = false;
413
414#ifdef G4DEBUG_FIELD
415 if( fNoZeroStep > fActionThreshold_NoZeroSteps )
416 {
417 if( fNoZeroStep > fSevereActionThreshold_NoZeroSteps )
418 G4cout <<
" Above 'Severe Action' threshold -- for Zero steps. ";
419 else
420 G4cout <<
" Above 'action' threshold -- for Zero steps. ";
421 G4cout <<
" Number of zero steps = " << fNoZeroStep <<
G4endl;
423 CurrentState, CurrentProposedStepLength,
424 NewSafety, do_loop_count, pPhysVol );
425 }
426 if( (fVerboseLevel > 1) && (do_loop_count > fMax_loop_count-10 ))
427 {
428 if( do_loop_count == fMax_loop_count-9 )
429 {
430 G4cout <<
" G4PropagatorInField::ComputeStep(): " <<
G4endl
431 <<
" Difficult track - taking many sub steps." <<
G4endl;
432 printStatus( SubStepStartState, SubStepStartState, CurrentProposedStepLength,
433 NewSafety, 0, pPhysVol );
434 }
435 printStatus( SubStepStartState, CurrentState, CurrentProposedStepLength,
436 NewSafety, do_loop_count, pPhysVol );
437 }
438#endif
439
440 ++do_loop_count;
441
442 } while( (!intersects )
443 && (!fParticleIsLooping)
444 && (StepTaken + kCarTolerance < CurrentProposedStepLength)
445 && ( do_loop_count < fMax_loop_count ) );
446
447 if( do_loop_count >= fMax_loop_count
448 && (StepTaken + kCarTolerance < CurrentProposedStepLength) )
449 {
450 fParticleIsLooping = true;
451 }
452 if ( ( fParticleIsLooping ) && (fVerboseLevel > 0) )
453 {
455 CurrentProposedStepLength, methodName,
456 CurrentState.GetMomentum(), pPhysVol );
457 }
458
459 if( !intersects )
460 {
461
462
463
464 End_PointAndTangent = CurrentState;
465 TruePathLength = StepTaken;
466
467
468
469
470 }
471 fLastStepInVolume = intersects;
472
473
474
475 pFieldTrack = End_PointAndTangent;
476
477#ifdef G4VERBOSE
478
479
481 - End_PointAndTangent.
GetCurveLength()) > 3.e-4 * TruePathLength )
482 {
483 std::ostringstream message;
484 message << "Curve length mis-match between original state "
485 <<
"and proposed endpoint of propagation." <<
G4endl
486 << " The curve length of the endpoint should be: "
488 << " and it is instead: "
490 << " A difference of: "
493 <<
" Original state = " << OriginalState <<
G4endl
494 << " Proposed state = " << End_PointAndTangent;
496 }
497#endif
498
499 if( TruePathLength+kCarTolerance >= CurrentProposedStepLength )
500 {
501 fNoZeroStep = 0;
502 }
503 else
504 {
505
506
507
508 if( TruePathLength < std::max( fZeroStepThreshold, 0.5*kCarTolerance ) )
509 {
510 ++fNoZeroStep;
511 }
512 else
513 {
514 fNoZeroStep = 0;
515 }
516 }
517 if( fNoZeroStep > fAbandonThreshold_NoZeroSteps )
518 {
519 fParticleIsLooping = true;
521 fFull_CurveLen_of_LastAttempt, pPhysVol );
522 fNoZeroStep = 0;
523 }
524
526 return TruePathLength;
527}
G4double epsilon(G4double density, G4double temperature)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
G4GLOB_DLL std::ostream G4cerr
G4double GetDeltaChord() const
void OnComputeStep(const G4FieldTrack *track)
void SetDeltaChord(G4double newval)
G4double AdvanceChordLimited(G4FieldTrack &yCurrent, G4double stepInitial, G4double epsStep_Relative, const G4ThreeVector &latestSafetyOrigin, G4double lasestSafetyRadius)
G4double GetMinimumEpsilonStep() const
G4double GetDeltaOneStep() const
G4double GetCurveLength() const
virtual void LocateGlobalPointWithinVolume(const G4ThreeVector &position)
G4VPhysicalVolume * GetWorldVolume() const
void printStatus(const G4FieldTrack &startFT, const G4FieldTrack ¤tFT, G4double requestStep, G4double safety, G4int step, G4VPhysicalVolume *startVolume)
G4FieldManager * FindAndSetFieldManager(G4VPhysicalVolume *pCurrentPhysVol)
void PrintStepLengthDiagnostic(G4double currentProposedStepLength, G4double decreaseFactor, G4double stepTrial, const G4FieldTrack &aFieldTrack)
void ReportStuckParticle(G4int noZeroSteps, G4double proposedStep, G4double lastTriedStep, G4VPhysicalVolume *physVol)
G4ChordFinder * GetChordFinder()
G4bool IntersectChord(const G4ThreeVector &StartPointA, const G4ThreeVector &EndPointB, G4double &NewSafety, G4double &LinearStepLength, G4ThreeVector &IntersectionPoint)
void ReportLoopingParticle(G4int count, G4double StepTaken, G4double stepRequest, const char *methodName, const G4ThreeVector &momentumVec, G4VPhysicalVolume *physVol)
void SetEpsilonStep(G4double newEps)
void CreateNewTrajectorySegment()
virtual void TakeIntermediatePoint(G4ThreeVector newPoint)=0
const G4String & GetName() const