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