Geant4_10
G4ITTransportation.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 // $Id: G4ITTransportation.cc 71125 2013-06-11 15:39:09Z gcosmo $
27 //
31 //
32 // Contact : Mathieu Karamitros (kara (AT) cenbg . in2p3 . fr)
33 //
34 // History :
35 // -----------
36 // =======================================================================
37 // Modified:
38 // 28 Oct 2011, P.Gumpl./J.Ap: Detect gravity field, use magnetic moment
39 // 20 Nov 2008, J.Apostolakis: Push safety to helper - after ComputeSafety
40 // 9 Nov 2007, J.Apostolakis: Flag for short steps, push safety to helper
41 // 19 Jan 2006, P.MoraDeFreitas: Fix for suspended tracks (StartTracking)
42 // 11 Aug 2004, M.Asai: Add G4VSensitiveDetector* for updating stepPoint.
43 // 21 June 2003, J.Apostolakis: Calling field manager with
44 // track, to enable it to configure its accuracy
45 // 13 May 2003, J.Apostolakis: Zero field areas now taken into
46 // account correclty in all cases (thanks to W Pokorski).
47 // 29 June 2001, J.Apostolakis, D.Cote-Ahern, P.Gumplinger:
48 // correction for spin tracking
49 // 20 Febr 2001, J.Apostolakis: update for new FieldTrack
50 // 22 Sept 2000, V.Grichine: update of Kinetic Energy
51 // ---------------------------------------------------
52 // 10 Oct 2011, M.Karamitros: G4ITTransportation created
53 // Created: 19 March 1997, J. Apostolakis
54 // =======================================================================
55 //
56 // -------------------------------------------------------------------
57 
58 #include "G4ITTransportation.hh"
59 #include "G4SystemOfUnits.hh"
62 #include "G4ProductionCutsTable.hh"
63 #include "G4ParticleTable.hh"
64 #include "G4ITNavigator.hh"
65 #include "G4PropagatorInField.hh"
66 #include "G4FieldManager.hh"
67 #include "G4ChordFinder.hh"
68 #include "G4SafetyHelper.hh"
69 #include "G4FieldManagerStore.hh"
70 
71 #include "G4UnitsTable.hh"
72 
74 
75 #ifndef State
76 #define State(theXInfo) (fTransportationState->theXInfo)
77 #endif
78 
79 //#define G4DEBUG_TRANSPORT 1
80 
83  InitProcessState(fTransportationState, fpState),
84  fThreshold_Warning_Energy( 100 * MeV ),
85  fThreshold_Important_Energy( 250 * MeV ),
86  fThresholdTrials( 10 ),
87  fUnimportant_Energy( 1 * MeV ), // Not used
88  fSumEnergyKilled( 0.0 ), fMaxEnergyKilled( 0.0 ),
89  fShortStepOptimisation(false), // Old default: true (=fast short steps)
90  fVerboseLevel( verbose )
91 {
93  G4TransportationManager* transportMgr ;
95  G4ITTransportationManager* ITtransportMgr ;
97  fLinearNavigator = ITtransportMgr->GetNavigatorForTracking() ;
98  fFieldPropagator = transportMgr->GetPropagatorInField() ;
99  fpSafetyHelper = 0; // transportMgr->GetSafetyHelper(); // New
100 
101  // Cannot determine whether a field exists here, as it would
102  // depend on the relative order of creating the detector's
103  // field and this process. That order is not guaranted.
104  // Instead later the method DoesGlobalFieldExist() is called
105 
106  enableAtRestDoIt = false;
107  enableAlongStepDoIt = true;
108  enablePostStepDoIt = true;
109  SetProcessSubType(60);
112  fInstantiateProcessState = true;
113 }
114 
115 
117  G4VITProcess(right),
118  InitProcessState(fTransportationState, fpState)
119 {
120  // Copy attributes
121  fVerboseLevel = right.fVerboseLevel ;
129 
130  // Setup Navigators
131  G4TransportationManager* transportMgr ;
133  G4ITTransportationManager* ITtransportMgr ;
135  fLinearNavigator = ITtransportMgr->GetNavigatorForTracking() ;
136  fFieldPropagator = transportMgr->GetPropagatorInField() ;
137  fpSafetyHelper = 0; //transportMgr->GetSafetyHelper(); // New
138 
139  // Cannot determine whether a field exists here, as it would
140  // depend on the relative order of creating the detector's
141  // field and this process. That order is not guaranted.
142  // Instead later the method DoesGlobalFieldExist() is called
143 
144  enableAtRestDoIt = false;
145  enableAlongStepDoIt = true;
146  enablePostStepDoIt = true;
147 
151  fInstantiateProcessState = right.fInstantiateProcessState;
152 }
153 
154 G4ITTransportation& G4ITTransportation::operator=(const G4ITTransportation& right)
155 {
156  if(this == &right) return *this;
157  return *this;
158 }
159 
164  fCurrentTouchableHandle(0)
165 {
170  fMomentumChanged = false;
171  fEnergyChanged = false;
172  fEndGlobalTimeComputed = false;
174  fParticleIsLooping = false;
175  static G4ThreadLocal G4TouchableHandle *nullTouchableHandle = 0;
176  if (!nullTouchableHandle) nullTouchableHandle = new G4TouchableHandle; // Points to (G4VTouchable*) 0
177  fCurrentTouchableHandle = *nullTouchableHandle;
178  fGeometryLimitedStep = false;
180  fPreviousSafety = 0.0;
181  fNoLooperTrials = false;
182  endpointDistance= -1;
183 }
184 
186 {
187  ;
188 }
189 
191 {
192 #ifdef G4VERBOSE
193  if( (fVerboseLevel > 0) && (fSumEnergyKilled > 0.0 ) )
194  {
195  G4cout << " G4ITTransportation: Statistics for looping particles " << G4endl;
196  G4cout << " Sum of energy of loopers killed: " << fSumEnergyKilled << G4endl;
197  G4cout << " Max energy of loopers killed: " << fMaxEnergyKilled << G4endl;
198  }
199 #endif
200 }
201 
203 {
204  G4TransportationManager* transportMgr;
206 
207  // fFieldExists= transportMgr->GetFieldManager()->DoesFieldExist();
208  // return fFieldExists;
209  return transportMgr->GetFieldManager()->DoesFieldExist();
210 }
211 
212 
214 //
215 // Responsibilities:
216 // Find whether the geometry limits the Step, and to what length
217 // Calculate the new value of the safety and return it.
218 // Store the final time, position and momentum.
220  const G4Track& track,
221  G4double ,
222  G4double currentMinimumStep,
223  G4double& currentSafety,
224  G4GPILSelection* selection)
225 {
226  G4double geometryStepLength(-1.0), newSafety(-1.0) ;
227 
228  State(fParticleIsLooping) = false ;
229  State(fEndGlobalTimeComputed) = false ;
230  State(fGeometryLimitedStep) = false ;
231 
232  // Initial actions moved to StartTrack()
233  // --------------------------------------
234  // Note: in case another process changes touchable handle
235  // it will be necessary to add here (for all steps)
236  State(fCurrentTouchableHandle) = track.GetTouchableHandle();
237 
238  // GPILSelection is set to defaule value of CandidateForSelection
239  // It is a return value
240  //
241  *selection = CandidateForSelection ;
242 
243  // Get initial Energy/Momentum of the track
244  //
245  const G4DynamicParticle* pParticle = track.GetDynamicParticle() ;
246 // const G4ParticleDefinition* pParticleDef = pParticle->GetDefinition() ;
247  G4ThreeVector startMomentumDir = pParticle->GetMomentumDirection() ;
248  G4ThreeVector startPosition = track.GetPosition() ;
249 
250  // G4double theTime = track.GetGlobalTime() ;
251 
252  // The Step Point safety can be limited by other geometries and/or the
253  // assumptions of any process - it's not always the geometrical safety.
254  // We calculate the starting point's isotropic safety here.
255  //
256  G4ThreeVector OriginShift = startPosition - State(fPreviousSftOrigin) ;
257  G4double MagSqShift = OriginShift.mag2() ;
258  if( MagSqShift >= sqr(State(fPreviousSafety)) )
259  {
260  currentSafety = 0.0 ;
261  }
262  else
263  {
264  currentSafety = State(fPreviousSafety) - std::sqrt(MagSqShift) ;
265  }
266 
267  // Is the particle charged ?
268  //
269  G4double particleCharge = pParticle->GetCharge() ;
270 
271 
272  // There is no need to locate the current volume. It is Done elsewhere:
273  // On track construction
274  // By the tracking, after all AlongStepDoIts, in "Relocation"
275 
276  // Check whether the particle have an (EM) field force exerting upon it
277  //
278  G4FieldManager* fieldMgr=0;
279  G4bool fieldExertsForce = false ;
280  if( (particleCharge != 0.0) )
281  {
282  fieldMgr= fFieldPropagator->FindAndSetFieldManager( track.GetVolume() );
283  if (fieldMgr != 0)
284  {
285  // Message the field Manager, to configure it for this track
286  fieldMgr->ConfigureForTrack( &track );
287  // Moved here, in order to allow a transition
288  // from a zero-field status (with fieldMgr->(field)0
289  // to a finite field status
290 
291  // If the field manager has no field, there is no field !
292  fieldExertsForce = (fieldMgr->GetDetectorField() != 0);
293  }
294  }
295 
296  // G4cout << " G4Transport: field exerts force= " << fieldExertsForce
297  // << " fieldMgr= " << fieldMgr << G4endl;
298 
299  // Choose the calculation of the transportation: Field or not
300  //
301  if( !fieldExertsForce )
302  {
303  G4double linearStepLength ;
304  if( fShortStepOptimisation && (currentMinimumStep <= currentSafety) )
305  {
306  // The Step is guaranteed to be taken
307  //
308  geometryStepLength = currentMinimumStep ;
309  State(fGeometryLimitedStep) = false ;
310  }
311  else
312  {
313  // Find whether the straight path intersects a volume
314  //
315  linearStepLength = fLinearNavigator->ComputeStep( startPosition,
316  startMomentumDir,
317  currentMinimumStep,
318  newSafety) ;
319  // Remember last safety origin & value.
320  //
321  State(fPreviousSftOrigin) = startPosition ;
322  State(fPreviousSafety) = newSafety ;
323  // fpSafetyHelper->SetCurrentSafety( newSafety, startPosition);
324 
325  // The safety at the initial point has been re-calculated:
326  //
327  currentSafety = newSafety ;
328 
329  State(fGeometryLimitedStep)= (linearStepLength <= currentMinimumStep);
330  if( State(fGeometryLimitedStep) )
331  {
332  // The geometry limits the Step size (an intersection was found.)
333  geometryStepLength = linearStepLength ;
334  }
335  else
336  {
337  // The full Step is taken.
338  geometryStepLength = currentMinimumStep ;
339  }
340  }
341  State(endpointDistance) = geometryStepLength ;
342 
343  // Calculate final position
344  //
345  State(fTransportEndPosition) = startPosition+geometryStepLength*startMomentumDir ;
346 
347  // Momentum direction, energy and polarisation are unchanged by transport
348  //
349  State(fTransportEndMomentumDir) = startMomentumDir ;
350  State(fTransportEndKineticEnergy) = track.GetKineticEnergy() ;
351  State(fTransportEndSpin) = track.GetPolarization();
352  State(fParticleIsLooping) = false ;
353  State(fMomentumChanged) = false ;
354  State(fEndGlobalTimeComputed) = true ;
355  State(theInteractionTimeLeft) = State(endpointDistance)/track.GetVelocity();
356  State(fCandidateEndGlobalTime) = State(theInteractionTimeLeft)+track.GetGlobalTime();
357 
358  // G4cout << "track.GetVelocity() : " << track.GetVelocity() << G4endl;
359  // G4cout << "State(endpointDistance) : " << G4BestUnit(State(endpointDistance),"Length") << G4endl;
360  // G4cout << "State(theInteractionTimeLeft) : " << G4BestUnit(State(theInteractionTimeLeft),"Time") << G4endl;
361  // G4cout << "track.GetGlobalTime() : " << G4BestUnit(track.GetGlobalTime(),"Time") << G4endl;
362 
363  }
364  else // A field exerts force
365  {
366 
367  G4ExceptionDescription exceptionDescription;
368  exceptionDescription << "ITTransportation does not support external fields.";
369  exceptionDescription << " If you are dealing with a tradiational MC simulation, ";
370  exceptionDescription << "please use G4Transportation.";
371 
372  G4Exception("G4ITTransportation::AlongStepGetPhysicalInteractionLength","NoExternalFieldSupport",
373  FatalException,exceptionDescription);
374  /*
375  G4double momentumMagnitude = pParticle->GetTotalMomentum() ;
376  // G4ThreeVector EndUnitMomentum ;
377  G4double lengthAlongCurve ;
378  G4double restMass = pParticleDef->GetPDGMass() ;
379 
380  fFieldPropagator->SetChargeMomentumMass( particleCharge, // in e+ units
381  momentumMagnitude, // in Mev/c
382  restMass ) ;
383 
384  G4ThreeVector spin = track.GetPolarization() ;
385  G4FieldTrack aFieldTrack = G4FieldTrack( startPosition,
386  track.GetMomentumDirection(),
387  0.0,
388  track.GetKineticEnergy(),
389  restMass,
390  track.GetVelocity(),
391  track.GetGlobalTime(), // Lab.
392  track.GetProperTime(), // Part.
393  &spin ) ;
394  if( currentMinimumStep > 0 )
395  {
396  // Do the Transport in the field (non recti-linear)
397  //
398  lengthAlongCurve = fFieldPropagator->ComputeStep( aFieldTrack,
399  currentMinimumStep,
400  currentSafety,
401  track.GetVolume() ) ;
402  State(fGeometryLimitedStep)= lengthAlongCurve < currentMinimumStep;
403  if( State(fGeometryLimitedStep) )
404  {
405  geometryStepLength = lengthAlongCurve ;
406  }
407  else
408  {
409  geometryStepLength = currentMinimumStep ;
410  }
411 
412  // Remember last safety origin & value.
413  //
414  State(fPreviousSftOrigin) = startPosition ;
415  State(fPreviousSafety) = currentSafety ;
416  fpSafetyHelper->SetCurrentSafety( newSafety, startPosition);
417  }
418  else
419  {
420  geometryStepLength = lengthAlongCurve= 0.0 ;
421  State(fGeometryLimitedStep) = false ;
422  }
423 
424  // Get the End-Position and End-Momentum (Dir-ection)
425  //
426  State(fTransportEndPosition) = aFieldTrack.GetPosition() ;
427 
428  // Momentum: Magnitude and direction can be changed too now ...
429  //
430  State(fMomentumChanged) = true ;
431  State(fTransportEndMomentumDir) = aFieldTrack.GetMomentumDir() ;
432 
433  State(fTransportEndKineticEnergy) = aFieldTrack.GetKineticEnergy() ;
434 
435  if( fFieldPropagator->GetCurrentFieldManager()->DoesFieldChangeEnergy() )
436  {
437  // If the field can change energy, then the time must be integrated
438  // - so this should have been updated
439  //
440  State(fCandidateEndGlobalTime) = aFieldTrack.GetLabTimeOfFlight();
441  State(fEndGlobalTimeComputed) = true;
442 
443  State(theInteractionTimeLeft) = State(fCandidateEndGlobalTime) - track.GetGlobalTime() ;
444 
445  // was ( State(fCandidateEndGlobalTime) != track.GetGlobalTime() );
446  // a cleaner way is to have FieldTrack knowing whether time is updated.
447  }
448  else
449  {
450  // The energy should be unchanged by field transport,
451  // - so the time changed will be calculated elsewhere
452  //
453  State(fEndGlobalTimeComputed) = false;
454 
455  // Check that the integration preserved the energy
456  // - and if not correct this!
457  G4double startEnergy= track.GetKineticEnergy();
458  G4double endEnergy= State(fTransportEndKineticEnergy);
459 
460  static G4int no_inexact_steps=0, no_large_ediff;
461  G4double absEdiff = std::fabs(startEnergy- endEnergy);
462  if( absEdiff > perMillion * endEnergy )
463  {
464  no_inexact_steps++;
465  // Possible statistics keeping here ...
466  }
467 #ifdef G4VERBOSE
468  if( fVerboseLevel > 1 )
469  {
470  if( std::fabs(startEnergy- endEnergy) > perThousand * endEnergy )
471  {
472  static G4int no_warnings= 0, warnModulo=1, moduloFactor= 10;
473  no_large_ediff ++;
474  if( (no_large_ediff% warnModulo) == 0 )
475  {
476  no_warnings++;
477  G4cout << "WARNING - G4Transportation::AlongStepGetPIL() "
478  << " Energy change in Step is above 1^-3 relative value. " << G4endl
479  << " Relative change in 'tracking' step = "
480  << std::setw(15) << (endEnergy-startEnergy)/startEnergy << G4endl
481  << " Starting E= " << std::setw(12) << startEnergy / MeV << " MeV " << G4endl
482  << " Ending E= " << std::setw(12) << endEnergy / MeV << " MeV " << G4endl;
483  G4cout << " Energy has been corrected -- however, review"
484  << " field propagation parameters for accuracy." << G4endl;
485  if( (fVerboseLevel > 2 ) || (no_warnings<4) || (no_large_ediff == warnModulo * moduloFactor) )
486  {
487  G4cout << " These include EpsilonStepMax(/Min) in G4FieldManager "
488  << " which determine fractional error per step for integrated quantities. " << G4endl
489  << " Note also the influence of the permitted number of integration steps."
490  << G4endl;
491  }
492  G4cerr << "ERROR - G4Transportation::AlongStepGetPIL()" << G4endl
493  << " Bad 'endpoint'. Energy change detected"
494  << " and corrected. "
495  << " Has occurred already "
496  << no_large_ediff << " times." << G4endl;
497  if( no_large_ediff == warnModulo * moduloFactor )
498  {
499  warnModulo *= moduloFactor;
500  }
501  }
502  }
503  } // end of if (fVerboseLevel)
504 #endif
505  // Correct the energy for fields that conserve it
506  // This - hides the integration error
507  // - but gives a better physical answer
508  State(fTransportEndKineticEnergy)= track.GetKineticEnergy();
509  }
510 
511  State(fTransportEndSpin) = aFieldTrack.GetSpin();
512  State(fParticleIsLooping) = fFieldPropagator->IsParticleLooping() ;
513  State(endpointDistance) = (State(fTransportEndPosition) - startPosition).mag() ;
514  // State(theInteractionTimeLeft) = track.GetVelocity()/State(endpointDistance) ;
515 */
516  }
517 
518  // If we are asked to go a step length of 0, and we are on a boundary
519  // then a boundary will also limit the step -> we must flag this.
520  //
521  if( currentMinimumStep == 0.0 )
522  {
523  if( currentSafety == 0.0 )
524  {
525  State(fGeometryLimitedStep) = true ;
526 // G4cout << "!!!! Safety is NULL, on the Boundary !!!!!" << G4endl;
527 // G4cout << " Track position : " << track.GetPosition() /nanometer << G4endl;
528  }
529  }
530 
531  // Update the safety starting from the end-point,
532  // if it will become negative at the end-point.
533  //
534  if( currentSafety < State(endpointDistance) )
535  {
536  // if( particleCharge == 0.0 )
537  // G4cout << " Avoiding call to ComputeSafety : charge = 0.0 " << G4endl;
538 
539  if( particleCharge != 0.0 )
540  {
541 
542  G4double endSafety =
543  fLinearNavigator->ComputeSafety( State(fTransportEndPosition)) ;
544  currentSafety = endSafety ;
545  State(fPreviousSftOrigin) = State(fTransportEndPosition) ;
546  State(fPreviousSafety) = currentSafety ;
547  fpSafetyHelper->SetCurrentSafety( currentSafety, State(fTransportEndPosition));
548 
549  // Because the Stepping Manager assumes it is from the start point,
550  // add the StepLength
551  //
552  currentSafety += State(endpointDistance) ;
553 
554 #ifdef G4DEBUG_TRANSPORT
555  G4cout.precision(12) ;
556  G4cout << "***G4Transportation::AlongStepGPIL ** " << G4endl ;
557  G4cout << " Called Navigator->ComputeSafety at " << State(fTransportEndPosition)
558  << " and it returned safety= " << endSafety << G4endl ;
559  G4cout << " Adding endpoint distance " << State(endpointDistance)
560  << " to obtain pseudo-safety= " << currentSafety << G4endl ;
561 #endif
562  }
563  }
564 
565  // fParticleChange.ProposeTrueStepLength(geometryStepLength) ;
566 
567  return geometryStepLength ;
568 }
569 
570 
572  const G4Step& /*step*/,
573  const double timeStep,
574  double& oPhysicalStep)
575 {
576  const G4DynamicParticle* pParticle = track.GetDynamicParticle() ;
577  G4ThreeVector startMomentumDir = pParticle->GetMomentumDirection() ;
578  G4ThreeVector startPosition = track.GetPosition() ;
579 
580  track.CalculateVelocity();
581  G4double initialVelocity = track.GetVelocity() ;
582 
583  State(fGeometryLimitedStep) = false;
584 
586  // !!! CASE NO FIELD !!!
588  State(fCandidateEndGlobalTime) = timeStep + track.GetGlobalTime();
589  State(fEndGlobalTimeComputed) = true ;
590 
591  // Choose the calculation of the transportation: Field or not
592  //
593  if( !State(fMomentumChanged) )
594  {
595  // G4cout << "Momentum has not changed" << G4endl;
596  fParticleChange.ProposeVelocity(initialVelocity);
597  oPhysicalStep = initialVelocity*timeStep ;
598 
599  // Calculate final position
600  //
601  State(fTransportEndPosition) = startPosition + oPhysicalStep*startMomentumDir ;
602  }
603 }
604 
605 
607 //
608 // Initialize ParticleChange (by setting all its members equal
609 // to corresponding members in G4Track)
610 #include "G4ParticleTable.hh"
612  const G4Step& stepData )
613 {
614 
615  // G4cout << "G4ITTransportation::AlongStepDoIt" << G4endl;
616  // set pdefOpticalPhoton
617  //Andrea Dotti: the following statement should be in a single line:
618  // G4-MT transformation tools get confused if statement spans two lines
619  // If needed contact: adotti@slac.stanford.edu
620  static G4ThreadLocal G4ParticleDefinition* pdefOpticalPhoton = 0 ; if (!pdefOpticalPhoton) pdefOpticalPhoton= G4ParticleTable::GetParticleTable()->FindParticle("opticalphoton");
621 
622  static G4ThreadLocal G4int noCalls=0;
623  noCalls++;
624 
625  fParticleChange.Initialize(track) ;
626 
627  // Code for specific process
628  //
629  fParticleChange.ProposePosition(State(fTransportEndPosition)) ;
630  fParticleChange.ProposeMomentumDirection(State(fTransportEndMomentumDir)) ;
631  fParticleChange.ProposeEnergy(State(fTransportEndKineticEnergy)) ;
632  fParticleChange.SetMomentumChanged(State(fMomentumChanged)) ;
633 
634  fParticleChange.ProposePolarization(State(fTransportEndSpin));
635 
636  G4double deltaTime = 0.0 ;
637 
638  // Calculate Lab Time of Flight (ONLY if field Equations used it!)
639  // G4double endTime = State(fCandidateEndGlobalTime);
640  // G4double delta_time = endTime - startTime;
641 
642  G4double startTime = track.GetGlobalTime() ;
646  if (State(fEndGlobalTimeComputed) == false)
647  {
648  // The time was not integrated .. make the best estimate possible
649  //
650  G4double initialVelocity = stepData.GetPreStepPoint()->GetVelocity() ;
651  G4double stepLength = track.GetStepLength() ;
652 
653  deltaTime= 0.0; // in case initialVelocity = 0
654  if (track.GetParticleDefinition() == pdefOpticalPhoton)
655  {
656  // For only Optical Photon, final velocity is used
657  double finalVelocity = track.CalculateVelocityForOpticalPhoton();
658  fParticleChange.ProposeVelocity(finalVelocity);
659  deltaTime = stepLength/finalVelocity ;
660  }
661  else if( initialVelocity > 0.0 )
662  {
663  deltaTime = stepLength/initialVelocity ;
664  }
665 
666  State(fCandidateEndGlobalTime) = startTime + deltaTime ;
667  }
668  else
669  {
670  deltaTime = State(fCandidateEndGlobalTime) - startTime ;
671  }
672 
673  fParticleChange.ProposeGlobalTime( State(fCandidateEndGlobalTime) ) ;
674  fParticleChange.ProposeLocalTime( track.GetLocalTime() + deltaTime) ;
675  /*
676  // Now Correct by Lorentz factor to get delta "proper" Time
677 
678  G4double restMass = track.GetDynamicParticle()->GetMass() ;
679  G4double deltaProperTime = deltaTime*( restMass/track.GetTotalEnergy() ) ;
680 
681  fParticleChange.ProposeProperTime(track.GetProperTime() + deltaProperTime) ;
682 */
683 
684  fParticleChange. ProposeTrueStepLength( track.GetStepLength() ) ;
685 
688 
689  // If the particle is caught looping or is stuck (in very difficult
690  // boundaries) in a magnetic field (doing many steps)
691  // THEN this kills it ...
692  //
693  if ( State(fParticleIsLooping) )
694  {
695  G4double endEnergy= State(fTransportEndKineticEnergy);
696 
697  if( (endEnergy < fThreshold_Important_Energy)
698  || (State(fNoLooperTrials) >= fThresholdTrials ) )
699  {
700  // Kill the looping particle
701  //
702  // G4cout << "G4ITTransportation will killed the molecule"<< G4endl;
704 
705  // 'Bare' statistics
706  fSumEnergyKilled += endEnergy;
707  if( endEnergy > fMaxEnergyKilled)
708  {
709  fMaxEnergyKilled= endEnergy;
710  }
711 
712 #ifdef G4VERBOSE
713  if( (fVerboseLevel > 1) ||
714  ( endEnergy > fThreshold_Warning_Energy ) )
715  {
716  G4cout << " G4ITTransportation is killing track that is looping or stuck "
717  << G4endl
718  << " This track has " << track.GetKineticEnergy() / MeV
719  << " MeV energy." << G4endl;
720  G4cout << " Number of trials = " << State(fNoLooperTrials)
721  << " No of calls to AlongStepDoIt = " << noCalls
722  << G4endl;
723  }
724 #endif
725  State(fNoLooperTrials)=0;
726  }
727  else
728  {
729  State(fNoLooperTrials) ++;
730 #ifdef G4VERBOSE
731  if( (fVerboseLevel > 2) )
732  {
733  G4cout << " G4ITTransportation::AlongStepDoIt(): Particle looping - "
734  << " Number of trials = " << State(fNoLooperTrials)
735  << " No of calls to = " << noCalls
736  << G4endl;
737  }
738 #endif
739  }
740  }
741  else
742  {
743  State(fNoLooperTrials)=0;
744  }
745 
746  // Another (sometimes better way) is to use a user-limit maximum Step size
747  // to alleviate this problem ..
748 
749  // Introduce smooth curved trajectories to particle-change
750  //
753 
754  return &fParticleChange ;
755 }
756 
758 //
759 // This ensures that the PostStep action is always called,
760 // so that it can do the relocation if it is needed.
761 //
762 
765  const G4Track& , // track
766  G4double , // previousStepSize
767  G4ForceCondition* pForceCond
768  )
769 {
770  *pForceCond = Forced ;
771  return DBL_MAX ; // was kInfinity ; but convention now is DBL_MAX
772 }
773 
775 //
776 
778  const G4Step& )
779 {
780  // G4cout << "G4ITTransportation::PostStepDoIt" << G4endl;
781  G4TouchableHandle retCurrentTouchable ; // The one to return
782  G4bool isLastStep= false;
783 
784  // Initialize ParticleChange (by setting all its members equal
785  // to corresponding members in G4Track)
786  fParticleChange.Initialize(track) ; // To initialise TouchableChange
787 
789 
790  // If the Step was determined by the volume boundary,
791  // logically relocate the particle
792 
793  if(State(fGeometryLimitedStep))
794  {
795 
796  // G4cout << "Step is limited by geometry " << "track ID : " << track.GetTrackID() << G4endl;
797 
798  // fCurrentTouchable will now become the previous touchable,
799  // and what was the previous will be freed.
800  // (Needed because the preStepPoint can point to the previous touchable)
801 
802  if( State(fCurrentTouchableHandle)->GetVolume() == 0 )
803  {
804  G4ExceptionDescription exceptionDescription ;
805  exceptionDescription << "No current touchable found " ;
806  G4Exception(" G4ITTransportation::PostStepDoIt","G4ITTransportation001",
807  FatalErrorInArgument,exceptionDescription);
808  }
809 
812  LocateGlobalPointAndUpdateTouchableHandle( track.GetPosition(),
813  track.GetMomentumDirection(),
814  State(fCurrentTouchableHandle),
815  true ) ;
816  // Check whether the particle is out of the world volume
817  // If so it has exited and must be killed.
818  //
819  if( State(fCurrentTouchableHandle)->GetVolume() == 0 )
820  {
821 #ifdef G4VERBOSE
822  if(fVerboseLevel > 0)
823  {
824  G4cout << "Track position : " << track.GetPosition() / nanometer << " [nm]"
825  << " Track ID : " << track.GetTrackID()<< G4endl;
826  G4cout << "G4ITTransportation will killed the track because State(fCurrentTouchableHandle)->GetVolume() == 0"<< G4endl;
827  }
828 #endif
830  }
831 
832  retCurrentTouchable = State(fCurrentTouchableHandle) ;
833 
834 // G4cout << "Current volume : " << track.GetVolume()->GetName()
835 // << " Next volume : "
836 // << (State(fCurrentTouchableHandle)->GetVolume() ? State(fCurrentTouchableHandle)->GetVolume()->GetName():"OutWorld")
837 // << " Position : " << track.GetPosition() / nanometer
838 // << " track ID : " << track.GetTrackID()
839 // << G4endl;
840 
841  fParticleChange.SetTouchableHandle( State(fCurrentTouchableHandle) ) ;
842 
843  // Update the Step flag which identifies the Last Step in a volume
844  isLastStep = fLinearNavigator->ExitedMotherVolume()
846 
847 #ifdef G4DEBUG_TRANSPORT
848  // Checking first implementation of flagging Last Step in Volume
851 
852  if( ! (exiting || entering) )
853  {
854  G4cout << " Transport> : Proposed isLastStep= " << isLastStep
855  << " Exiting " << fLinearNavigator->ExitedMotherVolume()
856  << " Entering " << fLinearNavigator->EnteredDaughterVolume()
857  << " Track position : " << track.GetPosition() /nanometer << " [nm]"
858  << G4endl;
859  G4cout << " Track position : " << track.GetPosition() /nanometer << G4endl;
860  }
861 #endif
862  }
863  else // fGeometryLimitedStep is false
864  {
865  // This serves only to move the Navigator's location
866  //
868 
869  // The value of the track's current Touchable is retained.
870  // (and it must be correct because we must use it below to
871  // overwrite the (unset) one in particle change)
872  // It must be fCurrentTouchable too ??
873  //
875  retCurrentTouchable = track.GetTouchableHandle() ;
876 
877  isLastStep= false;
878 #ifdef G4DEBUG_TRANSPORT
879  // Checking first implementation of flagging Last Step in Volume
880  G4cout << " Transport> Proposed isLastStep= " << isLastStep
881  << " Geometry did not limit step. Position : "
882  << track.GetPosition()/ nanometer << G4endl;
883 #endif
884  } // endif ( fGeometryLimitedStep )
885 
887 
888  const G4VPhysicalVolume* pNewVol = retCurrentTouchable->GetVolume() ;
889  const G4Material* pNewMaterial = 0 ;
890  const G4VSensitiveDetector* pNewSensitiveDetector = 0 ;
891 
892  if( pNewVol != 0 )
893  {
894  pNewMaterial= pNewVol->GetLogicalVolume()->GetMaterial();
895  pNewSensitiveDetector= pNewVol->GetLogicalVolume()->GetSensitiveDetector();
896  }
897 
898  // ( <const_cast> pNewMaterial ) ;
899  // ( <const_cast> pNewSensitiveDetector) ;
900 
903 
904  const G4MaterialCutsCouple* pNewMaterialCutsCouple = 0;
905  if( pNewVol != 0 )
906  {
907  pNewMaterialCutsCouple=pNewVol->GetLogicalVolume()->GetMaterialCutsCouple();
908  }
909 
910  if( pNewVol!=0 && pNewMaterialCutsCouple!=0 && pNewMaterialCutsCouple->GetMaterial()!=pNewMaterial )
911  {
912  // for parametrized volume
913  //
914  pNewMaterialCutsCouple =
916  ->GetMaterialCutsCouple(pNewMaterial,
917  pNewMaterialCutsCouple->GetProductionCuts());
918  }
919  fParticleChange.SetMaterialCutsCoupleInTouchable( pNewMaterialCutsCouple );
920 
921  // temporarily until Get/Set Material of ParticleChange,
922  // and StepPoint can be made const.
923  // Set the touchable in ParticleChange
924  // this must always be done because the particle change always
925  // uses this value to overwrite the current touchable pointer.
926  //
927  fParticleChange.SetTouchableHandle(retCurrentTouchable) ;
928 
929  return &fParticleChange ;
930 }
931 
932 // New method takes over the responsibility to reset the state of G4Transportation
933 // object at the start of a new track or the resumption of a suspended track.
934 
935 void
937 {
939  if(fInstantiateProcessState)
941  // Will set in the same time fTransportationState
942 
943  // The actions here are those that were taken in AlongStepGPIL
944  // when track.GetCurrentStepNumber()==1
945 
946  // reset safety value and center
947  //
948  // State(fPreviousSafety) = 0.0 ;
949  // State(fPreviousSftOrigin) = G4ThreeVector(0.,0.,0.) ;
950 
951  // reset looping counter -- for motion in field
952  // State(fNoLooperTrials)= 0;
953  // Must clear this state .. else it depends on last track's value
954  // --> a better solution would set this from state of suspended track TODO ?
955  // Was if( aTrack->GetCurrentStepNumber()==1 ) { .. }
956 
957  // ChordFinder reset internal state
958  //
959  if( DoesGlobalFieldExist() )
960  {
962  // Resets all state of field propagator class (ONLY)
963  // including safety values (in case of overlaps and to wipe for first track).
964 
965  // G4ChordFinder* chordF= fFieldPropagator->GetChordFinder();
966  // if( chordF ) chordF->ResetStepEstimate();
967  }
968 
969  // Make sure to clear the chord finders of all fields (ie managers)
970  static G4ThreadLocal G4FieldManagerStore* fieldMgrStore = 0 ; if (!fieldMgrStore) fieldMgrStore= G4FieldManagerStore::GetInstance();
971  fieldMgrStore->ClearAllChordFindersState();
972 
973  // Update the current touchable handle (from the track's)
974  //
975  State(fCurrentTouchableHandle) = track->GetTouchableHandle();
976 
978 }
979 
980 #undef State
void SetMaterialInTouchable(G4Material *fMaterial)
void SetTouchableHandle(const G4TouchableHandle &fTouchable)
const G4ThreeVector & GetPolarization() const
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
G4double GetLocalTime() const
G4ParticleChangeForTransport fParticleChange
void SetMaterialCutsCoupleInTouchable(const G4MaterialCutsCouple *fMaterialCutsCouple)
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
std::vector< G4ThreeVector > * GimmeTrajectoryVectorAndForgetIt() const
CLHEP::Hep3Vector G4ThreeVector
G4double GetVelocity() const
const G4DynamicParticle * GetDynamicParticle() const
G4Material * GetMaterial() const
virtual G4double ComputeSafety(const G4ThreeVector &globalpoint, const G4double pProposedMaxLength=DBL_MAX, const G4bool keepState=false)
const G4ThreeVector & GetPosition() const
virtual G4double ComputeStep(const G4ThreeVector &pGlobalPoint, const G4ThreeVector &pDirection, const G4double pCurrentProposedStepLength, G4double &pNewSafety)
G4TrackStatus GetTrackStatus() const
virtual G4VParticleChange * PostStepDoIt(const G4Track &track, const G4Step &)
void ProposePolarization(G4double Px, G4double Py, G4double Pz)
virtual G4VParticleChange * AlongStepDoIt(const G4Track &track, const G4Step &stepData)
G4ProcessState * fpState
virtual void LocateGlobalPointWithinVolume(const G4ThreeVector &position)
G4ITTransportation(const G4String &aName="ITTransportation", G4int verbosityLevel=1)
G4double GetVelocity() const
virtual G4double AlongStepGetPhysicalInteractionLength(const G4Track &track, G4double, G4double currentMinimumStep, G4double &currentSafety, G4GPILSelection *selection)
#define G4ThreadLocal
Definition: tls.hh:52
G4bool enableAtRestDoIt
Definition: G4VProcess.hh:350
#define State(theXInfo)
int G4int
Definition: G4Types.hh:78
G4bool ExitedMotherVolume() const
void SetInstantiateProcessState(G4bool flag)
void ProposePosition(G4double x, G4double y, G4double z)
G4ReferenceCountedHandle< G4VTouchable > G4TouchableHandle
G4ITNavigator * fLinearNavigator
int nanometer
Definition: hepunit.py:35
G4StepPoint * GetPreStepPoint() const
virtual void Initialize(const G4Track &)
virtual void StartTracking(G4Track *)
Definition: G4VProcess.cc:101
G4double CalculateVelocityForOpticalPhoton() const
Definition: G4Track.cc:247
virtual void ConfigureForTrack(const G4Track *)
virtual void StartTracking(G4Track *aTrack)
void SetGeometricallyLimitedStep()
G4double GetKineticEnergy() const
static G4ITTransportationManager * GetTransportationManager()
virtual G4double PostStepGetPhysicalInteractionLength(const G4Track &, G4double, G4ForceCondition *pForceCond)
G4GLOB_DLL std::ostream G4cout
virtual void ComputeStep(const G4Track &, const G4Step &, const double timeStep, double &spaceStep)
bool G4bool
Definition: G4Types.hh:79
const G4ThreeVector & GetMomentumDirection() const
G4PropagatorInField * fFieldPropagator
G4double GetCharge() const
G4bool DoesFieldExist() const
virtual void StartTracking(G4Track *)
Definition: G4VITProcess.cc:83
void SetInstantiateProcessState(G4bool flag)
const G4ParticleDefinition * GetParticleDefinition() const
void SetProcessSubType(G4int)
Definition: G4VProcess.hh:432
Definition: G4Step.hh:76
G4bool enablePostStepDoIt
Definition: G4VProcess.hh:352
G4int GetTrackID() const
G4double GetGlobalTime() const
G4double CalculateVelocity() const
Definition: G4Track.cc:214
G4FieldManager * FindAndSetFieldManager(G4VPhysicalVolume *pCurrentPhysVol)
const G4TouchableHandle & GetTouchableHandle() const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
static G4TransportationManager * GetTransportationManager()
G4FieldManager * GetFieldManager() const
static G4ProductionCutsTable * GetProductionCutsTable()
const G4ThreeVector & GetMomentumDirection() const
#define InitProcessState(destination, source)
Definition: G4VITProcess.hh:53
G4LogicalVolume * GetLogicalVolume() const
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
static G4FieldManagerStore * GetInstance()
void SetSensitiveDetectorInTouchable(G4VSensitiveDetector *fSensitiveDetector)
static G4ParticleTable * GetParticleTable()
virtual G4VPhysicalVolume * GetVolume(G4int depth=0) const
Definition: G4VTouchable.cc:44
void SetPointerToVectorOfAuxiliaryPoints(std::vector< G4ThreeVector > *theNewVectorPointer)
void ProposeGlobalTime(G4double t)
G4bool EnteredDaughterVolume() const
G4VParticleChange * pParticleChange
Definition: G4VProcess.hh:283
void ProposeEnergy(G4double finalEnergy)
G4SafetyHelper * fpSafetyHelper
double mag2() const
G4VPhysicalVolume * GetVolume() const
void SetCurrentSafety(G4double val, const G4ThreeVector &pos)
#define G4endl
Definition: G4ios.hh:61
G4bool enableAlongStepDoIt
Definition: G4VProcess.hh:351
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)
G4double fThreshold_Important_Energy
const G4Field * GetDetectorField() const
G4ForceCondition
void SetMomentumChanged(G4bool b)
G4PropagatorInField * GetPropagatorInField() const
void ProposeVelocity(G4double finalVelocity)
G4ProductionCuts * GetProductionCuts() const
#define DBL_MAX
Definition: templates.hh:83
G4VSensitiveDetector * GetSensitiveDetector() const
void ProposeLocalTime(G4double t)
G4double GetStepLength() const
G4GPILSelection
const G4Material * GetMaterial() const