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