159{
160 G4double geometryStepLength= -1.0, newSafety= -1.0;
161 fParticleIsLooping = false ;
162
163
164
165
166
167
168
169
170
171
173
174
175
180
181
182
183
184
185
186
187 G4ThreeVector OriginShift = startPosition - fPreviousSftOrigin ;
189 if( MagSqShift >=
sqr(fPreviousSafety) )
190 {
191 currentSafety = 0.0 ;
192 }
193 else
194 {
195 currentSafety = fPreviousSafety - std::sqrt(MagSqShift) ;
196 }
197
198
199
203
204 fGeometryLimitedStep = false ;
205
206
207
208
209
210
211
212
214 G4bool fieldExertsForce = false ;
215
217 G4bool fieldExists=
false;
218
220 if( fieldMgr != 0 )
221 {
222
224
225
226
227
228
230 fieldExists = (ptrField!=0) ;
231 if( fieldExists )
232 {
234
235 if( (particleCharge != 0.0)
236 || (fUseMagneticMoment && (magneticMoment != 0.0) )
237 || (gravityOn && (restMass != 0.0) )
238 )
239 {
240 fieldExertsForce = fieldExists;
241 }
242 }
243 }
244
245
246
247 if( !fieldExertsForce )
248 {
250 if( fShortStepOptimisation && (currentMinimumStep <= currentSafety) )
251 {
252
253
254 geometryStepLength = currentMinimumStep ;
255 fGeometryLimitedStep = false ;
256 }
257 else
258 {
259
260
261 linearStepLength = fLinearNavigator->
ComputeStep( startPosition,
262 startMomentumDir,
263 currentMinimumStep,
264 newSafety) ;
265
266
267 fPreviousSftOrigin = startPosition ;
268 fPreviousSafety = newSafety ;
270
271 currentSafety = newSafety ;
272
273 fGeometryLimitedStep= (linearStepLength <= currentMinimumStep);
274 if( fGeometryLimitedStep )
275 {
276
277 geometryStepLength = linearStepLength ;
278 }
279 else
280 {
281
282 geometryStepLength = currentMinimumStep ;
283 }
284 }
285 fEndPointDistance = geometryStepLength ;
286
287
288
289 fTransportEndPosition = startPosition+geometryStepLength*startMomentumDir ;
290
291
292
293 fTransportEndMomentumDir = startMomentumDir ;
296 fParticleIsLooping = false ;
297 fMomentumChanged = false ;
298 fEndGlobalTimeComputed = false ;
299 }
300 else
301 {
305
307 momentumMagnitude,
308 restMass ) ;
309
313 0.0,
315 restMass,
319 &spin ) ;
320 if( currentMinimumStep > 0 )
321 {
322
323
324 lengthAlongCurve = fFieldPropagator->
ComputeStep( aFieldTrack,
325 currentMinimumStep,
326 currentSafety,
328 fGeometryLimitedStep= lengthAlongCurve < currentMinimumStep;
329 if( fGeometryLimitedStep ) {
330 geometryStepLength = lengthAlongCurve ;
331 } else {
332 geometryStepLength = currentMinimumStep ;
333 }
334
335
336
337 fPreviousSftOrigin = startPosition ;
338 fPreviousSafety = currentSafety ;
340 }
341 else
342 {
343 geometryStepLength = lengthAlongCurve= 0.0 ;
344 fGeometryLimitedStep = false ;
345 }
346
347
348
349 fTransportEndPosition = aFieldTrack.
GetPosition() ;
350
351
352
353 fMomentumChanged = true ;
355
357
359 {
360
361
362
364 fEndGlobalTimeComputed = true;
365
366
367
368 }
369 else
370 {
371
372
373
374 fEndGlobalTimeComputed = false;
375
376
377
379 G4double endEnergy= fTransportEndKineticEnergy;
380
381 static G4int no_inexact_steps=0, no_large_ediff;
382 G4double absEdiff = std::fabs(startEnergy- endEnergy);
383 if( absEdiff > perMillion * endEnergy )
384 {
385 no_inexact_steps++;
386
387 }
388 if( fVerboseLevel > 1 )
389 {
390 if( std::fabs(startEnergy- endEnergy) > perThousand * endEnergy )
391 {
392 static G4int no_warnings= 0, warnModulo=1, moduloFactor= 10;
393 no_large_ediff ++;
394 if( (no_large_ediff% warnModulo) == 0 )
395 {
396 no_warnings++;
397 G4cout <<
"WARNING - G4Transportation::AlongStepGetPIL() "
398 <<
" Energy change in Step is above 1^-3 relative value. " <<
G4endl
399 << " Relative change in 'tracking' step = "
400 << std::setw(15) << (endEnergy-startEnergy)/startEnergy <<
G4endl
401 <<
" Starting E= " << std::setw(12) << startEnergy / MeV <<
" MeV " <<
G4endl
402 <<
" Ending E= " << std::setw(12) << endEnergy / MeV <<
" MeV " <<
G4endl;
403 G4cout <<
" Energy has been corrected -- however, review"
404 <<
" field propagation parameters for accuracy." <<
G4endl;
405 if( (fVerboseLevel > 2 ) || (no_warnings<4) || (no_large_ediff == warnModulo * moduloFactor) ){
406 G4cout <<
" These include EpsilonStepMax(/Min) in G4FieldManager "
407 <<
" which determine fractional error per step for integrated quantities. " <<
G4endl
408 << " Note also the influence of the permitted number of integration steps."
410 }
411 G4cerr <<
"ERROR - G4Transportation::AlongStepGetPIL()" <<
G4endl
412 << " Bad 'endpoint'. Energy change detected"
413 << " and corrected. "
414 << " Has occurred already "
415 << no_large_ediff <<
" times." <<
G4endl;
416 if( no_large_ediff == warnModulo * moduloFactor )
417 {
418 warnModulo *= moduloFactor;
419 }
420 }
421 }
422 }
423
424
425
426
428 }
429
430 fTransportEndSpin = aFieldTrack.
GetSpin();
432 fEndPointDistance = (fTransportEndPosition - startPosition).mag() ;
433 }
434
435
436
437
438 if( currentMinimumStep == 0.0 )
439 {
440 if( currentSafety == 0.0 ) fGeometryLimitedStep = true ;
441 }
442
443
444
445
446 if( currentSafety < fEndPointDistance )
447 {
448 if( particleCharge != 0.0 )
449 {
452 currentSafety = endSafety ;
453 fPreviousSftOrigin = fTransportEndPosition ;
454 fPreviousSafety = currentSafety ;
456
457
458
459
460 currentSafety += fEndPointDistance ;
461
462#ifdef G4DEBUG_TRANSPORT
464 G4cout <<
"***G4Transportation::AlongStepGPIL ** " <<
G4endl ;
465 G4cout <<
" Called Navigator->ComputeSafety at " << fTransportEndPosition
466 <<
" and it returned safety= " << endSafety <<
G4endl ;
467 G4cout <<
" Adding endpoint distance " << fEndPointDistance
468 <<
" to obtain pseudo-safety= " << currentSafety <<
G4endl ;
469 }else{
470 G4cout <<
"***G4Transportation::AlongStepGPIL ** " <<
G4endl ;
471 G4cout <<
" Avoiding call to ComputeSafety : " <<
G4endl;
474#endif
475 }
476 }
477
479
480 return geometryStepLength ;
481}
G4DLLIMPORT std::ostream G4cerr
const G4ThreeVector & GetMomentumDirection() const
G4double GetCharge() const
G4ParticleDefinition * GetDefinition() const
G4double GetMagneticMoment() const
G4double GetTotalMomentum() const
G4bool DoesFieldChangeEnergy() const
virtual void ConfigureForTrack(const G4Track *)
const G4Field * GetDetectorField() const
const G4ThreeVector & GetMomentumDir() const
G4double GetKineticEnergy() const
G4ThreeVector GetPosition() const
G4ThreeVector GetSpin() const
G4double GetLabTimeOfFlight() const
G4bool IsGravityActive() 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 GetPDGMass() const
G4double ComputeStep(G4FieldTrack &pFieldTrack, G4double pCurrentProposedStepLength, G4double &pNewSafety, G4VPhysicalVolume *pPhysVol=0)
G4FieldManager * FindAndSetFieldManager(G4VPhysicalVolume *pCurrentPhysVol)
G4FieldManager * GetCurrentFieldManager()
G4bool IsParticleLooping() const
void SetChargeMomentumMass(G4double charge, G4double momentum, G4double pMass)
void SetCurrentSafety(G4double val, const G4ThreeVector &pos)
G4double GetVelocity() const
G4VPhysicalVolume * GetVolume() const
const G4ThreeVector & GetPosition() const
const G4ThreeVector & GetMomentumDirection() const
const G4ThreeVector & GetPolarization() const
void ProposeTrueStepLength(G4double truePathLength)