Geant4  10.01.p02
G4VParticleChange.cc
Go to the documentation of this file.
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: G4VParticleChange.cc 78546 2014-01-07 09:35:55Z gcosmo $
28 //
29 //
30 // --------------------------------------------------------------
31 // GEANT 4 class implementation file
32 //
33 //
34 // ------------------------------------------------------------
35 // Implemented for the new scheme 23 Mar. 1998 H.Kurahige
36 // --------------------------------------------------------------
37 
38 #include "G4VParticleChange.hh"
39 #include "G4SystemOfUnits.hh"
40 #include "G4Track.hh"
41 #include "G4Step.hh"
42 #include "G4TrackFastVector.hh"
43 #include "G4ExceptionSeverity.hh"
44 
47 
49  :theListOfSecondaries(0),
50  theNumberOfSecondaries(0),
51  theSizeOftheListOfSecondaries(G4TrackFastVectorSize),
52  theStatusChange(fAlive),
53  theSteppingControlFlag(NormalCondition),
54  theLocalEnergyDeposit(0.0),
55  theNonIonizingEnergyDeposit(0.0),
56  theTrueStepLength(0.0),
57  theFirstStepInVolume(false),
58  theLastStepInVolume(false),
59  theParentWeight(1.0),
60  isParentWeightProposed(false),
61  fSetSecondaryWeightByProcess(false),
62  theParentGlobalTime(0.0),
63  verboseLevel(1),
64  debugFlag(false)
65 {
66 #ifdef G4VERBOSE
67  // activate CHeckIt if in VERBOSE mode
68  debugFlag = true;
69 #endif
71 }
72 
74  // check if tracks still exist in theListOfSecondaries
75  if (theNumberOfSecondaries>0) {
76 #ifdef G4VERBOSE
77  if (verboseLevel>0) {
78  G4cout << "G4VParticleChange::~G4VParticleChange() Warning ";
79  G4cout << "theListOfSecondaries is not empty ";
80  }
81 #endif
82  for (G4int index= 0; index<theNumberOfSecondaries; index++){
83  delete (*theListOfSecondaries)[index] ;
84  }
85  }
86  delete theListOfSecondaries;
87 }
88 
90  :theListOfSecondaries(0),
91  theNumberOfSecondaries(0),
92  theSizeOftheListOfSecondaries(G4TrackFastVectorSize),
93  theStatusChange( right.theStatusChange),
94  theSteppingControlFlag(right.theSteppingControlFlag),
95  theLocalEnergyDeposit(right.theLocalEnergyDeposit),
96  theNonIonizingEnergyDeposit(right.theNonIonizingEnergyDeposit),
97  theTrueStepLength(right.theTrueStepLength),
98  theFirstStepInVolume( right.theFirstStepInVolume),
99  theLastStepInVolume(right.theLastStepInVolume),
100  theParentWeight(right.theParentWeight),
101  isParentWeightProposed(false),
102  fSetSecondaryWeightByProcess(right.fSetSecondaryWeightByProcess),
103  theParentGlobalTime(0.0),
104  verboseLevel(right.verboseLevel),
105  debugFlag(right.debugFlag)
106 {
109  for (G4int index = 0; index<theNumberOfSecondaries; index++){
110  G4Track* newTrack = new G4Track(*((*right.theListOfSecondaries)[index] ));
111  theListOfSecondaries->SetElement(index, newTrack);
112  }
113 }
114 
115 
117 {
118  if (this != &right){
119  if (theNumberOfSecondaries>0) {
120 #ifdef G4VERBOSE
121  if (verboseLevel>0) {
122  G4cout << "G4VParticleChange: assignment operator Warning ";
123  G4cout << "theListOfSecondaries is not empty ";
124  }
125 #endif
126  for (G4int index = 0; index<theNumberOfSecondaries; index++){
127  if ( (*theListOfSecondaries)[index] ) delete ((*theListOfSecondaries)[index]) ;
128  }
129  }
130  delete theListOfSecondaries;
131 
134  for (G4int index = 0; index<theNumberOfSecondaries; index++){
135  G4Track* newTrack = new G4Track(*((*right.theListOfSecondaries)[index] ));
136  theListOfSecondaries->SetElement(index, newTrack);
137  }
143 
146 
150 
152 
153  verboseLevel = right.verboseLevel;
154  debugFlag = right.debugFlag;
155 
156  }
157  return *this;
158 }
159 
161 {
162  return (this == (G4VParticleChange *) &right);
163 }
164 
165 
167 {
168  return (this != (G4VParticleChange *) &right);
169 }
170 
172 {
173  if (debugFlag) CheckSecondary(*aTrack);
174 
175  // add a secondary after size check
177  // Set weight of secondary tracks
181  } else {
182  delete aTrack;
183 
184 #ifdef G4VERBOSE
185  if (verboseLevel>0) {
186  G4cout << "G4VParticleChange::AddSecondary() Warning ";
187  G4cout << "theListOfSecondaries is full !! " << G4endl;
188  G4cout << " The track is deleted " << G4endl;
189  }
190 #endif
191  G4Exception("G4VParticleChange::AddSecondary",
192  "TRACK101", JustWarning,
193  "Secondary Bug is full. The track is deleted");
194  }
195 }
196 
197 
198 
199 // Virtual methods for updating G4Step
200 //
201 
203 {
204  // Update the G4Step specific attributes
209 
210  if (theFirstStepInVolume) {pStep->SetFirstStepFlag();}
211  else {pStep->ClearFirstStepFlag();}
212  if (theLastStepInVolume) {pStep->SetLastStepFlag();}
213  else {pStep->ClearLastStepFlag();}
214 
215  return pStep;
216 }
217 
218 
220 {
223  }
224  return UpdateStepInfo(Step);
225 }
226 
227 
229 {
230  if (isParentWeightProposed ){
231  G4double initialWeight = Step->GetPreStepPoint()->GetWeight();
232  G4double currentWeight = Step->GetPostStepPoint()->GetWeight();
233  G4double finalWeight = (theParentWeight/initialWeight)*currentWeight;
234  Step->GetPostStepPoint()->SetWeight( finalWeight );
235  }
236  return UpdateStepInfo(Step);
237 }
238 
240 {
241  if (isParentWeightProposed ){
243  }
244  return UpdateStepInfo(Step);
245 }
246 
247 
248 //----------------------------------------------------------------
249 // methods for printing messages
250 //
251 
253 {
254 
255 // Show header
256  G4int olprc = G4cout.precision(3);
257  G4cout << " -----------------------------------------------"
258  << G4endl;
259  G4cout << " G4ParticleChange Information " << std::setw(20) << G4endl;
260  G4cout << " -----------------------------------------------"
261  << G4endl;
262 
263  G4cout << " # of 2ndaries : "
264  << std::setw(20) << theNumberOfSecondaries
265  << G4endl;
266 
267  if (theNumberOfSecondaries >0) {
268  G4cout << " Pointer to 2ndaries : "
269  << std::setw(20) << GetSecondary(0)
270  << G4endl;
271  G4cout << " (Showed only 1st one)"
272  << G4endl;
273  }
274  G4cout << " -----------------------------------------------"
275  << G4endl;
276 
277  G4cout << " Energy Deposit (MeV): "
278  << std::setw(20) << theLocalEnergyDeposit/MeV
279  << G4endl;
280 
281  G4cout << " Non-ionizing Energy Deposit (MeV): "
282  << std::setw(20) << theNonIonizingEnergyDeposit/MeV
283  << G4endl;
284 
285  G4cout << " Track Status : "
286  << std::setw(20);
287  if( theStatusChange == fAlive ){
288  G4cout << " Alive";
289  } else if( theStatusChange == fStopButAlive ){
290  G4cout << " StopButAlive";
291  } else if( theStatusChange == fStopAndKill ){
292  G4cout << " StopAndKill";
294  G4cout << " KillTrackAndSecondaries";
295  } else if( theStatusChange == fSuspend ){
296  G4cout << " Suspend";
297  } else if( theStatusChange == fPostponeToNextEvent ){
298  G4cout << " PostponeToNextEvent";
299  }
300  G4cout << G4endl;
301  G4cout << " True Path Length (mm) : "
302  << std::setw(20) << theTrueStepLength/mm
303  << G4endl;
304  G4cout << " Stepping Control : "
305  << std::setw(20) << theSteppingControlFlag
306  << G4endl;
307  if (theFirstStepInVolume) {
308  G4cout << " First Step In the voulme : " << G4endl;
309  }
310  if (theLastStepInVolume) {
311  G4cout << " Last Step In the voulme : " << G4endl;
312  }
313  G4cout.precision(olprc);
314 }
315 
317 #ifdef G4VERBOSE
318  aTrack
319 #endif
320 )
321 {
322 
323  G4bool exitWithError = false;
324  G4double accuracy;
325  static G4ThreadLocal G4int nError = 0;
326 #ifdef G4VERBOSE
327  const G4int maxError = 30;
328 #endif
329 
330  // Energy deposit should not be negative
331  G4bool itsOKforEnergy = true;
332  accuracy = -1.0*theLocalEnergyDeposit/MeV;
333  if (accuracy > accuracyForWarning) {
334  itsOKforEnergy = false;
335  nError += 1;
336  exitWithError = (accuracy > accuracyForException);
337 #ifdef G4VERBOSE
338  if (nError < maxError) {
339  G4cout << " G4VParticleChange::CheckIt : ";
340  G4cout << "the energy deposit is negative !!"
341  << " Difference: " << accuracy << "[MeV] " <<G4endl;
342  G4cout << aTrack.GetDefinition()->GetParticleName()
343  << " E=" << aTrack.GetKineticEnergy()/MeV
344  << " pos=" << aTrack.GetPosition().x()/m
345  << ", " << aTrack.GetPosition().y()/m
346  << ", " << aTrack.GetPosition().z()/m
347  <<G4endl;
348  }
349 #endif
350  }
351 
352  // true path length should not be negative
353  G4bool itsOKforStepLength = true;
354  accuracy = -1.0*theTrueStepLength/mm;
355  if (accuracy > accuracyForWarning) {
356  itsOKforStepLength = false;
357  nError += 1;
358  exitWithError = (accuracy > accuracyForException);
359 #ifdef G4VERBOSE
360  if (nError < maxError) {
361  G4cout << " G4VParticleChange::CheckIt : ";
362  G4cout << "the true step length is negative !!"
363  << " Difference: " << accuracy << "[MeV] " <<G4endl;
364  G4cout << aTrack.GetDefinition()->GetParticleName()
365  << " E=" << aTrack.GetKineticEnergy()/MeV
366  << " pos=" << aTrack.GetPosition().x()/m
367  << ", " << aTrack.GetPosition().y()/m
368  << ", " << aTrack.GetPosition().z()/m
369  <<G4endl;
370  }
371 #endif
372 
373  }
374 #ifdef G4VERBOSE
375  if (!itsOKforStepLength || !itsOKforEnergy ){
376  // dump out information of this particle change
377  DumpInfo();
378  }
379 #endif
380 
381  // Exit with error
382  if (exitWithError) {
383  G4Exception("G4VParticleChange::CheckIt",
384  "TRACK001", EventMustBeAborted,
385  "Step length and/or energy deposit was illegal");
386  }
387 
388  // correction
389  if ( !itsOKforStepLength ) {
390  theTrueStepLength = (1.e-12)*mm;
391  }
392  if ( !itsOKforEnergy ) {
393  theLocalEnergyDeposit = 0.0;
394  }
395  return (itsOKforStepLength && itsOKforEnergy );
396 }
397 
399 {
400  G4bool exitWithError = false;
401  G4double accuracy;
402  static G4ThreadLocal G4int nError = 0;
403 #ifdef G4VERBOSE
404  const G4int maxError = 30;
405 #endif
406 
407  // MomentumDirection should be unit vector
408  G4bool itsOKforMomentum = true;
409  if (aTrack.GetKineticEnergy()>0.){
410  accuracy = std::fabs((aTrack.GetMomentumDirection()).mag2()-1.0);
411  if (accuracy > accuracyForWarning) {
412  itsOKforMomentum = false;
413  nError += 1;
414  exitWithError = exitWithError || (accuracy > accuracyForException);
415 #ifdef G4VERBOSE
416  if (nError < maxError) {
417  G4cout << " G4VParticleChange::CheckSecondary : ";
418  G4cout << "the Momentum direction is not unit vector !! "
419  << " Difference: " << accuracy << G4endl;
420  G4cout << aTrack.GetDefinition()->GetParticleName()
421  << " E=" << aTrack.GetKineticEnergy()/MeV
422  << " pos=" << aTrack.GetPosition().x()/m
423  << ", " << aTrack.GetPosition().y()/m
424  << ", " << aTrack.GetPosition().z()/m
425  <<G4endl;
426  }
427 #endif
428  }
429  }
430 
431  // Kinetic Energy should not be negative
432  G4bool itsOKforEnergy = true;
433  accuracy = -1.0*(aTrack.GetKineticEnergy())/MeV;
434  if (accuracy > accuracyForWarning) {
435  itsOKforEnergy = false;
436  nError += 1;
437  exitWithError = exitWithError || (accuracy > accuracyForException);
438 #ifdef G4VERBOSE
439  if (nError < maxError) {
440  G4cout << " G4VParticleChange::CheckSecondary : ";
441  G4cout << "the kinetic energy is negative !!"
442  << " Difference: " << accuracy << "[MeV] " <<G4endl;
443  G4cout << " G4VParticleChange::CheckSecondary : ";
444  G4cout << "the global time of secondary is earlier than the parent !!"
445  << " Difference: " << accuracy << "[ns] " <<G4endl;
446  G4cout << aTrack.GetDefinition()->GetParticleName()
447  << " E=" << aTrack.GetKineticEnergy()/MeV
448  << " pos=" << aTrack.GetPosition().x()/m
449  << ", " << aTrack.GetPosition().y()/m
450  << ", " << aTrack.GetPosition().z()/m
451  <<G4endl;
452  }
453 #endif
454  }
455  // Check timing of secondaries
456  G4bool itsOKforTiming = true;
457 
458  accuracy = (theParentGlobalTime - aTrack.GetGlobalTime())/ns;
459  if (accuracy > accuracyForWarning){
460  itsOKforTiming = false;
461  nError += 1;
462  exitWithError = (accuracy > accuracyForException);
463 #ifdef G4VERBOSE
464  if (nError < maxError) {
465  G4cout << " G4VParticleChange::CheckSecondary : ";
466  G4cout << "the global time of secondary goes back comapared to the parent !!"
467  << " Difference: " << accuracy << "[ns] " <<G4endl;
468  G4cout << aTrack.GetDefinition()->GetParticleName()
469  << " E=" << aTrack.GetKineticEnergy()/MeV
470  << " pos=" << aTrack.GetPosition().x()/m
471  << ", " << aTrack.GetPosition().y()/m
472  << ", " << aTrack.GetPosition().z()/m
473  << " time=" << aTrack.GetGlobalTime()/ns
474  << " parent time=" << theParentGlobalTime/ns
475  <<G4endl;
476  }
477 #endif
478  }
479 
480  // Exit with error
481  if (exitWithError) {
482  G4Exception("G4VParticleChange::CheckSecondary",
483  "TRACK001", EventMustBeAborted,
484  "Secondary with illegal energy/momentum ");
485  }
486 
487  G4bool itsOK = itsOKforMomentum && itsOKforEnergy && itsOKforTiming;
488 
489  //correction
490  if (!itsOKforMomentum) {
491  G4double vmag = (aTrack.GetMomentumDirection()).mag();
492  aTrack.SetMomentumDirection((1./vmag)*aTrack.GetMomentumDirection());
493  }
494  if (!itsOKforEnergy) {
495  aTrack.SetKineticEnergy(0.0);
496  }
497 
498  if (!itsOK) {
499  this->DumpInfo();
500 
501  }
502 
503 
504  return itsOK;
505 }
506 
507 
509 {
510  return accuracyForWarning;
511 }
512 
514 {
515  return accuracyForException;
516 }
517 
518 
520 //Obsolete methods for parent weight
523 {
524 }
525 
526 
528 {
529  return true;
530 }
void SetFirstStepFlag()
G4ParticleDefinition * GetDefinition() const
void SetStepLength(G4double value)
void SetElement(G4int anIndex, Type *anElement)
Definition: G4FastVector.hh:76
static const double MeV
Definition: G4SIunits.hh:193
void AddNonIonizingEnergyDeposit(G4double value)
virtual G4Step * UpdateStepForAlongStep(G4Step *Step)
G4Track * GetSecondary(G4int anIndex) const
G4double GetWeight() const
G4VParticleChange & operator=(const G4VParticleChange &right)
void SetLastStepFlag()
const G4ThreeVector & GetPosition() const
G4TrackFastVector * theListOfSecondaries
void SetWeight(G4double aValue)
virtual void DumpInfo() const
G4bool operator==(const G4VParticleChange &right) const
#define G4ThreadLocal
Definition: tls.hh:89
void ClearLastStepFlag()
void ClearFirstStepFlag()
int G4int
Definition: G4Types.hh:78
const G4String & GetParticleName() const
static const G4double accuracyForException
void SetWeight(G4double aValue)
G4StepPoint * GetPreStepPoint() const
virtual G4bool CheckIt(const G4Track &)
void SetParentWeightByProcess(G4bool)
G4double GetKineticEnergy() const
virtual ~G4VParticleChange()
G4bool operator!=(const G4VParticleChange &right) const
G4GLOB_DLL std::ostream G4cout
void SetControlFlag(G4SteppingControl StepControlFlag)
bool G4bool
Definition: G4Types.hh:79
G4SteppingControl theSteppingControlFlag
G4bool IsParentWeightSetByProcess() const
Definition: G4Step.hh:76
G4double GetGlobalTime() const
void AddSecondary(G4Track *aSecondary)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
const G4ThreeVector & GetMomentumDirection() const
G4Step * UpdateStepInfo(G4Step *Step)
G4FastVector< G4Track, G4TrackFastVectorSize > G4TrackFastVector
G4StepPoint * GetPostStepPoint() const
void AddTotalEnergyDeposit(G4double value)
Definition: Step.hh:41
G4double GetAccuracyForException() const
#define G4endl
Definition: G4ios.hh:61
const G4int G4TrackFastVectorSize
static const double m
Definition: G4SIunits.hh:110
G4TrackStatus theStatusChange
void SetKineticEnergy(const G4double aValue)
double G4double
Definition: G4Types.hh:76
virtual G4Step * UpdateStepForPostStep(G4Step *Step)
static const G4double accuracyForWarning
G4double GetAccuracyForWarning() const
static const double mm
Definition: G4SIunits.hh:102
#define ns
Definition: xmlparse.cc:597
virtual G4Step * UpdateStepForAtRest(G4Step *Step)
G4bool CheckSecondary(G4Track &)
G4double theNonIonizingEnergyDeposit
void SetMomentumDirection(const G4ThreeVector &aValue)