ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4Transportation.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4Transportation.cc
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 //
28 // ------------------------------------------------------------
29 // GEANT 4 include file implementation
30 //
31 // ------------------------------------------------------------
32 //
33 // This class is a process responsible for the transportation of
34 // a particle, ie the geometrical propagation that encounters the
35 // geometrical sub-volumes of the detectors.
36 //
37 // It is also tasked with the key role of proposing the "isotropic safety",
38 // which will be used to update the post-step point's safety.
39 //
40 // =======================================================================
41 // Created: 19 March 1997, J. Apostolakis
42 // =======================================================================
43 
44 #include "G4Transportation.hh"
47 
48 #include "G4PhysicalConstants.hh"
49 #include "G4SystemOfUnits.hh"
50 #include "G4ProductionCutsTable.hh"
51 #include "G4ParticleTable.hh"
52 
53 #include "G4ChargeState.hh"
54 #include "G4EquationOfMotion.hh"
55 
56 #include "G4FieldManagerStore.hh"
58 
60 
64 
66 //
67 // Constructor
68 
70  : G4VProcess( G4String("Transportation"), fTransportation ),
72  fPreviousSftOrigin( 0.,0.,0. ),
73  fPreviousSafety( 0.0 ),
74  fEndPointDistance( -1.0 ),
75  fShortStepOptimisation( false ) // Old default: true (=fast short steps)
76 {
77  SetProcessSubType(static_cast<G4int>(TRANSPORTATION));
78  pParticleChange= &fParticleChange; // Required to conform to G4VProcess
79  SetVerboseLevel(verbosity);
80 
81  G4TransportationManager* transportMgr ;
82 
84 
85  fLinearNavigator = transportMgr->GetNavigatorForTracking() ;
86 
87  fFieldPropagator = transportMgr->GetPropagatorInField() ;
88 
89  fpSafetyHelper = transportMgr->GetSafetyHelper(); // New
90 
91  fpLogger = new G4TransportationLogger("G4Transportation", verbosity);
92 
94  // Use the old defaults: Warning = 100 MeV, Important = 250 MeV, No Trials = 10;
95 
97  // Should be done by Set methods in SetHighLooperThresholds -- making sure
98 
99  // Cannot determine whether a field exists here, as it would
100  // depend on the relative order of creating the detector's
101  // field and this process. That order is not guaranted.
103  // This value must be updated using DoesAnyFieldExist() at least at the
104  // start of each Run -- for now this is at the Start of every Track. TODO
105 
106  static G4ThreadLocal G4TouchableHandle* pNullTouchableHandle = 0;
107  if ( !pNullTouchableHandle)
108  {
109  pNullTouchableHandle = new G4TouchableHandle;
110  }
111  fCurrentTouchableHandle = *pNullTouchableHandle;
112  // Points to (G4VTouchable*) 0
113 
114 
115 #ifdef G4VERBOSE
116  if( verboseLevel > 0)
117  {
118  G4cout << " G4Transportation constructor> set fShortStepOptimisation to ";
119  if ( fShortStepOptimisation ) { G4cout << "true" << G4endl; }
120  else { G4cout << "false" << G4endl; }
121  }
122 #endif
123 }
124 
126 
128 {
129  if( fSumEnergyKilled > 0.0 )
130  {
132  }
133  delete fpLogger;
134 }
135 
137 
138 void
139 G4Transportation::PrintStatistics( std::ostream& outStr) const
140 {
141  outStr << " G4Transportation: Statistics for looping particles " << G4endl;
142  if( fSumEnergyKilled > 0.0 || fNumLoopersKilled > 0 )
143  {
144  outStr << " Sum of energy of looping tracks killed: "
145  << fSumEnergyKilled / CLHEP::MeV << " MeV "
146  << " from " << fNumLoopersKilled << " tracks " << G4endl
147  << " Sum of energy of non-electrons : "
149  << " from " << fNumLoopersKilled_NonElectron << " tracks "
150  << G4endl;
151  outStr << " Max energy of *any type* looper killed: " << fMaxEnergyKilled
152  << " its PDG was " << fMaxEnergyKilledPDG << G4endl;
153  if( fMaxEnergyKilled_NonElectron > 0.0 )
154  {
155  outStr << " Max energy of non-electron looper killed: "
157  << " its PDG was " << fMaxEnergyKilled_NonElecPDG << G4endl;
158  }
159  if( fMaxEnergySaved > 0.0 )
160  {
161  outStr << " Max energy of loopers 'saved': " << fMaxEnergySaved << G4endl;
162  outStr << " Sum of energy of loopers 'saved': "
163  << fSumEnergySaved << G4endl;
164  outStr << " Sum of energy of unstable loopers 'saved': "
166  }
167  }
168  else
169  {
170  outStr << " No looping tracks found or killed. " << G4endl;
171  }
172 }
173 
175 //
176 // Responsibilities:
177 // Find whether the geometry limits the Step, and to what length
178 // Calculate the new value of the safety and return it.
179 // Store the final time, position and momentum.
180 
183  G4double, // previousStepSize
184  G4double currentMinimumStep,
185  G4double& currentSafety,
186  G4GPILSelection* selection )
187 {
188  G4double geometryStepLength= -1.0, newSafety= -1.0;
190 
191  // Initial actions moved to StartTrack()
192  // --------------------------------------
193  // Note: in case another process changes touchable handle
194  // it will be necessary to add here (for all steps)
195  // fCurrentTouchableHandle = aTrack->GetTouchableHandle();
196 
197  // GPILSelection is set to defaule value of CandidateForSelection
198  // It is a return value
199  //
200  *selection = CandidateForSelection ;
201 
203  fLastStepInVolume= false;
204  fNewTrack = false;
205 
207 
208  // Get initial Energy/Momentum of the track
209  //
210  const G4DynamicParticle* pParticle = track.GetDynamicParticle() ;
211  const G4ParticleDefinition* pParticleDef = pParticle->GetDefinition() ;
212  G4ThreeVector startMomentumDir = pParticle->GetMomentumDirection() ;
213  G4ThreeVector startPosition = track.GetPosition() ;
214 
215  // The Step Point safety can be limited by other geometries and/or the
216  // assumptions of any process - it's not always the geometrical safety.
217  // We calculate the starting point's isotropic safety here.
218  //
219  G4ThreeVector OriginShift = startPosition - fPreviousSftOrigin ;
220  G4double MagSqShift = OriginShift.mag2() ;
221  if( MagSqShift >= sqr(fPreviousSafety) )
222  {
223  currentSafety = 0.0 ;
224  }
225  else
226  {
227  currentSafety = fPreviousSafety - std::sqrt(MagSqShift) ;
228  }
229 
230  // Is the particle charged or has it a magnetic moment?
231  //
232  G4double particleCharge = pParticle->GetCharge() ;
233  G4double magneticMoment = pParticle->GetMagneticMoment() ;
234  G4double restMass = pParticle->GetMass() ;
235 
237 
238  // There is no need to locate the current volume. It is Done elsewhere:
239  // On track construction
240  // By the tracking, after all AlongStepDoIts, in "Relocation"
241 
242  // Check if the particle has a force, EM or gravitational, exerted on it
243  //
244  G4FieldManager* fieldMgr=0;
245  G4bool fieldExertsForce = false ;
246 
247  fieldMgr = fFieldPropagator->FindAndSetFieldManager( track.GetVolume() );
248  G4bool eligibleEM = (particleCharge != 0.0)
249  || ( fUseMagneticMoment && (magneticMoment != 0.0) );
250  G4bool eligibleGrav = fUseGravity && (restMass != 0.0) ;
251 
252  if( (fieldMgr!=nullptr) && (eligibleEM||eligibleGrav) )
253  {
254  // User can configure the field Manager for this track
255  fieldMgr->ConfigureForTrack( &track );
256  // Called here to allow a transition from no-field pointer
257  // to finite field (non-zero pointer).
258 
259  // If the field manager has no field ptr, the field is zero
260  // by definition ( = there is no field ! )
261  const G4Field* ptrField= fieldMgr->GetDetectorField();
262  if( ptrField )
263  {
264  fieldExertsForce = eligibleEM
265  || ( eligibleGrav && ptrField->IsGravityActive() );
266  }
267  // || (gravityOn && (restMass != 0.0)) )
268 
269  }
270  fFieldExertedForce = fieldExertsForce;
271 
272  if( !fieldExertsForce )
273  {
274  G4double linearStepLength ;
275  if( fShortStepOptimisation && (currentMinimumStep <= currentSafety) )
276  {
277  // The Step is guaranteed to be taken
278  //
279  geometryStepLength = currentMinimumStep ;
281  }
282  else
283  {
284  // Find whether the straight path intersects a volume
285  //
286  linearStepLength = fLinearNavigator->ComputeStep( startPosition,
287  startMomentumDir,
288  currentMinimumStep,
289  newSafety) ;
290  // Remember last safety origin & value.
291  //
292  fPreviousSftOrigin = startPosition ;
293  fPreviousSafety = newSafety ;
294  fpSafetyHelper->SetCurrentSafety( newSafety, startPosition);
295 
296  currentSafety = newSafety ;
297 
298  fGeometryLimitedStep= (linearStepLength <= currentMinimumStep);
300  {
301  // The geometry limits the Step size (an intersection was found.)
302  geometryStepLength = linearStepLength ;
303  }
304  else
305  {
306  // The full Step is taken.
307  geometryStepLength = currentMinimumStep ;
308  }
309  }
310  fEndPointDistance = geometryStepLength ;
311 
312  // Calculate final position
313  //
314  fTransportEndPosition = startPosition+geometryStepLength*startMomentumDir;
315 
316  // Momentum direction, energy and polarisation are unchanged by transport
317  //
318  fTransportEndMomentumDir = startMomentumDir ;
324  }
325  else // A field exerts force
326  {
327  G4double momentumMagnitude = pParticle->GetTotalMomentum() ;
328  G4ThreeVector EndUnitMomentum ;
329  G4double lengthAlongCurve ;
330 
331  // The charge can change (dynamic)
332  //
333  G4ChargeState chargeState(particleCharge,
334  magneticMoment,
335  pParticleDef->GetPDGSpin() );
336 
337  auto equationOfMotion = fFieldPropagator->GetCurrentEquationOfMotion();
338 
339  equationOfMotion->SetChargeMomentumMass( chargeState,
340  momentumMagnitude,
341  restMass);
342 
343  G4FieldTrack aFieldTrack = G4FieldTrack( startPosition,
344  track.GetGlobalTime(), // Lab.
345  track.GetMomentumDirection(),
346  track.GetKineticEnergy(),
347  restMass,
348  particleCharge,
349  track.GetPolarization(),
350  pParticleDef->GetPDGMagneticMoment(),
351  0.0, // Length along track
352  pParticleDef->GetPDGSpin() );
353 
354  if( currentMinimumStep > 0 )
355  {
356  // Do the Transport in the field (non recti-linear)
357  //
358  lengthAlongCurve = fFieldPropagator->ComputeStep( aFieldTrack,
359  currentMinimumStep,
360  currentSafety,
361  track.GetVolume(),
363 
365  //
366  // It is possible that step was reduced in PropagatorInField due to
367  // previous zero steps. To cope with case that reduced step is taken
368  // in full, we must rely on PiF to obtain this value
369 
370  geometryStepLength = std::min( lengthAlongCurve, currentMinimumStep );
371 
372  // Remember last safety origin & value.
373  //
374  fPreviousSftOrigin = startPosition ;
375  fPreviousSafety = currentSafety ;
376  fpSafetyHelper->SetCurrentSafety( currentSafety, startPosition);
377  }
378  else
379  {
380  geometryStepLength = lengthAlongCurve= 0.0 ;
382  }
383 
384  // Get the End-Position and End-Momentum (Dir-ection)
385  //
386  fTransportEndPosition = aFieldTrack.GetPosition() ;
387 
388  fTransportEndSpin = aFieldTrack.GetSpin();
390  fEndPointDistance = (fTransportEndPosition - startPosition).mag() ;
391 
392  // Momentum: Magnitude and direction can be changed too now ...
393  //
395  fTransportEndMomentumDir = aFieldTrack.GetMomentumDir() ;
396 
398 
400  {
401  // If the field can change energy, then the time must be integrated
402  // - so this should have been updated
403  //
405  fEndGlobalTimeComputed = true;
406 
407  // was ( fCandidateEndGlobalTime != track.GetGlobalTime() );
408  // a cleaner way is to have FieldTrack knowing whether time is updated.
409  }
410  else
411  {
412  // The energy should be unchanged by field transport,
413  // - so the time changed will be calculated elsewhere
414  //
415  fEndGlobalTimeComputed = false;
416 
417  // Check that the integration preserved the energy
418  // - and if not correct this!
419  G4double startEnergy= track.GetKineticEnergy();
421 
422  static G4ThreadLocal G4int no_inexact_steps=0, no_large_ediff;
423  G4double absEdiff = std::fabs(startEnergy- endEnergy);
424  if( absEdiff > perMillion * endEnergy )
425  {
426  no_inexact_steps++;
427  // Possible statistics keeping here ...
428  }
429  if( verboseLevel > 1 )
430  {
431  if( std::fabs(startEnergy- endEnergy) > perThousand * endEnergy )
432  {
433  static G4ThreadLocal G4int no_warnings= 0, warnModulo=1,
434  moduloFactor= 10;
435  no_large_ediff ++;
436  if( (no_large_ediff% warnModulo) == 0 )
437  {
438  no_warnings++;
439  std::ostringstream message;
440  message << "Energy change in Step is above 1^-3 relative value. "
441  << G4endl
442  << " Relative change in 'tracking' step = "
443  << std::setw(15) << (endEnergy-startEnergy)/startEnergy
444  << G4endl
445  << " Starting E= " << std::setw(12) << startEnergy / MeV
446  << " MeV " << G4endl
447  << " Ending E= " << std::setw(12) << endEnergy / MeV
448  << " MeV " << G4endl
449  << "Energy has been corrected -- however, review"
450  << " field propagation parameters for accuracy." << G4endl;
451  if ( (verboseLevel > 2 ) || (no_warnings<4)
452  || (no_large_ediff == warnModulo * moduloFactor) )
453  {
454  message << "These include EpsilonStepMax(/Min) in G4FieldManager " << G4endl
455  << "which determine fractional error per step for integrated quantities. " << G4endl
456  << "Note also the influence of the permitted number of integration steps."
457  << G4endl;
458  }
459  message << "Bad 'endpoint'. Energy change detected and corrected."
460  << G4endl
461  << "Has occurred already " << no_large_ediff << " times.";
462  G4Exception("G4Transportation::AlongStepGetPIL()",
463  "EnergyChange", JustWarning, message);
464  if( no_large_ediff == warnModulo * moduloFactor )
465  {
466  warnModulo *= moduloFactor;
467  }
468  }
469  }
470  } // end of if (verboseLevel)
471 
472  // Correct the energy for fields that conserve it
473  // This - hides the integration error
474  // - but gives a better physical answer
475  //
477  }
478 
479  }
480 
481  // If we are asked to go a step length of 0, and we are on a boundary
482  // then a boundary will also limit the step -> we must flag this.
483  //
484  if( currentMinimumStep == 0.0 )
485  {
486  if( currentSafety == 0.0 ) { fGeometryLimitedStep = true; }
487  }
488 
489  // Update the safety starting from the end-point,
490  // if it will become negative at the end-point.
491  //
492  if( currentSafety < fEndPointDistance )
493  {
494  if( particleCharge != 0.0 )
495  {
496  G4double endSafety =
498  currentSafety = endSafety ;
499  fPreviousSftOrigin = fTransportEndPosition ;
500  fPreviousSafety = currentSafety ;
502 
503  // Because the Stepping Manager assumes it is from the start point,
504  // add the StepLength
505  //
506  currentSafety += fEndPointDistance ;
507 
508 #ifdef G4DEBUG_TRANSPORT
509  G4cout.precision(12) ;
510  G4cout << "***G4Transportation::AlongStepGPIL ** " << G4endl ;
511  G4cout << " Called Navigator->ComputeSafety at " << fTransportEndPosition
512  << " and it returned safety= " << endSafety << G4endl ;
513  G4cout << " Adding endpoint distance " << fEndPointDistance
514  << " to obtain pseudo-safety= " << currentSafety << G4endl ;
515  }
516  else
517  {
518  G4cout << "***G4Transportation::AlongStepGPIL ** " << G4endl ;
519  G4cout << " Avoiding call to ComputeSafety : " << G4endl;
520  G4cout << " charge = " << particleCharge << G4endl;
521  G4cout << " mag moment = " << magneticMoment << G4endl;
522 #endif
523  }
524  }
525 
526  fParticleChange.ProposeTrueStepLength(geometryStepLength) ;
527 
528  return geometryStepLength ;
529 }
530 
532 //
533 // Initialize ParticleChange (by setting all its members equal
534 // to corresponding members in G4Track)
535 
537  const G4Step& stepData )
538 {
539  static G4ThreadLocal G4long noCallsASDI=0;
540  const char *methodName= "AlongStepDoIt";
541  noCallsASDI++;
542 
543  fParticleChange.Initialize(track) ;
544 
545  // Code for specific process
546  //
551 
553 
554  G4double deltaTime = 0.0 ;
555 
556  // Calculate Lab Time of Flight (ONLY if field Equations used it!)
557  // G4double endTime = fCandidateEndGlobalTime;
558  // G4double delta_time = endTime - startTime;
559 
560  G4double startTime = track.GetGlobalTime() ;
561 
563  {
564  // The time was not integrated .. make the best estimate possible
565  //
566  G4double initialVelocity = stepData.GetPreStepPoint()->GetVelocity();
567  G4double stepLength = track.GetStepLength();
568 
569  deltaTime= 0.0; // in case initialVelocity = 0
570  if ( initialVelocity > 0.0 ) { deltaTime = stepLength/initialVelocity; }
571 
572  fCandidateEndGlobalTime = startTime + deltaTime ;
573  fParticleChange.ProposeLocalTime( track.GetLocalTime() + deltaTime) ;
574  }
575  else
576  {
577  deltaTime = fCandidateEndGlobalTime - startTime ;
579  }
580 
581 
582  // Now Correct by Lorentz factor to get delta "proper" Time
583 
584  G4double restMass = track.GetDynamicParticle()->GetMass() ;
585  G4double deltaProperTime = deltaTime*( restMass/track.GetTotalEnergy() ) ;
586 
587  fParticleChange.ProposeProperTime(track.GetProperTime() + deltaProperTime) ;
588  //fParticleChange.ProposeTrueStepLength( track.GetStepLength() ) ;
589 
590  // If the particle is caught looping or is stuck (in very difficult
591  // boundaries) in a magnetic field (doing many steps) THEN can kill it ...
592  //
593  if ( fParticleIsLooping )
594  {
596  fNoLooperTrials ++;
597  auto particleType= track.GetDynamicParticle()->GetParticleDefinition();
598 
599  G4bool stable = particleType->GetPDGStable();
600  G4bool candidateForEnd = (endEnergy < fThreshold_Important_Energy)
602  G4bool unstableAndKillable = !stable && ( fAbandonUnstableTrials != 0);
603  G4bool unstableForEnd = (endEnergy < fThreshold_Important_Energy)
605  if( (candidateForEnd && stable) || (unstableAndKillable && unstableForEnd) )
606  {
607  // Kill the looping particle
608  //
610  G4int particlePDG= particleType->GetPDGEncoding();
611  const G4int electronPDG= 11; // G4Electron::G4Electron()->GetPDGEncoding();
612 
613  // Simple statistics
614  fSumEnergyKilled += endEnergy;
615  fSumEnerSqKilled = endEnergy * endEnergy;
617 
618  if( endEnergy > fMaxEnergyKilled ) {
619  fMaxEnergyKilled = endEnergy;
620  fMaxEnergyKilledPDG = particlePDG;
621  }
622  if( particleType->GetPDGEncoding() != electronPDG )
623  {
624  fSumEnergyKilled_NonElectron += endEnergy;
625  fSumEnerSqKilled_NonElectron += endEnergy * endEnergy;
627 
628  if( endEnergy > fMaxEnergyKilled_NonElectron )
629  {
630  fMaxEnergyKilled_NonElectron = endEnergy;
631  fMaxEnergyKilled_NonElecPDG = particlePDG;
632  }
633  }
634 
636  {
637  fpLogger->ReportLoopingTrack( track, stepData, fNoLooperTrials,
638  noCallsASDI, methodName );
639  }
640  fNoLooperTrials=0;
641  }
642  else
643  {
645  if( fNoLooperTrials == 1 ) {
646  fSumEnergySaved += endEnergy;
647  if ( !stable )
648  fSumEnergyUnstableSaved += endEnergy;
649  }
650 #ifdef G4VERBOSE
651  if( verboseLevel > 2 && ! fSilenceLooperWarnings )
652  {
653  G4cout << " " << methodName
654  << " Particle is looping but is saved ..." << G4endl
655  << " Number of trials = " << fNoLooperTrials << G4endl
656  << " No of calls to = " << noCallsASDI << G4endl;
657  }
658 #endif
659  }
660  }
661  else
662  {
663  fNoLooperTrials=0;
664  }
665 
666  // Another (sometimes better way) is to use a user-limit maximum Step size
667  // to alleviate this problem ..
668 
669  // Introduce smooth curved trajectories to particle-change
670  //
673 
674  return &fParticleChange ;
675 }
676 
678 //
679 // This ensures that the PostStep action is always called,
680 // so that it can do the relocation if it is needed.
681 //
682 
685  G4double, // previousStepSize
686  G4ForceCondition* pForceCond )
687 {
688  fFieldExertedForce = false; // Not known
689  *pForceCond = Forced ;
690  return DBL_MAX ; // was kInfinity ; but convention now is DBL_MAX
691 }
692 
694 //
695 
697  const G4Step& )
698 {
699  G4TouchableHandle retCurrentTouchable ; // The one to return
700  G4bool isLastStep= false;
701 
702  // Initialize ParticleChange (by setting all its members equal
703  // to corresponding members in G4Track)
704  // fParticleChange.Initialize(track) ; // To initialise TouchableChange
705 
707 
708  // If the Step was determined by the volume boundary,
709  // logically relocate the particle
710 
712  {
713  // fCurrentTouchable will now become the previous touchable,
714  // and what was the previous will be freed.
715  // (Needed because the preStepPoint can point to the previous touchable)
716 
719  LocateGlobalPointAndUpdateTouchableHandle( track.GetPosition(),
720  track.GetMomentumDirection(),
722  true ) ;
723  // Check whether the particle is out of the world volume
724  // If so it has exited and must be killed.
725  //
726  if( fCurrentTouchableHandle->GetVolume() == 0 )
727  {
729  }
730  retCurrentTouchable = fCurrentTouchableHandle ;
732 
733  // Update the Step flag which identifies the Last Step in a volume
734  if( !fFieldExertedForce )
735  isLastStep = fLinearNavigator->ExitedMotherVolume()
737  else
738  isLastStep = fFieldPropagator->IsLastStepInVolume();
739  }
740  else // fGeometryLimitedStep is false
741  {
742  // This serves only to move the Navigator's location
743  //
745 
746  // The value of the track's current Touchable is retained.
747  // (and it must be correct because we must use it below to
748  // overwrite the (unset) one in particle change)
749  // It must be fCurrentTouchable too ??
750  //
752  retCurrentTouchable = track.GetTouchableHandle() ;
753 
754  isLastStep= false;
755  } // endif ( fGeometryLimitedStep )
756  fLastStepInVolume= isLastStep;
757 
760 
761  const G4VPhysicalVolume* pNewVol = retCurrentTouchable->GetVolume() ;
762  const G4Material* pNewMaterial = 0 ;
763  const G4VSensitiveDetector* pNewSensitiveDetector = 0 ;
764 
765  if( pNewVol != 0 )
766  {
767  pNewMaterial= pNewVol->GetLogicalVolume()->GetMaterial();
768  pNewSensitiveDetector= pNewVol->GetLogicalVolume()->GetSensitiveDetector();
769  }
770 
773 
774  const G4MaterialCutsCouple* pNewMaterialCutsCouple = 0;
775  if( pNewVol != 0 )
776  {
777  pNewMaterialCutsCouple=pNewVol->GetLogicalVolume()->GetMaterialCutsCouple();
778  }
779 
780  if ( pNewVol!=0 && pNewMaterialCutsCouple!=0
781  && pNewMaterialCutsCouple->GetMaterial()!=pNewMaterial )
782  {
783  // for parametrized volume
784  //
785  pNewMaterialCutsCouple =
787  ->GetMaterialCutsCouple(pNewMaterial,
788  pNewMaterialCutsCouple->GetProductionCuts());
789  }
790  fParticleChange.SetMaterialCutsCoupleInTouchable( pNewMaterialCutsCouple );
791 
792  // temporarily until Get/Set Material of ParticleChange,
793  // and StepPoint can be made const.
794  // Set the touchable in ParticleChange
795  // this must always be done because the particle change always
796  // uses this value to overwrite the current touchable pointer.
797  //
798  fParticleChange.SetTouchableHandle(retCurrentTouchable) ;
799 
800  return &fParticleChange ;
801 }
802 
804 // New method takes over the responsibility to reset the state of
805 // G4Transportation object at the start of a new track or the resumption
806 // of a suspended track.
807 //
808 
809 void
811 {
813  fNewTrack= true;
814  fFirstStepInVolume= true;
815  fLastStepInVolume= false;
816 
817  // The actions here are those that were taken in AlongStepGPIL
818  // when track.GetCurrentStepNumber()==1
819 
820  // Whether field exists should be determined at run level -- TODO
822 
823  // reset safety value and center
824  //
825  fPreviousSafety = 0.0 ;
826  fPreviousSftOrigin = G4ThreeVector(0.,0.,0.) ;
827 
828  // reset looping counter -- for motion in field
829  fNoLooperTrials= 0;
830  // Must clear this state .. else it depends on last track's value
831  // --> a better solution would set this from state of suspended track TODO ?
832  // Was if( aTrack->GetCurrentStepNumber()==1 ) { .. }
833 
834  // ChordFinder reset internal state
835  //
837  {
839  // Resets all state of field propagator class (ONLY) including safety
840  // values (in case of overlaps and to wipe for first track).
841  }
842 
843  // Make sure to clear the chord finders of all fields (i.e. managers)
844  //
846  fieldMgrStore->ClearAllChordFindersState();
847 
848  // Update the current touchable handle (from the track's)
849  //
851 
852  // Inform field propagator of new track
853  //
855 }
856 
858 //
859 
861 {
862  G4bool lastValue= fUseMagneticMoment;
863  fUseMagneticMoment= useMoment;
865  return lastValue;
866 }
867 
869 //
870 
872 {
873  G4bool lastValue= fUseGravity;
874  fUseGravity= useGravity;
876  return lastValue;
877 }
878 
880 //
881 // Supress (or not) warnings about 'looping' particles
882 
884 {
885  fSilenceLooperWarnings= val; // Flag to *Supress* all 'looper' warnings
886  // G4CoupledTransportation::fSilenceLooperWarnings= val;
887 }
888 
890 //
892 {
893  return fSilenceLooperWarnings;
894 }
895 
896 
898 //
900 {
901  // Restores the old high values -- potentially appropriate for energy-frontier
902  // HEP experiments.
903  // Caution: All tracks with E < 100 MeV that are found to loop are
904  SetThresholdWarningEnergy( 100.0 * CLHEP::MeV ); // Warn above this energy
905  SetThresholdImportantEnergy( 250.0 * CLHEP::MeV ); // Extra trial above this En
906 
907  G4int maxTrials = 10;
908  SetThresholdTrials( maxTrials );
909 
910  PushThresholdsToLogger(); // Again, to be sure
912 }
913 
915 void G4Transportation::SetLowLooperThresholds() // Values for low-E applications
916 {
917  // These values were the default in Geant4 10.5 - beta
918  SetThresholdWarningEnergy( 1.0 * CLHEP::keV ); // Warn above this En
919  SetThresholdImportantEnergy( 1.0 * CLHEP::MeV ); // Extra trials above it
920 
921  G4int maxTrials = 30; // A new value - was 10
922  SetThresholdTrials( maxTrials );
923 
924  PushThresholdsToLogger(); // Again, to be sure
926 }
927 
929 //
930 void
931 G4Transportation::ReportMissingLogger( const char* methodName )
932 {
933  const char* message= "Logger object missing from G4Transportation object";
934  G4String classAndMethod= G4String("G4Transportation") + G4String( methodName );
935  G4Exception(classAndMethod, "Missing Logger", JustWarning, message);
936 }
937 
938 
940 //
941 void
943 {
944  PushThresholdsToLogger(); // To be absolutely certain they are in sync
945  fpLogger->ReportLooperThresholds("G4Transportation");
946 }
947 
949 //
950 void G4Transportation::ProcessDescription(std::ostream& outStr) const
951 
952 // StreamInfo(std::ostream& out, const G4ParticleDefinition& part, G4bool rst) const
953 
954 {
955  G4String indent = " "; // : "");
956  G4int oldPrec= outStr.precision(6);
957  // outStr << std::setprecision(6);
958  outStr << G4endl << indent << GetProcessName() << ": ";
959 
960  outStr << " Parameters for looping particles: " << G4endl
961  << " warning-E = " << fThreshold_Warning_Energy / CLHEP::MeV << " MeV " << G4endl
962  << " important E = " << fThreshold_Important_Energy / CLHEP::MeV << " MeV " << G4endl
963  << " thresholdTrials " << fThresholdTrials << G4endl;
964  outStr.precision(oldPrec);
965 }