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