| Geant4 Cross Reference |
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 include file implementation 29 // GEANT 4 include file implementation
30 // 30 //
31 // ------------------------------------------- 31 // ------------------------------------------------------------
32 // 32 //
33 // This class is a process responsible for the 33 // This class is a process responsible for the transportation of
34 // a particle, ie the geometrical propagation 34 // a particle, ie the geometrical propagation that encounters the
35 // geometrical sub-volumes of the detectors. 35 // geometrical sub-volumes of the detectors.
36 // 36 //
37 // It is also tasked with the key role of prop 37 // It is also tasked with the key role of proposing the "isotropic safety",
38 // which will be used to update the post-ste 38 // which will be used to update the post-step point's safety.
39 // 39 //
40 // =========================================== 40 // =======================================================================
41 // Created: 19 March 1997, J. Apostolakis 41 // Created: 19 March 1997, J. Apostolakis
42 // =========================================== 42 // =======================================================================
43 43
44 #include "G4Transportation.hh" 44 #include "G4Transportation.hh"
45 #include "G4TransportationProcessType.hh" 45 #include "G4TransportationProcessType.hh"
46 #include "G4TransportationLogger.hh" 46 #include "G4TransportationLogger.hh"
47 47
48 #include "G4PhysicalConstants.hh" 48 #include "G4PhysicalConstants.hh"
49 #include "G4SystemOfUnits.hh" 49 #include "G4SystemOfUnits.hh"
50 #include "G4ProductionCutsTable.hh" 50 #include "G4ProductionCutsTable.hh"
51 #include "G4ParticleTable.hh" 51 #include "G4ParticleTable.hh"
52 52
53 #include "G4ChargeState.hh" 53 #include "G4ChargeState.hh"
54 #include "G4EquationOfMotion.hh" 54 #include "G4EquationOfMotion.hh"
55 55
56 #include "G4FieldManagerStore.hh" 56 #include "G4FieldManagerStore.hh"
57 57
58 #include "G4Navigator.hh" 58 #include "G4Navigator.hh"
59 #include "G4PropagatorInField.hh" 59 #include "G4PropagatorInField.hh"
60 #include "G4TransportationManager.hh" 60 #include "G4TransportationManager.hh"
61 61
62 #include "G4TransportationParameters.hh" 62 #include "G4TransportationParameters.hh"
63 63
64 class G4VSensitiveDetector; 64 class G4VSensitiveDetector;
65 65
66 G4bool G4Transportation::fUseMagneticMoment=fa 66 G4bool G4Transportation::fUseMagneticMoment=false;
67 G4bool G4Transportation::fUseGravity= false; 67 G4bool G4Transportation::fUseGravity= false;
68 G4bool G4Transportation::fSilenceLooperWarning 68 G4bool G4Transportation::fSilenceLooperWarnings= false;
69 69
70 ////////////////////////////////////////////// 70 //////////////////////////////////////////////////////////////////////////
71 // 71 //
72 // Constructor 72 // Constructor
73 73
74 G4Transportation::G4Transportation( G4int verb 74 G4Transportation::G4Transportation( G4int verbosity, const G4String& aName )
75 : G4VProcess( aName, fTransportation ), 75 : G4VProcess( aName, fTransportation ),
76 fFieldExertedForce( false ), 76 fFieldExertedForce( false ),
77 fPreviousSftOrigin( 0.,0.,0. ), 77 fPreviousSftOrigin( 0.,0.,0. ),
78 fPreviousSafety( 0.0 ), 78 fPreviousSafety( 0.0 ),
79 fEndPointDistance( -1.0 ), 79 fEndPointDistance( -1.0 ),
80 fShortStepOptimisation( false ) // Old def 80 fShortStepOptimisation( false ) // Old default: true (=fast short steps)
81 { 81 {
82 SetProcessSubType(static_cast<G4int>(TRANSPO 82 SetProcessSubType(static_cast<G4int>(TRANSPORTATION));
83 pParticleChange= &fParticleChange; // Requ 83 pParticleChange= &fParticleChange; // Required to conform to G4VProcess
84 SetVerboseLevel(verbosity); 84 SetVerboseLevel(verbosity);
85 85
86 G4TransportationManager* transportMgr ; 86 G4TransportationManager* transportMgr ;
87 87
88 transportMgr = G4TransportationManager::GetT 88 transportMgr = G4TransportationManager::GetTransportationManager() ;
89 89
90 fLinearNavigator = transportMgr->GetNavigato 90 fLinearNavigator = transportMgr->GetNavigatorForTracking() ;
91 91
92 fFieldPropagator = transportMgr->GetPropagat 92 fFieldPropagator = transportMgr->GetPropagatorInField() ;
93 93
94 fpSafetyHelper = transportMgr->GetSafetyHe 94 fpSafetyHelper = transportMgr->GetSafetyHelper(); // New
95 95
96 fpLogger = new G4TransportationLogger("G4Tra 96 fpLogger = new G4TransportationLogger("G4Transportation", verbosity);
97 97
98 if( G4TransportationParameters::Exists() ) 98 if( G4TransportationParameters::Exists() )
99 { 99 {
100 auto trParams= G4TransportationParameters: 100 auto trParams= G4TransportationParameters::Instance();
101 101
102 SetThresholdWarningEnergy( trParams->GetW 102 SetThresholdWarningEnergy( trParams->GetWarningEnergy() );
103 SetThresholdImportantEnergy( trParams->Get 103 SetThresholdImportantEnergy( trParams->GetImportantEnergy() );
104 SetThresholdTrials( trParams->GetNumberOfT 104 SetThresholdTrials( trParams->GetNumberOfTrials() );
105 G4Transportation::fSilenceLooperWarnings= 105 G4Transportation::fSilenceLooperWarnings= trParams->GetSilenceAllLooperWarnings();
106 } 106 }
107 else { 107 else {
108 SetHighLooperThresholds(); 108 SetHighLooperThresholds();
109 // Use the old defaults: Warning = 100 Me 109 // Use the old defaults: Warning = 100 MeV, Important = 250 MeV, No Trials = 10;
110 } 110 }
111 111
112 PushThresholdsToLogger(); 112 PushThresholdsToLogger();
113 // Should be done by Set methods in SetHighL 113 // Should be done by Set methods in SetHighLooperThresholds -- making sure
114 114
115 static G4ThreadLocal G4TouchableHandle* pNul 115 static G4ThreadLocal G4TouchableHandle* pNullTouchableHandle = 0;
116 if ( !pNullTouchableHandle) 116 if ( !pNullTouchableHandle)
117 { 117 {
118 pNullTouchableHandle = new G4TouchableHand 118 pNullTouchableHandle = new G4TouchableHandle;
119 } 119 }
120 fCurrentTouchableHandle = *pNullTouchableHan 120 fCurrentTouchableHandle = *pNullTouchableHandle;
121 // Points to (G4VTouchable*) 0 121 // Points to (G4VTouchable*) 0
122 122
123 123
124 #ifdef G4VERBOSE 124 #ifdef G4VERBOSE
125 if( verboseLevel > 0) 125 if( verboseLevel > 0)
126 { 126 {
127 G4cout << " G4Transportation constructor> 127 G4cout << " G4Transportation constructor> set fShortStepOptimisation to ";
128 if ( fShortStepOptimisation ) { G4cout < 128 if ( fShortStepOptimisation ) { G4cout << "true" << G4endl; }
129 else { G4cout < 129 else { G4cout << "false" << G4endl; }
130 } 130 }
131 #endif 131 #endif
132 } 132 }
133 133
134 ////////////////////////////////////////////// 134 //////////////////////////////////////////////////////////////////////////
135 135
136 G4Transportation::~G4Transportation() 136 G4Transportation::~G4Transportation()
137 { 137 {
138 if( fSumEnergyKilled > 0.0 ) 138 if( fSumEnergyKilled > 0.0 )
139 { 139 {
140 PrintStatistics( G4cout ); 140 PrintStatistics( G4cout );
141 } 141 }
142 delete fpLogger; 142 delete fpLogger;
143 } 143 }
144 144
145 ////////////////////////////////////////////// 145 //////////////////////////////////////////////////////////////////////////
146 146
147 void 147 void
148 G4Transportation::PrintStatistics( std::ostrea 148 G4Transportation::PrintStatistics( std::ostream& outStr) const
149 { 149 {
150 outStr << " G4Transportation: Statistics fo 150 outStr << " G4Transportation: Statistics for looping particles " << G4endl;
151 if( fSumEnergyKilled > 0.0 || fNumLoopersKi 151 if( fSumEnergyKilled > 0.0 || fNumLoopersKilled > 0 )
152 { 152 {
153 outStr << " Sum of energy of looping t 153 outStr << " Sum of energy of looping tracks killed: "
154 << fSumEnergyKilled / CLHEP::MeV 154 << fSumEnergyKilled / CLHEP::MeV << " MeV "
155 << " from " << fNumLoopersKilled 155 << " from " << fNumLoopersKilled << " tracks " << G4endl
156 << " Sum of energy of non-elect 156 << " Sum of energy of non-electrons : "
157 << fSumEnergyKilled_NonElectron / 157 << fSumEnergyKilled_NonElectron / CLHEP::MeV << " MeV "
158 << " from " << fNumLoopersKilled 158 << " from " << fNumLoopersKilled_NonElectron << " tracks "
159 << G4endl; 159 << G4endl;
160 outStr << " Max energy of *any type* 160 outStr << " Max energy of *any type* looper killed: " << fMaxEnergyKilled
161 << " its PDG was " << fMaxEner 161 << " its PDG was " << fMaxEnergyKilledPDG << G4endl;
162 if( fMaxEnergyKilled_NonElectron > 0.0 ) 162 if( fMaxEnergyKilled_NonElectron > 0.0 )
163 { 163 {
164 outStr << " Max energy of non-elect 164 outStr << " Max energy of non-electron looper killed: "
165 << fMaxEnergyKilled_NonElectro 165 << fMaxEnergyKilled_NonElectron
166 << " its PDG was " << fMaxE 166 << " its PDG was " << fMaxEnergyKilled_NonElecPDG << G4endl;
167 } 167 }
168 if( fMaxEnergySaved > 0.0 ) 168 if( fMaxEnergySaved > 0.0 )
169 { 169 {
170 outStr << " Max energy of loopers ' 170 outStr << " Max energy of loopers 'saved': " << fMaxEnergySaved << G4endl;
171 outStr << " Sum of energy of looper 171 outStr << " Sum of energy of loopers 'saved': "
172 << fSumEnergySaved << G4endl; 172 << fSumEnergySaved << G4endl;
173 outStr << " Sum of energy of unstab 173 outStr << " Sum of energy of unstable loopers 'saved': "
174 << fSumEnergyUnstableSaved << 174 << fSumEnergyUnstableSaved << G4endl;
175 } 175 }
176 } 176 }
177 else 177 else
178 { 178 {
179 outStr << " No looping tracks found or k 179 outStr << " No looping tracks found or killed. " << G4endl;
180 } 180 }
181 } 181 }
182 182
183 ////////////////////////////////////////////// 183 //////////////////////////////////////////////////////////////////////////
184 // 184 //
185 // Responsibilities: 185 // Responsibilities:
186 // Find whether the geometry limits the Ste 186 // Find whether the geometry limits the Step, and to what length
187 // Calculate the new value of the safety an 187 // Calculate the new value of the safety and return it.
188 // Store the final time, position and momen 188 // Store the final time, position and momentum.
189 189
190 G4double G4Transportation::AlongStepGetPhysica 190 G4double G4Transportation::AlongStepGetPhysicalInteractionLength(
191 const G4Track& track, 191 const G4Track& track,
192 G4double, // previousStepSize 192 G4double, // previousStepSize
193 G4double currentMinimumStep, G4double& curre 193 G4double currentMinimumStep, G4double& currentSafety,
194 G4GPILSelection* selection) 194 G4GPILSelection* selection)
195 { 195 {
196 // Initial actions moved to StartTrack() 196 // Initial actions moved to StartTrack()
197 // -------------------------------------- 197 // --------------------------------------
198 // Note: in case another process changes tou 198 // Note: in case another process changes touchable handle
199 // it will be necessary to add here (for 199 // it will be necessary to add here (for all steps)
200 // fCurrentTouchableHandle = aTrack->GetTouc 200 // fCurrentTouchableHandle = aTrack->GetTouchableHandle();
201 201
202 // GPILSelection is set to defaule value of 202 // GPILSelection is set to defaule value of CandidateForSelection
203 // It is a return value 203 // It is a return value
204 // 204 //
205 *selection = CandidateForSelection; 205 *selection = CandidateForSelection;
206 206
207 // Get initial Energy/Momentum of the track 207 // Get initial Energy/Momentum of the track
208 // 208 //
209 const G4ThreeVector startPosition = track 209 const G4ThreeVector startPosition = track.GetPosition();
210 const G4ThreeVector startMomentumDir = track 210 const G4ThreeVector startMomentumDir = track.GetMomentumDirection();
211 211
212 // The Step Point safety can be limited by o 212 // The Step Point safety can be limited by other geometries and/or the
213 // assumptions of any process - it's not alw 213 // assumptions of any process - it's not always the geometrical safety.
214 // We calculate the starting point's isotrop 214 // We calculate the starting point's isotropic safety here.
215 { 215 {
216 const G4double MagSqShift = (startPosition 216 const G4double MagSqShift = (startPosition - fPreviousSftOrigin).mag2();
217 217
218 if(MagSqShift >= sqr(fPreviousSafety)) 218 if(MagSqShift >= sqr(fPreviousSafety))
219 currentSafety = 0.0; 219 currentSafety = 0.0;
220 else 220 else
221 currentSafety = fPreviousSafety - std::s 221 currentSafety = fPreviousSafety - std::sqrt(MagSqShift);
222 } 222 }
223 223
224 // Is the particle charged or has it a magne 224 // Is the particle charged or has it a magnetic moment?
225 // 225 //
226 const G4DynamicParticle* pParticle = track.G 226 const G4DynamicParticle* pParticle = track.GetDynamicParticle();
227 227
228 const G4double particleMass = pParticle->G 228 const G4double particleMass = pParticle->GetMass();
229 const G4double particleCharge = pParticle->G 229 const G4double particleCharge = pParticle->GetCharge();
230 const G4double kineticEnergy = pParticle->Ge 230 const G4double kineticEnergy = pParticle->GetKineticEnergy();
231 231
232 const G4double magneticMoment = pParticle 232 const G4double magneticMoment = pParticle->GetMagneticMoment();
233 const G4ThreeVector particleSpin = pParticle 233 const G4ThreeVector particleSpin = pParticle->GetPolarization();
234 234
235 // There is no need to locate the current vo 235 // There is no need to locate the current volume. It is Done elsewhere:
236 // On track construction 236 // On track construction
237 // By the tracking, after all AlongStepDoI 237 // By the tracking, after all AlongStepDoIts, in "Relocation"
238 238
239 // Check if the particle has a force, EM or 239 // Check if the particle has a force, EM or gravitational, exerted on it
240 // 240 //
241 241
242 G4bool eligibleEM = 242 G4bool eligibleEM =
243 (particleCharge != 0.0) || ((magneticMomen 243 (particleCharge != 0.0) || ((magneticMoment != 0.0) && fUseMagneticMoment);
244 G4bool eligibleGrav = (particleMass != 0.0) 244 G4bool eligibleGrav = (particleMass != 0.0) && fUseGravity;
245 245
246 fFieldExertedForce = false; 246 fFieldExertedForce = false;
247 247
248 if(eligibleEM || eligibleGrav) 248 if(eligibleEM || eligibleGrav)
249 { 249 {
250 if(G4FieldManager* fieldMgr = 250 if(G4FieldManager* fieldMgr =
251 fFieldPropagator->FindAndSetFieldMana 251 fFieldPropagator->FindAndSetFieldManager(track.GetVolume()))
252 { 252 {
253 // User can configure the field Manager 253 // User can configure the field Manager for this track
254 fieldMgr->ConfigureForTrack(&track); 254 fieldMgr->ConfigureForTrack(&track);
255 // Called here to allow a transition fro 255 // Called here to allow a transition from no-field pointer
256 // to finite field (non-zero pointer). 256 // to finite field (non-zero pointer).
257 257
258 // If the field manager has no field ptr 258 // If the field manager has no field ptr, the field is zero
259 // by definition ( = there is no field 259 // by definition ( = there is no field ! )
260 if(const G4Field* ptrField = fieldMgr->G 260 if(const G4Field* ptrField = fieldMgr->GetDetectorField())
261 fFieldExertedForce = 261 fFieldExertedForce =
262 eligibleEM || (eligibleGrav && ptrFi 262 eligibleEM || (eligibleGrav && ptrField->IsGravityActive());
263 } 263 }
264 } 264 }
265 265
266 G4double geometryStepLength = currentMinimum 266 G4double geometryStepLength = currentMinimumStep;
267 267
268 if(currentMinimumStep == 0.0) 268 if(currentMinimumStep == 0.0)
269 { 269 {
270 fEndPointDistance = 0.0; 270 fEndPointDistance = 0.0;
271 fGeometryLimitedStep = false; // Old cod 271 fGeometryLimitedStep = false; // Old code: = (currentSafety == 0.0);
272 // Changed to avoid problems when setting 272 // Changed to avoid problems when setting the step status (now done in AlongStepDoIt)
273 fMomentumChanged = false; 273 fMomentumChanged = false;
274 fParticleIsLooping = false; 274 fParticleIsLooping = false;
275 fEndGlobalTimeComputed = false; 275 fEndGlobalTimeComputed = false;
276 fTransportEndPosition = startPosition 276 fTransportEndPosition = startPosition;
277 fTransportEndMomentumDir = startMomentum 277 fTransportEndMomentumDir = startMomentumDir;
278 fTransportEndKineticEnergy = kineticEnergy 278 fTransportEndKineticEnergy = kineticEnergy;
279 fTransportEndSpin = particleSpin; 279 fTransportEndSpin = particleSpin;
280 } 280 }
281 else if(!fFieldExertedForce) 281 else if(!fFieldExertedForce)
282 { 282 {
283 fGeometryLimitedStep = false; 283 fGeometryLimitedStep = false;
284 if(geometryStepLength > currentSafety || ! 284 if(geometryStepLength > currentSafety || !fShortStepOptimisation)
285 { 285 {
286 const G4double linearStepLength = fLinea 286 const G4double linearStepLength = fLinearNavigator->ComputeStep(
287 startPosition, startMomentumDir, curre 287 startPosition, startMomentumDir, currentMinimumStep, currentSafety);
288 288
289 if(linearStepLength <= currentMinimumSte 289 if(linearStepLength <= currentMinimumStep)
290 { 290 {
291 geometryStepLength = linearStepLength; 291 geometryStepLength = linearStepLength;
292 fGeometryLimitedStep = true; 292 fGeometryLimitedStep = true;
293 } 293 }
294 // Remember last safety origin & value. 294 // Remember last safety origin & value.
295 // 295 //
296 fPreviousSftOrigin = startPosition; 296 fPreviousSftOrigin = startPosition;
297 fPreviousSafety = currentSafety; 297 fPreviousSafety = currentSafety;
298 fpSafetyHelper->SetCurrentSafety(current 298 fpSafetyHelper->SetCurrentSafety(currentSafety, startPosition);
299 } 299 }
300 300
301 fEndPointDistance = geometryStepLength; 301 fEndPointDistance = geometryStepLength;
302 302
303 fMomentumChanged = false; 303 fMomentumChanged = false;
304 fParticleIsLooping = false; 304 fParticleIsLooping = false;
305 fEndGlobalTimeComputed = false; 305 fEndGlobalTimeComputed = false;
306 fTransportEndPosition = 306 fTransportEndPosition =
307 startPosition + geometryStepLength * sta 307 startPosition + geometryStepLength * startMomentumDir;
308 fTransportEndMomentumDir = startMomentum 308 fTransportEndMomentumDir = startMomentumDir;
309 fTransportEndKineticEnergy = kineticEnergy 309 fTransportEndKineticEnergy = kineticEnergy;
310 fTransportEndSpin = particleSpin; 310 fTransportEndSpin = particleSpin;
311 } 311 }
312 else // A field exerts force 312 else // A field exerts force
313 { 313 {
314 const auto pParticleDef = pParticle->Ge 314 const auto pParticleDef = pParticle->GetDefinition();
315 const auto particlePDGSpin = pParticleDef- 315 const auto particlePDGSpin = pParticleDef->GetPDGSpin();
316 const auto particlePDGMagM = pParticleDef- 316 const auto particlePDGMagM = pParticleDef->GetPDGMagneticMoment();
317 317
318 auto equationOfMotion = fFieldPropagator-> 318 auto equationOfMotion = fFieldPropagator->GetCurrentEquationOfMotion();
319 319
320 // The charge can change (dynamic), theref 320 // The charge can change (dynamic), therefore the use of G4ChargeState
321 // 321 //
322 equationOfMotion->SetChargeMomentumMass( 322 equationOfMotion->SetChargeMomentumMass(
323 G4ChargeState(particleCharge, magneticMo 323 G4ChargeState(particleCharge, magneticMoment, particlePDGSpin),
324 pParticle->GetTotalMomentum(), particleM 324 pParticle->GetTotalMomentum(), particleMass);
325 325
326 G4FieldTrack aFieldTrack(startPosition, 326 G4FieldTrack aFieldTrack(startPosition,
327 track.GetGlobalTi 327 track.GetGlobalTime(), // Lab.
328 startMomentumDir, 328 startMomentumDir, kineticEnergy, particleMass,
329 particleCharge, p 329 particleCharge, particleSpin, particlePDGMagM,
330 0.0, // Length a 330 0.0, // Length along track
331 particlePDGSpin); 331 particlePDGSpin);
332 332
333 // Do the Transport in the field (non rect 333 // Do the Transport in the field (non recti-linear)
334 // 334 //
335 const G4double lengthAlongCurve = fFieldPr 335 const G4double lengthAlongCurve = fFieldPropagator->ComputeStep(
336 aFieldTrack, currentMinimumStep, current 336 aFieldTrack, currentMinimumStep, currentSafety, track.GetVolume(),
337 kineticEnergy < fThreshold_Important_Ene 337 kineticEnergy < fThreshold_Important_Energy);
338 338
339 if(lengthAlongCurve < geometryStepLength) 339 if(lengthAlongCurve < geometryStepLength)
340 geometryStepLength = lengthAlongCurve; 340 geometryStepLength = lengthAlongCurve;
341 341
342 // Remember last safety origin & value. 342 // Remember last safety origin & value.
343 // 343 //
344 fPreviousSftOrigin = startPosition; 344 fPreviousSftOrigin = startPosition;
345 fPreviousSafety = currentSafety; 345 fPreviousSafety = currentSafety;
346 fpSafetyHelper->SetCurrentSafety(currentSa 346 fpSafetyHelper->SetCurrentSafety(currentSafety, startPosition);
347 347
348 fGeometryLimitedStep = fFieldPropagator->I 348 fGeometryLimitedStep = fFieldPropagator->IsLastStepInVolume();
349 // 349 //
350 // It is possible that step was reduced in 350 // It is possible that step was reduced in PropagatorInField due to
351 // previous zero steps. To cope with case 351 // previous zero steps. To cope with case that reduced step is taken
352 // in full, we must rely on PiF to obtain 352 // in full, we must rely on PiF to obtain this value
353 353
354 G4bool changesEnergy = 354 G4bool changesEnergy =
355 fFieldPropagator->GetCurrentFieldManager 355 fFieldPropagator->GetCurrentFieldManager()->DoesFieldChangeEnergy();
356 356
357 fMomentumChanged = true; 357 fMomentumChanged = true;
358 fParticleIsLooping = fFieldPropagator->IsP 358 fParticleIsLooping = fFieldPropagator->IsParticleLooping();
359 359
360 fEndGlobalTimeComputed = changesEnergy; 360 fEndGlobalTimeComputed = changesEnergy;
361 fTransportEndPosition = aFieldTrack.Get 361 fTransportEndPosition = aFieldTrack.GetPosition();
362 fTransportEndMomentumDir = aFieldTrack.Get 362 fTransportEndMomentumDir = aFieldTrack.GetMomentumDir();
363 363
364 // G4cout << " G4Transport: End of step pM 364 // G4cout << " G4Transport: End of step pMag = " << aFieldTrack.GetMomentum().mag() << G4endl;
365 365
366 fEndPointDistance = (fTransportEndPositio 366 fEndPointDistance = (fTransportEndPosition - startPosition).mag();
367 367
368 // Ignore change in energy for fields that 368 // Ignore change in energy for fields that conserve energy
369 // This hides the integration error, but g 369 // This hides the integration error, but gives a better physical answer
370 fTransportEndKineticEnergy = 370 fTransportEndKineticEnergy =
371 changesEnergy ? aFieldTrack.GetKineticEn 371 changesEnergy ? aFieldTrack.GetKineticEnergy() : kineticEnergy;
372 fTransportEndSpin = aFieldTrack.GetSpin(); 372 fTransportEndSpin = aFieldTrack.GetSpin();
373 373
374 if(fEndGlobalTimeComputed) 374 if(fEndGlobalTimeComputed)
375 { 375 {
376 // If the field can change energy, then 376 // If the field can change energy, then the time must be integrated
377 // - so this should have been updated 377 // - so this should have been updated
378 // 378 //
379 fCandidateEndGlobalTime = aFieldTrack.Ge 379 fCandidateEndGlobalTime = aFieldTrack.GetLabTimeOfFlight();
380 380
381 // was ( fCandidateEndGlobalTime != trac 381 // was ( fCandidateEndGlobalTime != track.GetGlobalTime() );
382 // a cleaner way is to have FieldTrack k 382 // a cleaner way is to have FieldTrack knowing whether time is updated.
383 } 383 }
384 #if defined(G4VERBOSE) || defined(G4DEBUG_TRAN 384 #if defined(G4VERBOSE) || defined(G4DEBUG_TRANSPORT)
385 else 385 else
386 { 386 {
387 // The energy should be unchanged by fie 387 // The energy should be unchanged by field transport,
388 // - so the time changed will be calc 388 // - so the time changed will be calculated elsewhere
389 // 389 //
390 // Check that the integration preserved 390 // Check that the integration preserved the energy
391 // - and if not correct this! 391 // - and if not correct this!
392 G4double startEnergy = kineticEnergy; 392 G4double startEnergy = kineticEnergy;
393 G4double endEnergy = fTransportEndKine 393 G4double endEnergy = fTransportEndKineticEnergy;
394 394
395 static G4ThreadLocal G4int no_large_edif 395 static G4ThreadLocal G4int no_large_ediff;
396 if(verboseLevel > 1) 396 if(verboseLevel > 1)
397 { 397 {
398 if(std::fabs(startEnergy - endEnergy) 398 if(std::fabs(startEnergy - endEnergy) > perThousand * endEnergy)
399 { 399 {
400 static G4ThreadLocal G4int no_warnin 400 static G4ThreadLocal G4int no_warnings = 0, warnModulo = 1,
401 moduloFac 401 moduloFactor = 10;
402 no_large_ediff++; 402 no_large_ediff++;
403 if((no_large_ediff % warnModulo) == 403 if((no_large_ediff % warnModulo) == 0)
404 { 404 {
405 no_warnings++; 405 no_warnings++;
406 std::ostringstream message; 406 std::ostringstream message;
407 message << "Energy change in Step 407 message << "Energy change in Step is above 1^-3 relative value. "
408 << G4endl << " Relativ 408 << G4endl << " Relative change in 'tracking' step = "
409 << std::setw(15) << (endEn 409 << std::setw(15) << (endEnergy - startEnergy) / startEnergy
410 << G4endl << " Startin 410 << G4endl << " Starting E= " << std::setw(12)
411 << startEnergy / MeV << " 411 << startEnergy / MeV << " MeV " << G4endl
412 << " Ending E= " << 412 << " Ending E= " << std::setw(12) << endEnergy / MeV
413 << " MeV " << G4endl 413 << " MeV " << G4endl
414 << "Energy has been correc 414 << "Energy has been corrected -- however, review"
415 << " field propagation par 415 << " field propagation parameters for accuracy." << G4endl;
416 if((verboseLevel > 2) || (no_warni 416 if((verboseLevel > 2) || (no_warnings < 4) ||
417 (no_large_ediff == warnModulo * 417 (no_large_ediff == warnModulo * moduloFactor))
418 { 418 {
419 message << "These include Epsilo 419 message << "These include EpsilonStepMax(/Min) in G4FieldManager "
420 << G4endl 420 << G4endl
421 << "which determine frac 421 << "which determine fractional error per step for "
422 "integrated quantitie 422 "integrated quantities. "
423 << G4endl 423 << G4endl
424 << "Note also the influe 424 << "Note also the influence of the permitted number of "
425 "integration steps." 425 "integration steps."
426 << G4endl; 426 << G4endl;
427 } 427 }
428 message << "Bad 'endpoint'. Energy 428 message << "Bad 'endpoint'. Energy change detected and corrected."
429 << G4endl << "Has occurred 429 << G4endl << "Has occurred already " << no_large_ediff
430 << " times."; 430 << " times.";
431 G4Exception("G4Transportation::Alo 431 G4Exception("G4Transportation::AlongStepGetPIL()", "EnergyChange",
432 JustWarning, message); 432 JustWarning, message);
433 if(no_large_ediff == warnModulo * 433 if(no_large_ediff == warnModulo * moduloFactor)
434 { 434 {
435 warnModulo *= moduloFactor; 435 warnModulo *= moduloFactor;
436 } 436 }
437 } 437 }
438 } 438 }
439 } // end of if (verboseLevel) 439 } // end of if (verboseLevel)
440 } 440 }
441 #endif 441 #endif
442 } 442 }
443 443
444 // Update the safety starting from the end-p 444 // Update the safety starting from the end-point,
445 // if it will become negative at the end-poi 445 // if it will become negative at the end-point.
446 // 446 //
447 if(currentSafety < fEndPointDistance) 447 if(currentSafety < fEndPointDistance)
448 { 448 {
449 if(particleCharge != 0.0) 449 if(particleCharge != 0.0)
450 { 450 {
451 G4double endSafety = 451 G4double endSafety =
452 fLinearNavigator->ComputeSafety(fTrans 452 fLinearNavigator->ComputeSafety(fTransportEndPosition);
453 currentSafety = endSafety; 453 currentSafety = endSafety;
454 fPreviousSftOrigin = fTransportEndPositi 454 fPreviousSftOrigin = fTransportEndPosition;
455 fPreviousSafety = currentSafety; 455 fPreviousSafety = currentSafety;
456 fpSafetyHelper->SetCurrentSafety(current 456 fpSafetyHelper->SetCurrentSafety(currentSafety, fTransportEndPosition);
457 457
458 // Because the Stepping Manager assumes 458 // Because the Stepping Manager assumes it is from the start point,
459 // add the StepLength 459 // add the StepLength
460 // 460 //
461 currentSafety += fEndPointDistance; 461 currentSafety += fEndPointDistance;
462 462
463 #ifdef G4DEBUG_TRANSPORT 463 #ifdef G4DEBUG_TRANSPORT
464 G4cout.precision(12); 464 G4cout.precision(12);
465 G4cout << "***G4Transportation::AlongSte 465 G4cout << "***G4Transportation::AlongStepGPIL ** " << G4endl;
466 G4cout << " Called Navigator->ComputeSa 466 G4cout << " Called Navigator->ComputeSafety at " << fTransportEndPosition
467 << " and it returned safety= 467 << " and it returned safety= " << endSafety << G4endl;
468 G4cout << " Adding endpoint distance 468 G4cout << " Adding endpoint distance " << fEndPointDistance
469 << " to obtain pseudo-safety= 469 << " to obtain pseudo-safety= " << currentSafety << G4endl;
470 } 470 }
471 else 471 else
472 { 472 {
473 G4cout << "***G4Transportation::AlongSte 473 G4cout << "***G4Transportation::AlongStepGPIL ** " << G4endl;
474 G4cout << " Avoiding call to ComputeSaf 474 G4cout << " Avoiding call to ComputeSafety : " << G4endl;
475 G4cout << " charge = " << particl 475 G4cout << " charge = " << particleCharge << G4endl;
476 G4cout << " mag moment = " << magneti 476 G4cout << " mag moment = " << magneticMoment << G4endl;
477 #endif 477 #endif
478 } 478 }
479 } 479 }
480 480
481 fFirstStepInVolume = fNewTrack || fLastStepI 481 fFirstStepInVolume = fNewTrack || fLastStepInVolume;
482 fLastStepInVolume = false; 482 fLastStepInVolume = false;
483 fNewTrack = false; 483 fNewTrack = false;
484 484
485 fParticleChange.ProposeFirstStepInVolume(fFi 485 fParticleChange.ProposeFirstStepInVolume(fFirstStepInVolume);
486 fParticleChange.ProposeTrueStepLength(geomet 486 fParticleChange.ProposeTrueStepLength(geometryStepLength);
487 487
488 return geometryStepLength; 488 return geometryStepLength;
489 } 489 }
490 490
491 ////////////////////////////////////////////// 491 //////////////////////////////////////////////////////////////////////////
492 // 492 //
493 // Initialize ParticleChange (by setting al 493 // Initialize ParticleChange (by setting all its members equal
494 // to correspond 494 // to corresponding members in G4Track)
495 495
496 G4VParticleChange* G4Transportation::AlongStep 496 G4VParticleChange* G4Transportation::AlongStepDoIt( const G4Track& track,
497 497 const G4Step& stepData )
498 { 498 {
499 #if defined(G4VERBOSE) || defined(G4DEBUG_TRAN 499 #if defined(G4VERBOSE) || defined(G4DEBUG_TRANSPORT)
500 static G4ThreadLocal G4long noCallsASDI=0; 500 static G4ThreadLocal G4long noCallsASDI=0;
501 noCallsASDI++; 501 noCallsASDI++;
502 #else 502 #else
503 #define noCallsASDI 0 503 #define noCallsASDI 0
504 #endif 504 #endif
505 505
506 if(fGeometryLimitedStep) 506 if(fGeometryLimitedStep)
507 { 507 {
508 stepData.GetPostStepPoint()->SetStepStatus 508 stepData.GetPostStepPoint()->SetStepStatus(fGeomBoundary);
509 } 509 }
510 510
511 fParticleChange.Initialize(track) ; 511 fParticleChange.Initialize(track) ;
512 512
513 // Code for specific process 513 // Code for specific process
514 // 514 //
515 fParticleChange.ProposePosition(fTransportEn 515 fParticleChange.ProposePosition(fTransportEndPosition) ;
516 fParticleChange.ProposeMomentumDirection(fTr 516 fParticleChange.ProposeMomentumDirection(fTransportEndMomentumDir) ;
517 fParticleChange.ProposeEnergy(fTransportEndK 517 fParticleChange.ProposeEnergy(fTransportEndKineticEnergy) ;
518 fParticleChange.SetMomentumChanged(fMomentum 518 fParticleChange.SetMomentumChanged(fMomentumChanged) ;
519 519
520 fParticleChange.ProposePolarization(fTranspo 520 fParticleChange.ProposePolarization(fTransportEndSpin);
521 521
522 G4double deltaTime = 0.0 ; 522 G4double deltaTime = 0.0 ;
523 523
524 // Calculate Lab Time of Flight (ONLY if fi 524 // Calculate Lab Time of Flight (ONLY if field Equations used it!)
525 // G4double endTime = fCandidateEndGlobalT 525 // G4double endTime = fCandidateEndGlobalTime;
526 // G4double delta_time = endTime - startTime 526 // G4double delta_time = endTime - startTime;
527 527
528 G4double startTime = track.GetGlobalTime() ; 528 G4double startTime = track.GetGlobalTime() ;
529 529
530 if (!fEndGlobalTimeComputed) 530 if (!fEndGlobalTimeComputed)
531 { 531 {
532 // The time was not integrated .. make th 532 // The time was not integrated .. make the best estimate possible
533 // 533 //
534 G4double initialVelocity = stepData.GetPr 534 G4double initialVelocity = stepData.GetPreStepPoint()->GetVelocity();
535 G4double stepLength = track.GetStepL 535 G4double stepLength = track.GetStepLength();
536 536
537 deltaTime= 0.0; // in case initialVeloci 537 deltaTime= 0.0; // in case initialVelocity = 0
538 if ( initialVelocity > 0.0 ) { deltaTime 538 if ( initialVelocity > 0.0 ) { deltaTime = stepLength/initialVelocity; }
539 539
540 fCandidateEndGlobalTime = startTime + d 540 fCandidateEndGlobalTime = startTime + deltaTime ;
541 fParticleChange.ProposeLocalTime( track. 541 fParticleChange.ProposeLocalTime( track.GetLocalTime() + deltaTime) ;
542 } 542 }
543 else 543 else
544 { 544 {
545 deltaTime = fCandidateEndGlobalTime - sta 545 deltaTime = fCandidateEndGlobalTime - startTime ;
546 fParticleChange.ProposeGlobalTime( fCandi 546 fParticleChange.ProposeGlobalTime( fCandidateEndGlobalTime ) ;
547 } 547 }
548 548
549 549
550 // Now Correct by Lorentz factor to get delt 550 // Now Correct by Lorentz factor to get delta "proper" Time
551 551
552 G4double restMass = track.GetDynamicP 552 G4double restMass = track.GetDynamicParticle()->GetMass() ;
553 G4double deltaProperTime = deltaTime*( restM 553 G4double deltaProperTime = deltaTime*( restMass/track.GetTotalEnergy() ) ;
554 554
555 fParticleChange.ProposeProperTime(track.GetP 555 fParticleChange.ProposeProperTime(track.GetProperTime() + deltaProperTime) ;
556 //fParticleChange.ProposeTrueStepLength( tra 556 //fParticleChange.ProposeTrueStepLength( track.GetStepLength() ) ;
557 557
558 // If the particle is caught looping or is s 558 // If the particle is caught looping or is stuck (in very difficult
559 // boundaries) in a magnetic field (doing ma 559 // boundaries) in a magnetic field (doing many steps) THEN can kill it ...
560 // 560 //
561 if ( fParticleIsLooping ) 561 if ( fParticleIsLooping )
562 { 562 {
563 G4double endEnergy= fTransportEndKinetic 563 G4double endEnergy= fTransportEndKineticEnergy;
564 fNoLooperTrials ++; 564 fNoLooperTrials ++;
565 auto particleType= track.GetDynamicParti 565 auto particleType= track.GetDynamicParticle()->GetParticleDefinition();
566 566
567 G4bool stable = particleType->GetPDGStab 567 G4bool stable = particleType->GetPDGStable();
568 G4bool candidateForEnd = (endEnergy < fT 568 G4bool candidateForEnd = (endEnergy < fThreshold_Important_Energy)
569 || (fNoLooperTrial 569 || (fNoLooperTrials >= fThresholdTrials) ;
570 G4bool unstableAndKillable = !stable && 570 G4bool unstableAndKillable = !stable && ( fAbandonUnstableTrials != 0);
571 G4bool unstableForEnd = (endEnergy < fTh 571 G4bool unstableForEnd = (endEnergy < fThreshold_Important_Energy)
572 && (fNoLooperTri 572 && (fNoLooperTrials >= fAbandonUnstableTrials) ;
573 if( (candidateForEnd && stable) || (unst 573 if( (candidateForEnd && stable) || (unstableAndKillable && unstableForEnd) )
574 { 574 {
575 // Kill the looping particle 575 // Kill the looping particle
576 // 576 //
577 fParticleChange.ProposeTrackStatus( fS 577 fParticleChange.ProposeTrackStatus( fStopAndKill ) ;
578 G4int particlePDG= particleType->GetPD 578 G4int particlePDG= particleType->GetPDGEncoding();
579 const G4int electronPDG= 11; // G4Elec 579 const G4int electronPDG= 11; // G4Electron::G4Electron()->GetPDGEncoding();
580 580
581 // Simple statistics 581 // Simple statistics
582 fSumEnergyKilled += endEnergy; 582 fSumEnergyKilled += endEnergy;
583 fSumEnerSqKilled += endEnergy * endEne 583 fSumEnerSqKilled += endEnergy * endEnergy;
584 fNumLoopersKilled++; 584 fNumLoopersKilled++;
585 585
586 if( endEnergy > fMaxEnergyKilled ) { 586 if( endEnergy > fMaxEnergyKilled ) {
587 fMaxEnergyKilled = endEnergy; 587 fMaxEnergyKilled = endEnergy;
588 fMaxEnergyKilledPDG = particlePDG; 588 fMaxEnergyKilledPDG = particlePDG;
589 } 589 }
590 if( particleType->GetPDGEncoding() != 590 if( particleType->GetPDGEncoding() != electronPDG )
591 { 591 {
592 fSumEnergyKilled_NonElectron += end 592 fSumEnergyKilled_NonElectron += endEnergy;
593 fSumEnerSqKilled_NonElectron += end 593 fSumEnerSqKilled_NonElectron += endEnergy * endEnergy;
594 fNumLoopersKilled_NonElectron++; 594 fNumLoopersKilled_NonElectron++;
595 595
596 if( endEnergy > fMaxEnergyKilled_No 596 if( endEnergy > fMaxEnergyKilled_NonElectron )
597 { 597 {
598 fMaxEnergyKilled_NonElectron = e 598 fMaxEnergyKilled_NonElectron = endEnergy;
599 fMaxEnergyKilled_NonElecPDG = p 599 fMaxEnergyKilled_NonElecPDG = particlePDG;
600 } 600 }
601 } 601 }
602 602
603 if( endEnergy > fThreshold_Warning_Ene 603 if( endEnergy > fThreshold_Warning_Energy && ! fSilenceLooperWarnings )
604 { 604 {
605 fpLogger->ReportLoopingTrack( track, 605 fpLogger->ReportLoopingTrack( track, stepData, fNoLooperTrials,
606 noCall 606 noCallsASDI, __func__ );
607 } 607 }
608 fNoLooperTrials=0; 608 fNoLooperTrials=0;
609 } 609 }
610 else 610 else
611 { 611 {
612 fMaxEnergySaved = std::max( endEnergy, 612 fMaxEnergySaved = std::max( endEnergy, fMaxEnergySaved);
613 if( fNoLooperTrials == 1 ) { 613 if( fNoLooperTrials == 1 ) {
614 fSumEnergySaved += endEnergy; 614 fSumEnergySaved += endEnergy;
615 if ( !stable ) 615 if ( !stable )
616 fSumEnergyUnstableSaved += endEne 616 fSumEnergyUnstableSaved += endEnergy;
617 } 617 }
618 #ifdef G4VERBOSE 618 #ifdef G4VERBOSE
619 if( verboseLevel > 2 && ! fSilenceLoop 619 if( verboseLevel > 2 && ! fSilenceLooperWarnings )
620 { 620 {
621 G4cout << " " << __func__ 621 G4cout << " " << __func__
622 << " Particle is looping but 622 << " Particle is looping but is saved ..." << G4endl
623 << " Number of trials = " < 623 << " Number of trials = " << fNoLooperTrials << G4endl
624 << " No of calls to = " << 624 << " No of calls to = " << noCallsASDI << G4endl;
625 } 625 }
626 #endif 626 #endif
627 } 627 }
628 } 628 }
629 else 629 else
630 { 630 {
631 fNoLooperTrials=0; 631 fNoLooperTrials=0;
632 } 632 }
633 633
634 // Another (sometimes better way) is to use 634 // Another (sometimes better way) is to use a user-limit maximum Step size
635 // to alleviate this problem .. 635 // to alleviate this problem ..
636 636
637 // Introduce smooth curved trajectories to p 637 // Introduce smooth curved trajectories to particle-change
638 // 638 //
639 fParticleChange.SetPointerToVectorOfAuxiliar 639 fParticleChange.SetPointerToVectorOfAuxiliaryPoints
640 (fFieldPropagator->GimmeTrajectoryVectorAn 640 (fFieldPropagator->GimmeTrajectoryVectorAndForgetIt() );
641 641
642 return &fParticleChange ; 642 return &fParticleChange ;
643 } 643 }
644 644
645 ////////////////////////////////////////////// 645 //////////////////////////////////////////////////////////////////////////
646 // 646 //
647 // This ensures that the PostStep action is a 647 // This ensures that the PostStep action is always called,
648 // so that it can do the relocation if it is 648 // so that it can do the relocation if it is needed.
649 // 649 //
650 650
651 G4double G4Transportation:: 651 G4double G4Transportation::
652 PostStepGetPhysicalInteractionLength( const G4 652 PostStepGetPhysicalInteractionLength( const G4Track&,
653 G4 653 G4double, // previousStepSize
654 G4 654 G4ForceCondition* pForceCond )
655 { 655 {
656 fFieldExertedForce = false; // Not known 656 fFieldExertedForce = false; // Not known
657 *pForceCond = Forced ; 657 *pForceCond = Forced ;
658 return DBL_MAX ; // was kInfinity ; but con 658 return DBL_MAX ; // was kInfinity ; but convention now is DBL_MAX
659 } 659 }
660 660
661 ////////////////////////////////////////////// 661 /////////////////////////////////////////////////////////////////////////////
662 662
663 void G4Transportation::SetTouchableInformation 663 void G4Transportation::SetTouchableInformation(const G4TouchableHandle& touchable)
664 { 664 {
665 const G4VPhysicalVolume* pNewVol = touchable 665 const G4VPhysicalVolume* pNewVol = touchable->GetVolume() ;
666 const G4Material* pNewMaterial = 0 ; 666 const G4Material* pNewMaterial = 0 ;
667 G4VSensitiveDetector* pNewSensitiveDetector 667 G4VSensitiveDetector* pNewSensitiveDetector = 0;
668 668
669 if( pNewVol != 0 ) 669 if( pNewVol != 0 )
670 { 670 {
671 pNewMaterial= pNewVol->GetLogicalVolume()- 671 pNewMaterial= pNewVol->GetLogicalVolume()->GetMaterial();
672 pNewSensitiveDetector= pNewVol->GetLogical 672 pNewSensitiveDetector= pNewVol->GetLogicalVolume()->GetSensitiveDetector();
673 } 673 }
674 674
675 fParticleChange.SetMaterialInTouchable( (G4M 675 fParticleChange.SetMaterialInTouchable( (G4Material *) pNewMaterial ) ;
676 fParticleChange.SetSensitiveDetectorInToucha 676 fParticleChange.SetSensitiveDetectorInTouchable( pNewSensitiveDetector ) ;
677 // temporarily until Get/Set Material of Par 677 // temporarily until Get/Set Material of ParticleChange,
678 // and StepPoint can be made const. 678 // and StepPoint can be made const.
679 679
680 const G4MaterialCutsCouple* pNewMaterialCuts 680 const G4MaterialCutsCouple* pNewMaterialCutsCouple = 0;
681 if( pNewVol != 0 ) 681 if( pNewVol != 0 )
682 { 682 {
683 pNewMaterialCutsCouple=pNewVol->GetLogical 683 pNewMaterialCutsCouple=pNewVol->GetLogicalVolume()->GetMaterialCutsCouple();
684 } 684 }
685 685
686 if ( pNewVol!=0 && pNewMaterialCutsCouple!=0 686 if ( pNewVol!=0 && pNewMaterialCutsCouple!=0
687 && pNewMaterialCutsCouple->GetMaterial()!= 687 && pNewMaterialCutsCouple->GetMaterial()!=pNewMaterial )
688 { 688 {
689 // for parametrized volume 689 // for parametrized volume
690 // 690 //
691 pNewMaterialCutsCouple = 691 pNewMaterialCutsCouple =
692 G4ProductionCutsTable::GetProductionCuts 692 G4ProductionCutsTable::GetProductionCutsTable()
693 ->GetMaterialCuts 693 ->GetMaterialCutsCouple(pNewMaterial,
694 pNewMaterialCut 694 pNewMaterialCutsCouple->GetProductionCuts());
695 } 695 }
696 fParticleChange.SetMaterialCutsCoupleInTouch 696 fParticleChange.SetMaterialCutsCoupleInTouchable( pNewMaterialCutsCouple );
697 697
698 // Set the touchable in ParticleChange 698 // Set the touchable in ParticleChange
699 // this must always be done because the part 699 // this must always be done because the particle change always
700 // uses this value to overwrite the current 700 // uses this value to overwrite the current touchable pointer.
701 // 701 //
702 fParticleChange.SetTouchableHandle(touchable 702 fParticleChange.SetTouchableHandle(touchable) ;
703 } 703 }
704 704
705 ////////////////////////////////////////////// 705 /////////////////////////////////////////////////////////////////////////////
706 // 706 //
707 707
708 G4VParticleChange* G4Transportation::PostStepD 708 G4VParticleChange* G4Transportation::PostStepDoIt( const G4Track& track,
709 709 const G4Step& )
710 { 710 {
711 G4TouchableHandle retCurrentTouchable ; / 711 G4TouchableHandle retCurrentTouchable ; // The one to return
712 G4bool isLastStep= false; 712 G4bool isLastStep= false;
713 713
714 // Initialize ParticleChange (by setting al 714 // Initialize ParticleChange (by setting all its members equal
715 // to correspond 715 // to corresponding members in G4Track)
716 // fParticleChange.Initialize(track) ; // T 716 // fParticleChange.Initialize(track) ; // To initialise TouchableChange
717 717
718 fParticleChange.ProposeTrackStatus(track.Get 718 fParticleChange.ProposeTrackStatus(track.GetTrackStatus()) ;
719 719
720 // If the Step was determined by the volume 720 // If the Step was determined by the volume boundary,
721 // logically relocate the particle 721 // logically relocate the particle
722 722
723 if(fGeometryLimitedStep) 723 if(fGeometryLimitedStep)
724 { 724 {
725 // fCurrentTouchable will now become the p 725 // fCurrentTouchable will now become the previous touchable,
726 // and what was the previous will be freed 726 // and what was the previous will be freed.
727 // (Needed because the preStepPoint can po 727 // (Needed because the preStepPoint can point to the previous touchable)
728 728
729 fLinearNavigator->SetGeometricallyLimitedS 729 fLinearNavigator->SetGeometricallyLimitedStep() ;
730 fLinearNavigator-> 730 fLinearNavigator->
731 LocateGlobalPointAndUpdateTouchableHandle( 731 LocateGlobalPointAndUpdateTouchableHandle( track.GetPosition(),
732 732 track.GetMomentumDirection(),
733 733 fCurrentTouchableHandle,
734 734 true ) ;
735 // Check whether the particle is out of th 735 // Check whether the particle is out of the world volume
736 // If so it has exited and must be killed. 736 // If so it has exited and must be killed.
737 // 737 //
738 if( fCurrentTouchableHandle->GetVolume() = 738 if( fCurrentTouchableHandle->GetVolume() == 0 )
739 { 739 {
740 fParticleChange.ProposeTrackStatus( fSt 740 fParticleChange.ProposeTrackStatus( fStopAndKill ) ;
741 } 741 }
742 retCurrentTouchable = fCurrentTouchableHan 742 retCurrentTouchable = fCurrentTouchableHandle ;
743 fParticleChange.SetTouchableHandle( fCurre 743 fParticleChange.SetTouchableHandle( fCurrentTouchableHandle ) ;
744 744
745 // Update the Step flag which identifies t 745 // Update the Step flag which identifies the Last Step in a volume
746 if( !fFieldExertedForce ) 746 if( !fFieldExertedForce )
747 isLastStep = fLinearNavigator->ExitedM 747 isLastStep = fLinearNavigator->ExitedMotherVolume()
748 || fLinearNavigator->Entered 748 || fLinearNavigator->EnteredDaughterVolume() ;
749 else 749 else
750 isLastStep = fFieldPropagator->IsLastSt 750 isLastStep = fFieldPropagator->IsLastStepInVolume();
751 } 751 }
752 else // fGeometryLimitedStep 752 else // fGeometryLimitedStep is false
753 { 753 {
754 // This serves only to move the Navigator' 754 // This serves only to move the Navigator's location
755 // 755 //
756 fLinearNavigator->LocateGlobalPointWithinV 756 fLinearNavigator->LocateGlobalPointWithinVolume( track.GetPosition() ) ;
757 757
758 // The value of the track's current Toucha 758 // The value of the track's current Touchable is retained.
759 // (and it must be correct because we must 759 // (and it must be correct because we must use it below to
760 // overwrite the (unset) one in particle c 760 // overwrite the (unset) one in particle change)
761 // It must be fCurrentTouchable too ?? 761 // It must be fCurrentTouchable too ??
762 // 762 //
763 fParticleChange.SetTouchableHandle( track. 763 fParticleChange.SetTouchableHandle( track.GetTouchableHandle() ) ;
764 retCurrentTouchable = track.GetTouchableHa 764 retCurrentTouchable = track.GetTouchableHandle() ;
765 765
766 isLastStep= false; 766 isLastStep= false;
767 } // endif ( fGeometryLimitedStep ) 767 } // endif ( fGeometryLimitedStep )
768 fLastStepInVolume= isLastStep; 768 fLastStepInVolume= isLastStep;
769 769
770 fParticleChange.ProposeFirstStepInVolume(fFi 770 fParticleChange.ProposeFirstStepInVolume(fFirstStepInVolume);
771 fParticleChange.ProposeLastStepInVolume(isLa 771 fParticleChange.ProposeLastStepInVolume(isLastStep);
772 772
773 SetTouchableInformation(retCurrentTouchable) 773 SetTouchableInformation(retCurrentTouchable);
774 774
775 return &fParticleChange ; 775 return &fParticleChange ;
776 } 776 }
777 777
778 ////////////////////////////////////////////// 778 /////////////////////////////////////////////////////////////////////////////
779 // New method takes over the responsibility to 779 // New method takes over the responsibility to reset the state of
780 // G4Transportation object at the start of a n 780 // G4Transportation object at the start of a new track or the resumption
781 // of a suspended track. 781 // of a suspended track.
782 // 782 //
783 783
784 void 784 void
785 G4Transportation::StartTracking(G4Track* aTrac 785 G4Transportation::StartTracking(G4Track* aTrack)
786 { 786 {
787 G4VProcess::StartTracking(aTrack); 787 G4VProcess::StartTracking(aTrack);
788 fNewTrack= true; 788 fNewTrack= true;
789 fFirstStepInVolume= true; 789 fFirstStepInVolume= true;
790 fLastStepInVolume= false; 790 fLastStepInVolume= false;
791 791
792 // The actions here are those that were take 792 // The actions here are those that were taken in AlongStepGPIL
793 // when track.GetCurrentStepNumber()==1 793 // when track.GetCurrentStepNumber()==1
794 794
795 // Whether field exists should be determined 795 // Whether field exists should be determined at run level -- TODO
796 G4FieldManagerStore* fieldMgrStore= G4FieldM 796 G4FieldManagerStore* fieldMgrStore= G4FieldManagerStore::GetInstance();
797 fAnyFieldExists = fieldMgrStore->size() > 0; 797 fAnyFieldExists = fieldMgrStore->size() > 0;
798 798
799 // reset safety value and center 799 // reset safety value and center
800 // 800 //
801 fPreviousSafety = 0.0 ; 801 fPreviousSafety = 0.0 ;
802 fPreviousSftOrigin = G4ThreeVector(0.,0.,0.) 802 fPreviousSftOrigin = G4ThreeVector(0.,0.,0.) ;
803 803
804 // reset looping counter -- for motion in fi 804 // reset looping counter -- for motion in field
805 fNoLooperTrials= 0; 805 fNoLooperTrials= 0;
806 // Must clear this state .. else it depends 806 // Must clear this state .. else it depends on last track's value
807 // --> a better solution would set this fro 807 // --> a better solution would set this from state of suspended track TODO ?
808 // Was if( aTrack->GetCurrentStepNumber()==1 808 // Was if( aTrack->GetCurrentStepNumber()==1 ) { .. }
809 809
810 // ChordFinder reset internal state 810 // ChordFinder reset internal state
811 // 811 //
812 if( fFieldPropagator && fAnyFieldExists ) 812 if( fFieldPropagator && fAnyFieldExists )
813 { 813 {
814 fFieldPropagator->ClearPropagatorState(); 814 fFieldPropagator->ClearPropagatorState();
815 // Resets all state of field propagator 815 // Resets all state of field propagator class (ONLY) including safety
816 // values (in case of overlaps and to w 816 // values (in case of overlaps and to wipe for first track).
817 } 817 }
818 818
819 // Make sure to clear the chord finders of a 819 // Make sure to clear the chord finders of all fields (i.e. managers)
820 // 820 //
821 fieldMgrStore->ClearAllChordFindersState(); 821 fieldMgrStore->ClearAllChordFindersState();
822 822
823 // Update the current touchable handle (fro 823 // Update the current touchable handle (from the track's)
824 // 824 //
825 fCurrentTouchableHandle = aTrack->GetTouchab 825 fCurrentTouchableHandle = aTrack->GetTouchableHandle();
826 826
827 // Inform field propagator of new track 827 // Inform field propagator of new track
828 // 828 //
829 fFieldPropagator->PrepareNewTrack(); 829 fFieldPropagator->PrepareNewTrack();
830 } 830 }
831 831
832 ////////////////////////////////////////////// 832 /////////////////////////////////////////////////////////////////////////////
833 // 833 //
834 834
835 G4bool G4Transportation::EnableMagneticMoment( 835 G4bool G4Transportation::EnableMagneticMoment(G4bool useMoment)
836 { 836 {
837 G4bool lastValue= fUseMagneticMoment; 837 G4bool lastValue= fUseMagneticMoment;
838 fUseMagneticMoment= useMoment; 838 fUseMagneticMoment= useMoment;
839 return lastValue; 839 return lastValue;
840 } 840 }
841 841
842 ////////////////////////////////////////////// 842 /////////////////////////////////////////////////////////////////////////////
843 // 843 //
844 844
845 G4bool G4Transportation::EnableGravity(G4bool 845 G4bool G4Transportation::EnableGravity(G4bool useGravity)
846 { 846 {
847 G4bool lastValue= fUseGravity; 847 G4bool lastValue= fUseGravity;
848 fUseGravity= useGravity; 848 fUseGravity= useGravity;
849 return lastValue; 849 return lastValue;
850 } 850 }
851 851
852 ////////////////////////////////////////////// 852 /////////////////////////////////////////////////////////////////////////////
853 // 853 //
854 // Supress (or not) warnings about 'looping' 854 // Supress (or not) warnings about 'looping' particles
855 855
856 void G4Transportation::SetSilenceLooperWarning 856 void G4Transportation::SetSilenceLooperWarnings( G4bool val)
857 { 857 {
858 fSilenceLooperWarnings= val; // Flag to *Su 858 fSilenceLooperWarnings= val; // Flag to *Supress* all 'looper' warnings
859 } 859 }
860 860
861 ////////////////////////////////////////////// 861 /////////////////////////////////////////////////////////////////////////////
862 // 862 //
863 G4bool G4Transportation::GetSilenceLooperWarni 863 G4bool G4Transportation::GetSilenceLooperWarnings()
864 { 864 {
865 return fSilenceLooperWarnings; 865 return fSilenceLooperWarnings;
866 } 866 }
867 867
868 868
869 ////////////////////////////////////////////// 869 /////////////////////////////////////////////////////////////////////////////
870 // 870 //
871 void G4Transportation::SetHighLooperThresholds 871 void G4Transportation::SetHighLooperThresholds()
872 { 872 {
873 // Restores the old high values -- potential 873 // Restores the old high values -- potentially appropriate for energy-frontier
874 // HEP experiments. 874 // HEP experiments.
875 // Caution: All tracks with E < 100 MeV tha 875 // Caution: All tracks with E < 100 MeV that are found to loop are
876 SetThresholdWarningEnergy( 100.0 * CLHEP: 876 SetThresholdWarningEnergy( 100.0 * CLHEP::MeV ); // Warn above this energy
877 SetThresholdImportantEnergy( 250.0 * CLHEP: 877 SetThresholdImportantEnergy( 250.0 * CLHEP::MeV ); // Extra trial above this En
878 878
879 G4int maxTrials = 10; 879 G4int maxTrials = 10;
880 SetThresholdTrials( maxTrials ); 880 SetThresholdTrials( maxTrials );
881 881
882 PushThresholdsToLogger(); // Again, to be s 882 PushThresholdsToLogger(); // Again, to be sure
883 if( verboseLevel ) ReportLooperThresholds() 883 if( verboseLevel ) ReportLooperThresholds();
884 } 884 }
885 885
886 ////////////////////////////////////////////// 886 /////////////////////////////////////////////////////////////////////////////
887 void G4Transportation::SetLowLooperThresholds( 887 void G4Transportation::SetLowLooperThresholds() // Values for low-E applications
888 { 888 {
889 // These values were the default in Geant4 1 889 // These values were the default in Geant4 10.5 - beta
890 SetThresholdWarningEnergy( 1.0 * CLHEP:: 890 SetThresholdWarningEnergy( 1.0 * CLHEP::keV ); // Warn above this En
891 SetThresholdImportantEnergy( 1.0 * CLHEP:: 891 SetThresholdImportantEnergy( 1.0 * CLHEP::MeV ); // Extra trials above it
892 892
893 G4int maxTrials = 30; // A new value - was 893 G4int maxTrials = 30; // A new value - was 10
894 SetThresholdTrials( maxTrials ); 894 SetThresholdTrials( maxTrials );
895 895
896 PushThresholdsToLogger(); // Again, to be s 896 PushThresholdsToLogger(); // Again, to be sure
897 if( verboseLevel ) ReportLooperThresholds() 897 if( verboseLevel ) ReportLooperThresholds();
898 } 898 }
899 899
900 ////////////////////////////////////////////// 900 /////////////////////////////////////////////////////////////////////////////
901 // 901 //
902 void 902 void
903 G4Transportation::ReportMissingLogger( const c 903 G4Transportation::ReportMissingLogger( const char* methodName )
904 { 904 {
905 const char* message= "Logger object missing 905 const char* message= "Logger object missing from G4Transportation object";
906 G4String classAndMethod= G4String("G4Transp 906 G4String classAndMethod= G4String("G4Transportation") + G4String( methodName );
907 G4Exception(classAndMethod, "Missing Logger 907 G4Exception(classAndMethod, "Missing Logger", JustWarning, message);
908 } 908 }
909 909
910 910
911 ////////////////////////////////////////////// 911 /////////////////////////////////////////////////////////////////////////////
912 // 912 //
913 void 913 void
914 G4Transportation::ReportLooperThresholds() 914 G4Transportation::ReportLooperThresholds()
915 { 915 {
916 PushThresholdsToLogger(); // To be absolut 916 PushThresholdsToLogger(); // To be absolutely certain they are in sync
917 fpLogger->ReportLooperThresholds("G4Transpo 917 fpLogger->ReportLooperThresholds("G4Transportation");
918 } 918 }
919 919
920 ////////////////////////////////////////////// 920 /////////////////////////////////////////////////////////////////////////////
921 // 921 //
922 void G4Transportation::ProcessDescription(std: 922 void G4Transportation::ProcessDescription(std::ostream& outStr) const
923 923
924 // StreamInfo(std::ostream& out, const G4Parti 924 // StreamInfo(std::ostream& out, const G4ParticleDefinition& part, G4bool rst) const
925 925
926 { 926 {
927 G4String indent = " "; // : ""); 927 G4String indent = " "; // : "");
928 G4long oldPrec= outStr.precision(6); 928 G4long oldPrec= outStr.precision(6);
929 // outStr << std::setprecision(6); 929 // outStr << std::setprecision(6);
930 outStr << G4endl << indent << GetProcessName 930 outStr << G4endl << indent << GetProcessName() << ": ";
931 931
932 outStr << " Parameters for looping particl 932 outStr << " Parameters for looping particles: " << G4endl
933 << " warning-E = " << fThreshold_ 933 << " warning-E = " << fThreshold_Warning_Energy / CLHEP::MeV << " MeV " << G4endl
934 << " important E = " << fThreshol 934 << " important E = " << fThreshold_Important_Energy / CLHEP::MeV << " MeV " << G4endl
935 << " thresholdTrials " << fThresh 935 << " thresholdTrials " << fThresholdTrials << G4endl;
936 outStr.precision(oldPrec); 936 outStr.precision(oldPrec);
937 } 937 }
938 938