Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/dna/management/src/G4ITTransportation.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  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 /// \brief This class is a slightly modified version of G4Transportation
 28 ///        initially written by John Apostolakis and colleagues
 29 ///        But it should use the exact same algorithm
 30 //
 31 // Contact : Mathieu Karamitros (kara (AT) cenbg . in2p3 . fr)
 32 //
 33 // History :
 34 // -----------
 35 // =======================================================================
 36 // Modified:
 37 //   28 Oct  2011, P.Gumpl./J.Ap: Detect gravity field, use magnetic moment
 38 //   20 Nov  2008, J.Apostolakis: Push safety to helper - after ComputeSafety
 39 //    9 Nov  2007, J.Apostolakis: Flag for short steps, push safety to helper
 40 //   19 Jan  2006, P.MoraDeFreitas: Fix for suspended tracks (StartTracking)
 41 //   11 Aug  2004, M.Asai: Add G4VSensitiveDetector* for updating stepPoint.
 42 //   21 June 2003, J.Apostolakis: Calling field manager with
 43 //                     track, to enable it to configure its accuracy
 44 //   13 May  2003, J.Apostolakis: Zero field areas now taken into
 45 //                     account correclty in all cases (thanks to W Pokorski).
 46 //   29 June 2001, J.Apostolakis, D.Cote-Ahern, P.Gumplinger:
 47 //                     correction for spin tracking
 48 //   20 Febr 2001, J.Apostolakis:  update for new FieldTrack
 49 //   22 Sept 2000, V.Grichine:     update of Kinetic Energy
 50 //  ---------------------------------------------------
 51 //   10 Oct  2011, M.Karamitros:   G4ITTransportation created
 52 // Created:  19 March 1997, J. Apostolakis
 53 // =======================================================================
 54 //
 55 // -------------------------------------------------------------------
 56 
 57 #include "G4ITTransportation.hh"
 58 
 59 #include "G4ChordFinder.hh"
 60 #include "G4FieldManager.hh"
 61 #include "G4FieldManagerStore.hh"
 62 #include "G4IT.hh"
 63 #include "G4ITNavigator.hh"
 64 #include "G4ITSafetyHelper.hh"
 65 #include "G4ITTransportationManager.hh"
 66 #include "G4LowEnergyEmProcessSubType.hh"
 67 #include "G4ParticleTable.hh"
 68 #include "G4ProductionCutsTable.hh"
 69 #include "G4PropagatorInField.hh"
 70 #include "G4ReferenceCast.hh"
 71 #include "G4SystemOfUnits.hh"
 72 #include "G4TrackingInformation.hh"
 73 #include "G4TransportationManager.hh"
 74 #include "G4UnitsTable.hh"
 75 
 76 #include <memory>
 77 
 78 class G4VSensitiveDetector;
 79 
 80 #ifndef PrepareState
 81 #  define PrepareState()                                                       \
 82     G4ITTransportationState* __state = this->GetState<G4ITTransportationState>()
 83 #endif
 84 
 85 #ifndef State
 86 #define State(theXInfo) (__state->theXInfo)
 87 #endif
 88 
 89 //#define DEBUG_MEM
 90 
 91 #ifdef DEBUG_MEM
 92 #include "G4MemStat.hh"
 93 using namespace G4MemStat;
 94 using G4MemStat::MemStat;
 95 #endif
 96 
 97 //#define G4DEBUG_TRANSPORT 1
 98 
 99 G4ITTransportation::G4ITTransportation(const G4String& aName, int verbose) :
