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