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