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