124{
127
128
129
130
131 const char* methodName = "G4PropagatorInField::ComputeStep";
132 if (CurrentProposedStepLength<kCarTolerance)
133 {
134 return kInfinity;
135 }
136
137
138
139 if (fpTrajectoryFilter)
140 {
142 }
143
144 fFirstStepInVolume = fNewTrack ? true : fLastStepInVolume;
145 fLastStepInVolume = false;
146 fNewTrack = false;
147
148 if( fVerboseLevel > 2 )
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
171 fParticleIsLooping = false;
172
173
174
175
176
177 if( !fSetFieldMgr )
178 {
180 }
181 fSetFieldMgr = false;
182
185
186
187
188
189 if( CurrentProposedStepLength >= fLargestAcceptableStep )
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,
199 fLargestAcceptableStep );
200 }
207
208
209
210
212
213
214
215 if( fNoZeroStep > fActionThreshold_NoZeroSteps )
216 {
218
219 stepTrial = fFull_CurveLen_of_LastAttempt;
220 if( (stepTrial <= 0.0) && (fLast_ProposedStepLength > 0.0) )
221 {
222 stepTrial = fLast_ProposedStepLength;
223 }
224
226 if( (fNoZeroStep < fSevereActionThreshold_NoZeroSteps)
227 && (stepTrial > 100.0*fZeroStepThreshold) )
228 {
229
230
231 decreaseFactor= 0.25;
232 }
233 else
234 {
235
236
237
238
239 if( stepTrial > 100.0*fZeroStepThreshold )
240 decreaseFactor = 0.35;
241 else if( stepTrial > 30.0*fZeroStepThreshold )
242 decreaseFactor= 0.5;
243 else if( stepTrial > 10.0*fZeroStepThreshold )
244 decreaseFactor= 0.75;
245 else
246 decreaseFactor= 0.9;
247 }
248 stepTrial *= decreaseFactor;
249
250#ifdef G4DEBUG_FIELD
251 if( fVerboseLevel > 2
252 || (fNoZeroStep >= fSevereActionThreshold_NoZeroSteps) )
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.";
276 fParticleIsLooping = true;
277 return 0;
278 }
279 if( stepTrial < CurrentProposedStepLength )
280 {
281 CurrentProposedStepLength = stepTrial;
282 }
283 }
284 fLast_ProposedStepLength = CurrentProposedStepLength;
285
286 G4int do_loop_count = 0;
287 do
288 {
291
292 if(!first_substep)
293 {
294 if( fVerboseLevel > 4 )
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 &&
307 fIncreaseChordDistanceThreshold > 0 &&
308 do_loop_count > fIncreaseChordDistanceThreshold &&
309 do_loop_count % fIncreaseChordDistanceThreshold == 0)
310 {
313 );
314 }
315
316
317
319 CurrentState,
320 h_TrialStepSize,
321 fEpsilonStep,
322 fPreviousSftOrigin,
323 fPreviousSafety );
324
325
326 fFull_CurveLen_of_LastAttempt = s_length_taken;
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
353 G4bool found_intersection = fIntersectionLocator->
354 EstimateIntersectionPoint( SubStepStartState, CurrentState,
355 InterSectionPointE, IntersectPointVelct_G,
356 recalculatedEndPt, fPreviousSafety,
357 fPreviousSftOrigin);
358 intersects = found_intersection;
359 if( found_intersection )
360 {
361 End_PointAndTangent= IntersectPointVelct_G;
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
377 < (endExpected*(1.0-CLHEP::perMillion));
378
381
382
383
384
385 CurrentState = IntersectPointVelct_G;
386 s_length_taken = stepAchieved;
387 if( shortEnd )
388 {
389 fParticleIsLooping = true;
390 }
391 }
392 }
393 }
394 if( !intersects )
395 {
396 StepTaken += s_length_taken;
397
398 if (fpTrajectoryFilter)
399 {
401 }
402 }
403 first_substep = false;
404
405#ifdef G4DEBUG_FIELD
406 if( fNoZeroStep > fActionThreshold_NoZeroSteps )
407 {
408 if( fNoZeroStep > fSevereActionThreshold_NoZeroSteps )
409 G4cout <<
" Above 'Severe Action' threshold -- for Zero steps. ";
410 else
411 G4cout <<
" Above 'action' threshold -- for Zero steps. ";
412 G4cout <<
" Number of zero steps = " << fNoZeroStep <<
G4endl;
414 CurrentState, CurrentProposedStepLength,
415 NewSafety, do_loop_count, pPhysVol );
416 }
417 if( (fVerboseLevel > 1) && (do_loop_count > fMax_loop_count-10 ))
418 {
419 if( do_loop_count == fMax_loop_count-9 )
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 )
434 && (!fParticleIsLooping)
435 && (StepTaken + kCarTolerance < CurrentProposedStepLength)
436 && ( do_loop_count < fMax_loop_count ) );
437
438 if( do_loop_count >= fMax_loop_count
439 && (StepTaken + kCarTolerance < CurrentProposedStepLength) )
440 {
441 fParticleIsLooping = true;
442 }
443 if ( ( fParticleIsLooping ) && (fVerboseLevel > 0) )
444 {
446 CurrentProposedStepLength, methodName,
447 CurrentState.GetMomentum(), pPhysVol );
448 }
449
450 if( !intersects )
451 {
452
453
454
455 End_PointAndTangent = CurrentState;
456 TruePathLength = StepTaken;
457
458
459
460
461 }
462 fLastStepInVolume = intersects;
463
464
465
466 pFieldTrack = End_PointAndTangent;
467
468#ifdef G4VERBOSE
469
470
472 - End_PointAndTangent.
GetCurveLength()) > 3.e-4 * TruePathLength )
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
485 << " Proposed state = " << End_PointAndTangent;
487 }
488#endif
489
490 if( TruePathLength+kCarTolerance >= CurrentProposedStepLength )
491 {
492 fNoZeroStep = 0;
493 }
494 else
495 {
496
497
498
499 if( TruePathLength < std::max( fZeroStepThreshold, 0.5*kCarTolerance ) )
500 {
501 ++fNoZeroStep;
502 }
503 else
504 {
505 fNoZeroStep = 0;
506 }
507 }
508 if( fNoZeroStep > fAbandonThreshold_NoZeroSteps )
509 {
510 fParticleIsLooping = true;
512 fFull_CurveLen_of_LastAttempt, pPhysVol );
513 fNoZeroStep = 0;
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)
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
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()
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)
void SetEpsilonStep(G4double newEps)
void CreateNewTrajectorySegment()
virtual void TakeIntermediatePoint(G4ThreeVector newPoint)=0
G4LogicalVolume * GetLogicalVolume() const
const G4String & GetName() const