Geant4  10.03.p01
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
G4Transportation.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // $Id: G4Transportation.cc 2011/06/10 16:19:46 japost Exp japost $
28 //
29 // ------------------------------------------------------------
30 // GEANT 4 include file implementation
31 //
32 // ------------------------------------------------------------
33 //
34 // This class is a process responsible for the transportation of
35 // a particle, ie the geometrical propagation that encounters the
36 // geometrical sub-volumes of the detectors.
37 //
38 // It is also tasked with the key role of proposing the "isotropic safety",
39 // which will be used to update the post-step point's safety.
40 //
41 // =======================================================================
42 // Modified:
43 // 10 Jan 2015, M.Kelsey: Use G4DynamicParticle mass, NOT PDGMass
44 // 28 Oct 2011, P.Gumpl./J.Ap: Detect gravity field, use magnetic moment
45 // 20 Nov 2008, J.Apostolakis: Push safety to helper - after ComputeSafety
46 // 9 Nov 2007, J.Apostolakis: Flag for short steps, push safety to helper
47 // 19 Jan 2006, P.MoraDeFreitas: Fix for suspended tracks (StartTracking)
48 // 11 Aug 2004, M.Asai: Add G4VSensitiveDetector* for updating stepPoint.
49 // 21 June 2003, J.Apostolakis: Calling field manager with
50 // track, to enable it to configure its accuracy
51 // 13 May 2003, J.Apostolakis: Zero field areas now taken into
52 // account correclty in all cases (thanks to W Pokorski).
53 // 29 June 2001, J.Apostolakis, D.Cote-Ahern, P.Gumplinger:
54 // correction for spin tracking
55 // 20 Febr 2001, J.Apostolakis: update for new FieldTrack
56 // 22 Sept 2000, V.Grichine: update of Kinetic Energy
57 // Created: 19 March 1997, J. Apostolakis
58 // =======================================================================
59 
60 #include "G4Transportation.hh"
62 
63 #include "G4PhysicalConstants.hh"
64 #include "G4SystemOfUnits.hh"
65 #include "G4ProductionCutsTable.hh"
66 #include "G4ParticleTable.hh"
67 
68 #include "G4ChargeState.hh"
69 #include "G4EquationOfMotion.hh"
70 
71 #include "G4FieldManagerStore.hh"
72 
74 
75 G4bool G4Transportation::fUseMagneticMoment=false;
76 
77 // #define G4DEBUG_TRANSPORT 1
78 
80 //
81 // Constructor
82 
84  : G4VProcess( G4String("Transportation"), fTransportation ),
85  fTransportEndPosition( 0.0, 0.0, 0.0 ),
86  fTransportEndMomentumDir( 0.0, 0.0, 0.0 ),
87  fTransportEndKineticEnergy( 0.0 ),
88  fTransportEndSpin( 0.0, 0.0, 0.0 ),
89  fMomentumChanged(true),
90  fEndGlobalTimeComputed(false),
91  fCandidateEndGlobalTime(0.0),
92  fParticleIsLooping( false ),
93  fNewTrack( true ),
94  fFirstStepInVolume( true ),
95  fLastStepInVolume( false ),
96  fGeometryLimitedStep(true),
97  fFieldExertedForce( false ),
98  fPreviousSftOrigin( 0.,0.,0. ),
99  fPreviousSafety( 0.0 ),
100  // fParticleChange(),
101  fEndPointDistance( -1.0 ),
102  fThreshold_Warning_Energy( 100 * MeV ),
103  fThreshold_Important_Energy( 250 * MeV ),
104  fThresholdTrials( 10 ),
105  fNoLooperTrials( 0 ),
106  fSumEnergyKilled( 0.0 ), fMaxEnergyKilled( 0.0 ),
107  fShortStepOptimisation( false ), // Old default: true (=fast short steps)
108  fVerboseLevel( verbosity )
109 {
110  // set Process Sub Type
111  SetProcessSubType(static_cast<G4int>(TRANSPORTATION));
112  pParticleChange= &fParticleChange; // Required to conform to G4VProcess
113 
114  G4TransportationManager* transportMgr ;
115 
117 
118  fLinearNavigator = transportMgr->GetNavigatorForTracking() ;
119 
120  fFieldPropagator = transportMgr->GetPropagatorInField() ;
121 
122  fpSafetyHelper = transportMgr->GetSafetyHelper(); // New
123 
124  // Cannot determine whether a field exists here, as it would
125  // depend on the relative order of creating the detector's
126  // field and this process. That order is not guaranted.
127  // Instead later the method DoesGlobalFieldExist() is called
128 
129  static G4ThreadLocal G4TouchableHandle* pNullTouchableHandle = 0;
130  if ( !pNullTouchableHandle) { pNullTouchableHandle = new G4TouchableHandle; }
131  fCurrentTouchableHandle = *pNullTouchableHandle;
132  // Points to (G4VTouchable*) 0
133 
134 #ifdef G4VERBOSE
135  if( fVerboseLevel > 0)
136  {
137  G4cout << " G4Transportation constructor> set fShortStepOptimisation to ";
138  if ( fShortStepOptimisation ) G4cout << "true" << G4endl;
139  else G4cout << "false" << G4endl;
140  }
141 #endif
142 }
143 
145 
147 {
148  if( (fVerboseLevel > 0) && (fSumEnergyKilled > 0.0 ) )
149  {
150  G4cout << " G4Transportation: Statistics for looping particles " << G4endl;
151  G4cout << " Sum of energy of loopers killed: " << fSumEnergyKilled << G4endl;
152  G4cout << " Max energy of loopers killed: " << fMaxEnergyKilled << G4endl;
153  }
154 }
155 
157 //
158 // Responsibilities:
159 // Find whether the geometry limits the Step, and to what length
160 // Calculate the new value of the safety and return it.
161 // Store the final time, position and momentum.
162 
165  G4double, // previousStepSize
166  G4double currentMinimumStep,
167  G4double& currentSafety,
168  G4GPILSelection* selection )
169 {
170  G4double geometryStepLength= -1.0, newSafety= -1.0;
171  fParticleIsLooping = false ;
172 
173  // Initial actions moved to StartTrack()
174  // --------------------------------------
175  // Note: in case another process changes touchable handle
176  // it will be necessary to add here (for all steps)
177  // fCurrentTouchableHandle = aTrack->GetTouchableHandle();
178 
179  // GPILSelection is set to defaule value of CandidateForSelection
180  // It is a return value
181  //
182  *selection = CandidateForSelection ;
183 
184  fFirstStepInVolume= fNewTrack || fLastStepInVolume;
185  // G4cout << " Transport::AlongStep GPIL: 1st-step= " << fFirstStepInVolume << " newTrack= " << fNewTrack << " fLastStep-in-Vol= " << fLastStepInVolume << G4endl;
186  fLastStepInVolume= false;
187  fNewTrack = false;
188 
189  fParticleChange.ProposeFirstStepInVolume(fFirstStepInVolume);
190 
191  // Get initial Energy/Momentum of the track
192  //
193  const G4DynamicParticle* pParticle = track.GetDynamicParticle() ;
194  const G4ParticleDefinition* pParticleDef = pParticle->GetDefinition() ;
195  G4ThreeVector startMomentumDir = pParticle->GetMomentumDirection() ;
196  G4ThreeVector startPosition = track.GetPosition() ;
197 
198  // G4double theTime = track.GetGlobalTime() ;
199 
200  // The Step Point safety can be limited by other geometries and/or the
201  // assumptions of any process - it's not always the geometrical safety.
202  // We calculate the starting point's isotropic safety here.
203  //
204  G4ThreeVector OriginShift = startPosition - fPreviousSftOrigin ;
205  G4double MagSqShift = OriginShift.mag2() ;
206  if( MagSqShift >= sqr(fPreviousSafety) )
207  {
208  currentSafety = 0.0 ;
209  }
210  else
211  {
212  currentSafety = fPreviousSafety - std::sqrt(MagSqShift) ;
213  }
214 
215  // Is the particle charged or has it a magnetic moment?
216  //
217  G4double particleCharge = pParticle->GetCharge() ;
218  G4double magneticMoment = pParticle->GetMagneticMoment() ;
219  G4double restMass = pParticle->GetMass() ;
220 
221  fGeometryLimitedStep = false ;
222  // fEndGlobalTimeComputed = false ;
223 
224  // There is no need to locate the current volume. It is Done elsewhere:
225  // On track construction
226  // By the tracking, after all AlongStepDoIts, in "Relocation"
227 
228  // Check if the particle has a force, EM or gravitational, exerted on it
229  //
230  G4FieldManager* fieldMgr=0;
231  G4bool fieldExertsForce = false ;
232 
233  G4bool gravityOn = false;
234  G4bool fieldExists= false; // Field is not 0 (null pointer)
235 
236  fieldMgr = fFieldPropagator->FindAndSetFieldManager( track.GetVolume() );
237  if( fieldMgr != 0 )
238  {
239  // Message the field Manager, to configure it for this track
240  fieldMgr->ConfigureForTrack( &track );
241  // Is here to allow a transition from no-field pointer
242  // to finite field (non-zero pointer).
243 
244  // If the field manager has no field ptr, the field is zero
245  // by definition ( = there is no field ! )
246  const G4Field* ptrField= fieldMgr->GetDetectorField();
247  fieldExists = (ptrField!=0) ;
248  if( fieldExists )
249  {
250  gravityOn= ptrField->IsGravityActive();
251 
252  if( (particleCharge != 0.0)
253  || (fUseMagneticMoment && (magneticMoment != 0.0) )
254  || (gravityOn && (restMass != 0.0) )
255  )
256  {
257  fieldExertsForce = fieldExists;
258  }
259  }
260  }
261  // G4cout << " G4Transport: field exerts force= " << fieldExertsForce
262  // << " fieldMgr= " << fieldMgr << G4endl;
263  fFieldExertedForce = fieldExertsForce;
264 
265  if( !fieldExertsForce )
266  {
267  G4double linearStepLength ;
268  if( fShortStepOptimisation && (currentMinimumStep <= currentSafety) )
269  {
270  // The Step is guaranteed to be taken
271  //
272  geometryStepLength = currentMinimumStep ;
273  fGeometryLimitedStep = false ;
274  }
275  else
276  {
277  // Find whether the straight path intersects a volume
278  //
279  linearStepLength = fLinearNavigator->ComputeStep( startPosition,
280  startMomentumDir,
281  currentMinimumStep,
282  newSafety) ;
283  // Remember last safety origin & value.
284  //
285  fPreviousSftOrigin = startPosition ;
286  fPreviousSafety = newSafety ;
287  fpSafetyHelper->SetCurrentSafety( newSafety, startPosition);
288 
289  currentSafety = newSafety ;
290 
291  fGeometryLimitedStep= (linearStepLength <= currentMinimumStep);
292  if( fGeometryLimitedStep )
293  {
294  // The geometry limits the Step size (an intersection was found.)
295  geometryStepLength = linearStepLength ;
296  }
297  else
298  {
299  // The full Step is taken.
300  geometryStepLength = currentMinimumStep ;
301  }
302  }
303  fEndPointDistance = geometryStepLength ;
304 
305  // Calculate final position
306  //
307  fTransportEndPosition = startPosition+geometryStepLength*startMomentumDir ;
308 
309  // Momentum direction, energy and polarisation are unchanged by transport
310  //
311  fTransportEndMomentumDir = startMomentumDir ;
312  fTransportEndKineticEnergy = track.GetKineticEnergy() ;
313  fTransportEndSpin = track.GetPolarization();
314  fParticleIsLooping = false ;
315  fMomentumChanged = false ;
316  fEndGlobalTimeComputed = false ;
317  }
318  else // A field exerts force
319  {
320  G4double momentumMagnitude = pParticle->GetTotalMomentum() ;
321  G4ThreeVector EndUnitMomentum ;
322  G4double lengthAlongCurve ;
323 
324  G4ChargeState chargeState(particleCharge, // The charge can change (dynamic)
325  magneticMoment,
326  pParticleDef->GetPDGSpin() );
327  // For insurance, could set it again
328  // chargeState.SetPDGSpin(pParticleDef->GetPDGSpin() ); // Provisionally in same object
329 
330  G4EquationOfMotion* equationOfMotion =
331  (fFieldPropagator->GetChordFinder()->GetIntegrationDriver()->GetStepper())
332  ->GetEquationOfMotion();
333 
334 // equationOfMotion->SetChargeMomentumMass( particleCharge,
335  equationOfMotion->SetChargeMomentumMass( chargeState,
336  momentumMagnitude,
337  restMass);
338 
339  G4FieldTrack aFieldTrack = G4FieldTrack( startPosition,
340  track.GetGlobalTime(), // Lab.
341  // track.GetProperTime(), // Particle rest frame
342  track.GetMomentumDirection(),
343  track.GetKineticEnergy(),
344  restMass,
345  particleCharge,
346  track.GetPolarization(),
347  pParticleDef->GetPDGMagneticMoment(),
348  0.0, // Length along track
349  pParticleDef->GetPDGSpin()
350  ) ;
351 
352  if( currentMinimumStep > 0 )
353  {
354  // Do the Transport in the field (non recti-linear)
355  //
356  lengthAlongCurve = fFieldPropagator->ComputeStep( aFieldTrack,
357  currentMinimumStep,
358  currentSafety,
359  track.GetVolume() ) ;
360 
361  fGeometryLimitedStep= fFieldPropagator->IsLastStepInVolume();
362  // It is possible that step was reduced in PropagatorInField due to previous zero steps
363  // To cope with case that reduced step is taken in full, we must rely on PiF to obtain this
364  // value.
365 
366  geometryStepLength = std::min( lengthAlongCurve, currentMinimumStep );
367 
368  // Remember last safety origin & value.
369  //
370  fPreviousSftOrigin = startPosition ;
371  fPreviousSafety = currentSafety ;
372  fpSafetyHelper->SetCurrentSafety( currentSafety, startPosition);
373  }
374  else
375  {
376  geometryStepLength = lengthAlongCurve= 0.0 ;
377  fGeometryLimitedStep = false ;
378  }
379 
380  // Get the End-Position and End-Momentum (Dir-ection)
381  //
382  fTransportEndPosition = aFieldTrack.GetPosition() ;
383 
384  // Momentum: Magnitude and direction can be changed too now ...
385  //
386  fMomentumChanged = true ;
387  fTransportEndMomentumDir = aFieldTrack.GetMomentumDir() ;
388 
389  fTransportEndKineticEnergy = aFieldTrack.GetKineticEnergy() ;
390 
391  if( fFieldPropagator->GetCurrentFieldManager()->DoesFieldChangeEnergy() )
392  {
393  // If the field can change energy, then the time must be integrated
394  // - so this should have been updated
395  //
396  fCandidateEndGlobalTime = aFieldTrack.GetLabTimeOfFlight();
397  fEndGlobalTimeComputed = true;
398 
399  // was ( fCandidateEndGlobalTime != track.GetGlobalTime() );
400  // a cleaner way is to have FieldTrack knowing whether time is updated.
401  }
402  else
403  {
404  // The energy should be unchanged by field transport,
405  // - so the time changed will be calculated elsewhere
406  //
407  fEndGlobalTimeComputed = false;
408 
409  // Check that the integration preserved the energy
410  // - and if not correct this!
411  G4double startEnergy= track.GetKineticEnergy();
412  G4double endEnergy= fTransportEndKineticEnergy;
413 
414  static G4ThreadLocal G4int no_inexact_steps=0, no_large_ediff;
415  G4double absEdiff = std::fabs(startEnergy- endEnergy);
416  if( absEdiff > perMillion * endEnergy )
417  {
418  no_inexact_steps++;
419  // Possible statistics keeping here ...
420  }
421  if( fVerboseLevel > 1 )
422  {
423  if( std::fabs(startEnergy- endEnergy) > perThousand * endEnergy )
424  {
425  static G4ThreadLocal G4int no_warnings= 0, warnModulo=1,
426  moduloFactor= 10;
427  no_large_ediff ++;
428  if( (no_large_ediff% warnModulo) == 0 )
429  {
430  no_warnings++;
431  G4cout << "WARNING - G4Transportation::AlongStepGetPIL() "
432  << " Energy change in Step is above 1^-3 relative value. " << G4endl
433  << " Relative change in 'tracking' step = "
434  << std::setw(15) << (endEnergy-startEnergy)/startEnergy << G4endl
435  << " Starting E= " << std::setw(12) << startEnergy / MeV << " MeV " << G4endl
436  << " Ending E= " << std::setw(12) << endEnergy / MeV << " MeV " << G4endl;
437  G4cout << " Energy has been corrected -- however, review"
438  << " field propagation parameters for accuracy." << G4endl;
439  if( (fVerboseLevel > 2 ) || (no_warnings<4) || (no_large_ediff == warnModulo * moduloFactor) )
440  {
441  G4cout << " These include EpsilonStepMax(/Min) in G4FieldManager "
442  << " which determine fractional error per step for integrated quantities. " << G4endl
443  << " Note also the influence of the permitted number of integration steps."
444  << G4endl;
445  }
446  G4cerr << "ERROR - G4Transportation::AlongStepGetPIL()" << G4endl
447  << " Bad 'endpoint'. Energy change detected"
448  << " and corrected. "
449  << " Has occurred already "
450  << no_large_ediff << " times." << G4endl;
451  if( no_large_ediff == warnModulo * moduloFactor )
452  {
453  warnModulo *= moduloFactor;
454  }
455  }
456  }
457  } // end of if (fVerboseLevel)
458 
459  // Correct the energy for fields that conserve it
460  // This - hides the integration error
461  // - but gives a better physical answer
462  fTransportEndKineticEnergy= track.GetKineticEnergy();
463  }
464 
465  fTransportEndSpin = aFieldTrack.GetSpin();
466  fParticleIsLooping = fFieldPropagator->IsParticleLooping() ;
467  fEndPointDistance = (fTransportEndPosition - startPosition).mag() ;
468  }
469 
470  // If we are asked to go a step length of 0, and we are on a boundary
471  // then a boundary will also limit the step -> we must flag this.
472  //
473  if( currentMinimumStep == 0.0 )
474  {
475  if( currentSafety == 0.0 ) { fGeometryLimitedStep = true; }
476  }
477 
478  // Update the safety starting from the end-point,
479  // if it will become negative at the end-point.
480  //
481  if( currentSafety < fEndPointDistance )
482  {
483  if( particleCharge != 0.0 )
484  {
485  G4double endSafety =
486  fLinearNavigator->ComputeSafety( fTransportEndPosition) ;
487  currentSafety = endSafety ;
488  fPreviousSftOrigin = fTransportEndPosition ;
489  fPreviousSafety = currentSafety ;
490  fpSafetyHelper->SetCurrentSafety( currentSafety, fTransportEndPosition);
491 
492  // Because the Stepping Manager assumes it is from the start point,
493  // add the StepLength
494  //
495  currentSafety += fEndPointDistance ;
496 
497 #ifdef G4DEBUG_TRANSPORT
498  G4cout.precision(12) ;
499  G4cout << "***G4Transportation::AlongStepGPIL ** " << G4endl ;
500  G4cout << " Called Navigator->ComputeSafety at " << fTransportEndPosition
501  << " and it returned safety= " << endSafety << G4endl ;
502  G4cout << " Adding endpoint distance " << fEndPointDistance
503  << " to obtain pseudo-safety= " << currentSafety << G4endl ;
504  }
505  else
506  {
507  G4cout << "***G4Transportation::AlongStepGPIL ** " << G4endl ;
508  G4cout << " Avoiding call to ComputeSafety : " << G4endl;
509  G4cout << " charge = " << particleCharge << G4endl;
510  G4cout << " mag moment = " << magneticMoment << G4endl;
511 #endif
512  }
513  }
514 
515  fParticleChange.ProposeTrueStepLength(geometryStepLength) ;
516 
517  return geometryStepLength ;
518 }
519 
521 //
522 // Initialize ParticleChange (by setting all its members equal
523 // to corresponding members in G4Track)
524 
526  const G4Step& stepData )
527 {
528  static G4ThreadLocal G4int noCalls=0;
529  noCalls++;
530 
531  fParticleChange.Initialize(track) ;
532 
533  // Code for specific process
534  //
535  fParticleChange.ProposePosition(fTransportEndPosition) ;
536  fParticleChange.ProposeMomentumDirection(fTransportEndMomentumDir) ;
537  fParticleChange.ProposeEnergy(fTransportEndKineticEnergy) ;
538  fParticleChange.SetMomentumChanged(fMomentumChanged) ;
539 
540  fParticleChange.ProposePolarization(fTransportEndSpin);
541 
542  G4double deltaTime = 0.0 ;
543 
544  // Calculate Lab Time of Flight (ONLY if field Equations used it!)
545  // G4double endTime = fCandidateEndGlobalTime;
546  // G4double delta_time = endTime - startTime;
547 
548  G4double startTime = track.GetGlobalTime() ;
549 
550  if (!fEndGlobalTimeComputed)
551  {
552  // The time was not integrated .. make the best estimate possible
553  //
554  G4double initialVelocity = stepData.GetPreStepPoint()->GetVelocity();
555  G4double stepLength = track.GetStepLength();
556 
557  deltaTime= 0.0; // in case initialVelocity = 0
558  if ( initialVelocity > 0.0 ) { deltaTime = stepLength/initialVelocity; }
559 
560  fCandidateEndGlobalTime = startTime + deltaTime ;
561  fParticleChange.ProposeLocalTime( track.GetLocalTime() + deltaTime) ;
562  }
563  else
564  {
565  deltaTime = fCandidateEndGlobalTime - startTime ;
566  fParticleChange.ProposeGlobalTime( fCandidateEndGlobalTime ) ;
567  }
568 
569 
570  // Now Correct by Lorentz factor to get delta "proper" Time
571 
572  G4double restMass = track.GetDynamicParticle()->GetMass() ;
573  G4double deltaProperTime = deltaTime*( restMass/track.GetTotalEnergy() ) ;
574 
575  fParticleChange.ProposeProperTime(track.GetProperTime() + deltaProperTime) ;
576  //fParticleChange. ProposeTrueStepLength( track.GetStepLength() ) ;
577 
578  // If the particle is caught looping or is stuck (in very difficult
579  // boundaries) in a magnetic field (doing many steps)
580  // THEN this kills it ...
581  //
582  if ( fParticleIsLooping )
583  {
584  G4double endEnergy= fTransportEndKineticEnergy;
585 
586  if( (endEnergy < fThreshold_Important_Energy)
587  || (fNoLooperTrials >= fThresholdTrials ) )
588  {
589  // Kill the looping particle
590  //
591  fParticleChange.ProposeTrackStatus( fStopAndKill ) ;
592 
593  // 'Bare' statistics
594  fSumEnergyKilled += endEnergy;
595  if( endEnergy > fMaxEnergyKilled) { fMaxEnergyKilled= endEnergy; }
596 
597 #ifdef G4VERBOSE
598  if( (fVerboseLevel > 1) &&
599  ( endEnergy > fThreshold_Warning_Energy ) )
600  {
601  G4cout << " G4Transportation is killing track that is looping or stuck "
602  << G4endl
603  << " This track has " << track.GetKineticEnergy() / MeV
604  << " MeV energy." << G4endl;
605  G4cout << " Number of trials = " << fNoLooperTrials
606  << " No of calls to AlongStepDoIt = " << noCalls
607  << G4endl;
608  }
609 #endif
610  fNoLooperTrials=0;
611  }
612  else
613  {
614  fNoLooperTrials ++;
615 #ifdef G4VERBOSE
616  if( (fVerboseLevel > 2) )
617  {
618  G4cout << " G4Transportation::AlongStepDoIt(): Particle looping - "
619  << " Number of trials = " << fNoLooperTrials
620  << " No of calls to = " << noCalls
621  << G4endl;
622  }
623 #endif
624  }
625  }
626  else
627  {
628  fNoLooperTrials=0;
629  }
630 
631  // Another (sometimes better way) is to use a user-limit maximum Step size
632  // to alleviate this problem ..
633 
634  // Introduce smooth curved trajectories to particle-change
635  //
637  (fFieldPropagator->GimmeTrajectoryVectorAndForgetIt() );
638 
639  return &fParticleChange ;
640 }
641 
643 //
644 // This ensures that the PostStep action is always called,
645 // so that it can do the relocation if it is needed.
646 //
647 
650  G4double, // previousStepSize
651  G4ForceCondition* pForceCond )
652 {
653  fFieldExertedForce = false; // Not known
654  *pForceCond = Forced ;
655  return DBL_MAX ; // was kInfinity ; but convention now is DBL_MAX
656 }
657 
659 //
660 
662  const G4Step& )
663 {
664  G4TouchableHandle retCurrentTouchable ; // The one to return
665  G4bool isLastStep= false;
666 
667  // Initialize ParticleChange (by setting all its members equal
668  // to corresponding members in G4Track)
669  // fParticleChange.Initialize(track) ; // To initialise TouchableChange
670 
671  fParticleChange.ProposeTrackStatus(track.GetTrackStatus()) ;
672 
673  // If the Step was determined by the volume boundary,
674  // logically relocate the particle
675 
676  if(fGeometryLimitedStep)
677  {
678  // fCurrentTouchable will now become the previous touchable,
679  // and what was the previous will be freed.
680  // (Needed because the preStepPoint can point to the previous touchable)
681 
682  fLinearNavigator->SetGeometricallyLimitedStep() ;
683  fLinearNavigator->
684  LocateGlobalPointAndUpdateTouchableHandle( track.GetPosition(),
685  track.GetMomentumDirection(),
686  fCurrentTouchableHandle,
687  true ) ;
688  // Check whether the particle is out of the world volume
689  // If so it has exited and must be killed.
690  //
691  if( fCurrentTouchableHandle->GetVolume() == 0 )
692  {
693  fParticleChange.ProposeTrackStatus( fStopAndKill ) ;
694  }
695  retCurrentTouchable = fCurrentTouchableHandle ;
696  fParticleChange.SetTouchableHandle( fCurrentTouchableHandle ) ;
697 
698  // Update the Step flag which identifies the Last Step in a volume
699  if( !fFieldExertedForce )
700  isLastStep = fLinearNavigator->ExitedMotherVolume()
701  | fLinearNavigator->EnteredDaughterVolume() ;
702  else
703  isLastStep = fFieldPropagator->IsLastStepInVolume();
704  }
705  else // fGeometryLimitedStep is false
706  {
707  // This serves only to move the Navigator's location
708  //
709  fLinearNavigator->LocateGlobalPointWithinVolume( track.GetPosition() ) ;
710 
711  // The value of the track's current Touchable is retained.
712  // (and it must be correct because we must use it below to
713  // overwrite the (unset) one in particle change)
714  // It must be fCurrentTouchable too ??
715  //
716  fParticleChange.SetTouchableHandle( track.GetTouchableHandle() ) ;
717  retCurrentTouchable = track.GetTouchableHandle() ;
718 
719  isLastStep= false;
720  } // endif ( fGeometryLimitedStep )
721  fLastStepInVolume= isLastStep;
722 
723  fParticleChange.ProposeFirstStepInVolume(fFirstStepInVolume);
724  fParticleChange.ProposeLastStepInVolume(isLastStep);
725 
726  const G4VPhysicalVolume* pNewVol = retCurrentTouchable->GetVolume() ;
727  const G4Material* pNewMaterial = 0 ;
728  const G4VSensitiveDetector* pNewSensitiveDetector = 0 ;
729 
730  if( pNewVol != 0 )
731  {
732  pNewMaterial= pNewVol->GetLogicalVolume()->GetMaterial();
733  pNewSensitiveDetector= pNewVol->GetLogicalVolume()->GetSensitiveDetector();
734  }
735 
736  // ( <const_cast> pNewMaterial ) ;
737  // ( <const_cast> pNewSensitiveDetector) ;
738 
739  fParticleChange.SetMaterialInTouchable( (G4Material *) pNewMaterial ) ;
740  fParticleChange.SetSensitiveDetectorInTouchable( (G4VSensitiveDetector *) pNewSensitiveDetector ) ;
741 
742  const G4MaterialCutsCouple* pNewMaterialCutsCouple = 0;
743  if( pNewVol != 0 )
744  {
745  pNewMaterialCutsCouple=pNewVol->GetLogicalVolume()->GetMaterialCutsCouple();
746  }
747 
748  if( pNewVol!=0 && pNewMaterialCutsCouple!=0 && pNewMaterialCutsCouple->GetMaterial()!=pNewMaterial )
749  {
750  // for parametrized volume
751  //
752  pNewMaterialCutsCouple =
754  ->GetMaterialCutsCouple(pNewMaterial,
755  pNewMaterialCutsCouple->GetProductionCuts());
756  }
757  fParticleChange.SetMaterialCutsCoupleInTouchable( pNewMaterialCutsCouple );
758 
759  // temporarily until Get/Set Material of ParticleChange,
760  // and StepPoint can be made const.
761  // Set the touchable in ParticleChange
762  // this must always be done because the particle change always
763  // uses this value to overwrite the current touchable pointer.
764  //
765  fParticleChange.SetTouchableHandle(retCurrentTouchable) ;
766 
767  return &fParticleChange ;
768 }
769 
770 // New method takes over the responsibility to reset the state of G4Transportation
771 // object at the start of a new track or the resumption of a suspended track.
772 
773 void
775 {
777  fNewTrack= true;
778  fFirstStepInVolume= true;
779  fLastStepInVolume= false;
780 
781  // The actions here are those that were taken in AlongStepGPIL
782  // when track.GetCurrentStepNumber()==1
783 
784  // reset safety value and center
785  //
786  fPreviousSafety = 0.0 ;
787  fPreviousSftOrigin = G4ThreeVector(0.,0.,0.) ;
788 
789  // reset looping counter -- for motion in field
790  fNoLooperTrials= 0;
791  // Must clear this state .. else it depends on last track's value
792  // --> a better solution would set this from state of suspended track TODO ?
793  // Was if( aTrack->GetCurrentStepNumber()==1 ) { .. }
794 
795  // ChordFinder reset internal state
796  //
797  if( DoesGlobalFieldExist() )
798  {
799  fFieldPropagator->ClearPropagatorState();
800  // Resets all state of field propagator class (ONLY)
801  // including safety values (in case of overlaps and to wipe for first track).
802 
803  // G4ChordFinder* chordF= fFieldPropagator->GetChordFinder();
804  // if( chordF ) chordF->ResetStepEstimate();
805  }
806 
807  // Make sure to clear the chord finders of all fields (ie managers)
808  //
810  fieldMgrStore->ClearAllChordFindersState();
811 
812  // Update the current touchable handle (from the track's)
813  //
814  fCurrentTouchableHandle = aTrack->GetTouchableHandle();
815 
816  // Inform field propagator of new track
817  fFieldPropagator->PrepareNewTrack();
818 }
819 
822 {
823  G4bool lastValue= fUseMagneticMoment;
824  fUseMagneticMoment= useMoment;
825  G4CoupledTransportation::fUseMagneticMoment= useMoment;
826  return lastValue;
827 }
void SetMaterialInTouchable(G4Material *fMaterial)
static constexpr double perMillion
Definition: G4SIunits.hh:334
void SetTouchableHandle(const G4TouchableHandle &fTouchable)
const G4ThreeVector & GetPolarization() const
G4VParticleChange * PostStepDoIt(const G4Track &track, const G4Step &stepData)
G4SafetyHelper * GetSafetyHelper() const
G4double GetLocalTime() const
G4Material * GetMaterial() const
void SetMaterialCutsCoupleInTouchable(const G4MaterialCutsCouple *fMaterialCutsCouple)
virtual void SetChargeMomentumMass(G4ChargeState particleCharge, G4double MomentumXc, G4double MassXc2)=0
G4double GetProperTime() const
std::vector< G4ThreeVector > * GimmeTrajectoryVectorAndForgetIt() const
CLHEP::Hep3Vector G4ThreeVector
#define fNewTrack
void StartTracking(G4Track *aTrack)
G4double GetKineticEnergy() const
const G4DynamicParticle * GetDynamicParticle() const
G4VParticleChange * AlongStepDoIt(const G4Track &track, const G4Step &stepData)
const G4MagIntegratorStepper * GetStepper() const
const G4ThreeVector & GetPosition() const
const G4ThreeVector & GetMomentumDir() const
G4TrackStatus GetTrackStatus() const
G4ThreeVector GetSpin() const
virtual G4double ComputeStep(const G4ThreeVector &pGlobalPoint, const G4ThreeVector &pDirection, const G4double pCurrentProposedStepLength, G4double &pNewSafety)
Definition: G4Navigator.cc:747
G4Navigator * GetNavigatorForTracking() const
void ProposePolarization(G4double Px, G4double Py, G4double Pz)
G4ParticleDefinition * GetDefinition() const
G4double GetVelocity() const
#define G4ThreadLocal
Definition: tls.hh:89
int G4int
Definition: G4Types.hh:78
#define fFieldExertedForce
void ProposePosition(G4double x, G4double y, G4double z)
G4ReferenceCountedHandle< G4VTouchable > G4TouchableHandle
G4double GetTotalMomentum() const
G4StepPoint * GetPreStepPoint() const
virtual void Initialize(const G4Track &)
virtual void StartTracking(G4Track *)
Definition: G4VProcess.cc:101
virtual void ConfigureForTrack(const G4Track *)
void SetGeometricallyLimitedStep()
G4double GetKineticEnergy() const
G4ThreeVector GetPosition() const
G4GLOB_DLL std::ostream G4cout
G4double GetMass() const
bool G4bool
Definition: G4Types.hh:79
const G4ThreeVector & GetMomentumDirection() const
G4double GetCharge() const
void ProposeTrueStepLength(G4double truePathLength)
G4bool DoesGlobalFieldExist()
void ProposeFirstStepInVolume(G4bool flag)
void SetProcessSubType(G4int)
Definition: G4VProcess.hh:432
G4bool EnteredDaughterVolume() const
G4bool DoesFieldChangeEnergy() const
Definition: G4Step.hh:76
G4double GetGlobalTime() const
G4FieldManager * FindAndSetFieldManager(G4VPhysicalVolume *pCurrentPhysVol)
G4Transportation(G4int verbosityLevel=1)
const G4TouchableHandle & GetTouchableHandle() const
static G4TransportationManager * GetTransportationManager()
G4double PostStepGetPhysicalInteractionLength(const G4Track &, G4double previousStepSize, G4ForceCondition *pForceCond)
G4double ComputeStep(G4FieldTrack &pFieldTrack, G4double pCurrentProposedStepLength, G4double &pNewSafety, G4VPhysicalVolume *pPhysVol=0)
static G4ProductionCutsTable * GetProductionCutsTable()
const G4ThreeVector & GetMomentumDirection() const
G4double AlongStepGetPhysicalInteractionLength(const G4Track &track, G4double previousStepSize, G4double currentMinimumStep, G4double &currentSafety, G4GPILSelection *selection)
G4LogicalVolume * GetLogicalVolume() const
void ProposeProperTime(G4double finalProperTime)
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
static G4FieldManagerStore * GetInstance()
void SetSensitiveDetectorInTouchable(G4VSensitiveDetector *fSensitiveDetector)
#define fPreviousSftOrigin
G4bool IsParticleLooping() const
G4double GetLabTimeOfFlight() const
virtual G4VPhysicalVolume * GetVolume(G4int depth=0) const
Definition: G4VTouchable.cc:44
void SetPointerToVectorOfAuxiliaryPoints(std::vector< G4ThreeVector > *theNewVectorPointer)
G4bool IsGravityActive() const
Definition: G4Field.hh:98
void ProposeGlobalTime(G4double t)
G4ChordFinder * GetChordFinder()
G4FieldManager * GetCurrentFieldManager()
#define fPreviousSafety
G4VParticleChange * pParticleChange
Definition: G4VProcess.hh:283
void ProposeEnergy(G4double finalEnergy)
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
double mag2() const
G4VPhysicalVolume * GetVolume() const
const G4MaterialCutsCouple * GetMaterialCutsCouple() const
G4double GetTotalEnergy() const
G4double GetPDGSpin() const
void SetCurrentSafety(G4double val, const G4ThreeVector &pos)
#define G4endl
Definition: G4ios.hh:61
G4bool ExitedMotherVolume() const
static constexpr double MeV
Definition: G4SIunits.hh:214
void ProposeLastStepInVolume(G4bool flag)
static constexpr double perThousand
Definition: G4SIunits.hh:333
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
T sqr(const T &x)
Definition: templates.hh:145
G4double GetPDGMagneticMoment() const
double G4double
Definition: G4Types.hh:76
void ProposeTrackStatus(G4TrackStatus status)
const G4Field * GetDetectorField() const
virtual G4double ComputeSafety(const G4ThreeVector &globalpoint, const G4double pProposedMaxLength=DBL_MAX, const G4bool keepState=true)
G4ForceCondition
G4MagInt_Driver * GetIntegrationDriver()
void SetMomentumChanged(G4bool b)
G4PropagatorInField * GetPropagatorInField() const
G4ProductionCuts * GetProductionCuts() const
#define DBL_MAX
Definition: templates.hh:83
G4VSensitiveDetector * GetSensitiveDetector() const
G4double GetMagneticMoment() const
void ProposeLocalTime(G4double t)
G4double GetStepLength() const
virtual void LocateGlobalPointWithinVolume(const G4ThreeVector &position)
Definition: G4Navigator.cc:582
G4GPILSelection
const G4Material * GetMaterial() const
G4GLOB_DLL std::ostream G4cerr
static G4bool EnableUseMagneticMoment(G4bool useMoment=true)