124{
127
128
129
130
131 const char* methodName = "G4PropagatorInField::ComputeStep";
133 {
135 }
136
137
138
140 {
142 }
143
147
149 {
151 G4cout <<
" Starting FT: " << pFieldTrack;
152 G4cout <<
" Requested length = " << CurrentProposedStepLength <<
G4endl;
154 if( pPhysVol )
156 else
159 }
160
161
162
164 G4double TruePathLength = CurrentProposedStepLength;
168 G4bool first_substep =
true;
169
172
173
174
175
176
178 {
180 }
182
185
186
187
188
190 {
192 StartPointA = pFieldTrack.GetPosition();
193 VelocityUnit = pFieldTrack.GetMomentumDir();
194
195 G4double trialProposedStep = 1.e2 * ( 10.0 *
cm +
197 GetSolid()->DistanceToOut(StartPointA, VelocityUnit) );
198 CurrentProposedStepLength =
std::min( trialProposedStep,
200 }
207
208
209
210
212
213
214
216 {
218
221 {
223 }
224
228 {
229
230
231 decreaseFactor= 0.25;
232 }
233 else
234 {
235
236
237
238
240 decreaseFactor = 0.35;
242 decreaseFactor= 0.5;
244 decreaseFactor= 0.75;
245 else
246 decreaseFactor= 0.9;
247 }
248 stepTrial *= decreaseFactor;
249
250#ifdef G4DEBUG_FIELD
253 {
254 G4cerr <<
" " << methodName
255 <<
" Decreasing step after " <<
fNoZeroStep <<
" zero steps "
256 << " - in volume " << pPhysVol;
257 if( pPhysVol )
258 G4cerr <<
" with name " << pPhysVol->GetName();
259 else
260 G4cerr <<
" i.e. *unknown* volume.";
263 stepTrial, pFieldTrack);
264 }
265#endif
266 if( stepTrial == 0.0 )
267 {
268 std::ostringstream message;
269 message << "Particle abandoned due to lack of progress in field."
271 <<
" Properties : " << pFieldTrack <<
G4endl
272 <<
" Attempting a zero step = " << stepTrial <<
G4endl
273 <<
" while attempting to progress after " <<
fNoZeroStep
274 << " trial steps. Will abandon step.";
277 return 0;
278 }
279 if( stepTrial < CurrentProposedStepLength )
280 {
281 CurrentProposedStepLength = stepTrial;
282 }
283 }
285
286 G4int do_loop_count = 0;
287 do
288 {
291
292 if(!first_substep)
293 {
295 {
296 G4cout <<
" PiF: Calling Nav/Locate Global Point within-Volume "
298 }
300 }
301
302
303
304 h_TrialStepSize = CurrentProposedStepLength - StepTaken;
305
306 if (canRelaxDeltaChord &&
310 {
313 );
314 }
315
316
317
319 CurrentState,
320 h_TrialStepSize,
324
325
327
331
332
333
335 NewSafety, LinearStepLength,
336 InterSectionPointE );
337
338
339
340 if( first_substep )
341 {
342 currentSafety = NewSafety;
343 }
344
345 if( intersects )
346 {
348
349
350
351 G4bool recalculatedEndPt =
false;
352
354 EstimateIntersectionPoint( SubStepStartState, CurrentState,
355 InterSectionPointE, IntersectPointVelct_G,
358 intersects = found_intersection;
359 if( found_intersection )
360 {
362 StepTaken = TruePathLength = IntersectPointVelct_G.
GetCurveLength()
364 }
365 else
366 {
367
368
369
370 if( recalculatedEndPt )
371 {
372 G4double endAchieved = IntersectPointVelct_G.GetCurveLength();
373 G4double endExpected = CurrentState.GetCurveLength();
374
375
376 G4bool shortEnd = endAchieved
378
381
382
383
384
385 CurrentState = IntersectPointVelct_G;
386 s_length_taken = stepAchieved;
387 if( shortEnd )
388 {
390 }
391 }
392 }
393 }
394 if( !intersects )
395 {
396 StepTaken += s_length_taken;
397
399 {
401 }
402 }
403 first_substep = false;
404
405#ifdef G4DEBUG_FIELD
407 {
409 G4cout <<
" Above 'Severe Action' threshold -- for Zero steps. ";
410 else
411 G4cout <<
" Above 'action' threshold -- for Zero steps. ";
414 CurrentState, CurrentProposedStepLength,
415 NewSafety, do_loop_count, pPhysVol );
416 }
418 {
420 {
421 G4cout <<
" G4PropagatorInField::ComputeStep(): " <<
G4endl
422 <<
" Difficult track - taking many sub steps." <<
G4endl;
423 printStatus( SubStepStartState, SubStepStartState, CurrentProposedStepLength,
424 NewSafety, 0, pPhysVol );
425 }
426 printStatus( SubStepStartState, CurrentState, CurrentProposedStepLength,
427 NewSafety, do_loop_count, pPhysVol );
428 }
429#endif
430
431 ++do_loop_count;
432
433 } while( (!intersects )
437
440 {
442 }
444 {
446 CurrentProposedStepLength, methodName,
447 CurrentState.GetMomentum(), pPhysVol );
448 }
449
450 if( !intersects )
451 {
452
453
454
456 TruePathLength = StepTaken;
457
458
459
460
461 }
463
464
465
467
468#ifdef G4VERBOSE
469
470
473 {
474 std::ostringstream message;
475 message << "Curve length mis-match between original state "
476 <<
"and proposed endpoint of propagation." <<
G4endl
477 << " The curve length of the endpoint should be: "
479 << " and it is instead: "
481 << " A difference of: "
484 <<
" Original state = " << OriginalState <<
G4endl
487 }
488#endif
489
490 if( TruePathLength+
kCarTolerance >= CurrentProposedStepLength )
491 {
493 }
494 else
495 {
496
497
498
500 {
502 }
503 else
504 {
506 }
507 }
509 {
514 }
515
517 return TruePathLength;
518}
G4double epsilon(G4double density, G4double temperature)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
static constexpr double cm
G4GLOB_DLL std::ostream G4cerr
G4double GetDeltaChord() const
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
G4VCurvedTrajectoryFilter * fpTrajectoryFilter
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)
G4int fAbandonThreshold_NoZeroSteps
void ReportStuckParticle(G4int noZeroSteps, G4double proposedStep, G4double lastTriedStep, G4VPhysicalVolume *physVol)
G4bool fFirstStepInVolume
G4int fSevereActionThreshold_NoZeroSteps
G4ChordFinder * GetChordFinder()
void ReportLoopingParticle(G4int count, G4double StepTaken, G4double stepRequest, const char *methodName, G4ThreeVector momentumVec, G4VPhysicalVolume *physVol)
G4bool IntersectChord(const G4ThreeVector &StartPointA, const G4ThreeVector &EndPointB, G4double &NewSafety, G4double &LinearStepLength, G4ThreeVector &IntersectionPoint)
G4int fIncreaseChordDistanceThreshold
G4int fActionThreshold_NoZeroSteps
void SetEpsilonStep(G4double newEps)
void CreateNewTrajectorySegment()
virtual void TakeIntermediatePoint(G4ThreeVector newPoint)=0
G4LogicalVolume * GetLogicalVolume() const
const G4String & GetName() const
static const G4double kInfinity
static constexpr double perMillion
T min(const T t1, const T t2)
brief Return the smallest of the two arguments