Bug Summary

File:scratch/gluex/scan-build-work/hdgeant4/src/G4fixes/G4CoupledTransportation.cc
Location:line 393, column 7
Description:Value stored to 'startMassSafety' is never read

Annotated Source Code

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