Geant4  10.03.p03
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 94289 2015-11-11 08:33:40Z 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 "G4IT.hh"
60 #include "G4TrackingInformation.hh"
61 #include "G4SystemOfUnits.hh"
64 #include "G4ProductionCutsTable.hh"
65 #include "G4ParticleTable.hh"
66 #include "G4ITNavigator.hh"
67 #include "G4PropagatorInField.hh"
68 #include "G4FieldManager.hh"
69 #include "G4ChordFinder.hh"
70 #include "G4ITSafetyHelper.hh"
71 #include "G4FieldManagerStore.hh"
72 
73 #include "G4UnitsTable.hh"
74 #include "G4ReferenceCast.hh"
75 
77 
78 #ifndef PrepareState
79 #define PrepareState() G4ITTransportationState* __state = this->GetState<G4ITTransportationState>();
80 #endif
81 
82 #ifndef State
83 #define State(theXInfo) (__state->theXInfo)
84 #endif
85 
86 //#define DEBUG_MEM
87 
88 #ifdef DEBUG_MEM
89 #include "G4MemStat.hh"
90 using namespace G4MemStat;
91 using G4MemStat::MemStat;
92 #endif
93 
94 //#define G4DEBUG_TRANSPORT 1
95 
98  fThreshold_Warning_Energy(100 * MeV),
99  fThreshold_Important_Energy(250 * MeV),
100  fThresholdTrials(10),
101  fUnimportant_Energy(1 * MeV), // Not used
102  fSumEnergyKilled(0.0),
103  fMaxEnergyKilled(0.0),
104  fShortStepOptimisation(false), // Old default: true (=fast short steps)
105  fVerboseLevel(verbose)
106 {
108  G4TransportationManager* transportMgr;
110  G4ITTransportationManager* ITtransportMgr;
112  fLinearNavigator = ITtransportMgr->GetNavigatorForTracking();
113  fFieldPropagator = transportMgr->GetPropagatorInField();
114  fpSafetyHelper = ITtransportMgr->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  enableAtRestDoIt = false;
122  enableAlongStepDoIt = true;
123  enablePostStepDoIt = true;
124  SetProcessSubType(60);
127  fInstantiateProcessState = true;
128 
130  /*
131  if(fTransportationState == 0)
132  {
133  G4cout << "KILL in G4ITTransportation::G4ITTransportation" << G4endl;
134  abort();
135  }
136  */
137 }
138 
140  G4VITProcess(right)
141 {
142  // Copy attributes
151 
152  // Setup Navigators
153  G4TransportationManager* transportMgr;
155  G4ITTransportationManager* ITtransportMgr;
157  fLinearNavigator = ITtransportMgr->GetNavigatorForTracking();
158  fFieldPropagator = transportMgr->GetPropagatorInField();
159  fpSafetyHelper = ITtransportMgr->GetSafetyHelper(); // New
160 
161  // Cannot determine whether a field exists here, as it would
162  // depend on the relative order of creating the detector's
163  // field and this process. That order is not guaranted.
164  // Instead later the method DoesGlobalFieldExist() is called
165 
166  enableAtRestDoIt = false;
167  enableAlongStepDoIt = true;
168  enablePostStepDoIt = true;
169 
173  fInstantiateProcessState = right.fInstantiateProcessState;
174 }
175 
176 G4ITTransportation& G4ITTransportation::operator=(const G4ITTransportation& /*right*/)
177 {
178 // if (this == &right) return *this;
179  return *this;
180 }
181 
186  G4ProcessState(), fCurrentTouchableHandle(0)
187 {
191  fTransportEndSpin = G4ThreeVector(0, 0, 0);
192  fMomentumChanged = false;
193  fEnergyChanged = false;
194  fEndGlobalTimeComputed = false;
196  fParticleIsLooping = false;
197 
198  static G4ThreadLocal G4TouchableHandle *nullTouchableHandle = 0;
199  if (!nullTouchableHandle) nullTouchableHandle = new G4TouchableHandle;
200  // Points to (G4VTouchable*) 0
201 
202  fCurrentTouchableHandle = *nullTouchableHandle;
203  fGeometryLimitedStep = false;
205  fPreviousSafety = 0.0;
206  fNoLooperTrials = false;
207  fEndPointDistance = -1;
208 }
209 
211 {
212  ;
213 }
214 
216 {
217 #ifdef G4VERBOSE
218  if ((fVerboseLevel > 0) && (fSumEnergyKilled > 0.0))
219  {
220  G4cout << " G4ITTransportation: Statistics for looping particles "
221  << G4endl;
222  G4cout << " Sum of energy of loopers killed: "
223  << fSumEnergyKilled << G4endl;
224  G4cout << " Max energy of loopers killed: "
225  << fMaxEnergyKilled << G4endl;
226  }
227 #endif
228 }
229 
231 {
233 }
234 
236 {
237  G4TransportationManager* transportMgr;
239 
240  // fFieldExists= transportMgr->GetFieldManager()->DoesFieldExist();
241  // return fFieldExists;
242  return transportMgr->GetFieldManager()->DoesFieldExist();
243 }
244 
246 //
247 // Responsibilities:
248 // Find whether the geometry limits the Step, and to what length
249 // Calculate the new value of the safety and return it.
250 // Store the final time, position and momentum.
251 G4double
254  G4double,
255  G4double currentMinimumStep,
256  G4double& currentSafety,
257  G4GPILSelection* selection)
258 {
259  PrepareState()
260  G4double geometryStepLength(-1.0), newSafety(-1.0);
261 
262  State(fParticleIsLooping) = false;
263  State(fEndGlobalTimeComputed) = false;
264  State(fGeometryLimitedStep) = false;
265 
266  // Initial actions moved to StartTrack()
267  // --------------------------------------
268  // Note: in case another process changes touchable handle
269  // it will be necessary to add here (for all steps)
270  // State(fCurrentTouchableHandle) = track.GetTouchableHandle();
271 
272  // GPILSelection is set to defaule value of CandidateForSelection
273  // It is a return value
274  //
275  *selection = CandidateForSelection;
276 
277  // Get initial Energy/Momentum of the track
278  //
279  const G4DynamicParticle* pParticle = track.GetDynamicParticle();
280 // const G4ParticleDefinition* pParticleDef = pParticle->GetDefinition() ;
281  G4ThreeVector startMomentumDir = pParticle->GetMomentumDirection();
282  G4ThreeVector startPosition = track.GetPosition();
283 
284  // G4double theTime = track.GetGlobalTime() ;
285 
286  // The Step Point safety can be limited by other geometries and/or the
287  // assumptions of any process - it's not always the geometrical safety.
288  // We calculate the starting point's isotropic safety here.
289  //
290  G4ThreeVector OriginShift = startPosition - State(fPreviousSftOrigin);
291  G4double MagSqShift = OriginShift.mag2();
292  if (MagSqShift >= sqr(State(fPreviousSafety)))
293  {
294  currentSafety = 0.0;
295  }
296  else
297  {
298  currentSafety = State(fPreviousSafety) - std::sqrt(MagSqShift);
299  }
300 
301  // Is the particle charged ?
302  //
303  G4double particleCharge = pParticle->GetCharge();
304 
305  // There is no need to locate the current volume. It is Done elsewhere:
306  // On track construction
307  // By the tracking, after all AlongStepDoIts, in "Relocation"
308 
309  // Check whether the particle have an (EM) field force exerting upon it
310  //
311  G4FieldManager* fieldMgr = 0;
312  G4bool fieldExertsForce = false;
313  if ((particleCharge != 0.0))
314  {
316  if (fieldMgr != 0)
317  {
318  // Message the field Manager, to configure it for this track
319  fieldMgr->ConfigureForTrack(&track);
320  // Moved here, in order to allow a transition
321  // from a zero-field status (with fieldMgr->(field)0
322  // to a finite field status
323 
324  // If the field manager has no field, there is no field !
325  fieldExertsForce = (fieldMgr->GetDetectorField() != 0);
326  }
327  }
328 
329  // G4cout << " G4Transport: field exerts force= " << fieldExertsForce
330  // << " fieldMgr= " << fieldMgr << G4endl;
331 
332  // Choose the calculation of the transportation: Field or not
333  //
334  if (!fieldExertsForce)
335  {
336  G4double linearStepLength;
337  if (fShortStepOptimisation && (currentMinimumStep <= currentSafety))
338  {
339  // The Step is guaranteed to be taken
340  //
341  geometryStepLength = currentMinimumStep;
342  State(fGeometryLimitedStep) = false;
343  }
344  else
345  {
346  // Find whether the straight path intersects a volume
347  //
348  // fLinearNavigator->SetNavigatorState(GetIT(track)->GetTrackingInfo()->GetNavigatorState());
349  linearStepLength = fLinearNavigator->ComputeStep(startPosition,
350  startMomentumDir,
351  currentMinimumStep,
352  newSafety);
353 
354 // G4cout << "linearStepLength : " << G4BestUnit(linearStepLength,"Length")
355 // << " | currentMinimumStep: " << currentMinimumStep
356 // << " | trackID: " << track.GetTrackID() << G4endl;
357 
358  // Remember last safety origin & value.
359  //
360  State(fPreviousSftOrigin) = startPosition;
361  State(fPreviousSafety) = newSafety;
362 
363  G4TrackStateManager& trackStateMan = GetIT(track)->GetTrackingInfo()
365  fpSafetyHelper->LoadTrackState(trackStateMan);
366  // fpSafetyHelper->SetTrackState(state);
367  fpSafetyHelper->SetCurrentSafety(newSafety,
368  State(fTransportEndPosition));
370 
371  // The safety at the initial point has been re-calculated:
372  //
373  currentSafety = newSafety;
374 
375  State(fGeometryLimitedStep) = (linearStepLength <= currentMinimumStep);
376  if (State(fGeometryLimitedStep))
377  {
378  // The geometry limits the Step size (an intersection was found.)
379  geometryStepLength = linearStepLength;
380  }
381  else
382  {
383  // The full Step is taken.
384  geometryStepLength = currentMinimumStep;
385  }
386  }
387  State(fEndPointDistance) = geometryStepLength;
388 
389  // Calculate final position
390  //
391  State(fTransportEndPosition) = startPosition
392  + geometryStepLength * startMomentumDir;
393 
394  // Momentum direction, energy and polarisation are unchanged by transport
395  //
396  State(fTransportEndMomentumDir) = startMomentumDir;
397  State(fTransportEndKineticEnergy) = track.GetKineticEnergy();
398  State(fTransportEndSpin) = track.GetPolarization();
399  State(fParticleIsLooping) = false;
400  State(fMomentumChanged) = false;
401  State(fEndGlobalTimeComputed) = true;
402  State(theInteractionTimeLeft) = State(fEndPointDistance)
403  / track.GetVelocity();
404  State(fCandidateEndGlobalTime) = State(theInteractionTimeLeft)
405  + track.GetGlobalTime();
406 /*
407  G4cout << "track.GetVelocity() : "
408  << track.GetVelocity() << G4endl;
409  G4cout << "State(endpointDistance) : "
410  << G4BestUnit(State(endpointDistance),"Length") << G4endl;
411  G4cout << "State(theInteractionTimeLeft) : "
412  << G4BestUnit(State(theInteractionTimeLeft),"Time") << G4endl;
413  G4cout << "track.GetGlobalTime() : "
414  << G4BestUnit(track.GetGlobalTime(),"Time") << G4endl;
415 */
416  }
417  else // A field exerts force
418  {
419 
420  G4ExceptionDescription exceptionDescription;
421  exceptionDescription
422  << "ITTransportation does not support external fields.";
423  exceptionDescription
424  << " If you are dealing with a tradiational MC simulation, ";
425  exceptionDescription << "please use G4Transportation.";
426 
427  G4Exception("G4ITTransportation::AlongStepGetPhysicalInteractionLength",
428  "NoExternalFieldSupport", FatalException, exceptionDescription);
429  /*
430  G4double momentumMagnitude = pParticle->GetTotalMomentum() ;
431  // G4ThreeVector EndUnitMomentum ;
432  G4double lengthAlongCurve ;
433  G4double restMass = pParticleDef->GetPDGMass() ;
434 
435  fFieldPropagator->SetChargeMomentumMass( particleCharge, // in e+ units
436  momentumMagnitude, // in Mev/c
437  restMass ) ;
438 
439  G4ThreeVector spin = track.GetPolarization() ;
440  G4FieldTrack aFieldTrack = G4FieldTrack( startPosition,
441  track.GetMomentumDirection(),
442  0.0,
443  track.GetKineticEnergy(),
444  restMass,
445  track.GetVelocity(),
446  track.GetGlobalTime(), // Lab.
447  track.GetProperTime(), // Part.
448  &spin ) ;
449  if( currentMinimumStep > 0 )
450  {
451  // Do the Transport in the field (non recti-linear)
452  //
453  lengthAlongCurve = fFieldPropagator->ComputeStep( aFieldTrack,
454  currentMinimumStep,
455  currentSafety,
456  track.GetVolume() ) ;
457  State(fGeometryLimitedStep)= lengthAlongCurve < currentMinimumStep;
458  if( State(fGeometryLimitedStep) )
459  {
460  geometryStepLength = lengthAlongCurve ;
461  }
462  else
463  {
464  geometryStepLength = currentMinimumStep ;
465  }
466 
467  // Remember last safety origin & value.
468  //
469  State(fPreviousSftOrigin) = startPosition ;
470  State(fPreviousSafety) = currentSafety ;
471  fpSafetyHelper->SetCurrentSafety( newSafety, startPosition);
472  }
473  else
474  {
475  geometryStepLength = lengthAlongCurve= 0.0 ;
476  State(fGeometryLimitedStep) = false ;
477  }
478 
479  // Get the End-Position and End-Momentum (Dir-ection)
480  //
481  State(fTransportEndPosition) = aFieldTrack.GetPosition() ;
482 
483  // Momentum: Magnitude and direction can be changed too now ...
484  //
485  State(fMomentumChanged) = true ;
486  State(fTransportEndMomentumDir) = aFieldTrack.GetMomentumDir() ;
487 
488  State(fTransportEndKineticEnergy) = aFieldTrack.GetKineticEnergy() ;
489 
490  if( fFieldPropagator->GetCurrentFieldManager()->DoesFieldChangeEnergy() )
491  {
492  // If the field can change energy, then the time must be integrated
493  // - so this should have been updated
494  //
495  State(fCandidateEndGlobalTime) = aFieldTrack.GetLabTimeOfFlight();
496  State(fEndGlobalTimeComputed) = true;
497 
498  State(theInteractionTimeLeft) = State(fCandidateEndGlobalTime) -
499  track.GetGlobalTime() ;
500 
501  // was ( State(fCandidateEndGlobalTime) != track.GetGlobalTime() );
502  // a cleaner way is to have FieldTrack knowing whether time is updated.
503  }
504  else
505  {
506  // The energy should be unchanged by field transport,
507  // - so the time changed will be calculated elsewhere
508  //
509  State(fEndGlobalTimeComputed) = false;
510 
511  // Check that the integration preserved the energy
512  // - and if not correct this!
513  G4double startEnergy= track.GetKineticEnergy();
514  G4double endEnergy= State(fTransportEndKineticEnergy);
515 
516  static G4int no_inexact_steps=0, no_large_ediff;
517  G4double absEdiff = std::fabs(startEnergy- endEnergy);
518  if( absEdiff > perMillion * endEnergy )
519  {
520  no_inexact_steps++;
521  // Possible statistics keeping here ...
522  }
523  #ifdef G4VERBOSE
524  if( fVerboseLevel > 1 )
525  {
526  if( std::fabs(startEnergy- endEnergy) > perThousand * endEnergy )
527  {
528  static G4int no_warnings= 0, warnModulo=1, moduloFactor= 10;
529  no_large_ediff ++;
530  if( (no_large_ediff% warnModulo) == 0 )
531  {
532  no_warnings++;
533  G4cout << "WARNING - G4Transportation::AlongStepGetPIL() "
534  << " Energy change in Step is above 1^-3 relative value. " << G4endl
535  << " Relative change in 'tracking' step = "
536  << std::setw(15) << (endEnergy-startEnergy)/startEnergy << G4endl
537  << " Starting E= " << std::setw(12) << startEnergy / MeV << " MeV "
538  << G4endl
539  << " Ending E= " << std::setw(12) << endEnergy / MeV << " MeV "
540  << G4endl;
541  G4cout << " Energy has been corrected -- however, review"
542  << " field propagation parameters for accuracy." << G4endl;
543  if( (fVerboseLevel > 2 ) || (no_warnings<4) ||
544  (no_large_ediff == warnModulo * moduloFactor) )
545  {
546  G4cout << " These include EpsilonStepMax(/Min) in G4FieldManager "
547  << " which determine fractional error per step for integrated quantities. "
548  << G4endl
549  << " Note also the influence of the permitted number of integration steps."
550  << G4endl;
551  }
552  G4cerr << "ERROR - G4Transportation::AlongStepGetPIL()" << G4endl
553  << " Bad 'endpoint'. Energy change detected"
554  << " and corrected. "
555  << " Has occurred already "
556  << no_large_ediff << " times." << G4endl;
557  if( no_large_ediff == warnModulo * moduloFactor )
558  {
559  warnModulo *= moduloFactor;
560  }
561  }
562  }
563  } // end of if (fVerboseLevel)
564  #endif
565  // Correct the energy for fields that conserve it
566  // This - hides the integration error
567  // - but gives a better physical answer
568  State(fTransportEndKineticEnergy)= track.GetKineticEnergy();
569  }
570 
571  State(fTransportEndSpin) = aFieldTrack.GetSpin();
572  State(fParticleIsLooping) = fFieldPropagator->IsParticleLooping() ;
573  State(endpointDistance) = (State(fTransportEndPosition) -
574  startPosition).mag() ;
575  // State(theInteractionTimeLeft) =
576  track.GetVelocity()/State(endpointDistance) ;
577  */
578  }
579 
580  // If we are asked to go a step length of 0, and we are on a boundary
581  // then a boundary will also limit the step -> we must flag this.
582  //
583  if (currentMinimumStep == 0.0)
584  {
585  if (currentSafety == 0.0)
586  {
587  State(fGeometryLimitedStep) = true;
588 // G4cout << "!!!! Safety is NULL, on the Boundary !!!!!" << G4endl;
589 // G4cout << " Track position : " << track.GetPosition() /nanometer
590 // << G4endl;
591  }
592  }
593 
594  // Update the safety starting from the end-point,
595  // if it will become negative at the end-point.
596  //
597  if (currentSafety < State(fEndPointDistance))
598  {
599  // if( particleCharge == 0.0 )
600  // G4cout << " Avoiding call to ComputeSafety : charge = 0.0 " << G4endl;
601 
602  if (particleCharge != 0.0)
603  {
604 
605  G4double endSafety = fLinearNavigator->ComputeSafety(
606  State(fTransportEndPosition));
607  currentSafety = endSafety;
608  State(fPreviousSftOrigin) = State(fTransportEndPosition);
609  State(fPreviousSafety) = currentSafety;
610 
611  /*
612  G4VTrackStateHandle state =
613  GetIT(track)->GetTrackingInfo()->GetTrackState(fpSafetyHelper);
614  */
615  G4TrackStateManager& trackStateMan = GetIT(track)->GetTrackingInfo()
617  fpSafetyHelper->LoadTrackState(trackStateMan);
618  // fpSafetyHelper->SetTrackState(state);
619  fpSafetyHelper->SetCurrentSafety(currentSafety,
620  State(fTransportEndPosition));
622 
623  // Because the Stepping Manager assumes it is from the start point,
624  // add the StepLength
625  //
626  currentSafety += State(fEndPointDistance);
627 
628 #ifdef G4DEBUG_TRANSPORT
629  G4cout.precision(12);
630  G4cout << "***G4Transportation::AlongStepGPIL ** " << G4endl;
631  G4cout << " Called Navigator->ComputeSafety at "
632  << State(fTransportEndPosition)
633  << " and it returned safety= " << endSafety << G4endl;
634  G4cout << " Adding endpoint distance " << State(fEndPointDistance)
635  << " to obtain pseudo-safety= " << currentSafety << G4endl;
636 #endif
637  }
638  }
639 
640  // fParticleChange.ProposeTrueStepLength(geometryStepLength) ;
641 
642 // G4cout << "G4ITTransportation::AlongStepGetPhysicalInteractionLength = "
643 // << G4BestUnit(geometryStepLength,"Length") << G4endl;
644 
645  return geometryStepLength;
646 }
647 
649  const G4Step& /*step*/,
650  const double timeStep,
651  double& oPhysicalStep)
652 {
653  PrepareState()
654  const G4DynamicParticle* pParticle = track.GetDynamicParticle();
655  G4ThreeVector startMomentumDir = pParticle->GetMomentumDirection();
656  G4ThreeVector startPosition = track.GetPosition();
657 
658  track.CalculateVelocity();
659  G4double initialVelocity = track.GetVelocity();
660 
661  State(fGeometryLimitedStep) = false;
662 
664  // !!! CASE NO FIELD !!!
666  State(fCandidateEndGlobalTime) = timeStep + track.GetGlobalTime();
667  State(fEndGlobalTimeComputed) = true;
668 
669  // Choose the calculation of the transportation: Field or not
670  //
671  if (!State(fMomentumChanged))
672  {
673  // G4cout << "Momentum has not changed" << G4endl;
674  fParticleChange.ProposeVelocity(initialVelocity);
675  oPhysicalStep = initialVelocity * timeStep;
676 
677  // Calculate final position
678  //
679  State(fTransportEndPosition) = startPosition
680  + oPhysicalStep * startMomentumDir;
681  }
682 }
683 
685 //
686 // Initialize ParticleChange (by setting all its members equal
687 // to corresponding members in G4Track)
688 #include "G4ParticleTable.hh"
690  const G4Step& stepData)
691 {
692 
693 #if defined (DEBUG_MEM)
694  MemStat mem_first, mem_second, mem_diff;
695 #endif
696 
697 #if defined (DEBUG_MEM)
698  mem_first = MemoryUsage();
699 #endif
700 
701  PrepareState()
702 
703  // G4cout << "G4ITTransportation::AlongStepDoIt" << G4endl;
704  // set pdefOpticalPhoton
705  // Andrea Dotti: the following statement should be in a single line:
706  // G4-MT transformation tools get confused if statement spans two lines
707  // If needed contact: adotti@slac.stanford.edu
708  static G4ThreadLocal G4ParticleDefinition* pdefOpticalPhoton = 0;
709  if (!pdefOpticalPhoton) pdefOpticalPhoton =
711 
712  static G4ThreadLocal G4int noCalls = 0;
713  noCalls++;
714 
716 
717  // Code for specific process
718  //
719  fParticleChange.ProposePosition(State(fTransportEndPosition));
720  fParticleChange.ProposeMomentumDirection(State(fTransportEndMomentumDir));
721  fParticleChange.ProposeEnergy(State(fTransportEndKineticEnergy));
722  fParticleChange.SetMomentumChanged(State(fMomentumChanged));
723 
724  fParticleChange.ProposePolarization(State(fTransportEndSpin));
725 
726  G4double deltaTime = 0.0;
727 
728  // Calculate Lab Time of Flight (ONLY if field Equations used it!)
729  // G4double endTime = State(fCandidateEndGlobalTime);
730  // G4double delta_time = endTime - startTime;
731 
732  G4double startTime = track.GetGlobalTime();
736  if (State(fEndGlobalTimeComputed) == false)
737  {
738  // The time was not integrated .. make the best estimate possible
739  //
740  G4double initialVelocity = stepData.GetPreStepPoint()->GetVelocity();
741  G4double stepLength = track.GetStepLength();
742 
743  deltaTime = 0.0; // in case initialVelocity = 0
744  if (track.GetParticleDefinition() == pdefOpticalPhoton)
745  {
746  // For only Optical Photon, final velocity is used
747  double finalVelocity = track.CalculateVelocityForOpticalPhoton();
748  fParticleChange.ProposeVelocity(finalVelocity);
749  deltaTime = stepLength / finalVelocity;
750  }
751  else if (initialVelocity > 0.0)
752  {
753  deltaTime = stepLength / initialVelocity;
754  }
755 
756  State(fCandidateEndGlobalTime) = startTime + deltaTime;
757  }
758  else
759  {
760  deltaTime = State(fCandidateEndGlobalTime) - startTime;
761  }
762 
763  fParticleChange.ProposeGlobalTime(State(fCandidateEndGlobalTime));
764  fParticleChange.ProposeLocalTime(track.GetLocalTime() + deltaTime);
765  /*
766  // Now Correct by Lorentz factor to get delta "proper" Time
767 
768  G4double restMass = track.GetDynamicParticle()->GetMass() ;
769  G4double deltaProperTime = deltaTime*( restMass/track.GetTotalEnergy() ) ;
770 
771  fParticleChange.ProposeProperTime(track.GetProperTime() + deltaProperTime) ;
772  */
773 
775 
778 
779  // If the particle is caught looping or is stuck (in very difficult
780  // boundaries) in a magnetic field (doing many steps)
781  // THEN this kills it ...
782  //
783  if (State(fParticleIsLooping))
784  {
785  G4double endEnergy = State(fTransportEndKineticEnergy);
786 
787  if ((endEnergy < fThreshold_Important_Energy) || (State(fNoLooperTrials)
788  >= fThresholdTrials))
789  {
790  // Kill the looping particle
791  //
792  // G4cout << "G4ITTransportation will killed the molecule"<< G4endl;
794 
795  // 'Bare' statistics
796  fSumEnergyKilled += endEnergy;
797  if (endEnergy > fMaxEnergyKilled)
798  {
799  fMaxEnergyKilled = endEnergy;
800  }
801 
802 #ifdef G4VERBOSE
803  if ((fVerboseLevel > 1) || (endEnergy > fThreshold_Warning_Energy))
804  {
805  G4cout
806  << " G4ITTransportation is killing track that is looping or stuck "
807  << G4endl<< " This track has " << track.GetKineticEnergy() / MeV
808  << " MeV energy." << G4endl;
809  G4cout << " Number of trials = " << State(fNoLooperTrials)
810  << " No of calls to AlongStepDoIt = " << noCalls
811  << G4endl;
812  }
813 #endif
814  State(fNoLooperTrials) = 0;
815  }
816  else
817  {
818  State(fNoLooperTrials)++;
819 #ifdef G4VERBOSE
820  if ((fVerboseLevel > 2))
821  {
822  G4cout << " G4ITTransportation::AlongStepDoIt(): Particle looping - "
823  << " Number of trials = " << State(fNoLooperTrials)
824  << " No of calls to = " << noCalls << G4endl;
825  }
826 #endif
827  }
828  }
829  else
830  {
831  State(fNoLooperTrials)=0;
832  }
833 
834  // Another (sometimes better way) is to use a user-limit maximum Step size
835  // to alleviate this problem ..
836 
837  // Introduce smooth curved trajectories to particle-change
838  //
841 
842 #if defined (DEBUG_MEM)
843  mem_second = MemoryUsage();
844  mem_diff = mem_second-mem_first;
845  G4cout << "\t || MEM || End of G4ITTransportation::AlongStepDoIt, diff is: "
846  << mem_diff << G4endl;
847 #endif
848 
849  return &fParticleChange;
850 }
851 
853 //
854 // This ensures that the PostStep action is always called,
855 // so that it can do the relocation if it is needed.
856 //
857 
858 G4double
861  G4double, // previousStepSize
862  G4ForceCondition* pForceCond)
863 {
864  *pForceCond = Forced;
865  return DBL_MAX; // was kInfinity ; but convention now is DBL_MAX
866 }
867 
869 //
870 
872  const G4Step&)
873 {
874  // G4cout << "G4ITTransportation::PostStepDoIt" << G4endl;
875 
876  PrepareState()
877  G4TouchableHandle retCurrentTouchable; // The one to return
878  G4bool isLastStep = false;
879 
880  // Initialize ParticleChange (by setting all its members equal
881  // to corresponding members in G4Track)
882  fParticleChange.Initialize(track); // To initialise TouchableChange
883 
885 
886  // If the Step was determined by the volume boundary,
887  // logically relocate the particle
888 
889  if (State(fGeometryLimitedStep))
890  {
891 
892  if(fVerboseLevel)
893  {
894  G4cout << "Step is limited by geometry "
895  << "track ID : " << track.GetTrackID() << G4endl;
896  }
897 
898  // fCurrentTouchable will now become the previous touchable,
899  // and what was the previous will be freed.
900  // (Needed because the preStepPoint can point to the previous touchable)
901 
902  if ( State(fCurrentTouchableHandle)->GetVolume() == 0)
903  {
904  G4ExceptionDescription exceptionDescription;
905  exceptionDescription << "No current touchable found ";
906  G4Exception(" G4ITTransportation::PostStepDoIt", "G4ITTransportation001",
907  FatalErrorInArgument, exceptionDescription);
908  }
909 
910  fLinearNavigator->SetGeometricallyLimitedStep();
911  fLinearNavigator->LocateGlobalPointAndUpdateTouchableHandle(
912  track.GetPosition(), track.GetMomentumDirection(),
913  State(fCurrentTouchableHandle), true);
914  // Check whether the particle is out of the world volume
915  // If so it has exited and must be killed.
916  //
917  if ( State(fCurrentTouchableHandle)->GetVolume() == 0)
918  {
919  // abort();
920 #ifdef G4VERBOSE
921  if (fVerboseLevel > 0)
922  {
923  G4cout << "Track position : " << track.GetPosition() / nanometer
924  << " [nm]" << " Track ID : " << track.GetTrackID() << G4endl;
925  G4cout << "G4ITTransportation will killed the track because "
926  "State(fCurrentTouchableHandle)->GetVolume() == 0"<< G4endl;
927  }
928 #endif
930  }
931 
932  retCurrentTouchable = State(fCurrentTouchableHandle);
933 
934 // G4cout << "Current volume : " << track.GetVolume()->GetName()
935 // << " Next volume : "
936 // << (State(fCurrentTouchableHandle)->GetVolume() ?
937 // State(fCurrentTouchableHandle)->GetVolume()->GetName():"OutWorld")
938 // << " Position : " << track.GetPosition() / nanometer
939 // << " track ID : " << track.GetTrackID()
940 // << G4endl;
941 
942  fParticleChange.SetTouchableHandle(State(fCurrentTouchableHandle));
943 
944  // Update the Step flag which identifies the Last Step in a volume
945  isLastStep = fLinearNavigator->ExitedMotherVolume()
946  | fLinearNavigator->EnteredDaughterVolume();
947 
948 #ifdef G4DEBUG_TRANSPORT
949  // Checking first implementation of flagging Last Step in Volume
950  G4bool exiting = fLinearNavigator->ExitedMotherVolume();
951  G4bool entering = fLinearNavigator->EnteredDaughterVolume();
952 
953  if( ! (exiting || entering) )
954  {
955  G4cout << " Transport> : Proposed isLastStep= " << isLastStep
956  << " Exiting " << fLinearNavigator->ExitedMotherVolume()
957  << " Entering " << fLinearNavigator->EnteredDaughterVolume()
958  << " Track position : " << track.GetPosition() /nanometer << " [nm]"
959  << G4endl;
960  G4cout << " Track position : " << track.GetPosition() /nanometer
961  << G4endl;
962  }
963 #endif
964  }
965  else // fGeometryLimitedStep is false
966  {
967  // This serves only to move the Navigator's location
968  //
969 // abort();
970  fLinearNavigator->LocateGlobalPointWithinVolume(track.GetPosition());
971 
972  // The value of the track's current Touchable is retained.
973  // (and it must be correct because we must use it below to
974  // overwrite the (unset) one in particle change)
975  // It must be fCurrentTouchable too ??
976  //
978  retCurrentTouchable = track.GetTouchableHandle();
979 
980  isLastStep = false;
981 #ifdef G4DEBUG_TRANSPORT
982  // Checking first implementation of flagging Last Step in Volume
983  G4cout << " Transport> Proposed isLastStep= " << isLastStep
984  << " Geometry did not limit step. Position : "
985  << track.GetPosition()/ nanometer << G4endl;
986 #endif
987  } // endif ( fGeometryLimitedStep )
988 
990 
991  const G4VPhysicalVolume* pNewVol = retCurrentTouchable->GetVolume();
992  const G4Material* pNewMaterial = 0;
993  const G4VSensitiveDetector* pNewSensitiveDetector = 0;
994 
995  if (pNewVol != 0)
996  {
997  pNewMaterial = pNewVol->GetLogicalVolume()->GetMaterial();
998  pNewSensitiveDetector = pNewVol->GetLogicalVolume()->GetSensitiveDetector();
999  }
1000 
1001  // ( <const_cast> pNewMaterial ) ;
1002  // ( <const_cast> pNewSensitiveDetector) ;
1003 
1006  (G4VSensitiveDetector *) pNewSensitiveDetector);
1007 
1008  const G4MaterialCutsCouple* pNewMaterialCutsCouple = 0;
1009  if (pNewVol != 0)
1010  {
1011  pNewMaterialCutsCouple =
1013  }
1014 
1015  if (pNewVol != 0 && pNewMaterialCutsCouple != 0
1016  && pNewMaterialCutsCouple->GetMaterial() != pNewMaterial)
1017  {
1018  // for parametrized volume
1019  //
1020  pNewMaterialCutsCouple = G4ProductionCutsTable::GetProductionCutsTable()
1021  ->GetMaterialCutsCouple(pNewMaterial,
1022  pNewMaterialCutsCouple->GetProductionCuts());
1023  }
1024  fParticleChange.SetMaterialCutsCoupleInTouchable(pNewMaterialCutsCouple);
1025 
1026  // temporarily until Get/Set Material of ParticleChange,
1027  // and StepPoint can be made const.
1028  // Set the touchable in ParticleChange
1029  // this must always be done because the particle change always
1030  // uses this value to overwrite the current touchable pointer.
1031  //
1032  fParticleChange.SetTouchableHandle(retCurrentTouchable);
1033 
1034  return &fParticleChange;
1035 }
1036 
1037 // New method takes over the responsibility to reset the state of
1038 // G4Transportation object at the start of a new track or the resumption of
1039 // a suspended track.
1040 
1042 {
1044  if (fInstantiateProcessState)
1045  {
1046 // G4VITProcess::fpState = new G4ITTransportationState();
1048  // Will set in the same time fTransportationState
1049  }
1050 
1053  GetIT(track)->GetTrackingInfo()->GetTrackStateManager());
1054 
1055  // The actions here are those that were taken in AlongStepGPIL
1056  // when track.GetCurrentStepNumber()==1
1057 
1058  // reset safety value and center
1059  //
1060  // State(fPreviousSafety) = 0.0 ;
1061  // State(fPreviousSftOrigin) = G4ThreeVector(0.,0.,0.) ;
1062 
1063  // reset looping counter -- for motion in field
1064  // State(fNoLooperTrials)= 0;
1065  // Must clear this state .. else it depends on last track's value
1066  // --> a better solution would set this from state of suspended track TODO ?
1067  // Was if( aTrack->GetCurrentStepNumber()==1 ) { .. }
1068 
1069  // ChordFinder reset internal state
1070  //
1071  if (DoesGlobalFieldExist())
1072  {
1074  // Resets all state of field propagator class (ONLY)
1075  // including safety values (in case of overlaps and to wipe for first track).
1076 
1077  // G4ChordFinder* chordF= fFieldPropagator->GetChordFinder();
1078  // if( chordF ) chordF->ResetStepEstimate();
1079  }
1080 
1081  // Make sure to clear the chord finders of all fields (ie managers)
1082  static G4ThreadLocal G4FieldManagerStore* fieldMgrStore = 0;
1083  if (!fieldMgrStore) fieldMgrStore = G4FieldManagerStore::GetInstance();
1084  fieldMgrStore->ClearAllChordFindersState();
1085 
1086  // Update the current touchable handle (from the track's)
1087  //
1088  PrepareState()
1089  State(fCurrentTouchableHandle) = track->GetTouchableHandle();
1090 
1092 }
1093 
1094 #undef State
1095 #undef PrepareState
void SetMaterialInTouchable(G4Material *fMaterial)
G4ITSafetyHelper * GetSafetyHelper() const
void SetTouchableHandle(const G4TouchableHandle &fTouchable)
const G4ThreeVector & GetPolarization() const
G4ParticleDefinition * FindParticle(G4int PDGEncoding)
G4double GetLocalTime() const
G4ParticleChangeForTransport fParticleChange
G4Material * GetMaterial() const
void SetMaterialCutsCoupleInTouchable(const G4MaterialCutsCouple *fMaterialCutsCouple)
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
std::vector< G4ThreeVector > * GimmeTrajectoryVectorAndForgetIt() const
CLHEP::Hep3Vector G4ThreeVector
G4double GetVelocity() const
virtual void BuildPhysicsTable(const G4ParticleDefinition &)
const G4DynamicParticle * GetDynamicParticle() const
const G4ThreeVector & GetPosition() const
G4TrackStatus GetTrackStatus() const
virtual G4VParticleChange * PostStepDoIt(const G4Track &track, const G4Step &)
void ProposePolarization(G4double Px, G4double Py, G4double Pz)
static constexpr double nanometer
Definition: G4SIunits.hh:101
virtual G4VParticleChange * AlongStepDoIt(const G4Track &track, const G4Step &stepData)
G4double GetVelocity() const
virtual G4double AlongStepGetPhysicalInteractionLength(const G4Track &track, G4double, G4double currentMinimumStep, G4double &currentSafety, G4GPILSelection *selection)
#define G4ThreadLocal
Definition: tls.hh:89
G4bool enableAtRestDoIt
Definition: G4VProcess.hh:350
#define State(theXInfo)
int G4int
Definition: G4Types.hh:78
MemStat MemoryUsage()
Definition: G4MemStat.cc:55
void SetInstantiateProcessState(G4bool flag)
void ProposePosition(G4double x, G4double y, G4double z)
G4ReferenceCountedHandle< G4VTouchable > G4TouchableHandle
G4ITNavigator * fLinearNavigator
#define PrepareState()
G4StepPoint * GetPreStepPoint() const
virtual void Initialize(const G4Track &)
virtual void StartTracking(G4Track *)
Definition: G4VProcess.cc:101
G4double CalculateVelocityForOpticalPhoton() const
Definition: G4Track.cc:255
virtual void ConfigureForTrack(const G4Track *)
virtual void StartTracking(G4Track *aTrack)
G4IT * GetIT(const G4Track *track)
Definition: G4IT.cc:49
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
virtual void ResetTrackState()
G4double GetCharge() const
void ProposeTrueStepLength(G4double truePathLength)
G4bool DoesFieldExist() const
virtual void StartTracking(G4Track *)
Definition: G4VITProcess.cc:86
void SetInstantiateProcessState(G4bool flag)
const G4ParticleDefinition * GetParticleDefinition() const
void SetProcessSubType(G4int)
Definition: G4VProcess.hh:432
G4shared_ptr< G4ProcessState > fpState
Definition: G4Step.hh:76
G4bool enablePostStepDoIt
Definition: G4VProcess.hh:352
G4int GetTrackID() const
G4double GetGlobalTime() const
G4double CalculateVelocity() const
Definition: G4Track.cc:222
G4FieldManager * FindAndSetFieldManager(G4VPhysicalVolume *pCurrentPhysVol)
const G4TouchableHandle & GetTouchableHandle() const
G4ITSafetyHelper * fpSafetyHelper
virtual void NewTrackState()
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
static G4TransportationManager * GetTransportationManager()
G4FieldManager * GetFieldManager() const
G4TrackingInformation * GetTrackingInfo()
Definition: G4IT.hh:144
static G4ProductionCutsTable * GetProductionCutsTable()
const G4ThreeVector & GetMomentumDirection() const
G4LogicalVolume * GetLogicalVolume() const
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
static G4FieldManagerStore * GetInstance()
void SetSensitiveDetectorInTouchable(G4VSensitiveDetector *fSensitiveDetector)
static G4ParticleTable * GetParticleTable()
#define fPreviousSftOrigin
virtual G4VPhysicalVolume * GetVolume(G4int depth=0) const
Definition: G4VTouchable.cc:44
void SetPointerToVectorOfAuxiliaryPoints(std::vector< G4ThreeVector > *theNewVectorPointer)
void ProposeGlobalTime(G4double t)
G4ITTransportation(const G4String &aName="ITTransportation", G4int verbosityLevel=0)
#define fPreviousSafety
G4ITNavigator * GetNavigatorForTracking() const
G4VParticleChange * pParticleChange
Definition: G4VProcess.hh:283
void ProposeEnergy(G4double finalEnergy)
double mag2() const
G4VPhysicalVolume * GetVolume() const
const G4MaterialCutsCouple * GetMaterialCutsCouple() const
#define G4endl
Definition: G4ios.hh:61
static constexpr double MeV
Definition: G4SIunits.hh:214
G4bool enableAlongStepDoIt
Definition: G4VProcess.hh:351
void ProposeLastStepInVolume(G4bool flag)
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
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
G4TrackStateManager & GetTrackStateManager()
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
virtual void SaveTrackState(G4TrackStateManager &manager)
void ProposeLocalTime(G4double t)
void SetCurrentSafety(G4double val, const G4ThreeVector &pos)
G4double GetStepLength() const
virtual void LoadTrackState(G4TrackStateManager &manager)
G4GPILSelection
const G4Material * GetMaterial() const