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