Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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$
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 
45 const G4double G4VParticleChange::accuracyForWarning = 1.0e-9;
46 const G4double G4VParticleChange::accuracyForException = 0.001;
47 
48 G4VParticleChange::G4VParticleChange()
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
70  theListOfSecondaries = new G4TrackFastVector();
71 }
72 
73 G4VParticleChange::~G4VParticleChange() {
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 
89 G4VParticleChange::G4VParticleChange(const G4VParticleChange &right)
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 {
107  theListOfSecondaries = new G4TrackFastVector();
108  theNumberOfSecondaries = right.theNumberOfSecondaries;
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 
116 G4VParticleChange & G4VParticleChange::operator=(const G4VParticleChange &right)
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 
132  theListOfSecondaries = new G4TrackFastVector();
133  theNumberOfSecondaries = right.theNumberOfSecondaries;
134  for (G4int index = 0; index<theNumberOfSecondaries; index++){
135  G4Track* newTrack = new G4Track(*((*right.theListOfSecondaries)[index] ));
136  theListOfSecondaries->SetElement(index, newTrack);
137  }
138  theStatusChange = right.theStatusChange;
139  theSteppingControlFlag = right.theSteppingControlFlag;
140  theLocalEnergyDeposit = right.theLocalEnergyDeposit;
141  theNonIonizingEnergyDeposit = right.theNonIonizingEnergyDeposit;
142  theTrueStepLength = right.theTrueStepLength;
143 
144  theFirstStepInVolume = right.theFirstStepInVolume;
145  theLastStepInVolume = right.theLastStepInVolume;
146 
147  theParentWeight = right.theParentWeight;
148  isParentWeightProposed = right.isParentWeightProposed;
149  fSetSecondaryWeightByProcess = right.fSetSecondaryWeightByProcess;
150 
151  theParentGlobalTime = right.theParentGlobalTime;
152 
153  verboseLevel = right.verboseLevel;
154  debugFlag = right.debugFlag;
155 
156  }
157  return *this;
158 }
159 
160 G4bool G4VParticleChange::operator==(const G4VParticleChange &right) const
161 {
162  return (this == (G4VParticleChange *) &right);
163 }
164 
165 
166 G4bool G4VParticleChange::operator!=(const G4VParticleChange &right) const
167 {
168  return (this != (G4VParticleChange *) &right);
169 }
170 
171 void G4VParticleChange::AddSecondary(G4Track *aTrack)
172 {
173  if (debugFlag) CheckSecondary(*aTrack);
174 
175  // add a secondary after size check
176  if (theSizeOftheListOfSecondaries > theNumberOfSecondaries) {
177  // Set weight of secondary tracks
178  if (!fSetSecondaryWeightByProcess) aTrack->SetWeight(theParentWeight);
179  theListOfSecondaries->SetElement(theNumberOfSecondaries, aTrack);
180  theNumberOfSecondaries++;
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 
202 G4Step* G4VParticleChange::UpdateStepInfo(G4Step* pStep)
203 {
204  // Update the G4Step specific attributes
205  pStep->SetStepLength( theTrueStepLength );
206  pStep->AddTotalEnergyDeposit( theLocalEnergyDeposit );
207  pStep->AddNonIonizingEnergyDeposit( theNonIonizingEnergyDeposit );
208  pStep->SetControlFlag( theSteppingControlFlag );
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 
219 G4Step* G4VParticleChange::UpdateStepForAtRest(G4Step* Step)
220 {
221  if (isParentWeightProposed ){
222  Step->GetPostStepPoint()->SetWeight( theParentWeight );
223  }
224  return UpdateStepInfo(Step);
225 }
226 
227 
228 G4Step* G4VParticleChange::UpdateStepForAlongStep(G4Step* Step)
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 
239 G4Step* G4VParticleChange::UpdateStepForPostStep(G4Step* Step)
240 {
241  if (isParentWeightProposed ){
242  Step->GetPostStepPoint()->SetWeight( theParentWeight );
243  }
244  return UpdateStepInfo(Step);
245 }
246 
247 
248 //----------------------------------------------------------------
249 // methods for printing messages
250 //
251 
252 void G4VParticleChange::DumpInfo() const
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";
293  } else if( theStatusChange == fKillTrackAndSecondaries ){
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 
316 G4bool G4VParticleChange::CheckIt(const G4Track&
317 #ifdef G4VERBOSE
318  aTrack
319 #endif
320 )
321 {
322 
323  G4bool exitWithError = false;
324  G4double accuracy;
325  static 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 
398 G4bool G4VParticleChange::CheckSecondary(G4Track& aTrack)
399 {
400  G4bool exitWithError = false;
401  G4double accuracy;
402  static 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  return itsOK;
499 }
500 
501 
502 G4double G4VParticleChange::GetAccuracyForWarning() const
503 {
504  return accuracyForWarning;
505 }
506 
507 G4double G4VParticleChange::GetAccuracyForException() const
508 {
509  return accuracyForException;
510 }
511 
512 
514 //Obsolete methods for parent weight
516 void G4VParticleChange::SetParentWeightByProcess(G4bool )
517 {
518 }
519 
520 
521 G4bool G4VParticleChange::IsParentWeightSetByProcess() const
522 {
523  return true;
524 }
Generated on Sat May 25 2013 14:34:12 for Geant4 by   doxygen 1.8.4