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