186{
187
188
189
190
191
192
193
194
195
197
198
199
202
203
204
205
206 {
207 const G4double MagSqShift = (startPosition - fPreviousSftOrigin).mag2();
208
209 if(MagSqShift >=
sqr(fPreviousSafety))
210 currentSafety = 0.0;
211 else
212 currentSafety = fPreviousSafety - std::sqrt(MagSqShift);
213 }
214
215
216
218
222
225
226
227
228
229
230
231
232
234 (particleCharge != 0.0) || ((magneticMoment != 0.0) && fUseMagneticMoment);
235 G4bool eligibleGrav = (particleMass != 0.0) && fUseGravity;
236
237 fFieldExertedForce = false;
238
239 if(eligibleEM || eligibleGrav)
240 {
243 {
244
245 fieldMgr->ConfigureForTrack(&track);
246
247
248
249
250
251 if(
const G4Field* ptrField = fieldMgr->GetDetectorField())
252 fFieldExertedForce =
253 eligibleEM || (eligibleGrav && ptrField->IsGravityActive());
254 }
255 }
256
257 G4double geometryStepLength = currentMinimumStep;
258
259 if(currentMinimumStep == 0.0)
260 {
261 fEndPointDistance = 0.0;
262
263 fGeometryLimitedStep = (currentSafety == 0.0);
264
265 fMomentumChanged = false;
266 fParticleIsLooping = false;
267 fEndGlobalTimeComputed = false;
268 fTransportEndPosition = startPosition;
269 fTransportEndMomentumDir = startMomentumDir;
270 fTransportEndKineticEnergy = kineticEnergy;
271 fTransportEndSpin = particleSpin;
272 }
273 else if(!fFieldExertedForce)
274 {
275 fGeometryLimitedStep = false;
276 if(geometryStepLength > currentSafety || !fShortStepOptimisation)
277 {
279 startPosition, startMomentumDir, currentMinimumStep, currentSafety);
280
281 if(linearStepLength <= currentMinimumStep)
282 {
283 geometryStepLength = linearStepLength;
284 fGeometryLimitedStep = true;
285 }
286
287
288 fPreviousSftOrigin = startPosition;
289 fPreviousSafety = currentSafety;
291 }
292
293 fEndPointDistance = geometryStepLength;
294
295 fMomentumChanged = false;
296 fParticleIsLooping = false;
297 fEndGlobalTimeComputed = false;
298 fTransportEndPosition =
299 startPosition + geometryStepLength * startMomentumDir;
300 fTransportEndMomentumDir = startMomentumDir;
301 fTransportEndKineticEnergy = kineticEnergy;
302 fTransportEndSpin = particleSpin;
303 }
304 else
305 {
307 const auto particlePDGSpin = pParticleDef->
GetPDGSpin();
308 const auto particlePDGMagM = pParticleDef->GetPDGMagneticMoment();
309
311
312
313
315 G4ChargeState(particleCharge, magneticMoment, particlePDGSpin),
317
320 startMomentumDir, kineticEnergy, particleMass,
321 particleCharge, particleSpin, particlePDGMagM,
322 0.0,
323 particlePDGSpin);
324
325
326
328 aFieldTrack, currentMinimumStep, currentSafety, track.
GetVolume(),
329 kineticEnergy < fThreshold_Important_Energy);
330
331 if(lengthAlongCurve < geometryStepLength)
332 geometryStepLength = lengthAlongCurve;
333
334
335
336 fPreviousSftOrigin = startPosition;
337 fPreviousSafety = currentSafety;
339
341
342
343
344
345
348
349 fMomentumChanged = true;
351
352 fEndGlobalTimeComputed = changesEnergy;
353 fTransportEndPosition = aFieldTrack.GetPosition();
354 fTransportEndMomentumDir = aFieldTrack.GetMomentumDir();
355
356 fEndPointDistance = (fTransportEndPosition - startPosition).mag();
357
358
359
360 fTransportEndKineticEnergy =
361 changesEnergy ? aFieldTrack.GetKineticEnergy() : kineticEnergy;
362 fTransportEndSpin = aFieldTrack.GetSpin();
363
364 if(fEndGlobalTimeComputed)
365 {
366
367
368
369 fCandidateEndGlobalTime = aFieldTrack.GetLabTimeOfFlight();
370
371
372
373 }
374#if defined(G4VERBOSE) || defined(G4DEBUG_TRANSPORT)
375 else
376 {
377
378
379
380
381
382 G4double startEnergy = kineticEnergy;
383 G4double endEnergy = fTransportEndKineticEnergy;
384
386 G4double absEdiff = std::fabs(startEnergy - endEnergy);
387 if(absEdiff > perMillion * endEnergy)
388 {
389 no_inexact_steps++;
390
391 }
393 {
394 if(std::fabs(startEnergy - endEnergy) > perThousand * endEnergy)
395 {
397 moduloFactor = 10;
398 no_large_ediff++;
399 if((no_large_ediff % warnModulo) == 0)
400 {
401 no_warnings++;
402 std::ostringstream message;
403 message << "Energy change in Step is above 1^-3 relative value. "
404 <<
G4endl <<
" Relative change in 'tracking' step = "
405 << std::setw(15) << (endEnergy - startEnergy) / startEnergy
406 <<
G4endl <<
" Starting E= " << std::setw(12)
407 << startEnergy / MeV <<
" MeV " <<
G4endl
408 << " Ending E= " << std::setw(12) << endEnergy / MeV
410 << "Energy has been corrected -- however, review"
411 <<
" field propagation parameters for accuracy." <<
G4endl;
413 (no_large_ediff == warnModulo * moduloFactor))
414 {
415 message << "These include EpsilonStepMax(/Min) in G4FieldManager "
417 << "which determine fractional error per step for "
418 "integrated quantities. "
420 << "Note also the influence of the permitted number of "
421 "integration steps."
423 }
424 message << "Bad 'endpoint'. Energy change detected and corrected."
425 <<
G4endl <<
"Has occurred already " << no_large_ediff
426 << " times.";
427 G4Exception(
"G4Transportation::AlongStepGetPIL()",
"EnergyChange",
429 if(no_large_ediff == warnModulo * moduloFactor)
430 {
431 warnModulo *= moduloFactor;
432 }
433 }
434 }
435 }
436 }
437#endif
438 }
439
440
441
442
443 if(currentSafety < fEndPointDistance)
444 {
445 if(particleCharge != 0.0)
446 {
449 currentSafety = endSafety;
450 fPreviousSftOrigin = fTransportEndPosition;
451 fPreviousSafety = currentSafety;
453
454
455
456
457 currentSafety += fEndPointDistance;
458
459#ifdef G4DEBUG_TRANSPORT
461 G4cout <<
"***G4Transportation::AlongStepGPIL ** " <<
G4endl;
462 G4cout <<
" Called Navigator->ComputeSafety at " << fTransportEndPosition
463 <<
" and it returned safety= " << endSafety <<
G4endl;
464 G4cout <<
" Adding endpoint distance " << fEndPointDistance
465 <<
" to obtain pseudo-safety= " << currentSafety <<
G4endl;
466 }
467 else
468 {
469 G4cout <<
"***G4Transportation::AlongStepGPIL ** " <<
G4endl;
470 G4cout <<
" Avoiding call to ComputeSafety : " <<
G4endl;
473#endif
474 }
475 }
476
477 fFirstStepInVolume = fNewTrack || fLastStepInVolume;
478 fLastStepInVolume = false;
479 fNewTrack = false;
480
483
484 return geometryStepLength;
485}
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
G4double GetCharge() const
G4ParticleDefinition * GetDefinition() const
G4double GetKineticEnergy() const
G4double GetMagneticMoment() const
G4double GetTotalMomentum() const
const G4ThreeVector & GetPolarization() const
virtual void SetChargeMomentumMass(G4ChargeState particleCharge, G4double MomentumXc, G4double MassXc2)=0
G4bool DoesFieldChangeEnergy() const
virtual G4double ComputeSafety(const G4ThreeVector &globalpoint, const G4double pProposedMaxLength=DBL_MAX, const G4bool keepState=true)
virtual G4double ComputeStep(const G4ThreeVector &pGlobalPoint, const G4ThreeVector &pDirection, const G4double pCurrentProposedStepLength, G4double &pNewSafety)
G4double GetPDGSpin() const
G4FieldManager * FindAndSetFieldManager(G4VPhysicalVolume *pCurrentPhysVol)
G4FieldManager * GetCurrentFieldManager()
G4bool IsLastStepInVolume()
G4bool IsParticleLooping() const
G4EquationOfMotion * GetCurrentEquationOfMotion()
G4double ComputeStep(G4FieldTrack &pFieldTrack, G4double pCurrentProposedStepLength, G4double &pNewSafety, G4VPhysicalVolume *pPhysVol=nullptr, G4bool canRelaxDeltaChord=false)
void SetCurrentSafety(G4double val, const G4ThreeVector &pos)
G4VPhysicalVolume * GetVolume() const
const G4ThreeVector & GetPosition() const
const G4ThreeVector & GetMomentumDirection() const
void ProposeFirstStepInVolume(G4bool flag)
void ProposeTrueStepLength(G4double truePathLength)