Bug Summary

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