Geant4  10.01.p02
G4RadioactiveDecay.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 // MODULE: G4RadioactiveDecay.cc
27 //
28 // Author: F Lei & P R Truscott
29 // Organisation: DERA UK
30 // Customer: ESA/ESTEC, NOORDWIJK
31 // Contract: 12115/96/JG/NL Work Order No. 3
32 //
33 // Documentation avaialable at http://www.space.dera.gov.uk/space_env/rdm.html
34 // These include:
35 // User Requirement Document (URD)
36 // Software Specification Documents (SSD)
37 // Software User Manual (SUM)
38 // Technical Note (TN) on the physics and algorithms
39 //
40 // The test and example programs are not included in the public release of
41 // G4 but they can be downloaded from
42 // http://www.space.qinetiq.com/space_env/rdm.html
43 //
44 // CHANGE HISTORY
45 // --------------
46 // 06 Aug 2014, L.G. Sarmiento Proton decay mode added mimicking the alpha decay
47 //
48 // 03 Oct 2012, V. Ivanchenko removed internal table for mean free path
49 // similar to what is done for as G4Decay
50 // 10 July 2012, L. Desorgher
51 // -In LoadDecayTable:
52 // Add LoadedNuclei.push_back(theParentNucleus.GetParticleName());
53 // also for the case where user data files are used. Correction for bug
54 // 1324. Changes proposed by Joa L.
55 //
56 //
57 // 01 May 2012, L. Desorgher
58 // -Force the reading of user file to theIsotopeTable
59 // -Merge the development by Fan Lei for activation computation
60 //
61 // 17 Oct 2011, L. Desorgher
62 // -Add possibility for the user to load its own decay file.
63 // -Set halflifethreshold negative by default to allow the tracking of all
64 // excited nuclei resulting from a radioactive decay
65 //
66 // 01 June 2011, M. Kelsey -- Add directional biasing interface to allow for
67 // "collimation" of decay daughters.
68 // 16 February 2006, V.Ivanchenko fix problem in IsApplicable connected with
69 // 8.0 particle design
70 // 18 October 2002, F. Lei
71 // in the case of beta decay, added a check of the end-energy
72 // to ensure it is > 0.
73 // ENSDF occationally have beta decay entries with zero energies
74 //
75 // 27 Sepetember 2001, F. Lei
76 // verboselevel(0) used in constructor
77 //
78 // 01 November 2000, F.Lei
79 // added " ee = e0 +1. ;" as line 763
80 // tagged as "radiative_decay-V02-00-02"
81 // 28 October 2000, F Lei
82 // added fast beta decay mode. Many files have been changed.
83 // tagged as "radiative_decay-V02-00-01"
84 //
85 // 25 October 2000, F Lei, DERA UK
86 // 1) line 1185 added 'const' to work with tag "Track-V02-00-00"
87 // tagged as "radiative_decay-V02-00-00"
88 // 14 April 2000, F Lei, DERA UK
89 // 0.b.4 release. Changes are:
90 // 1) Use PhotonEvaporation instead of DiscreteGammaDeexcitation
91 // 2) VR: Significant efficiency improvement
92 //
93 // 29 February 2000, P R Truscott, DERA UK
94 // 0.b.3 release.
95 //
97 //
98 #include "G4RadioactiveDecay.hh"
99 #include "G4RadioactiveDecaymessenger.hh"
100 
101 #include "G4SystemOfUnits.hh"
102 #include "G4DynamicParticle.hh"
103 #include "G4DecayProducts.hh"
104 #include "G4DecayTable.hh"
105 #include "G4ParticleChangeForRadDecay.hh"
106 #include "G4ITDecayChannel.hh"
107 #include "G4BetaMinusDecayChannel.hh"
108 #include "G4BetaPlusDecayChannel.hh"
109 #include "G4KshellECDecayChannel.hh"
110 #include "G4LshellECDecayChannel.hh"
111 #include "G4MshellECDecayChannel.hh"
112 // #include "G4AlphaDecayChannel.hh"
113 #include "G4AlphaDecay.hh"
114 #include "G4ProtonDecayChannel.hh"
115 #include "G4VDecayChannel.hh"
116 #include "G4RadioactiveDecayMode.hh"
117 #include "G4Ions.hh"
118 #include "G4IonTable.hh"
119 #include "G4RIsotopeTable.hh"
120 #include "G4BetaDecayType.hh"
121 #include "G4BetaDecayCorrections.hh"
122 #include "Randomize.hh"
123 #include "G4LogicalVolumeStore.hh"
124 #include "G4NuclearLevelManager.hh"
125 #include "G4NuclearLevelStore.hh"
126 #include "G4ThreeVector.hh"
127 #include "G4Electron.hh"
128 #include "G4Positron.hh"
129 #include "G4Neutron.hh"
130 #include "G4Gamma.hh"
131 #include "G4Alpha.hh"
132 #include "G4Proton.hh"
133 
134 #include "G4HadronicProcessType.hh"
135 #include "G4LossTableManager.hh"
136 #include "G4VAtomDeexcitation.hh"
137 
138 #include <vector>
139 #include <sstream>
140 #include <algorithm>
141 #include <fstream>
142 
143 using namespace CLHEP;
144 
145 const G4double G4RadioactiveDecay::levelTolerance = 2.0*keV;
146 const G4ThreeVector G4RadioactiveDecay::origin(0.,0.,0.);
147 #ifdef G4MULTITHREADED
148 DecayTableMap* G4RadioactiveDecay::master_dkmap = 0;
149 #endif
150 
151 G4RadioactiveDecay::G4RadioactiveDecay(const G4String& processName)
152  : G4VRestDiscreteProcess(processName, fDecay), isInitialised(false),
153  forceDecayDirection(0.,0.,0.), forceDecayHalfAngle(0.*deg), verboseLevel(0)
154 {
155 #ifdef G4VERBOSE
156  if (GetVerboseLevel() > 1) {
157  G4cout << "G4RadioactiveDecay constructor: processName = " << processName
158  << G4endl;
159  }
160 #endif
161 
162  SetProcessSubType(fRadioactiveDecay);
163 
164  theRadioactiveDecaymessenger = new G4RadioactiveDecaymessenger(this);
165  pParticleChange = &fParticleChangeForRadDecay;
166  theIsotopeTable = new G4RIsotopeTable();
167 
168  // Regsiter the isotope table to the ion table.
169  // Although we are touching the ion table, which is shared, we are not
170  // adding particles in this operation. We can therefore do the registration
171  // for each instance of the RDM process and do not need to restrict it to the
172  // master process. It's possible that future, more optimized versions of
173  // G4IonTable will require to test for master.
174 
175  G4IonTable* theIonTable = G4ParticleTable::GetParticleTable()->GetIonTable();
176  theIonTable->RegisterIsotopeTable(theIsotopeTable);
177 
178  // Reset the list of user defined data files
179  theUserRadioactiveDataFiles.clear();
180 
181  // Instantiate the map of decay tables
182 #ifdef G4MULTITHREADED
183  if(!master_dkmap) master_dkmap = new DecayTableMap;
184 #endif
185  dkmap = new DecayTableMap;
186 
187  // Apply default values.
188  NSourceBin = 1;
189  SBin[0] = 0.* s;
190  SBin[1] = 1.* s;
191  SProfile[0] = 1.;
192  SProfile[1] = 0.;
193  NDecayBin = 1;
194  DBin[0] = 0. * s ;
195  DBin[1] = 1. * s;
196  DProfile[0] = 1.;
197  DProfile[1] = 0.;
198  decayWindows[0] = 0;
199  G4RadioactivityTable* rTable = new G4RadioactivityTable() ;
200  theRadioactivityTables.push_back(rTable);
201  NSplit = 1;
202  AnalogueMC = true ;
203  FBeta = false ;
204  BRBias = true ;
205  applyICM = true ;
206  applyARM = true ;
207  halflifethreshold = nanosecond;
208 
209  // RDM applies to xall logical volumes as default
210  isAllVolumesMode = true;
211  SelectAllVolumes();
212 }
213 
214 
215 G4RadioactiveDecay::~G4RadioactiveDecay()
216 {
217  delete theRadioactiveDecaymessenger;
218  delete theIsotopeTable;
219  for (DecayTableMap::iterator i = dkmap->begin(); i != dkmap->end(); i++) {
220  delete i->second;
221  }
222  dkmap->clear();
223  delete dkmap;
224 }
225 
226 
227 G4bool G4RadioactiveDecay::IsApplicable(const G4ParticleDefinition& aParticle)
228 {
229  // All particles other than G4Ions, are rejected by default
230  if (((const G4Ions*)(&aParticle))->GetExcitationEnergy() > 0.) {return true;}
231  if (aParticle.GetParticleName() == "GenericIon") {
232  return true;
233  } else if (!(aParticle.GetParticleType() == "nucleus")
234  || aParticle.GetPDGLifeTime() < 0. ) {
235  return false;
236  }
237 
238  // Determine whether the nuclide falls into the correct A and Z range
239  G4int A = ((const G4Ions*) (&aParticle))->GetAtomicMass();
240  G4int Z = ((const G4Ions*) (&aParticle))->GetAtomicNumber();
241  if (A > theNucleusLimits.GetAMax() || A < theNucleusLimits.GetAMin())
242  {return false;}
243  else if (Z > theNucleusLimits.GetZMax() || Z < theNucleusLimits.GetZMin())
244  {return false;}
245  return true;
246 }
247 
248 G4DecayTable* G4RadioactiveDecay::GetDecayTable(const G4ParticleDefinition* aNucleus)
249 {
250  G4String key = aNucleus->GetParticleName();
251  DecayTableMap::iterator table_ptr = dkmap->find(key);
252 
253  G4DecayTable* theDecayTable = 0;
254  if (table_ptr == dkmap->end() ) { // If table not there,
255  theDecayTable = LoadDecayTable(*aNucleus); // load from file and
256  (*dkmap)[key] = theDecayTable; // store in library
257  } else {
258  theDecayTable = table_ptr->second;
259  }
260 
261  return theDecayTable;
262 }
263 
264 
265 void G4RadioactiveDecay::SelectAVolume(const G4String aVolume)
266 {
267  G4LogicalVolumeStore* theLogicalVolumes;
268  G4LogicalVolume* volume;
269  theLogicalVolumes = G4LogicalVolumeStore::GetInstance();
270  for (size_t i = 0; i < theLogicalVolumes->size(); i++) {
271  volume=(*theLogicalVolumes)[i];
272  if (volume->GetName() == aVolume) {
273  ValidVolumes.push_back(aVolume);
274  std::sort(ValidVolumes.begin(), ValidVolumes.end());
275  // sort need for performing binary_search
276 #ifdef G4VERBOSE
277  if (GetVerboseLevel()>0)
278  G4cout << " RDM Applies to : " << aVolume << G4endl;
279 #endif
280  } else if(i == theLogicalVolumes->size()) {
281  G4cerr << "SelectAVolume: "<< aVolume
282  << " is not a valid logical volume name" << G4endl;
283  }
284  }
285 }
286 
287 
288 void G4RadioactiveDecay::DeselectAVolume(const G4String aVolume)
289 {
290  G4LogicalVolumeStore* theLogicalVolumes;
291  G4LogicalVolume* volume;
292  theLogicalVolumes = G4LogicalVolumeStore::GetInstance();
293  for (size_t i = 0; i < theLogicalVolumes->size(); i++){
294  volume=(*theLogicalVolumes)[i];
295  if (volume->GetName() == aVolume) {
296  std::vector<G4String>::iterator location;
297  location = std::find(ValidVolumes.begin(),ValidVolumes.end(),aVolume);
298  if (location != ValidVolumes.end()) {
299  ValidVolumes.erase(location);
300  std::sort(ValidVolumes.begin(), ValidVolumes.end());
301  isAllVolumesMode =false;
302  } else {
303  G4cerr << " DeselectVolume:" << aVolume << " is not in the list "
304  << G4endl;
305  }
306 #ifdef G4VERBOSE
307  if (GetVerboseLevel() > 0)
308  G4cout << " DeselectVolume: " << aVolume << " is removed from list "
309  << G4endl;
310 #endif
311  } else if (i == theLogicalVolumes->size()) {
312  G4cerr << " DeselectVolume:" << aVolume
313  << "is not a valid logical volume name" << G4endl;
314  }
315  }
316 }
317 
318 
319 void G4RadioactiveDecay::SelectAllVolumes()
320 {
321  G4LogicalVolumeStore* theLogicalVolumes;
322  G4LogicalVolume* volume;
323  theLogicalVolumes = G4LogicalVolumeStore::GetInstance();
324  ValidVolumes.clear();
325 #ifdef G4VERBOSE
326  if (GetVerboseLevel()>0)
327  G4cout << " RDM Applies to all Volumes" << G4endl;
328 #endif
329  for (size_t i = 0; i < theLogicalVolumes->size(); i++){
330  volume = (*theLogicalVolumes)[i];
331  ValidVolumes.push_back(volume->GetName());
332 #ifdef G4VERBOSE
333  if (GetVerboseLevel()>0)
334  G4cout << " RDM Applies to Volume " << volume->GetName() << G4endl;
335 #endif
336  }
337  std::sort(ValidVolumes.begin(), ValidVolumes.end());
338  // sort needed in order to allow binary_search
339  isAllVolumesMode=true;
340 }
341 
342 
343 void G4RadioactiveDecay::DeselectAllVolumes()
344 {
345  ValidVolumes.clear();
346  isAllVolumesMode=false;
347 #ifdef G4VERBOSE
348  if (GetVerboseLevel() > 0) G4cout << "RDM removed from all volumes" << G4endl;
349 #endif
350 }
351 
352 
353 G4bool
354 G4RadioactiveDecay::IsRateTableReady(const G4ParticleDefinition& aParticle)
355 {
356  // Check whether the radioactive decay rates table for the ion has already
357  // been calculated.
358  G4String aParticleName = aParticle.GetParticleName();
359  for (size_t i = 0; i < theDecayRateTableVector.size(); i++) {
360  if (theDecayRateTableVector[i].GetIonName() == aParticleName) return true;
361  }
362  return false;
363 }
364 
365 // GetDecayRateTable
366 // retrieve the decayratetable for the specified aParticle
367 
368 void
369 G4RadioactiveDecay::GetDecayRateTable(const G4ParticleDefinition& aParticle)
370 {
371  G4String aParticleName = aParticle.GetParticleName();
372 
373  for (size_t i = 0; i < theDecayRateTableVector.size(); i++) {
374  if (theDecayRateTableVector[i].GetIonName() == aParticleName) {
375  theDecayRateVector = theDecayRateTableVector[i].GetItsRates();
376  }
377  }
378 #ifdef G4VERBOSE
379  if (GetVerboseLevel() > 0) {
380  G4cout << "The DecayRate Table for " << aParticleName << " is selected."
381  << G4endl;
382  }
383 #endif
384 }
385 
386 // GetTaoTime performs the convolution of the source time profile function
387 // with the decay constants in the decay chain.
388 G4double G4RadioactiveDecay::GetTaoTime(const G4double t, const G4double tao)
389 {
390  long double taotime = 0.L;
391  G4int nbin;
392  if ( t > SBin[NSourceBin]) {
393  nbin = NSourceBin;}
394  else {
395  nbin = 0;
396  while (t > SBin[nbin]) nbin++;
397  nbin--;}
398  long double lt = t ;
399  long double ltao = tao;
400 
401  if (nbin > 0) {
402  for (G4int i = 0; i < nbin; i++) {
403  taotime += (long double)SProfile[i] *
404  (std::exp(-(lt-(long double)SBin[i+1])/ltao)-std::exp(-(lt-(long double)SBin[i])/ltao));
405  }
406  }
407  taotime += (long double)SProfile[nbin] * (1.L-std::exp(-(lt-(long double)SBin[nbin])/ltao));
408  if (taotime < 0.) {
409  G4cout <<" Tao time =: " <<taotime << " reset to zero!"<<G4endl;
410  G4cout <<" t = " << t <<" tao = " <<tao <<G4endl;
411  G4cout << SBin[nbin] << " " <<SBin[0] << G4endl;
412  taotime = 0.;
413  }
414 #ifdef G4VERBOSE
415  if (GetVerboseLevel()>1)
416  {G4cout <<" Tao time: " <<taotime <<G4endl;}
417 #endif
418  return (G4double)taotime ;
419 }
420 
421 /*
422 // Other implementation tests to avoid use of long double
423 G4double G4RadioactiveDecay::GetTaoTime(const G4double t, const G4double tao)
424 {
425  long double taotime =0.L;
426  G4int nbin;
427  if ( t > SBin[NSourceBin]) {
428  nbin = NSourceBin;}
429  else {
430  nbin = 0;
431  while (t > SBin[nbin]) nbin++;
432  nbin--;}
433  long double lt = t ;
434  long double ltao = tao;
435  long double factor,factor1,dt1,dt;
436  if (nbin > 0) {
437  for (G4int i = 0; i < nbin; i++)
438  { long double s1=SBin[i];
439  long double s2=SBin[i+1];
440  dt1=(s2-s1)/ltao;
441  if (dt1 <50.) {
442  factor1=std::exp(dt1)-1.;
443  if (factor1<dt1) factor1 =dt1;
444  dt=(lt-s1)/ltao;
445  factor=std::exp(-dt);
446  }
447  else {
448  factor1=1.-std::exp(-dt1);
449  dt=(lt-s2)/ltao;
450  factor=std::exp(-dt);
451  }
452  G4cout<<(long double) SProfile[i] *factor*factor1<<'\t'<<std::endl;
453  long double test = (long double)SProfile[i] * (std::exp(-(lt-(long double)SBin[i+1])/ltao)-std::exp(-(lt-(long double)SBin[i])/ltao));
454  G4cout<<test<<std::endl;
455  taotime += (long double) SProfile[i] *factor*factor1;
456  }
457  }
458  long double s=SBin[nbin];
459  dt1=(lt-s)/ltao;
460  factor=1.-std::exp(-dt1);
461  taotime += (long double) SProfile[nbin] *factor;
462  if (taotime < 0.) {
463  G4cout <<" Tao time =: " <<taotime << " reset to zero!"<<G4endl;
464  G4cout <<" t = " << t <<" tao = " <<tao <<G4endl;
465  G4cout << SBin[nbin] << " " <<SBin[0] << G4endl;
466  taotime = 0.;
467  }
468 #ifdef G4VERBOSE
469  if (GetVerboseLevel()>1)
470  {G4cout <<" Tao time: " <<taotime <<G4endl;}
471 #endif
472  return (G4double)taotime ;
473 }
474 
475 
476 
477 G4double G4RadioactiveDecay::GetTaoTime(const G4double t, const G4double tao)
478 {
479  G4double taotime =0.;
480  G4int nbin;
481  if ( t > SBin[NSourceBin]) {
482  nbin = NSourceBin;}
483  else {
484  nbin = 0;
485  while (t > SBin[nbin]) nbin++;
486  nbin--;}
487  G4double lt = t ;
488  G4double ltao = tao;
489  G4double factor,factor1,dt1,dt;
490 
491  if (nbin > 0) {
492  for (G4int i = 0; i < nbin; i++)
493  { dt1=(SBin[i+1]-SBin[i])/ltao;
494  if (dt1 <50.) {
495  factor1=std::exp(dt1)-1.;
496  if (factor1<dt1) factor1 =dt1;
497  dt=(lt-SBin[i])/ltao;
498  factor=std::exp(-(lt-SBin[i])/ltao);
499  G4cout<<factor<<'\t'<<factor1<<std::endl;
500  }
501  else {
502  factor1=1.-std::exp(-dt1);
503  factor=std::exp(-(lt-SBin[i+1])/ltao);
504  }
505  G4cout<<factor<<'\t'<<factor1<<std::endl;
506  taotime += SProfile[i] *factor*factor1;
507  G4cout<<taotime<<std::endl;
508  }
509  }
510  dt1=(lt-SBin[nbin])/ltao;
511  factor=1.-std::exp(-dt1);
512  if (factor<(dt1-0.5*dt1*dt1)) factor =dt1-0.5*dt1*dt1;
513 
514 
515  taotime += SProfile[nbin] *factor;
516  G4cout<<factor<<'\t'<<taotime<<std::endl;
517  if (taotime < 0.) {
518  G4cout <<" Tao time =: " <<taotime << " reset to zero!"<<G4endl;
519  G4cout <<" t = " << t <<" tao = " <<tao <<G4endl;
520  G4cout << SBin[nbin] << " " <<SBin[0] << G4endl;
521  taotime = 0.;
522  }
523 
524 #ifdef G4VERBOSE
525  if (GetVerboseLevel()>1)
526  {G4cout <<" Tao time: " <<taotime <<G4endl;}
527 #endif
528  return (G4double)taotime ;
529 }
530 */
532 // //
533 // GetDecayTime //
534 // Randomly select a decay time for the decay process, following the //
535 // supplied decay time bias scheme. //
536 // //
538 
539 G4double G4RadioactiveDecay::GetDecayTime()
540 {
541  G4double decaytime = 0.;
542  G4double rand = G4UniformRand();
543  G4int i = 0;
544  while ( DProfile[i] < rand) i++;
545  rand = G4UniformRand();
546  decaytime = DBin[i] + rand*(DBin[i+1]-DBin[i]);
547 #ifdef G4VERBOSE
548  if (GetVerboseLevel()>1)
549  {G4cout <<" Decay time: " <<decaytime/s <<"[s]" <<G4endl;}
550 #endif
551  return decaytime;
552 }
553 
554 
555 G4int G4RadioactiveDecay::GetDecayTimeBin(const G4double aDecayTime)
556 {
557  G4int i = 0;
558  while ( aDecayTime > DBin[i] ) i++;
559  return i;
560 }
561 
563 // //
564 // GetMeanLifeTime (required by the base class) //
565 // //
567 
568 G4double G4RadioactiveDecay::GetMeanLifeTime(const G4Track& theTrack,
569  G4ForceCondition*)
570 {
571  // For varience reduction implementation the time is set to 0 so as to
572  // force the particle to decay immediately.
573  // In analogueMC mode it return the particle's mean-life.
574 
575  G4double meanlife = 0.;
576  if (AnalogueMC) {
577  const G4DynamicParticle* theParticle = theTrack.GetDynamicParticle();
578  const G4ParticleDefinition* theParticleDef = theParticle->GetDefinition();
579  G4double theLife = theParticleDef->GetPDGLifeTime();
580 
581 #ifdef G4VERBOSE
582  if (GetVerboseLevel() > 2) {
583  G4cout << "G4RadioactiveDecay::GetMeanLifeTime() " << G4endl;
584  G4cout << "KineticEnergy: " << theParticle->GetKineticEnergy()/GeV
585  << " GeV, Mass: " << theParticle->GetMass()/GeV
586  << " GeV, Life time: " << theLife/ns << " ns " << G4endl;
587  }
588 #endif
589  if (theParticleDef->GetPDGStable()) {meanlife = DBL_MAX;}
590  else if (theLife < 0.0) {meanlife = DBL_MAX;}
591  else {meanlife = theLife;}
592  // Set meanlife to zero for excited istopes which are not in the
593  // RDM database
594  if (((const G4Ions*)(theParticleDef))->GetExcitationEnergy() > 0. &&
595  meanlife == DBL_MAX) {meanlife = 0.;}
596  }
597 #ifdef G4VERBOSE
598  if (GetVerboseLevel() > 1)
599  G4cout << " mean life time: " << meanlife/s << " s " << G4endl;
600 #endif
601 
602  return meanlife;
603 }
604 
606 // //
607 // GetMeanFreePath for decay in flight //
608 // //
610 
611 G4double G4RadioactiveDecay::GetMeanFreePath (const G4Track& aTrack, G4double,
612  G4ForceCondition*)
613 {
614  const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
615  const G4ParticleDefinition* aParticleDef = aParticle->GetDefinition();
616  G4double tau = aParticleDef->GetPDGLifeTime();
617  G4double aMass = aParticle->GetMass();
618 
619 #ifdef G4VERBOSE
620  if (GetVerboseLevel() > 2) {
621  G4cout << "G4RadioactiveDecay::GetMeanFreePath() " << G4endl;
622  G4cout << " KineticEnergy: " << aParticle->GetKineticEnergy()/GeV
623  << " GeV, Mass: " << aMass/GeV << " GeV, tau: " << tau << " ns "
624  << G4endl;
625  }
626 #endif
627  G4double pathlength = DBL_MAX;
628  if (tau != -1) {
629  // Ion can decay
630 
631  if (tau < -1000.0) {
632  pathlength = DBL_MIN; // nuclide had very short lifetime or wasn't in table
633 
634  } else if (tau < 0.0) {
635  G4cout << aParticleDef->GetParticleName() << " has lifetime " << tau << G4endl;
636  G4ExceptionDescription ed;
637  ed << "Ion has negative lifetime " << tau
638  << " but is not stable. Setting mean free path to DBL_MAX" << G4endl;
639  G4Exception("G4RadioactiveDecay::GetMeanFreePath()", "HAD_RDM_011",
640  JustWarning, ed);
641  pathlength = DBL_MAX;
642 
643  } else {
644  // Calculate mean free path
645  G4double betaGamma = aParticle->GetTotalMomentum()/aMass;
646  pathlength = c_light*tau*betaGamma;
647 
648  if (pathlength < DBL_MIN) {
649  pathlength = DBL_MIN;
650 #ifdef G4VERBOSE
651  if (GetVerboseLevel() > 2) {
652  G4cout << "G4Decay::GetMeanFreePath: "
653  << aParticleDef->GetParticleName()
654  << " stops, kinetic energy = "
655  << aParticle->GetKineticEnergy()/keV <<" keV " << G4endl;
656  }
657 #endif
658  }
659  }
660  }
661 
662 #ifdef G4VERBOSE
663  if (GetVerboseLevel() > 1) {
664  G4cout << "mean free path: "<< pathlength/m << " m" << G4endl;
665  }
666 #endif
667  return pathlength;
668 }
669 
671 // //
672 // BuildPhysicsTable - initialisation of atomic de-excitation //
673 // //
675 
676 void G4RadioactiveDecay::BuildPhysicsTable(const G4ParticleDefinition&)
677 {
678  if (!isInitialised) {
679  isInitialised = true;
680  G4LossTableManager* theManager = G4LossTableManager::Instance();
681  G4VAtomDeexcitation* p = theManager->AtomDeexcitation();
682  if(p) { p->InitialiseAtomicDeexcitation(); }
683  }
684 }
685 
687 // //
688 // LoadDecayTable loads the decay scheme from the RadioactiveDecay database //
689 // for the parent nucleus. //
690 // //
692 
693 #ifdef G4MULTITHREADED
694 #include "G4AutoLock.hh"
695 G4Mutex G4RadioactiveDecay::radioactiveDecayMutex = G4MUTEX_INITIALIZER;
696 #endif
697 
698 G4DecayTable*
699 G4RadioactiveDecay::LoadDecayTable(const G4ParticleDefinition& theParentNucleus)
700 {
701  // Generate input data file name using Z and A of the parent nucleus
702  // file containing radioactive decay data.
703  G4int A = ((const G4Ions*)(&theParentNucleus))->GetAtomicMass();
704  G4int Z = ((const G4Ions*)(&theParentNucleus))->GetAtomicNumber();
705  G4int lvl = ((const G4Ions*)(&theParentNucleus))->GetIsomerLevel();
706  G4DecayTable* theDecayTable = 0;
707 
708 #ifdef G4MULTITHREADED
709  G4AutoLock lk(&G4RadioactiveDecay::radioactiveDecayMutex);
710 
711  G4String key = theParentNucleus.GetParticleName();
712  DecayTableMap::iterator master_table_ptr = master_dkmap->find(key);
713 
714  if (master_table_ptr != master_dkmap->end() ) { // If table is there
715  return master_table_ptr->second;
716  }
717 #endif
718 
719  // Create and initialise variables used in the method.
720  theDecayTable = new G4DecayTable();
721 
722  //Check if data have been provided by the user
723  G4String file= theUserRadioactiveDataFiles[1000*A+Z];
724 
725  if (file =="") {
726  if (!getenv("G4RADIOACTIVEDATA") ) {
727  G4cout << "Please setenv G4RADIOACTIVEDATA to point to the radioactive decay data files."
728  << G4endl;
729  throw G4HadronicException(__FILE__, __LINE__, " Please setenv G4RADIOACTIVEDATA to point to the radioactive decay data files.");
730  }
731  G4String dirName = getenv("G4RADIOACTIVEDATA");
732 
733  std::ostringstream os;
734  os <<dirName <<"/z" <<Z <<".a" <<A <<'\0';
735  file = os.str();
736  }
737 
738  std::ifstream DecaySchemeFile(file);
739 
740  G4bool found(false);
741  if (DecaySchemeFile) {
742  // Initialise variables used for reading in radioactive decay data.
743  G4int nMode = 8;
744  G4bool modeFirstRecord[8];
745  G4double modeTotalBR[8] = {0.0};
746  G4double modeSumBR[8];
747  for (G4int i = 0; i < nMode; i++) {
748  modeFirstRecord[i] = true;
749  modeSumBR[i] = 0.0;
750  }
751 
752  G4bool complete(false);
753  char inputChars[100]={' '};
754  G4String inputLine;
755  G4String recordType("");
756  G4RadioactiveDecayMode theDecayMode;
757  G4double a(0.0);
758  G4double b(0.0);
759  G4double c(0.0);
760  G4int levelCounter = 0;
761  G4BetaDecayType betaType(allowed);
762  G4double e0;
763 
764  // Loop through each data file record until you identify the decay
765  // data relating to the nuclide of concern.
766 
767  while (!complete && !DecaySchemeFile.getline(inputChars, 100).eof()) {
768  inputLine = inputChars;
769  inputLine = inputLine.strip(1);
770  if (inputChars[0] != '#' && inputLine.length() != 0) {
771  std::istringstream tmpStream(inputLine);
772 
773  if (inputChars[0] == 'P') {
774  // Nucleus is a parent type. Check excitation level to see if it
775  // matches that of theParentNucleus
776  tmpStream >> recordType >> a >> b;
777  if (found) {
778  complete = true;
779 // else {found = (std::abs(a*keV - E) < levelTolerance);}
780  } else {
781  found = (levelCounter == lvl);
782  }
783  levelCounter++;
784 
785  } else if (found) {
786  // The right part of the radioactive decay data file has been found. Search
787  // through it to determine the mode of decay of the subsequent records.
788  if (inputChars[0] == 'W') {
789 #ifdef G4VERBOSE
790  if (GetVerboseLevel() > 0) {
791  // a comment line identified and print out the message
792  G4cout << " Warning in G4RadioactiveDecay::LoadDecayTable " << G4endl;
793  G4cout << " In data file " << file << G4endl;
794  G4cout << " " << inputLine << G4endl;
795  }
796 #endif
797  } else {
798  tmpStream >> theDecayMode >> a >> b >> c >> betaType;
799 
800  // Allowed transitions are the default. Forbidden transitions are
801  // indicated in the last column.
802  if (inputLine.length() < 80) betaType = allowed;
803  a /= 1000.;
804  c /= 1000.;
805 
806  switch (theDecayMode) {
807 
808  case IT: // Isomeric transition
809  {
810  G4ITDecayChannel* anITChannel =
811  new G4ITDecayChannel(GetVerboseLevel(),
812  (const G4Ions*)& theParentNucleus, b);
813  anITChannel->SetICM(applyICM);
814  anITChannel->SetARM(applyARM);
815  anITChannel->SetHLThreshold(halflifethreshold);
816  theDecayTable->Insert(anITChannel);
817  }
818  break;
819 
820  case BetaMinus:
821  {
822  if (modeFirstRecord[1]) {
823  modeFirstRecord[1] = false;
824  modeTotalBR[1] = b;
825  } else {
826  if (c > 0.) {
827  e0 = c*MeV/0.511;
828  G4BetaDecayCorrections corrections(Z+1, A);
829 
830  // array to store function shape
831  G4int npti = 100;
832  G4double* pdf = new G4double[npti];
833 
834  G4double e; // Total electron energy in units of electron mass
835  G4double p; // Electron momentum in units of electron mass
836  G4double f; // Spectral shape function value
837  for (G4int ptn = 0; ptn < npti; ptn++) {
838  // Calculate simple phase space spectrum
839  e = 1. + e0*(ptn+0.5)/100.;
840  p = std::sqrt(e*e - 1.);
841  f = p*e*(e0-e+1)*(e0-e+1);
842 
843  // Apply Fermi factor to get allowed shape
844  f *= corrections.FermiFunction(e);
845 
846  // Apply shape factor for forbidden transitions
847  f *= corrections.ShapeFactor(betaType, p, e0-e+1.);
848  pdf[ptn] = f;
849  }
850 
851  G4RandGeneral* aRandomEnergy = new G4RandGeneral( pdf, npti);
852  G4BetaMinusDecayChannel *aBetaMinusChannel = new
853  G4BetaMinusDecayChannel(GetVerboseLevel(), &theParentNucleus,
854  b, c*MeV, a*MeV, 0, FBeta, aRandomEnergy);
855  aBetaMinusChannel->SetICM(applyICM);
856  aBetaMinusChannel->SetARM(applyARM);
857  aBetaMinusChannel->SetHLThreshold(halflifethreshold);
858  theDecayTable->Insert(aBetaMinusChannel);
859  modeSumBR[1] += b;
860  delete[] pdf;
861  } // c > 0
862  } // if not first record
863  }
864  break;
865 
866  case BetaPlus:
867  {
868  if (modeFirstRecord[2]) {
869  modeFirstRecord[2] = false;
870  modeTotalBR[2] = b;
871  } else {
872  e0 = c*MeV/0.510999 - 2.;
873  // Need to test e0 for nuclei which have Q < 2Me in their
874  // data files (e.g. z67.a162)
875  if (e0 > 0.) {
876  G4BetaDecayCorrections corrections(1-Z, A);
877 
878  // array to store function shape
879  G4int npti = 100;
880  G4double* pdf = new G4double[npti];
881 
882  G4double e; // Total positron energy in units of electron mass
883  G4double p; // Positron momentum in units of electron mass
884  G4double f; // Spectral shape function value
885  for (G4int ptn = 0; ptn < npti; ptn++) {
886  // Calculate simple phase space spectrum
887  e = 1. + e0*(ptn+0.5)/100.;
888  p = std::sqrt(e*e - 1.);
889  f = p*e*(e0-e+1)*(e0-e+1);
890 
891  // Apply Fermi factor to get allowed shape
892  f *= corrections.FermiFunction(e);
893 
894  // Apply shape factor for forbidden transitions
895  f *= corrections.ShapeFactor(betaType, p, e0-e+1.);
896  pdf[ptn] = f;
897  }
898  G4RandGeneral* aRandomEnergy = new G4RandGeneral(pdf, npti);
899  G4BetaPlusDecayChannel* aBetaPlusChannel = new
900  G4BetaPlusDecayChannel(GetVerboseLevel(),
901  &theParentNucleus, b,
902  (c-1.021998)*MeV, a*MeV, 0,
903  FBeta, aRandomEnergy);
904  aBetaPlusChannel->SetICM(applyICM);
905  aBetaPlusChannel->SetARM(applyARM);
906  aBetaPlusChannel->SetHLThreshold(halflifethreshold);
907  theDecayTable->Insert(aBetaPlusChannel);
908  modeSumBR[2] += b;
909  delete[] pdf;
910  } // if e0 > 0
911  } // if not first record
912  }
913  break;
914 
915  case KshellEC: // K-shell electron capture
916 
917  if (modeFirstRecord[3]) {
918  modeFirstRecord[3] = false;
919  modeTotalBR[3] = b;
920  } else {
921  G4KshellECDecayChannel* aKECChannel =
922  new G4KshellECDecayChannel(GetVerboseLevel(),
923  &theParentNucleus,
924  b, c*MeV, a*MeV);
925  aKECChannel->SetICM(applyICM);
926  aKECChannel->SetARM(applyARM);
927  aKECChannel->SetHLThreshold(halflifethreshold);
928  theDecayTable->Insert(aKECChannel);
929  modeSumBR[3] += b;
930  }
931  break;
932 
933  case LshellEC: // L-shell electron capture
934 
935  if (modeFirstRecord[4]) {
936  modeFirstRecord[4] = false;
937  modeTotalBR[4] = b;
938  } else {
939  G4LshellECDecayChannel *aLECChannel = new
940  G4LshellECDecayChannel (GetVerboseLevel(), &theParentNucleus,
941  b, c*MeV, a*MeV);
942  aLECChannel->SetICM(applyICM);
943  aLECChannel->SetARM(applyARM);
944  aLECChannel->SetHLThreshold(halflifethreshold);
945  theDecayTable->Insert(aLECChannel);
946  modeSumBR[4] += b;
947  }
948  break;
949 
950  case MshellEC: // M-shell electron capture
951  // In this implementation it is added to L-shell case
952  if (modeFirstRecord[5]) {
953  modeFirstRecord[5] = false;
954  modeTotalBR[5] = b;
955  } else {
956  G4MshellECDecayChannel* aMECChannel =
957  new G4MshellECDecayChannel(GetVerboseLevel(),
958  &theParentNucleus,
959  b, c*MeV, a*MeV);
960  aMECChannel->SetICM(applyICM);
961  aMECChannel->SetARM(applyARM);
962  aMECChannel->SetHLThreshold(halflifethreshold);
963  theDecayTable->Insert(aMECChannel);
964  modeSumBR[5] += b;
965  }
966  break;
967 
968  case Alpha:
969  if (modeFirstRecord[6]) {
970  modeFirstRecord[6] = false;
971  modeTotalBR[6] = b;
972  } else {
973  G4AlphaDecay* anAlphaChannel =
974  new G4AlphaDecay(&theParentNucleus, b, c*MeV, a*MeV);
975 // anAlphaChannel->DumpInfo();
976 /*
977  G4AlphaDecayChannel* anAlphaChannel =
978  new G4AlphaDecayChannel(GetVerboseLevel(),
979  &theParentNucleus,
980  b, c*MeV, a*MeV);
981  anAlphaChannel->SetICM(applyICM);
982  anAlphaChannel->SetARM(applyARM);
983 */
984  anAlphaChannel->SetHLThreshold(halflifethreshold);
985  theDecayTable->Insert(anAlphaChannel);
986  modeSumBR[6] += b;
987  }
988  break;
989 
990  case Proton:
991  if (modeFirstRecord[7]) {
992  modeFirstRecord[7] = false;
993  modeTotalBR[7] = b;
994  } else {
995  G4ProtonDecayChannel* aProtonChannel =
996  new G4ProtonDecayChannel(GetVerboseLevel(),
997  &theParentNucleus,
998  b, c*MeV, a*MeV);
999  aProtonChannel->SetICM(applyICM);
1000  aProtonChannel->SetARM(applyARM);
1001  aProtonChannel->SetHLThreshold(halflifethreshold);
1002  theDecayTable->Insert(aProtonChannel);
1003  modeSumBR[7] += b;
1004  }
1005  break;
1006 
1007  case Beta2Minus:
1008  // Not yet implemented
1009  // G4cout << " Double beta- decay, a = " << a << ", b = " << b << ", c = " << c << G4endl;
1010  break;
1011 
1012  case Beta2Plus:
1013  // Not yet implemented
1014  // G4cout << " Double beta+ decay, a = " << a << ", b = " << b << ", c = " << c << G4endl;
1015  break;
1016 
1017  case Proton2:
1018  // Not yet implemented
1019  // G4cout << " Double proton decay, a = " << a << ", b = " << b << ", c = " << c << G4endl;
1020  break;
1021 
1022  case SpFission:
1023  // Not yet implemented
1024  //G4cout<<"Sp fission channel"<<a<<'\t'<<b<<'\t'<<c<<std::endl;
1025  break;
1026  case RDM_ERROR:
1027 
1028  default:
1029  G4Exception("G4RadioactiveDecay::LoadDecayTable()", "HAD_RDM_000",
1030  FatalException, "Selected decay mode does not exist");
1031  } // switch
1032  } // if char == W
1033  } // if char == P
1034  } // if char != #
1035  } // While
1036 
1037  // Go through the decay table and make sure that the branching ratios are
1038  // correctly normalised.
1039 
1040  G4VDecayChannel* theChannel = 0;
1041  G4NuclearDecayChannel* theNuclearDecayChannel = 0;
1042  G4String mode = "";
1043 
1044  G4double theBR = 0.0;
1045  for (G4int i = 0; i < theDecayTable->entries(); i++) {
1046  theChannel = theDecayTable->GetDecayChannel(i);
1047  if (theChannel->GetKinematicsName() == "alpha decay") {
1048  theDecayMode = Alpha;
1049  } else {
1050  theNuclearDecayChannel = static_cast<G4NuclearDecayChannel*>(theChannel);
1051  theDecayMode = theNuclearDecayChannel->GetDecayMode();
1052  }
1053  if (theDecayMode != IT) {
1054  theBR = theChannel->GetBR();
1055  theChannel->SetBR(theBR*modeTotalBR[theDecayMode]/modeSumBR[theDecayMode]);
1056  }
1057  }
1058  } // if (DecaySchemeFile)
1059  DecaySchemeFile.close();
1060 
1061  if (!found && lvl > 0) {
1062  // Case where IT cascade for excited isotopes has no entries in RDM database
1063  // Decay mode is isomeric transition.
1064  G4ITDecayChannel* anITChannel = new G4ITDecayChannel
1065  (GetVerboseLevel(), (const G4Ions*) &theParentNucleus, 1);
1066  anITChannel->SetICM(applyICM);
1067  anITChannel->SetARM(applyARM);
1068  anITChannel->SetHLThreshold(halflifethreshold);
1069  theDecayTable->Insert(anITChannel);
1070  }
1071  if (!theDecayTable) {
1072 
1073  // There is no radioactive decay data for this nucleus. Return a null
1074  // decay table.
1075  G4cerr << "G4RadoactiveDecay::LoadDecayTable() : cannot find ion radioactive decay file "
1076  << G4endl;
1077  theDecayTable = 0;
1078  return theDecayTable;
1079  }
1080 
1081  if (theDecayTable && GetVerboseLevel() > 1) {
1082  theDecayTable->DumpInfo();
1083  }
1084 
1085 #ifdef G4MULTITHREADED
1086  (*master_dkmap)[key] = theDecayTable; // store in master library
1087 #endif
1088  return theDecayTable;
1089 }
1090 
1091 
1092 void G4RadioactiveDecay::AddUserDecayDataFile(G4int Z, G4int A,G4String filename)
1093 {
1094  if (Z < 1 || A < 2) G4cout << "Z and A not valid!" << G4endl;
1095 
1096  std::ifstream DecaySchemeFile(filename);
1097  if (DecaySchemeFile) {
1098  G4int ID_ion = A*1000 + Z;
1099  theUserRadioactiveDataFiles[ID_ion] = filename;
1100  theIsotopeTable->AddUserDecayDataFile(Z,A,filename);
1101  } else {
1102  G4cout << "The file " << filename << " does not exist!" << G4endl;
1103  }
1104 }
1105 
1106 
1107 void
1108 G4RadioactiveDecay::SetDecayRate(G4int theZ, G4int theA, G4double theE,
1109  G4int theG, std::vector<G4double> theRates,
1110  std::vector<G4double> theTaos)
1111 {
1112  //fill the decay rate vector
1113  theDecayRate.SetZ(theZ);
1114  theDecayRate.SetA(theA);
1115  theDecayRate.SetE(theE);
1116  theDecayRate.SetGeneration(theG);
1117  theDecayRate.SetDecayRateC(theRates);
1118  theDecayRate.SetTaos(theTaos);
1119 }
1120 
1121 
1122 void
1123 G4RadioactiveDecay::AddDecayRateTable(const G4ParticleDefinition& theParentNucleus)
1124 {
1125  // 1) To calculate all the coefficiecies required to derive the
1126  // radioactivities for all progeny of theParentNucleus
1127  //
1128  // 2) Add the coefficiencies to the decay rate table vector
1129  //
1130 
1131  //
1132  // Create and initialise variables used in the method.
1133  //
1134  theDecayRateVector.clear();
1135 
1136  G4int nGeneration = 0;
1137  std::vector<G4double> rates;
1138  std::vector<G4double> taos;
1139 
1140  // start rate is -1.
1141  // Eq.4.26 of the Technical Note
1142  rates.push_back(-1.);
1143  //
1144  //
1145  G4int A = ((const G4Ions*)(&theParentNucleus))->GetAtomicMass();
1146  G4int Z = ((const G4Ions*)(&theParentNucleus))->GetAtomicNumber();
1147  G4double E = ((const G4Ions*)(&theParentNucleus))->GetExcitationEnergy();
1148  G4double tao = theParentNucleus.GetPDGLifeTime();
1149  if (tao < 0.) tao = 1e-100;
1150  taos.push_back(tao);
1151  G4int nEntry = 0;
1152 
1153  //fill the decay rate with the intial isotope data
1154  SetDecayRate(Z,A,E,nGeneration,rates,taos);
1155 
1156  // store the decay rate in decay rate vector
1157  theDecayRateVector.push_back(theDecayRate);
1158  nEntry++;
1159 
1160  // now start treating the sencondary generations..
1161 
1162  G4bool stable = false;
1163  G4int i;
1164  G4int j;
1165  G4VDecayChannel* theChannel = 0;
1166  G4NuclearDecayChannel* theNuclearDecayChannel = 0;
1167  G4ITDecayChannel* theITChannel = 0;
1168  G4BetaMinusDecayChannel *theBetaMinusChannel = 0;
1169  G4BetaPlusDecayChannel *theBetaPlusChannel = 0;
1170 // G4AlphaDecayChannel *theAlphaChannel = 0;
1171  G4AlphaDecay* theAlphaChannel = 0;
1172  G4ProtonDecayChannel *theProtonChannel = 0;
1173  G4RadioactiveDecayMode theDecayMode;
1174  G4double theBR = 0.0;
1175  G4int AP = 0;
1176  G4int ZP = 0;
1177  G4int AD = 0;
1178  G4int ZD = 0;
1179  G4double EP = 0.;
1180  std::vector<G4double> TP;
1181  std::vector<G4double> RP;
1182  G4ParticleDefinition *theDaughterNucleus;
1183  G4double daughterExcitation;
1184  G4ParticleDefinition *aParentNucleus;
1185  G4IonTable* theIonTable;
1186  G4DecayTable *aTempDecayTable;
1187  G4double theRate;
1188  G4double TaoPlus;
1189  G4int nS = 0;
1190  G4int nT = nEntry;
1191  G4double brs[8];
1192  //
1193  theIonTable =
1194  (G4IonTable*)(G4ParticleTable::GetParticleTable()->GetIonTable());
1195 
1196  while (!stable) {
1197  nGeneration++;
1198  for (j = nS; j < nT; j++) {
1199  ZP = theDecayRateVector[j].GetZ();
1200  AP = theDecayRateVector[j].GetA();
1201  EP = theDecayRateVector[j].GetE();
1202  RP = theDecayRateVector[j].GetDecayRateC();
1203  TP = theDecayRateVector[j].GetTaos();
1204  if (GetVerboseLevel() > 0) {
1205  G4cout << "G4RadioactiveDecay::AddDecayRateTable : daughters of ("
1206  << ZP << ", " << AP << ", " << EP
1207  << ") are being calculated, generation = " << nGeneration
1208  << G4endl;
1209  }
1210 
1211  aParentNucleus = theIonTable->GetIon(ZP,AP,EP);
1212  aTempDecayTable = GetDecayTable(aParentNucleus);
1213 
1214  G4DecayTable* theDecayTable = new G4DecayTable();
1215  for (G4int k = 0; k < 8; k++) brs[k] = 0.0;
1216 
1217  // Go through the decay table and to combine the same decay channels
1218  for (i = 0; i < aTempDecayTable->entries(); i++) {
1219  theChannel = aTempDecayTable->GetDecayChannel(i);
1220  theNuclearDecayChannel = static_cast<G4NuclearDecayChannel*>(theChannel);
1221  theDecayMode = theNuclearDecayChannel->GetDecayMode();
1222  daughterExcitation = theNuclearDecayChannel->GetDaughterExcitation();
1223  theDaughterNucleus = theNuclearDecayChannel->GetDaughterNucleus() ;
1224  AD = ((const G4Ions*)(theDaughterNucleus))->GetAtomicMass();
1225  ZD = ((const G4Ions*)(theDaughterNucleus))->GetAtomicNumber();
1226  G4NuclearLevelManager* levelManager =
1227  G4NuclearLevelStore::GetInstance()->GetManager(ZD, AD);
1228  if (levelManager->NumberOfLevels() ) {
1229  const G4NuclearLevel* level =
1230  levelManager->NearestLevel (daughterExcitation);
1231 
1232  if (std::abs(daughterExcitation - level->Energy()) < levelTolerance) {
1233  // Level half-life is in ns and the threshold is set to 1 micros
1234  // by default, user can set it via the UI command
1235  if (level->HalfLife()*ns >= halflifethreshold){
1236  // save the metastable nucleus
1237  theDecayTable->Insert(theChannel);
1238  } else {
1239  brs[theDecayMode] += theChannel->GetBR();
1240  }
1241  } else {
1242  brs[theDecayMode] += theChannel->GetBR();
1243  }
1244  } else {
1245  brs[theDecayMode] += theChannel->GetBR();
1246  }
1247  }
1248  brs[2] = brs[2]+brs[3]+brs[4]+brs[5];
1249  brs[3] = brs[4] =brs[5] = 0.0;
1250  for (i= 0; i<8; i++){
1251  if (brs[i] > 0.) {
1252  switch ( i ) {
1253  case 0:
1254  // Decay mode is isomeric transition
1255  theITChannel = new G4ITDecayChannel(0,
1256  (const G4Ions*) aParentNucleus, brs[0]);
1257  theDecayTable->Insert(theITChannel);
1258  break;
1259 
1260  case 1:
1261  // Decay mode is beta-
1262  theBetaMinusChannel = new G4BetaMinusDecayChannel(0, aParentNucleus,
1263  brs[1], 0.*MeV, 0.*MeV, 1, false, 0);
1264  theDecayTable->Insert(theBetaMinusChannel);
1265  break;
1266 
1267  case 2:
1268  // Decay mode is beta+ + EC.
1269  theBetaPlusChannel = new G4BetaPlusDecayChannel(GetVerboseLevel(),
1270  aParentNucleus, brs[2], 0.*MeV, 0.*MeV, 1, false, 0);
1271  theDecayTable->Insert(theBetaPlusChannel);
1272  break;
1273 
1274  case 6:
1275  // Decay mode is alpha.
1276 /*
1277  theAlphaChannel = new G4AlphaDecayChannel(GetVerboseLevel(),
1278  aParentNucleus,
1279  brs[6], 0.*MeV, 0.*MeV);
1280 */
1281  theAlphaChannel = new G4AlphaDecay(aParentNucleus, brs[6], 0.*MeV,
1282  0.*MeV);
1283 
1284  theDecayTable->Insert(theAlphaChannel);
1285  break;
1286 
1287  case 7:
1288  // Decay mode is proton.
1289  theProtonChannel = new G4ProtonDecayChannel(GetVerboseLevel(),
1290  aParentNucleus,
1291  brs[7], 0.*MeV, 0.*MeV);
1292  theDecayTable->Insert(theProtonChannel);
1293  break;
1294 
1295  default:
1296  break;
1297  }
1298  }
1299  }
1300 
1301  // loop over all branches in theDecayTable
1302  //
1303  for (i = 0; i < theDecayTable->entries(); i++){
1304  theChannel = theDecayTable->GetDecayChannel(i);
1305  theNuclearDecayChannel = static_cast<G4NuclearDecayChannel*>(theChannel);
1306  theBR = theChannel->GetBR();
1307  theDaughterNucleus = theNuclearDecayChannel->GetDaughterNucleus();
1308  // First check if the decay of the original nucleus is an IT channel,
1309  // if true create a new groud-level nucleus
1310  if (theNuclearDecayChannel->GetDecayMode() == IT && nGeneration == 1) {
1311  A = ((const G4Ions*)(theDaughterNucleus))->GetAtomicMass();
1312  Z = ((const G4Ions*)(theDaughterNucleus))->GetAtomicNumber();
1313  theDaughterNucleus=theIonTable->GetIon(Z,A,0.);
1314  }
1315  if (IsApplicable(*theDaughterNucleus) && theBR &&
1316  aParentNucleus != theDaughterNucleus) {
1317  // need to make sure daughter has decay table
1318  aTempDecayTable = GetDecayTable(theDaughterNucleus);
1319 
1320  if (aTempDecayTable->entries() ) {
1321  A = ((const G4Ions*)(theDaughterNucleus))->GetAtomicMass();
1322  Z = ((const G4Ions*)(theDaughterNucleus))->GetAtomicNumber();
1323  E = ((const G4Ions*)(theDaughterNucleus))->GetExcitationEnergy();
1324 
1325  TaoPlus = theDaughterNucleus->GetPDGLifeTime();
1326  if (TaoPlus <= 0.) TaoPlus = 1e-100;
1327 
1328  // first set the taos, one simply need to add to the parent ones
1329  taos.clear();
1330  taos = TP;
1331  size_t k;
1332  //check that TaoPlus differs from other taos from at least 1.e5 relative difference
1333  //for (k = 0; k < TP.size(); k++){
1334  //if (std::abs((TaoPlus-TP[k])/TP[k])<1.e-5 ) TaoPlus=1.00001*TP[k];
1335  //}
1336  taos.push_back(TaoPlus);
1337  // now calculate the coefficiencies
1338  //
1339  // they are in two parts, first the less than n ones
1340  // Eq 4.24 of the TN
1341  rates.clear();
1342  long double ta1,ta2;
1343  ta2 = (long double)TaoPlus;
1344  for (k = 0; k < RP.size(); k++){
1345  ta1 = (long double)TP[k];
1346  if (ta1 == ta2) {
1347  theRate = 1.e100;
1348  } else {
1349  theRate = ta1/(ta1-ta2);
1350  }
1351  theRate = theRate * theBR * RP[k];
1352  rates.push_back(theRate);
1353  }
1354 
1355  // the sencond part: the n:n coefficiency
1356  // Eq 4.25 of the TN. Note Yn+1 is zero apart from Y1 which is -1
1357  // as treated at line 1013
1358  theRate = 0.;
1359  long double aRate, aRate1;
1360  aRate1 = 0.L;
1361  for (k = 0; k < RP.size(); k++){
1362  ta1 = (long double)TP[k];
1363  if (ta1 == ta2 ) {
1364  aRate = 1.e100;
1365  } else {
1366  aRate = ta2/(ta1-ta2);
1367  }
1368  aRate = aRate * (long double)(theBR * RP[k]);
1369  aRate1 += aRate;
1370  }
1371  theRate = -aRate1;
1372  rates.push_back(theRate);
1373  SetDecayRate (Z,A,E,nGeneration,rates,taos);
1374  theDecayRateVector.push_back(theDecayRate);
1375  nEntry++;
1376  }
1377  } // end of testing daughter nucleus
1378  } // end of i loop( the branches)
1379  // delete theDecayTable;
1380 
1381  } //end of for j loop
1382  nS = nT;
1383  nT = nEntry;
1384  if (nS == nT) stable = true;
1385  }
1386 
1387  // end of while loop
1388  // the calculation completed here
1389 
1390 
1391  // fill the first part of the decay rate table
1392  // which is the name of the original particle (isotope)
1393  theDecayRateTable.SetIonName(theParentNucleus.GetParticleName());
1394 
1395  // now fill the decay table with the newly completed decay rate vector
1396  theDecayRateTable.SetItsRates(theDecayRateVector);
1397 
1398  // finally add the decayratetable to the tablevector
1399  theDecayRateTableVector.push_back(theDecayRateTable);
1400 }
1401 
1403 // //
1404 // SetSourceTimeProfile //
1405 // read in the source time profile function (histogram) //
1406 // //
1408 
1409 void G4RadioactiveDecay::SetSourceTimeProfile(G4String filename)
1410 {
1411  std::ifstream infile ( filename, std::ios::in );
1412  if (!infile) {
1413  G4ExceptionDescription ed;
1414  ed << " Could not open file " << filename << G4endl;
1415  G4Exception("G4RadioactiveDecay::SetSourceTimeProfile()", "HAD_RDM_001",
1416  FatalException, ed);
1417  }
1418 
1419  G4double bin, flux;
1420  NSourceBin = -1;
1421  while (infile >> bin >> flux ) {
1422  NSourceBin++;
1423  if (NSourceBin > 99) {
1424  G4Exception("G4RadioactiveDecay::SetSourceTimeProfile()", "HAD_RDM_002",
1425  FatalException, "Input source time file too big (>100 rows)");
1426 
1427  } else {
1428  SBin[NSourceBin] = bin * s;
1429  SProfile[NSourceBin] = flux;
1430  }
1431  }
1432  SetAnalogueMonteCarlo(0);
1433  infile.close();
1434 
1435 #ifdef G4VERBOSE
1436  if (GetVerboseLevel()>1)
1437  {G4cout <<" Source Timeprofile Nbin = " << NSourceBin <<G4endl;}
1438 #endif
1439 }
1440 
1442 // //
1443 // SetDecayBiasProfile //
1444 // read in the decay bias scheme function (histogram) //
1445 // //
1447 
1448 void G4RadioactiveDecay::SetDecayBias(G4String filename)
1449 {
1450 
1451  std::ifstream infile(filename, std::ios::in);
1452  if (!infile) G4Exception("G4RadioactiveDecay::SetDecayBias()", "HAD_RDM_003",
1453  FatalException, "Unable to open bias data file" );
1454 
1455  G4double bin, flux;
1456  G4int dWindows = 0;
1457  G4int i ;
1458 
1459  theRadioactivityTables.clear();
1460 
1461  NDecayBin = -1;
1462  while (infile >> bin >> flux ) {
1463  NDecayBin++;
1464  if (NDecayBin > 99) {
1465  G4Exception("G4RadioactiveDecay::SetDecayBias()", "HAD_RDM_004",
1466  FatalException, "Input bias file too big (>100 rows)" );
1467  } else {
1468  DBin[NDecayBin] = bin * s;
1469  DProfile[NDecayBin] = flux;
1470  if (flux > 0.) {
1471  decayWindows[NDecayBin] = dWindows;
1472  dWindows++;
1473  G4RadioactivityTable *rTable = new G4RadioactivityTable() ;
1474  theRadioactivityTables.push_back(rTable);
1475  }
1476  }
1477  }
1478  for ( i = 1; i<= NDecayBin; i++) DProfile[i] += DProfile[i-1];
1479  for ( i = 0; i<= NDecayBin; i++) DProfile[i] /= DProfile[NDecayBin];
1480  // converted to accumulated probabilities
1481 
1482  SetAnalogueMonteCarlo(0);
1483  infile.close();
1484 
1485 #ifdef G4VERBOSE
1486  if (GetVerboseLevel()>1)
1487  {G4cout <<" Decay Bias Profile Nbin = " << NDecayBin <<G4endl;}
1488 #endif
1489 }
1490 
1492 // //
1493 // DecayIt //
1494 // //
1496 
1497 G4VParticleChange*
1498 G4RadioactiveDecay::DecayIt(const G4Track& theTrack, const G4Step&)
1499 {
1500  // Initialize G4ParticleChange object, get particle details and decay table
1501 
1502  fParticleChangeForRadDecay.Initialize(theTrack);
1503  const G4DynamicParticle* theParticle = theTrack.GetDynamicParticle();
1504  const G4ParticleDefinition* theParticleDef = theParticle->GetDefinition();
1505 
1506  // First check whether RDM applies to the current logical volume
1507  if (!isAllVolumesMode) {
1508  if (!std::binary_search(ValidVolumes.begin(), ValidVolumes.end(),
1509  theTrack.GetVolume()->GetLogicalVolume()->GetName())) {
1510 #ifdef G4VERBOSE
1511  if (GetVerboseLevel()>0) {
1512  G4cout <<"G4RadioactiveDecay::DecayIt : "
1513  << theTrack.GetVolume()->GetLogicalVolume()->GetName()
1514  << " is not selected for the RDM"<< G4endl;
1515  G4cout << " There are " << ValidVolumes.size() << " volumes" << G4endl;
1516  G4cout << " The Valid volumes are " << G4endl;
1517  for (size_t i = 0; i< ValidVolumes.size(); i++)
1518  G4cout << ValidVolumes[i] << G4endl;
1519  }
1520 #endif
1521  fParticleChangeForRadDecay.SetNumberOfSecondaries(0);
1522 
1523  // Kill the parent particle.
1524  fParticleChangeForRadDecay.ProposeTrackStatus(fStopAndKill) ;
1525  fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0);
1526  ClearNumberOfInteractionLengthLeft();
1527  return &fParticleChangeForRadDecay;
1528  }
1529  }
1530 
1531  // Now check if particle is valid for RDM
1532  if (!(IsApplicable(*theParticleDef) ) ) {
1533  // Particle is not an ion or is outside the nucleuslimits for decay
1534 
1535 #ifdef G4VERBOSE
1536  if (GetVerboseLevel()>0) {
1537  G4cerr << "G4RadioactiveDecay::DecayIt : "
1538  << theParticleDef->GetParticleName()
1539  << " is not a valid nucleus for the RDM"<< G4endl;
1540  }
1541 #endif
1542  fParticleChangeForRadDecay.SetNumberOfSecondaries(0);
1543 
1544  // Kill the parent particle
1545  fParticleChangeForRadDecay.ProposeTrackStatus(fStopAndKill) ;
1546  fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0);
1547  ClearNumberOfInteractionLengthLeft();
1548  return &fParticleChangeForRadDecay;
1549  }
1550 
1551  G4DecayTable* theDecayTable = GetDecayTable(theParticleDef);
1552 
1553  if (theDecayTable == 0 || theDecayTable->entries() == 0) {
1554  // No data in the decay table. Set particle change parameters
1555  // to indicate this.
1556 #ifdef G4VERBOSE
1557  if (GetVerboseLevel() > 0) {
1558  G4cerr <<"G4RadioactiveDecay::DecayIt : decay table not defined for ";
1559  G4cerr <<theParticleDef->GetParticleName() <<G4endl;
1560  }
1561 #endif
1562  fParticleChangeForRadDecay.SetNumberOfSecondaries(0);
1563 
1564  // Kill the parent particle.
1565  fParticleChangeForRadDecay.ProposeTrackStatus(fStopAndKill) ;
1566  fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0);
1567  ClearNumberOfInteractionLengthLeft();
1568  return &fParticleChangeForRadDecay;
1569 
1570  } else {
1571  // Data found. Try to decay nucleus
1572  G4double energyDeposit = 0.0;
1573  G4double finalGlobalTime = theTrack.GetGlobalTime();
1574  G4double finalLocalTime = theTrack.GetLocalTime();
1575  G4int index;
1576  G4ThreeVector currentPosition;
1577  currentPosition = theTrack.GetPosition();
1578 
1579  // Check whether use Analogue or VR implementation
1580  if (AnalogueMC) {
1581 #ifdef G4VERBOSE
1582  if (GetVerboseLevel() > 0)
1583  G4cout <<"DecayIt: Analogue MC version " << G4endl;
1584 #endif
1585 
1586  G4DecayProducts* products = DoDecay(*theParticleDef);
1587 
1588  // Check if the product is the same as input and kill the track if
1589  // necessary to prevent infinite loop (11/05/10, F.Lei)
1590  if ( products->entries() == 1) {
1591  fParticleChangeForRadDecay.SetNumberOfSecondaries(0);
1592  fParticleChangeForRadDecay.ProposeTrackStatus(fStopAndKill);
1593  fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0);
1594  ClearNumberOfInteractionLengthLeft();
1595  return &fParticleChangeForRadDecay;
1596  }
1597 
1598  // Get parent particle information and boost the decay products to the
1599  // laboratory frame based on this information.
1600 
1601  //The Parent Energy used for the boost should be the total energy of
1602  // the nucleus of the parent ion without the energy of the shell electrons
1603  // (correction for bug 1359 by L. Desorgher)
1604  G4double ParentEnergy = theParticle->GetKineticEnergy()
1605  + theParticle->GetParticleDefinition()->GetPDGMass();
1606  G4ThreeVector ParentDirection(theParticle->GetMomentumDirection());
1607 
1608  if (theTrack.GetTrackStatus() == fStopButAlive) {
1609  //this condition seems to be always True, further investigation is needed (L.Desorgher)
1610 
1611  // The particle is decayed at rest.
1612  // since the time is still for rest particle in G4 we need to add the
1613  // additional time lapsed between the particle come to rest and the
1614  // actual decay. This time is simply sampled with the mean-life of
1615  // the particle. But we need to protect the case PDGTime < 0.
1616  // (F.Lei 11/05/10)
1617  G4double temptime = -std::log( G4UniformRand())
1618  *theParticleDef->GetPDGLifeTime();
1619  if (temptime < 0.) temptime = 0.;
1620  finalGlobalTime += temptime;
1621  finalLocalTime += temptime;
1622  energyDeposit += theParticle->GetKineticEnergy();
1623  }
1624  products->Boost(ParentEnergy, ParentDirection);
1625 
1626  // Add products in theParticleChangeForRadDecay.
1627  G4int numberOfSecondaries = products->entries();
1628  fParticleChangeForRadDecay.SetNumberOfSecondaries(numberOfSecondaries);
1629 #ifdef G4VERBOSE
1630  if (GetVerboseLevel()>1) {
1631  G4cout <<"G4RadioactiveDecay::DecayIt : Decay vertex :";
1632  G4cout <<" Time: " <<finalGlobalTime/ns <<"[ns]";
1633  G4cout <<" X:" <<(theTrack.GetPosition()).x() /cm <<"[cm]";
1634  G4cout <<" Y:" <<(theTrack.GetPosition()).y() /cm <<"[cm]";
1635  G4cout <<" Z:" <<(theTrack.GetPosition()).z() /cm <<"[cm]";
1636  G4cout << G4endl;
1637  G4cout <<"G4Decay::DecayIt : decay products in Lab. Frame" <<G4endl;
1638  products->DumpInfo();
1639  products->IsChecked();
1640  }
1641 #endif
1642  for (index=0; index < numberOfSecondaries; index++) {
1643  G4Track* secondary = new G4Track(products->PopProducts(),
1644  finalGlobalTime, currentPosition);
1645  secondary->SetGoodForTrackingFlag();
1646  secondary->SetTouchableHandle(theTrack.GetTouchableHandle());
1647  fParticleChangeForRadDecay.AddSecondary(secondary);
1648  }
1649  delete products;
1650  // end of analogue MC algorithm
1651 
1652  } else {
1653  // Variance Reduction Method
1654 #ifdef G4VERBOSE
1655  if (GetVerboseLevel()>0)
1656  G4cout << "DecayIt: Variance Reduction version " << G4endl;
1657 #endif
1658  if (!IsRateTableReady(*theParticleDef)) {
1659  // if the decayrates are not ready, calculate them and
1660  // add to the rate table vector
1661  AddDecayRateTable(*theParticleDef);
1662  }
1663  //retrieve the rates
1664  GetDecayRateTable(*theParticleDef);
1665 
1666  // declare some of the variables required in the implementation
1667  G4ParticleDefinition* parentNucleus;
1668  G4IonTable* theIonTable;
1669  G4int PZ;
1670  G4int PA;
1671  G4double PE;
1672  G4String keyName;
1673  std::vector<G4double> PT;
1674  std::vector<G4double> PR;
1675  G4double taotime;
1676  long double decayRate;
1677 
1678  size_t i;
1679  size_t j;
1680  G4int numberOfSecondaries;
1681  G4int totalNumberOfSecondaries = 0;
1682  G4double currentTime = 0.;
1683  G4int ndecaych;
1684  G4DynamicParticle* asecondaryparticle;
1685  std::vector<G4DynamicParticle*> secondaryparticles;
1686  std::vector<G4double> pw;
1687  std::vector<G4double> ptime;
1688  pw.clear();
1689  ptime.clear();
1690 
1691  //now apply the nucleus splitting
1692  for (G4int n = 0; n < NSplit; n++) {
1693  // Get the decay time following the decay probability function
1694  // suppllied by user
1695  G4double theDecayTime = GetDecayTime();
1696  G4int nbin = GetDecayTimeBin(theDecayTime);
1697 
1698  // calculate the first part of the weight function
1699  G4double weight1 = 1.;
1700  if (nbin == 1) {
1701  weight1 = 1./DProfile[nbin-1]
1702  *(DBin[nbin]-DBin[nbin-1])/NSplit;
1703  } else if (nbin > 1) {
1704  weight1 = 1./(DProfile[nbin]-DProfile[nbin-2])
1705  *(DBin[nbin]-DBin[nbin-1])/NSplit;
1706  }
1707 
1708  // it should be calculated in seconds
1709  weight1 /= s ;
1710 
1711  // loop over all the possible secondaries of the nucleus
1712  // the first one is itself.
1713  for (i = 0; i < theDecayRateVector.size(); i++) {
1714  PZ = theDecayRateVector[i].GetZ();
1715  PA = theDecayRateVector[i].GetA();
1716  PE = theDecayRateVector[i].GetE();
1717  PT = theDecayRateVector[i].GetTaos();
1718  PR = theDecayRateVector[i].GetDecayRateC();
1719 
1720  // Calculate the decay rate of the isotope
1721  // decayRate is the radioactivity of isotope (PZ,PA,PE) at the
1722  // time 'theDecayTime'
1723  // it will be used to calculate the statistical weight of the
1724  // decay products of this isotope
1725 
1726  // G4cout <<"PA= "<< PA << " PZ= " << PZ << " PE= "<< PE <<G4endl;
1727  decayRate = 0.L;
1728  for (j = 0; j < PT.size(); j++) {
1729  taotime = GetTaoTime(theDecayTime,PT[j]);
1730  decayRate -= PR[j] * (long double)taotime;
1731  // Eq.4.23 of of the TN
1732  // note the negative here is required as the rate in the
1733  // equation is defined to be negative,
1734  // i.e. decay away, but we need positive value here.
1735 
1736  // G4cout << j << "\t"<< PT[j]/s <<"\t"<<PR[j]<< "\t"
1737  // << decayRate << G4endl;
1738  }
1739 
1740  // add the isotope to the radioactivity tables
1741  // G4cout <<theDecayTime/s <<"\t"<<nbin<<G4endl;
1742  // G4cout << theTrack.GetWeight() <<"\t"<<weight1<<"\t"<<decayRate<< G4endl;
1743  theRadioactivityTables[decayWindows[nbin-1]]->AddIsotope(PZ,PA,PE,weight1*decayRate,theTrack.GetWeight());
1744 
1745  // Now calculate the statistical weight
1746  // One needs to fold the source bias function with the decaytime
1747  // also need to include the track weight! (F.Lei, 28/10/10)
1748  G4double weight = weight1*decayRate*theTrack.GetWeight();
1749 
1750  // decay the isotope
1751  theIonTable = (G4IonTable *)(G4ParticleTable::GetParticleTable()->GetIonTable());
1752  parentNucleus = theIonTable->GetIon(PZ,PA,PE);
1753 
1754  // Create a temprary products buffer.
1755  // Its contents to be transfered to the products at the end of the loop
1756  G4DecayProducts* tempprods = 0;
1757 
1758  // Decide whether to apply branching ratio bias or not
1759  if (BRBias) {
1760  G4DecayTable* decayTable = GetDecayTable(parentNucleus);
1761 
1762  ndecaych = G4int(decayTable->entries()*G4UniformRand());
1763  G4VDecayChannel* theDecayChannel = decayTable->GetDecayChannel(ndecaych);
1764  if (theDecayChannel == 0) {
1765  // Decay channel not found.
1766 #ifdef G4VERBOSE
1767  if (GetVerboseLevel()>0) {
1768  G4cerr << " G4RadioactiveDecay::DoIt : cannot determine decay channel ";
1769  G4cerr << " for this nucleus; decay as if no biasing active ";
1770  G4cerr << G4endl;
1771  decayTable ->DumpInfo();
1772  }
1773 #endif
1774  tempprods = DoDecay(*parentNucleus); // DHW 6 Dec 2010 - do decay as if no biasing
1775  // to avoid deref of temppprods = 0
1776  } else {
1777  // A decay channel has been identified, so execute the DecayIt.
1778  G4double tempmass = parentNucleus->GetPDGMass();
1779  tempprods = theDecayChannel->DecayIt(tempmass);
1780  weight *= (theDecayChannel->GetBR())*(decayTable->entries());
1781  }
1782  } else {
1783  tempprods = DoDecay(*parentNucleus);
1784  }
1785 
1786  // save the secondaries for buffers
1787  numberOfSecondaries = tempprods->entries();
1788  currentTime = finalGlobalTime + theDecayTime;
1789  for (index = 0; index < numberOfSecondaries; index++) {
1790  asecondaryparticle = tempprods->PopProducts();
1791  if (asecondaryparticle->GetDefinition()->GetBaryonNumber() < 5) {
1792  pw.push_back(weight);
1793  ptime.push_back(currentTime);
1794  secondaryparticles.push_back(asecondaryparticle);
1795  }
1796  }
1797  delete tempprods;
1798 
1799  } // end of i loop
1800  } // end of n loop
1801 
1802  // now deal with the secondaries in the two stl containers
1803  // and submmit them back to the tracking manager
1804  totalNumberOfSecondaries = pw.size();
1805  fParticleChangeForRadDecay.SetNumberOfSecondaries(totalNumberOfSecondaries);
1806  for (index=0; index < totalNumberOfSecondaries; index++) {
1807  G4Track* secondary = new G4Track(secondaryparticles[index],
1808  ptime[index], currentPosition);
1809  secondary->SetGoodForTrackingFlag();
1810  secondary->SetTouchableHandle(theTrack.GetTouchableHandle());
1811  secondary->SetWeight(pw[index]);
1812  fParticleChangeForRadDecay.AddSecondary(secondary);
1813  }
1814  // make sure the original track is set to stop and its kinematic energy collected
1815  //
1816  //theTrack.SetTrackStatus(fStopButAlive);
1817  //energyDeposit += theParticle->GetKineticEnergy();
1818 
1819  } // End of Variance Reduction
1820 
1821  // Kill the parent particle
1822  fParticleChangeForRadDecay.ProposeTrackStatus(fStopAndKill) ;
1823  fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(energyDeposit);
1824  fParticleChangeForRadDecay.ProposeLocalTime(finalLocalTime);
1825  // Reset NumberOfInteractionLengthLeft.
1826  ClearNumberOfInteractionLengthLeft();
1827 
1828  return &fParticleChangeForRadDecay ;
1829  }
1830 }
1831 
1832 
1833 G4DecayProducts*
1834 G4RadioactiveDecay::DoDecay(const G4ParticleDefinition& theParticleDef)
1835 {
1836  G4DecayProducts* products = 0;
1837  G4DecayTable* theDecayTable = GetDecayTable(&theParticleDef);
1838 
1839  // Choose a decay channel.
1840 #ifdef G4VERBOSE
1841  if (GetVerboseLevel() > 0) G4cout << "Select a channel..." << G4endl;
1842 #endif
1843 
1844  G4VDecayChannel* theDecayChannel = theDecayTable->SelectADecayChannel();
1845 
1846  if (theDecayChannel == 0) {
1847  // Decay channel not found.
1848  G4cerr << "G4RadioactiveDecay::DoIt : can not determine decay channel";
1849  G4cerr << G4endl;
1850  } else {
1851  // A decay channel has been identified, so execute the DecayIt.
1852 #ifdef G4VERBOSE
1853  if (GetVerboseLevel() > 1) {
1854  G4cerr << "G4RadioactiveDecay::DoIt : selected decay channel addr:";
1855  G4cerr << theDecayChannel << G4endl;
1856  }
1857 #endif
1858  products = theDecayChannel->DecayIt(theParticleDef.GetPDGMass() );
1859 
1860  // Apply directional bias if requested by user
1861  CollimateDecay(products);
1862  }
1863 
1864  return products;
1865 }
1866 
1867 
1868 // Apply directional bias for "visible" daughters (e+-, gamma, n, p, alpha)
1869 
1870 void G4RadioactiveDecay::CollimateDecay(G4DecayProducts* products) {
1871  if (origin == forceDecayDirection) return; // No collimation requested
1872  if (180.*deg == forceDecayHalfAngle) return;
1873  if (0 == products || 0 == products->entries()) return;
1874 
1875 #ifdef G4VERBOSE
1876  if (GetVerboseLevel() > 0) G4cout << "Begin of CollimateDecay..." << G4endl;
1877 #endif
1878 
1879  // Particles suitable for directional biasing (for if-blocks below)
1880  static const G4ParticleDefinition* electron = G4Electron::Definition();
1881  static const G4ParticleDefinition* positron = G4Positron::Definition();
1882  static const G4ParticleDefinition* neutron = G4Neutron::Definition();
1883  static const G4ParticleDefinition* gamma = G4Gamma::Definition();
1884  static const G4ParticleDefinition* alpha = G4Alpha::Definition();
1885  static const G4ParticleDefinition* proton = G4Proton::Definition();
1886 
1887  G4ThreeVector newDirection; // Re-use to avoid memory churn
1888  for (G4int i=0; i<products->entries(); i++) {
1889  G4DynamicParticle* daughter = (*products)[i];
1890  const G4ParticleDefinition* daughterType =
1891  daughter->GetParticleDefinition();
1892  if (daughterType == electron || daughterType == positron ||
1893  daughterType == neutron || daughterType == gamma ||
1894  daughterType == alpha || daughterType == proton) CollimateDecayProduct(daughter);
1895  }
1896 }
1897 
1898 void G4RadioactiveDecay::CollimateDecayProduct(G4DynamicParticle* daughter) {
1899 #ifdef G4VERBOSE
1900  if (GetVerboseLevel() > 1) {
1901  G4cout << "CollimateDecayProduct for daughter "
1902  << daughter->GetParticleDefinition()->GetParticleName() << G4endl;
1903  }
1904 #endif
1905 
1906  G4ThreeVector collimate = ChooseCollimationDirection();
1907  if (origin != collimate) daughter->SetMomentumDirection(collimate);
1908 }
1909 
1910 
1911 // Choose random direction within collimation cone
1912 
1913 G4ThreeVector G4RadioactiveDecay::ChooseCollimationDirection() const {
1914  if (origin == forceDecayDirection) return origin; // Don't do collimation
1915  if (forceDecayHalfAngle == 180.*deg) return origin;
1916 
1917  G4ThreeVector dir = forceDecayDirection;
1918 
1919  // Return direction offset by random throw
1920  if (forceDecayHalfAngle > 0.) {
1921  // Generate uniform direction around central axis
1922  G4double phi = 2.*pi*G4UniformRand();
1923  G4double cosMin = std::cos(forceDecayHalfAngle);
1924  G4double cosTheta = (1.-cosMin)*G4UniformRand() + cosMin; // [cosMin,1.)
1925 
1926  dir.setPhi(dir.phi()+phi);
1927  dir.setTheta(dir.theta()+std::acos(cosTheta));
1928  }
1929 
1930 #ifdef G4VERBOSE
1931  if (GetVerboseLevel()>1)
1932  G4cout << " ChooseCollimationDirection returns " << dir << G4endl;
1933 #endif
1934 
1935  return dir;
1936 }
G4ParticleChangeForRadDecay.hh
cm
static const double cm
Definition: G4SIunits.hh:106
allowed
Definition: G4BetaDecayType.hh:41
G4BetaMinusDecayChannel.hh
G4IonTable::RegisterIsotopeTable
void RegisterIsotopeTable(G4VIsotopeTable *table)
Definition: G4IonTable.cc:1457
G4VDecayChannel::GetBR
G4double GetBR() const
Definition: G4VDecayChannel.hh:345
G4RadioactiveDecay::theDecayRateTable
G4RadioactiveDecayRateVector theDecayRateTable
Definition: G4RadioactiveDecay.hh:299
Proton2
Definition: G4RadioactiveDecayMode.hh:66
MeV
static const double MeV
Definition: G4SIunits.hh:193
G4RadioactiveDecay::DoDecay
G4DecayProducts * DoDecay(const G4ParticleDefinition &theParticleDef)
Definition: G4RadioactiveDecay.cc:1834
G4VDecayChannel::SetBR
void SetBR(G4double value)
Definition: G4VDecayChannel.hh:342
G4LossTableManager::Instance
static G4LossTableManager * Instance()
Definition: G4LossTableManager.cc:114
G4Track::GetLocalTime
G4double GetLocalTime() const
G4RadioactiveDecay::SelectAVolume
void SelectAVolume(const G4String aVolume)
Definition: G4RadioactiveDecay.cc:265
DecayTableMap
std::map< G4String, G4DecayTable * > DecayTableMap
Definition: G4RadioactiveDecay.hh:79
G4ExceptionDescription
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
G4LogicalVolume::GetName
G4String GetName() const
G4DynamicParticle::GetKineticEnergy
G4double GetKineticEnergy() const
G4ThreeVector
CLHEP::Hep3Vector G4ThreeVector
Definition: G4ThreeVector.hh:42
G4RadioactiveDecayMode.hh
G4RadioactiveDecay::IsRateTableReady
G4bool IsRateTableReady(const G4ParticleDefinition &)
Definition: G4RadioactiveDecay.cc:354
G4NuclearLevelManager.hh
fStopButAlive
Definition: G4TrackStatus.hh:54
G4Track::GetDynamicParticle
const G4DynamicParticle * GetDynamicParticle() const
Alpha
Definition: G4RadioactiveDecayMode.hh:65
G4RadioactiveDecay::ValidVolumes
std::vector< G4String > ValidVolumes
Definition: G4RadioactiveDecay.hh:294
LshellEC
Definition: G4RadioactiveDecayMode.hh:65
G4NuclearLevel::HalfLife
G4double HalfLife() const
Definition: G4NuclearLevel.cc:222
G4RadioactiveDecayRate::SetTaos
void SetTaos(std::vector< G4double > value)
Definition: G4RadioactiveDecayRate.hh:106
z
G4double z
Definition: TRTMaterials.hh:39
G4String::strip
G4String strip(G4int strip_Type=trailing, char c=' ')
G4RadioactiveDecayRate::SetDecayRateC
void SetDecayRateC(std::vector< G4double > value)
Definition: G4RadioactiveDecayRate.hh:104
G4RIsotopeTable::AddUserDecayDataFile
void AddUserDecayDataFile(G4int Z, G4int A, G4String filename)
Definition: G4RIsotopeTable.cc:272
G4ITDecayChannel.hh
G4ParticleDefinition::GetPDGStable
G4bool GetPDGStable() const
G4IonTable::GetIon
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
Definition: G4IonTable.cc:463
G4NuclearLevel::Energy
G4double Energy() const
Definition: G4NuclearLevel.cc:212
G4Alpha::Definition
static G4Alpha * Definition()
Definition: G4Alpha.cc:49
G4INCL::Math::pi
const G4double pi
Definition: G4INCLGlobals.hh:68
G4Track::GetPosition
const G4ThreeVector & GetPosition() const
G4NucleusLimits::GetZMax
G4int GetZMax() const
Definition: G4NucleusLimits.hh:100
G4ParticleChangeForRadDecay::AddSecondary
void AddSecondary(G4Track *aSecondary)
Definition: G4ParticleChangeForRadDecay.hh:71
G4Track::GetTrackStatus
G4TrackStatus GetTrackStatus() const
G4VDecayChannel::GetKinematicsName
const G4String & GetKinematicsName() const
Definition: G4VDecayChannel.hh:339
G4RadioactiveDecay::theDecayRateVector
G4RadioactiveDecayRates theDecayRateVector
Definition: G4RadioactiveDecay.hh:298
G4DynamicParticle.hh
G4Electron::Definition
static G4Electron * Definition()
Definition: G4Electron.cc:49
G4Track::SetTouchableHandle
void SetTouchableHandle(const G4TouchableHandle &apValue)
G4DecayTable::GetDecayChannel
G4VDecayChannel * GetDecayChannel(G4int index) const
Definition: G4DecayTable.hh:113
G4RadioactiveDecay::GetMeanLifeTime
G4double GetMeanLifeTime(const G4Track &theTrack, G4ForceCondition *condition)
Definition: G4RadioactiveDecay.cc:568
G4RadioactiveDecay::GetVerboseLevel
G4int GetVerboseLevel() const
Definition: G4RadioactiveDecay.hh:171
G4LossTableManager.hh
G4DynamicParticle::SetMomentumDirection
void SetMomentumDirection(const G4ThreeVector &aDirection)
G4NuclearLevelStore::GetInstance
static G4NuclearLevelStore * GetInstance()
Definition: G4NuclearLevelStore.cc:47
neutron
Definition: G4DataQuestionaire.hh:35
G4RadioactiveDecaymessenger.hh
G4DynamicParticle::GetDefinition
G4ParticleDefinition * GetDefinition() const
G4Ions.hh
a
G4double a
Definition: TRTMaterials.hh:39
G4RadioactiveDecay::SetSourceTimeProfile
void SetSourceTimeProfile(G4String filename)
Definition: G4RadioactiveDecay.cc:1409
G4RadioactiveDecay::IsApplicable
G4bool IsApplicable(const G4ParticleDefinition &)
Definition: G4RadioactiveDecay.cc:227
G4DecayProducts::Boost
void Boost(G4double totalEnergy, const G4ThreeVector &momentumDirection)
Definition: G4DecayProducts.cc:176
G4Positron::Definition
static G4Positron * Definition()
Definition: G4Positron.cc:49
G4VProcess::ClearNumberOfInteractionLengthLeft
void ClearNumberOfInteractionLengthLeft()
Definition: G4VProcess.hh:447
BetaPlus
Definition: G4RadioactiveDecayMode.hh:65
JustWarning
Definition: G4ExceptionSeverity.hh:64
G4int
int G4int
Definition: G4Types.hh:78
G4NucleusLimits::GetAMin
G4int GetAMin() const
Definition: G4NucleusLimits.hh:91
G4MUTEX_INITIALIZER
#define G4MUTEX_INITIALIZER
Definition: G4Threading.hh:175
G4ParticleDefinition::GetParticleName
const G4String & GetParticleName() const
Definition: G4ParticleDefinition.hh:159
G4NuclearLevelStore.hh
G4IonTable.hh
G4Track::SetWeight
void SetWeight(G4double aValue)
G4VParticleChange::ProposeLocalEnergyDeposit
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
BetaMinus
Definition: G4RadioactiveDecayMode.hh:65
G4RadioactiveDecayRateVector::SetItsRates
void SetItsRates(G4RadioactiveDecayRates arate)
Definition: G4RadioactiveDecayRateVector.hh:102
G4RadioactiveDecay::fParticleChangeForRadDecay
G4ParticleChangeForRadDecay fParticleChangeForRadDecay
Definition: G4RadioactiveDecay.hh:322
G4DynamicParticle::GetTotalMomentum
G4double GetTotalMomentum() const
G4Proton::Definition
static G4Proton * Definition()
Definition: G4Proton.cc:49
L
static const G4int L[nN]
Definition: G4ElectroNuclearCrossSection.cc:113
s
static const double s
Definition: G4SIunits.hh:150
G4DecayTable::entries
G4int entries() const
Definition: G4DecayTable.hh:100
G4BetaDecayCorrections.hh
G4ParticleChangeForDecay::Initialize
virtual void Initialize(const G4Track &)
Definition: G4ParticleChangeForDecay.cc:124
G4ParticleTable::GetIonTable
G4IonTable * GetIonTable() const
Definition: G4ParticleTable.cc:636
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:93
G4cout
G4GLOB_DLL std::ostream G4cout
G4Ions
Definition: G4Ions.hh:51
G4RadioactiveDecay::origin
static const G4ThreeVector origin
Definition: G4RadioactiveDecay.hh:285
G4RadioactiveDecay::SetDecayRate
void SetDecayRate(G4int, G4int, G4double, G4int, std::vector< G4double >, std::vector< G4double >)
Definition: G4RadioactiveDecay.cc:1108
deg
static const double deg
Definition: G4SIunits.hh:133
G4DynamicParticle::GetMass
G4double GetMass() const
G4Neutron.hh
G4NucleusLimits::GetAMax
G4int GetAMax() const
Definition: G4NucleusLimits.hh:94
G4VDecayChannel.hh
G4bool
bool G4bool
Definition: G4Types.hh:79
KshellEC
Definition: G4RadioactiveDecayMode.hh:65
G4DecayProducts::DumpInfo
void DumpInfo() const
Definition: G4DecayProducts.cc:311
G4RandGeneral
#define G4RandGeneral
Definition: Randomize.hh:90
G4DynamicParticle::GetMomentumDirection
const G4ThreeVector & GetMomentumDirection() const
G4DecayTable::DumpInfo
void DumpInfo() const
Definition: G4DecayTable.cc:101
G4RadioactiveDecay::CollimateDecayProduct
void CollimateDecayProduct(G4DynamicParticle *product)
Definition: G4RadioactiveDecay.cc:1898
G4Proton.hh
G4RadioactiveDecayMode
G4RadioactiveDecayMode
Definition: G4RadioactiveDecayMode.hh:64
G4Alpha.hh
G4LogicalVolumeStore::GetInstance
static G4LogicalVolumeStore * GetInstance()
Definition: G4LogicalVolumeStore.cc:181
G4Positron.hh
G4RadioactiveDecay::ChooseCollimationDirection
G4ThreeVector ChooseCollimationDirection() const
Definition: G4RadioactiveDecay.cc:1913
G4RadioactiveDecay::GetDecayTimeBin
G4int GetDecayTimeBin(const G4double aDecayTime)
Definition: G4RadioactiveDecay.cc:555
G4BetaDecayType.hh
G4VProcess::SetProcessSubType
void SetProcessSubType(G4int)
Definition: G4VProcess.hh:432
GeV
static const double GeV
Definition: G4SIunits.hh:196
G4ParticleDefinition::GetParticleType
const G4String & GetParticleType() const
Definition: G4ParticleDefinition.hh:180
G4Step
Definition: G4Step.hh:76
Randomize.hh
G4VAtomDeexcitation.hh
Beta2Plus
Definition: G4RadioactiveDecayMode.hh:66
G4DynamicParticle::GetParticleDefinition
const G4ParticleDefinition * GetParticleDefinition() const
n
const G4int n
Definition: G4UrQMD1_3Interface.hh:144
G4Track::GetGlobalTime
G4double GetGlobalTime() const
A
static const G4double A[nN]
Definition: G4ElectroNuclearCrossSection.cc:115
G4RadioactiveDecay::GetTaoTime
G4double GetTaoTime(const G4double, const G4double)
Definition: G4RadioactiveDecay.cc:388
G4RadioactiveDecay::theNucleusLimits
G4NucleusLimits theNucleusLimits
Definition: G4RadioactiveDecay.hh:270
G4Track::GetTouchableHandle
const G4TouchableHandle & GetTouchableHandle() const
G4BetaDecayType
G4BetaDecayType
Definition: G4BetaDecayType.hh:39
G4DecayTable::SelectADecayChannel
G4VDecayChannel * SelectADecayChannel()
Definition: G4DecayTable.cc:81
G4RadioactiveDecay::GetDecayTable
G4DecayTable * GetDecayTable(const G4ParticleDefinition *)
Definition: G4RadioactiveDecay.cc:248
G4NuclearLevelStore::GetManager
G4NuclearLevelManager * GetManager(G4int Z, G4int A)
Definition: G4NuclearLevelStore.cc:125
G4ProtonDecayChannel.hh
G4Exception
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4NuclearDecayChannel::GetDaughterNucleus
G4ParticleDefinition * GetDaughterNucleus()
Definition: G4NuclearDecayChannel.hh:106
G4DecayTable::Insert
void Insert(G4VDecayChannel *aChannel)
Definition: G4DecayTable.cc:60
G4DecayProducts.hh
G4RadioactiveDecay::GetMeanFreePath
G4double GetMeanFreePath(const G4Track &theTrack, G4double previousStepSize, G4ForceCondition *condition)
Definition: G4RadioactiveDecay.cc:611
G4RadioactiveDecay::theDecayRate
G4RadioactiveDecayRate theDecayRate
Definition: G4RadioactiveDecay.hh:297
G4RadioactiveDecay::levelTolerance
static const G4double levelTolerance
Definition: G4RadioactiveDecay.hh:305
G4Mutex
G4int G4Mutex
Definition: G4Threading.hh:173
G4Electron.hh
G4RadioactiveDecay::theDecayRateTableVector
G4RadioactiveDecayRateTable theDecayRateTableVector
Definition: G4RadioactiveDecay.hh:300
G4RadioactiveDecay::theUserRadioactiveDataFiles
std::map< G4int, G4String > theUserRadioactiveDataFiles
Definition: G4RadioactiveDecay.hh:308
nanosecond
static const double nanosecond
Definition: G4SIunits.hh:139
Proton
Definition: G4RadioactiveDecayMode.hh:66
G4VPhysicalVolume::GetLogicalVolume
G4LogicalVolume * GetLogicalVolume() const
G4ParticleDefinition::GetPDGMass
G4double GetPDGMass() const
Definition: G4ParticleDefinition.hh:161
fRadioactiveDecay
Definition: G4HadronicProcessType.hh:51
G4ParticleTable::GetParticleTable
static G4ParticleTable * GetParticleTable()
Definition: G4ParticleTable.cc:96
G4RadioactiveDecay::GetDecayRateTable
void GetDecayRateTable(const G4ParticleDefinition &)
Definition: G4RadioactiveDecay.cc:369
G4RadioactiveDecay::theRadioactiveDecaymessenger
G4RadioactiveDecaymessenger * theRadioactiveDecaymessenger
Definition: G4RadioactiveDecay.hh:267
G4VParticleChange::SetNumberOfSecondaries
void SetNumberOfSecondaries(G4int totSecondaries)
G4DecayProducts::PopProducts
G4DynamicParticle * PopProducts()
Definition: G4DecayProducts.cc:141
G4VProcess::pParticleChange
G4VParticleChange * pParticleChange
Definition: G4VProcess.hh:283
G4NuclearDecayChannel::GetDecayMode
G4RadioactiveDecayMode GetDecayMode()
Definition: G4NuclearDecayChannel.hh:100
G4RadioactiveDecay::AddDecayRateTable
void AddDecayRateTable(const G4ParticleDefinition &)
Definition: G4RadioactiveDecay.cc:1123
FatalException
Definition: G4ExceptionSeverity.hh:60
G4HadronicProcessType.hh
DBL_MIN
#define DBL_MIN
Definition: templates.hh:75
G4RadioactiveDecay::LoadDecayTable
G4DecayTable * LoadDecayTable(const G4ParticleDefinition &theParentNucleus)
Definition: G4RadioactiveDecay.cc:699
G4Neutron::Definition
static G4Neutron * Definition()
Definition: G4Neutron.cc:54
fStopAndKill
Definition: G4TrackStatus.hh:56
G4Track::GetVolume
G4VPhysicalVolume * GetVolume() const
G4RadioactiveDecay::theRadioactivityTables
std::vector< G4RadioactivityTable * > theRadioactivityTables
Definition: G4RadioactiveDecay.hh:303
G4Track::GetWeight
G4double GetWeight() const
G4ParticleDefinition::GetPDGLifeTime
G4double GetPDGLifeTime() const
G4RadioactiveDecay::G4RadioactiveDecay
G4RadioactiveDecay(const G4String &processName="RadioactiveDecay")
Definition: G4RadioactiveDecay.cc:151
G4endl
#define G4endl
Definition: G4ios.hh:61
m
static const double m
Definition: G4SIunits.hh:110
G4RadioactiveDecay::forceDecayDirection
G4ThreeVector forceDecayDirection
Definition: G4RadioactiveDecay.hh:283
G4DecayProducts::IsChecked
G4bool IsChecked() const
Definition: G4DecayProducts.cc:241
G4LossTableManager::AtomDeexcitation
G4VAtomDeexcitation * AtomDeexcitation()
Definition: G4LossTableManager.cc:1198
keV
static const double keV
Definition: G4SIunits.hh:195
G4AlphaDecay::SetHLThreshold
void SetHLThreshold(G4double HLT)
Definition: G4AlphaDecay.hh:55
G4DecayProducts::entries
G4int entries() const
Definition: G4DecayProducts.hh:86
G4VDecayChannel::DecayIt
virtual G4DecayProducts * DecayIt(G4double parentMass=-1.0)=0
G4LogicalVolumeStore.hh
G4double
double G4double
Definition: G4Types.hh:76
G4VParticleChange::ProposeTrackStatus
void ProposeTrackStatus(G4TrackStatus status)
IT
Definition: G4RadioactiveDecayMode.hh:65
SpFission
Definition: G4RadioactiveDecayMode.hh:66
G4RadioactiveDecay::BuildPhysicsTable
void BuildPhysicsTable(const G4ParticleDefinition &)
Definition: G4RadioactiveDecay.cc:676
Beta2Minus
Definition: G4RadioactiveDecayMode.hh:66
G4BetaDecayCorrections::ShapeFactor
G4double ShapeFactor(const G4BetaDecayType &, const G4double &p_e, const G4double &e_nu)
Definition: G4BetaDecayCorrections.cc:126
G4SystemOfUnits.hh
G4ThreeVector.hh
G4BetaPlusDecayChannel.hh
G4ForceCondition
G4ForceCondition
Definition: G4ForceCondition.hh:49
G4NuclearDecayChannel::SetHLThreshold
void SetHLThreshold(G4double hl)
Definition: G4NuclearDecayChannel.hh:92
G4RadioactiveDecay::SetAnalogueMonteCarlo
void SetAnalogueMonteCarlo(G4bool r)
Definition: G4RadioactiveDecay.hh:184
alpha
static const G4double alpha
Definition: G4LivermoreBremsstrahlungModel.cc:79
G4KshellECDecayChannel.hh
G4BetaDecayCorrections::FermiFunction
G4double FermiFunction(const G4double &W)
Definition: G4BetaDecayCorrections.cc:57
G4RadioactiveDecay::DeselectAVolume
void DeselectAVolume(const G4String aVolume)
Definition: G4RadioactiveDecay.cc:288
DBL_MAX
#define DBL_MAX
Definition: templates.hh:83
fDecay
Definition: G4ProcessType.hh:51
RDM_ERROR
Definition: G4RadioactiveDecayMode.hh:65
G4DecayTable.hh
G4LshellECDecayChannel.hh
ns
#define ns
Definition: xmlparse.cc:597
G4NuclearLevelManager::NearestLevel
const G4NuclearLevel * NearestLevel(G4double energy, G4double eDiffMax=1.e+8) const
Definition: G4NuclearLevelManager.cc:118
G4AutoLock.hh
G4RadioactiveDecay::theIsotopeTable
G4RIsotopeTable * theIsotopeTable
Definition: G4RadioactiveDecay.hh:268
G4RadioactiveDecay::AddUserDecayDataFile
void AddUserDecayDataFile(G4int Z, G4int A, G4String filename)
Definition: G4RadioactiveDecay.cc:1092
G4MshellECDecayChannel.hh
MshellEC
Definition: G4RadioactiveDecayMode.hh:65
G4RadioactiveDecay::CollimateDecay
void CollimateDecay(G4DecayProducts *products)
Definition: G4RadioactiveDecay.cc:1870
G4RIsotopeTable.hh
G4RadioactiveDecay::SetDecayBias
void SetDecayBias(G4String filename)
Definition: G4RadioactiveDecay.cc:1448
G4NucleusLimits::GetZMin
G4int GetZMin() const
Definition: G4NucleusLimits.hh:97
G4Gamma::Definition
static G4Gamma * Definition()
Definition: G4Gamma.cc:49
G4Track::SetGoodForTrackingFlag
void SetGoodForTrackingFlag(G4bool value=true)
G4RadioactiveDecay::DecayIt
G4VParticleChange * DecayIt(const G4Track &theTrack, const G4Step &theStep)
Definition: G4RadioactiveDecay.cc:1498
G4Gamma.hh
G4AlphaDecay.hh
G4cerr
G4GLOB_DLL std::ostream G4cerr
G4RadioactiveDecay.hh