Bug Summary

File:scratch/gluex/scan-build-work/hdgeant4^hdr4251/src/G4fixes/G4Transportation.cc
Location:line 817, column 3
Description:Called C++ object pointer is null

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: G4Transportation.cc 2011/06/10 16:19:46 japost Exp japost $
28//
29// ------------------------------------------------------------
30// GEANT 4 include file implementation
31//
32// ------------------------------------------------------------
33//
34// This class is a process responsible for the transportation of
35// a particle, ie the geometrical propagation that encounters the
36// geometrical sub-volumes of the detectors.
37//
38// It is also tasked with the key role of proposing the "isotropic safety",
39// which will be used to update the post-step point's safety.
40//
41// =======================================================================
42// Modified:
43// 10 Jan 2015, M.Kelsey: Use G4DynamicParticle mass, NOT PDGMass
44// 28 Oct 2011, P.Gumpl./J.Ap: Detect gravity field, use magnetic moment
45// 20 Nov 2008, J.Apostolakis: Push safety to helper - after ComputeSafety
46// 9 Nov 2007, J.Apostolakis: Flag for short steps, push safety to helper
47// 19 Jan 2006, P.MoraDeFreitas: Fix for suspended tracks (StartTracking)
48// 11 Aug 2004, M.Asai: Add G4VSensitiveDetector* for updating stepPoint.
49// 21 June 2003, J.Apostolakis: Calling field manager with
50// track, to enable it to configure its accuracy
51// 13 May 2003, J.Apostolakis: Zero field areas now taken into
52// account correclty in all cases (thanks to W Pokorski).
53// 29 June 2001, J.Apostolakis, D.Cote-Ahern, P.Gumplinger:
54// correction for spin tracking
55// 20 Febr 2001, J.Apostolakis: update for new FieldTrack
56// 22 Sept 2000, V.Grichine: update of Kinetic Energy
57// Created: 19 March 1997, J. Apostolakis
58// =======================================================================
59
60#include "G4Transportation.hh"
61#include "G4TransportationProcessType.hh"
62
63#include "G4PhysicalConstants.hh"
64#include "G4SystemOfUnits.hh"
65#include "G4ProductionCutsTable.hh"
66#include "G4ParticleTable.hh"
67
68#include "G4ChargeState.hh"
69#include "G4EquationOfMotion.hh"
70
71#include "G4FieldManagerStore.hh"
72
73class G4VSensitiveDetector;
74
75G4bool G4Transportation::fUseMagneticMoment=false;
76
77// #define G4DEBUG_TRANSPORT 1
78
79//////////////////////////////////////////////////////////////////////////
80//
81// Constructor
82
83G4Transportation::G4Transportation( G4int verbosity )
84 : G4VProcess( G4String("Transportation"), fTransportation ),
85 fTransportEndPosition( 0.0, 0.0, 0.0 ),
86 fTransportEndMomentumDir( 0.0, 0.0, 0.0 ),
87 fTransportEndKineticEnergy( 0.0 ),
88 fTransportEndSpin( 0.0, 0.0, 0.0 ),
89 fMomentumChanged(true),
90 fEndGlobalTimeComputed(false),
91 fCandidateEndGlobalTime(0.0),
92 fParticleIsLooping( false ),
93 fNewTrack( true ),
94 fFirstStepInVolume( true ),
95 fLastStepInVolume( false ),
96 fGeometryLimitedStep(true),
97 fFieldExertedForce( false ),
98 fPreviousSftOrigin( 0.,0.,0. ),
99 fPreviousSafety( 0.0 ),
100 // fParticleChange(),
101 fEndPointDistance( -1.0 ),
102 fThreshold_Warning_Energy( 100 * MeV ),
103 fThreshold_Important_Energy( 250 * MeV ),
104 fThresholdTrials( 10 ),
105 fNoLooperTrials( 0 ),
106 fSumEnergyKilled( 0.0 ), fMaxEnergyKilled( 0.0 ),
107 fShortStepOptimisation( false ), // Old default: true (=fast short steps)
108 fVerboseLevel( verbosity )
109{
110 // set Process Sub Type
111 SetProcessSubType(static_cast<G4int>(TRANSPORTATION));
112 pParticleChange= &fParticleChange; // Required to conform to G4VProcess
113
114 G4TransportationManager* transportMgr ;
115
116 transportMgr = G4TransportationManager::GetTransportationManager() ;
117
118 fLinearNavigator = transportMgr->GetNavigatorForTracking() ;
119
120 fFieldPropagator = transportMgr->GetPropagatorInField() ;
121
122 fpSafetyHelper = transportMgr->GetSafetyHelper(); // New
123
124 // Cannot determine whether a field exists here, as it would
125 // depend on the relative order of creating the detector's
126 // field and this process. That order is not guaranted.
127 // Instead later the method DoesGlobalFieldExist() is called
128
129 static G4ThreadLocalthread_local G4TouchableHandle* pNullTouchableHandle = 0;
130 if ( !pNullTouchableHandle) { pNullTouchableHandle = new G4TouchableHandle; }
131 fCurrentTouchableHandle = *pNullTouchableHandle;
132 // Points to (G4VTouchable*) 0
133
134#ifdef G4VERBOSE1
135 if( fVerboseLevel > 0)
136 {
137 G4cout(*G4cout_p) << " G4Transportation constructor> set fShortStepOptimisation to ";
138 if ( fShortStepOptimisation ) G4cout(*G4cout_p) << "true" << G4endlstd::endl;
139 else G4cout(*G4cout_p) << "false" << G4endlstd::endl;
140 }
141#endif
142}
143
144//////////////////////////////////////////////////////////////////////////
145
146G4Transportation::~G4Transportation()
147{
148 if( (fVerboseLevel > 0) && (fSumEnergyKilled > 0.0 ) )
149 {
150 G4cout(*G4cout_p) << " G4Transportation: Statistics for looping particles " << G4endlstd::endl;
151 G4cout(*G4cout_p) << " Sum of energy of loopers killed: " << fSumEnergyKilled << G4endlstd::endl;
152 G4cout(*G4cout_p) << " Max energy of loopers killed: " << fMaxEnergyKilled << G4endlstd::endl;
153 }
154}
155
156//////////////////////////////////////////////////////////////////////////
157//
158// Responsibilities:
159// Find whether the geometry limits the Step, and to what length
160// Calculate the new value of the safety and return it.
161// Store the final time, position and momentum.
162
163G4double G4Transportation::
164AlongStepGetPhysicalInteractionLength( const G4Track& track,
165 G4double, // previousStepSize
166 G4double currentMinimumStep,
167 G4double& currentSafety,
168 G4GPILSelection* selection )
169{
170 G4double geometryStepLength= -1.0, newSafety= -1.0;
171 fParticleIsLooping = false ;
172
173 // Initial actions moved to StartTrack()
174 // --------------------------------------
175 // Note: in case another process changes touchable handle
176 // it will be necessary to add here (for all steps)
177 // fCurrentTouchableHandle = aTrack->GetTouchableHandle();
178
179 // GPILSelection is set to defaule value of CandidateForSelection
180 // It is a return value
181 //
182 *selection = CandidateForSelection ;
183
184 fFirstStepInVolume= fNewTrack || fLastStepInVolume;
185 // G4cout << " Transport::AlongStep GPIL: 1st-step= " << fFirstStepInVolume << " newTrack= " << fNewTrack << " fLastStep-in-Vol= " << fLastStepInVolume << G4endl;
186 fLastStepInVolume= false;
187 fNewTrack = false;
188
189 fParticleChange.ProposeFirstStepInVolume(fFirstStepInVolume);
190
191 // Get initial Energy/Momentum of the track
192 //
193 const G4DynamicParticle* pParticle = track.GetDynamicParticle() ;
194 const G4ParticleDefinition* pParticleDef = pParticle->GetDefinition() ;
195 G4ThreeVector startMomentumDir = pParticle->GetMomentumDirection() ;
196 G4ThreeVector startPosition = track.GetPosition() ;
197
198 // G4double theTime = track.GetGlobalTime() ;
199
200 // The Step Point safety can be limited by other geometries and/or the
201 // assumptions of any process - it's not always the geometrical safety.
202 // We calculate the starting point's isotropic safety here.
203 //
204 G4ThreeVector OriginShift = startPosition - fPreviousSftOrigin ;
205 G4double MagSqShift = OriginShift.mag2() ;
206 if( MagSqShift >= sqr(fPreviousSafety) )
207 {
208 currentSafety = 0.0 ;
209 }
210 else
211 {
212 currentSafety = fPreviousSafety - std::sqrt(MagSqShift) ;
213 }
214
215 // Is the particle charged or has it a magnetic moment?
216 //
217 G4double particleCharge = pParticle->GetCharge() ;
218 G4double magneticMoment = pParticle->GetMagneticMoment() ;
219 G4double restMass = pParticle->GetMass() ;
220
221 fGeometryLimitedStep = false ;
222 // fEndGlobalTimeComputed = false ;
223
224 // There is no need to locate the current volume. It is Done elsewhere:
225 // On track construction
226 // By the tracking, after all AlongStepDoIts, in "Relocation"
227
228 // Check if the particle has a force, EM or gravitational, exerted on it
229 //
230 G4FieldManager* fieldMgr=0;
231 G4bool fieldExertsForce = false ;
232
233 G4bool gravityOn = false;
234 G4bool fieldExists= false; // Field is not 0 (null pointer)
235
236 fieldMgr = fFieldPropagator->FindAndSetFieldManager( track.GetVolume() );
237 if( fieldMgr != 0 )
238 {
239 // Message the field Manager, to configure it for this track
240 fieldMgr->ConfigureForTrack( &track );
241 // Is here to allow a transition from no-field pointer
242 // to finite field (non-zero pointer).
243
244 // If the field manager has no field ptr, the field is zero
245 // by definition ( = there is no field ! )
246 const G4Field* ptrField= fieldMgr->GetDetectorField();
247 fieldExists = (ptrField!=0) ;
248 if( fieldExists )
249 {
250 gravityOn= ptrField->IsGravityActive();
251
252 if( (particleCharge != 0.0)
253 || (fUseMagneticMoment && (magneticMoment != 0.0) )
254 || (gravityOn && (restMass != 0.0) )
255 )
256 {
257 fieldExertsForce = fieldExists;
258 }
259 }
260 }
261 // G4cout << " G4Transport: field exerts force= " << fieldExertsForce
262 // << " fieldMgr= " << fieldMgr << G4endl;
263 fFieldExertedForce = fieldExertsForce;
264
265 if( !fieldExertsForce )
266 {
267 G4double linearStepLength ;
268 if( fShortStepOptimisation && (currentMinimumStep <= currentSafety) )
269 {
270 // The Step is guaranteed to be taken
271 //
272 geometryStepLength = currentMinimumStep ;
273 fGeometryLimitedStep = false ;
274 }
275 else
276 {
277 // Find whether the straight path intersects a volume
278 //
279 linearStepLength = fLinearNavigator->ComputeStep( startPosition,
280 startMomentumDir,
281 currentMinimumStep,
282 newSafety) ;
283 // Remember last safety origin & value.
284 //
285 fPreviousSftOrigin = startPosition ;
286 fPreviousSafety = newSafety ;
287 fpSafetyHelper->SetCurrentSafety( newSafety, startPosition);
288
289 currentSafety = newSafety ;
290
291 fGeometryLimitedStep= (linearStepLength <= currentMinimumStep);
292 if( fGeometryLimitedStep )
293 {
294 // The geometry limits the Step size (an intersection was found.)
295 geometryStepLength = linearStepLength ;
296 }
297 else
298 {
299 // The full Step is taken.
300 geometryStepLength = currentMinimumStep ;
301 }
302 }
303 fEndPointDistance = geometryStepLength ;
304
305 // Calculate final position
306 //
307 fTransportEndPosition = startPosition+geometryStepLength*startMomentumDir ;
308
309 // Momentum direction, energy and polarisation are unchanged by transport
310 //
311 fTransportEndMomentumDir = startMomentumDir ;
312 fTransportEndKineticEnergy = track.GetKineticEnergy() ;
313 fTransportEndSpin = track.GetPolarization();
314 fParticleIsLooping = false ;
315 fMomentumChanged = false ;
316 fEndGlobalTimeComputed = false ;
317 }
318 else // A field exerts force
319 {
320 G4double momentumMagnitude = pParticle->GetTotalMomentum() ;
321 G4ThreeVector EndUnitMomentum ;
322 G4double lengthAlongCurve ;
323
324 G4ChargeState chargeState(particleCharge, // The charge can change (dynamic)
325 magneticMoment,
326 pParticleDef->GetPDGSpin() );
327 // For insurance, could set it again
328 // chargeState.SetPDGSpin(pParticleDef->GetPDGSpin() ); // Provisionally in same object
329
330 G4EquationOfMotion* equationOfMotion =
331 (fFieldPropagator->GetChordFinder()->GetIntegrationDriver()->GetStepper())
332 ->GetEquationOfMotion();
333
334// equationOfMotion->SetChargeMomentumMass( particleCharge,
335 equationOfMotion->SetChargeMomentumMass( chargeState,
336 momentumMagnitude,
337 restMass);
338
339 G4FieldTrack aFieldTrack = G4FieldTrack( startPosition,
340 track.GetGlobalTime(), // Lab.
341 // track.GetProperTime(), // Particle rest frame
342 track.GetMomentumDirection(),
343 track.GetKineticEnergy(),
344 restMass,
345 particleCharge,
346 track.GetPolarization(),
347 pParticleDef->GetPDGMagneticMoment(),
348 0.0, // Length along track
349 pParticleDef->GetPDGSpin()
350 ) ;
351
352 if( currentMinimumStep > 0 )
353 {
354 // Do the Transport in the field (non recti-linear)
355 //
356 lengthAlongCurve = fFieldPropagator->ComputeStep( aFieldTrack,
357 currentMinimumStep,
358 currentSafety,
359 track.GetVolume() ) ;
360
361 fGeometryLimitedStep= fFieldPropagator->IsLastStepInVolume();
362 // It is possible that step was reduced in PropagatorInField due to previous zero steps
363 // To cope with case that reduced step is taken in full, we must rely on PiF to obtain this
364 // value.
365
366 geometryStepLength = std::min( lengthAlongCurve, currentMinimumStep );
367
368 // Remember last safety origin & value.
369 //
370 fPreviousSftOrigin = startPosition ;
371 fPreviousSafety = currentSafety ;
372 fpSafetyHelper->SetCurrentSafety( currentSafety, startPosition);
373 }
374 else
375 {
376 geometryStepLength = lengthAlongCurve= 0.0 ;
377 fGeometryLimitedStep = false ;
378 }
379
380 // Get the End-Position and End-Momentum (Dir-ection)
381 //
382 fTransportEndPosition = aFieldTrack.GetPosition() ;
383
384 // Momentum: Magnitude and direction can be changed too now ...
385 //
386 fMomentumChanged = true ;
387 fTransportEndMomentumDir = aFieldTrack.GetMomentumDir() ;
388
389 fTransportEndKineticEnergy = aFieldTrack.GetKineticEnergy() ;
390
391 if( fFieldPropagator->GetCurrentFieldManager()->DoesFieldChangeEnergy() )
392 {
393 // If the field can change energy, then the time must be integrated
394 // - so this should have been updated
395 //
396 fCandidateEndGlobalTime = aFieldTrack.GetLabTimeOfFlight();
397 fEndGlobalTimeComputed = true;
398
399 // was ( fCandidateEndGlobalTime != track.GetGlobalTime() );
400 // a cleaner way is to have FieldTrack knowing whether time is updated.
401 }
402 else
403 {
404 // The energy should be unchanged by field transport,
405 // - so the time changed will be calculated elsewhere
406 //
407 fEndGlobalTimeComputed = false;
408
409 // Check that the integration preserved the energy
410 // - and if not correct this!
411 G4double startEnergy= track.GetKineticEnergy();
412 G4double endEnergy= fTransportEndKineticEnergy;
413
414 static G4ThreadLocalthread_local G4int no_inexact_steps=0, no_large_ediff;
415 G4double absEdiff = std::fabs(startEnergy- endEnergy);
416 if( absEdiff > perMillion * endEnergy )
417 {
418 no_inexact_steps++;
419 // Possible statistics keeping here ...
420 }
421 if( fVerboseLevel > 1 )
422 {
423 if( std::fabs(startEnergy- endEnergy) > perThousand * endEnergy )
424 {
425 static G4ThreadLocalthread_local G4int no_warnings= 0, warnModulo=1,
426 moduloFactor= 10;
427 no_large_ediff ++;
428 if( (no_large_ediff% warnModulo) == 0 )
429 {
430 no_warnings++;
431 G4cout(*G4cout_p) << "WARNING - G4Transportation::AlongStepGetPIL() "
432 << " Energy change in Step is above 1^-3 relative value. " << G4endlstd::endl
433 << " Relative change in 'tracking' step = "
434 << std::setw(15) << (endEnergy-startEnergy)/startEnergy << G4endlstd::endl
435 << " Starting E= " << std::setw(12) << startEnergy / MeV << " MeV " << G4endlstd::endl
436 << " Ending E= " << std::setw(12) << endEnergy / MeV << " MeV " << G4endlstd::endl;
437 G4cout(*G4cout_p) << " Energy has been corrected -- however, review"
438 << " field propagation parameters for accuracy." << G4endlstd::endl;
439 if( (fVerboseLevel > 2 ) || (no_warnings<4) || (no_large_ediff == warnModulo * moduloFactor) )
440 {
441 G4cout(*G4cout_p) << " These include EpsilonStepMax(/Min) in G4FieldManager "
442 << " which determine fractional error per step for integrated quantities. " << G4endlstd::endl
443 << " Note also the influence of the permitted number of integration steps."
444 << G4endlstd::endl;
445 }
446 G4cerr(*G4cerr_p) << "ERROR - G4Transportation::AlongStepGetPIL()" << G4endlstd::endl
447 << " Bad 'endpoint'. Energy change detected"
448 << " and corrected. "
449 << " Has occurred already "
450 << no_large_ediff << " times." << G4endlstd::endl;
451 if( no_large_ediff == warnModulo * moduloFactor )
452 {
453 warnModulo *= moduloFactor;
454 }
455 }
456 }
457 } // end of if (fVerboseLevel)
458
459 // Correct the energy for fields that conserve it
460 // This - hides the integration error
461 // - but gives a better physical answer
462 fTransportEndKineticEnergy= track.GetKineticEnergy();
463 }
464
465 fTransportEndSpin = aFieldTrack.GetSpin();
466 fParticleIsLooping = fFieldPropagator->IsParticleLooping() ;
467 fEndPointDistance = (fTransportEndPosition - startPosition).mag() ;
468 }
469
470 // If we are asked to go a step length of 0, and we are on a boundary
471 // then a boundary will also limit the step -> we must flag this.
472 //
473 if( currentMinimumStep == 0.0 )
474 {
475 if( currentSafety == 0.0 ) { fGeometryLimitedStep = true; }
476 }
477
478 // Update the safety starting from the end-point,
479 // if it will become negative at the end-point.
480 //
481 if( currentSafety < fEndPointDistance )
482 {
483 if( particleCharge != 0.0 )
484 {
485 G4double endSafety =
486 fLinearNavigator->ComputeSafety( fTransportEndPosition) ;
487 currentSafety = endSafety ;
488 fPreviousSftOrigin = fTransportEndPosition ;
489 fPreviousSafety = currentSafety ;
490 fpSafetyHelper->SetCurrentSafety( currentSafety, fTransportEndPosition);
491
492 // Because the Stepping Manager assumes it is from the start point,
493 // add the StepLength
494 //
495 currentSafety += fEndPointDistance ;
496
497#ifdef G4DEBUG_TRANSPORT
498 G4cout(*G4cout_p).precision(12) ;
499 G4cout(*G4cout_p) << "***G4Transportation::AlongStepGPIL ** " << G4endlstd::endl ;
500 G4cout(*G4cout_p) << " Called Navigator->ComputeSafety at " << fTransportEndPosition
501 << " and it returned safety= " << endSafety << G4endlstd::endl ;
502 G4cout(*G4cout_p) << " Adding endpoint distance " << fEndPointDistance
503 << " to obtain pseudo-safety= " << currentSafety << G4endlstd::endl ;
504 }
505 else
506 {
507 G4cout(*G4cout_p) << "***G4Transportation::AlongStepGPIL ** " << G4endlstd::endl ;
508 G4cout(*G4cout_p) << " Avoiding call to ComputeSafety : " << G4endlstd::endl;
509 G4cout(*G4cout_p) << " charge = " << particleCharge << G4endlstd::endl;
510 G4cout(*G4cout_p) << " mag moment = " << magneticMoment << G4endlstd::endl;
511#endif
512 }
513 }
514
515 fParticleChange.ProposeTrueStepLength(geometryStepLength) ;
516
517 return geometryStepLength ;
518}
519
520//////////////////////////////////////////////////////////////////////////
521//
522// Initialize ParticleChange (by setting all its members equal
523// to corresponding members in G4Track)
524
525G4VParticleChange* G4Transportation::AlongStepDoIt( const G4Track& track,
526 const G4Step& stepData )
527{
528 static G4ThreadLocalthread_local G4int noCalls=0;
529 noCalls++;
530
531 fParticleChange.Initialize(track) ;
532
533 // Code for specific process
534 //
535 fParticleChange.ProposePosition(fTransportEndPosition) ;
536 fParticleChange.ProposeMomentumDirection(fTransportEndMomentumDir) ;
537 fParticleChange.ProposeEnergy(fTransportEndKineticEnergy) ;
538 fParticleChange.SetMomentumChanged(fMomentumChanged) ;
539
540 fParticleChange.ProposePolarization(fTransportEndSpin);
541
542 G4double deltaTime = 0.0 ;
543
544 // Calculate Lab Time of Flight (ONLY if field Equations used it!)
545 // G4double endTime = fCandidateEndGlobalTime;
546 // G4double delta_time = endTime - startTime;
547
548 G4double startTime = track.GetGlobalTime() ;
549
550 if (!fEndGlobalTimeComputed)
551 {
552 // The time was not integrated .. make the best estimate possible
553 //
554 G4double initialVelocity = stepData.GetPreStepPoint()->GetVelocity();
555 G4double stepLength = track.GetStepLength();
556
557 deltaTime= 0.0; // in case initialVelocity = 0
558 if ( initialVelocity > 0.0 ) { deltaTime = stepLength/initialVelocity; }
559
560 fCandidateEndGlobalTime = startTime + deltaTime ;
561 fParticleChange.ProposeLocalTime( track.GetLocalTime() + deltaTime) ;
562 }
563 else
564 {
565 deltaTime = fCandidateEndGlobalTime - startTime ;
566 fParticleChange.ProposeGlobalTime( fCandidateEndGlobalTime ) ;
567 }
568
569
570 // Now Correct by Lorentz factor to get delta "proper" Time
571
572 G4double restMass = track.GetDynamicParticle()->GetMass() ;
573 G4double deltaProperTime = deltaTime*( restMass/track.GetTotalEnergy() ) ;
574
575 fParticleChange.ProposeProperTime(track.GetProperTime() + deltaProperTime) ;
576 //fParticleChange. ProposeTrueStepLength( track.GetStepLength() ) ;
577
578 // If the particle is caught looping or is stuck (in very difficult
579 // boundaries) in a magnetic field (doing many steps)
580 // THEN this kills it ...
581 //
582 if ( fParticleIsLooping )
583 {
584 G4double endEnergy= fTransportEndKineticEnergy;
585
586 if( (endEnergy < fThreshold_Important_Energy)
587 || (fNoLooperTrials >= fThresholdTrials ) )
588 {
589 // Kill the looping particle
590 //
591 fParticleChange.ProposeTrackStatus( fStopAndKill ) ;
592
593 // 'Bare' statistics
594 fSumEnergyKilled += endEnergy;
595 if( endEnergy > fMaxEnergyKilled) { fMaxEnergyKilled= endEnergy; }
596
597#ifdef G4VERBOSE1
598 if( (fVerboseLevel > 1) &&
599 ( endEnergy > fThreshold_Warning_Energy ) )
600 {
601 G4cout(*G4cout_p) << " G4Transportation is killing track that is looping or stuck "
602 << G4endlstd::endl
603 << " This track has " << track.GetKineticEnergy() / MeV
604 << " MeV energy." << G4endlstd::endl;
605 G4cout(*G4cout_p) << " Number of trials = " << fNoLooperTrials
606 << " No of calls to AlongStepDoIt = " << noCalls
607 << G4endlstd::endl;
608 }
609#endif
610 fNoLooperTrials=0;
611 }
612 else
613 {
614 fNoLooperTrials ++;
615#ifdef G4VERBOSE1
616 if( (fVerboseLevel > 2) )
617 {
618 G4cout(*G4cout_p) << " G4Transportation::AlongStepDoIt(): Particle looping - "
619 << " Number of trials = " << fNoLooperTrials
620 << " No of calls to = " << noCalls
621 << G4endlstd::endl;
622 }
623#endif
624 }
625 }
626 else
627 {
628 fNoLooperTrials=0;
629 }
630
631 // Another (sometimes better way) is to use a user-limit maximum Step size
632 // to alleviate this problem ..
633
634 // Introduce smooth curved trajectories to particle-change
635 //
636 fParticleChange.SetPointerToVectorOfAuxiliaryPoints
637 (fFieldPropagator->GimmeTrajectoryVectorAndForgetIt() );
638
639 return &fParticleChange ;
640}
641
642//////////////////////////////////////////////////////////////////////////
643//
644// This ensures that the PostStep action is always called,
645// so that it can do the relocation if it is needed.
646//
647
648G4double G4Transportation::
649PostStepGetPhysicalInteractionLength( const G4Track&,
650 G4double, // previousStepSize
651 G4ForceCondition* pForceCond )
652{
653 fFieldExertedForce = false; // Not known
654 *pForceCond = Forced ;
655 return DBL_MAXstd::numeric_limits<double>::max() ; // was kInfinity ; but convention now is DBL_MAX
656}
657
658/////////////////////////////////////////////////////////////////////////////
659//
660
661G4VParticleChange* G4Transportation::PostStepDoIt( const G4Track& track,
662 const G4Step& )
663{
664 G4TouchableHandle retCurrentTouchable ; // The one to return
665 G4bool isLastStep= false;
666
667 // Initialize ParticleChange (by setting all its members equal
668 // to corresponding members in G4Track)
669 // fParticleChange.Initialize(track) ; // To initialise TouchableChange
670
671 fParticleChange.ProposeTrackStatus(track.GetTrackStatus()) ;
672
673 // If the Step was determined by the volume boundary,
674 // logically relocate the particle
675
676 if(fGeometryLimitedStep)
677 {
678 // fCurrentTouchable will now become the previous touchable,
679 // and what was the previous will be freed.
680 // (Needed because the preStepPoint can point to the previous touchable)
681
682 fLinearNavigator->SetGeometricallyLimitedStep() ;
683 fLinearNavigator->
684 LocateGlobalPointAndUpdateTouchableHandle( track.GetPosition(),
685 track.GetMomentumDirection(),
686 fCurrentTouchableHandle,
687 true ) ;
688 // Check whether the particle is out of the world volume
689 // If so it has exited and must be killed.
690 //
691 if( fCurrentTouchableHandle->GetVolume() == 0 )
692 {
693 fParticleChange.ProposeTrackStatus( fStopAndKill ) ;
694 }
695 retCurrentTouchable = fCurrentTouchableHandle ;
696 fParticleChange.SetTouchableHandle( fCurrentTouchableHandle ) ;
697
698 // Update the Step flag which identifies the Last Step in a volume
699 if( !fFieldExertedForce )
700 isLastStep = fLinearNavigator->ExitedMotherVolume()
701 | fLinearNavigator->EnteredDaughterVolume() ;
702 else
703 isLastStep = fFieldPropagator->IsLastStepInVolume();
704 }
705 else // fGeometryLimitedStep is false
706 {
707 // This serves only to move the Navigator's location
708 //
709 fLinearNavigator->LocateGlobalPointWithinVolume( track.GetPosition() ) ;
710
711 // The value of the track's current Touchable is retained.
712 // (and it must be correct because we must use it below to
713 // overwrite the (unset) one in particle change)
714 // It must be fCurrentTouchable too ??
715 //
716 fParticleChange.SetTouchableHandle( track.GetTouchableHandle() ) ;
717 retCurrentTouchable = track.GetTouchableHandle() ;
718
719 isLastStep= false;
720 } // endif ( fGeometryLimitedStep )
721 fLastStepInVolume= isLastStep;
722
723 fParticleChange.ProposeFirstStepInVolume(fFirstStepInVolume);
724 fParticleChange.ProposeLastStepInVolume(isLastStep);
725
726 const G4VPhysicalVolume* pNewVol = retCurrentTouchable->GetVolume() ;
727 const G4Material* pNewMaterial = 0 ;
728 const G4VSensitiveDetector* pNewSensitiveDetector = 0 ;
729
730 if( pNewVol != 0 )
731 {
732 pNewMaterial= pNewVol->GetLogicalVolume()->GetMaterial();
733 pNewSensitiveDetector= pNewVol->GetLogicalVolume()->GetSensitiveDetector();
734 }
735
736 // ( <const_cast> pNewMaterial ) ;
737 // ( <const_cast> pNewSensitiveDetector) ;
738
739 fParticleChange.SetMaterialInTouchable( (G4Material *) pNewMaterial ) ;
740 fParticleChange.SetSensitiveDetectorInTouchable( (G4VSensitiveDetector *) pNewSensitiveDetector ) ;
741
742 const G4MaterialCutsCouple* pNewMaterialCutsCouple = 0;
743 if( pNewVol != 0 )
744 {
745 pNewMaterialCutsCouple=pNewVol->GetLogicalVolume()->GetMaterialCutsCouple();
746 }
747
748 if( pNewVol!=0 && pNewMaterialCutsCouple!=0 && pNewMaterialCutsCouple->GetMaterial()!=pNewMaterial )
749 {
750 // for parametrized volume
751 //
752 pNewMaterialCutsCouple =
753 G4ProductionCutsTable::GetProductionCutsTable()
754 ->GetMaterialCutsCouple(pNewMaterial,
755 pNewMaterialCutsCouple->GetProductionCuts());
756 }
757 fParticleChange.SetMaterialCutsCoupleInTouchable( pNewMaterialCutsCouple );
758
759 // temporarily until Get/Set Material of ParticleChange,
760 // and StepPoint can be made const.
761 // Set the touchable in ParticleChange
762 // this must always be done because the particle change always
763 // uses this value to overwrite the current touchable pointer.
764 //
765 fParticleChange.SetTouchableHandle(retCurrentTouchable) ;
766
767 return &fParticleChange ;
768}
769
770// New method takes over the responsibility to reset the state of G4Transportation
771// object at the start of a new track or the resumption of a suspended track.
772
773void
774G4Transportation::StartTracking(G4Track* aTrack)
775{
776 G4VProcess::StartTracking(aTrack);
1
Value assigned to field 'fFieldPropagator'
→
777 fNewTrack= true;
778 fFirstStepInVolume= true;
779 fLastStepInVolume= false;
780
781 // The actions here are those that were taken in AlongStepGPIL
782 // when track.GetCurrentStepNumber()==1
783
784 // reset safety value and center
785 //
786 fPreviousSafety = 0.0 ;
787 fPreviousSftOrigin = G4ThreeVector(0.,0.,0.) ;
788
789 // reset looping counter -- for motion in field
790 fNoLooperTrials= 0;
791 // Must clear this state .. else it depends on last track's value
792 // --> a better solution would set this from state of suspended track TODO ?
793 // Was if( aTrack->GetCurrentStepNumber()==1 ) { .. }
794
795 // ChordFinder reset internal state
796 //
797 if( fFieldPropagator )
2
←
Assuming pointer value is null
→
3
←
Taking false branch
→
798 {
799 fFieldPropagator->ClearPropagatorState();
800 // Resets all state of field propagator class (ONLY)
801 // including safety values (in case of overlaps and to wipe for first track).
802
803 // G4ChordFinder* chordF= fFieldPropagator->GetChordFinder();
804 // if( chordF ) chordF->ResetStepEstimate();
805 }
806
807 // Make sure to clear the chord finders of all fields (ie managers)
808 //
809 G4FieldManagerStore* fieldMgrStore = G4FieldManagerStore::GetInstance();
810 fieldMgrStore->ClearAllChordFindersState();
811
812 // Update the current touchable handle (from the track's)
813 //
814 fCurrentTouchableHandle = aTrack->GetTouchableHandle();
815
816 // Inform field propagator of new track
817 fFieldPropagator->PrepareNewTrack();
4
←
Called C++ object pointer is null
818}
819
820#include "G4CoupledTransportation.hh"
821G4bool G4Transportation::EnableUseMagneticMoment(G4bool useMoment)
822{
823 G4bool lastValue= fUseMagneticMoment;
824 fUseMagneticMoment= useMoment;
825 G4CoupledTransportation::fUseMagneticMoment= useMoment;
826 return lastValue;
827}