Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/transportation/src/G4CoupledTransportation.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 ]

Diff markup

Differences between /processes/transportation/src/G4CoupledTransportation.cc (Version 11.3.0) and /processes/transportation/src/G4CoupledTransportation.cc (Version 10.7.p2)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
 26 //                                                 26 //
 27 //                                                 27 //
 28 // -------------------------------------------     28 // ------------------------------------------------------------
 29 //  GEANT 4 class implementation                   29 //  GEANT 4 class implementation
 30 //                                                 30 //
 31 // ===========================================     31 // =======================================================================
 32 // Created:  19 March 1997, J. Apostolakis         32 // Created:  19 March 1997, J. Apostolakis
 33 // ===========================================     33 // =======================================================================
 34                                                    34 
 35 #include "G4CoupledTransportation.hh"              35 #include "G4CoupledTransportation.hh"
 36 #include "G4TransportationProcessType.hh"          36 #include "G4TransportationProcessType.hh"
 37 #include "G4TransportationLogger.hh"               37 #include "G4TransportationLogger.hh"
 38                                                    38 
 39 #include "G4PhysicalConstants.hh"                  39 #include "G4PhysicalConstants.hh"
 40 #include "G4SystemOfUnits.hh"                      40 #include "G4SystemOfUnits.hh"
 41 #include "G4ProductionCutsTable.hh"                41 #include "G4ProductionCutsTable.hh"
 42 #include "G4ParticleTable.hh"                      42 #include "G4ParticleTable.hh"
 43 #include "G4ChordFinder.hh"                        43 #include "G4ChordFinder.hh"
 44 #include "G4Field.hh"                              44 #include "G4Field.hh"
 45 #include "G4FieldTrack.hh"                         45 #include "G4FieldTrack.hh"
 46 #include "G4FieldManagerStore.hh"                  46 #include "G4FieldManagerStore.hh"
 47 #include "G4PathFinder.hh"                     << 
 48                                                    47 
 49 #include "G4PropagatorInField.hh"              <<  48 #include "G4Transportation.hh"    //  In order to use fUseMagneticMoment
 50 #include "G4TransportationManager.hh"          << 
 51                                                    49 
 52 class G4VSensitiveDetector;                        50 class G4VSensitiveDetector;
 53                                                    51 
 54 G4bool G4CoupledTransportation::fSignifyStepIn     52 G4bool G4CoupledTransportation::fSignifyStepInAnyVolume= false;
 55 // This mode must apply to all threads             53 // This mode must apply to all threads 
 56                                                    54 
                                                   >>  55 G4bool G4CoupledTransportation::fUseMagneticMoment=false;
                                                   >>  56 G4bool G4CoupledTransportation::fUseGravity= false;
                                                   >>  57 G4bool G4CoupledTransportation::fSilenceLooperWarnings= false;
 57 //////////////////////////////////////////////     58 //////////////////////////////////////////////////////////////////////////
 58 //                                                 59 //
 59 // Constructor                                     60 // Constructor
 60                                                    61 
 61 G4CoupledTransportation::G4CoupledTransportati     62 G4CoupledTransportation::G4CoupledTransportation( G4int verbosity )
 62   : G4Transportation( verbosity, "CoupledTrans <<  63   : G4VProcess( G4String("CoupledTransportation"), fTransportation ),
                                                   >>  64     fTransportEndPosition(0.0, 0.0, 0.0),
                                                   >>  65     fTransportEndMomentumDir(0.0, 0.0, 0.0),
                                                   >>  66     fTransportEndKineticEnergy(0.0), 
                                                   >>  67     fTransportEndSpin(0.0, 0.0, 0.0),
                                                   >>  68     fMomentumChanged(false), 
                                                   >>  69     fEndGlobalTimeComputed(false),
                                                   >>  70     fCandidateEndGlobalTime(0.0),
                                                   >>  71     fParticleIsLooping( false ),
                                                   >>  72     fNewTrack( true ),
                                                   >>  73     fPreviousSftOrigin( 0.,0.,0. ),
 63     fPreviousMassSafety( 0.0 ),                    74     fPreviousMassSafety( 0.0 ),
 64     fPreviousFullSafety( 0.0 ),                    75     fPreviousFullSafety( 0.0 ),
 65     fMassGeometryLimitedStep( false ),             76     fMassGeometryLimitedStep( false ), 
 66     fFirstStepInMassVolume( true )             <<  77     fAnyGeometryLimitedStep( false ), 
                                                   >>  78     fEndpointDistance( -1.0 ), 
                                                   >>  79     fSumEnergyKilled( 0.0 ), fMaxEnergyKilled( 0.0 ), 
                                                   >>  80     fFirstStepInMassVolume( true ),
                                                   >>  81     fFirstStepInAnyVolume( true )
 67 {                                                  82 {
 68   SetProcessSubType(static_cast<G4int>(COUPLED     83   SetProcessSubType(static_cast<G4int>(COUPLED_TRANSPORTATION));
 69   // SetVerboseLevel is called in the construc <<  84   SetVerboseLevel(verbosity);
 70                                                    85 
                                                   >>  86   G4TransportationManager* transportMgr ; 
                                                   >>  87 
                                                   >>  88   transportMgr = G4TransportationManager::GetTransportationManager() ; 
                                                   >>  89 
                                                   >>  90   fMassNavigator = transportMgr->GetNavigatorForTracking() ; 
                                                   >>  91   fFieldPropagator = transportMgr->GetPropagatorInField() ;
                                                   >>  92   // fGlobalFieldMgr = transportMgr->GetFieldManager() ;
                                                   >>  93   fNavigatorId= transportMgr->ActivateNavigator( fMassNavigator ); 
 71   if( verboseLevel > 0 )                           94   if( verboseLevel > 0 )
 72   {                                                95   {
 73     G4cout << " G4CoupledTransportation constr     96     G4cout << " G4CoupledTransportation constructor: ----- " << G4endl;
 74     G4cout << " Verbose level is " << verboseL     97     G4cout << " Verbose level is " << verboseLevel << G4endl;
                                                   >>  98     G4cout << " Navigator Id obtained in G4CoupledTransportation constructor " 
                                                   >>  99            << fNavigatorId << G4endl;
 75     G4cout << " Reports First/Last in "           100     G4cout << " Reports First/Last in " 
 76            << (fSignifyStepInAnyVolume ? " any    101            << (fSignifyStepInAnyVolume ? " any " : " mass " )
 77            << " geometry " << G4endl;             102            << " geometry " << G4endl;
 78   }                                               103   }
 79   fPathFinder=  G4PathFinder::GetInstance();      104   fPathFinder=  G4PathFinder::GetInstance(); 
                                                   >> 105   fpSafetyHelper = transportMgr->GetSafetyHelper();  // New 
                                                   >> 106 
                                                   >> 107   fpLogger = new G4TransportationLogger("G4Transportation", verbosity);
                                                   >> 108 
                                                   >> 109   SetHighLooperThresholds();
                                                   >> 110   // Use the old defaults: Warning = 100 MeV, Important = 250 MeV, No Trials = 10;
                                                   >> 111   
                                                   >> 112   PushThresholdsToLogger();
                                                   >> 113   // Should be done by Set methods in SetHighLooperThresholds -- making sure
                                                   >> 114   
                                                   >> 115   static G4ThreadLocal G4TouchableHandle* pNullTouchableHandle = 0;
                                                   >> 116   if ( !pNullTouchableHandle)  { pNullTouchableHandle = new G4TouchableHandle; }
                                                   >> 117   fCurrentTouchableHandle = *pNullTouchableHandle;
                                                   >> 118     // Points to (G4VTouchable*) 0
                                                   >> 119 
                                                   >> 120   fAnyFieldExists= DoesAnyFieldExist();
 80 }                                                 121 }
 81                                                   122 
 82 //////////////////////////////////////////////    123 //////////////////////////////////////////////////////////////////////////
 83                                                   124 
 84 G4CoupledTransportation::~G4CoupledTransportat    125 G4CoupledTransportation::~G4CoupledTransportation()
 85 {                                                 126 {
                                                   >> 127   if( fSumEnergyKilled > 0.0 )
                                                   >> 128   {
                                                   >> 129      PrintStatistics( G4cout );
                                                   >> 130   }
                                                   >> 131   delete fpLogger;
                                                   >> 132 }
                                                   >> 133 
                                                   >> 134 //////////////////////////////////////////////////////////////////////////
                                                   >> 135 
                                                   >> 136 void
                                                   >> 137 G4CoupledTransportation::PrintStatistics( std::ostream& outStr) const
                                                   >> 138 { 
                                                   >> 139   if( fSumEnergyKilled > 0.0 )
                                                   >> 140   { 
                                                   >> 141     outStr << " G4CoupledTransportation: Statistics for looping particles "
                                                   >> 142            << G4endl;
                                                   >> 143     outStr << "   Sum of energy of loopers killed: "
                                                   >> 144            <<  fSumEnergyKilled / MeV << " MeV " << G4endl;
                                                   >> 145     outStr << "   Max energy of loopers killed: "
                                                   >> 146            <<  fMaxEnergyKilled / MeV << " MeV " << G4endl;
                                                   >> 147 
                                                   >> 148 
                                                   >> 149     outStr << "   Max energy of loopers 'saved':  " << fMaxEnergySaved << G4endl;
                                                   >> 150     outStr << "   Sum of energy of loopers 'saved': "
                                                   >> 151            <<  fSumEnergySaved << G4endl;
                                                   >> 152     outStr << "   Sum of energy of unstable loopers 'saved': "
                                                   >> 153            << fSumEnergyUnstableSaved << G4endl;
                                                   >> 154   } 
 86 }                                                 155 }
 87                                                   156 
 88 //////////////////////////////////////////////    157 //////////////////////////////////////////////////////////////////////////
 89 //                                                158 //
 90 // Responsibilities:                              159 // Responsibilities:
 91 //    Find whether the geometry limits the Ste    160 //    Find whether the geometry limits the Step, and to what length
 92 //    Calculate the new value of the safety an    161 //    Calculate the new value of the safety and return it.
 93 //    Store the final time, position and momen    162 //    Store the final time, position and momentum.
 94                                                   163 
 95 G4double G4CoupledTransportation::                164 G4double G4CoupledTransportation::
 96 AlongStepGetPhysicalInteractionLength( const G    165 AlongStepGetPhysicalInteractionLength( const G4Track&  track,
 97                                              G    166                                              G4double, //  previousStepSize
 98                                              G    167                                              G4double  currentMinimumStep,
 99                                              G    168                                              G4double& proposedSafetyForStart,
100                                              G    169                                              G4GPILSelection* selection )
101 {                                                 170 {
102   G4double geometryStepLength;                    171   G4double geometryStepLength; 
                                                   >> 172   G4double startMassSafety= 0.0; // estimated safety for start point (mass geometry)
103   G4double startFullSafety= 0.0; // estimated     173   G4double startFullSafety= 0.0; // estimated safety for start point (all geometries)
104   G4double safetyProposal= -1.0; // local copy    174   G4double safetyProposal= -1.0; // local copy of proposal 
105                                                   175 
106   G4ThreeVector  EndUnitMomentum ;                176   G4ThreeVector  EndUnitMomentum ;
107   G4double       lengthAlongCurve = 0.0 ;         177   G4double       lengthAlongCurve = 0.0 ;
108                                                   178  
109   fParticleIsLooping = false ;                    179   fParticleIsLooping = false ;
110                                                   180 
111   // Initial actions moved to  StartTrack()       181   // Initial actions moved to  StartTrack()   
112   // --------------------------------------       182   // --------------------------------------
113   // Note: in case another process changes tou    183   // Note: in case another process changes touchable handle
114   //    it will be necessary to add here (for     184   //    it will be necessary to add here (for all steps)   
115   // fCurrentTouchableHandle = aTrack->GetTouc    185   // fCurrentTouchableHandle = aTrack->GetTouchableHandle();
116                                                   186 
117   // GPILSelection is set to defaule value of     187   // GPILSelection is set to defaule value of CandidateForSelection
118   // It is a return value                         188   // It is a return value
119   //                                              189   //
120   *selection = CandidateForSelection ;            190   *selection = CandidateForSelection ;
121                                                   191 
122   fFirstStepInMassVolume = fNewTrack || fMassG << 192   fFirstStepInMassVolume = fNewTrack | fMassGeometryLimitedStep ; 
123   fFirstStepInVolume     = fNewTrack || fGeome << 193   fFirstStepInAnyVolume =  fNewTrack | fAnyGeometryLimitedStep ;
124                                                   194 
125 #ifdef G4DEBUG_TRANSPORT                          195 #ifdef G4DEBUG_TRANSPORT
126   G4cout << "  CoupledTransport::AlongStep GPI    196   G4cout << "  CoupledTransport::AlongStep GPIL:  "
127          << "  1st-step:  any= "  <<fFirstStep << 197          << "  1st-step:  any= "  <<fFirstStepInAnyVolume  << " ( geom= "
128          << fGeometryLimitedStep << " ) "      << 198          << fAnyGeometryLimitedStep << " ) "
129          <<           " mass= " << fFirstStepI    199          <<           " mass= " << fFirstStepInMassVolume << " ( geom= "
130          << fMassGeometryLimitedStep << " ) "     200          << fMassGeometryLimitedStep << " ) " 
131          << "  newTrack= " << fNewTrack << G4e    201          << "  newTrack= " << fNewTrack << G4endl;
132 #endif                                            202 #endif
133                                                   203   
134   // fLastStepInVolume= false;                    204   // fLastStepInVolume= false;
135   fNewTrack = false;                              205   fNewTrack = false;
136                                                   206 
137   // Get initial Energy/Momentum of the track     207   // Get initial Energy/Momentum of the track
138   //                                              208   //
139   const G4DynamicParticle*    pParticle  = tra    209   const G4DynamicParticle*    pParticle  = track.GetDynamicParticle() ;
140   const G4ParticleDefinition* pParticleDef   =    210   const G4ParticleDefinition* pParticleDef   = pParticle->GetDefinition() ;
141   G4ThreeVector startMomentumDir       = pPart    211   G4ThreeVector startMomentumDir       = pParticle->GetMomentumDirection() ;
142   G4ThreeVector startPosition          = track    212   G4ThreeVector startPosition          = track.GetPosition() ;
143   G4VPhysicalVolume* currentVolume= track.GetV    213   G4VPhysicalVolume* currentVolume= track.GetVolume(); 
144                                                   214 
145 #ifdef G4DEBUG_TRANSPORT                          215 #ifdef G4DEBUG_TRANSPORT
146   if( verboseLevel > 1 )                          216   if( verboseLevel > 1 )
147   {                                               217   {
148     G4cout << "G4CoupledTransportation::AlongS    218     G4cout << "G4CoupledTransportation::AlongStepGPIL> called in volume " 
149            << currentVolume->GetName() << G4en    219            << currentVolume->GetName() << G4endl; 
150   }                                               220   }
151 #endif                                            221 #endif
152   // G4double   theTime        = track.GetGlob    222   // G4double   theTime        = track.GetGlobalTime() ;
153                                                   223 
154   // The Step Point safety can be limited by o    224   // The Step Point safety can be limited by other geometries and/or the 
155   // assumptions of any process - it's not alw    225   // assumptions of any process - it's not always the geometrical safety.
156   // We calculate the starting point's isotrop    226   // We calculate the starting point's isotropic safety here.
157   //                                              227   //
158   G4ThreeVector OriginShift = startPosition -     228   G4ThreeVector OriginShift = startPosition - fPreviousSftOrigin ;
159   G4double      MagSqShift  = OriginShift.mag2    229   G4double      MagSqShift  = OriginShift.mag2() ;
                                                   >> 230   startMassSafety = 0.0; 
160   startFullSafety= 0.0;                           231   startFullSafety= 0.0; 
161                                                   232 
162   // Recall that FullSafety <= MassSafety         233   // Recall that FullSafety <= MassSafety 
163   // Original: if( MagSqShift < sqr(fPreviousM    234   // Original: if( MagSqShift < sqr(fPreviousMassSafety) ) {
164   if( MagSqShift < sqr(fPreviousFullSafety) )     235   if( MagSqShift < sqr(fPreviousFullSafety) )
165   {                                               236   {
166      G4double mag_shift= std::sqrt(MagSqShift)    237      G4double mag_shift= std::sqrt(MagSqShift); 
                                                   >> 238      startMassSafety = std::max( (fPreviousMassSafety - mag_shift), 0.0); 
167      startFullSafety = std::max( (fPreviousFul    239      startFullSafety = std::max( (fPreviousFullSafety - mag_shift), 0.0);
168        // Need to be consistent between full s    240        // Need to be consistent between full safety with Mass safety
169        // in order reproduce results in simple    241        // in order reproduce results in simple case
170        // --> use same calculation method         242        // --> use same calculation method
171                                                   243 
172      // Only compute full safety if massSafety    244      // Only compute full safety if massSafety > 0.  Else it remains 0
173      // startFullSafety = fPathFinder->Compute    245      // startFullSafety = fPathFinder->ComputeSafety( startPosition ); 
174   }                                               246   }
175                                                   247 
176   // Is the particle charged or has it a magne    248   // Is the particle charged or has it a magnetic moment?
177   //                                              249   //
178   G4double particleCharge = pParticle->GetChar    250   G4double particleCharge = pParticle->GetCharge() ;
179   G4double magneticMoment = pParticle->GetMagn    251   G4double magneticMoment = pParticle->GetMagneticMoment() ;
180   G4double       restMass = pParticle->GetMass    252   G4double       restMass = pParticle->GetMass() ; 
181                                                   253 
182   fMassGeometryLimitedStep = false ; //  Set d    254   fMassGeometryLimitedStep = false ; //  Set default - alex
183   fGeometryLimitedStep     = false;            << 255   fAnyGeometryLimitedStep = false; 
184                                                   256 
185   // There is no need to locate the current vo    257   // There is no need to locate the current volume. It is Done elsewhere:
186   //   On track construction                      258   //   On track construction 
187   //   By the tracking, after all AlongStepDoI    259   //   By the tracking, after all AlongStepDoIts, in "Relocation"
188                                                   260 
189   // Check if the particle has a force, EM or     261   // Check if the particle has a force, EM or gravitational, exerted on it
190   //                                              262   //
191   G4FieldManager* fieldMgr= nullptr;              263   G4FieldManager* fieldMgr= nullptr;
192   G4bool          fieldExertsForce = false ;      264   G4bool          fieldExertsForce = false ;
193                                                   265 
194   const G4Field* ptrField= nullptr;               266   const G4Field* ptrField= nullptr;
195                                                   267 
196   fieldMgr = fFieldPropagator->FindAndSetField    268   fieldMgr = fFieldPropagator->FindAndSetFieldManager( track.GetVolume() );
197   G4bool eligibleEM = (particleCharge != 0.0)     269   G4bool eligibleEM = (particleCharge != 0.0)
198                    || ( fUseMagneticMoment &&     270                    || ( fUseMagneticMoment && (magneticMoment != 0.0) );
199   G4bool eligibleGrav =  fUseGravity && (restM    271   G4bool eligibleGrav =  fUseGravity && (restMass != 0.0) ;
200                                                   272 
201   if( (fieldMgr!=nullptr) && (eligibleEM||elig    273   if( (fieldMgr!=nullptr) && (eligibleEM||eligibleGrav) )
202   {                                               274   {
203      // Message the field Manager, to configur    275      // Message the field Manager, to configure it for this track
204      //                                           276      //
205      fieldMgr->ConfigureForTrack( &track );       277      fieldMgr->ConfigureForTrack( &track );
206                                                   278 
207      // The above call can transition from a n    279      // The above call can transition from a null field-ptr oto a finite field.
208      // If the field manager has no field ptr,    280      // If the field manager has no field ptr, the field is zero 
209      // by definition ( = there is no field !     281      // by definition ( = there is no field ! )
210      //                                           282      //
211      ptrField= fieldMgr->GetDetectorField();      283      ptrField= fieldMgr->GetDetectorField();
212                                                   284  
213      if( ptrField != nullptr)                     285      if( ptrField != nullptr)
214      {                                            286      { 
215         fieldExertsForce = eligibleEM             287         fieldExertsForce = eligibleEM
216               || ( eligibleGrav && ptrField->I    288               || ( eligibleGrav && ptrField->IsGravityActive() );
217      }                                            289      }
218   }                                               290   }
219   G4double momentumMagnitude = pParticle->GetT    291   G4double momentumMagnitude = pParticle->GetTotalMomentum() ;
220                                                   292 
221   if( fieldExertsForce )                          293   if( fieldExertsForce )
222   {                                               294   {
223      auto equationOfMotion= fFieldPropagator->    295      auto equationOfMotion= fFieldPropagator->GetCurrentEquationOfMotion();
224                                                   296  
225      G4ChargeState chargeState(particleCharge,    297      G4ChargeState chargeState(particleCharge, // Charge can change (dynamic)
226                                magneticMoment,    298                                magneticMoment,
227                                pParticleDef->G    299                                pParticleDef->GetPDGSpin() );
228      if( equationOfMotion )                       300      if( equationOfMotion )
229      {                                            301      {
230         equationOfMotion->SetChargeMomentumMas    302         equationOfMotion->SetChargeMomentumMass( chargeState,
231                                                   303                                                  momentumMagnitude,
232                                                   304                                                  restMass );
233      }                                            305      }
234 #ifdef G4DEBUG_TRANSPORT                          306 #ifdef G4DEBUG_TRANSPORT
235      else                                         307      else
236      {                                            308      {
237         G4cerr << " ERROR in G4CoupledTranspor    309         G4cerr << " ERROR in G4CoupledTransportation> "
238                << "Cannot find valid Equation     310                << "Cannot find valid Equation of motion: "      << G4endl;
239                << " Unable to pass Charge, Mom    311                << " Unable to pass Charge, Momentum and Mass "  << G4endl;
240      }                                            312      }
241 #endif                                            313 #endif     
242   }                                               314   }
243                                                   315 
244   G4ThreeVector polarizationVec  = track.GetPo    316   G4ThreeVector polarizationVec  = track.GetPolarization() ;
245   G4FieldTrack  aFieldTrack = G4FieldTrack(sta    317   G4FieldTrack  aFieldTrack = G4FieldTrack(startPosition, 
246                                            tra    318                                            track.GetGlobalTime(), // Lab.
247                                            tra    319                                            track.GetMomentumDirection(),
248                                            tra    320                                            track.GetKineticEnergy(),
249                                            res    321                                            restMass,
250                                            par    322                                            particleCharge, 
251                                            pol    323                                            polarizationVec, 
252                                            pPa    324                                            pParticleDef->GetPDGMagneticMoment(),
253                                            0.0    325                                            0.0,  // Length along track
254                                            pPa    326                                            pParticleDef->GetPDGSpin() ) ;
255   G4int stepNo= track.GetCurrentStepNumber();     327   G4int stepNo= track.GetCurrentStepNumber(); 
256                                                   328 
257   ELimited limitedStep;                           329   ELimited limitedStep; 
258   G4FieldTrack endTrackState('a');  //  Defaul    330   G4FieldTrack endTrackState('a');  //  Default values
259                                                   331 
260   fMassGeometryLimitedStep = false ;    //  de << 332   fMassGeometryLimitedStep = false ;    //  default 
261   fGeometryLimitedStep     = false;            << 333   fAnyGeometryLimitedStep  = false ;
262   if( currentMinimumStep > 0 )                    334   if( currentMinimumStep > 0 )
263   {                                               335   {
264       G4double newMassSafety= 0.0;     //  tem    336       G4double newMassSafety= 0.0;     //  temp. for recalculation
265                                                   337 
266       // Do the Transport in the field (non re    338       // Do the Transport in the field (non recti-linear)
267       //                                          339       //
268       lengthAlongCurve = fPathFinder->ComputeS    340       lengthAlongCurve = fPathFinder->ComputeStep( aFieldTrack,
269                                                   341                                                    currentMinimumStep, 
270                                                << 342                                                    fNavigatorId,
271                                                   343                                                    stepNo,
272                                                   344                                                    newMassSafety,
273                                                   345                                                    limitedStep,
274                                                   346                                                    endTrackState,
275                                                   347                                                    currentVolume ) ;
276                                                   348 
277       G4double newFullSafety= fPathFinder->Get    349       G4double newFullSafety= fPathFinder->GetCurrentSafety();  
278         // this was estimated already in step     350         // this was estimated already in step above
279                                                   351 
280       if( limitedStep == kUnique || limitedSte    352       if( limitedStep == kUnique || limitedStep == kSharedTransport )
281       {                                           353       {
282         fMassGeometryLimitedStep = true ;         354         fMassGeometryLimitedStep = true ;
283       }                                           355       }
284                                                   356       
285       fGeometryLimitedStep = (fPathFinder->Get << 357       fAnyGeometryLimitedStep = (fPathFinder->GetNumberGeometriesLimitingStep() != 0);
286                                                   358 
287 #ifdef G4DEBUG_TRANSPORT                          359 #ifdef G4DEBUG_TRANSPORT
288       if( fMassGeometryLimitedStep && !fGeomet << 360       if( fMassGeometryLimitedStep && !fAnyGeometryLimitedStep )
289       {                                           361       {
290         std::ostringstream message;               362         std::ostringstream message;
291         message << " ERROR in determining geom    363         message << " ERROR in determining geometries limiting the step" << G4endl;
292         message << "  Limiting:  mass=" << fMa    364         message << "  Limiting:  mass=" << fMassGeometryLimitedStep
293                 << " any= " << fGeometryLimite << 365                 << " any= " << fAnyGeometryLimitedStep << G4endl;
294         message << "Incompatible conditions -     366         message << "Incompatible conditions - by which geometries was it limited ?";
295         G4Exception("G4CoupledTransportation::    367         G4Exception("G4CoupledTransportation::AlongStepGetPhysicalInteractionLength()", 
296                     "PathFinderConfused", Fata    368                     "PathFinderConfused", FatalException, message); 
297       }                                           369       }
298 #endif                                            370 #endif
299                                                   371 
300       geometryStepLength = std::min( lengthAlo    372       geometryStepLength = std::min( lengthAlongCurve, currentMinimumStep); 
301                                                   373 
302       // Momentum:  Magnitude and direction ca    374       // Momentum:  Magnitude and direction can be changed too now ...
303       //                                          375       //
304       fMomentumChanged         = true ;           376       fMomentumChanged         = true ; 
305       fTransportEndMomentumDir = endTrackState    377       fTransportEndMomentumDir = endTrackState.GetMomentumDir() ;
306                                                   378 
307       // Remember last safety origin & value.     379       // Remember last safety origin & value.
308       fPreviousSftOrigin  = startPosition ;       380       fPreviousSftOrigin  = startPosition ;
309       fPreviousMassSafety = newMassSafety ;       381       fPreviousMassSafety = newMassSafety ;         
310       fPreviousFullSafety = newFullSafety ;       382       fPreviousFullSafety = newFullSafety ; 
311       // fpSafetyHelper->SetCurrentSafety( new    383       // fpSafetyHelper->SetCurrentSafety( newFullSafety, startPosition);
312                                                   384 
313 #ifdef G4DEBUG_TRANSPORT                          385 #ifdef G4DEBUG_TRANSPORT
314       if( verboseLevel > 1 )                      386       if( verboseLevel > 1 )
315       {                                           387       {
316         G4cout << "G4Transport:CompStep> "        388         G4cout << "G4Transport:CompStep> " 
317                << " called the pathfinder for     389                << " called the pathfinder for a new step at " << startPosition
318                << " and obtained step = " << l    390                << " and obtained step = " << lengthAlongCurve << G4endl;
319         G4cout << "  New safety (preStep) = "  << 391         G4cout << "  New safety (preStep) = " << newMassSafety 
                                                   >> 392                << " versus precalculated = "  << startMassSafety << G4endl; 
320       }                                           393       }
321 #endif                                            394 #endif
322                                                   395 
323       // Store as best estimate value             396       // Store as best estimate value
                                                   >> 397       startMassSafety    = newMassSafety ; 
324       startFullSafety    = newFullSafety ;        398       startFullSafety    = newFullSafety ; 
325                                                   399 
326       // Get the End-Position and End-Momentum    400       // Get the End-Position and End-Momentum (Dir-ection)
327       fTransportEndPosition = endTrackState.Ge    401       fTransportEndPosition = endTrackState.GetPosition() ;
328       fTransportEndKineticEnergy  = endTrackSt    402       fTransportEndKineticEnergy  = endTrackState.GetKineticEnergy() ; 
329   }                                               403   }
330   else                                            404   else
331   {                                               405   {
332       geometryStepLength   = lengthAlongCurve=    406       geometryStepLength   = lengthAlongCurve= 0.0 ;
333       fMomentumChanged         = false ;          407       fMomentumChanged         = false ; 
334       // fMassGeometryLimitedStep = false ;       408       // fMassGeometryLimitedStep = false ;   //  --- ???
335       // fGeometryLimitedStep = true;          << 409       // fAnyGeometryLimitedStep = true;
336       fTransportEndMomentumDir = track.GetMome    410       fTransportEndMomentumDir = track.GetMomentumDirection();
337       fTransportEndKineticEnergy  = track.GetK    411       fTransportEndKineticEnergy  = track.GetKineticEnergy();
338                                                   412 
339       fTransportEndPosition = startPosition;      413       fTransportEndPosition = startPosition;
340                                                   414 
341       endTrackState= aFieldTrack;  // Ensures     415       endTrackState= aFieldTrack;  // Ensures that time is updated
                                                   >> 416 
                                                   >> 417       // If the step length requested is 0, and we are on a boundary
                                                   >> 418       //   then a boundary will also limit the step.
                                                   >> 419       if( startMassSafety == 0.0 )
                                                   >> 420       {
                                                   >> 421          fMassGeometryLimitedStep = true ;
                                                   >> 422          fAnyGeometryLimitedStep = true;
                                                   >> 423       }
                                                   >> 424       //   TODO:  Add explicit logical status for being at a boundary
342   }                                               425   }
343   // G4FieldTrack aTrackState(endTrackState);     426   // G4FieldTrack aTrackState(endTrackState);  
344                                                   427 
345   if( !fieldExertsForce )                         428   if( !fieldExertsForce ) 
346   {                                               429   { 
347       fParticleIsLooping         = false ;        430       fParticleIsLooping         = false ; 
348       fMomentumChanged           = false ;        431       fMomentumChanged           = false ; 
349       fEndGlobalTimeComputed     = false ;        432       fEndGlobalTimeComputed     = false ; 
350   }                                               433   } 
351   else                                            434   else 
352   {                                               435   { 
353       fParticleIsLooping = fFieldPropagator->I    436       fParticleIsLooping = fFieldPropagator->IsParticleLooping() ;
354                                                   437   
355 #ifdef G4DEBUG_TRANSPORT                          438 #ifdef G4DEBUG_TRANSPORT
356       if( verboseLevel > 1 )                      439       if( verboseLevel > 1 )
357       {                                           440       {
358         G4cout << " G4CT::CS End Position = "     441         G4cout << " G4CT::CS End Position = "
359                << fTransportEndPosition << G4e    442                << fTransportEndPosition << G4endl; 
360         G4cout << " G4CT::CS End Direction = "    443         G4cout << " G4CT::CS End Direction = "
361                << fTransportEndMomentumDir <<     444                << fTransportEndMomentumDir << G4endl; 
362       }                                           445       }
363 #endif                                            446 #endif
364       if( fFieldPropagator->GetCurrentFieldMan    447       if( fFieldPropagator->GetCurrentFieldManager()->DoesFieldChangeEnergy() )
365       {                                           448       {
366           // If the field can change energy, t    449           // If the field can change energy, then the time must be integrated
367           //    - so this should have been upd    450           //    - so this should have been updated
368           //                                      451           //
369           fCandidateEndGlobalTime   = endTrack    452           fCandidateEndGlobalTime   = endTrackState.GetLabTimeOfFlight();
370           fEndGlobalTimeComputed    = true;       453           fEndGlobalTimeComputed    = true;
371                                                   454   
372           // was ( fCandidateEndGlobalTime !=     455           // was ( fCandidateEndGlobalTime != track.GetGlobalTime() );
373           // a cleaner way is to have FieldTra    456           // a cleaner way is to have FieldTrack knowing whether time
374           // is updated                           457           // is updated
375       }                                           458       }
376       else                                        459       else
377       {                                           460       {
378           // The energy should be unchanged by    461           // The energy should be unchanged by field transport,
379           //    - so the time changed will be     462           //    - so the time changed will be calculated elsewhere
380           //                                      463           //
381           fEndGlobalTimeComputed = false;         464           fEndGlobalTimeComputed = false;
382                                                   465   
383 #ifdef G4VERBOSE                               << 
384           // Check that the integration preser    466           // Check that the integration preserved the energy 
385           //     -  and if not correct this!      467           //     -  and if not correct this!
386           G4double  startEnergy= track.GetKine    468           G4double  startEnergy= track.GetKineticEnergy();
387           G4double  endEnergy= fTransportEndKi    469           G4double  endEnergy= fTransportEndKineticEnergy; 
388                                                   470       
                                                   >> 471           static G4ThreadLocal G4int no_inexact_steps=0; // , no_large_ediff;
389           G4double absEdiff = std::fabs(startE    472           G4double absEdiff = std::fabs(startEnergy- endEnergy);
                                                   >> 473           if( absEdiff > perMillion * endEnergy )
                                                   >> 474           {
                                                   >> 475             no_inexact_steps++;
                                                   >> 476             // Possible statistics keeping here ...
                                                   >> 477           }
                                                   >> 478 #ifdef G4VERBOSE
390           if( (verboseLevel > 1) && ( absEdiff    479           if( (verboseLevel > 1) && ( absEdiff > perThousand * endEnergy) )
391           {                                       480           {
392             ReportInexactEnergy(startEnergy, e    481             ReportInexactEnergy(startEnergy, endEnergy); 
393           }  // end of if (verboseLevel)          482           }  // end of if (verboseLevel)
394 #endif                                            483 #endif
395           // Correct the energy for fields tha    484           // Correct the energy for fields that conserve it
396           //  This - hides the integration err    485           //  This - hides the integration error
397           //       - but gives a better physic    486           //       - but gives a better physical answer
398           fTransportEndKineticEnergy= track.Ge    487           fTransportEndKineticEnergy= track.GetKineticEnergy();
399       }                                           488       }
400   }                                               489   }
401                                                   490 
402   fEndPointDistance   = (fTransportEndPosition << 491   fEndpointDistance   = (fTransportEndPosition - startPosition).mag() ;
403   fTransportEndSpin = endTrackState.GetSpin();    492   fTransportEndSpin = endTrackState.GetSpin();
404                                                   493 
405   // Calculate the safety                         494   // Calculate the safety
406                                                   495 
407   safetyProposal= startFullSafety;   // used t    496   safetyProposal= startFullSafety;   // used to be startMassSafety
408     // Changed to accomodate processes that ca    497     // Changed to accomodate processes that cannot update the safety
409                                                   498 
410   // Update safety for the end-point, if becom    499   // Update safety for the end-point, if becomes negative at the end-point.
411                                                   500 
412   if(   (startFullSafety < fEndPointDistance ) << 501   if(   (startFullSafety < fEndpointDistance ) 
413         && ( particleCharge != 0.0 ) )  // Onl    502         && ( particleCharge != 0.0 ) )  // Only needed to prepare for MSC
414    //   && !fGeometryLimitedStep ) // To-Try:  << 503    //   && !fAnyGeometryLimitedStep ) // To-Try: No safety update if at boundary
415   {                                               504   {
416       G4double endFullSafety =                    505       G4double endFullSafety =
417         fPathFinder->ComputeSafety( fTransport    506         fPathFinder->ComputeSafety( fTransportEndPosition); 
418         // Expected mission -- only mass geome    507         // Expected mission -- only mass geometry's safety
419         //   fLinearNavigator->ComputeSafety(  << 508         //   fMassNavigator->ComputeSafety( fTransportEndPosition) ;
420         // Yet discrete processes only have po    509         // Yet discrete processes only have poststep -- and this cannot 
421         //   currently revise the safety          510         //   currently revise the safety  
422         //   ==> so we use the all-geometry sa    511         //   ==> so we use the all-geometry safety as a precaution
423                                                   512 
424       fpSafetyHelper->SetCurrentSafety( endFul    513       fpSafetyHelper->SetCurrentSafety( endFullSafety, fTransportEndPosition);
425         // Pushing safety to Helper avoids rec    514         // Pushing safety to Helper avoids recalculation at this point
426                                                   515 
427       G4ThreeVector centerPt= G4ThreeVector(0.    516       G4ThreeVector centerPt= G4ThreeVector(0.0, 0.0, 0.0);  // Used for return value
428       G4double endMassSafety= fPathFinder->Obt << 517       G4double endMassSafety= fPathFinder->ObtainSafety( fNavigatorId, centerPt); 
429         //  Retrieves the mass value from Path    518         //  Retrieves the mass value from PathFinder (it calculated it)
430                                                   519 
431       fPreviousMassSafety = endMassSafety ;       520       fPreviousMassSafety = endMassSafety ; 
432       fPreviousFullSafety = endFullSafety;        521       fPreviousFullSafety = endFullSafety; 
433       fPreviousSftOrigin = fTransportEndPositi    522       fPreviousSftOrigin = fTransportEndPosition ;
434                                                   523 
435       // The convention (Stepping Manager's) i    524       // The convention (Stepping Manager's) is safety from the start point
436       //                                          525       //
437       safetyProposal = endFullSafety + fEndPoi << 526       safetyProposal = endFullSafety + fEndpointDistance;
438           //  --> was endMassSafety               527           //  --> was endMassSafety
439       // Changed to accomodate processes that     528       // Changed to accomodate processes that cannot update the safety
440                                                   529 
441 #ifdef G4DEBUG_TRANSPORT                          530 #ifdef G4DEBUG_TRANSPORT 
442       G4int prec= G4cout.precision(12) ;          531       G4int prec= G4cout.precision(12) ;
443       G4cout << "***CoupledTransportation::Alo    532       G4cout << "***CoupledTransportation::AlongStepGPIL ** " << G4endl  ;
444       G4cout << "  Revised Safety at endpoint     533       G4cout << "  Revised Safety at endpoint "  << fTransportEndPosition
445              << "   give safety values: Mass=     534              << "   give safety values: Mass= " << endMassSafety 
446              << "  All= " << endFullSafety <<     535              << "  All= " << endFullSafety << G4endl ; 
447       G4cout << "  Adding endpoint distance "  << 536       G4cout << "  Adding endpoint distance " << fEndpointDistance 
448              << "   to obtain pseudo-safety= "    537              << "   to obtain pseudo-safety= " << safetyProposal << G4endl ; 
449       G4cout.precision(prec);                     538       G4cout.precision(prec); 
450   }                                               539   }  
451   else                                            540   else
452   {                                               541   {
453       G4int prec= G4cout.precision(12) ;          542       G4int prec= G4cout.precision(12) ;
454       G4cout << "***CoupledTransportation::Alo    543       G4cout << "***CoupledTransportation::AlongStepGPIL ** " << G4endl  ;
455       G4cout << "  Quick Safety estimate at en    544       G4cout << "  Quick Safety estimate at endpoint "
456              << fTransportEndPosition             545              << fTransportEndPosition
457              << "   gives safety endpoint valu    546              << "   gives safety endpoint value = "
458              << startFullSafety - fEndPointDis << 547              << startFullSafety - fEndpointDistance
459              << "  using start-point value " <    548              << "  using start-point value " << startFullSafety 
460              << "  and endpointDistance " << f << 549              << "  and endpointDistance " << fEndpointDistance << G4endl; 
461       G4cout.precision(prec);                     550       G4cout.precision(prec); 
462 #endif                                            551 #endif
463   }                                               552   }          
464                                                   553 
465   proposedSafetyForStart= safetyProposal;         554   proposedSafetyForStart= safetyProposal; 
466   fParticleChange.ProposeTrueStepLength(geomet    555   fParticleChange.ProposeTrueStepLength(geometryStepLength) ;
467                                                   556 
468   return geometryStepLength ;                     557   return geometryStepLength ;
469 }                                                 558 }
470                                                   559 
                                                   >> 560 //////////////////////////////////////////////////////////////////////////
                                                   >> 561 
                                                   >> 562 G4VParticleChange*
                                                   >> 563 G4CoupledTransportation::AlongStepDoIt( const G4Track& track,
                                                   >> 564                                         const G4Step&  stepData )
                                                   >> 565 {
                                                   >> 566   static G4ThreadLocal G4long noCallsCT_ASDI=0;
                                                   >> 567   const char *methodName= "AlongStepDoIt";
                                                   >> 568   
                                                   >> 569   noCallsCT_ASDI++;
                                                   >> 570 
                                                   >> 571   fParticleChange.Initialize(track) ;
                                                   >> 572     // sets all its members to the value of corresponding members in G4Track
                                                   >> 573 
                                                   >> 574   //  Code specific for Transport
                                                   >> 575   //
                                                   >> 576   fParticleChange.ProposePosition(fTransportEndPosition) ;
                                                   >> 577   fParticleChange.ProposeMomentumDirection(fTransportEndMomentumDir) ;
                                                   >> 578   fParticleChange.ProposeEnergy(fTransportEndKineticEnergy) ;
                                                   >> 579   fParticleChange.SetMomentumChanged(fMomentumChanged) ;
                                                   >> 580 
                                                   >> 581   fParticleChange.ProposePolarization(fTransportEndSpin);
                                                   >> 582   
                                                   >> 583   G4double deltaTime = 0.0 ;
                                                   >> 584 
                                                   >> 585   // Calculate  Lab Time of Flight (ONLY if field Equations used it!)
                                                   >> 586      // G4double endTime   = fCandidateEndGlobalTime;
                                                   >> 587      // G4double delta_time = endTime - startTime;
                                                   >> 588 
                                                   >> 589   G4double startTime = track.GetGlobalTime() ;
                                                   >> 590   
                                                   >> 591   if (!fEndGlobalTimeComputed)
                                                   >> 592   {
                                                   >> 593      G4double finalInverseVel= DBL_MAX, initialInverseVel=DBL_MAX; 
                                                   >> 594 
                                                   >> 595      // The time was not integrated .. make the best estimate possible
                                                   >> 596      //
                                                   >> 597      G4double finalVelocity   = track.GetVelocity() ;
                                                   >> 598      if( finalVelocity > 0.0 ) { finalInverseVel= 1.0 / finalVelocity; }
                                                   >> 599      G4double initialVelocity = stepData.GetPreStepPoint()->GetVelocity() ;
                                                   >> 600      if( initialVelocity > 0.0 ) { initialInverseVel= 1.0 / initialVelocity; }
                                                   >> 601      G4double stepLength      = track.GetStepLength() ;
                                                   >> 602 
                                                   >> 603      if (finalVelocity > 0.0)
                                                   >> 604      {
                                                   >> 605         // deltaTime = stepLength/finalVelocity ;
                                                   >> 606         G4double meanInverseVelocity = 0.5
                                                   >> 607                                      * (initialInverseVel + finalInverseVel);
                                                   >> 608         deltaTime = stepLength * meanInverseVelocity ;
                                                   >> 609      }
                                                   >> 610      else
                                                   >> 611      {
                                                   >> 612         deltaTime = stepLength * initialInverseVel ;
                                                   >> 613       }  //  Could do with better estimate for final step (finalVelocity = 0) ?
                                                   >> 614 
                                                   >> 615      fCandidateEndGlobalTime   = startTime + deltaTime ;
                                                   >> 616      fParticleChange.ProposeLocalTime(  track.GetLocalTime() + deltaTime) ;
                                                   >> 617   }
                                                   >> 618   else
                                                   >> 619   {
                                                   >> 620      deltaTime = fCandidateEndGlobalTime - startTime ;
                                                   >> 621      fParticleChange.ProposeGlobalTime( fCandidateEndGlobalTime ) ;
                                                   >> 622   }
                                                   >> 623 
                                                   >> 624   // Now Correct by Lorentz factor to get "proper" deltaTime
                                                   >> 625   
                                                   >> 626   G4double  restMass       = track.GetDynamicParticle()->GetMass() ;
                                                   >> 627   G4double deltaProperTime = deltaTime*( restMass/track.GetTotalEnergy() ) ;
                                                   >> 628 
                                                   >> 629   fParticleChange.ProposeProperTime(track.GetProperTime() + deltaProperTime) ;
                                                   >> 630   // fParticleChange. ProposeTrueStepLength( track.GetStepLength() ) ;
                                                   >> 631 
                                                   >> 632   // If the particle is caught looping or is stuck (in very difficult
                                                   >> 633   // boundaries) in a magnetic field (doing many steps) THEN this kills it ...
                                                   >> 634   //
                                                   >> 635   if ( fParticleIsLooping )
                                                   >> 636   {
                                                   >> 637      G4double endEnergy= fTransportEndKineticEnergy;
                                                   >> 638 
                                                   >> 639      const G4ParticleDefinition* particleType= 
                                                   >> 640          track.GetDynamicParticle() -> GetParticleDefinition();
                                                   >> 641      G4bool stable = particleType->GetPDGStable();
                                                   >> 642      
                                                   >> 643      G4bool candidateForEnd = (endEnergy < fThreshold_Important_Energy) 
                                                   >> 644                            || (fNoLooperTrials >= fThresholdTrials) ;
                                                   >> 645                                  
                                                   >> 646      if( candidateForEnd && stable )
                                                   >> 647      {
                                                   >> 648         const G4int electronPDG= 11; // G4Electron::G4Electron()->GetPDGEncoding();
                                                   >> 649         G4int particlePDG= particleType->GetPDGEncoding();
                                                   >> 650 
                                                   >> 651         // Kill the looping particle 
                                                   >> 652         //
                                                   >> 653         fParticleChange.ProposeTrackStatus( fStopAndKill )  ;
                                                   >> 654 
                                                   >> 655         // Simple statistics
                                                   >> 656         fSumEnergyKilled += endEnergy;
                                                   >> 657         fSumEnerSqKilled = endEnergy * endEnergy;
                                                   >> 658         fNumLoopersKilled++;
                                                   >> 659  
                                                   >> 660         if( endEnergy > fMaxEnergyKilled ) {
                                                   >> 661           fMaxEnergyKilled = endEnergy;
                                                   >> 662           fMaxEnergyKilledPDG = particlePDG; 
                                                   >> 663         }
                                                   >> 664         
                                                   >> 665         if(  particleType->GetPDGEncoding() != electronPDG )
                                                   >> 666         {
                                                   >> 667            fSumEnergyKilled_NonElectron += endEnergy;
                                                   >> 668            fSumEnerSqKilled_NonElectron += endEnergy * endEnergy;
                                                   >> 669            fNumLoopersKilled_NonElectron++;
                                                   >> 670            
                                                   >> 671            if( endEnergy > fMaxEnergyKilled_NonElectron )
                                                   >> 672            {
                                                   >> 673               fMaxEnergyKilled_NonElectron = endEnergy;
                                                   >> 674               fMaxEnergyKilled_NonElecPDG =  particlePDG;
                                                   >> 675            }
                                                   >> 676         }
                                                   >> 677 
                                                   >> 678         if( endEnergy > fThreshold_Warning_Energy && ! fSilenceLooperWarnings )           
                                                   >> 679         {
                                                   >> 680           fpLogger->ReportLoopingTrack( track, stepData, fNoLooperTrials,
                                                   >> 681                                         noCallsCT_ASDI, methodName );
                                                   >> 682         }
                                                   >> 683 
                                                   >> 684         fNoLooperTrials=0; 
                                                   >> 685       }
                                                   >> 686       else
                                                   >> 687       { 
                                                   >> 688         fNoLooperTrials ++;
                                                   >> 689 
                                                   >> 690         fMaxEnergySaved = std::max( endEnergy, fMaxEnergySaved);
                                                   >> 691         if( fNoLooperTrials == 1 ) {
                                                   >> 692           fSumEnergySaved += endEnergy;
                                                   >> 693           if ( !stable )
                                                   >> 694              fSumEnergyUnstableSaved += endEnergy;
                                                   >> 695         }
                                                   >> 696 #ifdef G4VERBOSE
                                                   >> 697         if( verboseLevel > 2 && ! fSilenceLooperWarnings )           
                                                   >> 698         {
                                                   >> 699           G4cout << "  ** G4CoupledTransportation::AlongStepDoIt():"
                                                   >> 700                  << " Particle is looping but is saved ..."  << G4endl
                                                   >> 701                  << "   Number of trials (this track) = " << fNoLooperTrials
                                                   >> 702                  << G4endl
                                                   >> 703                  << "   Steps by this track: " << track.GetCurrentStepNumber()
                                                   >> 704                  << G4endl
                                                   >> 705                  << "   Total no of calls to this method (all tracks) = "
                                                   >> 706                  << noCallsCT_ASDI << G4endl;
                                                   >> 707         }
                                                   >> 708 #endif
                                                   >> 709       }
                                                   >> 710   }
                                                   >> 711   else
                                                   >> 712   { 
                                                   >> 713       fNoLooperTrials=0; 
                                                   >> 714   }
                                                   >> 715 
                                                   >> 716   // Another (sometimes better way) is to use a user-limit maximum Step size
                                                   >> 717   // to alleviate this problem ..
                                                   >> 718   // Add smooth curved trajectories to particle-change
                                                   >> 719   //
                                                   >> 720   // fParticleChange.SetPointerToVectorOfAuxiliaryPoints
                                                   >> 721   //   (fFieldPropagator->GimmeTrajectoryVectorAndForgetIt() );
                                                   >> 722 
                                                   >> 723   return &fParticleChange ;
                                                   >> 724 }
                                                   >> 725 
                                                   >> 726 //////////////////////////////////////////////////////////////////////////
                                                   >> 727 //
                                                   >> 728 //  This ensures that the PostStep action is always called,
                                                   >> 729 //  so that it can do the relocation if it is needed.
                                                   >> 730 // 
                                                   >> 731 
                                                   >> 732 G4double G4CoupledTransportation::
                                                   >> 733 PostStepGetPhysicalInteractionLength( const G4Track&,
                                                   >> 734                                             G4double, // previousStepSize
                                                   >> 735                                             G4ForceCondition* pForceCond )
                                                   >> 736 { 
                                                   >> 737   // Must act as PostStep action -- to relocate particle
                                                   >> 738   *pForceCond = Forced ;    
                                                   >> 739   return DBL_MAX ;
                                                   >> 740 }
                                                   >> 741 
471 //////////////////////////////////////////////    742 /////////////////////////////////////////////////////////////////////////////
472                                                   743 
473 void G4CoupledTransportation::                    744 void G4CoupledTransportation::
474 ReportMove( G4ThreeVector OldVector, G4ThreeVe    745 ReportMove( G4ThreeVector OldVector, G4ThreeVector NewVector,
475             const G4String& Quantity )            746             const G4String& Quantity )
476 {                                                 747 {
477     G4ThreeVector moveVec = ( NewVector - OldV    748     G4ThreeVector moveVec = ( NewVector - OldVector );
478                                                   749 
479     G4cerr << G4endl                              750     G4cerr << G4endl
480            << "*******************************    751            << "**************************************************************"
481            << G4endl;                             752            << G4endl;
482     G4cerr << "Endpoint has moved between valu    753     G4cerr << "Endpoint has moved between value expected from TransportEndPosition "
483            << " and value from Track in PostSt    754            << " and value from Track in PostStepDoIt. " << G4endl
484            << "Change of " << Quantity << " is    755            << "Change of " << Quantity << " is " << moveVec.mag() / mm
485            << " mm long, "                        756            << " mm long, "
486            << " and its vector is " << (1.0/mm    757            << " and its vector is " << (1.0/mm) * moveVec << " mm " << G4endl
487            << "Endpoint of ComputeStep was " <    758            << "Endpoint of ComputeStep was " << OldVector
488            << " and current position to locate    759            << " and current position to locate is " << NewVector << G4endl;
489 }                                                 760 }
490                                                   761 
491 //////////////////////////////////////////////    762 /////////////////////////////////////////////////////////////////////////////
492                                                   763 
493 G4VParticleChange* G4CoupledTransportation::Po    764 G4VParticleChange* G4CoupledTransportation::PostStepDoIt( const G4Track& track,
494                                                   765                                                           const G4Step& )
495 {                                                 766 {
496   G4TouchableHandle retCurrentTouchable ;   //    767   G4TouchableHandle retCurrentTouchable ;   // The one to return
497                                                   768 
498   // Initialize ParticleChange  (by setting al    769   // Initialize ParticleChange  (by setting all its members equal
499   //                             to correspond    770   //                             to corresponding members in G4Track)
500   // fParticleChange.Initialize(track) ;  // T    771   // fParticleChange.Initialize(track) ;  // To initialise TouchableChange
501                                                   772 
502   fParticleChange.ProposeTrackStatus(track.Get    773   fParticleChange.ProposeTrackStatus(track.GetTrackStatus()) ;
503                                                   774 
504   if( fSignifyStepInAnyVolume )                   775   if( fSignifyStepInAnyVolume )
505   {                                               776   {
506      fParticleChange.ProposeFirstStepInVolume( << 777      fParticleChange.ProposeFirstStepInVolume( fFirstStepInAnyVolume );
507   }                                               778   }
508   else                                            779   else
509   {                                               780   {
510      fParticleChange.ProposeFirstStepInVolume(    781      fParticleChange.ProposeFirstStepInVolume( fFirstStepInMassVolume );
511   }                                               782   }
512                                                   783   
513   // Check that the end position and direction    784   // Check that the end position and direction are preserved 
514   // since call to AlongStepDoIt                  785   // since call to AlongStepDoIt
515                                                   786 
516 #ifdef G4DEBUG_TRANSPORT                          787 #ifdef G4DEBUG_TRANSPORT
517   if( ( verboseLevel > 0 )                        788   if( ( verboseLevel > 0 )
518      && ((fTransportEndPosition - track.GetPos    789      && ((fTransportEndPosition - track.GetPosition()).mag2() >= 1.0e-16) )
519   {                                               790   {
520      ReportMove( track.GetPosition(), fTranspo    791      ReportMove( track.GetPosition(), fTransportEndPosition,
521                  "End of Step Position" );        792                  "End of Step Position" ); 
522      G4cerr << " Problem in G4CoupledTransport    793      G4cerr << " Problem in G4CoupledTransportation::PostStepDoIt " << G4endl; 
523   }                                               794   }
524                                                   795 
525   // If the Step was determined by the volume     796   // If the Step was determined by the volume boundary, relocate the particle
526   // The pathFinder will know that the geometr    797   // The pathFinder will know that the geometry limited the step (!?)
527                                                   798 
528   if( verboseLevel > 0 )                          799   if( verboseLevel > 0 )
529   {                                               800   {
530      G4cout << " Calling PathFinder::Locate()     801      G4cout << " Calling PathFinder::Locate() from " 
531             << " G4CoupledTransportation::Post    802             << " G4CoupledTransportation::PostStepDoIt " << G4endl;
532      G4cout << "   fGeometryLimitedStep is " < << 803      G4cout << "  fAnyGeometryLimitedStep is " << fAnyGeometryLimitedStep
                                                   >> 804             << G4endl;
                                                   >> 805 
533   }                                               806   }
534 #endif                                            807 #endif
535                                                   808 
536   if(fGeometryLimitedStep)                     << 809   if(fAnyGeometryLimitedStep)
537   {                                               810   {  
538     fPathFinder->Locate( track.GetPosition(),     811     fPathFinder->Locate( track.GetPosition(), 
539                          track.GetMomentumDire    812                          track.GetMomentumDirection(),
540                          true);                   813                          true); 
541                                                   814 
542     // fCurrentTouchable will now become the p    815     // fCurrentTouchable will now become the previous touchable, 
543     // and what was the previous will be freed    816     // and what was the previous will be freed.
544     // (Needed because the preStepPoint can po    817     // (Needed because the preStepPoint can point to the previous touchable)
545                                                   818 
546     fCurrentTouchableHandle=                      819     fCurrentTouchableHandle= 
547       fPathFinder->CreateTouchableHandle( G4Tr << 820       fPathFinder->CreateTouchableHandle( fNavigatorId );
548                                                   821 
549 #ifdef G4DEBUG_TRANSPORT                          822 #ifdef G4DEBUG_TRANSPORT
                                                   >> 823     if( verboseLevel > 0 )
                                                   >> 824     {
                                                   >> 825       G4cout << "G4CoupledTransportation::PostStepDoIt --- fNavigatorId = " 
                                                   >> 826              << fNavigatorId << G4endl;
                                                   >> 827     }
550     if( verboseLevel > 1 )                        828     if( verboseLevel > 1 )
551     {                                             829     {
552        G4VPhysicalVolume* vol= fCurrentTouchab    830        G4VPhysicalVolume* vol= fCurrentTouchableHandle->GetVolume(); 
553        G4cout << "CHECK !!!!!!!!!!! fCurrentTo    831        G4cout << "CHECK !!!!!!!!!!! fCurrentTouchableHandle->GetVolume() = "
554               << vol;                             832               << vol;
555        if( vol ) { G4cout << "Name=" << vol->G    833        if( vol ) { G4cout << "Name=" << vol->GetName(); }
556        G4cout << G4endl;                          834        G4cout << G4endl;
557     }                                             835     }
558 #endif                                            836 #endif
559                                                   837 
560     // Check whether the particle is out of th    838     // Check whether the particle is out of the world volume 
561     // If so it has exited and must be killed.    839     // If so it has exited and must be killed.
562     //                                            840     //
563     if( fCurrentTouchableHandle->GetVolume() =    841     if( fCurrentTouchableHandle->GetVolume() == 0 )
564     {                                             842     {
565        fParticleChange.ProposeTrackStatus( fSt    843        fParticleChange.ProposeTrackStatus( fStopAndKill ) ;
566     }                                             844     }
567     retCurrentTouchable = fCurrentTouchableHan    845     retCurrentTouchable = fCurrentTouchableHandle ;
568     // fParticleChange.SetTouchableHandle( fCu    846     // fParticleChange.SetTouchableHandle( fCurrentTouchableHandle ) ;
569   }                                               847   }
570   else  // fGeometryLimitedStep  is false      << 848   else                 // fAnyGeometryLimitedStep  is false
571   {                                               849   { 
572 #ifdef G4DEBUG_TRANSPORT                          850 #ifdef G4DEBUG_TRANSPORT
573     if( verboseLevel > 1 )                        851     if( verboseLevel > 1 )
574     {                                             852     {
575       G4cout << "G4CoupledTransportation::Post << 853        G4cout << "G4CoupledTransportation::PostStepDoIt -- "
576              << " fGeometryLimitedStep  = " << << 854               << " fAnyGeometryLimitedStep  = " << fAnyGeometryLimitedStep  
577              << " must be false " << G4endl;   << 855               << " must be false " << G4endl;
578     }                                             856     }
579 #endif                                            857 #endif
580     // This serves only to move each of the Na    858     // This serves only to move each of the Navigator's location
581     //                                            859     //
582     // fLinearNavigator->LocateGlobalPointWith    860     // fLinearNavigator->LocateGlobalPointWithinVolume( track.GetPosition() ) ;
583                                                   861 
584     fPathFinder->ReLocate( track.GetPosition()    862     fPathFinder->ReLocate( track.GetPosition() );
585                            // track.GetMomentu    863                            // track.GetMomentumDirection() ); 
586                                                   864 
587     // Keep the value of the track's current T    865     // Keep the value of the track's current Touchable is retained,
588     //  and use it to overwrite the (unset) on    866     //  and use it to overwrite the (unset) one in particle change.
589     // Expect this must be fCurrentTouchable t    867     // Expect this must be fCurrentTouchable too
590     //   - could it be different, eg at the st    868     //   - could it be different, eg at the start of a step ?
591     //                                            869     //
592     retCurrentTouchable = track.GetTouchableHa    870     retCurrentTouchable = track.GetTouchableHandle() ;
593     // fParticleChange.SetTouchableHandle( tra    871     // fParticleChange.SetTouchableHandle( track.GetTouchableHandle() ) ;
594   }  // endif ( fGeometryLimitedStep )         << 872   }         // endif ( fAnyGeometryLimitedStep ) 
595                                                   873 
596 #ifdef G4DEBUG_NAVIGATION                         874 #ifdef G4DEBUG_NAVIGATION  
597   G4cout << "  CoupledTransport::AlongStep GPI    875   G4cout << "  CoupledTransport::AlongStep GPIL:  "
598          << " last-step:  any= " << fGeometryL << 876          << " last-step:  any= "  << fAnyGeometryLimitedStep
599          << " mass= " << fMassGeometryLimitedS << 877          << " . ..... x . " 
                                                   >> 878          <<            " mass= " <<  fMassGeometryLimitedStep
                                                   >> 879          << G4endl;
600 #endif                                            880 #endif
601                                                   881   
602   if( fSignifyStepInAnyVolume )                   882   if( fSignifyStepInAnyVolume )
603     fParticleChange.ProposeLastStepInVolume(fG << 883      fParticleChange.ProposeLastStepInVolume(fAnyGeometryLimitedStep);
604   else                                            884   else
605      fParticleChange.ProposeLastStepInVolume(f    885      fParticleChange.ProposeLastStepInVolume(fMassGeometryLimitedStep);
606                                                   886   
607   SetTouchableInformation(retCurrentTouchable) << 887   const G4VPhysicalVolume* pNewVol = retCurrentTouchable->GetVolume() ;
                                                   >> 888   const G4Material* pNewMaterial   = 0 ;
                                                   >> 889   const G4VSensitiveDetector* pNewSensitiveDetector   = 0 ;
                                                   >> 890                                                                                        
                                                   >> 891   if( pNewVol != 0 )
                                                   >> 892   {
                                                   >> 893     pNewMaterial= pNewVol->GetLogicalVolume()->GetMaterial();
                                                   >> 894     pNewSensitiveDetector= pNewVol->GetLogicalVolume()->GetSensitiveDetector();
                                                   >> 895   }
                                                   >> 896 
                                                   >> 897   // ( const_cast<G4Material *> pNewMaterial ) ;
                                                   >> 898   // ( const_cast<G4VSensitiveDetetor *> pNewSensitiveDetector) ;
                                                   >> 899 
                                                   >> 900   fParticleChange.SetMaterialInTouchable( (G4Material *) pNewMaterial ) ;
                                                   >> 901   fParticleChange.SetSensitiveDetectorInTouchable( (G4VSensitiveDetector *) pNewSensitiveDetector ) ;
                                                   >> 902     // "temporarily" until Get/Set Material of ParticleChange, 
                                                   >> 903     // and StepPoint can be made const. 
                                                   >> 904 
                                                   >> 905   const G4MaterialCutsCouple* pNewMaterialCutsCouple = 0;
                                                   >> 906   if( pNewVol != 0 )
                                                   >> 907   {
                                                   >> 908     pNewMaterialCutsCouple=pNewVol->GetLogicalVolume()->GetMaterialCutsCouple();
                                                   >> 909     if( pNewMaterialCutsCouple!=0 
                                                   >> 910         && pNewMaterialCutsCouple->GetMaterial()!=pNewMaterial )
                                                   >> 911       {
                                                   >> 912         // for parametrized volume
                                                   >> 913         //
                                                   >> 914         pNewMaterialCutsCouple = G4ProductionCutsTable::GetProductionCutsTable()
                                                   >> 915           ->GetMaterialCutsCouple(pNewMaterial,
                                                   >> 916                                   pNewMaterialCutsCouple->GetProductionCuts());
                                                   >> 917       }
                                                   >> 918   }
                                                   >> 919   fParticleChange.SetMaterialCutsCoupleInTouchable( pNewMaterialCutsCouple );
                                                   >> 920 
                                                   >> 921   // Must always set the touchable in ParticleChange, whether relocated or not
                                                   >> 922   //
                                                   >> 923   fParticleChange.SetTouchableHandle(retCurrentTouchable) ;
608                                                   924 
609   return &fParticleChange ;                       925   return &fParticleChange ;
610 }                                                 926 }
611                                                   927 
612 //////////////////////////////////////////////    928 /////////////////////////////////////////////////////////////////////////////
613 // New method takes over the responsibility to    929 // New method takes over the responsibility to reset the state of 
614 // G4CoupledTransportation object:                930 // G4CoupledTransportation object:
615 //      - at the start of a new track,  and       931 //      - at the start of a new track,  and
616 //      - on the resumption of a suspended tra    932 //      - on the resumption of a suspended track. 
617 //                                                933 //
618 void                                              934 void 
619 G4CoupledTransportation::StartTracking(G4Track    935 G4CoupledTransportation::StartTracking(G4Track* aTrack)
620 {                                                 936 {
621   G4Transportation::StartTracking(aTrack);     << 937 
                                                   >> 938   G4TransportationManager* transportMgr =
                                                   >> 939     G4TransportationManager::GetTransportationManager();
                                                   >> 940 
                                                   >> 941   // G4VProcess::StartTracking(aTrack);
                                                   >> 942   fNewTrack= true;
622                                                   943   
                                                   >> 944   //  The 'initialising' actions
                                                   >> 945   //     once taken in AlongStepGPIL -- if ( track.GetCurrentStepNumber()==1 )
                                                   >> 946 
                                                   >> 947   // fStartedNewTrack= true; 
                                                   >> 948 
                                                   >> 949   fMassNavigator = transportMgr->GetNavigatorForTracking() ; 
                                                   >> 950   fNavigatorId= transportMgr->ActivateNavigator( fMassNavigator );
                                                   >> 951 
623   G4ThreeVector position = aTrack->GetPosition    952   G4ThreeVector position = aTrack->GetPosition(); 
624   G4ThreeVector direction = aTrack->GetMomentu    953   G4ThreeVector direction = aTrack->GetMomentumDirection();
625                                                   954 
626   fPathFinder->PrepareNewTrack( position, dire    955   fPathFinder->PrepareNewTrack( position, direction); 
627   // This implies a call to fPathFinder->Locat    956   // This implies a call to fPathFinder->Locate( position, direction ); 
628                                                   957 
                                                   >> 958   // Whether field exists should be determined at run level -- TODO
                                                   >> 959   fAnyFieldExists= DoesAnyFieldExist(); 
                                                   >> 960 
629   // reset safety value and center                961   // reset safety value and center
630   //                                              962   //
631   fPreviousMassSafety  = 0.0 ;                    963   fPreviousMassSafety  = 0.0 ; 
632   fPreviousFullSafety  = 0.0 ;                    964   fPreviousFullSafety  = 0.0 ; 
633   fPreviousSftOrigin = G4ThreeVector(0.,0.,0.)    965   fPreviousSftOrigin = G4ThreeVector(0.,0.,0.) ;
                                                   >> 966   
                                                   >> 967   // reset looping counter -- for motion in field  
                                                   >> 968   fNoLooperTrials= 0; 
                                                   >> 969 
                                                   >> 970   // Must clear this state .. else it depends on last track's value
                                                   >> 971   //  --> a better solution would set this from state of suspended track TODO ? 
                                                   >> 972   // Was if( aTrack->GetCurrentStepNumber()==1 ) { .. }
                                                   >> 973 
                                                   >> 974   // ChordFinder reset internal state
                                                   >> 975   //
                                                   >> 976   if( fFieldPropagator && fAnyFieldExists )
                                                   >> 977   {
                                                   >> 978      fFieldPropagator->ClearPropagatorState();   
                                                   >> 979        // Resets safety values, in case of overlaps.  
                                                   >> 980 
                                                   >> 981      G4ChordFinder* chordF= fFieldPropagator->GetChordFinder();
                                                   >> 982      if( chordF )  { chordF->ResetStepEstimate(); }
                                                   >> 983   }
                                                   >> 984 
                                                   >> 985   // Clear the chord finders of all fields (ie managers) derived objects
                                                   >> 986   //
                                                   >> 987   G4FieldManagerStore* fieldMgrStore = G4FieldManagerStore::GetInstance();
                                                   >> 988   fieldMgrStore->ClearAllChordFindersState(); 
                                                   >> 989 
                                                   >> 990 #ifdef G4DEBUG_TRANSPORT
                                                   >> 991   if( verboseLevel > 1 )
                                                   >> 992   {
                                                   >> 993     G4cout << " Returning touchable handle " << fCurrentTouchableHandle
                                                   >> 994            << G4endl;
                                                   >> 995   }
                                                   >> 996 #endif
                                                   >> 997 
                                                   >> 998   // Update the current touchable handle  (from the track's)
                                                   >> 999   //
                                                   >> 1000   fCurrentTouchableHandle = aTrack->GetTouchableHandle();  
634 }                                                 1001 }
635                                                   1002 
636 //////////////////////////////////////////////    1003 /////////////////////////////////////////////////////////////////////////////
637                                                   1004 
638 void                                              1005 void 
639 G4CoupledTransportation::EndTracking()            1006 G4CoupledTransportation::EndTracking()
640 {                                                 1007 {
641   G4TransportationManager::GetTransportationMa    1008   G4TransportationManager::GetTransportationManager()->InactivateAll();
642   fPathFinder->EndTrack();                        1009   fPathFinder->EndTrack(); 
643     // Resets TransportationManager to use ord    1010     // Resets TransportationManager to use ordinary Navigator
644 }                                                 1011 }
645                                                   1012 
646 //////////////////////////////////////////////    1013 /////////////////////////////////////////////////////////////////////////////
647                                                   1014 
648 void                                              1015 void
649 G4CoupledTransportation::                         1016 G4CoupledTransportation::
650 ReportInexactEnergy(G4double startEnergy, G4do    1017 ReportInexactEnergy(G4double startEnergy, G4double endEnergy)
651 {                                                 1018 {
652   static G4ThreadLocal G4int no_warnings= 0, w    1019   static G4ThreadLocal G4int no_warnings= 0, warnModulo=1,
653                              moduloFactor= 10,    1020                              moduloFactor= 10, no_large_ediff= 0; 
654                                                   1021 
655   if( std::fabs(startEnergy- endEnergy) > perT    1022   if( std::fabs(startEnergy- endEnergy) > perThousand * endEnergy )
656   {                                               1023   {
657     no_large_ediff ++;                            1024     no_large_ediff ++;
658     if( (no_large_ediff% warnModulo) == 0 )       1025     if( (no_large_ediff% warnModulo) == 0 )
659     {                                             1026     {
660       no_warnings++;                              1027       no_warnings++;
661       std::ostringstream message;                 1028       std::ostringstream message;
662       message << "Energy change in Step is abo    1029       message << "Energy change in Step is above 1^-3 relative value. "
663               << G4endl                           1030               << G4endl
664               << "   Relative change in 'track    1031               << "   Relative change in 'tracking' step = " 
665               << std::setw(15) << (endEnergy-s    1032               << std::setw(15) << (endEnergy-startEnergy)/startEnergy
666               << G4endl                           1033               << G4endl
667               << "   Starting E= " << std::set    1034               << "   Starting E= " << std::setw(12) << startEnergy / MeV
668               << " MeV " << G4endl                1035               << " MeV " << G4endl
669               << "   Ending   E= " << std::set    1036               << "   Ending   E= " << std::setw(12) << endEnergy   / MeV
670               << " MeV " << G4endl                1037               << " MeV " << G4endl
671               << "Energy has been corrected --    1038               << "Energy has been corrected -- however, review"
672               << " field propagation parameter    1039               << " field propagation parameters for accuracy." << G4endl;
673       if ( (verboseLevel > 2 ) || (no_warnings    1040       if ( (verboseLevel > 2 ) || (no_warnings<4)
674         || (no_large_ediff == warnModulo * mod    1041         || (no_large_ediff == warnModulo * moduloFactor) )
675       {                                           1042       {
676         message << "These include EpsilonStepM    1043         message << "These include EpsilonStepMax(/Min) in G4FieldManager,"
677                 << G4endl                         1044                 << G4endl
678                 << "which determine fractional    1045                 << "which determine fractional error per step for integrated quantities."
679                 << G4endl                         1046                 << G4endl
680                 << "Note also the influence of    1047                 << "Note also the influence of the permitted number of integration steps."
681                 << G4endl;                        1048                 << G4endl;
682       }                                           1049       }
683       message << "Bad 'endpoint'. Energy chang    1050       message << "Bad 'endpoint'. Energy change detected and corrected."
684               << G4endl                           1051               << G4endl
685               << "Has occurred already " << no    1052               << "Has occurred already " << no_large_ediff << " times.";
686       G4Exception("G4CoupledTransportation::Al    1053       G4Exception("G4CoupledTransportation::AlongStepGetPIL()", 
687                   "EnergyChange", JustWarning,    1054                   "EnergyChange", JustWarning, message);
688       if( no_large_ediff == warnModulo * modul    1055       if( no_large_ediff == warnModulo * moduloFactor )
689       {                                           1056       {
690         warnModulo *= moduloFactor;               1057         warnModulo *= moduloFactor;
691       }                                           1058       }
692     }                                             1059     }
693   }                                               1060   }
                                                   >> 1061 }
                                                   >> 1062 
                                                   >> 1063 /////////////////////////////////////////////////////////////////////////////
                                                   >> 1064 
                                                   >> 1065 G4bool G4CoupledTransportation::EnableMagneticMoment(G4bool useMoment)
                                                   >> 1066 {
                                                   >> 1067   G4bool lastValue= fUseMagneticMoment;
                                                   >> 1068   fUseMagneticMoment= useMoment;
                                                   >> 1069   G4Transportation::fUseMagneticMoment= useMoment;
                                                   >> 1070   return lastValue;
                                                   >> 1071 }
                                                   >> 1072 
                                                   >> 1073 /////////////////////////////////////////////////////////////////////////////
                                                   >> 1074 //
                                                   >> 1075 
                                                   >> 1076 G4bool G4CoupledTransportation::EnableGravity(G4bool useGravity)
                                                   >> 1077 {
                                                   >> 1078   G4bool lastValue= fUseGravity;
                                                   >> 1079   fUseGravity= useGravity;
                                                   >> 1080   G4Transportation::fUseGravity= useGravity;
                                                   >> 1081   return lastValue;
                                                   >> 1082 }
                                                   >> 1083 
                                                   >> 1084 /////////////////////////////////////////////////////////////////////////////
                                                   >> 1085 //
                                                   >> 1086 //  Supress (or not) warnings about 'looping' particles
                                                   >> 1087 
                                                   >> 1088 void G4CoupledTransportation::SetSilenceLooperWarnings( G4bool val)
                                                   >> 1089 {
                                                   >> 1090   fSilenceLooperWarnings= val;  // Flag to *Supress* all 'looper' warnings  
                                                   >> 1091   // G4CoupledTransportation::fSilenceLooperWarnings= val;
                                                   >> 1092 }
                                                   >> 1093 
                                                   >> 1094 /////////////////////////////////////////////////////////////////////////////
                                                   >> 1095 //
                                                   >> 1096 G4bool G4CoupledTransportation::GetSilenceLooperWarnings()
                                                   >> 1097 {
                                                   >> 1098   return fSilenceLooperWarnings; 
                                                   >> 1099 }
                                                   >> 1100 
                                                   >> 1101 /////////////////////////////////////////////////////////////////////////////
                                                   >> 1102 //
                                                   >> 1103 void G4CoupledTransportation::SetHighLooperThresholds()
                                                   >> 1104 {
                                                   >> 1105   // Setting old 'high' values for thresholds - old values, potentially appropriate
                                                   >> 1106   //   for the energy frontier HEP experiments
                                                   >> 1107   SetThresholdWarningEnergy(    100.0 * CLHEP::MeV );  //  Warn above this energy
                                                   >> 1108   SetThresholdImportantEnergy(  250.0 * CLHEP::MeV );  //  Give a few trial above this E); 
                                                   >> 1109 
                                                   >> 1110   G4int maxTrials = 10;
                                                   >> 1111   SetThresholdTrials( maxTrials );
                                                   >> 1112 
                                                   >> 1113   if( verboseLevel )  ReportLooperThresholds();
                                                   >> 1114 }
                                                   >> 1115 
                                                   >> 1116 /////////////////////////////////////////////////////////////////////////////
                                                   >> 1117 void G4CoupledTransportation::SetLowLooperThresholds() // Values for low-E applications
                                                   >> 1118 {
                                                   >> 1119   // These values were the default in Geant4 10.5 - beta
                                                   >> 1120   SetThresholdWarningEnergy(     1.0 * CLHEP::keV ); // Warn above this En
                                                   >> 1121   SetThresholdImportantEnergy(   1.0 * CLHEP::MeV ); // Extra trials above it
                                                   >> 1122 
                                                   >> 1123   G4int maxTrials = 30; //  A new value - was 10
                                                   >> 1124   SetThresholdTrials( maxTrials );
                                                   >> 1125   
                                                   >> 1126   if( verboseLevel )  ReportLooperThresholds();  
                                                   >> 1127 }
                                                   >> 1128 
                                                   >> 1129 /////////////////////////////////////////////////////////////////////////////
                                                   >> 1130 //
                                                   >> 1131 void
                                                   >> 1132 G4CoupledTransportation::ReportMissingLogger( const char* methodName )
                                                   >> 1133 {
                                                   >> 1134    const char* message= "Logger object missing from G4CoupledTransportation";
                                                   >> 1135    G4String classAndMethod= G4String("G4CoupledTransportation") + G4String( methodName );
                                                   >> 1136    G4Exception(classAndMethod, "Missing Logger", JustWarning, message);
                                                   >> 1137 
                                                   >> 1138    if( verboseLevel )  ReportLooperThresholds();  
                                                   >> 1139 }
                                                   >> 1140 
                                                   >> 1141 /////////////////////////////////////////////////////////////////////////////
                                                   >> 1142 //
                                                   >> 1143 void
                                                   >> 1144 G4CoupledTransportation::ReportLooperThresholds()
                                                   >> 1145 {
                                                   >> 1146    PushThresholdsToLogger();  // To be absolutely certain they are in sync
                                                   >> 1147    fpLogger->ReportLooperThresholds("G4CoupledTransportation");
694 }                                                 1148 }
695                                                   1149