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