1 | // |
---|
2 | // ******************************************************************** |
---|
3 | // * License and Disclaimer * |
---|
4 | // * * |
---|
5 | // * The Geant4 software is copyright of the Copyright Holders of * |
---|
6 | // * the Geant4 Collaboration. It is provided under the terms and * |
---|
7 | // * conditions of the Geant4 Software License, included in the file * |
---|
8 | // * LICENSE and available at http://cern.ch/geant4/license . These * |
---|
9 | // * include a list of copyright holders. * |
---|
10 | // * * |
---|
11 | // * Neither the authors of this software system, nor their employing * |
---|
12 | // * institutes,nor the agencies providing financial support for this * |
---|
13 | // * work make any representation or warranty, express or implied, * |
---|
14 | // * regarding this software system or assume any liability for its * |
---|
15 | // * use. Please see the license in the file LICENSE and URL above * |
---|
16 | // * for the full disclaimer and the limitation of liability. * |
---|
17 | // * * |
---|
18 | // * This code implementation is the result of the scientific and * |
---|
19 | // * technical work of the GEANT4 collaboration. * |
---|
20 | // * By using, copying, modifying or distributing the software (or * |
---|
21 | // * any work based on the software) you agree to acknowledge its * |
---|
22 | // * use in resulting scientific publications, and indicate your * |
---|
23 | // * acceptance of all terms of the Geant4 Software license. * |
---|
24 | // ******************************************************************** |
---|
25 | // |
---|
26 | // |
---|
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-04-beta-01 $ |
---|
30 | // ------------------------------------------------------------ |
---|
31 | // GEANT 4 class implementation |
---|
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 |
---|
45 | // Created: 19 March 1997, J. Apostolakis |
---|
46 | // ======================================================================= |
---|
47 | |
---|
48 | #include "G4CoupledTransportation.hh" |
---|
49 | #include "G4ProductionCutsTable.hh" |
---|
50 | #include "G4ParticleTable.hh" |
---|
51 | #include "G4ChordFinder.hh" |
---|
52 | #include "G4FieldManagerStore.hh" |
---|
53 | class G4VSensitiveDetector; |
---|
54 | |
---|
55 | ////////////////////////////////////////////////////////////////////////// |
---|
56 | // |
---|
57 | // Constructor |
---|
58 | |
---|
59 | G4CoupledTransportation::G4CoupledTransportation( G4int verboseLevel ) |
---|
60 | : G4VProcess( G4String("CoupledTransportation"), fTransportation ), |
---|
61 | fParticleIsLooping( false ), |
---|
62 | fPreviousSftOrigin (0.,0.,0.), |
---|
63 | fPreviousMassSafety ( 0.0 ), |
---|
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 |
---|
72 | { |
---|
73 | G4TransportationManager* transportMgr ; |
---|
74 | |
---|
75 | transportMgr = G4TransportationManager::GetTransportationManager() ; |
---|
76 | |
---|
77 | fMassNavigator = transportMgr->GetNavigatorForTracking() ; |
---|
78 | fFieldPropagator = transportMgr->GetPropagatorInField() ; |
---|
79 | // fGlobalFieldMgr = transportMgr->GetFieldManager() ; |
---|
80 | fNavigatorId= transportMgr->ActivateNavigator( fMassNavigator ); |
---|
81 | if( fVerboseLevel > 0 ){ |
---|
82 | G4cout << " G4CoupledTransportation constructor: ----- " << G4endl; |
---|
83 | G4cout << " Verbose level is " << fVerboseLevel << G4endl; |
---|
84 | G4cout << " Navigator Id obtained in G4CoupledTransportation constructor " |
---|
85 | << fNavigatorId << G4endl; |
---|
86 | } |
---|
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; |
---|
97 | } |
---|
98 | |
---|
99 | ////////////////////////////////////////////////////////////////////////// |
---|
100 | |
---|
101 | G4CoupledTransportation::~G4CoupledTransportation() |
---|
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 | } |
---|
110 | } |
---|
111 | |
---|
112 | ////////////////////////////////////////////////////////////////////////// |
---|
113 | // |
---|
114 | // Responsibilities: |
---|
115 | // Find whether the geometry limits the Step, and to what length |
---|
116 | // Calculate the new value of the safety and return it. |
---|
117 | // Store the final time, position and momentum. |
---|
118 | |
---|
119 | G4double G4CoupledTransportation:: |
---|
120 | AlongStepGetPhysicalInteractionLength( const G4Track& track, |
---|
121 | G4double, // previousStepSize |
---|
122 | G4double currentMinimumStep, |
---|
123 | G4double& proposedSafetyForStart, |
---|
124 | G4GPILSelection* selection ) |
---|
125 | { |
---|
126 | G4double geometryStepLength; |
---|
127 | G4double startMassSafety= 0.0; // estimated safety for start point (mass geometry) |
---|
128 | G4double startFullSafety= 0.0; // estimated safety for start point (all geometries) |
---|
129 | G4double safetyProposal= -1.0; // local copy of proposal |
---|
130 | |
---|
131 | G4ThreeVector EndUnitMomentum ; |
---|
132 | G4double lengthAlongCurve=0.0 ; |
---|
133 | |
---|
134 | fParticleIsLooping = false ; |
---|
135 | |
---|
136 | // Initial actions moved to StartTrack() |
---|
137 | // -------------------------------------- |
---|
138 | // Note: in case another process changes touchable handle |
---|
139 | // it will be necessary to add here (for all steps) |
---|
140 | // fCurrentTouchableHandle = aTrack->GetTouchableHandle(); |
---|
141 | |
---|
142 | // GPILSelection is set to defaule value of CandidateForSelection |
---|
143 | // It is a return value |
---|
144 | // |
---|
145 | *selection = CandidateForSelection ; |
---|
146 | |
---|
147 | // Get initial Energy/Momentum of the track |
---|
148 | // |
---|
149 | const G4DynamicParticle* pParticle = track.GetDynamicParticle() ; |
---|
150 | const G4ParticleDefinition* pParticleDef = pParticle->GetDefinition() ; |
---|
151 | G4ThreeVector startMomentumDir = pParticle->GetMomentumDirection() ; |
---|
152 | G4ThreeVector startPosition = track.GetPosition() ; |
---|
153 | G4VPhysicalVolume* currentVolume= track.GetVolume(); |
---|
154 | |
---|
155 | #ifdef G4DEBUG_TRANSPORT |
---|
156 | if( fVerboseLevel > 1 ) { |
---|
157 | G4cout << "G4CoupledTransportation::AlongStepGPIL> called in volume " |
---|
158 | << currentVolume->GetName() << G4endl; |
---|
159 | } |
---|
160 | #endif |
---|
161 | // G4double theTime = track.GetGlobalTime() ; |
---|
162 | |
---|
163 | // The Step Point safety can be limited by other geometries and/or the |
---|
164 | // assumptions of any process - it's not always the geometrical safety. |
---|
165 | // We calculate the starting point's isotropic safety here. |
---|
166 | // |
---|
167 | G4ThreeVector OriginShift = startPosition - fPreviousSftOrigin ; |
---|
168 | G4double MagSqShift = OriginShift.mag2() ; |
---|
169 | startMassSafety = 0.0; |
---|
170 | startFullSafety= 0.0; |
---|
171 | |
---|
172 | // Recall that FullSafety <= MassSafety |
---|
173 | // Original: if( MagSqShift < sqr(fPreviousMassSafety) ) { |
---|
174 | if( MagSqShift < sqr(fPreviousFullSafety) ) { // Revision proposed by Alex H, 2 Oct 07 |
---|
175 | G4double mag_shift= std::sqrt(MagSqShift); |
---|
176 | startMassSafety = std::max( (fPreviousMassSafety - mag_shift), 0.0); |
---|
177 | startFullSafety = std::max( (fPreviousFullSafety - mag_shift), 0.0); |
---|
178 | // Need to be consistent between full safety with Mass safety |
---|
179 | // in order reproduce results in simple case --> use same calculation method |
---|
180 | |
---|
181 | // Only compute full safety if massSafety > 0. Else it remains 0 |
---|
182 | // startFullSafety = fPathFinder->ComputeSafety( startPosition ); |
---|
183 | } |
---|
184 | |
---|
185 | // Is the particle charged ? |
---|
186 | // |
---|
187 | G4double particleCharge = pParticle->GetCharge() ; |
---|
188 | |
---|
189 | fMassGeometryLimitedStep = false ; // Set default - alex |
---|
190 | fAnyGeometryLimitedStep = false; |
---|
191 | |
---|
192 | // fEndGlobalTimeComputed = false ; |
---|
193 | |
---|
194 | // There is no need to locate the current volume. It is Done elsewhere: |
---|
195 | // On track construction |
---|
196 | // By the tracking, after all AlongStepDoIts, in "Relocation" |
---|
197 | // Check whether the particle has an (EM) field force exerted upon it |
---|
198 | // |
---|
199 | G4FieldManager* fieldMgr=0; |
---|
200 | G4bool fieldExertsForce = false; |
---|
201 | if( (particleCharge != 0.0) ) // || (magneticMoment != 0.0 ) ) |
---|
202 | { |
---|
203 | fieldMgr= fFieldPropagator->FindAndSetFieldManager( currentVolume ); |
---|
204 | if (fieldMgr != 0) { |
---|
205 | // Message the field Manager, to configure it for this track |
---|
206 | fieldMgr->ConfigureForTrack( &track ); |
---|
207 | fieldExertsForce = (fieldMgr->GetDetectorField() != 0); |
---|
208 | } |
---|
209 | // the PathFinder will recognise whether the field exerts force |
---|
210 | } |
---|
211 | G4double momentumMagnitude = pParticle->GetTotalMomentum() ; |
---|
212 | G4double restMass = pParticleDef->GetPDGMass() ; |
---|
213 | |
---|
214 | fFieldPropagator->SetChargeMomentumMass( particleCharge, // in e+ units |
---|
215 | momentumMagnitude, // in Mev/c |
---|
216 | restMass ) ; |
---|
217 | // This should be obsolete - the call to SetChargeAndMoments below should do the work |
---|
218 | |
---|
219 | G4ThreeVector spin = track.GetPolarization() ; |
---|
220 | G4FieldTrack theFieldTrack = G4FieldTrack( startPosition, |
---|
221 | track.GetMomentumDirection(), |
---|
222 | 0.0, |
---|
223 | track.GetKineticEnergy(), |
---|
224 | restMass, |
---|
225 | 0.0, // UNUSED: track.GetVelocity(), |
---|
226 | track.GetGlobalTime(), // Lab. |
---|
227 | track.GetProperTime(), // Part. |
---|
228 | &spin ) ; |
---|
229 | theFieldTrack.SetChargeAndMoments( particleCharge ); // EM moments -- future extension |
---|
230 | |
---|
231 | G4int stepNo= track.GetCurrentStepNumber(); |
---|
232 | |
---|
233 | ELimited limitedStep; |
---|
234 | G4FieldTrack endTrackState('a'); // Default values |
---|
235 | |
---|
236 | fMassGeometryLimitedStep = false ; // default |
---|
237 | fAnyGeometryLimitedStep = false ; |
---|
238 | if( currentMinimumStep > 0 ) { |
---|
239 | G4double newMassSafety= 0.0; // temp. for recalculation |
---|
240 | |
---|
241 | // Do the Transport in the field (non recti-linear) |
---|
242 | // |
---|
243 | lengthAlongCurve = fPathFinder->ComputeStep( theFieldTrack, |
---|
244 | currentMinimumStep, |
---|
245 | fNavigatorId, |
---|
246 | stepNo, |
---|
247 | newMassSafety, |
---|
248 | limitedStep, |
---|
249 | endTrackState, |
---|
250 | currentVolume ) ; |
---|
251 | // G4cout << " PathFinder ComputeStep returns " << lengthAlongCurve << G4endl; |
---|
252 | |
---|
253 | G4double newFullSafety= fPathFinder->GetCurrentSafety(); |
---|
254 | // this was estimated already in step above |
---|
255 | // G4double newFullStep= fPathFinder->GetMinimumStep(); |
---|
256 | |
---|
257 | if( limitedStep == kUnique || limitedStep == kSharedTransport ) { |
---|
258 | fMassGeometryLimitedStep = true ; |
---|
259 | } |
---|
260 | |
---|
261 | fAnyGeometryLimitedStep = (fPathFinder->GetNumberGeometriesLimitingStep() != 0); |
---|
262 | |
---|
263 | // #ifdef G4DEBUG |
---|
264 | if( fMassGeometryLimitedStep && !fAnyGeometryLimitedStep ){ |
---|
265 | G4cerr << " Error in determining geometries limiting the step" << G4endl; |
---|
266 | G4cerr << " Limiting: mass=" << fMassGeometryLimitedStep |
---|
267 | << " any= " << fAnyGeometryLimitedStep << G4endl; |
---|
268 | G4Exception("G4CoupledTransportation::AlongStepGetPhysicalInteractionLength()", |
---|
269 | "PathFinderConfused", |
---|
270 | FatalException, |
---|
271 | "Incompatible conditions - was limited by a geometry?"); |
---|
272 | } |
---|
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!) |
---|
280 | |
---|
281 | geometryStepLength = std::min( lengthAlongCurve, currentMinimumStep); |
---|
282 | |
---|
283 | // Momentum: Magnitude and direction can be changed too now ... |
---|
284 | // |
---|
285 | fMomentumChanged = true ; |
---|
286 | fTransportEndMomentumDir = endTrackState.GetMomentumDir() ; |
---|
287 | |
---|
288 | // Remember last safety origin & value. |
---|
289 | fPreviousSftOrigin = startPosition ; |
---|
290 | fPreviousMassSafety = newMassSafety ; |
---|
291 | fPreviousFullSafety = newFullSafety ; |
---|
292 | // fpSafetyHelper->SetCurrentSafety( newFullSafety, startPosition); |
---|
293 | |
---|
294 | #ifdef G4DEBUG_TRANSPORT |
---|
295 | if( fVerboseLevel > 1 ){ |
---|
296 | G4cout << "G4Transport:CompStep> " |
---|
297 | << " called the pathfinder for a new step at " << startPosition |
---|
298 | << " and obtained step = " << lengthAlongCurve << G4endl; |
---|
299 | G4cout << " New safety (preStep) = " << newMassSafety |
---|
300 | << " versus precalculated = " << startMassSafety << G4endl; |
---|
301 | } |
---|
302 | #endif |
---|
303 | |
---|
304 | // Store as best estimate value |
---|
305 | startMassSafety = newMassSafety ; |
---|
306 | startFullSafety = newFullSafety ; |
---|
307 | |
---|
308 | // Get the End-Position and End-Momentum (Dir-ection) |
---|
309 | fTransportEndPosition = endTrackState.GetPosition() ; |
---|
310 | fTransportEndKineticEnergy = endTrackState.GetKineticEnergy() ; |
---|
311 | } else { |
---|
312 | geometryStepLength = lengthAlongCurve= 0.0 ; |
---|
313 | fMomentumChanged = false ; |
---|
314 | // fMassGeometryLimitedStep = false ; // --- ??? |
---|
315 | // fAnyGeometryLimitedStep = true; |
---|
316 | fTransportEndMomentumDir = track.GetMomentumDirection(); |
---|
317 | fTransportEndKineticEnergy = track.GetKineticEnergy(); |
---|
318 | |
---|
319 | fTransportEndPosition = startPosition; |
---|
320 | // If the step length requested is 0, and we are on a boundary |
---|
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 |
---|
327 | } |
---|
328 | // G4FieldTrack aTrackState(endTrackState); |
---|
329 | |
---|
330 | if( !fieldExertsForce ) |
---|
331 | { |
---|
332 | fParticleIsLooping = false ; |
---|
333 | fMomentumChanged = false ; |
---|
334 | fEndGlobalTimeComputed = false ; |
---|
335 | // G4cout << " global time is false " << G4endl; |
---|
336 | } |
---|
337 | else |
---|
338 | { |
---|
339 | |
---|
340 | #ifdef G4DEBUG_TRANSPORT |
---|
341 | if( fVerboseLevel > 1 ){ |
---|
342 | G4cout << " G4CT::CS End Position = " << fTransportEndPosition << G4endl; |
---|
343 | G4cout << " G4CT::CS End Direction = " << fTransportEndMomentumDir << G4endl; |
---|
344 | } |
---|
345 | #endif |
---|
346 | // G4cout << " G4CoupledTransportation Before if change energy statement: " << fFieldPropagator->GetCurrentFieldManager()->DoesFieldChangeEnergy() << G4endl; |
---|
347 | if( fFieldPropagator->GetCurrentFieldManager()->DoesFieldChangeEnergy() ) |
---|
348 | { |
---|
349 | // If the field can change energy, then the time must be integrated |
---|
350 | // - so this should have been updated |
---|
351 | // |
---|
352 | fCandidateEndGlobalTime = endTrackState.GetLabTimeOfFlight(); |
---|
353 | fEndGlobalTimeComputed = true; |
---|
354 | // G4cout << " setting global time to true - why? " << G4endl; |
---|
355 | |
---|
356 | // was ( fCandidateEndGlobalTime != track.GetGlobalTime() ); |
---|
357 | // a cleaner way is to have FieldTrack knowing whether time is updated. |
---|
358 | } |
---|
359 | else |
---|
360 | { |
---|
361 | // The energy should be unchanged by field transport, |
---|
362 | // - so the time changed will be calculated elsewhere |
---|
363 | // |
---|
364 | fEndGlobalTimeComputed = false; |
---|
365 | |
---|
366 | // Check that the integration preserved the energy |
---|
367 | // - and if not correct this! |
---|
368 | G4double startEnergy= track.GetKineticEnergy(); |
---|
369 | G4double endEnergy= fTransportEndKineticEnergy; |
---|
370 | |
---|
371 | static G4int no_inexact_steps=0; // , no_large_ediff; |
---|
372 | G4double absEdiff = std::fabs(startEnergy- endEnergy); |
---|
373 | if( absEdiff > perMillion * endEnergy ) { |
---|
374 | no_inexact_steps++; |
---|
375 | // Possible statistics keeping here ... |
---|
376 | } |
---|
377 | #ifdef G4VERBOSE |
---|
378 | if( (fVerboseLevel > 1) && ( absEdiff > perThousand * endEnergy) ){ |
---|
379 | ReportInexactEnergy(startEnergy, endEnergy); |
---|
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(); |
---|
387 | } |
---|
388 | } |
---|
389 | |
---|
390 | endpointDistance = (fTransportEndPosition - startPosition).mag() ; |
---|
391 | fParticleIsLooping = fFieldPropagator->IsParticleLooping() ; |
---|
392 | |
---|
393 | fTransportEndSpin = endTrackState.GetSpin(); |
---|
394 | |
---|
395 | // Calculate the safety |
---|
396 | safetyProposal= startFullSafety; // used to be startMassSafety |
---|
397 | // Changed to accomodate processes that cannot update the safety -- JA 22 Nov 06 |
---|
398 | |
---|
399 | // Update safety for the end-point, if becomes negative at the end-point. |
---|
400 | |
---|
401 | if( (startFullSafety < endpointDistance ) |
---|
402 | && ( particleCharge != 0.0 ) ) // Only needed to prepare for Mult Scat. |
---|
403 | // && !fAnyGeometryLimitedStep ) // To-Try: No safety update if at a boundary |
---|
404 | { |
---|
405 | G4double endFullSafety = |
---|
406 | fPathFinder->ComputeSafety( fTransportEndPosition); |
---|
407 | // Expected mission -- only mass geometry's safety |
---|
408 | // fMassNavigator->ComputeSafety( fTransportEndPosition) ; |
---|
409 | // Yet discrete processes only have poststep -- and this cannot |
---|
410 | // currently revise the safety |
---|
411 | // ==> so we use the all-geometry safety as a precaution |
---|
412 | |
---|
413 | fpSafetyHelper->SetCurrentSafety( endFullSafety, fTransportEndPosition); |
---|
414 | // Pushing safety to Helper avoids recalculation at this point |
---|
415 | |
---|
416 | G4ThreeVector centerPt= G4ThreeVector(0.0, 0.0, 0.0); // Used for return value |
---|
417 | G4double endMassSafety= fPathFinder->ObtainSafety( fNavigatorId, centerPt); |
---|
418 | // Retrieves the mass value from PathFinder (it calculated it) |
---|
419 | |
---|
420 | fPreviousMassSafety = endMassSafety ; |
---|
421 | fPreviousFullSafety = endFullSafety; |
---|
422 | fPreviousSftOrigin = fTransportEndPosition ; |
---|
423 | |
---|
424 | // The convention (Stepping Manager's) is safety from the start point |
---|
425 | // |
---|
426 | safetyProposal = endFullSafety + endpointDistance; |
---|
427 | // --> was endMassSafety |
---|
428 | // Changed to accomodate processes that cannot update the safety -- JA 22 Nov 06 |
---|
429 | |
---|
430 | // #define G4DEBUG_TRANSPORT 1 |
---|
431 | |
---|
432 | #ifdef G4DEBUG_TRANSPORT |
---|
433 | int prec= G4cout.precision(12) ; |
---|
434 | G4cout << "***Transportation::AlongStepGPIL ** " << G4endl ; |
---|
435 | G4cout << " Revised Safety at endpoint " << fTransportEndPosition |
---|
436 | << " give safety values: Mass= " << endMassSafety |
---|
437 | << " All= " << endFullSafety << G4endl ; |
---|
438 | G4cout << " Adding endpoint distance " << endpointDistance |
---|
439 | << " to obtain pseudo-safety= " << safetyProposal << G4endl ; |
---|
440 | G4cout.precision(prec); |
---|
441 | } |
---|
442 | else{ |
---|
443 | int prec= G4cout.precision(12) ; |
---|
444 | G4cout << "***Transportation::AlongStepGPIL ** " << G4endl ; |
---|
445 | G4cout << " Quick Safety estimate at endpoint " << fTransportEndPosition |
---|
446 | << " gives safety endpoint value = " << startFullSafety - endpointDistance |
---|
447 | << " using start-point value " << startFullSafety |
---|
448 | << " and endpointDistance " << endpointDistance << G4endl; |
---|
449 | G4cout.precision(prec); |
---|
450 | #endif |
---|
451 | } |
---|
452 | |
---|
453 | proposedSafetyForStart= safetyProposal; |
---|
454 | fParticleChange.ProposeTrueStepLength(geometryStepLength) ; |
---|
455 | |
---|
456 | return geometryStepLength ; |
---|
457 | } |
---|
458 | |
---|
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 | } |
---|
609 | |
---|
610 | ////////////////////////////////////////////////////////////////////////// |
---|
611 | // |
---|
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 ) |
---|
628 | { |
---|
629 | G4ThreeVector moveVec = ( NewVector - OldVector ); |
---|
630 | |
---|
631 | G4cerr << G4endl |
---|
632 | << "**************************************************************" << G4endl; |
---|
633 | G4cerr << "Endpoint has moved between value expected from TransportEndPosition " |
---|
634 | << " and value from Track in PostStepDoIt. " << G4endl |
---|
635 | << "Change of " << Quantity << " is " << moveVec.mag() / mm << " mm long, " |
---|
636 | << " and its vector is " << (1.0/mm) * moveVec << " mm " << G4endl |
---|
637 | << "Endpoint of ComputeStep was " << OldVector |
---|
638 | << " and current position to locate is " << NewVector << G4endl; |
---|
639 | } |
---|
640 | |
---|
641 | ///////////////////////////////////////////////////////////////////////////// |
---|
642 | |
---|
643 | G4VParticleChange* G4CoupledTransportation::PostStepDoIt( const G4Track& track, |
---|
644 | const G4Step& ) |
---|
645 | { |
---|
646 | G4TouchableHandle retCurrentTouchable ; // The one to return |
---|
647 | |
---|
648 | // Initialize ParticleChange (by setting all its members equal |
---|
649 | // to corresponding members in G4Track) |
---|
650 | // fParticleChange.Initialize(track) ; // To initialise TouchableChange |
---|
651 | |
---|
652 | fParticleChange.ProposeTrackStatus(track.GetTrackStatus()) ; |
---|
653 | |
---|
654 | // Check that the end position and direction are preserved |
---|
655 | // since call to AlongStepDoIt |
---|
656 | if( (fTransportEndPosition - track.GetPosition()).mag2() >= 1.0e-16 ){ |
---|
657 | static G4String EndLabelString("End of Step Position"); |
---|
658 | ReportMove( track.GetPosition(), fTransportEndPosition, EndLabelString ); |
---|
659 | G4cerr << " Problem in G4CoupledTransportation::PostStepDoIt " << G4endl; |
---|
660 | } |
---|
661 | |
---|
662 | // If the Step was determined by the volume boundary, relocate the particle |
---|
663 | // The pathFinder will know that the geometry limited the step (!?) |
---|
664 | |
---|
665 | if( fVerboseLevel > 0 ){ |
---|
666 | G4cout << " Calling PathFinder::Locate() from " |
---|
667 | << " G4CoupledTransportation::PostStepDoIt " << G4endl; |
---|
668 | G4cout << " fAnyGeometryLimitedStep is " << fAnyGeometryLimitedStep << G4endl; |
---|
669 | |
---|
670 | } |
---|
671 | if(fAnyGeometryLimitedStep) |
---|
672 | { |
---|
673 | fPathFinder->Locate( track.GetPosition(), |
---|
674 | track.GetMomentumDirection(), |
---|
675 | true); |
---|
676 | |
---|
677 | // fCurrentTouchable will now become the previous touchable, |
---|
678 | // and what was the previous will be freed. |
---|
679 | // (Needed because the preStepPoint can point to the previous touchable) |
---|
680 | if( fVerboseLevel > 0 ) |
---|
681 | G4cout << "G4CoupledTransportation::PostStepDoIt --- fNavigatorId = " |
---|
682 | << fNavigatorId << G4endl; |
---|
683 | |
---|
684 | fCurrentTouchableHandle= |
---|
685 | fPathFinder->CreateTouchableHandle( fNavigatorId ); |
---|
686 | |
---|
687 | // Check whether the particle is out of the world volume |
---|
688 | // If so it has exited and must be killed. |
---|
689 | // |
---|
690 | #ifdef G4DEBUG_TRANSPORT |
---|
691 | if( fVerboseLevel > 1 ){ |
---|
692 | G4VPhysicalVolume* vol= fCurrentTouchableHandle->GetVolume(); |
---|
693 | G4cout << "CHECK !!!!!!!!!!! fCurrentTouchableHandle->GetVolume() = " << vol; |
---|
694 | if( vol ) { G4cout << "Name=" << vol->GetName(); } |
---|
695 | G4cout << G4endl; |
---|
696 | } |
---|
697 | #endif |
---|
698 | if( fCurrentTouchableHandle->GetVolume() == 0 ) |
---|
699 | { |
---|
700 | fParticleChange.ProposeTrackStatus( fStopAndKill ) ; |
---|
701 | } |
---|
702 | retCurrentTouchable = fCurrentTouchableHandle ; |
---|
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 | |
---|
710 | } |
---|
711 | else // fAnyGeometryLimitedStep is false |
---|
712 | { |
---|
713 | #ifdef G4DEBUG_TRANSPORT |
---|
714 | if( fVerboseLevel > 1 ){ |
---|
715 | G4cout << "G4CoupledTransportation::PostStepDoIt -- " |
---|
716 | << " fAnyGeometryLimitedStep = " << fAnyGeometryLimitedStep |
---|
717 | << " must be false " << G4endl; |
---|
718 | } |
---|
719 | #endif |
---|
720 | // This serves only to move each of the Navigator's location |
---|
721 | // |
---|
722 | // fLinearNavigator->LocateGlobalPointWithinVolume( track.GetPosition() ) ; |
---|
723 | |
---|
724 | // G4cout << "G4CoupledTransportation calling PathFinder::ReLocate() " << G4endl; |
---|
725 | fPathFinder->ReLocate( track.GetPosition() ); |
---|
726 | // track.GetMomentumDirection() ); |
---|
727 | |
---|
728 | // Keep the value of the track's current Touchable is retained, |
---|
729 | // and use it to overwrite the (unset) one in particle change. |
---|
730 | // Expect this must be fCurrentTouchable too |
---|
731 | // - could it be different, eg at the start of a step ? |
---|
732 | // |
---|
733 | retCurrentTouchable = track.GetTouchableHandle() ; |
---|
734 | // fParticleChange.SetTouchableHandle( track.GetTouchableHandle() ) ; |
---|
735 | |
---|
736 | // Have not reached a boundary |
---|
737 | fParticleChange.ProposeLastStepInVolume(false); |
---|
738 | } // endif ( fAnyGeometryLimitedStep ) |
---|
739 | |
---|
740 | const G4VPhysicalVolume* pNewVol = retCurrentTouchable->GetVolume() ; |
---|
741 | const G4Material* pNewMaterial = 0 ; |
---|
742 | const G4VSensitiveDetector* pNewSensitiveDetector = 0 ; |
---|
743 | |
---|
744 | if( pNewVol != 0 ) |
---|
745 | { |
---|
746 | pNewMaterial= pNewVol->GetLogicalVolume()->GetMaterial(); |
---|
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) ; |
---|
777 | |
---|
778 | return &fParticleChange ; |
---|
779 | } |
---|
780 | |
---|
781 | // New method takes over the responsibility to reset the state of |
---|
782 | // G4CoupledTransportation object: |
---|
783 | // - at the start of a new track, and |
---|
784 | // - on the resumption of a suspended track. |
---|
785 | |
---|
786 | void |
---|
787 | G4CoupledTransportation::StartTracking(G4Track* aTrack) |
---|
788 | { |
---|
789 | |
---|
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 | // } |
---|
806 | G4ThreeVector position = aTrack->GetPosition(); |
---|
807 | G4ThreeVector direction = aTrack->GetMomentumDirection(); |
---|
808 | |
---|
809 | // if( fVerboseLevel > 1 ){ |
---|
810 | // G4cout << " Calling PathFinder::PrepareNewTrack from " |
---|
811 | // << " G4CoupledTransportation::StartTracking -- which calls Locate()" << G4endl; |
---|
812 | // } |
---|
813 | fPathFinder->PrepareNewTrack( position, direction); |
---|
814 | // This implies a call to fPathFinder->Locate( position, direction ); |
---|
815 | |
---|
816 | // Global field, if any, must exist before tracking is started |
---|
817 | fGlobalFieldExists= DoesGlobalFieldExist(); |
---|
818 | // reset safety value and center |
---|
819 | // |
---|
820 | fPreviousMassSafety = 0.0 ; |
---|
821 | fPreviousFullSafety = 0.0 ; |
---|
822 | fPreviousSftOrigin = G4ThreeVector(0.,0.,0.) ; |
---|
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 ) { .. } |
---|
829 | |
---|
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(); |
---|
842 | |
---|
843 | #ifdef G4DEBUG_TRANSPORT |
---|
844 | if( fVerboseLevel > 1 ){ |
---|
845 | G4cout << " Returning touchable handle " << fCurrentTouchableHandle << G4endl; |
---|
846 | } |
---|
847 | #endif |
---|
848 | |
---|
849 | // Update the current touchable handle (from the track's) |
---|
850 | // |
---|
851 | fCurrentTouchableHandle = aTrack->GetTouchableHandle(); |
---|
852 | } |
---|
853 | |
---|
854 | void |
---|
855 | G4CoupledTransportation:: |
---|
856 | ReportInexactEnergy(G4double startEnergy, G4double endEnergy) |
---|
857 | { |
---|
858 | static G4int no_warnings= 0, warnModulo=1, moduloFactor= 10, no_large_ediff= 0; |
---|
859 | |
---|
860 | if( std::fabs(startEnergy- endEnergy) > perThousand * endEnergy ) |
---|
861 | { |
---|
862 | no_large_ediff ++; |
---|
863 | if( (no_large_ediff% warnModulo) == 0 ) |
---|
864 | { |
---|
865 | no_warnings++; |
---|
866 | G4cout << "WARNING - G4CoupledTransportation::AlongStepGetPIL() " |
---|
867 | << " Energy change in Step is above 1^-3 relative value. " << G4endl |
---|
868 | << " Relative change in 'tracking' step = " |
---|
869 | << std::setw(15) << (endEnergy-startEnergy)/startEnergy << G4endl |
---|
870 | << " Starting E= " << std::setw(12) << startEnergy / MeV << " MeV " << G4endl |
---|
871 | << " Ending E= " << std::setw(12) << endEnergy / MeV << " MeV " << G4endl; |
---|
872 | G4cout << " Energy has been corrected -- however, review" |
---|
873 | << " field propagation parameters for accuracy." << G4endl; |
---|
874 | if( (fVerboseLevel > 2 ) || (no_warnings<4) || (no_large_ediff == warnModulo * moduloFactor) ){ |
---|
875 | G4cout << " These include EpsilonStepMax(/Min) in G4FieldManager " |
---|
876 | << " which determine fractional error per step for integrated quantities. " << G4endl |
---|
877 | << " Note also the influence of the permitted number of integration steps." |
---|
878 | << G4endl; |
---|
879 | } |
---|
880 | G4cerr << "ERROR - G4CoupledTransportation::AlongStepGetPIL()" << G4endl |
---|
881 | << " Bad 'endpoint'. Energy change detected" |
---|
882 | << " and corrected. " |
---|
883 | << " Has occurred already " |
---|
884 | << no_large_ediff << " times." << G4endl; |
---|
885 | if( no_large_ediff == warnModulo * moduloFactor ) |
---|
886 | { |
---|
887 | warnModulo *= moduloFactor; |
---|
888 | } |
---|
889 | } |
---|
890 | } |
---|
891 | } |
---|