| 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 // $Id: G4Transportation.cc,v 1.55.2.1 2006/06/29 21:12:28 gunter Exp $
>> 28 // GEANT4 tag $Name: geant4-08-01 $
27 // 29 //
28 // ------------------------------------------- 30 // ------------------------------------------------------------
29 // GEANT 4 include file implementation 31 // GEANT 4 include file implementation
30 // 32 //
31 // ------------------------------------------- 33 // ------------------------------------------------------------
32 // 34 //
33 // This class is a process responsible for the 35 // This class is a process responsible for the transportation of
34 // a particle, ie the geometrical propagation 36 // a particle, ie the geometrical propagation that encounters the
35 // geometrical sub-volumes of the detectors. 37 // geometrical sub-volumes of the detectors.
36 // 38 //
37 // It is also tasked with the key role of prop << 39 // It is also tasked with part of updating the "safety".
38 // which will be used to update the post-ste <<
39 // 40 //
40 // =========================================== 41 // =======================================================================
>> 42 // Modified:
>> 43 // 19 Jan 2006, P.MoraDeFreitas: Fix for suspended tracks (StartTracking)
>> 44 // 11 Aug 2004, M.Asai: Add G4VSensitiveDetector* for updating stepPoint.
>> 45 // 21 June 2003, J.Apostolakis: Calling field manager with
>> 46 // track, to enable it to configure its accuracy
>> 47 // 13 May 2003, J.Apostolakis: Zero field areas now taken into
>> 48 // account correclty in all cases (thanks to W Pokorski).
>> 49 // 29 June 2001, J.Apostolakis, D.Cote-Ahern, P.Gumplinger:
>> 50 // correction for spin tracking
>> 51 // 20 Febr 2001, J.Apostolakis: update for new FieldTrack
>> 52 // 22 Sept 2000, V.Grichine: update of Kinetic Energy
>> 53 // 9 June 1999, J.Apostolakis & S.Giani: protect full relocation
>> 54 // in DEBUG for track starting on surface that
>> 55 // goes step < tolerance.
41 // Created: 19 March 1997, J. Apostolakis 56 // Created: 19 March 1997, J. Apostolakis
42 // =========================================== 57 // =======================================================================
43 58
44 #include "G4Transportation.hh" 59 #include "G4Transportation.hh"
45 #include "G4TransportationProcessType.hh" <<
46 #include "G4TransportationLogger.hh" <<
47 <<
48 #include "G4PhysicalConstants.hh" <<
49 #include "G4SystemOfUnits.hh" <<
50 #include "G4ProductionCutsTable.hh" 60 #include "G4ProductionCutsTable.hh"
51 #include "G4ParticleTable.hh" 61 #include "G4ParticleTable.hh"
52 << 62 #include "G4ChordFinder.hh"
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; 63 class G4VSensitiveDetector;
65 64
66 G4bool G4Transportation::fUseMagneticMoment=fa <<
67 G4bool G4Transportation::fUseGravity= false; <<
68 G4bool G4Transportation::fSilenceLooperWarning <<
69 <<
70 ////////////////////////////////////////////// 65 //////////////////////////////////////////////////////////////////////////
71 // 66 //
72 // Constructor 67 // Constructor
73 68
74 G4Transportation::G4Transportation( G4int verb << 69 G4Transportation::G4Transportation( G4int verboseLevel )
75 : G4VProcess( aName, fTransportation ), << 70 : G4VProcess( G4String("Transportation"), fTransportation ),
76 fFieldExertedForce( false ), << 71 fParticleIsLooping( false ),
77 fPreviousSftOrigin( 0.,0.,0. ), << 72 fPreviousSftOrigin (0.,0.,0.),
78 fPreviousSafety( 0.0 ), << 73 fPreviousSafety ( 0.0 ),
79 fEndPointDistance( -1.0 ), << 74 fThreshold_Warning_Energy( 100 * MeV ),
80 fShortStepOptimisation( false ) // Old def << 75 fThreshold_Important_Energy( 250 * MeV ),
>> 76 fThresholdTrials( 10 ),
>> 77 fUnimportant_Energy( 1 * MeV ),
>> 78 fNoLooperTrials(0),
>> 79 fSumEnergyKilled( 0.0 ), fMaxEnergyKilled( 0.0 ),
>> 80 fVerboseLevel( verboseLevel )
81 { 81 {
82 SetProcessSubType(static_cast<G4int>(TRANSPO <<
83 pParticleChange= &fParticleChange; // Requ <<
84 SetVerboseLevel(verbosity); <<
85 <<
86 G4TransportationManager* transportMgr ; 82 G4TransportationManager* transportMgr ;
87 83
88 transportMgr = G4TransportationManager::GetT 84 transportMgr = G4TransportationManager::GetTransportationManager() ;
89 85
90 fLinearNavigator = transportMgr->GetNavigato 86 fLinearNavigator = transportMgr->GetNavigatorForTracking() ;
91 87
92 fFieldPropagator = transportMgr->GetPropagat << 88 // fGlobalFieldMgr = transportMgr->GetFieldManager() ;
93 <<
94 fpSafetyHelper = transportMgr->GetSafetyHe <<
95 <<
96 fpLogger = new G4TransportationLogger("G4Tra <<
97 89
98 if( G4TransportationParameters::Exists() ) << 90 fFieldPropagator = transportMgr->GetPropagatorInField() ;
99 { <<
100 auto trParams= G4TransportationParameters: <<
101 <<
102 SetThresholdWarningEnergy( trParams->GetW <<
103 SetThresholdImportantEnergy( trParams->Get <<
104 SetThresholdTrials( trParams->GetNumberOfT <<
105 G4Transportation::fSilenceLooperWarnings= <<
106 } <<
107 else { <<
108 SetHighLooperThresholds(); <<
109 // Use the old defaults: Warning = 100 Me <<
110 } <<
111 91
112 PushThresholdsToLogger(); << 92 // Cannot determine whether a field exists here,
113 // Should be done by Set methods in SetHighL << 93 // because it would only work if the field manager has informed
114 << 94 // about the detector's field before this transportation process
115 static G4ThreadLocal G4TouchableHandle* pNul << 95 // is constructed.
116 if ( !pNullTouchableHandle) << 96 // Instead later the method DoesGlobalFieldExist() is called
117 { <<
118 pNullTouchableHandle = new G4TouchableHand <<
119 } <<
120 fCurrentTouchableHandle = *pNullTouchableHan <<
121 // Points to (G4VTouchable*) 0 <<
122 97
>> 98 fCurrentTouchableHandle = new G4TouchableHistory();
123 99
124 #ifdef G4VERBOSE << 100 fEndGlobalTimeComputed = false;
125 if( verboseLevel > 0) << 101 fCandidateEndGlobalTime = 0;
126 { <<
127 G4cout << " G4Transportation constructor> <<
128 if ( fShortStepOptimisation ) { G4cout < <<
129 else { G4cout < <<
130 } <<
131 #endif <<
132 } 102 }
133 103
134 ////////////////////////////////////////////// 104 //////////////////////////////////////////////////////////////////////////
135 105
136 G4Transportation::~G4Transportation() 106 G4Transportation::~G4Transportation()
137 { 107 {
138 if( fSumEnergyKilled > 0.0 ) << 108 if( (fVerboseLevel > 0) && (fSumEnergyKilled > 0.0 ) ){
139 { << 109 G4cout << " G4Transportation: Statistics for looping particles " << G4endl;
140 PrintStatistics( G4cout ); << 110 G4cout << " Sum of energy of loopers killed: " << fSumEnergyKilled << G4endl;
141 } << 111 G4cout << " Max energy of loopers killed: " << fMaxEnergyKilled << G4endl;
142 delete fpLogger; << 112 }
143 } <<
144 <<
145 ////////////////////////////////////////////// <<
146 <<
147 void <<
148 G4Transportation::PrintStatistics( std::ostrea <<
149 { <<
150 outStr << " G4Transportation: Statistics fo <<
151 if( fSumEnergyKilled > 0.0 || fNumLoopersKi <<
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 } 113 }
182 114
183 ////////////////////////////////////////////// 115 //////////////////////////////////////////////////////////////////////////
184 // 116 //
185 // Responsibilities: 117 // Responsibilities:
186 // Find whether the geometry limits the Ste 118 // Find whether the geometry limits the Step, and to what length
187 // Calculate the new value of the safety an 119 // Calculate the new value of the safety and return it.
188 // Store the final time, position and momen 120 // Store the final time, position and momentum.
189 121
190 G4double G4Transportation::AlongStepGetPhysica << 122 G4double G4Transportation::
191 const G4Track& track, << 123 AlongStepGetPhysicalInteractionLength( const G4Track& track,
192 G4double, // previousStepSize << 124 G4double, // previousStepSize
193 G4double currentMinimumStep, G4double& curre << 125 G4double currentMinimumStep,
194 G4GPILSelection* selection) << 126 G4double& currentSafety,
>> 127 G4GPILSelection* selection )
195 { 128 {
196 // Initial actions moved to StartTrack() << 129 G4double geometryStepLength, newSafety ;
>> 130 fParticleIsLooping = false ;
>> 131
>> 132 // Initial actions moved to StartTrack()
197 // -------------------------------------- 133 // --------------------------------------
198 // Note: in case another process changes tou 134 // Note: in case another process changes touchable handle
199 // it will be necessary to add here (for << 135 // it will be necessary to add here (for all steps)
200 // fCurrentTouchableHandle = aTrack->GetTouc 136 // fCurrentTouchableHandle = aTrack->GetTouchableHandle();
201 137
202 // GPILSelection is set to defaule value of 138 // GPILSelection is set to defaule value of CandidateForSelection
203 // It is a return value 139 // It is a return value
204 // 140 //
205 *selection = CandidateForSelection; << 141 *selection = CandidateForSelection ;
206 142
207 // Get initial Energy/Momentum of the track 143 // Get initial Energy/Momentum of the track
208 // 144 //
209 const G4ThreeVector startPosition = track << 145 const G4DynamicParticle* pParticle = track.GetDynamicParticle() ;
210 const G4ThreeVector startMomentumDir = track << 146 const G4ParticleDefinition* pParticleDef = pParticle->GetDefinition() ;
>> 147 G4ThreeVector startMomentumDir = pParticle->GetMomentumDirection() ;
>> 148 G4ThreeVector startPosition = track.GetPosition() ;
>> 149
>> 150 // G4double theTime = track.GetGlobalTime() ;
211 151
212 // The Step Point safety can be limited by o << 152 // The Step Point safety can be limited by other geometries and/or the
213 // assumptions of any process - it's not alw 153 // assumptions of any process - it's not always the geometrical safety.
214 // We calculate the starting point's isotrop 154 // We calculate the starting point's isotropic safety here.
>> 155 //
>> 156 G4ThreeVector OriginShift = startPosition - fPreviousSftOrigin ;
>> 157 G4double MagSqShift = OriginShift.mag2() ;
>> 158 if( MagSqShift >= sqr(fPreviousSafety) )
215 { 159 {
216 const G4double MagSqShift = (startPosition << 160 currentSafety = 0.0 ;
217 << 161 }
218 if(MagSqShift >= sqr(fPreviousSafety)) << 162 else
219 currentSafety = 0.0; << 163 {
220 else << 164 currentSafety = fPreviousSafety - std::sqrt(MagSqShift) ;
221 currentSafety = fPreviousSafety - std::s <<
222 } 165 }
223 166
224 // Is the particle charged or has it a magne << 167 // Is the particle charged ?
225 // 168 //
226 const G4DynamicParticle* pParticle = track.G << 169 G4double particleCharge = pParticle->GetCharge() ;
227 <<
228 const G4double particleMass = pParticle->G <<
229 const G4double particleCharge = pParticle->G <<
230 const G4double kineticEnergy = pParticle->Ge <<
231 170
232 const G4double magneticMoment = pParticle << 171 fGeometryLimitedStep = false ;
233 const G4ThreeVector particleSpin = pParticle << 172 // fEndGlobalTimeComputed = false ;
234 173
235 // There is no need to locate the current vo 174 // There is no need to locate the current volume. It is Done elsewhere:
236 // On track construction << 175 // On track construction
237 // By the tracking, after all AlongStepDoI 176 // By the tracking, after all AlongStepDoIts, in "Relocation"
238 177
239 // Check if the particle has a force, EM or << 178 // Check whether the particle have an (EM) field force exerting upon it
240 // 179 //
>> 180 G4FieldManager* fieldMgr=0;
>> 181 G4bool fieldExertsForce = false ;
>> 182 if( (particleCharge != 0.0) )
>> 183 {
>> 184 fieldMgr= fFieldPropagator->FindAndSetFieldManager( track.GetVolume() );
>> 185 if (fieldMgr != 0) {
>> 186 // If the field manager has no field, there is no field !
>> 187 fieldExertsForce = (fieldMgr->GetDetectorField() != 0);
>> 188 }
>> 189 }
>> 190
>> 191 // Choose the calculation of the transportation: Field or not
>> 192 //
>> 193 if( !fieldExertsForce )
>> 194 {
>> 195 G4double linearStepLength ;
>> 196 if( currentMinimumStep <= currentSafety )
>> 197 {
>> 198 // The Step is guaranteed to be taken
>> 199 //
>> 200 geometryStepLength = currentMinimumStep ;
>> 201 fGeometryLimitedStep = false ;
>> 202 }
>> 203 else
>> 204 {
>> 205 // Find whether the straight path intersects a volume
>> 206 //
>> 207 linearStepLength = fLinearNavigator->ComputeStep( startPosition,
>> 208 startMomentumDir,
>> 209 currentMinimumStep,
>> 210 newSafety) ;
>> 211 // Remember last safety origin & value.
>> 212 //
>> 213 fPreviousSftOrigin = startPosition ;
>> 214 fPreviousSafety = newSafety ;
>> 215
>> 216 // The safety at the initial point has been re-calculated:
>> 217 //
>> 218 currentSafety = newSafety ;
>> 219
>> 220 if( linearStepLength <= currentMinimumStep)
>> 221 {
>> 222 // The geometry limits the Step size (an intersection was found.)
>> 223 //
>> 224 geometryStepLength = linearStepLength ;
>> 225 fGeometryLimitedStep = true ;
>> 226 }
>> 227 else
>> 228 {
>> 229 // The full Step is taken.
>> 230 //
>> 231 geometryStepLength = currentMinimumStep ;
>> 232 fGeometryLimitedStep = false ;
>> 233 }
>> 234 }
>> 235 endpointDistance = geometryStepLength ;
241 236
242 G4bool eligibleEM = << 237 // Calculate final position
243 (particleCharge != 0.0) || ((magneticMomen << 238 //
244 G4bool eligibleGrav = (particleMass != 0.0) << 239 fTransportEndPosition = startPosition+geometryStepLength*startMomentumDir ;
245 <<
246 fFieldExertedForce = false; <<
247 <<
248 if(eligibleEM || eligibleGrav) <<
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 <<
266 G4double geometryStepLength = currentMinimum <<
267 <<
268 if(currentMinimumStep == 0.0) <<
269 { <<
270 fEndPointDistance = 0.0; <<
271 fGeometryLimitedStep = false; // Old cod <<
272 // Changed to avoid problems when setting <<
273 fMomentumChanged = false; <<
274 fParticleIsLooping = false; <<
275 fEndGlobalTimeComputed = false; <<
276 fTransportEndPosition = startPosition <<
277 fTransportEndMomentumDir = startMomentum <<
278 fTransportEndKineticEnergy = kineticEnergy <<
279 fTransportEndSpin = particleSpin; <<
280 } <<
281 else if(!fFieldExertedForce) <<
282 { <<
283 fGeometryLimitedStep = false; <<
284 if(geometryStepLength > currentSafety || ! <<
285 { <<
286 const G4double linearStepLength = fLinea <<
287 startPosition, startMomentumDir, curre <<
288 <<
289 if(linearStepLength <= currentMinimumSte <<
290 { <<
291 geometryStepLength = linearStepLength; <<
292 fGeometryLimitedStep = true; <<
293 } <<
294 // Remember last safety origin & value. <<
295 // <<
296 fPreviousSftOrigin = startPosition; <<
297 fPreviousSafety = currentSafety; <<
298 fpSafetyHelper->SetCurrentSafety(current <<
299 } <<
300 <<
301 fEndPointDistance = geometryStepLength; <<
302 <<
303 fMomentumChanged = false; <<
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 240
348 fGeometryLimitedStep = fFieldPropagator->I << 241 // Momentum direction, energy and polarisation are unchanged by transport
349 // << 242 //
350 // It is possible that step was reduced in << 243 fTransportEndMomentumDir = startMomentumDir ;
351 // previous zero steps. To cope with case << 244 fTransportEndKineticEnergy = track.GetKineticEnergy() ;
352 // in full, we must rely on PiF to obtain << 245 fTransportEndSpin = track.GetPolarization();
>> 246 fParticleIsLooping = false ;
>> 247 fMomentumChanged = false ;
>> 248 fEndGlobalTimeComputed = false ;
>> 249 }
>> 250 else // A field exerts force
>> 251 {
>> 252 G4double momentumMagnitude = pParticle->GetTotalMomentum() ;
>> 253 G4ThreeVector EndUnitMomentum ;
>> 254 G4double lengthAlongCurve ;
>> 255 G4double restMass = pParticleDef->GetPDGMass() ;
>> 256
>> 257 fFieldPropagator->SetChargeMomentumMass( particleCharge, // in e+ units
>> 258 momentumMagnitude, // in Mev/c
>> 259 restMass ) ;
353 260
354 G4bool changesEnergy = << 261 // Message the field Manager, to configure it for this track
355 fFieldPropagator->GetCurrentFieldManager << 262 fieldMgr->ConfigureForTrack( &track );
356 263
357 fMomentumChanged = true; << 264 G4ThreeVector spin = track.GetPolarization() ;
358 fParticleIsLooping = fFieldPropagator->IsP << 265 G4FieldTrack aFieldTrack = G4FieldTrack( startPosition,
>> 266 track.GetMomentumDirection(),
>> 267 0.0,
>> 268 track.GetKineticEnergy(),
>> 269 restMass,
>> 270 track.GetVelocity(),
>> 271 track.GetGlobalTime(), // Lab.
>> 272 track.GetProperTime(), // Part.
>> 273 &spin ) ;
>> 274 if( currentMinimumStep > 0 )
>> 275 {
>> 276 // Do the Transport in the field (non recti-linear)
>> 277 //
>> 278 lengthAlongCurve = fFieldPropagator->ComputeStep( aFieldTrack,
>> 279 currentMinimumStep,
>> 280 currentSafety,
>> 281 track.GetVolume() ) ;
>> 282 if( lengthAlongCurve < currentMinimumStep)
>> 283 {
>> 284 geometryStepLength = lengthAlongCurve ;
>> 285 fGeometryLimitedStep = true ;
>> 286 }
>> 287 else
>> 288 {
>> 289 geometryStepLength = currentMinimumStep ;
>> 290 fGeometryLimitedStep = false ;
>> 291 }
>> 292 }
>> 293 else
>> 294 {
>> 295 geometryStepLength = lengthAlongCurve= 0.0 ;
>> 296 fGeometryLimitedStep = false ;
>> 297 }
359 298
360 fEndGlobalTimeComputed = changesEnergy; << 299 // Remember last safety origin & value.
361 fTransportEndPosition = aFieldTrack.Get << 300 //
362 fTransportEndMomentumDir = aFieldTrack.Get << 301 fPreviousSftOrigin = startPosition ;
>> 302 fPreviousSafety = currentSafety ;
>> 303
>> 304 // Get the End-Position and End-Momentum (Dir-ection)
>> 305 //
>> 306 fTransportEndPosition = aFieldTrack.GetPosition() ;
363 307
364 // G4cout << " G4Transport: End of step pM << 308 // Momentum: Magnitude and direction can be changed too now ...
>> 309 //
>> 310 fMomentumChanged = true ;
>> 311 fTransportEndMomentumDir = aFieldTrack.GetMomentumDir() ;
365 312
366 fEndPointDistance = (fTransportEndPositio << 313 fTransportEndKineticEnergy = aFieldTrack.GetKineticEnergy() ;
367 314
368 // Ignore change in energy for fields that << 315 if( fFieldPropagator->GetCurrentFieldManager()->DoesFieldChangeEnergy() )
369 // This hides the integration error, but g << 316 {
370 fTransportEndKineticEnergy = << 317 // If the field can change energy, then the time must be integrated
371 changesEnergy ? aFieldTrack.GetKineticEn << 318 // - so this should have been updated
372 fTransportEndSpin = aFieldTrack.GetSpin(); << 319 //
>> 320 fCandidateEndGlobalTime = aFieldTrack.GetLabTimeOfFlight();
>> 321 fEndGlobalTimeComputed = true;
373 322
374 if(fEndGlobalTimeComputed) << 323 // was ( fCandidateEndGlobalTime != track.GetGlobalTime() );
375 { << 324 // a cleaner way is to have FieldTrack knowing whether time is updated.
376 // If the field can change energy, then << 325 }
377 // - so this should have been updated << 326 else
378 // << 327 {
379 fCandidateEndGlobalTime = aFieldTrack.Ge << 328 // The energy should be unchanged by field transport,
>> 329 // - so the time changed will be calculated elsewhere
>> 330 //
>> 331 fEndGlobalTimeComputed = false;
380 332
381 // was ( fCandidateEndGlobalTime != trac << 333 // Check that the integration preserved the energy
382 // a cleaner way is to have FieldTrack k << 334 // - and if not correct this!
383 } << 335 G4double startEnergy= track.GetKineticEnergy();
384 #if defined(G4VERBOSE) || defined(G4DEBUG_TRAN << 336 G4double endEnergy= fTransportEndKineticEnergy;
385 else << 337
386 { << 338 static G4int no_inexact_steps=0, no_large_ediff;
387 // The energy should be unchanged by fie << 339 G4double absEdiff = std::fabs(startEnergy- endEnergy);
388 // - so the time changed will be calc << 340 if( absEdiff > perMillion * endEnergy )
389 // << 341 {
390 // Check that the integration preserved << 342 no_inexact_steps++;
391 // - and if not correct this! << 343 // Possible statistics keeping here ...
392 G4double startEnergy = kineticEnergy; << 344 }
393 G4double endEnergy = fTransportEndKine << 345 if( fVerboseLevel > 1 )
394 <<
395 static G4ThreadLocal G4int no_large_edif <<
396 if(verboseLevel > 1) <<
397 { <<
398 if(std::fabs(startEnergy - endEnergy) <<
399 { 346 {
400 static G4ThreadLocal G4int no_warnin << 347 if( std::fabs(startEnergy- endEnergy) > perThousand * endEnergy )
401 moduloFac <<
402 no_large_ediff++; <<
403 if((no_large_ediff % warnModulo) == <<
404 { 348 {
405 no_warnings++; << 349 static G4int no_warnings= 0, warnModulo=1, moduloFactor= 10;
406 std::ostringstream message; << 350 no_large_ediff ++;
407 message << "Energy change in Step << 351 if( (no_large_ediff% warnModulo) == 0 )
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 { 352 {
435 warnModulo *= moduloFactor; << 353 no_warnings++;
>> 354 G4cout << "WARNING - G4Transportation::AlongStepGetPIL()" << G4endl
>> 355 << " Energy changed in Step, more than 1/1000: " << G4endl
>> 356 << " Start= " << startEnergy << G4endl
>> 357 << " End= " << endEnergy << G4endl
>> 358 << " Relative change= "
>> 359 << (startEnergy-endEnergy)/startEnergy << G4endl;
>> 360 G4cout << " Energy has been corrected -- however, review"
>> 361 << " field propagation parameters for accuracy." << G4endl;
>> 362 G4cerr << "ERROR - G4Transportation::AlongStepGetPIL()" << G4endl
>> 363 << " Bad 'endpoint'. Energy change detected"
>> 364 << " and corrected," << G4endl
>> 365 << " occurred already "
>> 366 << no_large_ediff << " times." << G4endl;
>> 367 if( no_large_ediff == warnModulo * moduloFactor )
>> 368 {
>> 369 warnModulo *= moduloFactor;
>> 370 }
436 } 371 }
437 } 372 }
438 } << 373 } // end of if (fVerboseLevel)
439 } // end of if (verboseLevel) << 374
440 } << 375 // Correct the energy for fields that conserve it
441 #endif << 376 // This - hides the integration error
>> 377 // - but gives a better physical answer
>> 378 fTransportEndKineticEnergy= track.GetKineticEnergy();
>> 379 }
>> 380
>> 381 fTransportEndSpin = aFieldTrack.GetSpin();
>> 382 fParticleIsLooping = fFieldPropagator->IsParticleLooping() ;
>> 383 endpointDistance = (fTransportEndPosition - startPosition).mag() ;
>> 384 }
>> 385
>> 386 // If we are asked to go a step length of 0, and we are on a boundary
>> 387 // then a boundary will also limit the step -> we must flag this.
>> 388 //
>> 389 if( currentMinimumStep == 0.0 )
>> 390 {
>> 391 if( currentSafety == 0.0 ) fGeometryLimitedStep = true ;
442 } 392 }
443 393
444 // Update the safety starting from the end-p 394 // Update the safety starting from the end-point,
445 // if it will become negative at the end-poi 395 // if it will become negative at the end-point.
446 // 396 //
447 if(currentSafety < fEndPointDistance) << 397 if( currentSafety < endpointDistance )
448 { 398 {
449 if(particleCharge != 0.0) <<
450 { <<
451 G4double endSafety = 399 G4double endSafety =
452 fLinearNavigator->ComputeSafety(fTrans << 400 fLinearNavigator->ComputeSafety( fTransportEndPosition) ;
453 currentSafety = endSafety; << 401 currentSafety = endSafety ;
454 fPreviousSftOrigin = fTransportEndPositi << 402 fPreviousSftOrigin = fTransportEndPosition ;
455 fPreviousSafety = currentSafety; << 403 fPreviousSafety = currentSafety ;
456 fpSafetyHelper->SetCurrentSafety(current <<
457 404
458 // Because the Stepping Manager assumes << 405 // Because the Stepping Manager assumes it is from the start point,
459 // add the StepLength 406 // add the StepLength
460 // 407 //
461 currentSafety += fEndPointDistance; << 408 currentSafety += endpointDistance ;
462 409
463 #ifdef G4DEBUG_TRANSPORT << 410 #ifdef G4DEBUG_TRANSPORT
464 G4cout.precision(12); << 411 G4cout.precision(16) ;
465 G4cout << "***G4Transportation::AlongSte << 412 G4cout << "***Transportation::AlongStepGPIL ** " << G4endl ;
466 G4cout << " Called Navigator->ComputeSa << 413 G4cout << " Called Navigator->ComputeSafety at "
467 << " and it returned safety= << 414 << fTransportEndPosition
468 G4cout << " Adding endpoint distance << 415 << " and it returned safety= " << endSafety << G4endl ;
469 << " to obtain pseudo-safety= << 416 G4cout << " Adding endpoint distance " << endpointDistance
470 } << 417 << " we obtain pseudo-safety= " << currentSafety << G4endl ;
471 else <<
472 { <<
473 G4cout << "***G4Transportation::AlongSte <<
474 G4cout << " Avoiding call to ComputeSaf <<
475 G4cout << " charge = " << particl <<
476 G4cout << " mag moment = " << magneti <<
477 #endif 418 #endif
478 } << 419 }
479 } <<
480 420
481 fFirstStepInVolume = fNewTrack || fLastStepI << 421 fParticleChange.ProposeTrueStepLength(geometryStepLength) ;
482 fLastStepInVolume = false; <<
483 fNewTrack = false; <<
484 422
485 fParticleChange.ProposeFirstStepInVolume(fFi << 423 return geometryStepLength ;
486 fParticleChange.ProposeTrueStepLength(geomet <<
487 <<
488 return geometryStepLength; <<
489 } 424 }
490 425
491 ////////////////////////////////////////////// 426 //////////////////////////////////////////////////////////////////////////
492 // 427 //
493 // Initialize ParticleChange (by setting al 428 // Initialize ParticleChange (by setting all its members equal
494 // to correspond 429 // to corresponding members in G4Track)
495 430
496 G4VParticleChange* G4Transportation::AlongStep 431 G4VParticleChange* G4Transportation::AlongStepDoIt( const G4Track& track,
497 432 const G4Step& stepData )
498 { 433 {
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) ; 434 fParticleChange.Initialize(track) ;
512 435
513 // Code for specific process 436 // Code for specific process
514 // 437 //
515 fParticleChange.ProposePosition(fTransportEn 438 fParticleChange.ProposePosition(fTransportEndPosition) ;
516 fParticleChange.ProposeMomentumDirection(fTr 439 fParticleChange.ProposeMomentumDirection(fTransportEndMomentumDir) ;
517 fParticleChange.ProposeEnergy(fTransportEndK 440 fParticleChange.ProposeEnergy(fTransportEndKineticEnergy) ;
518 fParticleChange.SetMomentumChanged(fMomentum 441 fParticleChange.SetMomentumChanged(fMomentumChanged) ;
519 442
520 fParticleChange.ProposePolarization(fTranspo 443 fParticleChange.ProposePolarization(fTransportEndSpin);
521 << 444
522 G4double deltaTime = 0.0 ; 445 G4double deltaTime = 0.0 ;
523 446
524 // Calculate Lab Time of Flight (ONLY if fi 447 // Calculate Lab Time of Flight (ONLY if field Equations used it!)
525 // G4double endTime = fCandidateEndGlobalT << 448 // G4double endTime = fCandidateEndGlobalTime;
526 // G4double delta_time = endTime - startTime << 449 // G4double delta_time = endTime - startTime;
527 450
528 G4double startTime = track.GetGlobalTime() ; 451 G4double startTime = track.GetGlobalTime() ;
529 452
530 if (!fEndGlobalTimeComputed) 453 if (!fEndGlobalTimeComputed)
531 { 454 {
532 // The time was not integrated .. make th 455 // The time was not integrated .. make the best estimate possible
533 // 456 //
534 G4double initialVelocity = stepData.GetPr << 457 G4double finalVelocity = track.GetVelocity() ;
535 G4double stepLength = track.GetStepL << 458 G4double initialVelocity = stepData.GetPreStepPoint()->GetVelocity() ;
536 << 459 G4double stepLength = track.GetStepLength() ;
537 deltaTime= 0.0; // in case initialVeloci << 460
538 if ( initialVelocity > 0.0 ) { deltaTime << 461 static const G4ParticleDefinition* fOpticalPhoton =
539 << 462 G4ParticleTable::GetParticleTable()->FindParticle("opticalphoton");
>> 463 const G4DynamicParticle* fpDynamicParticle = track.GetDynamicParticle();
>> 464 if (fpDynamicParticle->GetDefinition()== fOpticalPhoton)
>> 465 {
>> 466 // A photon is in the medium of the final point
>> 467 // during the step, so it has the final velocity.
>> 468 deltaTime = stepLength/finalVelocity ;
>> 469 }
>> 470 else if (finalVelocity > 0.0)
>> 471 {
>> 472 G4double meanInverseVelocity ;
>> 473 // deltaTime = stepLength/finalVelocity ;
>> 474 meanInverseVelocity = 0.5
>> 475 * ( 1.0 / initialVelocity + 1.0 / finalVelocity ) ;
>> 476 deltaTime = stepLength * meanInverseVelocity ;
>> 477 }
>> 478 else
>> 479 {
>> 480 deltaTime = stepLength/initialVelocity ;
>> 481 }
540 fCandidateEndGlobalTime = startTime + d 482 fCandidateEndGlobalTime = startTime + deltaTime ;
541 fParticleChange.ProposeLocalTime( track. <<
542 } 483 }
543 else 484 else
544 { 485 {
545 deltaTime = fCandidateEndGlobalTime - sta 486 deltaTime = fCandidateEndGlobalTime - startTime ;
546 fParticleChange.ProposeGlobalTime( fCandi <<
547 } 487 }
548 488
>> 489 fParticleChange.ProposeGlobalTime( fCandidateEndGlobalTime ) ;
549 490
550 // Now Correct by Lorentz factor to get delt << 491 // Now Correct by Lorentz factor to get "proper" deltaTime
551 492
552 G4double restMass = track.GetDynamicP 493 G4double restMass = track.GetDynamicParticle()->GetMass() ;
553 G4double deltaProperTime = deltaTime*( restM 494 G4double deltaProperTime = deltaTime*( restMass/track.GetTotalEnergy() ) ;
554 495
555 fParticleChange.ProposeProperTime(track.GetP 496 fParticleChange.ProposeProperTime(track.GetProperTime() + deltaProperTime) ;
556 //fParticleChange.ProposeTrueStepLength( tra << 497 //fParticleChange. ProposeTrueStepLength( track.GetStepLength() ) ;
557 498
558 // If the particle is caught looping or is s 499 // If the particle is caught looping or is stuck (in very difficult
559 // boundaries) in a magnetic field (doing ma << 500 // boundaries) in a magnetic field (doing many steps)
>> 501 // THEN this kills it ...
560 // 502 //
561 if ( fParticleIsLooping ) 503 if ( fParticleIsLooping )
562 { 504 {
563 G4double endEnergy= fTransportEndKinetic 505 G4double endEnergy= fTransportEndKineticEnergy;
564 fNoLooperTrials ++; <<
565 auto particleType= track.GetDynamicParti <<
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 506
603 if( endEnergy > fThreshold_Warning_Ene << 507 if( (endEnergy < fThreshold_Important_Energy)
604 { << 508 || (fNoLooperTrials >= fThresholdTrials ) ){
605 fpLogger->ReportLoopingTrack( track, << 509 // Kill the looping particle
606 noCall << 510 //
607 } << 511 fParticleChange.ProposeTrackStatus( fStopAndKill ) ;
608 fNoLooperTrials=0; << 512
609 } << 513 // 'Bare' statistics
610 else << 514 fSumEnergyKilled += endEnergy;
611 { << 515 if( endEnergy > fMaxEnergyKilled) { fMaxEnergyKilled= endEnergy; }
612 fMaxEnergySaved = std::max( endEnergy, << 516
613 if( fNoLooperTrials == 1 ) { <<
614 fSumEnergySaved += endEnergy; <<
615 if ( !stable ) <<
616 fSumEnergyUnstableSaved += endEne <<
617 } <<
618 #ifdef G4VERBOSE 517 #ifdef G4VERBOSE
619 if( verboseLevel > 2 && ! fSilenceLoop << 518 if( (fVerboseLevel > 1) ||
620 { << 519 ( endEnergy > fThreshold_Warning_Energy ) ) {
621 G4cout << " " << __func__ << 520 G4cout << " G4Transportation is killing track that is looping or stuck "
622 << " Particle is looping but << 521 << G4endl
623 << " Number of trials = " < << 522 << " This track has " << track.GetKineticEnergy() / MeV
624 << " No of calls to = " << << 523 << " MeV energy." << G4endl;
625 } << 524 }
626 #endif 525 #endif
>> 526 fNoLooperTrials=0;
627 } 527 }
628 } << 528 else{
629 else << 529 fNoLooperTrials ++;
630 { << 530 }
>> 531 }else{
631 fNoLooperTrials=0; 532 fNoLooperTrials=0;
632 } 533 }
633 534
634 // Another (sometimes better way) is to use 535 // Another (sometimes better way) is to use a user-limit maximum Step size
635 // to alleviate this problem .. 536 // to alleviate this problem ..
636 537
637 // Introduce smooth curved trajectories to p 538 // Introduce smooth curved trajectories to particle-change
638 // 539 //
639 fParticleChange.SetPointerToVectorOfAuxiliar 540 fParticleChange.SetPointerToVectorOfAuxiliaryPoints
640 (fFieldPropagator->GimmeTrajectoryVectorAn 541 (fFieldPropagator->GimmeTrajectoryVectorAndForgetIt() );
641 542
642 return &fParticleChange ; 543 return &fParticleChange ;
643 } 544 }
644 545
645 ////////////////////////////////////////////// 546 //////////////////////////////////////////////////////////////////////////
646 // 547 //
647 // This ensures that the PostStep action is a 548 // This ensures that the PostStep action is always called,
648 // so that it can do the relocation if it is 549 // so that it can do the relocation if it is needed.
649 // 550 //
650 551
651 G4double G4Transportation:: 552 G4double G4Transportation::
652 PostStepGetPhysicalInteractionLength( const G4 553 PostStepGetPhysicalInteractionLength( const G4Track&,
653 G4 554 G4double, // previousStepSize
654 G4 555 G4ForceCondition* pForceCond )
655 { << 556 {
656 fFieldExertedForce = false; // Not known <<
657 *pForceCond = Forced ; 557 *pForceCond = Forced ;
658 return DBL_MAX ; // was kInfinity ; but con 558 return DBL_MAX ; // was kInfinity ; but convention now is DBL_MAX
659 } 559 }
660 560
661 ////////////////////////////////////////////// 561 /////////////////////////////////////////////////////////////////////////////
662 <<
663 void G4Transportation::SetTouchableInformation <<
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 } <<
704 <<
705 ////////////////////////////////////////////// <<
706 // 562 //
707 563
708 G4VParticleChange* G4Transportation::PostStepD 564 G4VParticleChange* G4Transportation::PostStepDoIt( const G4Track& track,
709 565 const G4Step& )
710 { 566 {
711 G4TouchableHandle retCurrentTouchable ; / << 567 G4TouchableHandle retCurrentTouchable ; // The one to return
712 G4bool isLastStep= false; <<
713 568
714 // Initialize ParticleChange (by setting al 569 // Initialize ParticleChange (by setting all its members equal
715 // to correspond 570 // to corresponding members in G4Track)
716 // fParticleChange.Initialize(track) ; // T 571 // fParticleChange.Initialize(track) ; // To initialise TouchableChange
717 572
718 fParticleChange.ProposeTrackStatus(track.Get 573 fParticleChange.ProposeTrackStatus(track.GetTrackStatus()) ;
719 574
720 // If the Step was determined by the volume 575 // If the Step was determined by the volume boundary,
721 // logically relocate the particle 576 // logically relocate the particle
722 << 577
723 if(fGeometryLimitedStep) 578 if(fGeometryLimitedStep)
724 { 579 {
725 // fCurrentTouchable will now become the p 580 // fCurrentTouchable will now become the previous touchable,
726 // and what was the previous will be freed 581 // and what was the previous will be freed.
727 // (Needed because the preStepPoint can po 582 // (Needed because the preStepPoint can point to the previous touchable)
728 583
729 fLinearNavigator->SetGeometricallyLimitedS 584 fLinearNavigator->SetGeometricallyLimitedStep() ;
730 fLinearNavigator-> 585 fLinearNavigator->
731 LocateGlobalPointAndUpdateTouchableHandle( 586 LocateGlobalPointAndUpdateTouchableHandle( track.GetPosition(),
732 587 track.GetMomentumDirection(),
733 588 fCurrentTouchableHandle,
734 589 true ) ;
735 // Check whether the particle is out of th 590 // Check whether the particle is out of the world volume
736 // If so it has exited and must be killed. 591 // If so it has exited and must be killed.
737 // 592 //
738 if( fCurrentTouchableHandle->GetVolume() = 593 if( fCurrentTouchableHandle->GetVolume() == 0 )
739 { 594 {
740 fParticleChange.ProposeTrackStatus( fSt 595 fParticleChange.ProposeTrackStatus( fStopAndKill ) ;
741 } 596 }
742 retCurrentTouchable = fCurrentTouchableHan 597 retCurrentTouchable = fCurrentTouchableHandle ;
743 fParticleChange.SetTouchableHandle( fCurre 598 fParticleChange.SetTouchableHandle( fCurrentTouchableHandle ) ;
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 } 599 }
752 else // fGeometryLimitedStep 600 else // fGeometryLimitedStep is false
753 { 601 {
>> 602 #ifdef G4DEBUG_POSTSTEP_TRANSPORT
>> 603
>> 604 // Although the location is changed, we know that the physical
>> 605 // volume remains constant.
>> 606 // In order to help in checking the user geometry
>> 607 // we perform a full-relocation and check its result
>> 608 // *except* if we have made a very small step from a boundary
>> 609 // (i.e. remaining inside the tolerance)
>> 610
>> 611 G4bool startAtSurface_And_MoveEpsilon ;
>> 612 startAtSurface_And_MoveEpsilon =
>> 613 (stepData.GetPreStepPoint()->GetSafety() == 0.0) &&
>> 614 (stepData.GetStepLength() < kCarTolerance ) ;
>> 615
>> 616 if( startAtSurface_And_MoveEpsilon )
>> 617 {
>> 618 fLinearNavigator->
>> 619 LocateGlobalPointAndUpdateTouchableHandle( track.GetPosition(),
>> 620 track.GetMomentumDirection(),
>> 621 fCurrentTouchableHandle,
>> 622 true );
>> 623 if( fCurrentTouchableHandle->GetVolume() != track.GetVolume() )
>> 624 {
>> 625 G4cerr << " ERROR: A relocation within safety has"
>> 626 << " caused a volume change! " << G4endl ;
>> 627 G4cerr << " The old volume is called "
>> 628 << track.GetVolume()->GetName() << G4endl ;
>> 629 G4cerr << " The new volume is called " ;
>> 630
>> 631 if ( fCurrentTouchableHandle->GetVolume() != 0 )
>> 632 {
>> 633 G4cerr << fCurrentTouchableHandle->GetVolume()->GetName()
>> 634 << G4endl ;
>> 635 }
>> 636 else
>> 637 {
>> 638 G4cerr << "Out of World" << G4endl ;
>> 639 }
>> 640 G4cerr.precision(7) ;
>> 641 G4cerr << " The position is " << track.GetPosition() << G4endl ;
>> 642
>> 643 // Let us relocate again, for debuging
>> 644 //
>> 645 fLinearNavigator->
>> 646 LocateGlobalPointAndUpdateTouchableHandle(track.GetPosition(),
>> 647 track.GetMomentumDirection(),
>> 648 fCurrentTouchableHandle,
>> 649 true ) ;
>> 650 G4cerr << " The newer volume is called " ;
>> 651
>> 652 if ( fCurrentTouchableHandle->GetVolume() != 0 )
>> 653 {
>> 654 G4cerr << fCurrentTouchableHandle->GetVolume()->GetName()
>> 655 << G4endl ;
>> 656 }
>> 657 else
>> 658 {
>> 659 G4cerr << "Out of World" << G4endl ;
>> 660 }
>> 661 }
>> 662
>> 663 assert( fCurrentTouchableHandle->GetVolume()->GetName() ==
>> 664 track.GetVolume()->GetName() ) ;
>> 665
>> 666 retCurrentTouchable = fCurrentTouchableHandle ;
>> 667 fParticleChange.SetTouchableHandle( fCurrentTouchableHandle ) ;
>> 668
>> 669 }
>> 670 else
>> 671 {
>> 672 retCurrentTouchable = track.GetTouchableHandle() ;
>> 673 fParticleChange.SetTouchableHandle( track.GetTouchableHandle() ) ;
>> 674 }
>> 675
>> 676 // This must be done in the above if ( AtSur ) fails
>> 677 // We also do it for if (true) in order to get debug/opt to
>> 678 // behave as exactly the same way as possible.
>> 679 //
>> 680 fLinearNavigator->LocateGlobalPointWithinVolume( track.GetPosition() ) ;
>> 681
>> 682 #else // ie #ifndef G4DEBUG_POSTSTEP_TRANSPORT does a quick relocation
>> 683
754 // This serves only to move the Navigator' 684 // This serves only to move the Navigator's location
755 // 685 //
756 fLinearNavigator->LocateGlobalPointWithinV 686 fLinearNavigator->LocateGlobalPointWithinVolume( track.GetPosition() ) ;
757 687
758 // The value of the track's current Toucha 688 // The value of the track's current Touchable is retained.
759 // (and it must be correct because we must 689 // (and it must be correct because we must use it below to
760 // overwrite the (unset) one in particle c 690 // overwrite the (unset) one in particle change)
761 // It must be fCurrentTouchable too ?? << 691 // Although in general this is fCurrentTouchable, at the start of
>> 692 // a step it could be different ... ??
762 // 693 //
763 fParticleChange.SetTouchableHandle( track. 694 fParticleChange.SetTouchableHandle( track.GetTouchableHandle() ) ;
764 retCurrentTouchable = track.GetTouchableHa 695 retCurrentTouchable = track.GetTouchableHandle() ;
765 696
766 isLastStep= false; << 697 #endif
767 } // endif ( fGeometryLimitedStep ) 698 } // endif ( fGeometryLimitedStep )
768 fLastStepInVolume= isLastStep; <<
769 <<
770 fParticleChange.ProposeFirstStepInVolume(fFi <<
771 fParticleChange.ProposeLastStepInVolume(isLa <<
772 699
773 SetTouchableInformation(retCurrentTouchable) << 700 const G4VPhysicalVolume* pNewVol = retCurrentTouchable->GetVolume() ;
>> 701 const G4Material* pNewMaterial = 0 ;
>> 702 const G4VSensitiveDetector* pNewSensitiveDetector = 0 ;
>> 703
>> 704 if( pNewVol != 0 )
>> 705 {
>> 706 pNewMaterial= pNewVol->GetLogicalVolume()->GetMaterial();
>> 707 pNewSensitiveDetector= pNewVol->GetLogicalVolume()->GetSensitiveDetector();
>> 708 }
>> 709
>> 710 // ( <const_cast> pNewMaterial ) ;
>> 711 // ( <const_cast> pNewSensitiveDetector) ;
>> 712
>> 713 fParticleChange.SetMaterialInTouchable( (G4Material *) pNewMaterial ) ;
>> 714 fParticleChange.SetSensitiveDetectorInTouchable( (G4VSensitiveDetector *) pNewSensitiveDetector ) ;
>> 715
>> 716 const G4MaterialCutsCouple* pNewMaterialCutsCouple = 0;
>> 717 if( pNewVol != 0 )
>> 718 {
>> 719 pNewMaterialCutsCouple=pNewVol->GetLogicalVolume()->GetMaterialCutsCouple();
>> 720 }
>> 721
>> 722 if( pNewVol!=0 && pNewMaterialCutsCouple!=0 && pNewMaterialCutsCouple->GetMaterial()!=pNewMaterial )
>> 723 {
>> 724 // for parametrized volume
>> 725 //
>> 726 pNewMaterialCutsCouple =
>> 727 G4ProductionCutsTable::GetProductionCutsTable()
>> 728 ->GetMaterialCutsCouple(pNewMaterial,
>> 729 pNewMaterialCutsCouple->GetProductionCuts());
>> 730 }
>> 731 fParticleChange.SetMaterialCutsCoupleInTouchable( pNewMaterialCutsCouple );
>> 732
>> 733 // temporarily until Get/Set Material of ParticleChange,
>> 734 // and StepPoint can be made const.
>> 735 // Set the touchable in ParticleChange
>> 736 // this must always be done because the particle change always
>> 737 // uses this value to overwrite the current touchable pointer.
>> 738 //
>> 739 fParticleChange.SetTouchableHandle(retCurrentTouchable) ;
774 740
775 return &fParticleChange ; 741 return &fParticleChange ;
776 } 742 }
777 743
778 ////////////////////////////////////////////// << 744 // New method takes over the responsibility to reset the state of G4Transportation
779 // New method takes over the responsibility to << 745 // object at the start of a new track or the resumption of a suspended track.
780 // G4Transportation object at the start of a n <<
781 // of a suspended track. <<
782 // <<
783 746
784 void 747 void
785 G4Transportation::StartTracking(G4Track* aTrac 748 G4Transportation::StartTracking(G4Track* aTrack)
786 { 749 {
787 G4VProcess::StartTracking(aTrack); 750 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 751
795 // Whether field exists should be determined << 752 // The actions here are those that were taken in AlongStepGPIL
796 G4FieldManagerStore* fieldMgrStore= G4FieldM << 753 // when track.GetCurrentStepNumber()==1
797 fAnyFieldExists = fieldMgrStore->size() > 0; << 754
798 <<
799 // reset safety value and center 755 // reset safety value and center
800 // 756 //
801 fPreviousSafety = 0.0 ; 757 fPreviousSafety = 0.0 ;
802 fPreviousSftOrigin = G4ThreeVector(0.,0.,0.) 758 fPreviousSftOrigin = G4ThreeVector(0.,0.,0.) ;
803 759
804 // reset looping counter -- for motion in fi 760 // reset looping counter -- for motion in field
805 fNoLooperTrials= 0; << 761 if( aTrack->GetCurrentStepNumber()==1 ) {
806 // Must clear this state .. else it depends << 762 fNoLooperTrials= 0;
807 // --> a better solution would set this fro << 763 }
808 // Was if( aTrack->GetCurrentStepNumber()==1 <<
809 764
810 // ChordFinder reset internal state 765 // ChordFinder reset internal state
811 // 766 //
812 if( fFieldPropagator && fAnyFieldExists ) << 767 if( DoesGlobalFieldExist() ) {
813 { << 768 G4ChordFinder* chordF= fFieldPropagator->GetChordFinder();
814 fFieldPropagator->ClearPropagatorState(); << 769 if( chordF ) chordF->ResetStepEstimate();
815 // Resets all state of field propagator <<
816 // values (in case of overlaps and to w <<
817 } 770 }
818 << 771
819 // Make sure to clear the chord finders of a <<
820 // <<
821 fieldMgrStore->ClearAllChordFindersState(); <<
822 <<
823 // Update the current touchable handle (fro 772 // Update the current touchable handle (from the track's)
824 // 773 //
825 fCurrentTouchableHandle = aTrack->GetTouchab 774 fCurrentTouchableHandle = aTrack->GetTouchableHandle();
826 <<
827 // Inform field propagator of new track <<
828 // <<
829 fFieldPropagator->PrepareNewTrack(); <<
830 } <<
831 <<
832 ////////////////////////////////////////////// <<
833 // <<
834 <<
835 G4bool G4Transportation::EnableMagneticMoment( <<
836 { <<
837 G4bool lastValue= fUseMagneticMoment; <<
838 fUseMagneticMoment= useMoment; <<
839 return lastValue; <<
840 } <<
841 <<
842 ////////////////////////////////////////////// <<
843 // <<
844 <<
845 G4bool G4Transportation::EnableGravity(G4bool <<
846 { <<
847 G4bool lastValue= fUseGravity; <<
848 fUseGravity= useGravity; <<
849 return lastValue; <<
850 } <<
851 <<
852 ////////////////////////////////////////////// <<
853 // <<
854 // Supress (or not) warnings about 'looping' <<
855 <<
856 void G4Transportation::SetSilenceLooperWarning <<
857 { <<
858 fSilenceLooperWarnings= val; // Flag to *Su <<
859 } 775 }
860 776
861 ////////////////////////////////////////////// <<
862 // <<
863 G4bool G4Transportation::GetSilenceLooperWarni <<
864 { <<
865 return fSilenceLooperWarnings; <<
866 } <<
867 <<
868 <<
869 ////////////////////////////////////////////// <<
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 <<
879 G4int maxTrials = 10; <<
880 SetThresholdTrials( maxTrials ); <<
881 <<
882 PushThresholdsToLogger(); // Again, to be s <<
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 <<
910 <<
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 } <<
938 777