100     G4VITProcess(aName, fTransportation),
101     fThreshold_Warning_Energy(100 * MeV),
102     fThreshold_Important_Energy(250 * MeV),
103     
104     fUnimportant_Energy(1 * MeV),  // Old default: true (=fast short steps)
105     fVerboseLevel(verbose)
106 {
107   pParticleChange = &fParticleChange;
108   G4TransportationManager* transportMgr;
109   transportMgr = G4TransportationManager::GetTransportationManager();
110   G4ITTransportationManager* ITtransportMgr;
111   ITtransportMgr = G4ITTransportationManager::GetTransportationManager();
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(fLowEnergyTransportation);
125   SetInstantiateProcessState(true);
126   G4VITProcess::SetInstantiateProcessState(false);
127   fInstantiateProcessState = true;
128 
129   G4VITProcess::fpState = std::make_shared<G4ITTransportationState>();
130   /*
131    if(fTransportationState == 0)
132    {
133    G4cout << "KILL in G4ITTransportation::G4ITTransportation" << G4endl;
134    abort();
135    }
136    */
137 }
138 
139 G4ITTransportation::G4ITTransportation(const G4ITTransportation& right) :
140     G4VITProcess(right)
141 {
142   // Copy attributes
143   fVerboseLevel = right.fVerboseLevel;
144   fThreshold_Warning_Energy = right.fThreshold_Warning_Energy;
145   fThreshold_Important_Energy = right.fThreshold_Important_Energy;
146   fThresholdTrials = right.fThresholdTrials;
147   fUnimportant_Energy = right.fUnimportant_Energy;
148   fSumEnergyKilled = right.fSumEnergyKilled;
149   fMaxEnergyKilled = right.fMaxEnergyKilled;
150   fShortStepOptimisation = right.fShortStepOptimisation;
151 
152   // Setup Navigators
153   G4TransportationManager* transportMgr;
154   transportMgr = G4TransportationManager::GetTransportationManager();
155   G4ITTransportationManager* ITtransportMgr;
156   ITtransportMgr = G4ITTransportationManager::GetTransportationManager();
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 
170   pParticleChange = &fParticleChange;
171   SetInstantiateProcessState(true);
172   G4VITProcess::SetInstantiateProcessState(false);
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 
182 //////////////////////////////////////////////////////////////////////////////
183 /// Process State
184 //////////////////////////////////////////////////////////////////////////////
185 G4ITTransportation::G4ITTransportationState::G4ITTransportationState() :
186      fCurrentTouchableHandle(nullptr)
187 {
188   fTransportEndPosition = G4ThreeVector(0, 0, 0);
189   fTransportEndMomentumDir = G4ThreeVector(0, 0, 0);
190   fTransportEndKineticEnergy = -1;
191   fTransportEndSpin = G4ThreeVector(0, 0, 0);
192   fMomentumChanged = false;
193   fEnergyChanged = false;
194   fEndGlobalTimeComputed = false;
195   fCandidateEndGlobalTime = -1;
196   fParticleIsLooping = false;
197 
198   static G4ThreadLocal G4TouchableHandle *nullTouchableHandle = nullptr;
199   if (nullTouchableHandle == nullptr) nullTouchableHandle = new G4TouchableHandle;
200   // Points to (G4VTouchable*) 0
201 
202   fCurrentTouchableHandle = *nullTouchableHandle;
203   fGeometryLimitedStep = false;
204   fPreviousSftOrigin = G4ThreeVector(0, 0, 0);
205   fPreviousSafety = 0.0;
206   fNoLooperTrials = 0;
207   fEndPointDistance = -1;
208 }
209 
210 G4ITTransportation::G4ITTransportationState::~G4ITTransportationState()
211 {
212   ;
213 }
214 
215 G4ITTransportation::~G4ITTransportation()
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 
230 void G4ITTransportation::BuildPhysicsTable(const G4ParticleDefinition&)
231 {
232   fpSafetyHelper->InitialiseHelper();
233 }
234 
235 G4bool G4ITTransportation::DoesGlobalFieldExist()
236 {
237   G4TransportationManager* transportMgr;
238   transportMgr = G4TransportationManager::GetTransportationManager();
239 
240   // fFieldExists= transportMgr->GetFieldManager()->DoesFieldExist();
241   // return fFieldExists;
242   return transportMgr->GetFieldManager()->DoesFieldExist();
243 }
244 
245 //////////////////////////////////////////////////////////////////////////
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
252 G4ITTransportation::
253 AlongStepGetPhysicalInteractionLength(const G4Track& track,
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 = nullptr;
312   G4bool fieldExertsForce = false;
313   if ((particleCharge != 0.0))
314   {
315     fieldMgr = fFieldPropagator->FindAndSetFieldManager(track.GetVolume());
316     if (fieldMgr != nullptr)
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() != nullptr);
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()
364           ->GetTrackStateManager();
365       fpSafetyHelper->LoadTrackState(trackStateMan);
366       // fpSafetyHelper->SetTrackState(state);
367       fpSafetyHelper->SetCurrentSafety(newSafety,
368                                        State(fTransportEndPosition));
369       fpSafetyHelper->ResetTrackState();
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()
616           ->GetTrackStateManager();
617       fpSafetyHelper->LoadTrackState(trackStateMan);
618       // fpSafetyHelper->SetTrackState(state);
619       fpSafetyHelper->SetCurrentSafety(currentSafety,
620                                        State(fTransportEndPosition));
621       fpSafetyHelper->ResetTrackState();
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 
648 void G4ITTransportation::ComputeStep(const G4Track& track,
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 
663   /////////////////////////
664   // !!! CASE NO FIELD !!!
665   /////////////////////////
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 
684 //////////////////////////////////////////////////////////////////////////
685 //
686 //   Initialize ParticleChange  (by setting all its members equal
687 //                               to corresponding members in G4Track)
688 #include "G4ParticleTable.hh"
689 G4VParticleChange* G4ITTransportation::AlongStepDoIt(const G4Track& track,
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 = nullptr;
709   if (pdefOpticalPhoton == nullptr) pdefOpticalPhoton =
710       G4ParticleTable::GetParticleTable()->FindParticle("opticalphoton");
711 
712   static G4ThreadLocal G4int noCalls = 0;
713   noCalls++;
714 
715   fParticleChange.Initialize(track);
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();
733   ///___________________________________________________________________________
734   /// !!!!!!!
735   /// A REVOIR !!!!
736   if (State(fEndGlobalTimeComputed) == 0)
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 
774   fParticleChange.ProposeTrueStepLength(track.GetStepLength());
775 
776   ///___________________________________________________________________________
777   ///
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;
793       fParticleChange.ProposeTrackStatus(fStopAndKill);
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   //
839   fParticleChange.SetPointerToVectorOfAuxiliaryPoints(
840       fFieldPropagator->GimmeTrajectoryVectorAndForgetIt());
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 
852 //////////////////////////////////////////////////////////////////////////
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
859 G4ITTransportation::
860 PostStepGetPhysicalInteractionLength(const G4Track&, // track
861                                      G4double, // previousStepSize
862                                      G4ForceCondition* pForceCond)
863 {
864   *pForceCond = Forced;
865   return DBL_MAX; // was kInfinity ; but convention now is DBL_MAX
866 }
867 
868 /////////////////////////////////////////////////////////////////////////////
869 //
870 
871 G4VParticleChange* G4ITTransportation::PostStepDoIt(const G4Track& track,
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 
884   fParticleChange.ProposeTrackStatus(track.GetTrackStatus());
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 != 0)
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() == nullptr)
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() == nullptr)
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
929       fParticleChange.ProposeTrackStatus( fStopAndKill );
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     //
977     fParticleChange.SetTouchableHandle(track.GetTouchableHandle());
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 
989   fParticleChange.ProposeLastStepInVolume(isLastStep);
990 
991   const G4VPhysicalVolume* pNewVol = retCurrentTouchable->GetVolume();
992   const G4Material* pNewMaterial = nullptr;
993   G4VSensitiveDetector* pNewSensitiveDetector = nullptr;
994 
995   if (pNewVol != nullptr)
996   {
997     pNewMaterial = pNewVol->GetLogicalVolume()->GetMaterial();
998     pNewSensitiveDetector = pNewVol->GetLogicalVolume()->GetSensitiveDetector();
999   }
1000 
1001   // ( <const_cast> pNewMaterial ) ;
1002 
1003   fParticleChange.SetMaterialInTouchable((G4Material *) pNewMaterial);
1004   fParticleChange.SetSensitiveDetectorInTouchable(pNewSensitiveDetector);
1005 
1006   const G4MaterialCutsCouple* pNewMaterialCutsCouple = nullptr;
1007   if (pNewVol != nullptr)
1008   {
1009     pNewMaterialCutsCouple =
1010         pNewVol->GetLogicalVolume()->GetMaterialCutsCouple();
1011   }
1012 
1013   if (pNewVol != nullptr && pNewMaterialCutsCouple != nullptr
1014       && pNewMaterialCutsCouple->GetMaterial() != pNewMaterial)
1015   {
1016     // for parametrized volume
1017     //
1018     pNewMaterialCutsCouple = G4ProductionCutsTable::GetProductionCutsTable()
1019         ->GetMaterialCutsCouple(pNewMaterial,
1020                                 pNewMaterialCutsCouple->GetProductionCuts());
1021   }
1022   fParticleChange.SetMaterialCutsCoupleInTouchable(pNewMaterialCutsCouple);
1023 
1024   // temporarily until Get/Set Material of ParticleChange,
1025   // and StepPoint can be made const.
1026   // Set the touchable in ParticleChange
1027   // this must always be done because the particle change always
1028   // uses this value to overwrite the current touchable pointer.
1029   //
1030   fParticleChange.SetTouchableHandle(retCurrentTouchable);
1031 
1032   return &fParticleChange;
1033 }
1034 
1035 // New method takes over the responsibility to reset the state of
1036 // G4Transportation object at the start of a new track or the resumption of
1037 // a suspended track.
1038 
1039 void G4ITTransportation::StartTracking(G4Track* track)
1040 {
1041   G4VProcess::StartTracking(track);
1042   if (fInstantiateProcessState)
1043   {
1044 //        G4VITProcess::fpState = new G4ITTransportationState();
1045     G4VITProcess::fpState = std::make_shared<G4ITTransportationState>();
1046     // Will set in the same time fTransportationState
1047   }
1048 
1049   fpSafetyHelper->NewTrackState();
1050   fpSafetyHelper->SaveTrackState(
1051       GetIT(track)->GetTrackingInfo()->GetTrackStateManager());
1052 
1053   // The actions here are those that were taken in AlongStepGPIL
1054   //   when track.GetCurrentStepNumber()==1
1055 
1056   // reset safety value and center
1057   //
1058   //    State(fPreviousSafety)    = 0.0 ;
1059   //    State(fPreviousSftOrigin) = G4ThreeVector(0.,0.,0.) ;
1060 
1061   // reset looping counter -- for motion in field
1062   //    State(fNoLooperTrials)= 0;
1063   // Must clear this state .. else it depends on last track's value
1064   //  --> a better solution would set this from state of suspended track TODO ?
1065   // Was if( aTrack->GetCurrentStepNumber()==1 ) { .. }
1066 
1067   // ChordFinder reset internal state
1068   //
1069   if (DoesGlobalFieldExist())
1070   {
1071     fFieldPropagator->ClearPropagatorState();
1072     // Resets all state of field propagator class (ONLY)
1073     // including safety values (in case of overlaps and to wipe for first track).
1074 
1075     // G4ChordFinder* chordF= fFieldPropagator->GetChordFinder();
1076     // if( chordF ) chordF->ResetStepEstimate();
1077   }
1078 
1079   // Make sure to clear the chord finders of all fields (ie managers)
1080   static G4ThreadLocal G4FieldManagerStore* fieldMgrStore = nullptr;
1081   if (fieldMgrStore == nullptr) fieldMgrStore = G4FieldManagerStore::GetInstance();
1082   fieldMgrStore->ClearAllChordFindersState();
1083 
1084   // Update the current touchable handle  (from the track's)
1085   //
1086   PrepareState();
1087   State(fCurrentTouchableHandle) = track->GetTouchableHandle();
1088 
1089   G4VITProcess::StartTracking(track);
1090 }
1091 
1092 #undef State
1093 #undef PrepareState
1094