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