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