| 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" <<
48 47
49 #include "G4PropagatorInField.hh" << 48 #include "G4Transportation.hh" // In order to use fUseMagneticMoment
50 #include "G4TransportationManager.hh" <<
51 49
52 class G4VSensitiveDetector; 50 class G4VSensitiveDetector;
53 51
54 G4bool G4CoupledTransportation::fSignifyStepIn 52 G4bool G4CoupledTransportation::fSignifyStepInAnyVolume= false;
55 // This mode must apply to all threads 53 // This mode must apply to all threads
56 54
>> 55 G4bool G4CoupledTransportation::fUseMagneticMoment=false;
>> 56 G4bool G4CoupledTransportation::fUseGravity= false;
>> 57 G4bool G4CoupledTransportation::fSilenceLooperWarnings= false;
57 ////////////////////////////////////////////// 58 //////////////////////////////////////////////////////////////////////////
58 // 59 //
59 // Constructor 60 // Constructor
60 61
61 G4CoupledTransportation::G4CoupledTransportati 62 G4CoupledTransportation::G4CoupledTransportation( G4int verbosity )
62 : G4Transportation( verbosity, "CoupledTrans << 63 : G4VProcess( G4String("CoupledTransportation"), fTransportation ),
>> 64 fTransportEndPosition(0.0, 0.0, 0.0),
>> 65 fTransportEndMomentumDir(0.0, 0.0, 0.0),
>> 66 fTransportEndKineticEnergy(0.0),
>> 67 fTransportEndSpin(0.0, 0.0, 0.0),
>> 68 fMomentumChanged(false),
>> 69 fEndGlobalTimeComputed(false),
>> 70 fCandidateEndGlobalTime(0.0),
>> 71 fParticleIsLooping( false ),
>> 72 fNewTrack( true ),
>> 73 fPreviousSftOrigin( 0.,0.,0. ),
63 fPreviousMassSafety( 0.0 ), 74 fPreviousMassSafety( 0.0 ),
64 fPreviousFullSafety( 0.0 ), 75 fPreviousFullSafety( 0.0 ),
65 fMassGeometryLimitedStep( false ), 76 fMassGeometryLimitedStep( false ),
66 fFirstStepInMassVolume( true ) << 77 fAnyGeometryLimitedStep( false ),
>> 78 fEndpointDistance( -1.0 ),
>> 79 fSumEnergyKilled( 0.0 ), fMaxEnergyKilled( 0.0 ),
>> 80 fFirstStepInMassVolume( true ),
>> 81 fFirstStepInAnyVolume( true )
67 { 82 {
68 SetProcessSubType(static_cast<G4int>(COUPLED 83 SetProcessSubType(static_cast<G4int>(COUPLED_TRANSPORTATION));
69 // SetVerboseLevel is called in the construc << 84 SetVerboseLevel(verbosity);
70 85
>> 86 G4TransportationManager* transportMgr ;
>> 87
>> 88 transportMgr = G4TransportationManager::GetTransportationManager() ;
>> 89
>> 90 fMassNavigator = transportMgr->GetNavigatorForTracking() ;
>> 91 fFieldPropagator = transportMgr->GetPropagatorInField() ;
>> 92 // fGlobalFieldMgr = transportMgr->GetFieldManager() ;
>> 93 fNavigatorId= transportMgr->ActivateNavigator( fMassNavigator );
71 if( verboseLevel > 0 ) 94 if( verboseLevel > 0 )
72 { 95 {
73 G4cout << " G4CoupledTransportation constr 96 G4cout << " G4CoupledTransportation constructor: ----- " << G4endl;
74 G4cout << " Verbose level is " << verboseL 97 G4cout << " Verbose level is " << verboseLevel << G4endl;
>> 98 G4cout << " Navigator Id obtained in G4CoupledTransportation constructor "
>> 99 << fNavigatorId << G4endl;
75 G4cout << " Reports First/Last in " 100 G4cout << " Reports First/Last in "
76 << (fSignifyStepInAnyVolume ? " any 101 << (fSignifyStepInAnyVolume ? " any " : " mass " )
77 << " geometry " << G4endl; 102 << " geometry " << G4endl;
78 } 103 }
79 fPathFinder= G4PathFinder::GetInstance(); 104 fPathFinder= G4PathFinder::GetInstance();
>> 105 fpSafetyHelper = transportMgr->GetSafetyHelper(); // New
>> 106
>> 107 fpLogger = new G4TransportationLogger("G4Transportation", verbosity);
>> 108
>> 109 SetHighLooperThresholds();
>> 110 // Use the old defaults: Warning = 100 MeV, Important = 250 MeV, No Trials = 10;
>> 111
>> 112 PushThresholdsToLogger();
>> 113 // Should be done by Set methods in SetHighLooperThresholds -- making sure
>> 114
>> 115 static G4ThreadLocal G4TouchableHandle* pNullTouchableHandle = 0;
>> 116 if ( !pNullTouchableHandle) { pNullTouchableHandle = new G4TouchableHandle; }
>> 117 fCurrentTouchableHandle = *pNullTouchableHandle;
>> 118 // Points to (G4VTouchable*) 0
>> 119
>> 120 fAnyFieldExists= DoesAnyFieldExist();
80 } 121 }
81 122
82 ////////////////////////////////////////////// 123 //////////////////////////////////////////////////////////////////////////
83 124
84 G4CoupledTransportation::~G4CoupledTransportat 125 G4CoupledTransportation::~G4CoupledTransportation()
85 { 126 {
>> 127 if( fSumEnergyKilled > 0.0 )
>> 128 {
>> 129 PrintStatistics( G4cout );
>> 130 }
>> 131 delete fpLogger;
>> 132 }
>> 133
>> 134 //////////////////////////////////////////////////////////////////////////
>> 135
>> 136 void
>> 137 G4CoupledTransportation::PrintStatistics( std::ostream& outStr) const
>> 138 {
>> 139 if( fSumEnergyKilled > 0.0 )
>> 140 {
>> 141 outStr << " G4CoupledTransportation: Statistics for looping particles "
>> 142 << G4endl;
>> 143 outStr << " Sum of energy of loopers killed: "
>> 144 << fSumEnergyKilled / MeV << " MeV " << G4endl;
>> 145 outStr << " Max energy of loopers killed: "
>> 146 << fMaxEnergyKilled / MeV << " MeV " << G4endl;
>> 147
>> 148
>> 149 outStr << " Max energy of loopers 'saved': " << fMaxEnergySaved << G4endl;
>> 150 outStr << " Sum of energy of loopers 'saved': "
>> 151 << fSumEnergySaved << G4endl;
>> 152 outStr << " Sum of energy of unstable loopers 'saved': "
>> 153 << fSumEnergyUnstableSaved << G4endl;
>> 154 }
86 } 155 }
87 156
88 ////////////////////////////////////////////// 157 //////////////////////////////////////////////////////////////////////////
89 // 158 //
90 // Responsibilities: 159 // Responsibilities:
91 // Find whether the geometry limits the Ste 160 // Find whether the geometry limits the Step, and to what length
92 // Calculate the new value of the safety an 161 // Calculate the new value of the safety and return it.
93 // Store the final time, position and momen 162 // Store the final time, position and momentum.
94 163
95 G4double G4CoupledTransportation:: 164 G4double G4CoupledTransportation::
96 AlongStepGetPhysicalInteractionLength( const G 165 AlongStepGetPhysicalInteractionLength( const G4Track& track,
97 G 166 G4double, // previousStepSize
98 G 167 G4double currentMinimumStep,
99 G 168 G4double& proposedSafetyForStart,
100 G 169 G4GPILSelection* selection )
101 { 170 {
102 G4double geometryStepLength; 171 G4double geometryStepLength;
>> 172 G4double startMassSafety= 0.0; // estimated safety for start point (mass geometry)
103 G4double startFullSafety= 0.0; // estimated 173 G4double startFullSafety= 0.0; // estimated safety for start point (all geometries)
104 G4double safetyProposal= -1.0; // local copy 174 G4double safetyProposal= -1.0; // local copy of proposal
105 175
106 G4ThreeVector EndUnitMomentum ; 176 G4ThreeVector EndUnitMomentum ;
107 G4double lengthAlongCurve = 0.0 ; 177 G4double lengthAlongCurve = 0.0 ;
108 178
109 fParticleIsLooping = false ; 179 fParticleIsLooping = false ;
110 180
111 // Initial actions moved to StartTrack() 181 // Initial actions moved to StartTrack()
112 // -------------------------------------- 182 // --------------------------------------
113 // Note: in case another process changes tou 183 // Note: in case another process changes touchable handle
114 // it will be necessary to add here (for 184 // it will be necessary to add here (for all steps)
115 // fCurrentTouchableHandle = aTrack->GetTouc 185 // fCurrentTouchableHandle = aTrack->GetTouchableHandle();
116 186
117 // GPILSelection is set to defaule value of 187 // GPILSelection is set to defaule value of CandidateForSelection
118 // It is a return value 188 // It is a return value
119 // 189 //
120 *selection = CandidateForSelection ; 190 *selection = CandidateForSelection ;
121 191
122 fFirstStepInMassVolume = fNewTrack || fMassG 192 fFirstStepInMassVolume = fNewTrack || fMassGeometryLimitedStep ;
123 fFirstStepInVolume = fNewTrack || fGeome << 193 fFirstStepInAnyVolume = fNewTrack || fAnyGeometryLimitedStep ;
124 194
125 #ifdef G4DEBUG_TRANSPORT 195 #ifdef G4DEBUG_TRANSPORT
126 G4cout << " CoupledTransport::AlongStep GPI 196 G4cout << " CoupledTransport::AlongStep GPIL: "
127 << " 1st-step: any= " <<fFirstStep << 197 << " 1st-step: any= " <<fFirstStepInAnyVolume << " ( geom= "
128 << fGeometryLimitedStep << " ) " << 198 << fAnyGeometryLimitedStep << " ) "
129 << " mass= " << fFirstStepI 199 << " mass= " << fFirstStepInMassVolume << " ( geom= "
130 << fMassGeometryLimitedStep << " ) " 200 << fMassGeometryLimitedStep << " ) "
131 << " newTrack= " << fNewTrack << G4e 201 << " newTrack= " << fNewTrack << G4endl;
132 #endif 202 #endif
133 203
134 // fLastStepInVolume= false; 204 // fLastStepInVolume= false;
135 fNewTrack = false; 205 fNewTrack = false;
136 206
137 // Get initial Energy/Momentum of the track 207 // Get initial Energy/Momentum of the track
138 // 208 //
139 const G4DynamicParticle* pParticle = tra 209 const G4DynamicParticle* pParticle = track.GetDynamicParticle() ;
140 const G4ParticleDefinition* pParticleDef = 210 const G4ParticleDefinition* pParticleDef = pParticle->GetDefinition() ;
141 G4ThreeVector startMomentumDir = pPart 211 G4ThreeVector startMomentumDir = pParticle->GetMomentumDirection() ;
142 G4ThreeVector startPosition = track 212 G4ThreeVector startPosition = track.GetPosition() ;
143 G4VPhysicalVolume* currentVolume= track.GetV 213 G4VPhysicalVolume* currentVolume= track.GetVolume();
144 214
145 #ifdef G4DEBUG_TRANSPORT 215 #ifdef G4DEBUG_TRANSPORT
146 if( verboseLevel > 1 ) 216 if( verboseLevel > 1 )
147 { 217 {
148 G4cout << "G4CoupledTransportation::AlongS 218 G4cout << "G4CoupledTransportation::AlongStepGPIL> called in volume "
149 << currentVolume->GetName() << G4en 219 << currentVolume->GetName() << G4endl;
150 } 220 }
151 #endif 221 #endif
152 // G4double theTime = track.GetGlob 222 // G4double theTime = track.GetGlobalTime() ;
153 223
154 // The Step Point safety can be limited by o 224 // The Step Point safety can be limited by other geometries and/or the
155 // assumptions of any process - it's not alw 225 // assumptions of any process - it's not always the geometrical safety.
156 // We calculate the starting point's isotrop 226 // We calculate the starting point's isotropic safety here.
157 // 227 //
158 G4ThreeVector OriginShift = startPosition - 228 G4ThreeVector OriginShift = startPosition - fPreviousSftOrigin ;
159 G4double MagSqShift = OriginShift.mag2 229 G4double MagSqShift = OriginShift.mag2() ;
>> 230 startMassSafety = 0.0;
160 startFullSafety= 0.0; 231 startFullSafety= 0.0;
161 232
162 // Recall that FullSafety <= MassSafety 233 // Recall that FullSafety <= MassSafety
163 // Original: if( MagSqShift < sqr(fPreviousM 234 // Original: if( MagSqShift < sqr(fPreviousMassSafety) ) {
164 if( MagSqShift < sqr(fPreviousFullSafety) ) 235 if( MagSqShift < sqr(fPreviousFullSafety) )
165 { 236 {
166 G4double mag_shift= std::sqrt(MagSqShift) 237 G4double mag_shift= std::sqrt(MagSqShift);
>> 238 startMassSafety = std::max( (fPreviousMassSafety - mag_shift), 0.0);
167 startFullSafety = std::max( (fPreviousFul 239 startFullSafety = std::max( (fPreviousFullSafety - mag_shift), 0.0);
168 // Need to be consistent between full s 240 // Need to be consistent between full safety with Mass safety
169 // in order reproduce results in simple 241 // in order reproduce results in simple case
170 // --> use same calculation method 242 // --> use same calculation method
171 243
172 // Only compute full safety if massSafety 244 // Only compute full safety if massSafety > 0. Else it remains 0
173 // startFullSafety = fPathFinder->Compute 245 // startFullSafety = fPathFinder->ComputeSafety( startPosition );
174 } 246 }
175 247
176 // Is the particle charged or has it a magne 248 // Is the particle charged or has it a magnetic moment?
177 // 249 //
178 G4double particleCharge = pParticle->GetChar 250 G4double particleCharge = pParticle->GetCharge() ;
179 G4double magneticMoment = pParticle->GetMagn 251 G4double magneticMoment = pParticle->GetMagneticMoment() ;
180 G4double restMass = pParticle->GetMass 252 G4double restMass = pParticle->GetMass() ;
181 253
182 fMassGeometryLimitedStep = false ; // Set d 254 fMassGeometryLimitedStep = false ; // Set default - alex
183 fGeometryLimitedStep = false; << 255 fAnyGeometryLimitedStep = false;
184 256
185 // There is no need to locate the current vo 257 // There is no need to locate the current volume. It is Done elsewhere:
186 // On track construction 258 // On track construction
187 // By the tracking, after all AlongStepDoI 259 // By the tracking, after all AlongStepDoIts, in "Relocation"
188 260
189 // Check if the particle has a force, EM or 261 // Check if the particle has a force, EM or gravitational, exerted on it
190 // 262 //
191 G4FieldManager* fieldMgr= nullptr; 263 G4FieldManager* fieldMgr= nullptr;
192 G4bool fieldExertsForce = false ; 264 G4bool fieldExertsForce = false ;
193 265
194 const G4Field* ptrField= nullptr; 266 const G4Field* ptrField= nullptr;
195 267
196 fieldMgr = fFieldPropagator->FindAndSetField 268 fieldMgr = fFieldPropagator->FindAndSetFieldManager( track.GetVolume() );
197 G4bool eligibleEM = (particleCharge != 0.0) 269 G4bool eligibleEM = (particleCharge != 0.0)
198 || ( fUseMagneticMoment && 270 || ( fUseMagneticMoment && (magneticMoment != 0.0) );
199 G4bool eligibleGrav = fUseGravity && (restM 271 G4bool eligibleGrav = fUseGravity && (restMass != 0.0) ;
200 272
201 if( (fieldMgr!=nullptr) && (eligibleEM||elig 273 if( (fieldMgr!=nullptr) && (eligibleEM||eligibleGrav) )
202 { 274 {
203 // Message the field Manager, to configur 275 // Message the field Manager, to configure it for this track
204 // 276 //
205 fieldMgr->ConfigureForTrack( &track ); 277 fieldMgr->ConfigureForTrack( &track );
206 278
207 // The above call can transition from a n 279 // The above call can transition from a null field-ptr oto a finite field.
208 // If the field manager has no field ptr, 280 // If the field manager has no field ptr, the field is zero
209 // by definition ( = there is no field ! 281 // by definition ( = there is no field ! )
210 // 282 //
211 ptrField= fieldMgr->GetDetectorField(); 283 ptrField= fieldMgr->GetDetectorField();
212 284
213 if( ptrField != nullptr) 285 if( ptrField != nullptr)
214 { 286 {
215 fieldExertsForce = eligibleEM 287 fieldExertsForce = eligibleEM
216 || ( eligibleGrav && ptrField->I 288 || ( eligibleGrav && ptrField->IsGravityActive() );
217 } 289 }
218 } 290 }
219 G4double momentumMagnitude = pParticle->GetT 291 G4double momentumMagnitude = pParticle->GetTotalMomentum() ;
220 292
221 if( fieldExertsForce ) 293 if( fieldExertsForce )
222 { 294 {
223 auto equationOfMotion= fFieldPropagator-> 295 auto equationOfMotion= fFieldPropagator->GetCurrentEquationOfMotion();
224 296
225 G4ChargeState chargeState(particleCharge, 297 G4ChargeState chargeState(particleCharge, // Charge can change (dynamic)
226 magneticMoment, 298 magneticMoment,
227 pParticleDef->G 299 pParticleDef->GetPDGSpin() );
228 if( equationOfMotion ) 300 if( equationOfMotion )
229 { 301 {
230 equationOfMotion->SetChargeMomentumMas 302 equationOfMotion->SetChargeMomentumMass( chargeState,
231 303 momentumMagnitude,
232 304 restMass );
233 } 305 }
234 #ifdef G4DEBUG_TRANSPORT 306 #ifdef G4DEBUG_TRANSPORT
235 else 307 else
236 { 308 {
237 G4cerr << " ERROR in G4CoupledTranspor 309 G4cerr << " ERROR in G4CoupledTransportation> "
238 << "Cannot find valid Equation 310 << "Cannot find valid Equation of motion: " << G4endl;
239 << " Unable to pass Charge, Mom 311 << " Unable to pass Charge, Momentum and Mass " << G4endl;
240 } 312 }
241 #endif 313 #endif
242 } 314 }
243 315
244 G4ThreeVector polarizationVec = track.GetPo 316 G4ThreeVector polarizationVec = track.GetPolarization() ;
245 G4FieldTrack aFieldTrack = G4FieldTrack(sta 317 G4FieldTrack aFieldTrack = G4FieldTrack(startPosition,
246 tra 318 track.GetGlobalTime(), // Lab.
247 tra 319 track.GetMomentumDirection(),
248 tra 320 track.GetKineticEnergy(),
249 res 321 restMass,
250 par 322 particleCharge,
251 pol 323 polarizationVec,
252 pPa 324 pParticleDef->GetPDGMagneticMoment(),
253 0.0 325 0.0, // Length along track
254 pPa 326 pParticleDef->GetPDGSpin() ) ;
255 G4int stepNo= track.GetCurrentStepNumber(); 327 G4int stepNo= track.GetCurrentStepNumber();
256 328
257 ELimited limitedStep; 329 ELimited limitedStep;
258 G4FieldTrack endTrackState('a'); // Defaul 330 G4FieldTrack endTrackState('a'); // Default values
259 331
260 fMassGeometryLimitedStep = false ; // de << 332 fMassGeometryLimitedStep = false ; // default
261 fGeometryLimitedStep = false; << 333 fAnyGeometryLimitedStep = false ;
262 if( currentMinimumStep > 0 ) 334 if( currentMinimumStep > 0 )
263 { 335 {
264 G4double newMassSafety= 0.0; // tem 336 G4double newMassSafety= 0.0; // temp. for recalculation
265 337
266 // Do the Transport in the field (non re 338 // Do the Transport in the field (non recti-linear)
267 // 339 //
268 lengthAlongCurve = fPathFinder->ComputeS 340 lengthAlongCurve = fPathFinder->ComputeStep( aFieldTrack,
269 341 currentMinimumStep,
270 << 342 fNavigatorId,
271 343 stepNo,
272 344 newMassSafety,
273 345 limitedStep,
274 346 endTrackState,
275 347 currentVolume ) ;
276 348
277 G4double newFullSafety= fPathFinder->Get 349 G4double newFullSafety= fPathFinder->GetCurrentSafety();
278 // this was estimated already in step 350 // this was estimated already in step above
279 351
280 if( limitedStep == kUnique || limitedSte 352 if( limitedStep == kUnique || limitedStep == kSharedTransport )
281 { 353 {
282 fMassGeometryLimitedStep = true ; 354 fMassGeometryLimitedStep = true ;
283 } 355 }
284 356
285 fGeometryLimitedStep = (fPathFinder->Get << 357 fAnyGeometryLimitedStep = (fPathFinder->GetNumberGeometriesLimitingStep() != 0);
286 358
287 #ifdef G4DEBUG_TRANSPORT 359 #ifdef G4DEBUG_TRANSPORT
288 if( fMassGeometryLimitedStep && !fGeomet << 360 if( fMassGeometryLimitedStep && !fAnyGeometryLimitedStep )
289 { 361 {
290 std::ostringstream message; 362 std::ostringstream message;
291 message << " ERROR in determining geom 363 message << " ERROR in determining geometries limiting the step" << G4endl;
292 message << " Limiting: mass=" << fMa 364 message << " Limiting: mass=" << fMassGeometryLimitedStep
293 << " any= " << fGeometryLimite << 365 << " any= " << fAnyGeometryLimitedStep << G4endl;
294 message << "Incompatible conditions - 366 message << "Incompatible conditions - by which geometries was it limited ?";
295 G4Exception("G4CoupledTransportation:: 367 G4Exception("G4CoupledTransportation::AlongStepGetPhysicalInteractionLength()",
296 "PathFinderConfused", Fata 368 "PathFinderConfused", FatalException, message);
297 } 369 }
298 #endif 370 #endif
299 371
300 geometryStepLength = std::min( lengthAlo 372 geometryStepLength = std::min( lengthAlongCurve, currentMinimumStep);
301 373
302 // Momentum: Magnitude and direction ca 374 // Momentum: Magnitude and direction can be changed too now ...
303 // 375 //
304 fMomentumChanged = true ; 376 fMomentumChanged = true ;
305 fTransportEndMomentumDir = endTrackState 377 fTransportEndMomentumDir = endTrackState.GetMomentumDir() ;
306 378
307 // Remember last safety origin & value. 379 // Remember last safety origin & value.
308 fPreviousSftOrigin = startPosition ; 380 fPreviousSftOrigin = startPosition ;
309 fPreviousMassSafety = newMassSafety ; 381 fPreviousMassSafety = newMassSafety ;
310 fPreviousFullSafety = newFullSafety ; 382 fPreviousFullSafety = newFullSafety ;
311 // fpSafetyHelper->SetCurrentSafety( new 383 // fpSafetyHelper->SetCurrentSafety( newFullSafety, startPosition);
312 384
313 #ifdef G4DEBUG_TRANSPORT 385 #ifdef G4DEBUG_TRANSPORT
314 if( verboseLevel > 1 ) 386 if( verboseLevel > 1 )
315 { 387 {
316 G4cout << "G4Transport:CompStep> " 388 G4cout << "G4Transport:CompStep> "
317 << " called the pathfinder for 389 << " called the pathfinder for a new step at " << startPosition
318 << " and obtained step = " << l 390 << " and obtained step = " << lengthAlongCurve << G4endl;
319 G4cout << " New safety (preStep) = " << 391 G4cout << " New safety (preStep) = " << newMassSafety
>> 392 << " versus precalculated = " << startMassSafety << G4endl;
320 } 393 }
321 #endif 394 #endif
322 395
323 // Store as best estimate value 396 // Store as best estimate value
>> 397 startMassSafety = newMassSafety ;
324 startFullSafety = newFullSafety ; 398 startFullSafety = newFullSafety ;
325 399
326 // Get the End-Position and End-Momentum 400 // Get the End-Position and End-Momentum (Dir-ection)
327 fTransportEndPosition = endTrackState.Ge 401 fTransportEndPosition = endTrackState.GetPosition() ;
328 fTransportEndKineticEnergy = endTrackSt 402 fTransportEndKineticEnergy = endTrackState.GetKineticEnergy() ;
329 } 403 }
330 else 404 else
331 { 405 {
332 geometryStepLength = lengthAlongCurve= 406 geometryStepLength = lengthAlongCurve= 0.0 ;
333 fMomentumChanged = false ; 407 fMomentumChanged = false ;
334 // fMassGeometryLimitedStep = false ; 408 // fMassGeometryLimitedStep = false ; // --- ???
335 // fGeometryLimitedStep = true; << 409 // fAnyGeometryLimitedStep = true;
336 fTransportEndMomentumDir = track.GetMome 410 fTransportEndMomentumDir = track.GetMomentumDirection();
337 fTransportEndKineticEnergy = track.GetK 411 fTransportEndKineticEnergy = track.GetKineticEnergy();
338 412
339 fTransportEndPosition = startPosition; 413 fTransportEndPosition = startPosition;
340 414
341 endTrackState= aFieldTrack; // Ensures 415 endTrackState= aFieldTrack; // Ensures that time is updated
>> 416
>> 417 // If the step length requested is 0, and we are on a boundary
>> 418 // then a boundary will also limit the step.
>> 419 if( startMassSafety == 0.0 )
>> 420 {
>> 421 fMassGeometryLimitedStep = true ;
>> 422 fAnyGeometryLimitedStep = true;
>> 423 }
>> 424 // TODO: Add explicit logical status for being at a boundary
342 } 425 }
343 // G4FieldTrack aTrackState(endTrackState); 426 // G4FieldTrack aTrackState(endTrackState);
344 427
345 if( !fieldExertsForce ) 428 if( !fieldExertsForce )
346 { 429 {
347 fParticleIsLooping = false ; 430 fParticleIsLooping = false ;
348 fMomentumChanged = false ; 431 fMomentumChanged = false ;
349 fEndGlobalTimeComputed = false ; 432 fEndGlobalTimeComputed = false ;
350 } 433 }
351 else 434 else
352 { 435 {
353 fParticleIsLooping = fFieldPropagator->I 436 fParticleIsLooping = fFieldPropagator->IsParticleLooping() ;
354 437
355 #ifdef G4DEBUG_TRANSPORT 438 #ifdef G4DEBUG_TRANSPORT
356 if( verboseLevel > 1 ) 439 if( verboseLevel > 1 )
357 { 440 {
358 G4cout << " G4CT::CS End Position = " 441 G4cout << " G4CT::CS End Position = "
359 << fTransportEndPosition << G4e 442 << fTransportEndPosition << G4endl;
360 G4cout << " G4CT::CS End Direction = " 443 G4cout << " G4CT::CS End Direction = "
361 << fTransportEndMomentumDir << 444 << fTransportEndMomentumDir << G4endl;
362 } 445 }
363 #endif 446 #endif
364 if( fFieldPropagator->GetCurrentFieldMan 447 if( fFieldPropagator->GetCurrentFieldManager()->DoesFieldChangeEnergy() )
365 { 448 {
366 // If the field can change energy, t 449 // If the field can change energy, then the time must be integrated
367 // - so this should have been upd 450 // - so this should have been updated
368 // 451 //
369 fCandidateEndGlobalTime = endTrack 452 fCandidateEndGlobalTime = endTrackState.GetLabTimeOfFlight();
370 fEndGlobalTimeComputed = true; 453 fEndGlobalTimeComputed = true;
371 454
372 // was ( fCandidateEndGlobalTime != 455 // was ( fCandidateEndGlobalTime != track.GetGlobalTime() );
373 // a cleaner way is to have FieldTra 456 // a cleaner way is to have FieldTrack knowing whether time
374 // is updated 457 // is updated
375 } 458 }
376 else 459 else
377 { 460 {
378 // The energy should be unchanged by 461 // The energy should be unchanged by field transport,
379 // - so the time changed will be 462 // - so the time changed will be calculated elsewhere
380 // 463 //
381 fEndGlobalTimeComputed = false; 464 fEndGlobalTimeComputed = false;
382 465
383 #ifdef G4VERBOSE <<
384 // Check that the integration preser 466 // Check that the integration preserved the energy
385 // - and if not correct this! 467 // - and if not correct this!
386 G4double startEnergy= track.GetKine 468 G4double startEnergy= track.GetKineticEnergy();
387 G4double endEnergy= fTransportEndKi 469 G4double endEnergy= fTransportEndKineticEnergy;
388 470
>> 471 static G4ThreadLocal G4int no_inexact_steps=0; // , no_large_ediff;
389 G4double absEdiff = std::fabs(startE 472 G4double absEdiff = std::fabs(startEnergy- endEnergy);
>> 473 if( absEdiff > perMillion * endEnergy )
>> 474 {
>> 475 no_inexact_steps++;
>> 476 // Possible statistics keeping here ...
>> 477 }
>> 478 #ifdef G4VERBOSE
390 if( (verboseLevel > 1) && ( absEdiff 479 if( (verboseLevel > 1) && ( absEdiff > perThousand * endEnergy) )
391 { 480 {
392 ReportInexactEnergy(startEnergy, e 481 ReportInexactEnergy(startEnergy, endEnergy);
393 } // end of if (verboseLevel) 482 } // end of if (verboseLevel)
394 #endif 483 #endif
395 // Correct the energy for fields tha 484 // Correct the energy for fields that conserve it
396 // This - hides the integration err 485 // This - hides the integration error
397 // - but gives a better physic 486 // - but gives a better physical answer
398 fTransportEndKineticEnergy= track.Ge 487 fTransportEndKineticEnergy= track.GetKineticEnergy();
399 } 488 }
400 } 489 }
401 490
402 fEndPointDistance = (fTransportEndPosition << 491 fEndpointDistance = (fTransportEndPosition - startPosition).mag() ;
403 fTransportEndSpin = endTrackState.GetSpin(); 492 fTransportEndSpin = endTrackState.GetSpin();
404 493
405 // Calculate the safety 494 // Calculate the safety
406 495
407 safetyProposal= startFullSafety; // used t 496 safetyProposal= startFullSafety; // used to be startMassSafety
408 // Changed to accomodate processes that ca 497 // Changed to accomodate processes that cannot update the safety
409 498
410 // Update safety for the end-point, if becom 499 // Update safety for the end-point, if becomes negative at the end-point.
411 500
412 if( (startFullSafety < fEndPointDistance ) << 501 if( (startFullSafety < fEndpointDistance )
413 && ( particleCharge != 0.0 ) ) // Onl 502 && ( particleCharge != 0.0 ) ) // Only needed to prepare for MSC
414 // && !fGeometryLimitedStep ) // To-Try: << 503 // && !fAnyGeometryLimitedStep ) // To-Try: No safety update if at boundary
415 { 504 {
416 G4double endFullSafety = 505 G4double endFullSafety =
417 fPathFinder->ComputeSafety( fTransport 506 fPathFinder->ComputeSafety( fTransportEndPosition);
418 // Expected mission -- only mass geome 507 // Expected mission -- only mass geometry's safety
419 // fLinearNavigator->ComputeSafety( << 508 // fMassNavigator->ComputeSafety( fTransportEndPosition) ;
420 // Yet discrete processes only have po 509 // Yet discrete processes only have poststep -- and this cannot
421 // currently revise the safety 510 // currently revise the safety
422 // ==> so we use the all-geometry sa 511 // ==> so we use the all-geometry safety as a precaution
423 512
424 fpSafetyHelper->SetCurrentSafety( endFul 513 fpSafetyHelper->SetCurrentSafety( endFullSafety, fTransportEndPosition);
425 // Pushing safety to Helper avoids rec 514 // Pushing safety to Helper avoids recalculation at this point
426 515
427 G4ThreeVector centerPt= G4ThreeVector(0. 516 G4ThreeVector centerPt= G4ThreeVector(0.0, 0.0, 0.0); // Used for return value
428 G4double endMassSafety= fPathFinder->Obt << 517 G4double endMassSafety= fPathFinder->ObtainSafety( fNavigatorId, centerPt);
429 // Retrieves the mass value from Path 518 // Retrieves the mass value from PathFinder (it calculated it)
430 519
431 fPreviousMassSafety = endMassSafety ; 520 fPreviousMassSafety = endMassSafety ;
432 fPreviousFullSafety = endFullSafety; 521 fPreviousFullSafety = endFullSafety;
433 fPreviousSftOrigin = fTransportEndPositi 522 fPreviousSftOrigin = fTransportEndPosition ;
434 523
435 // The convention (Stepping Manager's) i 524 // The convention (Stepping Manager's) is safety from the start point
436 // 525 //
437 safetyProposal = endFullSafety + fEndPoi << 526 safetyProposal = endFullSafety + fEndpointDistance;
438 // --> was endMassSafety 527 // --> was endMassSafety
439 // Changed to accomodate processes that 528 // Changed to accomodate processes that cannot update the safety
440 529
441 #ifdef G4DEBUG_TRANSPORT 530 #ifdef G4DEBUG_TRANSPORT
442 G4int prec= G4cout.precision(12) ; 531 G4int prec= G4cout.precision(12) ;
443 G4cout << "***CoupledTransportation::Alo 532 G4cout << "***CoupledTransportation::AlongStepGPIL ** " << G4endl ;
444 G4cout << " Revised Safety at endpoint 533 G4cout << " Revised Safety at endpoint " << fTransportEndPosition
445 << " give safety values: Mass= 534 << " give safety values: Mass= " << endMassSafety
446 << " All= " << endFullSafety << 535 << " All= " << endFullSafety << G4endl ;
447 G4cout << " Adding endpoint distance " << 536 G4cout << " Adding endpoint distance " << fEndpointDistance
448 << " to obtain pseudo-safety= " 537 << " to obtain pseudo-safety= " << safetyProposal << G4endl ;
449 G4cout.precision(prec); 538 G4cout.precision(prec);
450 } 539 }
451 else 540 else
452 { 541 {
453 G4int prec= G4cout.precision(12) ; 542 G4int prec= G4cout.precision(12) ;
454 G4cout << "***CoupledTransportation::Alo 543 G4cout << "***CoupledTransportation::AlongStepGPIL ** " << G4endl ;
455 G4cout << " Quick Safety estimate at en 544 G4cout << " Quick Safety estimate at endpoint "
456 << fTransportEndPosition 545 << fTransportEndPosition
457 << " gives safety endpoint valu 546 << " gives safety endpoint value = "
458 << startFullSafety - fEndPointDis << 547 << startFullSafety - fEndpointDistance
459 << " using start-point value " < 548 << " using start-point value " << startFullSafety
460 << " and endpointDistance " << f << 549 << " and endpointDistance " << fEndpointDistance << G4endl;
461 G4cout.precision(prec); 550 G4cout.precision(prec);
462 #endif 551 #endif
463 } 552 }
464 553
465 proposedSafetyForStart= safetyProposal; 554 proposedSafetyForStart= safetyProposal;
466 fParticleChange.ProposeTrueStepLength(geomet 555 fParticleChange.ProposeTrueStepLength(geometryStepLength) ;
467 556
468 return geometryStepLength ; 557 return geometryStepLength ;
469 } 558 }
470 559
>> 560 //////////////////////////////////////////////////////////////////////////
>> 561
>> 562 G4VParticleChange*
>> 563 G4CoupledTransportation::AlongStepDoIt( const G4Track& track,
>> 564 const G4Step& stepData )
>> 565 {
>> 566 static G4ThreadLocal G4long noCallsCT_ASDI=0;
>> 567 const char *methodName= "AlongStepDoIt";
>> 568
>> 569 noCallsCT_ASDI++;
>> 570
>> 571 fParticleChange.Initialize(track) ;
>> 572 // sets all its members to the value of corresponding members in G4Track
>> 573
>> 574 // Code specific for Transport
>> 575 //
>> 576 fParticleChange.ProposePosition(fTransportEndPosition) ;
>> 577 fParticleChange.ProposeMomentumDirection(fTransportEndMomentumDir) ;
>> 578 fParticleChange.ProposeEnergy(fTransportEndKineticEnergy) ;
>> 579 fParticleChange.SetMomentumChanged(fMomentumChanged) ;
>> 580
>> 581 fParticleChange.ProposePolarization(fTransportEndSpin);
>> 582
>> 583 G4double deltaTime = 0.0 ;
>> 584
>> 585 // Calculate Lab Time of Flight (ONLY if field Equations used it!)
>> 586 // G4double endTime = fCandidateEndGlobalTime;
>> 587 // G4double delta_time = endTime - startTime;
>> 588
>> 589 G4double startTime = track.GetGlobalTime() ;
>> 590
>> 591 if (!fEndGlobalTimeComputed)
>> 592 {
>> 593 G4double finalInverseVel= DBL_MAX, initialInverseVel=DBL_MAX;
>> 594
>> 595 // The time was not integrated .. make the best estimate possible
>> 596 //
>> 597 G4double finalVelocity = track.GetVelocity() ;
>> 598 if( finalVelocity > 0.0 ) { finalInverseVel= 1.0 / finalVelocity; }
>> 599 G4double initialVelocity = stepData.GetPreStepPoint()->GetVelocity() ;
>> 600 if( initialVelocity > 0.0 ) { initialInverseVel= 1.0 / initialVelocity; }
>> 601 G4double stepLength = track.GetStepLength() ;
>> 602
>> 603 if (finalVelocity > 0.0)
>> 604 {
>> 605 // deltaTime = stepLength/finalVelocity ;
>> 606 G4double meanInverseVelocity = 0.5
>> 607 * (initialInverseVel + finalInverseVel);
>> 608 deltaTime = stepLength * meanInverseVelocity ;
>> 609 }
>> 610 else
>> 611 {
>> 612 deltaTime = stepLength * initialInverseVel ;
>> 613 } // Could do with better estimate for final step (finalVelocity = 0) ?
>> 614
>> 615 fCandidateEndGlobalTime = startTime + deltaTime ;
>> 616 fParticleChange.ProposeLocalTime( track.GetLocalTime() + deltaTime) ;
>> 617 }
>> 618 else
>> 619 {
>> 620 deltaTime = fCandidateEndGlobalTime - startTime ;
>> 621 fParticleChange.ProposeGlobalTime( fCandidateEndGlobalTime ) ;
>> 622 }
>> 623
>> 624 // Now Correct by Lorentz factor to get "proper" deltaTime
>> 625
>> 626 G4double restMass = track.GetDynamicParticle()->GetMass() ;
>> 627 G4double deltaProperTime = deltaTime*( restMass/track.GetTotalEnergy() ) ;
>> 628
>> 629 fParticleChange.ProposeProperTime(track.GetProperTime() + deltaProperTime) ;
>> 630 // fParticleChange. ProposeTrueStepLength( track.GetStepLength() ) ;
>> 631
>> 632 // If the particle is caught looping or is stuck (in very difficult
>> 633 // boundaries) in a magnetic field (doing many steps) THEN this kills it ...
>> 634 //
>> 635 if ( fParticleIsLooping )
>> 636 {
>> 637 G4double endEnergy= fTransportEndKineticEnergy;
>> 638
>> 639 const G4ParticleDefinition* particleType=
>> 640 track.GetDynamicParticle() -> GetParticleDefinition();
>> 641 G4bool stable = particleType->GetPDGStable();
>> 642
>> 643 G4bool candidateForEnd = (endEnergy < fThreshold_Important_Energy)
>> 644 || (fNoLooperTrials >= fThresholdTrials) ;
>> 645
>> 646 if( candidateForEnd && stable )
>> 647 {
>> 648 const G4int electronPDG= 11; // G4Electron::G4Electron()->GetPDGEncoding();
>> 649 G4int particlePDG= particleType->GetPDGEncoding();
>> 650
>> 651 // Kill the looping particle
>> 652 //
>> 653 fParticleChange.ProposeTrackStatus( fStopAndKill ) ;
>> 654
>> 655 // Simple statistics
>> 656 fSumEnergyKilled += endEnergy;
>> 657 fSumEnerSqKilled = endEnergy * endEnergy;
>> 658 fNumLoopersKilled++;
>> 659
>> 660 if( endEnergy > fMaxEnergyKilled ) {
>> 661 fMaxEnergyKilled = endEnergy;
>> 662 fMaxEnergyKilledPDG = particlePDG;
>> 663 }
>> 664
>> 665 if( particleType->GetPDGEncoding() != electronPDG )
>> 666 {
>> 667 fSumEnergyKilled_NonElectron += endEnergy;
>> 668 fSumEnerSqKilled_NonElectron += endEnergy * endEnergy;
>> 669 fNumLoopersKilled_NonElectron++;
>> 670
>> 671 if( endEnergy > fMaxEnergyKilled_NonElectron )
>> 672 {
>> 673 fMaxEnergyKilled_NonElectron = endEnergy;
>> 674 fMaxEnergyKilled_NonElecPDG = particlePDG;
>> 675 }
>> 676 }
>> 677
>> 678 if( endEnergy > fThreshold_Warning_Energy && ! fSilenceLooperWarnings )
>> 679 {
>> 680 fpLogger->ReportLoopingTrack( track, stepData, fNoLooperTrials,
>> 681 noCallsCT_ASDI, methodName );
>> 682 }
>> 683
>> 684 fNoLooperTrials=0;
>> 685 }
>> 686 else
>> 687 {
>> 688 fNoLooperTrials ++;
>> 689
>> 690 fMaxEnergySaved = std::max( endEnergy, fMaxEnergySaved);
>> 691 if( fNoLooperTrials == 1 ) {
>> 692 fSumEnergySaved += endEnergy;
>> 693 if ( !stable )
>> 694 fSumEnergyUnstableSaved += endEnergy;
>> 695 }
>> 696 #ifdef G4VERBOSE
>> 697 if( verboseLevel > 2 && ! fSilenceLooperWarnings )
>> 698 {
>> 699 G4cout << " ** G4CoupledTransportation::AlongStepDoIt():"
>> 700 << " Particle is looping but is saved ..." << G4endl
>> 701 << " Number of trials (this track) = " << fNoLooperTrials
>> 702 << G4endl
>> 703 << " Steps by this track: " << track.GetCurrentStepNumber()
>> 704 << G4endl
>> 705 << " Total no of calls to this method (all tracks) = "
>> 706 << noCallsCT_ASDI << G4endl;
>> 707 }
>> 708 #endif
>> 709 }
>> 710 }
>> 711 else
>> 712 {
>> 713 fNoLooperTrials=0;
>> 714 }
>> 715
>> 716 // Another (sometimes better way) is to use a user-limit maximum Step size
>> 717 // to alleviate this problem ..
>> 718
>> 719 // Introduce smooth curved trajectories to particle-change
>> 720 //
>> 721 fParticleChange.SetPointerToVectorOfAuxiliaryPoints
>> 722 (fFieldPropagator->GimmeTrajectoryVectorAndForgetIt() );
>> 723
>> 724 return &fParticleChange ;
>> 725 }
>> 726
>> 727 //////////////////////////////////////////////////////////////////////////
>> 728 //
>> 729 // This ensures that the PostStep action is always called,
>> 730 // so that it can do the relocation if it is needed.
>> 731 //
>> 732
>> 733 G4double G4CoupledTransportation::
>> 734 PostStepGetPhysicalInteractionLength( const G4Track&,
>> 735 G4double, // previousStepSize
>> 736 G4ForceCondition* pForceCond )
>> 737 {
>> 738 // Must act as PostStep action -- to relocate particle
>> 739 *pForceCond = Forced ;
>> 740 return DBL_MAX ;
>> 741 }
>> 742
471 ////////////////////////////////////////////// 743 /////////////////////////////////////////////////////////////////////////////
472 744
473 void G4CoupledTransportation:: 745 void G4CoupledTransportation::
474 ReportMove( G4ThreeVector OldVector, G4ThreeVe 746 ReportMove( G4ThreeVector OldVector, G4ThreeVector NewVector,
475 const G4String& Quantity ) 747 const G4String& Quantity )
476 { 748 {
477 G4ThreeVector moveVec = ( NewVector - OldV 749 G4ThreeVector moveVec = ( NewVector - OldVector );
478 750
479 G4cerr << G4endl 751 G4cerr << G4endl
480 << "******************************* 752 << "**************************************************************"
481 << G4endl; 753 << G4endl;
482 G4cerr << "Endpoint has moved between valu 754 G4cerr << "Endpoint has moved between value expected from TransportEndPosition "
483 << " and value from Track in PostSt 755 << " and value from Track in PostStepDoIt. " << G4endl
484 << "Change of " << Quantity << " is 756 << "Change of " << Quantity << " is " << moveVec.mag() / mm
485 << " mm long, " 757 << " mm long, "
486 << " and its vector is " << (1.0/mm 758 << " and its vector is " << (1.0/mm) * moveVec << " mm " << G4endl
487 << "Endpoint of ComputeStep was " < 759 << "Endpoint of ComputeStep was " << OldVector
488 << " and current position to locate 760 << " and current position to locate is " << NewVector << G4endl;
489 } 761 }
490 762
491 ////////////////////////////////////////////// 763 /////////////////////////////////////////////////////////////////////////////
492 764
493 G4VParticleChange* G4CoupledTransportation::Po 765 G4VParticleChange* G4CoupledTransportation::PostStepDoIt( const G4Track& track,
494 766 const G4Step& )
495 { 767 {
496 G4TouchableHandle retCurrentTouchable ; // 768 G4TouchableHandle retCurrentTouchable ; // The one to return
497 769
498 // Initialize ParticleChange (by setting al 770 // Initialize ParticleChange (by setting all its members equal
499 // to correspond 771 // to corresponding members in G4Track)
500 // fParticleChange.Initialize(track) ; // T 772 // fParticleChange.Initialize(track) ; // To initialise TouchableChange
501 773
502 fParticleChange.ProposeTrackStatus(track.Get 774 fParticleChange.ProposeTrackStatus(track.GetTrackStatus()) ;
503 775
504 if( fSignifyStepInAnyVolume ) 776 if( fSignifyStepInAnyVolume )
505 { 777 {
506 fParticleChange.ProposeFirstStepInVolume( << 778 fParticleChange.ProposeFirstStepInVolume( fFirstStepInAnyVolume );
507 } 779 }
508 else 780 else
509 { 781 {
510 fParticleChange.ProposeFirstStepInVolume( 782 fParticleChange.ProposeFirstStepInVolume( fFirstStepInMassVolume );
511 } 783 }
512 784
513 // Check that the end position and direction 785 // Check that the end position and direction are preserved
514 // since call to AlongStepDoIt 786 // since call to AlongStepDoIt
515 787
516 #ifdef G4DEBUG_TRANSPORT 788 #ifdef G4DEBUG_TRANSPORT
517 if( ( verboseLevel > 0 ) 789 if( ( verboseLevel > 0 )
518 && ((fTransportEndPosition - track.GetPos 790 && ((fTransportEndPosition - track.GetPosition()).mag2() >= 1.0e-16) )
519 { 791 {
520 ReportMove( track.GetPosition(), fTranspo 792 ReportMove( track.GetPosition(), fTransportEndPosition,
521 "End of Step Position" ); 793 "End of Step Position" );
522 G4cerr << " Problem in G4CoupledTransport 794 G4cerr << " Problem in G4CoupledTransportation::PostStepDoIt " << G4endl;
523 } 795 }
524 796
525 // If the Step was determined by the volume 797 // If the Step was determined by the volume boundary, relocate the particle
526 // The pathFinder will know that the geometr 798 // The pathFinder will know that the geometry limited the step (!?)
527 799
528 if( verboseLevel > 0 ) 800 if( verboseLevel > 0 )
529 { 801 {
530 G4cout << " Calling PathFinder::Locate() 802 G4cout << " Calling PathFinder::Locate() from "
531 << " G4CoupledTransportation::Post 803 << " G4CoupledTransportation::PostStepDoIt " << G4endl;
532 G4cout << " fGeometryLimitedStep is " < << 804 G4cout << " fAnyGeometryLimitedStep is " << fAnyGeometryLimitedStep
>> 805 << G4endl;
>> 806
533 } 807 }
534 #endif 808 #endif
535 809
536 if(fGeometryLimitedStep) << 810 if(fAnyGeometryLimitedStep)
537 { 811 {
538 fPathFinder->Locate( track.GetPosition(), 812 fPathFinder->Locate( track.GetPosition(),
539 track.GetMomentumDire 813 track.GetMomentumDirection(),
540 true); 814 true);
541 815
542 // fCurrentTouchable will now become the p 816 // fCurrentTouchable will now become the previous touchable,
543 // and what was the previous will be freed 817 // and what was the previous will be freed.
544 // (Needed because the preStepPoint can po 818 // (Needed because the preStepPoint can point to the previous touchable)
545 819
546 fCurrentTouchableHandle= 820 fCurrentTouchableHandle=
547 fPathFinder->CreateTouchableHandle( G4Tr << 821 fPathFinder->CreateTouchableHandle( fNavigatorId );
548 822
549 #ifdef G4DEBUG_TRANSPORT 823 #ifdef G4DEBUG_TRANSPORT
>> 824 if( verboseLevel > 0 )
>> 825 {
>> 826 G4cout << "G4CoupledTransportation::PostStepDoIt --- fNavigatorId = "
>> 827 << fNavigatorId << G4endl;
>> 828 }
550 if( verboseLevel > 1 ) 829 if( verboseLevel > 1 )
551 { 830 {
552 G4VPhysicalVolume* vol= fCurrentTouchab 831 G4VPhysicalVolume* vol= fCurrentTouchableHandle->GetVolume();
553 G4cout << "CHECK !!!!!!!!!!! fCurrentTo 832 G4cout << "CHECK !!!!!!!!!!! fCurrentTouchableHandle->GetVolume() = "
554 << vol; 833 << vol;
555 if( vol ) { G4cout << "Name=" << vol->G 834 if( vol ) { G4cout << "Name=" << vol->GetName(); }
556 G4cout << G4endl; 835 G4cout << G4endl;
557 } 836 }
558 #endif 837 #endif
559 838
560 // Check whether the particle is out of th 839 // Check whether the particle is out of the world volume
561 // If so it has exited and must be killed. 840 // If so it has exited and must be killed.
562 // 841 //
563 if( fCurrentTouchableHandle->GetVolume() = 842 if( fCurrentTouchableHandle->GetVolume() == 0 )
564 { 843 {
565 fParticleChange.ProposeTrackStatus( fSt 844 fParticleChange.ProposeTrackStatus( fStopAndKill ) ;
566 } 845 }
567 retCurrentTouchable = fCurrentTouchableHan 846 retCurrentTouchable = fCurrentTouchableHandle ;
568 // fParticleChange.SetTouchableHandle( fCu 847 // fParticleChange.SetTouchableHandle( fCurrentTouchableHandle ) ;
569 } 848 }
570 else // fGeometryLimitedStep is false << 849 else // fAnyGeometryLimitedStep is false
571 { 850 {
572 #ifdef G4DEBUG_TRANSPORT 851 #ifdef G4DEBUG_TRANSPORT
573 if( verboseLevel > 1 ) 852 if( verboseLevel > 1 )
574 { 853 {
575 G4cout << "G4CoupledTransportation::Post << 854 G4cout << "G4CoupledTransportation::PostStepDoIt -- "
576 << " fGeometryLimitedStep = " << << 855 << " fAnyGeometryLimitedStep = " << fAnyGeometryLimitedStep
577 << " must be false " << G4endl; << 856 << " must be false " << G4endl;
578 } 857 }
579 #endif 858 #endif
580 // This serves only to move each of the Na 859 // This serves only to move each of the Navigator's location
581 // 860 //
582 // fLinearNavigator->LocateGlobalPointWith 861 // fLinearNavigator->LocateGlobalPointWithinVolume( track.GetPosition() ) ;
583 862
584 fPathFinder->ReLocate( track.GetPosition() 863 fPathFinder->ReLocate( track.GetPosition() );
585 // track.GetMomentu 864 // track.GetMomentumDirection() );
586 865
587 // Keep the value of the track's current T 866 // Keep the value of the track's current Touchable is retained,
588 // and use it to overwrite the (unset) on 867 // and use it to overwrite the (unset) one in particle change.
589 // Expect this must be fCurrentTouchable t 868 // Expect this must be fCurrentTouchable too
590 // - could it be different, eg at the st 869 // - could it be different, eg at the start of a step ?
591 // 870 //
592 retCurrentTouchable = track.GetTouchableHa 871 retCurrentTouchable = track.GetTouchableHandle() ;
593 // fParticleChange.SetTouchableHandle( tra 872 // fParticleChange.SetTouchableHandle( track.GetTouchableHandle() ) ;
594 } // endif ( fGeometryLimitedStep ) << 873 } // endif ( fAnyGeometryLimitedStep )
595 874
596 #ifdef G4DEBUG_NAVIGATION 875 #ifdef G4DEBUG_NAVIGATION
597 G4cout << " CoupledTransport::AlongStep GPI 876 G4cout << " CoupledTransport::AlongStep GPIL: "
598 << " last-step: any= " << fGeometryL << 877 << " last-step: any= " << fAnyGeometryLimitedStep
599 << " mass= " << fMassGeometryLimitedS << 878 << " . ..... x . "
>> 879 << " mass= " << fMassGeometryLimitedStep
>> 880 << G4endl;
600 #endif 881 #endif
601 882
602 if( fSignifyStepInAnyVolume ) 883 if( fSignifyStepInAnyVolume )
603 fParticleChange.ProposeLastStepInVolume(fG << 884 fParticleChange.ProposeLastStepInVolume(fAnyGeometryLimitedStep);
604 else 885 else
605 fParticleChange.ProposeLastStepInVolume(f 886 fParticleChange.ProposeLastStepInVolume(fMassGeometryLimitedStep);
606 887
607 SetTouchableInformation(retCurrentTouchable) << 888 const G4VPhysicalVolume* pNewVol = retCurrentTouchable->GetVolume() ;
>> 889 const G4Material* pNewMaterial = 0 ;
>> 890 const G4VSensitiveDetector* pNewSensitiveDetector = 0 ;
>> 891
>> 892 if( pNewVol != 0 )
>> 893 {
>> 894 pNewMaterial= pNewVol->GetLogicalVolume()->GetMaterial();
>> 895 pNewSensitiveDetector= pNewVol->GetLogicalVolume()->GetSensitiveDetector();
>> 896 }
>> 897
>> 898 // ( const_cast<G4Material *> pNewMaterial ) ;
>> 899 // ( const_cast<G4VSensitiveDetetor *> pNewSensitiveDetector) ;
>> 900
>> 901 fParticleChange.SetMaterialInTouchable( (G4Material *) pNewMaterial ) ;
>> 902 fParticleChange.SetSensitiveDetectorInTouchable( (G4VSensitiveDetector *) pNewSensitiveDetector ) ;
>> 903 // "temporarily" until Get/Set Material of ParticleChange,
>> 904 // and StepPoint can be made const.
>> 905
>> 906 const G4MaterialCutsCouple* pNewMaterialCutsCouple = 0;
>> 907 if( pNewVol != 0 )
>> 908 {
>> 909 pNewMaterialCutsCouple=pNewVol->GetLogicalVolume()->GetMaterialCutsCouple();
>> 910 if( pNewMaterialCutsCouple!=0
>> 911 && pNewMaterialCutsCouple->GetMaterial()!=pNewMaterial )
>> 912 {
>> 913 // for parametrized volume
>> 914 //
>> 915 pNewMaterialCutsCouple = G4ProductionCutsTable::GetProductionCutsTable()
>> 916 ->GetMaterialCutsCouple(pNewMaterial,
>> 917 pNewMaterialCutsCouple->GetProductionCuts());
>> 918 }
>> 919 }
>> 920 fParticleChange.SetMaterialCutsCoupleInTouchable( pNewMaterialCutsCouple );
>> 921
>> 922 // Must always set the touchable in ParticleChange, whether relocated or not
>> 923 //
>> 924 fParticleChange.SetTouchableHandle(retCurrentTouchable) ;
608 925
609 return &fParticleChange ; 926 return &fParticleChange ;
610 } 927 }
611 928
612 ////////////////////////////////////////////// 929 /////////////////////////////////////////////////////////////////////////////
613 // New method takes over the responsibility to 930 // New method takes over the responsibility to reset the state of
614 // G4CoupledTransportation object: 931 // G4CoupledTransportation object:
615 // - at the start of a new track, and 932 // - at the start of a new track, and
616 // - on the resumption of a suspended tra 933 // - on the resumption of a suspended track.
617 // 934 //
618 void 935 void
619 G4CoupledTransportation::StartTracking(G4Track 936 G4CoupledTransportation::StartTracking(G4Track* aTrack)
620 { 937 {
621 G4Transportation::StartTracking(aTrack); << 938
>> 939 G4TransportationManager* transportMgr =
>> 940 G4TransportationManager::GetTransportationManager();
>> 941
>> 942 // G4VProcess::StartTracking(aTrack);
>> 943 fNewTrack= true;
622 944
>> 945 // The 'initialising' actions
>> 946 // once taken in AlongStepGPIL -- if ( track.GetCurrentStepNumber()==1 )
>> 947
>> 948 // fStartedNewTrack= true;
>> 949
>> 950 fMassNavigator = transportMgr->GetNavigatorForTracking() ;
>> 951 fNavigatorId= transportMgr->ActivateNavigator( fMassNavigator );
>> 952
623 G4ThreeVector position = aTrack->GetPosition 953 G4ThreeVector position = aTrack->GetPosition();
624 G4ThreeVector direction = aTrack->GetMomentu 954 G4ThreeVector direction = aTrack->GetMomentumDirection();
625 955
626 fPathFinder->PrepareNewTrack( position, dire 956 fPathFinder->PrepareNewTrack( position, direction);
627 // This implies a call to fPathFinder->Locat 957 // This implies a call to fPathFinder->Locate( position, direction );
628 958
>> 959 // Whether field exists should be determined at run level -- TODO
>> 960 fAnyFieldExists= DoesAnyFieldExist();
>> 961
629 // reset safety value and center 962 // reset safety value and center
630 // 963 //
631 fPreviousMassSafety = 0.0 ; 964 fPreviousMassSafety = 0.0 ;
632 fPreviousFullSafety = 0.0 ; 965 fPreviousFullSafety = 0.0 ;
633 fPreviousSftOrigin = G4ThreeVector(0.,0.,0.) 966 fPreviousSftOrigin = G4ThreeVector(0.,0.,0.) ;
>> 967
>> 968 // reset looping counter -- for motion in field
>> 969 fNoLooperTrials= 0;
>> 970
>> 971 // Must clear this state .. else it depends on last track's value
>> 972 // --> a better solution would set this from state of suspended track TODO ?
>> 973 // Was if( aTrack->GetCurrentStepNumber()==1 ) { .. }
>> 974
>> 975 // ChordFinder reset internal state
>> 976 //
>> 977 if( fFieldPropagator && fAnyFieldExists )
>> 978 {
>> 979 fFieldPropagator->ClearPropagatorState();
>> 980 // Resets safety values, in case of overlaps.
>> 981
>> 982 G4ChordFinder* chordF= fFieldPropagator->GetChordFinder();
>> 983 if( chordF ) { chordF->ResetStepEstimate(); }
>> 984 }
>> 985
>> 986 // Clear the chord finders of all fields (ie managers) derived objects
>> 987 //
>> 988 G4FieldManagerStore* fieldMgrStore = G4FieldManagerStore::GetInstance();
>> 989 fieldMgrStore->ClearAllChordFindersState();
>> 990
>> 991 #ifdef G4DEBUG_TRANSPORT
>> 992 if( verboseLevel > 1 )
>> 993 {
>> 994 G4cout << " Returning touchable handle " << fCurrentTouchableHandle
>> 995 << G4endl;
>> 996 }
>> 997 #endif
>> 998
>> 999 // Update the current touchable handle (from the track's)
>> 1000 //
>> 1001 fCurrentTouchableHandle = aTrack->GetTouchableHandle();
634 } 1002 }
635 1003
636 ////////////////////////////////////////////// 1004 /////////////////////////////////////////////////////////////////////////////
637 1005
638 void 1006 void
639 G4CoupledTransportation::EndTracking() 1007 G4CoupledTransportation::EndTracking()
640 { 1008 {
641 G4TransportationManager::GetTransportationMa 1009 G4TransportationManager::GetTransportationManager()->InactivateAll();
642 fPathFinder->EndTrack(); 1010 fPathFinder->EndTrack();
643 // Resets TransportationManager to use ord 1011 // Resets TransportationManager to use ordinary Navigator
644 } 1012 }
645 1013
646 ////////////////////////////////////////////// 1014 /////////////////////////////////////////////////////////////////////////////
647 1015
648 void 1016 void
649 G4CoupledTransportation:: 1017 G4CoupledTransportation::
650 ReportInexactEnergy(G4double startEnergy, G4do 1018 ReportInexactEnergy(G4double startEnergy, G4double endEnergy)
651 { 1019 {
652 static G4ThreadLocal G4int no_warnings= 0, w 1020 static G4ThreadLocal G4int no_warnings= 0, warnModulo=1,
653 moduloFactor= 10, 1021 moduloFactor= 10, no_large_ediff= 0;
654 1022
655 if( std::fabs(startEnergy- endEnergy) > perT 1023 if( std::fabs(startEnergy- endEnergy) > perThousand * endEnergy )
656 { 1024 {
657 no_large_ediff ++; 1025 no_large_ediff ++;
658 if( (no_large_ediff% warnModulo) == 0 ) 1026 if( (no_large_ediff% warnModulo) == 0 )
659 { 1027 {
660 no_warnings++; 1028 no_warnings++;
661 std::ostringstream message; 1029 std::ostringstream message;
662 message << "Energy change in Step is abo 1030 message << "Energy change in Step is above 1^-3 relative value. "
663 << G4endl 1031 << G4endl
664 << " Relative change in 'track 1032 << " Relative change in 'tracking' step = "
665 << std::setw(15) << (endEnergy-s 1033 << std::setw(15) << (endEnergy-startEnergy)/startEnergy
666 << G4endl 1034 << G4endl
667 << " Starting E= " << std::set 1035 << " Starting E= " << std::setw(12) << startEnergy / MeV
668 << " MeV " << G4endl 1036 << " MeV " << G4endl
669 << " Ending E= " << std::set 1037 << " Ending E= " << std::setw(12) << endEnergy / MeV
670 << " MeV " << G4endl 1038 << " MeV " << G4endl
671 << "Energy has been corrected -- 1039 << "Energy has been corrected -- however, review"
672 << " field propagation parameter 1040 << " field propagation parameters for accuracy." << G4endl;
673 if ( (verboseLevel > 2 ) || (no_warnings 1041 if ( (verboseLevel > 2 ) || (no_warnings<4)
674 || (no_large_ediff == warnModulo * mod 1042 || (no_large_ediff == warnModulo * moduloFactor) )
675 { 1043 {
676 message << "These include EpsilonStepM 1044 message << "These include EpsilonStepMax(/Min) in G4FieldManager,"
677 << G4endl 1045 << G4endl
678 << "which determine fractional 1046 << "which determine fractional error per step for integrated quantities."
679 << G4endl 1047 << G4endl
680 << "Note also the influence of 1048 << "Note also the influence of the permitted number of integration steps."
681 << G4endl; 1049 << G4endl;
682 } 1050 }
683 message << "Bad 'endpoint'. Energy chang 1051 message << "Bad 'endpoint'. Energy change detected and corrected."
684 << G4endl 1052 << G4endl
685 << "Has occurred already " << no 1053 << "Has occurred already " << no_large_ediff << " times.";
686 G4Exception("G4CoupledTransportation::Al 1054 G4Exception("G4CoupledTransportation::AlongStepGetPIL()",
687 "EnergyChange", JustWarning, 1055 "EnergyChange", JustWarning, message);
688 if( no_large_ediff == warnModulo * modul 1056 if( no_large_ediff == warnModulo * moduloFactor )
689 { 1057 {
690 warnModulo *= moduloFactor; 1058 warnModulo *= moduloFactor;
691 } 1059 }
692 } 1060 }
693 } 1061 }
>> 1062 }
>> 1063
>> 1064 /////////////////////////////////////////////////////////////////////////////
>> 1065
>> 1066 G4bool G4CoupledTransportation::EnableMagneticMoment(G4bool useMoment)
>> 1067 {
>> 1068 G4bool lastValue= fUseMagneticMoment;
>> 1069 fUseMagneticMoment= useMoment;
>> 1070 G4Transportation::fUseMagneticMoment= useMoment;
>> 1071 return lastValue;
>> 1072 }
>> 1073
>> 1074 /////////////////////////////////////////////////////////////////////////////
>> 1075 //
>> 1076
>> 1077 G4bool G4CoupledTransportation::EnableGravity(G4bool useGravity)
>> 1078 {
>> 1079 G4bool lastValue= fUseGravity;
>> 1080 fUseGravity= useGravity;
>> 1081 G4Transportation::fUseGravity= useGravity;
>> 1082 return lastValue;
>> 1083 }
>> 1084
>> 1085 /////////////////////////////////////////////////////////////////////////////
>> 1086 //
>> 1087 // Supress (or not) warnings about 'looping' particles
>> 1088
>> 1089 void G4CoupledTransportation::SetSilenceLooperWarnings( G4bool val)
>> 1090 {
>> 1091 fSilenceLooperWarnings= val; // Flag to *Supress* all 'looper' warnings
>> 1092 // G4CoupledTransportation::fSilenceLooperWarnings= val;
>> 1093 }
>> 1094
>> 1095 /////////////////////////////////////////////////////////////////////////////
>> 1096 //
>> 1097 G4bool G4CoupledTransportation::GetSilenceLooperWarnings()
>> 1098 {
>> 1099 return fSilenceLooperWarnings;
>> 1100 }
>> 1101
>> 1102 /////////////////////////////////////////////////////////////////////////////
>> 1103 //
>> 1104 void G4CoupledTransportation::SetHighLooperThresholds()
>> 1105 {
>> 1106 // Setting old 'high' values for thresholds - old values, potentially appropriate
>> 1107 // for the energy frontier HEP experiments
>> 1108 SetThresholdWarningEnergy( 100.0 * CLHEP::MeV ); // Warn above this energy
>> 1109 SetThresholdImportantEnergy( 250.0 * CLHEP::MeV ); // Give a few trial above this E);
>> 1110
>> 1111 G4int maxTrials = 10;
>> 1112 SetThresholdTrials( maxTrials );
>> 1113
>> 1114 if( verboseLevel ) ReportLooperThresholds();
>> 1115 }
>> 1116
>> 1117 /////////////////////////////////////////////////////////////////////////////
>> 1118 void G4CoupledTransportation::SetLowLooperThresholds() // Values for low-E applications
>> 1119 {
>> 1120 // These values were the default in Geant4 10.5 - beta
>> 1121 SetThresholdWarningEnergy( 1.0 * CLHEP::keV ); // Warn above this En
>> 1122 SetThresholdImportantEnergy( 1.0 * CLHEP::MeV ); // Extra trials above it
>> 1123
>> 1124 G4int maxTrials = 30; // A new value - was 10
>> 1125 SetThresholdTrials( maxTrials );
>> 1126
>> 1127 if( verboseLevel ) ReportLooperThresholds();
>> 1128 }
>> 1129
>> 1130 /////////////////////////////////////////////////////////////////////////////
>> 1131 //
>> 1132 void
>> 1133 G4CoupledTransportation::ReportMissingLogger( const char* methodName )
>> 1134 {
>> 1135 const char* message= "Logger object missing from G4CoupledTransportation";
>> 1136 G4String classAndMethod= G4String("G4CoupledTransportation") + G4String( methodName );
>> 1137 G4Exception(classAndMethod, "Missing Logger", JustWarning, message);
>> 1138
>> 1139 if( verboseLevel ) ReportLooperThresholds();
>> 1140 }
>> 1141
>> 1142 /////////////////////////////////////////////////////////////////////////////
>> 1143 //
>> 1144 void
>> 1145 G4CoupledTransportation::ReportLooperThresholds()
>> 1146 {
>> 1147 PushThresholdsToLogger(); // To be absolutely certain they are in sync
>> 1148 fpLogger->ReportLooperThresholds("G4CoupledTransportation");
694 } 1149 }
695 1150