| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // 26 //
27 // 27 //
28 // ------------------------------------------- 28 // ------------------------------------------------------------
29 // GEANT 4 class implementation 29 // GEANT 4 class implementation
30 // 30 //
31 // =========================================== 31 // =======================================================================
32 // Created: 19 March 1997, J. Apostolakis 32 // Created: 19 March 1997, J. Apostolakis
33 // =========================================== 33 // =======================================================================
34 34
35 #include "G4CoupledTransportation.hh" 35 #include "G4CoupledTransportation.hh"
36 #include "G4TransportationProcessType.hh" 36 #include "G4TransportationProcessType.hh"
37 #include "G4TransportationLogger.hh" 37 #include "G4TransportationLogger.hh"
38 38
39 #include "G4PhysicalConstants.hh" 39 #include "G4PhysicalConstants.hh"
40 #include "G4SystemOfUnits.hh" 40 #include "G4SystemOfUnits.hh"
41 #include "G4ProductionCutsTable.hh" 41 #include "G4ProductionCutsTable.hh"
42 #include "G4ParticleTable.hh" 42 #include "G4ParticleTable.hh"
43 #include "G4ChordFinder.hh" 43 #include "G4ChordFinder.hh"
44 #include "G4Field.hh" 44 #include "G4Field.hh"
45 #include "G4FieldTrack.hh" 45 #include "G4FieldTrack.hh"
46 #include "G4FieldManagerStore.hh" 46 #include "G4FieldManagerStore.hh"
47 #include "G4PathFinder.hh" 47 #include "G4PathFinder.hh"
48 48
49 #include "G4PropagatorInField.hh" 49 #include "G4PropagatorInField.hh"
50 #include "G4TransportationManager.hh" 50 #include "G4TransportationManager.hh"
51 51
52 class G4VSensitiveDetector; 52 class G4VSensitiveDetector;
53 53
54 G4bool G4CoupledTransportation::fSignifyStepIn 54 G4bool G4CoupledTransportation::fSignifyStepInAnyVolume= false;
55 // This mode must apply to all threads 55 // This mode must apply to all threads
56 56
57 ////////////////////////////////////////////// 57 //////////////////////////////////////////////////////////////////////////
58 // 58 //
59 // Constructor 59 // Constructor
60 60
61 G4CoupledTransportation::G4CoupledTransportati 61 G4CoupledTransportation::G4CoupledTransportation( G4int verbosity )
62 : G4Transportation( verbosity, "CoupledTrans 62 : G4Transportation( verbosity, "CoupledTransportation" ),
63 fPreviousMassSafety( 0.0 ), 63 fPreviousMassSafety( 0.0 ),
64 fPreviousFullSafety( 0.0 ), 64 fPreviousFullSafety( 0.0 ),
65 fMassGeometryLimitedStep( false ), 65 fMassGeometryLimitedStep( false ),
66 fFirstStepInMassVolume( true ) 66 fFirstStepInMassVolume( true )
67 { 67 {
68 SetProcessSubType(static_cast<G4int>(COUPLED 68 SetProcessSubType(static_cast<G4int>(COUPLED_TRANSPORTATION));
69 // SetVerboseLevel is called in the construc 69 // SetVerboseLevel is called in the constructor of G4Transportation
70 70
71 if( verboseLevel > 0 ) 71 if( verboseLevel > 0 )
72 { 72 {
73 G4cout << " G4CoupledTransportation constr 73 G4cout << " G4CoupledTransportation constructor: ----- " << G4endl;
74 G4cout << " Verbose level is " << verboseL 74 G4cout << " Verbose level is " << verboseLevel << G4endl;
75 G4cout << " Reports First/Last in " 75 G4cout << " Reports First/Last in "
76 << (fSignifyStepInAnyVolume ? " any 76 << (fSignifyStepInAnyVolume ? " any " : " mass " )
77 << " geometry " << G4endl; 77 << " geometry " << G4endl;
78 } 78 }
79 fPathFinder= G4PathFinder::GetInstance(); 79 fPathFinder= G4PathFinder::GetInstance();
80 } 80 }
81 81
82 ////////////////////////////////////////////// 82 //////////////////////////////////////////////////////////////////////////
83 83
84 G4CoupledTransportation::~G4CoupledTransportat 84 G4CoupledTransportation::~G4CoupledTransportation()
85 { 85 {
86 } 86 }
87 87
88 ////////////////////////////////////////////// 88 //////////////////////////////////////////////////////////////////////////
89 // 89 //
90 // Responsibilities: 90 // Responsibilities:
91 // Find whether the geometry limits the Ste 91 // Find whether the geometry limits the Step, and to what length
92 // Calculate the new value of the safety an 92 // Calculate the new value of the safety and return it.
93 // Store the final time, position and momen 93 // Store the final time, position and momentum.
94 94
95 G4double G4CoupledTransportation:: 95 G4double G4CoupledTransportation::
96 AlongStepGetPhysicalInteractionLength( const G 96 AlongStepGetPhysicalInteractionLength( const G4Track& track,
97 G 97 G4double, // previousStepSize
98 G 98 G4double currentMinimumStep,
99 G 99 G4double& proposedSafetyForStart,
100 G 100 G4GPILSelection* selection )
101 { 101 {
102 G4double geometryStepLength; 102 G4double geometryStepLength;
103 G4double startFullSafety= 0.0; // estimated 103 G4double startFullSafety= 0.0; // estimated safety for start point (all geometries)
104 G4double safetyProposal= -1.0; // local copy 104 G4double safetyProposal= -1.0; // local copy of proposal
105 105
106 G4ThreeVector EndUnitMomentum ; 106 G4ThreeVector EndUnitMomentum ;
107 G4double lengthAlongCurve = 0.0 ; 107 G4double lengthAlongCurve = 0.0 ;
108 108
109 fParticleIsLooping = false ; 109 fParticleIsLooping = false ;
110 110
111 // Initial actions moved to StartTrack() 111 // Initial actions moved to StartTrack()
112 // -------------------------------------- 112 // --------------------------------------
113 // Note: in case another process changes tou 113 // Note: in case another process changes touchable handle
114 // it will be necessary to add here (for 114 // it will be necessary to add here (for all steps)
115 // fCurrentTouchableHandle = aTrack->GetTouc 115 // fCurrentTouchableHandle = aTrack->GetTouchableHandle();
116 116
117 // GPILSelection is set to defaule value of 117 // GPILSelection is set to defaule value of CandidateForSelection
118 // It is a return value 118 // It is a return value
119 // 119 //
120 *selection = CandidateForSelection ; 120 *selection = CandidateForSelection ;
121 121
122 fFirstStepInMassVolume = fNewTrack || fMassG 122 fFirstStepInMassVolume = fNewTrack || fMassGeometryLimitedStep ;
123 fFirstStepInVolume = fNewTrack || fGeome 123 fFirstStepInVolume = fNewTrack || fGeometryLimitedStep ;
124 124
125 #ifdef G4DEBUG_TRANSPORT 125 #ifdef G4DEBUG_TRANSPORT
126 G4cout << " CoupledTransport::AlongStep GPI 126 G4cout << " CoupledTransport::AlongStep GPIL: "
127 << " 1st-step: any= " <<fFirstStep 127 << " 1st-step: any= " <<fFirstStepInVolume << " ( geom= "
128 << fGeometryLimitedStep << " ) " 128 << fGeometryLimitedStep << " ) "
129 << " mass= " << fFirstStepI 129 << " mass= " << fFirstStepInMassVolume << " ( geom= "
130 << fMassGeometryLimitedStep << " ) " 130 << fMassGeometryLimitedStep << " ) "
131 << " newTrack= " << fNewTrack << G4e 131 << " newTrack= " << fNewTrack << G4endl;
132 #endif 132 #endif
133 133
134 // fLastStepInVolume= false; 134 // fLastStepInVolume= false;
135 fNewTrack = false; 135 fNewTrack = false;
136 136
137 // Get initial Energy/Momentum of the track 137 // Get initial Energy/Momentum of the track
138 // 138 //
139 const G4DynamicParticle* pParticle = tra 139 const G4DynamicParticle* pParticle = track.GetDynamicParticle() ;
140 const G4ParticleDefinition* pParticleDef = 140 const G4ParticleDefinition* pParticleDef = pParticle->GetDefinition() ;
141 G4ThreeVector startMomentumDir = pPart 141 G4ThreeVector startMomentumDir = pParticle->GetMomentumDirection() ;
142 G4ThreeVector startPosition = track 142 G4ThreeVector startPosition = track.GetPosition() ;
143 G4VPhysicalVolume* currentVolume= track.GetV 143 G4VPhysicalVolume* currentVolume= track.GetVolume();
144 144
145 #ifdef G4DEBUG_TRANSPORT 145 #ifdef G4DEBUG_TRANSPORT
146 if( verboseLevel > 1 ) 146 if( verboseLevel > 1 )
147 { 147 {
148 G4cout << "G4CoupledTransportation::AlongS 148 G4cout << "G4CoupledTransportation::AlongStepGPIL> called in volume "
149 << currentVolume->GetName() << G4en 149 << currentVolume->GetName() << G4endl;
150 } 150 }
151 #endif 151 #endif
152 // G4double theTime = track.GetGlob 152 // G4double theTime = track.GetGlobalTime() ;
153 153
154 // The Step Point safety can be limited by o 154 // The Step Point safety can be limited by other geometries and/or the
155 // assumptions of any process - it's not alw 155 // assumptions of any process - it's not always the geometrical safety.
156 // We calculate the starting point's isotrop 156 // We calculate the starting point's isotropic safety here.
157 // 157 //
158 G4ThreeVector OriginShift = startPosition - 158 G4ThreeVector OriginShift = startPosition - fPreviousSftOrigin ;
159 G4double MagSqShift = OriginShift.mag2 159 G4double MagSqShift = OriginShift.mag2() ;
160 startFullSafety= 0.0; 160 startFullSafety= 0.0;
161 161
162 // Recall that FullSafety <= MassSafety 162 // Recall that FullSafety <= MassSafety
163 // Original: if( MagSqShift < sqr(fPreviousM 163 // Original: if( MagSqShift < sqr(fPreviousMassSafety) ) {
164 if( MagSqShift < sqr(fPreviousFullSafety) ) 164 if( MagSqShift < sqr(fPreviousFullSafety) )
165 { 165 {
166 G4double mag_shift= std::sqrt(MagSqShift) 166 G4double mag_shift= std::sqrt(MagSqShift);
167 startFullSafety = std::max( (fPreviousFul 167 startFullSafety = std::max( (fPreviousFullSafety - mag_shift), 0.0);
168 // Need to be consistent between full s 168 // Need to be consistent between full safety with Mass safety
169 // in order reproduce results in simple 169 // in order reproduce results in simple case
170 // --> use same calculation method 170 // --> use same calculation method
171 171
172 // Only compute full safety if massSafety 172 // Only compute full safety if massSafety > 0. Else it remains 0
173 // startFullSafety = fPathFinder->Compute 173 // startFullSafety = fPathFinder->ComputeSafety( startPosition );
174 } 174 }
175 175
176 // Is the particle charged or has it a magne 176 // Is the particle charged or has it a magnetic moment?
177 // 177 //
178 G4double particleCharge = pParticle->GetChar 178 G4double particleCharge = pParticle->GetCharge() ;
179 G4double magneticMoment = pParticle->GetMagn 179 G4double magneticMoment = pParticle->GetMagneticMoment() ;
180 G4double restMass = pParticle->GetMass 180 G4double restMass = pParticle->GetMass() ;
181 181
182 fMassGeometryLimitedStep = false ; // Set d 182 fMassGeometryLimitedStep = false ; // Set default - alex
183 fGeometryLimitedStep = false; 183 fGeometryLimitedStep = false;
184 184
185 // There is no need to locate the current vo 185 // There is no need to locate the current volume. It is Done elsewhere:
186 // On track construction 186 // On track construction
187 // By the tracking, after all AlongStepDoI 187 // By the tracking, after all AlongStepDoIts, in "Relocation"
188 188
189 // Check if the particle has a force, EM or 189 // Check if the particle has a force, EM or gravitational, exerted on it
190 // 190 //
191 G4FieldManager* fieldMgr= nullptr; 191 G4FieldManager* fieldMgr= nullptr;
192 G4bool fieldExertsForce = false ; 192 G4bool fieldExertsForce = false ;
193 193
194 const G4Field* ptrField= nullptr; 194 const G4Field* ptrField= nullptr;
195 195
196 fieldMgr = fFieldPropagator->FindAndSetField 196 fieldMgr = fFieldPropagator->FindAndSetFieldManager( track.GetVolume() );
197 G4bool eligibleEM = (particleCharge != 0.0) 197 G4bool eligibleEM = (particleCharge != 0.0)
198 || ( fUseMagneticMoment && 198 || ( fUseMagneticMoment && (magneticMoment != 0.0) );
199 G4bool eligibleGrav = fUseGravity && (restM 199 G4bool eligibleGrav = fUseGravity && (restMass != 0.0) ;
200 200
201 if( (fieldMgr!=nullptr) && (eligibleEM||elig 201 if( (fieldMgr!=nullptr) && (eligibleEM||eligibleGrav) )
202 { 202 {
203 // Message the field Manager, to configur 203 // Message the field Manager, to configure it for this track
204 // 204 //
205 fieldMgr->ConfigureForTrack( &track ); 205 fieldMgr->ConfigureForTrack( &track );
206 206
207 // The above call can transition from a n 207 // The above call can transition from a null field-ptr oto a finite field.
208 // If the field manager has no field ptr, 208 // If the field manager has no field ptr, the field is zero
209 // by definition ( = there is no field ! 209 // by definition ( = there is no field ! )
210 // 210 //
211 ptrField= fieldMgr->GetDetectorField(); 211 ptrField= fieldMgr->GetDetectorField();
212 212
213 if( ptrField != nullptr) 213 if( ptrField != nullptr)
214 { 214 {
215 fieldExertsForce = eligibleEM 215 fieldExertsForce = eligibleEM
216 || ( eligibleGrav && ptrField->I 216 || ( eligibleGrav && ptrField->IsGravityActive() );
217 } 217 }
218 } 218 }
219 G4double momentumMagnitude = pParticle->GetT 219 G4double momentumMagnitude = pParticle->GetTotalMomentum() ;
220 220
221 if( fieldExertsForce ) 221 if( fieldExertsForce )
222 { 222 {
223 auto equationOfMotion= fFieldPropagator-> 223 auto equationOfMotion= fFieldPropagator->GetCurrentEquationOfMotion();
224 224
225 G4ChargeState chargeState(particleCharge, 225 G4ChargeState chargeState(particleCharge, // Charge can change (dynamic)
226 magneticMoment, 226 magneticMoment,
227 pParticleDef->G 227 pParticleDef->GetPDGSpin() );
228 if( equationOfMotion ) 228 if( equationOfMotion )
229 { 229 {
230 equationOfMotion->SetChargeMomentumMas 230 equationOfMotion->SetChargeMomentumMass( chargeState,
231 231 momentumMagnitude,
232 232 restMass );
233 } 233 }
234 #ifdef G4DEBUG_TRANSPORT 234 #ifdef G4DEBUG_TRANSPORT
235 else 235 else
236 { 236 {
237 G4cerr << " ERROR in G4CoupledTranspor 237 G4cerr << " ERROR in G4CoupledTransportation> "
238 << "Cannot find valid Equation 238 << "Cannot find valid Equation of motion: " << G4endl;
239 << " Unable to pass Charge, Mom 239 << " Unable to pass Charge, Momentum and Mass " << G4endl;
240 } 240 }
241 #endif 241 #endif
242 } 242 }
243 243
244 G4ThreeVector polarizationVec = track.GetPo 244 G4ThreeVector polarizationVec = track.GetPolarization() ;
245 G4FieldTrack aFieldTrack = G4FieldTrack(sta 245 G4FieldTrack aFieldTrack = G4FieldTrack(startPosition,
246 tra 246 track.GetGlobalTime(), // Lab.
247 tra 247 track.GetMomentumDirection(),
248 tra 248 track.GetKineticEnergy(),
249 res 249 restMass,
250 par 250 particleCharge,
251 pol 251 polarizationVec,
252 pPa 252 pParticleDef->GetPDGMagneticMoment(),
253 0.0 253 0.0, // Length along track
254 pPa 254 pParticleDef->GetPDGSpin() ) ;
255 G4int stepNo= track.GetCurrentStepNumber(); 255 G4int stepNo= track.GetCurrentStepNumber();
256 256
257 ELimited limitedStep; 257 ELimited limitedStep;
258 G4FieldTrack endTrackState('a'); // Defaul 258 G4FieldTrack endTrackState('a'); // Default values
259 259
260 fMassGeometryLimitedStep = false ; // de 260 fMassGeometryLimitedStep = false ; // default
261 fGeometryLimitedStep = false; 261 fGeometryLimitedStep = false;
262 if( currentMinimumStep > 0 ) 262 if( currentMinimumStep > 0 )
263 { 263 {
264 G4double newMassSafety= 0.0; // tem 264 G4double newMassSafety= 0.0; // temp. for recalculation
265 265
266 // Do the Transport in the field (non re 266 // Do the Transport in the field (non recti-linear)
267 // 267 //
268 lengthAlongCurve = fPathFinder->ComputeS 268 lengthAlongCurve = fPathFinder->ComputeStep( aFieldTrack,
269 269 currentMinimumStep,
270 270 G4TransportationManager::kMassNavigatorId,
271 271 stepNo,
272 272 newMassSafety,
273 273 limitedStep,
274 274 endTrackState,
275 275 currentVolume ) ;
276 276
277 G4double newFullSafety= fPathFinder->Get 277 G4double newFullSafety= fPathFinder->GetCurrentSafety();
278 // this was estimated already in step 278 // this was estimated already in step above
279 279
280 if( limitedStep == kUnique || limitedSte 280 if( limitedStep == kUnique || limitedStep == kSharedTransport )
281 { 281 {
282 fMassGeometryLimitedStep = true ; 282 fMassGeometryLimitedStep = true ;
283 } 283 }
284 284
285 fGeometryLimitedStep = (fPathFinder->Get 285 fGeometryLimitedStep = (fPathFinder->GetNumberGeometriesLimitingStep() != 0);
286 286
287 #ifdef G4DEBUG_TRANSPORT 287 #ifdef G4DEBUG_TRANSPORT
288 if( fMassGeometryLimitedStep && !fGeomet 288 if( fMassGeometryLimitedStep && !fGeometryLimitedStep )
289 { 289 {
290 std::ostringstream message; 290 std::ostringstream message;
291 message << " ERROR in determining geom 291 message << " ERROR in determining geometries limiting the step" << G4endl;
292 message << " Limiting: mass=" << fMa 292 message << " Limiting: mass=" << fMassGeometryLimitedStep
293 << " any= " << fGeometryLimite 293 << " any= " << fGeometryLimitedStep << G4endl;
294 message << "Incompatible conditions - 294 message << "Incompatible conditions - by which geometries was it limited ?";
295 G4Exception("G4CoupledTransportation:: 295 G4Exception("G4CoupledTransportation::AlongStepGetPhysicalInteractionLength()",
296 "PathFinderConfused", Fata 296 "PathFinderConfused", FatalException, message);
297 } 297 }
298 #endif 298 #endif
299 299
300 geometryStepLength = std::min( lengthAlo 300 geometryStepLength = std::min( lengthAlongCurve, currentMinimumStep);
301 301
302 // Momentum: Magnitude and direction ca 302 // Momentum: Magnitude and direction can be changed too now ...
303 // 303 //
304 fMomentumChanged = true ; 304 fMomentumChanged = true ;
305 fTransportEndMomentumDir = endTrackState 305 fTransportEndMomentumDir = endTrackState.GetMomentumDir() ;
306 306
307 // Remember last safety origin & value. 307 // Remember last safety origin & value.
308 fPreviousSftOrigin = startPosition ; 308 fPreviousSftOrigin = startPosition ;
309 fPreviousMassSafety = newMassSafety ; 309 fPreviousMassSafety = newMassSafety ;
310 fPreviousFullSafety = newFullSafety ; 310 fPreviousFullSafety = newFullSafety ;
311 // fpSafetyHelper->SetCurrentSafety( new 311 // fpSafetyHelper->SetCurrentSafety( newFullSafety, startPosition);
312 312
313 #ifdef G4DEBUG_TRANSPORT 313 #ifdef G4DEBUG_TRANSPORT
314 if( verboseLevel > 1 ) 314 if( verboseLevel > 1 )
315 { 315 {
316 G4cout << "G4Transport:CompStep> " 316 G4cout << "G4Transport:CompStep> "
317 << " called the pathfinder for 317 << " called the pathfinder for a new step at " << startPosition
318 << " and obtained step = " << l 318 << " and obtained step = " << lengthAlongCurve << G4endl;
319 G4cout << " New safety (preStep) = " 319 G4cout << " New safety (preStep) = " << newMassSafety << G4endl;
320 } 320 }
321 #endif 321 #endif
322 322
323 // Store as best estimate value 323 // Store as best estimate value
324 startFullSafety = newFullSafety ; 324 startFullSafety = newFullSafety ;
325 325
326 // Get the End-Position and End-Momentum 326 // Get the End-Position and End-Momentum (Dir-ection)
327 fTransportEndPosition = endTrackState.Ge 327 fTransportEndPosition = endTrackState.GetPosition() ;
328 fTransportEndKineticEnergy = endTrackSt 328 fTransportEndKineticEnergy = endTrackState.GetKineticEnergy() ;
329 } 329 }
330 else 330 else
331 { 331 {
332 geometryStepLength = lengthAlongCurve= 332 geometryStepLength = lengthAlongCurve= 0.0 ;
333 fMomentumChanged = false ; 333 fMomentumChanged = false ;
334 // fMassGeometryLimitedStep = false ; 334 // fMassGeometryLimitedStep = false ; // --- ???
335 // fGeometryLimitedStep = true; 335 // fGeometryLimitedStep = true;
336 fTransportEndMomentumDir = track.GetMome 336 fTransportEndMomentumDir = track.GetMomentumDirection();
337 fTransportEndKineticEnergy = track.GetK 337 fTransportEndKineticEnergy = track.GetKineticEnergy();
338 338
339 fTransportEndPosition = startPosition; 339 fTransportEndPosition = startPosition;
340 340
341 endTrackState= aFieldTrack; // Ensures 341 endTrackState= aFieldTrack; // Ensures that time is updated
342 } 342 }
343 // G4FieldTrack aTrackState(endTrackState); 343 // G4FieldTrack aTrackState(endTrackState);
344 344
345 if( !fieldExertsForce ) 345 if( !fieldExertsForce )
346 { 346 {
347 fParticleIsLooping = false ; 347 fParticleIsLooping = false ;
348 fMomentumChanged = false ; 348 fMomentumChanged = false ;
349 fEndGlobalTimeComputed = false ; 349 fEndGlobalTimeComputed = false ;
350 } 350 }
351 else 351 else
352 { 352 {
353 fParticleIsLooping = fFieldPropagator->I 353 fParticleIsLooping = fFieldPropagator->IsParticleLooping() ;
354 354
355 #ifdef G4DEBUG_TRANSPORT 355 #ifdef G4DEBUG_TRANSPORT
356 if( verboseLevel > 1 ) 356 if( verboseLevel > 1 )
357 { 357 {
358 G4cout << " G4CT::CS End Position = " 358 G4cout << " G4CT::CS End Position = "
359 << fTransportEndPosition << G4e 359 << fTransportEndPosition << G4endl;
360 G4cout << " G4CT::CS End Direction = " 360 G4cout << " G4CT::CS End Direction = "
361 << fTransportEndMomentumDir << 361 << fTransportEndMomentumDir << G4endl;
362 } 362 }
363 #endif 363 #endif
364 if( fFieldPropagator->GetCurrentFieldMan 364 if( fFieldPropagator->GetCurrentFieldManager()->DoesFieldChangeEnergy() )
365 { 365 {
366 // If the field can change energy, t 366 // If the field can change energy, then the time must be integrated
367 // - so this should have been upd 367 // - so this should have been updated
368 // 368 //
369 fCandidateEndGlobalTime = endTrack 369 fCandidateEndGlobalTime = endTrackState.GetLabTimeOfFlight();
370 fEndGlobalTimeComputed = true; 370 fEndGlobalTimeComputed = true;
371 371
372 // was ( fCandidateEndGlobalTime != 372 // was ( fCandidateEndGlobalTime != track.GetGlobalTime() );
373 // a cleaner way is to have FieldTra 373 // a cleaner way is to have FieldTrack knowing whether time
374 // is updated 374 // is updated
375 } 375 }
376 else 376 else
377 { 377 {
378 // The energy should be unchanged by 378 // The energy should be unchanged by field transport,
379 // - so the time changed will be 379 // - so the time changed will be calculated elsewhere
380 // 380 //
381 fEndGlobalTimeComputed = false; 381 fEndGlobalTimeComputed = false;
382 382
383 #ifdef G4VERBOSE 383 #ifdef G4VERBOSE
384 // Check that the integration preser 384 // Check that the integration preserved the energy
385 // - and if not correct this! 385 // - and if not correct this!
386 G4double startEnergy= track.GetKine 386 G4double startEnergy= track.GetKineticEnergy();
387 G4double endEnergy= fTransportEndKi 387 G4double endEnergy= fTransportEndKineticEnergy;
388 388
389 G4double absEdiff = std::fabs(startE 389 G4double absEdiff = std::fabs(startEnergy- endEnergy);
390 if( (verboseLevel > 1) && ( absEdiff 390 if( (verboseLevel > 1) && ( absEdiff > perThousand * endEnergy) )
391 { 391 {
392 ReportInexactEnergy(startEnergy, e 392 ReportInexactEnergy(startEnergy, endEnergy);
393 } // end of if (verboseLevel) 393 } // end of if (verboseLevel)
394 #endif 394 #endif
395 // Correct the energy for fields tha 395 // Correct the energy for fields that conserve it
396 // This - hides the integration err 396 // This - hides the integration error
397 // - but gives a better physic 397 // - but gives a better physical answer
398 fTransportEndKineticEnergy= track.Ge 398 fTransportEndKineticEnergy= track.GetKineticEnergy();
399 } 399 }
400 } 400 }
401 401
402 fEndPointDistance = (fTransportEndPosition 402 fEndPointDistance = (fTransportEndPosition - startPosition).mag() ;
403 fTransportEndSpin = endTrackState.GetSpin(); 403 fTransportEndSpin = endTrackState.GetSpin();
404 404
405 // Calculate the safety 405 // Calculate the safety
406 406
407 safetyProposal= startFullSafety; // used t 407 safetyProposal= startFullSafety; // used to be startMassSafety
408 // Changed to accomodate processes that ca 408 // Changed to accomodate processes that cannot update the safety
409 409
410 // Update safety for the end-point, if becom 410 // Update safety for the end-point, if becomes negative at the end-point.
411 411
412 if( (startFullSafety < fEndPointDistance ) 412 if( (startFullSafety < fEndPointDistance )
413 && ( particleCharge != 0.0 ) ) // Onl 413 && ( particleCharge != 0.0 ) ) // Only needed to prepare for MSC
414 // && !fGeometryLimitedStep ) // To-Try: 414 // && !fGeometryLimitedStep ) // To-Try: No safety update if at boundary
415 { 415 {
416 G4double endFullSafety = 416 G4double endFullSafety =
417 fPathFinder->ComputeSafety( fTransport 417 fPathFinder->ComputeSafety( fTransportEndPosition);
418 // Expected mission -- only mass geome 418 // Expected mission -- only mass geometry's safety
419 // fLinearNavigator->ComputeSafety( 419 // fLinearNavigator->ComputeSafety( fTransportEndPosition) ;
420 // Yet discrete processes only have po 420 // Yet discrete processes only have poststep -- and this cannot
421 // currently revise the safety 421 // currently revise the safety
422 // ==> so we use the all-geometry sa 422 // ==> so we use the all-geometry safety as a precaution
423 423
424 fpSafetyHelper->SetCurrentSafety( endFul 424 fpSafetyHelper->SetCurrentSafety( endFullSafety, fTransportEndPosition);
425 // Pushing safety to Helper avoids rec 425 // Pushing safety to Helper avoids recalculation at this point
426 426
427 G4ThreeVector centerPt= G4ThreeVector(0. 427 G4ThreeVector centerPt= G4ThreeVector(0.0, 0.0, 0.0); // Used for return value
428 G4double endMassSafety= fPathFinder->Obt 428 G4double endMassSafety= fPathFinder->ObtainSafety( G4TransportationManager::kMassNavigatorId, centerPt);
429 // Retrieves the mass value from Path 429 // Retrieves the mass value from PathFinder (it calculated it)
430 430
431 fPreviousMassSafety = endMassSafety ; 431 fPreviousMassSafety = endMassSafety ;
432 fPreviousFullSafety = endFullSafety; 432 fPreviousFullSafety = endFullSafety;
433 fPreviousSftOrigin = fTransportEndPositi 433 fPreviousSftOrigin = fTransportEndPosition ;
434 434
435 // The convention (Stepping Manager's) i 435 // The convention (Stepping Manager's) is safety from the start point
436 // 436 //
437 safetyProposal = endFullSafety + fEndPoi 437 safetyProposal = endFullSafety + fEndPointDistance;
438 // --> was endMassSafety 438 // --> was endMassSafety
439 // Changed to accomodate processes that 439 // Changed to accomodate processes that cannot update the safety
440 440
441 #ifdef G4DEBUG_TRANSPORT 441 #ifdef G4DEBUG_TRANSPORT
442 G4int prec= G4cout.precision(12) ; 442 G4int prec= G4cout.precision(12) ;
443 G4cout << "***CoupledTransportation::Alo 443 G4cout << "***CoupledTransportation::AlongStepGPIL ** " << G4endl ;
444 G4cout << " Revised Safety at endpoint 444 G4cout << " Revised Safety at endpoint " << fTransportEndPosition
445 << " give safety values: Mass= 445 << " give safety values: Mass= " << endMassSafety
446 << " All= " << endFullSafety << 446 << " All= " << endFullSafety << G4endl ;
447 G4cout << " Adding endpoint distance " 447 G4cout << " Adding endpoint distance " << fEndPointDistance
448 << " to obtain pseudo-safety= " 448 << " to obtain pseudo-safety= " << safetyProposal << G4endl ;
449 G4cout.precision(prec); 449 G4cout.precision(prec);
450 } 450 }
451 else 451 else
452 { 452 {
453 G4int prec= G4cout.precision(12) ; 453 G4int prec= G4cout.precision(12) ;
454 G4cout << "***CoupledTransportation::Alo 454 G4cout << "***CoupledTransportation::AlongStepGPIL ** " << G4endl ;
455 G4cout << " Quick Safety estimate at en 455 G4cout << " Quick Safety estimate at endpoint "
456 << fTransportEndPosition 456 << fTransportEndPosition
457 << " gives safety endpoint valu 457 << " gives safety endpoint value = "
458 << startFullSafety - fEndPointDis 458 << startFullSafety - fEndPointDistance
459 << " using start-point value " < 459 << " using start-point value " << startFullSafety
460 << " and endpointDistance " << f 460 << " and endpointDistance " << fEndPointDistance << G4endl;
461 G4cout.precision(prec); 461 G4cout.precision(prec);
462 #endif 462 #endif
463 } 463 }
464 464
465 proposedSafetyForStart= safetyProposal; 465 proposedSafetyForStart= safetyProposal;
466 fParticleChange.ProposeTrueStepLength(geomet 466 fParticleChange.ProposeTrueStepLength(geometryStepLength) ;
467 467
468 return geometryStepLength ; 468 return geometryStepLength ;
469 } 469 }
470 470
471 ////////////////////////////////////////////// 471 /////////////////////////////////////////////////////////////////////////////
472 472
473 void G4CoupledTransportation:: 473 void G4CoupledTransportation::
474 ReportMove( G4ThreeVector OldVector, G4ThreeVe 474 ReportMove( G4ThreeVector OldVector, G4ThreeVector NewVector,
475 const G4String& Quantity ) 475 const G4String& Quantity )
476 { 476 {
477 G4ThreeVector moveVec = ( NewVector - OldV 477 G4ThreeVector moveVec = ( NewVector - OldVector );
478 478
479 G4cerr << G4endl 479 G4cerr << G4endl
480 << "******************************* 480 << "**************************************************************"
481 << G4endl; 481 << G4endl;
482 G4cerr << "Endpoint has moved between valu 482 G4cerr << "Endpoint has moved between value expected from TransportEndPosition "
483 << " and value from Track in PostSt 483 << " and value from Track in PostStepDoIt. " << G4endl
484 << "Change of " << Quantity << " is 484 << "Change of " << Quantity << " is " << moveVec.mag() / mm
485 << " mm long, " 485 << " mm long, "
486 << " and its vector is " << (1.0/mm 486 << " and its vector is " << (1.0/mm) * moveVec << " mm " << G4endl
487 << "Endpoint of ComputeStep was " < 487 << "Endpoint of ComputeStep was " << OldVector
488 << " and current position to locate 488 << " and current position to locate is " << NewVector << G4endl;
489 } 489 }
490 490
491 ////////////////////////////////////////////// 491 /////////////////////////////////////////////////////////////////////////////
492 492
493 G4VParticleChange* G4CoupledTransportation::Po 493 G4VParticleChange* G4CoupledTransportation::PostStepDoIt( const G4Track& track,
494 494 const G4Step& )
495 { 495 {
496 G4TouchableHandle retCurrentTouchable ; // 496 G4TouchableHandle retCurrentTouchable ; // The one to return
497 497
498 // Initialize ParticleChange (by setting al 498 // Initialize ParticleChange (by setting all its members equal
499 // to correspond 499 // to corresponding members in G4Track)
500 // fParticleChange.Initialize(track) ; // T 500 // fParticleChange.Initialize(track) ; // To initialise TouchableChange
501 501
502 fParticleChange.ProposeTrackStatus(track.Get 502 fParticleChange.ProposeTrackStatus(track.GetTrackStatus()) ;
503 503
504 if( fSignifyStepInAnyVolume ) 504 if( fSignifyStepInAnyVolume )
505 { 505 {
506 fParticleChange.ProposeFirstStepInVolume( 506 fParticleChange.ProposeFirstStepInVolume( fFirstStepInVolume );
507 } 507 }
508 else 508 else
509 { 509 {
510 fParticleChange.ProposeFirstStepInVolume( 510 fParticleChange.ProposeFirstStepInVolume( fFirstStepInMassVolume );
511 } 511 }
512 512
513 // Check that the end position and direction 513 // Check that the end position and direction are preserved
514 // since call to AlongStepDoIt 514 // since call to AlongStepDoIt
515 515
516 #ifdef G4DEBUG_TRANSPORT 516 #ifdef G4DEBUG_TRANSPORT
517 if( ( verboseLevel > 0 ) 517 if( ( verboseLevel > 0 )
518 && ((fTransportEndPosition - track.GetPos 518 && ((fTransportEndPosition - track.GetPosition()).mag2() >= 1.0e-16) )
519 { 519 {
520 ReportMove( track.GetPosition(), fTranspo 520 ReportMove( track.GetPosition(), fTransportEndPosition,
521 "End of Step Position" ); 521 "End of Step Position" );
522 G4cerr << " Problem in G4CoupledTransport 522 G4cerr << " Problem in G4CoupledTransportation::PostStepDoIt " << G4endl;
523 } 523 }
524 524
525 // If the Step was determined by the volume 525 // If the Step was determined by the volume boundary, relocate the particle
526 // The pathFinder will know that the geometr 526 // The pathFinder will know that the geometry limited the step (!?)
527 527
528 if( verboseLevel > 0 ) 528 if( verboseLevel > 0 )
529 { 529 {
530 G4cout << " Calling PathFinder::Locate() 530 G4cout << " Calling PathFinder::Locate() from "
531 << " G4CoupledTransportation::Post 531 << " G4CoupledTransportation::PostStepDoIt " << G4endl;
532 G4cout << " fGeometryLimitedStep is " < 532 G4cout << " fGeometryLimitedStep is " << fGeometryLimitedStep << G4endl;
533 } 533 }
534 #endif 534 #endif
535 535
536 if(fGeometryLimitedStep) 536 if(fGeometryLimitedStep)
537 { 537 {
538 fPathFinder->Locate( track.GetPosition(), 538 fPathFinder->Locate( track.GetPosition(),
539 track.GetMomentumDire 539 track.GetMomentumDirection(),
540 true); 540 true);
541 541
542 // fCurrentTouchable will now become the p 542 // fCurrentTouchable will now become the previous touchable,
543 // and what was the previous will be freed 543 // and what was the previous will be freed.
544 // (Needed because the preStepPoint can po 544 // (Needed because the preStepPoint can point to the previous touchable)
545 545
546 fCurrentTouchableHandle= 546 fCurrentTouchableHandle=
547 fPathFinder->CreateTouchableHandle( G4Tr 547 fPathFinder->CreateTouchableHandle( G4TransportationManager::kMassNavigatorId );
548 548
549 #ifdef G4DEBUG_TRANSPORT 549 #ifdef G4DEBUG_TRANSPORT
550 if( verboseLevel > 1 ) 550 if( verboseLevel > 1 )
551 { 551 {
552 G4VPhysicalVolume* vol= fCurrentTouchab 552 G4VPhysicalVolume* vol= fCurrentTouchableHandle->GetVolume();
553 G4cout << "CHECK !!!!!!!!!!! fCurrentTo 553 G4cout << "CHECK !!!!!!!!!!! fCurrentTouchableHandle->GetVolume() = "
554 << vol; 554 << vol;
555 if( vol ) { G4cout << "Name=" << vol->G 555 if( vol ) { G4cout << "Name=" << vol->GetName(); }
556 G4cout << G4endl; 556 G4cout << G4endl;
557 } 557 }
558 #endif 558 #endif
559 559
560 // Check whether the particle is out of th 560 // Check whether the particle is out of the world volume
561 // If so it has exited and must be killed. 561 // If so it has exited and must be killed.
562 // 562 //
563 if( fCurrentTouchableHandle->GetVolume() = 563 if( fCurrentTouchableHandle->GetVolume() == 0 )
564 { 564 {
565 fParticleChange.ProposeTrackStatus( fSt 565 fParticleChange.ProposeTrackStatus( fStopAndKill ) ;
566 } 566 }
567 retCurrentTouchable = fCurrentTouchableHan 567 retCurrentTouchable = fCurrentTouchableHandle ;
568 // fParticleChange.SetTouchableHandle( fCu 568 // fParticleChange.SetTouchableHandle( fCurrentTouchableHandle ) ;
569 } 569 }
570 else // fGeometryLimitedStep is false 570 else // fGeometryLimitedStep is false
571 { 571 {
572 #ifdef G4DEBUG_TRANSPORT 572 #ifdef G4DEBUG_TRANSPORT
573 if( verboseLevel > 1 ) 573 if( verboseLevel > 1 )
574 { 574 {
575 G4cout << "G4CoupledTransportation::Post 575 G4cout << "G4CoupledTransportation::PostStepDoIt -- "
576 << " fGeometryLimitedStep = " << 576 << " fGeometryLimitedStep = " << fGeometryLimitedStep
577 << " must be false " << G4endl; 577 << " must be false " << G4endl;
578 } 578 }
579 #endif 579 #endif
580 // This serves only to move each of the Na 580 // This serves only to move each of the Navigator's location
581 // 581 //
582 // fLinearNavigator->LocateGlobalPointWith 582 // fLinearNavigator->LocateGlobalPointWithinVolume( track.GetPosition() ) ;
583 583
584 fPathFinder->ReLocate( track.GetPosition() 584 fPathFinder->ReLocate( track.GetPosition() );
585 // track.GetMomentu 585 // track.GetMomentumDirection() );
586 586
587 // Keep the value of the track's current T 587 // Keep the value of the track's current Touchable is retained,
588 // and use it to overwrite the (unset) on 588 // and use it to overwrite the (unset) one in particle change.
589 // Expect this must be fCurrentTouchable t 589 // Expect this must be fCurrentTouchable too
590 // - could it be different, eg at the st 590 // - could it be different, eg at the start of a step ?
591 // 591 //
592 retCurrentTouchable = track.GetTouchableHa 592 retCurrentTouchable = track.GetTouchableHandle() ;
593 // fParticleChange.SetTouchableHandle( tra 593 // fParticleChange.SetTouchableHandle( track.GetTouchableHandle() ) ;
594 } // endif ( fGeometryLimitedStep ) 594 } // endif ( fGeometryLimitedStep )
595 595
596 #ifdef G4DEBUG_NAVIGATION 596 #ifdef G4DEBUG_NAVIGATION
597 G4cout << " CoupledTransport::AlongStep GPI 597 G4cout << " CoupledTransport::AlongStep GPIL: "
598 << " last-step: any= " << fGeometryL 598 << " last-step: any= " << fGeometryLimitedStep << " . ..... x . "
599 << " mass= " << fMassGeometryLimitedS 599 << " mass= " << fMassGeometryLimitedStep << G4endl;
600 #endif 600 #endif
601 601
602 if( fSignifyStepInAnyVolume ) 602 if( fSignifyStepInAnyVolume )
603 fParticleChange.ProposeLastStepInVolume(fG 603 fParticleChange.ProposeLastStepInVolume(fGeometryLimitedStep);
604 else 604 else
605 fParticleChange.ProposeLastStepInVolume(f 605 fParticleChange.ProposeLastStepInVolume(fMassGeometryLimitedStep);
606 606
607 SetTouchableInformation(retCurrentTouchable) 607 SetTouchableInformation(retCurrentTouchable);
608 608
609 return &fParticleChange ; 609 return &fParticleChange ;
610 } 610 }
611 611
612 ////////////////////////////////////////////// 612 /////////////////////////////////////////////////////////////////////////////
613 // New method takes over the responsibility to 613 // New method takes over the responsibility to reset the state of
614 // G4CoupledTransportation object: 614 // G4CoupledTransportation object:
615 // - at the start of a new track, and 615 // - at the start of a new track, and
616 // - on the resumption of a suspended tra 616 // - on the resumption of a suspended track.
617 // 617 //
618 void 618 void
619 G4CoupledTransportation::StartTracking(G4Track 619 G4CoupledTransportation::StartTracking(G4Track* aTrack)
620 { 620 {
621 G4Transportation::StartTracking(aTrack); 621 G4Transportation::StartTracking(aTrack);
622 622
623 G4ThreeVector position = aTrack->GetPosition 623 G4ThreeVector position = aTrack->GetPosition();
624 G4ThreeVector direction = aTrack->GetMomentu 624 G4ThreeVector direction = aTrack->GetMomentumDirection();
625 625
626 fPathFinder->PrepareNewTrack( position, dire 626 fPathFinder->PrepareNewTrack( position, direction);
627 // This implies a call to fPathFinder->Locat 627 // This implies a call to fPathFinder->Locate( position, direction );
628 628
629 // reset safety value and center 629 // reset safety value and center
630 // 630 //
631 fPreviousMassSafety = 0.0 ; 631 fPreviousMassSafety = 0.0 ;
632 fPreviousFullSafety = 0.0 ; 632 fPreviousFullSafety = 0.0 ;
633 fPreviousSftOrigin = G4ThreeVector(0.,0.,0.) 633 fPreviousSftOrigin = G4ThreeVector(0.,0.,0.) ;
634 } 634 }
635 635
636 ////////////////////////////////////////////// 636 /////////////////////////////////////////////////////////////////////////////
637 637
638 void 638 void
639 G4CoupledTransportation::EndTracking() 639 G4CoupledTransportation::EndTracking()
640 { 640 {
641 G4TransportationManager::GetTransportationMa 641 G4TransportationManager::GetTransportationManager()->InactivateAll();
642 fPathFinder->EndTrack(); 642 fPathFinder->EndTrack();
643 // Resets TransportationManager to use ord 643 // Resets TransportationManager to use ordinary Navigator
644 } 644 }
645 645
646 ////////////////////////////////////////////// 646 /////////////////////////////////////////////////////////////////////////////
647 647
648 void 648 void
649 G4CoupledTransportation:: 649 G4CoupledTransportation::
650 ReportInexactEnergy(G4double startEnergy, G4do 650 ReportInexactEnergy(G4double startEnergy, G4double endEnergy)
651 { 651 {
652 static G4ThreadLocal G4int no_warnings= 0, w 652 static G4ThreadLocal G4int no_warnings= 0, warnModulo=1,
653 moduloFactor= 10, 653 moduloFactor= 10, no_large_ediff= 0;
654 654
655 if( std::fabs(startEnergy- endEnergy) > perT 655 if( std::fabs(startEnergy- endEnergy) > perThousand * endEnergy )
656 { 656 {
657 no_large_ediff ++; 657 no_large_ediff ++;
658 if( (no_large_ediff% warnModulo) == 0 ) 658 if( (no_large_ediff% warnModulo) == 0 )
659 { 659 {
660 no_warnings++; 660 no_warnings++;
661 std::ostringstream message; 661 std::ostringstream message;
662 message << "Energy change in Step is abo 662 message << "Energy change in Step is above 1^-3 relative value. "
663 << G4endl 663 << G4endl
664 << " Relative change in 'track 664 << " Relative change in 'tracking' step = "
665 << std::setw(15) << (endEnergy-s 665 << std::setw(15) << (endEnergy-startEnergy)/startEnergy
666 << G4endl 666 << G4endl
667 << " Starting E= " << std::set 667 << " Starting E= " << std::setw(12) << startEnergy / MeV
668 << " MeV " << G4endl 668 << " MeV " << G4endl
669 << " Ending E= " << std::set 669 << " Ending E= " << std::setw(12) << endEnergy / MeV
670 << " MeV " << G4endl 670 << " MeV " << G4endl
671 << "Energy has been corrected -- 671 << "Energy has been corrected -- however, review"
672 << " field propagation parameter 672 << " field propagation parameters for accuracy." << G4endl;
673 if ( (verboseLevel > 2 ) || (no_warnings 673 if ( (verboseLevel > 2 ) || (no_warnings<4)
674 || (no_large_ediff == warnModulo * mod 674 || (no_large_ediff == warnModulo * moduloFactor) )
675 { 675 {
676 message << "These include EpsilonStepM 676 message << "These include EpsilonStepMax(/Min) in G4FieldManager,"
677 << G4endl 677 << G4endl
678 << "which determine fractional 678 << "which determine fractional error per step for integrated quantities."
679 << G4endl 679 << G4endl
680 << "Note also the influence of 680 << "Note also the influence of the permitted number of integration steps."
681 << G4endl; 681 << G4endl;
682 } 682 }
683 message << "Bad 'endpoint'. Energy chang 683 message << "Bad 'endpoint'. Energy change detected and corrected."
684 << G4endl 684 << G4endl
685 << "Has occurred already " << no 685 << "Has occurred already " << no_large_ediff << " times.";
686 G4Exception("G4CoupledTransportation::Al 686 G4Exception("G4CoupledTransportation::AlongStepGetPIL()",
687 "EnergyChange", JustWarning, 687 "EnergyChange", JustWarning, message);
688 if( no_large_ediff == warnModulo * modul 688 if( no_large_ediff == warnModulo * moduloFactor )
689 { 689 {
690 warnModulo *= moduloFactor; 690 warnModulo *= moduloFactor;
691 } 691 }
692 } 692 }
693 } 693 }
694 } 694 }
695 695