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