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