Bug Summary

File:scratch/gluex/scan-build-work/hdgeant4/src/G4fixes/G4Transportation.cc
Location:line 811, 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{
109 // set Process Sub Type
110 SetProcessSubType(static_cast<G4int>(TRANSPORTATION));
111 pParticleChange= &fParticleChange; // Required to conform to G4VProcess
112 SetVerboseLevel(verbosity);
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( verboseLevel > 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( (verboseLevel > 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 fLastStepInVolume= false;
186 fNewTrack = false;
187
188 // Get initial Energy/Momentum of the track
189 //
190 const G4DynamicParticle* pParticle = track.GetDynamicParticle() ;
191 const G4ParticleDefinition* pParticleDef = pParticle->GetDefinition() ;
192 G4ThreeVector startMomentumDir = pParticle->GetMomentumDirection() ;
193 G4ThreeVector startPosition = track.GetPosition() ;
194
195 // G4double theTime = track.GetGlobalTime() ;
196
197 // The Step Point safety can be limited by other geometries and/or the
198 // assumptions of any process - it's not always the geometrical safety.
199 // We calculate the starting point's isotropic safety here.
200 //
201 G4ThreeVector OriginShift = startPosition - fPreviousSftOrigin ;
202 G4double MagSqShift = OriginShift.mag2() ;
203 if( MagSqShift >= sqr(fPreviousSafety) )
204 {
205 currentSafety = 0.0 ;
206 }
207 else
208 {
209 currentSafety = fPreviousSafety - std::sqrt(MagSqShift) ;
210 }
211
212 // Is the particle charged or has it a magnetic moment?
213 //
214 G4double particleCharge = pParticle->GetCharge() ;
215 G4double magneticMoment = pParticle->GetMagneticMoment() ;
216 G4double restMass = pParticle->GetMass() ;
217
218 fGeometryLimitedStep = false ;
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 G4bool fieldExists= false; // Field is not 0 (null pointer)
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 // Is here to allow a transition from no-field pointer
239 // to finite field (non-zero pointer).
240
241 // If the field manager has no field ptr, the field is zero
242 // by definition ( = there is no field ! )
243 const G4Field* ptrField= fieldMgr->GetDetectorField();
244 fieldExists = (ptrField!=0) ;
245 if( fieldExists )
246 {
247 gravityOn= ptrField->IsGravityActive();
248
249 if( (particleCharge != 0.0)
250 || (fUseMagneticMoment && (magneticMoment != 0.0) )
251 || (gravityOn && (restMass != 0.0) )
252 )
253 {
254 fieldExertsForce = fieldExists;
255 }
256 }
257 }
258 // G4cout << " G4Transport: field exerts force= " << fieldExertsForce
259 // << " fieldMgr= " << fieldMgr << G4endl;
260 fFieldExertedForce = fieldExertsForce;
261
262 if( !fieldExertsForce )
263 {
264 G4double linearStepLength ;
265 if( fShortStepOptimisation && (currentMinimumStep <= currentSafety) )
266 {
267 // The Step is guaranteed to be taken
268 //
269 geometryStepLength = currentMinimumStep ;
270 fGeometryLimitedStep = false ;
271 }
272 else
273 {
274 // Find whether the straight path intersects a volume
275 //
276 linearStepLength = fLinearNavigator->ComputeStep( startPosition,
277 startMomentumDir,
278 currentMinimumStep,
279 newSafety) ;
280 // Remember last safety origin & value.
281 //
282 fPreviousSftOrigin = startPosition ;
283 fPreviousSafety = newSafety ;
284 fpSafetyHelper->SetCurrentSafety( newSafety, startPosition);
285
286 currentSafety = newSafety ;
287
288 fGeometryLimitedStep= (linearStepLength <= currentMinimumStep);
289 if( fGeometryLimitedStep )
290 {
291 // The geometry limits the Step size (an intersection was found.)
292 geometryStepLength = linearStepLength ;
293 }
294 else
295 {
296 // The full Step is taken.
297 geometryStepLength = currentMinimumStep ;
298 }
299 }
300 fEndPointDistance = geometryStepLength ;
301
302 // Calculate final position
303 //
304 fTransportEndPosition = startPosition+geometryStepLength*startMomentumDir ;
305
306 // Momentum direction, energy and polarisation are unchanged by transport
307 //
308 fTransportEndMomentumDir = startMomentumDir ;
309 fTransportEndKineticEnergy = track.GetKineticEnergy() ;
310 fTransportEndSpin = track.GetPolarization();
311 fParticleIsLooping = false ;
312 fMomentumChanged = false ;
313 fEndGlobalTimeComputed = false ;
314 }
315 else // A field exerts force
316 {
317 G4double momentumMagnitude = pParticle->GetTotalMomentum() ;
318 G4ThreeVector EndUnitMomentum ;
319 G4double lengthAlongCurve ;
320
321 G4ChargeState chargeState(particleCharge, // The charge can change (dynamic)
322 magneticMoment,
323 pParticleDef->GetPDGSpin() );
324 // For insurance, could set it again
325 // chargeState.SetPDGSpin(pParticleDef->GetPDGSpin() ); // Provisionally in same object
326
327 G4EquationOfMotion* equationOfMotion = fFieldPropagator->GetCurrentEquationOfMotion();
328
329 equationOfMotion->SetChargeMomentumMass( chargeState,
330 momentumMagnitude,
331 restMass);
332
333 G4FieldTrack aFieldTrack = G4FieldTrack( startPosition,
334 track.GetGlobalTime(), // Lab.
335 // track.GetProperTime(), // Particle rest frame
336 track.GetMomentumDirection(),
337 track.GetKineticEnergy(),
338 restMass,
339 particleCharge,
340 track.GetPolarization(),
341 pParticleDef->GetPDGMagneticMoment(),
342 0.0, // Length along track
343 pParticleDef->GetPDGSpin()
344 ) ;
345
346 if( currentMinimumStep > 0 )
347 {
348 // Do the Transport in the field (non recti-linear)
349 //
350 lengthAlongCurve = fFieldPropagator->ComputeStep( aFieldTrack,
351 currentMinimumStep,
352 currentSafety,
353 track.GetVolume() ) ;
354
355 fGeometryLimitedStep= fFieldPropagator->IsLastStepInVolume();
356 // It is possible that step was reduced in PropagatorInField due to previous zero steps
357 // To cope with case that reduced step is taken in full, we must rely on PiF to obtain this
358 // value.
359
360 geometryStepLength = std::min( lengthAlongCurve, currentMinimumStep );
361
362 // Remember last safety origin & value.
363 //
364 fPreviousSftOrigin = startPosition ;
365 fPreviousSafety = currentSafety ;
366 fpSafetyHelper->SetCurrentSafety( currentSafety, startPosition);
367 }
368 else
369 {
370 geometryStepLength = lengthAlongCurve= 0.0 ;
371 fGeometryLimitedStep = false ;
372 }
373
374 // Get the End-Position and End-Momentum (Dir-ection)
375 //
376 fTransportEndPosition = aFieldTrack.GetPosition() ;
377
378 // Momentum: Magnitude and direction can be changed too now ...
379 //
380 fMomentumChanged = true ;
381 fTransportEndMomentumDir = aFieldTrack.GetMomentumDir() ;
382
383 fTransportEndKineticEnergy = aFieldTrack.GetKineticEnergy() ;
384
385 if( fFieldPropagator->GetCurrentFieldManager()->DoesFieldChangeEnergy() )
386 {
387 // If the field can change energy, then the time must be integrated
388 // - so this should have been updated
389 //
390 fCandidateEndGlobalTime = aFieldTrack.GetLabTimeOfFlight();
391 fEndGlobalTimeComputed = true;
392
393 // was ( fCandidateEndGlobalTime != track.GetGlobalTime() );
394 // a cleaner way is to have FieldTrack knowing whether time is updated.
395 }
396 else
397 {
398 // The energy should be unchanged by field transport,
399 // - so the time changed will be calculated elsewhere
400 //
401 fEndGlobalTimeComputed = false;
402
403 // Check that the integration preserved the energy
404 // - and if not correct this!
405 G4double startEnergy= track.GetKineticEnergy();
406 G4double endEnergy= fTransportEndKineticEnergy;
407
408 static G4ThreadLocalthread_local G4int no_inexact_steps=0, no_large_ediff;
409 G4double absEdiff = std::fabs(startEnergy- endEnergy);
410 if( absEdiff > perMillion * endEnergy )
411 {
412 no_inexact_steps++;
413 // Possible statistics keeping here ...
414 }
415 if( verboseLevel > 1 )
416 {
417 if( std::fabs(startEnergy- endEnergy) > perThousand * endEnergy )
418 {
419 static G4ThreadLocalthread_local G4int no_warnings= 0, warnModulo=1,
420 moduloFactor= 10;
421 no_large_ediff ++;
422 if( (no_large_ediff% warnModulo) == 0 )
423 {
424 no_warnings++;
425 G4cout(*G4cout_p) << "WARNING - G4Transportation::AlongStepGetPIL() "
426 << " Energy change in Step is above 1^-3 relative value. " << G4endlstd::endl
427 << " Relative change in 'tracking' step = "
428 << std::setw(15) << (endEnergy-startEnergy)/startEnergy << G4endlstd::endl
429 << " Starting E= " << std::setw(12) << startEnergy / MeV << " MeV " << G4endlstd::endl
430 << " Ending E= " << std::setw(12) << endEnergy / MeV << " MeV " << G4endlstd::endl;
431 G4cout(*G4cout_p) << " Energy has been corrected -- however, review"
432 << " field propagation parameters for accuracy." << G4endlstd::endl;
433 if( (verboseLevel > 2 ) || (no_warnings<4) || (no_large_ediff == warnModulo * moduloFactor) )
434 {
435 G4cout(*G4cout_p) << " These include EpsilonStepMax(/Min) in G4FieldManager "
436 << " which determine fractional error per step for integrated quantities. " << G4endlstd::endl
437 << " Note also the influence of the permitted number of integration steps."
438 << G4endlstd::endl;
439 }
440 G4cerr(*G4cerr_p) << "ERROR - G4Transportation::AlongStepGetPIL()" << G4endlstd::endl
441 << " Bad 'endpoint'. Energy change detected"
442 << " and corrected. "
443 << " Has occurred already "
444 << no_large_ediff << " times." << G4endlstd::endl;
445 if( no_large_ediff == warnModulo * moduloFactor )
446 {
447 warnModulo *= moduloFactor;
448 }
449 }
450 }
451 } // end of if (verboseLevel)
452
453 // Correct the energy for fields that conserve it
454 // This - hides the integration error
455 // - but gives a better physical answer
456 fTransportEndKineticEnergy= track.GetKineticEnergy();
457 }
458
459 fTransportEndSpin = aFieldTrack.GetSpin();
460 fParticleIsLooping = fFieldPropagator->IsParticleLooping() ;
461 fEndPointDistance = (fTransportEndPosition - startPosition).mag() ;
462 }
463
464 // If we are asked to go a step length of 0, and we are on a boundary
465 // then a boundary will also limit the step -> we must flag this.
466 //
467 if( currentMinimumStep == 0.0 )
468 {
469 if( currentSafety == 0.0 ) { fGeometryLimitedStep = true; }
470 }
471
472 // Update the safety starting from the end-point,
473 // if it will become negative at the end-point.
474 //
475 if( currentSafety < fEndPointDistance )
476 {
477 if( particleCharge != 0.0 )
478 {
479 G4double endSafety =
480 fLinearNavigator->ComputeSafety( fTransportEndPosition) ;
481 currentSafety = endSafety ;
482 fPreviousSftOrigin = fTransportEndPosition ;
483 fPreviousSafety = currentSafety ;
484 fpSafetyHelper->SetCurrentSafety( currentSafety, fTransportEndPosition);
485
486 // Because the Stepping Manager assumes it is from the start point,
487 // add the StepLength
488 //
489 currentSafety += fEndPointDistance ;
490
491#ifdef G4DEBUG_TRANSPORT
492 G4cout(*G4cout_p).precision(12) ;
493 G4cout(*G4cout_p) << "***G4Transportation::AlongStepGPIL ** " << G4endlstd::endl ;
494 G4cout(*G4cout_p) << " Called Navigator->ComputeSafety at " << fTransportEndPosition
495 << " and it returned safety= " << endSafety << G4endlstd::endl ;
496 G4cout(*G4cout_p) << " Adding endpoint distance " << fEndPointDistance
497 << " to obtain pseudo-safety= " << currentSafety << G4endlstd::endl ;
498 }
499 else
500 {
501 G4cout(*G4cout_p) << "***G4Transportation::AlongStepGPIL ** " << G4endlstd::endl ;
502 G4cout(*G4cout_p) << " Avoiding call to ComputeSafety : " << G4endlstd::endl;
503 G4cout(*G4cout_p) << " charge = " << particleCharge << G4endlstd::endl;
504 G4cout(*G4cout_p) << " mag moment = " << magneticMoment << G4endlstd::endl;
505#endif
506 }
507 }
508
509 fParticleChange.ProposeTrueStepLength(geometryStepLength) ;
510
511 return geometryStepLength ;
512}
513
514//////////////////////////////////////////////////////////////////////////
515//
516// Initialize ParticleChange (by setting all its members equal
517// to corresponding members in G4Track)
518
519G4VParticleChange* G4Transportation::AlongStepDoIt( const G4Track& track,
520 const G4Step& stepData )
521{
522 static G4ThreadLocalthread_local G4int noCalls=0;
523 noCalls++;
524
525 fParticleChange.Initialize(track) ;
526
527 // Code for specific process
528 //
529 fParticleChange.ProposePosition(fTransportEndPosition) ;
530 fParticleChange.ProposeMomentumDirection(fTransportEndMomentumDir) ;
531 fParticleChange.ProposeEnergy(fTransportEndKineticEnergy) ;
532 fParticleChange.SetMomentumChanged(fMomentumChanged) ;
533
534 fParticleChange.ProposePolarization(fTransportEndSpin);
535
536 G4double deltaTime = 0.0 ;
537
538 // Calculate Lab Time of Flight (ONLY if field Equations used it!)
539 // G4double endTime = fCandidateEndGlobalTime;
540 // G4double delta_time = endTime - startTime;
541
542 G4double startTime = track.GetGlobalTime() ;
543
544 if (!fEndGlobalTimeComputed)
545 {
546 // The time was not integrated .. make the best estimate possible
547 //
548 G4double initialVelocity = stepData.GetPreStepPoint()->GetVelocity();
549 G4double stepLength = track.GetStepLength();
550
551 deltaTime= 0.0; // in case initialVelocity = 0
552 if ( initialVelocity > 0.0 ) { deltaTime = stepLength/initialVelocity; }
553
554 fCandidateEndGlobalTime = startTime + deltaTime ;
555 fParticleChange.ProposeLocalTime( track.GetLocalTime() + deltaTime) ;
556 }
557 else
558 {
559 deltaTime = fCandidateEndGlobalTime - startTime ;
560 fParticleChange.ProposeGlobalTime( fCandidateEndGlobalTime ) ;
561 }
562
563
564 // Now Correct by Lorentz factor to get delta "proper" Time
565
566 G4double restMass = track.GetDynamicParticle()->GetMass() ;
567 G4double deltaProperTime = deltaTime*( restMass/track.GetTotalEnergy() ) ;
568
569 fParticleChange.ProposeProperTime(track.GetProperTime() + deltaProperTime) ;
570 //fParticleChange. ProposeTrueStepLength( track.GetStepLength() ) ;
571
572 // If the particle is caught looping or is stuck (in very difficult
573 // boundaries) in a magnetic field (doing many steps)
574 // THEN this kills it ...
575 //
576 if ( fParticleIsLooping )
577 {
578 G4double endEnergy= fTransportEndKineticEnergy;
579
580 if( (endEnergy < fThreshold_Important_Energy)
581 || (fNoLooperTrials >= fThresholdTrials ) )
582 {
583 // Kill the looping particle
584 //
585 fParticleChange.ProposeTrackStatus( fStopAndKill ) ;
586
587 // 'Bare' statistics
588 fSumEnergyKilled += endEnergy;
589 if( endEnergy > fMaxEnergyKilled) { fMaxEnergyKilled= endEnergy; }
590
591#ifdef G4VERBOSE1
592 if( (verboseLevel > 1) &&
593 ( endEnergy > fThreshold_Warning_Energy ) )
594 {
595 G4cout(*G4cout_p) << " G4Transportation is killing track that is looping or stuck "
596 << G4endlstd::endl
597 << " This track has " << track.GetKineticEnergy() / MeV
598 << " MeV energy." << G4endlstd::endl;
599 G4cout(*G4cout_p) << " Number of trials = " << fNoLooperTrials
600 << " No of calls to AlongStepDoIt = " << noCalls
601 << G4endlstd::endl;
602 }
603#endif
604 fNoLooperTrials=0;
605 }
606 else
607 {
608 fNoLooperTrials ++;
609#ifdef G4VERBOSE1
610 if( (verboseLevel > 2) )
611 {
612 G4cout(*G4cout_p) << " G4Transportation::AlongStepDoIt(): Particle looping - "
613 << " Number of trials = " << fNoLooperTrials
614 << " No of calls to = " << noCalls
615 << G4endlstd::endl;
616 }
617#endif
618 }
619 }
620 else
621 {
622 fNoLooperTrials=0;
623 }
624
625 // Another (sometimes better way) is to use a user-limit maximum Step size
626 // to alleviate this problem ..
627
628 // Introduce smooth curved trajectories to particle-change
629 //
630 fParticleChange.SetPointerToVectorOfAuxiliaryPoints
631 (fFieldPropagator->GimmeTrajectoryVectorAndForgetIt() );
632
633 return &fParticleChange ;
634}
635
636//////////////////////////////////////////////////////////////////////////
637//
638// This ensures that the PostStep action is always called,
639// so that it can do the relocation if it is needed.
640//
641
642G4double G4Transportation::
643PostStepGetPhysicalInteractionLength( const G4Track&,
644 G4double, // previousStepSize
645 G4ForceCondition* pForceCond )
646{
647 fFieldExertedForce = false; // Not known
648 *pForceCond = Forced ;
649 return DBL_MAXstd::numeric_limits<double>::max() ; // was kInfinity ; but convention now is DBL_MAX
650}
651
652/////////////////////////////////////////////////////////////////////////////
653//
654
655G4VParticleChange* G4Transportation::PostStepDoIt( const G4Track& track,
656 const G4Step& )
657{
658 G4TouchableHandle retCurrentTouchable ; // The one to return
659 G4bool isLastStep= false;
660
661 // Initialize ParticleChange (by setting all its members equal
662 // to corresponding members in G4Track)
663 // fParticleChange.Initialize(track) ; // To initialise TouchableChange
664
665 fParticleChange.ProposeTrackStatus(track.GetTrackStatus()) ;
666
667 // If the Step was determined by the volume boundary,
668 // logically relocate the particle
669
670 if(fGeometryLimitedStep)
671 {
672 // fCurrentTouchable will now become the previous touchable,
673 // and what was the previous will be freed.
674 // (Needed because the preStepPoint can point to the previous touchable)
675
676 fLinearNavigator->SetGeometricallyLimitedStep() ;
677 fLinearNavigator->
678 LocateGlobalPointAndUpdateTouchableHandle( track.GetPosition(),
679 track.GetMomentumDirection(),
680 fCurrentTouchableHandle,
681 true ) ;
682 // Check whether the particle is out of the world volume
683 // If so it has exited and must be killed.
684 //
685 if( fCurrentTouchableHandle->GetVolume() == 0 )
686 {
687 fParticleChange.ProposeTrackStatus( fStopAndKill ) ;
688 }
689 retCurrentTouchable = fCurrentTouchableHandle ;
690 fParticleChange.SetTouchableHandle( fCurrentTouchableHandle ) ;
691
692 // Update the Step flag which identifies the Last Step in a volume
693 if( !fFieldExertedForce )
694 isLastStep = fLinearNavigator->ExitedMotherVolume()
695 | fLinearNavigator->EnteredDaughterVolume() ;
696 else
697 isLastStep = fFieldPropagator->IsLastStepInVolume();
698 }
699 else // fGeometryLimitedStep is false
700 {
701 // This serves only to move the Navigator's location
702 //
703 fLinearNavigator->LocateGlobalPointWithinVolume( track.GetPosition() ) ;
704
705 // The value of the track's current Touchable is retained.
706 // (and it must be correct because we must use it below to
707 // overwrite the (unset) one in particle change)
708 // It must be fCurrentTouchable too ??
709 //
710 fParticleChange.SetTouchableHandle( track.GetTouchableHandle() ) ;
711 retCurrentTouchable = track.GetTouchableHandle() ;
712
713 isLastStep= false;
714 } // endif ( fGeometryLimitedStep )
715 fLastStepInVolume= isLastStep;
716
717 fParticleChange.ProposeFirstStepInVolume(fFirstStepInVolume);
718 fParticleChange.ProposeLastStepInVolume(isLastStep);
719
720 const G4VPhysicalVolume* pNewVol = retCurrentTouchable->GetVolume() ;
721 const G4Material* pNewMaterial = 0 ;
722 const G4VSensitiveDetector* pNewSensitiveDetector = 0 ;
723
724 if( pNewVol != 0 )
725 {
726 pNewMaterial= pNewVol->GetLogicalVolume()->GetMaterial();
727 pNewSensitiveDetector= pNewVol->GetLogicalVolume()->GetSensitiveDetector();
728 }
729
730 // ( <const_cast> pNewMaterial ) ;
731 // ( <const_cast> pNewSensitiveDetector) ;
732
733 fParticleChange.SetMaterialInTouchable( (G4Material *) pNewMaterial ) ;
734 fParticleChange.SetSensitiveDetectorInTouchable( (G4VSensitiveDetector *) pNewSensitiveDetector ) ;
735
736 const G4MaterialCutsCouple* pNewMaterialCutsCouple = 0;
737 if( pNewVol != 0 )
738 {
739 pNewMaterialCutsCouple=pNewVol->GetLogicalVolume()->GetMaterialCutsCouple();
740 }
741
742 if( pNewVol!=0 && pNewMaterialCutsCouple!=0 && pNewMaterialCutsCouple->GetMaterial()!=pNewMaterial )
743 {
744 // for parametrized volume
745 //
746 pNewMaterialCutsCouple =
747 G4ProductionCutsTable::GetProductionCutsTable()
748 ->GetMaterialCutsCouple(pNewMaterial,
749 pNewMaterialCutsCouple->GetProductionCuts());
750 }
751 fParticleChange.SetMaterialCutsCoupleInTouchable( pNewMaterialCutsCouple );
752
753 // temporarily until Get/Set Material of ParticleChange,
754 // and StepPoint can be made const.
755 // Set the touchable in ParticleChange
756 // this must always be done because the particle change always
757 // uses this value to overwrite the current touchable pointer.
758 //
759 fParticleChange.SetTouchableHandle(retCurrentTouchable) ;
760
761 return &fParticleChange ;
762}
763
764// New method takes over the responsibility to reset the state of G4Transportation
765// object at the start of a new track or the resumption of a suspended track.
766
767void
768G4Transportation::StartTracking(G4Track* aTrack)
769{
770 G4VProcess::StartTracking(aTrack);
1
Value assigned to field 'fFieldPropagator'
→
771 fNewTrack= true;
772 fFirstStepInVolume= true;
773 fLastStepInVolume= false;
774
775 // The actions here are those that were taken in AlongStepGPIL
776 // when track.GetCurrentStepNumber()==1
777
778 // reset safety value and center
779 //
780 fPreviousSafety = 0.0 ;
781 fPreviousSftOrigin = G4ThreeVector(0.,0.,0.) ;
782
783 // reset looping counter -- for motion in field
784 fNoLooperTrials= 0;
785 // Must clear this state .. else it depends on last track's value
786 // --> a better solution would set this from state of suspended track TODO ?
787 // Was if( aTrack->GetCurrentStepNumber()==1 ) { .. }
788
789 // ChordFinder reset internal state
790 //
791 if( fFieldPropagator )
2
←
Assuming pointer value is null
→
3
←
Taking false branch
→
792 {
793 fFieldPropagator->ClearPropagatorState();
794 // Resets all state of field propagator class (ONLY)
795 // including safety values (in case of overlaps and to wipe for first track).
796
797 // G4ChordFinder* chordF= fFieldPropagator->GetChordFinder();
798 // if( chordF ) chordF->ResetStepEstimate();
799 }
800
801 // Make sure to clear the chord finders of all fields (ie managers)
802 //
803 G4FieldManagerStore* fieldMgrStore = G4FieldManagerStore::GetInstance();
804 fieldMgrStore->ClearAllChordFindersState();
805
806 // Update the current touchable handle (from the track's)
807 //
808 fCurrentTouchableHandle = aTrack->GetTouchableHandle();
809
810 // Inform field propagator of new track
811 fFieldPropagator->PrepareNewTrack();
4
←
Called C++ object pointer is null
812}
813
814#include "G4CoupledTransportation.hh"
815G4bool G4Transportation::EnableUseMagneticMoment(G4bool useMoment)
816{
817 G4bool lastValue= fUseMagneticMoment;
818 fUseMagneticMoment= useMoment;
819 G4CoupledTransportation::fUseMagneticMoment= useMoment;
820 return lastValue;
821}