ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4Radioactivation.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4Radioactivation.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 //
27 // //
28 // File: G4Radioactivation.cc //
29 // Author: D.H. Wright (SLAC) //
30 // Date: 29 August 2017 //
31 // Description: activation process derived from the original //
32 // G4RadioactiveDecay of F. Lei and P.R. Truscott in which //
33 // biasing and activation calculations are separated from the //
34 // unbiased decay chain calculation performed in the base //
35 // class. //
36 // //
38 
39 #include "G4Radioactivation.hh"
41 
42 #include "G4SystemOfUnits.hh"
43 #include "G4DynamicParticle.hh"
44 #include "G4DecayProducts.hh"
45 #include "G4DecayTable.hh"
47 #include "G4ITDecay.hh"
48 #include "G4BetaDecayType.hh"
49 #include "G4BetaMinusDecay.hh"
50 #include "G4BetaPlusDecay.hh"
51 #include "G4ECDecay.hh"
52 #include "G4AlphaDecay.hh"
53 #include "G4TritonDecay.hh"
54 #include "G4ProtonDecay.hh"
55 #include "G4NeutronDecay.hh"
56 #include "G4SFDecay.hh"
57 #include "G4VDecayChannel.hh"
58 #include "G4NuclearDecay.hh"
60 #include "G4Fragment.hh"
61 #include "G4Ions.hh"
62 #include "G4IonTable.hh"
63 #include "G4BetaDecayType.hh"
64 #include "Randomize.hh"
65 #include "G4LogicalVolumeStore.hh"
66 #include "G4NuclearLevelData.hh"
67 #include "G4DeexPrecoParameters.hh"
68 #include "G4LevelManager.hh"
69 #include "G4ThreeVector.hh"
70 #include "G4Electron.hh"
71 #include "G4Positron.hh"
72 #include "G4Neutron.hh"
73 #include "G4Gamma.hh"
74 #include "G4Alpha.hh"
75 #include "G4Triton.hh"
76 #include "G4Proton.hh"
77 
78 #include "G4HadronicProcessType.hh"
80 #include "G4HadronicException.hh"
81 #include "G4LossTableManager.hh"
82 #include "G4VAtomDeexcitation.hh"
83 #include "G4UAtomicDeexcitation.hh"
84 #include "G4PhotonEvaporation.hh"
85 
86 #include <vector>
87 #include <sstream>
88 #include <algorithm>
89 #include <fstream>
90 
91 using namespace CLHEP;
92 
94  : G4RadioactiveDecayBase(processName)
95 {
96 #ifdef G4VERBOSE
97  if (GetVerboseLevel() > 1) {
98  G4cout << "G4Radioactivation constructor: processName = " << processName
99  << G4endl;
100  }
101 #endif
102 
103 // DHW SetProcessSubType(fRadioactiveDecay);
105 
106  // Apply default values.
107  NSourceBin = 1;
108  SBin[0] = 0.* s;
109  SBin[1] = 1.* s; // Convert to ns
110  SProfile[0] = 1.;
111  SProfile[1] = 0.;
112  NDecayBin = 1;
113  DBin[0] = 0. * s ;
114  DBin[1] = 1. * s;
115  DProfile[0] = 1.;
116  DProfile[1] = 0.;
117  decayWindows[0] = 0;
119  theRadioactivityTables.push_back(rTable);
120  NSplit = 1;
121  AnalogueMC = true;
122  BRBias = true;
124 }
125 
126 
127 void G4Radioactivation::ProcessDescription(std::ostream& outFile) const
128 {
129  outFile << "The G4Radioactivation process performs radioactive decay of\n"
130  << "nuclides (G4GenericIon) in biased mode which includes nucleus\n"
131  << "duplication, branching ratio biasing, source time convolution\n"
132  << "and detector time convolution. It is designed for use in\n"
133  << "activation physics.\n"
134  << "The required half-lives and decay schemes are retrieved from\n"
135  << "the RadioactiveDecay database which was derived from ENSDF.\n";
136 }
137 
138 
140 {
142 }
143 
145 {
146  G4String key = aNucleus->GetParticleName();
147  DecayTableMap::iterator table_ptr = dkmap->find(key);
148 
149  G4DecayTable* theDecayTable = 0;
150  if (table_ptr == dkmap->end() ) { // If table not there,
151  theDecayTable = LoadDecayTable(*aNucleus); // load from file and
152  if(theDecayTable) (*dkmap)[key] = theDecayTable; // store in library
153  } else {
154  theDecayTable = table_ptr->second;
155  }
156  return theDecayTable;
157 }
158 
159 G4bool
161 {
162  // Check whether the radioactive decay rates table for the ion has already
163  // been calculated.
164  G4String aParticleName = aParticle.GetParticleName();
165  for (size_t i = 0; i < theParentChainTable.size(); i++) {
166  if (theParentChainTable[i].GetIonName() == aParticleName) return true;
167  }
168  return false;
169 }
170 
171 
172 void
174 {
175  // Retrieve the decay rate table for the specified aParticle
176  G4String aParticleName = aParticle.GetParticleName();
177 
178  for (size_t i = 0; i < theParentChainTable.size(); i++) {
179  if (theParentChainTable[i].GetIonName() == aParticleName) {
180  theDecayRateVector = theParentChainTable[i].GetItsRates();
181  }
182  }
183 #ifdef G4VERBOSE
184  if (GetVerboseLevel() > 1) {
185  G4cout << "The DecayRate Table for " << aParticleName << " is selected."
186  << G4endl;
187  }
188 #endif
189 }
190 
191 // ConvolveSourceTimeProfile performs the convolution of the source time profile
192 // function with a single exponential characterized by a decay constant in the
193 // decay chain. The time profile is treated as a step function so that the
194 // convolution integral can be done bin-by-bin.
195 // This implements Eq. 4.13 of DERA technical note, with SProfile[i] = F(t')
196 
197 G4double
199 {
200  G4double convolvedTime = 0.0;
201  G4int nbin;
202  if ( t > SBin[NSourceBin]) {
203  nbin = NSourceBin;
204  } else {
205  nbin = 0;
206 
207  G4int loop = 0;
208  while (t > SBin[nbin]) { // Loop checking, 01.09.2015, D.Wright
209  loop++;
210  if (loop > 1000) {
211  G4Exception("G4Radioactivation::ConvolveSourceTimeProfile()",
212  "HAD_RDM_100", JustWarning, "While loop count exceeded");
213  break;
214  }
215  nbin++;
216  }
217  nbin--;
218  }
219 
220  // Use expm1 wherever possible to avoid large cancellation errors in
221  // 1 - exp(x) for small x
222  G4double earg = 0.0;
223  if (nbin > 0) {
224  for (G4int i = 0; i < nbin; i++) {
225  earg = (SBin[i+1] - SBin[i])/tau;
226  if (earg < 100.) {
227  convolvedTime += SProfile[i] * std::exp((SBin[i] - t)/tau) *
228  std::expm1(earg);
229  } else {
230  convolvedTime += SProfile[i] *
231  (std::exp(-(t-SBin[i+1])/tau)-std::exp(-(t-SBin[i])/tau));
232  }
233  }
234  }
235  convolvedTime -= SProfile[nbin] * std::expm1((SBin[nbin] - t)/tau);
236  // tau divided out of final result to provide probability of decay in window
237 
238  if (convolvedTime < 0.) {
239  G4cout << " Convolved time =: " << convolvedTime << " reset to zero! " << G4endl;
240  G4cout << " t = " << t << " tau = " << tau << G4endl;
241  G4cout << SBin[nbin] << " " << SBin[0] << G4endl;
242  convolvedTime = 0.;
243  }
244 #ifdef G4VERBOSE
245  if (GetVerboseLevel() > 2)
246  G4cout << " Convolved time: " << convolvedTime << G4endl;
247 #endif
248  return convolvedTime;
249 }
250 
251 
253 // //
254 // GetDecayTime //
255 // Randomly select a decay time for the decay process, following the //
256 // supplied decay time bias scheme. //
257 // //
259 
261 {
262  G4double decaytime = 0.;
263  G4double rand = G4UniformRand();
264  G4int i = 0;
265 
266  G4int loop = 0;
267  while (DProfile[i] < rand) { /* Loop checking, 01.09.2015, D.Wright */
268  // Entries in DProfile[i] are all between 0 and 1 and arranged in inreaseing order
269  // Comparison with rand chooses which time bin to sample
270  i++;
271  loop++;
272  if (loop > 100000) {
273  G4Exception("G4Radioactivation::GetDecayTime()", "HAD_RDM_100",
274  JustWarning, "While loop count exceeded");
275  break;
276  }
277  }
278 
279  rand = G4UniformRand();
280  decaytime = DBin[i] + rand*(DBin[i+1]-DBin[i]);
281 #ifdef G4VERBOSE
282  if (GetVerboseLevel() > 2)
283  G4cout <<" Decay time: " <<decaytime/s <<"[s]" <<G4endl;
284 #endif
285  return decaytime;
286 }
287 
288 
290 {
291  G4int i = 0;
292 
293  G4int loop = 0;
294  while (aDecayTime > DBin[i] ) { /* Loop checking, 01.09.2015, D.Wright */
295  i++;
296  loop++;
297  if (loop > 100000) {
298  G4Exception("G4Radioactivation::GetDecayTimeBin()", "HAD_RDM_100",
299  JustWarning, "While loop count exceeded");
300  break;
301  }
302  }
303 
304  return i;
305 }
306 
308 // //
309 // GetMeanLifeTime (required by the base class) //
310 // //
312 
315 {
316  // For variance reduction time is set to 0 so as to force the particle
317  // to decay immediately.
318  // In analogue mode it returns the particle's mean-life.
319  G4double meanlife = 0.;
320  if (AnalogueMC) meanlife = G4RadioactiveDecayBase::GetMeanLifeTime(theTrack, 0);
321  return meanlife;
322 }
323 
324 
325 void
327  G4int theG, std::vector<G4double> theCoefficients,
328  std::vector<G4double> theTaos)
329 // Why not make this a method of G4RadioactiveDecayRate? (e.g. SetParameters)
330 {
331  //fill the decay rate vector
332  ratesToDaughter.SetZ(theZ);
333  ratesToDaughter.SetA(theA);
334  ratesToDaughter.SetE(theE);
336  ratesToDaughter.SetDecayRateC(theCoefficients);
337  ratesToDaughter.SetTaos(theTaos);
338 }
339 
340 
343 {
344  // Use extended Bateman equation to calculate the radioactivities of all
345  // progeny of theParentNucleus. The coefficients required to do this are
346  // calculated using the method of P. Truscott (Ph.D. thesis and
347  // DERA Technical Note DERA/CIS/CIS2/7/36/4/10) 11 January 2000.
348  // Coefficients are then added to the decay rate table vector
349 
350  // Create and initialise variables used in the method.
351  theDecayRateVector.clear();
352 
353  G4int nGeneration = 0;
354 
355  std::vector<G4double> taos;
356 
357  // Dimensionless A coefficients of Eqs. 4.24 and 4.25 of the TN
358  std::vector<G4double> Acoeffs;
359 
360  // According to Eq. 4.26 the first coefficient (A_1:1) is -1
361  Acoeffs.push_back(-1.);
362 
363  G4int A = ((const G4Ions*)(&theParentNucleus))->GetAtomicMass();
364  G4int Z = ((const G4Ions*)(&theParentNucleus))->GetAtomicNumber();
365  G4double E = ((const G4Ions*)(&theParentNucleus))->GetExcitationEnergy();
366  G4double tao = theParentNucleus.GetPDGLifeTime();
367  if (tao < 0.) tao = 1e-100;
368  taos.push_back(tao);
369  G4int nEntry = 0;
370 
371  // Fill the decay rate container (G4RadioactiveDecayRate) with the parent
372  // isotope data
373  SetDecayRate(Z,A,E,nGeneration,Acoeffs,taos); // Fill TP with parent lifetime
374 
375  // store the decay rate in decay rate vector
377  nEntry++;
378 
379  // Now start treating the secondary generations.
380  G4bool stable = false;
381 // G4int i;
382  G4int j;
383  G4VDecayChannel* theChannel = 0;
384  G4NuclearDecay* theNuclearDecayChannel = 0;
385 
386  G4ITDecay* theITChannel = 0;
387  G4BetaMinusDecay* theBetaMinusChannel = 0;
388  G4BetaPlusDecay* theBetaPlusChannel = 0;
389  G4AlphaDecay* theAlphaChannel = 0;
390  G4ProtonDecay* theProtonChannel = 0;
391  G4TritonDecay* theTritonChannel = 0;
392  G4NeutronDecay* theNeutronChannel = 0;
393  G4SFDecay* theFissionChannel = 0;
394 
395  G4RadioactiveDecayMode theDecayMode;
396  G4double theBR = 0.0;
397  G4int AP = 0;
398  G4int ZP = 0;
399  G4int AD = 0;
400  G4int ZD = 0;
401  G4double EP = 0.;
402  std::vector<G4double> TP;
403  std::vector<G4double> RP; // A coefficients of the previous generation
404  G4ParticleDefinition *theDaughterNucleus;
405  G4double daughterExcitation;
406  G4double nearestEnergy = 0.0;
407  G4int nearestLevelIndex = 0;
408  G4ParticleDefinition *aParentNucleus;
409  G4IonTable* theIonTable;
410  G4DecayTable* parentDecayTable;
411  G4double theRate;
412  G4double TaoPlus;
413  G4int nS = 0; // Running index of first decay in a given generation
414  G4int nT = nEntry; // Total number of decays accumulated over entire history
415  const G4int nMode = 12;
416  G4double brs[nMode];
417  //
418  theIonTable =
420 
421  G4int loop = 0;
422  while (!stable) { /* Loop checking, 01.09.2015, D.Wright */
423  loop++;
424  if (loop > 10000) {
425  G4Exception("G4Radioactivation::CalculateChainsFromParent()", "HAD_RDM_100",
426  JustWarning, "While loop count exceeded");
427  break;
428  }
429  nGeneration++;
430  for (j = nS; j < nT; j++) {
431  // First time through, get data for parent nuclide
432  ZP = theDecayRateVector[j].GetZ();
433  AP = theDecayRateVector[j].GetA();
434  EP = theDecayRateVector[j].GetE();
435  RP = theDecayRateVector[j].GetDecayRateC();
436  TP = theDecayRateVector[j].GetTaos();
437  if (GetVerboseLevel() > 1) {
438  G4cout << "G4RadioactiveDecay::CalculateChainsFromParent: daughters of ("
439  << ZP << ", " << AP << ", " << EP
440  << ") are being calculated, generation = " << nGeneration
441  << G4endl;
442  }
443 // G4cout << " Taus = " << G4endl;
444 // for (G4int ii = 0; ii < TP.size(); ii++) G4cout << TP[ii] << ", " ;
445 // G4cout << G4endl;
446 
447  aParentNucleus = theIonTable->GetIon(ZP,AP,EP);
448  parentDecayTable = GetDecayTable1(aParentNucleus);
449 
450  G4DecayTable* summedDecayTable = new G4DecayTable();
451  // This instance of G4DecayTable is for accumulating BRs and decay
452  // channels. It will contain one decay channel per type of decay
453  // (alpha, beta, etc.); its branching ratio will be the sum of all
454  // branching ratios for that type of decay of the parent. If the
455  // halflife of a particular channel is longer than some threshold,
456  // that channel will be inserted specifically and its branching
457  // ratio will not be included in the above sums.
458  // This instance is not used to perform actual decays.
459 
460  for (G4int k = 0; k < nMode; k++) brs[k] = 0.0;
461 
462  // Go through the decay table and sum all channels having the same decay mode
463  for (G4int i = 0; i < parentDecayTable->entries(); i++) {
464  theChannel = parentDecayTable->GetDecayChannel(i);
465  theNuclearDecayChannel = static_cast<G4NuclearDecay*>(theChannel);
466  theDecayMode = theNuclearDecayChannel->GetDecayMode();
467  daughterExcitation = theNuclearDecayChannel->GetDaughterExcitation();
468  theDaughterNucleus = theNuclearDecayChannel->GetDaughterNucleus() ;
469  AD = ((const G4Ions*)(theDaughterNucleus))->GetAtomicMass();
470  ZD = ((const G4Ions*)(theDaughterNucleus))->GetAtomicNumber();
471  const G4LevelManager* levelManager =
473 
474  // Check each nuclide to see if it is metastable (lifetime > 1 usec)
475  // If so, add it to the decay chain by inserting its decay channel in
476  // summedDecayTable. If not, just add its BR to sum for that decay mode.
477  if (levelManager->NumberOfTransitions() ) {
478  nearestEnergy = levelManager->NearestLevelEnergy(daughterExcitation);
479  if (std::abs(daughterExcitation - nearestEnergy) < levelTolerance) {
480  // Level half-life is in ns and the threshold is set to 1 micros
481  // by default, user can set it via the UI command
482  nearestLevelIndex = levelManager->NearestLevelIndex(daughterExcitation);
483  if (levelManager->LifeTime(nearestLevelIndex)*ns >= halflifethreshold){
484  // save the metastable decay channel
485  summedDecayTable->Insert(theChannel);
486  } else {
487  brs[theDecayMode] += theChannel->GetBR();
488  }
489  } else {
490  brs[theDecayMode] += theChannel->GetBR();
491  }
492  } else {
493  brs[theDecayMode] += theChannel->GetBR();
494  }
495 
496  } // Combine decay channels (loop i)
497 
498  brs[2] = brs[2]+brs[3]+brs[4]+brs[5]+brs[6]; // Combine beta+ and EC
499  brs[3] = brs[4] = brs[5] = brs[6] = 0.0;
500  for (G4int i = 0; i < nMode; i++) { // loop over decay modes
501  if (brs[i] > 0.) {
502  switch (i) {
503  case 0:
504  // Decay mode is isomeric transition
505  theITChannel = new G4ITDecay(aParentNucleus, brs[0], 0.0, 0.0,
507 
508  summedDecayTable->Insert(theITChannel);
509  break;
510 
511  case 1:
512  // Decay mode is beta-
513  theBetaMinusChannel = new G4BetaMinusDecay(aParentNucleus, brs[1],
514  0.*MeV, 0.*MeV,
515  noFloat, allowed);
516  summedDecayTable->Insert(theBetaMinusChannel);
517  break;
518 
519  case 2:
520  // Decay mode is beta+ + EC.
521  theBetaPlusChannel = new G4BetaPlusDecay(aParentNucleus, brs[2],
522  0.*MeV, 0.*MeV,
523  noFloat, allowed);
524  summedDecayTable->Insert(theBetaPlusChannel);
525  break;
526 
527  case 7:
528  // Decay mode is alpha.
529  theAlphaChannel = new G4AlphaDecay(aParentNucleus, brs[7], 0.*MeV,
530  0.*MeV, noFloat);
531  summedDecayTable->Insert(theAlphaChannel);
532  break;
533 
534  case 8:
535  // Decay mode is proton.
536  theProtonChannel = new G4ProtonDecay(aParentNucleus, brs[8], 0.*MeV,
537  0.*MeV, noFloat);
538  summedDecayTable->Insert(theProtonChannel);
539  break;
540 
541  case 9:
542  // Decay mode is neutron.
543  theNeutronChannel = new G4NeutronDecay(aParentNucleus, brs[9], 0.*MeV,
544  0.*MeV, noFloat);
545  summedDecayTable->Insert(theNeutronChannel);
546  break;
547 
548  case 10:
549  // Decay mode is spontaneous fission
550  theFissionChannel = new G4SFDecay(aParentNucleus, brs[10], 0.*MeV,
551  0.*MeV, noFloat);
552  summedDecayTable->Insert(theFissionChannel);
553  break;
554  case 11:
555  // Decay mode is Triton.
556  theTritonChannel = new G4TritonDecay(aParentNucleus, brs[9], 0.*MeV,
557  0.*MeV, noFloat);
558  summedDecayTable->Insert(theTritonChannel);
559  break;
560 
561  default:
562  break;
563  }
564  }
565  }
566 
567  // loop over all branches in summedDecayTable
568  //
569  for (G4int i = 0; i < summedDecayTable->entries(); i++){
570  theChannel = summedDecayTable->GetDecayChannel(i);
571  theNuclearDecayChannel = static_cast<G4NuclearDecay*>(theChannel);
572  theBR = theChannel->GetBR();
573  theDaughterNucleus = theNuclearDecayChannel->GetDaughterNucleus();
574 
575  // First check if the decay of the original nucleus is an IT channel,
576  // if true create a new ground-state nucleus
577  if (theNuclearDecayChannel->GetDecayMode() == IT && nGeneration == 1) {
578 
579  A = ((const G4Ions*)(theDaughterNucleus))->GetAtomicMass();
580  Z = ((const G4Ions*)(theDaughterNucleus))->GetAtomicNumber();
581  theDaughterNucleus=theIonTable->GetIon(Z,A,0.);
582  }
583  if (IsApplicable(*theDaughterNucleus) && theBR > 0.0 &&
584  aParentNucleus != theDaughterNucleus) {
585  // need to make sure daughter has decay table
586  parentDecayTable = GetDecayTable1(theDaughterNucleus);
587  if (parentDecayTable->entries() ) {
588  A = ((const G4Ions*)(theDaughterNucleus))->GetAtomicMass();
589  Z = ((const G4Ions*)(theDaughterNucleus))->GetAtomicNumber();
590  E = ((const G4Ions*)(theDaughterNucleus))->GetExcitationEnergy();
591 
592  TaoPlus = theDaughterNucleus->GetPDGLifeTime();
593  if (TaoPlus <= 0.) TaoPlus = 1e-100;
594 
595  // first set the taos, one simply need to add to the parent ones
596  taos.clear();
597  taos = TP; // load lifetimes of all previous generations
598  size_t k;
599  //check that TaoPlus differs from other taos from at least 1.e5 relative difference
600  //for (k = 0; k < TP.size(); k++){
601  //if (std::abs((TaoPlus-TP[k])/TP[k])<1.e-5 ) TaoPlus=1.00001*TP[k];
602  //}
603  taos.push_back(TaoPlus); // add daughter lifetime to list
604  // now calculate the coefficiencies
605  //
606  // they are in two parts, first the less than n ones
607  // Eq 4.24 of the TN
608  Acoeffs.clear();
609  long double ta1,ta2;
610  ta2 = (long double)TaoPlus;
611  for (k = 0; k < RP.size(); k++){
612  ta1 = (long double)TP[k]; // loop over lifetimes of all previous generations
613  if (ta1 == ta2) {
614  theRate = 1.e100;
615  } else {
616  theRate = ta1/(ta1-ta2);
617  }
618  theRate = theRate * theBR * RP[k];
619  Acoeffs.push_back(theRate);
620  }
621 
622  // the second part: the n:n coefficiency
623  // Eq 4.25 of the TN. Note Yn+1 is zero apart from Y1 which is -1
624  // as treated at line 1013
625  theRate = 0.;
626  long double aRate, aRate1;
627  aRate1 = 0.L;
628  for (k = 0; k < RP.size(); k++){
629  ta1 = (long double)TP[k];
630  if (ta1 == ta2 ) {
631  aRate = 1.e100;
632  } else {
633  aRate = ta2/(ta1-ta2);
634  }
635  aRate = aRate * (long double)(theBR * RP[k]);
636  aRate1 += aRate;
637  }
638  theRate = -aRate1;
639  Acoeffs.push_back(theRate);
640  SetDecayRate (Z,A,E,nGeneration,Acoeffs,taos);
642  nEntry++;
643  } // there are entries in the table
644  } // nuclide is OK to decay
645  } // end of loop (i) over decay table branches
646 
647  delete summedDecayTable;
648 
649  } // Getting contents of decay rate vector (end loop on j)
650  nS = nT;
651  nT = nEntry;
652  if (nS == nT) stable = true;
653  } // while nuclide is not stable
654 
655  // end of while loop
656  // the calculation completed here
657 
658 
659  // fill the first part of the decay rate table
660  // which is the name of the original particle (isotope)
661  chainsFromParent.SetIonName(theParentNucleus.GetParticleName());
662 
663  // now fill the decay table with the newly completed decay rate vector
665 
666  // finally add the decayratetable to the tablevector
668 }
669 
671 // //
672 // SetSourceTimeProfile //
673 // read in the source time profile function (histogram) //
674 // //
676 
678 {
679  std::ifstream infile ( filename, std::ios::in );
680  if (!infile) {
682  ed << " Could not open file " << filename << G4endl;
683  G4Exception("G4Radioactivation::SetSourceTimeProfile()", "HAD_RDM_001",
684  FatalException, ed);
685  }
686 
687  G4double bin, flux;
688  NSourceBin = -1;
689 
690  G4int loop = 0;
691  while (infile >> bin >> flux) { /* Loop checking, 01.09.2015, D.Wright */
692  loop++;
693  if (loop > 10000) {
694  G4Exception("G4Radioactivation::SetSourceTimeProfile()", "HAD_RDM_100",
695  JustWarning, "While loop count exceeded");
696  break;
697  }
698 
699  NSourceBin++;
700  if (NSourceBin > 99) {
701  G4Exception("G4Radioactivation::SetSourceTimeProfile()", "HAD_RDM_002",
702  FatalException, "Input source time file too big (>100 rows)");
703 
704  } else {
705  SBin[NSourceBin] = bin * s; // Convert read-in time to ns
706  SProfile[NSourceBin] = flux; // Dimensionless
707  }
708  }
709 
710  AnalogueMC = false;
711  infile.close();
712 
713 #ifdef G4VERBOSE
714  if (GetVerboseLevel() > 2)
715  G4cout <<" Source Timeprofile Nbin = " << NSourceBin <<G4endl;
716 #endif
717 }
718 
720 // //
721 // SetDecayBiasProfile //
722 // read in the decay bias scheme function (histogram) //
723 // //
725 
727 {
728  std::ifstream infile(filename, std::ios::in);
729  if (!infile) G4Exception("G4Radioactivation::SetDecayBias()", "HAD_RDM_001",
730  FatalException, "Unable to open bias data file" );
731 
732  G4double bin, flux;
733  G4int dWindows = 0;
734  G4int i ;
735 
736  theRadioactivityTables.clear();
737 
738  NDecayBin = -1;
739 
740  G4int loop = 0;
741  while (infile >> bin >> flux ) { /* Loop checking, 01.09.2015, D.Wright */
742  NDecayBin++;
743  loop++;
744  if (loop > 10000) {
745  G4Exception("G4Radioactivation::SetDecayBias()", "HAD_RDM_100",
746  JustWarning, "While loop count exceeded");
747  break;
748  }
749 
750  if (NDecayBin > 99) {
751  G4Exception("G4Radioactivation::SetDecayBias()", "HAD_RDM_002",
752  FatalException, "Input bias file too big (>100 rows)" );
753  } else {
754  DBin[NDecayBin] = bin * s; // Convert read-in time to ns
755  DProfile[NDecayBin] = flux; // Dimensionless
756  if (flux > 0.) {
757  decayWindows[NDecayBin] = dWindows;
758  dWindows++;
760  theRadioactivityTables.push_back(rTable);
761  }
762  }
763  }
764  for ( i = 1; i<= NDecayBin; i++) DProfile[i] += DProfile[i-1]; // Cumulative flux vs i
765  for ( i = 0; i<= NDecayBin; i++) DProfile[i] /= DProfile[NDecayBin];
766  // Normalize so entries increase from 0 to 1
767  // converted to accumulated probabilities
768 
769  AnalogueMC = false;
770  infile.close();
771 
772 #ifdef G4VERBOSE
773  if (GetVerboseLevel() > 2)
774  G4cout <<" Decay Bias Profile Nbin = " << NDecayBin <<G4endl;
775 #endif
776 }
777 
779 // //
780 // DecayIt //
781 // //
783 
785 G4Radioactivation::DecayIt(const G4Track& theTrack, const G4Step&)
786 {
787  // Initialize G4ParticleChange object, get particle details and decay table
790  const G4DynamicParticle* theParticle = theTrack.GetDynamicParticle();
791  const G4ParticleDefinition* theParticleDef = theParticle->GetDefinition();
792 
793  // First check whether RDM applies to the current logical volume
794  if (!isAllVolumesMode) {
795  if (!std::binary_search(ValidVolumes.begin(), ValidVolumes.end(),
796  theTrack.GetVolume()->GetLogicalVolume()->GetName())) {
797 #ifdef G4VERBOSE
798  if (GetVerboseLevel()>0) {
799  G4cout <<"G4RadioactiveDecay::DecayIt : "
800  << theTrack.GetVolume()->GetLogicalVolume()->GetName()
801  << " is not selected for the RDM"<< G4endl;
802  G4cout << " There are " << ValidVolumes.size() << " volumes" << G4endl;
803  G4cout << " The Valid volumes are " << G4endl;
804  for (size_t i = 0; i< ValidVolumes.size(); i++)
805  G4cout << ValidVolumes[i] << G4endl;
806  }
807 #endif
809 
810  // Kill the parent particle.
815  }
816  }
817 
818  // Now check if particle is valid for RDM
819  if (!(IsApplicable(*theParticleDef) ) ) {
820  // Particle is not an ion or is outside the nucleuslimits for decay
821 #ifdef G4VERBOSE
822  if (GetVerboseLevel() > 1) {
823  G4cout << "G4RadioactiveDecay::DecayIt : "
824  << theParticleDef->GetParticleName()
825  << " is not an ion or is outside (Z,A) limits set for the decay. "
826  << " Set particle change accordingly. "
827  << G4endl;
828  }
829 #endif
831 
832  // Kill the parent particle
837  }
838 
839  G4DecayTable* theDecayTable = GetDecayTable1(theParticleDef);
840 
841  if (theDecayTable == 0 || theDecayTable->entries() == 0) {
842  // No data in the decay table. Set particle change parameters
843  // to indicate this.
844 #ifdef G4VERBOSE
845  if (GetVerboseLevel() > 1) {
846  G4cout << "G4RadioactiveDecay::DecayIt : "
847  << "decay table not defined for "
848  << theParticleDef->GetParticleName()
849  << ". Set particle change accordingly. "
850  << G4endl;
851  }
852 #endif
854 
855  // Kill the parent particle.
860 
861  } else {
862  // Data found. Try to decay nucleus
863  if (AnalogueMC) {
865 
866  } else {
867  // Proceed with decay using variance reduction
868  G4double energyDeposit = 0.0;
869  G4double finalGlobalTime = theTrack.GetGlobalTime();
870  G4double finalLocalTime = theTrack.GetLocalTime();
871  G4int index;
872  G4ThreeVector currentPosition;
873  currentPosition = theTrack.GetPosition();
874 
875  G4IonTable* theIonTable;
876  G4ParticleDefinition* parentNucleus;
877 
878  // Get decay chains for the given nuclide
879  if (!IsRateTableReady(*theParticleDef)) CalculateChainsFromParent(*theParticleDef);
880  GetChainsFromParent(*theParticleDef);
881 
882  // Declare some of the variables required in the implementation
883  G4int PZ;
884  G4int PA;
885  G4double PE;
886  G4String keyName;
887  std::vector<G4double> PT;
888  std::vector<G4double> PR;
889  G4double taotime;
890  long double decayRate;
891 
892  size_t i;
893  G4int numberOfSecondaries;
894  G4int totalNumberOfSecondaries = 0;
895  G4double currentTime = 0.;
896  G4int ndecaych;
897  G4DynamicParticle* asecondaryparticle;
898  std::vector<G4DynamicParticle*> secondaryparticles;
899  std::vector<G4double> pw;
900  std::vector<G4double> ptime;
901  pw.clear();
902  ptime.clear();
903 
904  // Now apply the nucleus splitting
905  for (G4int n = 0; n < NSplit; n++) {
906  // Get the decay time following the decay probability function
907  // supplied by user
908  G4double theDecayTime = GetDecayTime();
909  G4int nbin = GetDecayTimeBin(theDecayTime);
910 
911  // calculate the first part of the weight function
912  G4double weight1 = 1.;
913  if (nbin == 1) {
914  weight1 = 1./DProfile[nbin-1]
915  *(DBin[nbin]-DBin[nbin-1])/NSplit; // width of window in ns
916  } else if (nbin > 1) {
917  // Go from nbin to nbin-2 because flux entries in file alternate between 0 and 1
918  weight1 = 1./(DProfile[nbin]-DProfile[nbin-2])
919  *(DBin[nbin]-DBin[nbin-1])/NSplit;
920  // weight1 = (probability of choosing one of the bins)*(time width of bin)/NSplit
921  }
922  // it should be calculated in seconds
923  weight1 /= s ;
924 
925  // loop over all the possible secondaries of the nucleus
926  // the first one is itself.
927  for (i = 0; i < theDecayRateVector.size(); i++) {
928  PZ = theDecayRateVector[i].GetZ();
929  PA = theDecayRateVector[i].GetA();
930  PE = theDecayRateVector[i].GetE();
931  PT = theDecayRateVector[i].GetTaos();
932  PR = theDecayRateVector[i].GetDecayRateC();
933 
934  // The array of arrays theDecayRateVector contains all possible decay
935  // chains of a given parent nucleus (ZP,AP,EP) to a given descendant
936  // nuclide (Z,A,E).
937  //
938  // theDecayRateVector[0] contains the decay parameters of the parent
939  // nucleus
940  // PZ = ZP
941  // PA = AP
942  // PE = EP
943  // PT[] = {TP}
944  // PR[] = {RP}
945  //
946  // theDecayRateVector[1] contains the decay of the parent to the first
947  // generation daughter (Z1,A1,E1).
948  // PZ = Z1
949  // PA = A1
950  // PE = E1
951  // PT[] = {TP, T1}
952  // PR[] = {RP, R1}
953  //
954  // theDecayRateVector[2] contains the decay of the parent to the first
955  // generation daughter (Z1,A1,E1) and the decay of the first
956  // generation daughter to the second generation daughter (Z2,A2,E2).
957  // PZ = Z2
958  // PA = A2
959  // PE = E2
960  // PT[] = {TP, T1, T2}
961  // PR[] = {RP, R1, R2}
962  //
963  // theDecayRateVector[3] may contain a branch chain
964  // PZ = Z2a
965  // PA = A2a
966  // PE = E2a
967  // PT[] = {TP, T1, T2a}
968  // PR[] = {RP, R1, R2a}
969  //
970  // and so on.
971 
972  // Calculate the decay rate of the isotope. decayRate is the
973  // radioactivity of isotope (PZ,PA,PE) at 'theDecayTime'
974  // it will be used to calculate the statistical weight of the
975  // decay products of this isotope
976 
977  // For each nuclide, calculate all the decay chains which can reach
978  // the parent nuclide
979  decayRate = 0.L;
980  for (G4int j = 0; j < G4int(PT.size() ); j++) {
981  taotime = ConvolveSourceTimeProfile(theDecayTime,PT[j]);
982  decayRate -= PR[j] * (long double)taotime; // PRs are Acoeffs, taotime is inverse time
983  // Eq.4.23 of of the TN
984  // note the negative here is required as the rate in the
985  // equation is defined to be negative,
986  // i.e. decay away, but we need positive value here.
987 
988  // G4cout << j << "\t"<< PT[j]/s << "\t" << PR[j] << "\t" << decayRate << G4endl;
989  }
990 
991  // At this point any negative decay rates are probably small enough
992  // (order 10**-30) that negative values are likely due to cancellation
993  // errors. Set them to zero.
994  if (decayRate < 0.0) decayRate = 0.0;
995 
996  // G4cout <<theDecayTime/s <<"\t"<<nbin<<G4endl;
997  // G4cout << theTrack.GetWeight() <<"\t"<<weight1<<"\t"<<decayRate<< G4endl;
998 
999  // Add isotope to the radioactivity tables
1000  // One table for each observation time window specifed in
1001  // SetDecayBias(G4String filename)
1002 
1004  ->AddIsotope(PZ,PA,PE,weight1*decayRate,theTrack.GetWeight());
1005 
1006  // Now calculate the statistical weight
1007  // One needs to fold the source bias function with the decaytime
1008  // also need to include the track weight! (F.Lei, 28/10/10)
1009  G4double weight = weight1*decayRate*theTrack.GetWeight();
1010 
1011  // decay the isotope
1013  parentNucleus = theIonTable->GetIon(PZ,PA,PE);
1014 
1015  // Create a temprary products buffer.
1016  // Its contents to be transfered to the products at the end of the loop
1017  G4DecayProducts* tempprods = 0;
1018 
1019  // Decide whether to apply branching ratio bias or not
1020  if (BRBias) {
1021  G4DecayTable* decayTable = GetDecayTable1(parentNucleus);
1022  ndecaych = G4int(decayTable->entries()*G4UniformRand());
1023  G4VDecayChannel* theDecayChannel = decayTable->GetDecayChannel(ndecaych);
1024 
1025  if (theDecayChannel == 0) {
1026  // Decay channel not found.
1027 
1028  if (GetVerboseLevel() > 0) {
1029  G4cout << " G4RadioactiveDecay::DoIt : cannot determine decay channel ";
1030  G4cout << " for this nucleus; decay as if no biasing active. ";
1031  G4cout << G4endl;
1032  decayTable ->DumpInfo();
1033  }
1034 
1035  tempprods = DoDecay(*parentNucleus); // DHW 6 Dec 2010 - do decay as if no biasing
1036  // to avoid deref of temppprods = 0
1037  } else {
1038  // A decay channel has been identified, so execute the DecayIt.
1039  G4double tempmass = parentNucleus->GetPDGMass();
1040  tempprods = theDecayChannel->DecayIt(tempmass);
1041  weight *= (theDecayChannel->GetBR())*(decayTable->entries());
1042  }
1043  } else {
1044  tempprods = DoDecay(*parentNucleus);
1045  }
1046 
1047  // save the secondaries for buffers
1048  numberOfSecondaries = tempprods->entries();
1049  currentTime = finalGlobalTime + theDecayTime;
1050  for (index = 0; index < numberOfSecondaries; index++) {
1051  asecondaryparticle = tempprods->PopProducts();
1052  if (asecondaryparticle->GetDefinition()->GetPDGStable() ) {
1053  pw.push_back(weight);
1054  ptime.push_back(currentTime);
1055  secondaryparticles.push_back(asecondaryparticle);
1056  }
1057  // Generate gammas and Xrays from excited nucleus, added by L.Desorgher
1058  else if (((const G4Ions*)(asecondaryparticle->GetDefinition()))
1059  ->GetExcitationEnergy() > 0. && weight > 0.) { //Compute the gamma
1060  G4ParticleDefinition* apartDef = asecondaryparticle->GetDefinition();
1061  AddDeexcitationSpectrumForBiasMode(apartDef,weight,currentTime,pw,
1062  ptime,secondaryparticles);
1063  }
1064  }
1065 
1066  delete tempprods;
1067  } // end of i loop
1068  } // end of n loop
1069 
1070  // now deal with the secondaries in the two stl containers
1071  // and submmit them back to the tracking manager
1072  totalNumberOfSecondaries = pw.size();
1073  fParticleChangeForRadDecay.SetNumberOfSecondaries(totalNumberOfSecondaries);
1074  for (index=0; index < totalNumberOfSecondaries; index++) {
1075  G4Track* secondary = new G4Track(secondaryparticles[index],
1076  ptime[index], currentPosition);
1077  secondary->SetGoodForTrackingFlag();
1078  secondary->SetTouchableHandle(theTrack.GetTouchableHandle());
1079  secondary->SetWeight(pw[index]);
1081  }
1082 
1083  // Kill the parent particle
1087  // Reset NumberOfInteractionLengthLeft.
1089  } // end VR decay
1090 
1092  } // end of data found branch
1093 }
1094 
1095 
1096 // Add gamma, X-ray, conversion and auger electrons for bias mode
1097 void
1099  G4double weight,G4double currentTime,
1100  std::vector<double>& weights_v,
1101  std::vector<double>& times_v,
1102  std::vector<G4DynamicParticle*>& secondaries_v)
1103 {
1104  G4double elevel=((const G4Ions*)(apartDef))->GetExcitationEnergy();
1105  G4double life_time=apartDef->GetPDGLifeTime();
1106  G4ITDecay* anITChannel = 0;
1107 
1108  while (life_time < halflifethreshold && elevel > 0.) {
1109  anITChannel = new G4ITDecay(apartDef, 100., elevel, elevel, photonEvaporation);
1110  G4DecayProducts* pevap_products = anITChannel->DecayIt(0.);
1111  G4int nb_pevapSecondaries = pevap_products->entries();
1112 
1113  G4DynamicParticle* a_pevap_secondary = 0;
1114  G4ParticleDefinition* secDef = 0;
1115  for (G4int ind = 0; ind < nb_pevapSecondaries; ind++) {
1116  a_pevap_secondary= pevap_products->PopProducts();
1117  secDef = a_pevap_secondary->GetDefinition();
1118 
1119  if (secDef->GetBaryonNumber() > 4) {
1120  elevel = ((const G4Ions*)(secDef))->GetExcitationEnergy();
1121  life_time = secDef->GetPDGLifeTime();
1122  apartDef = secDef;
1123  if (secDef->GetPDGStable() ) {
1124  weights_v.push_back(weight);
1125  times_v.push_back(currentTime);
1126  secondaries_v.push_back(a_pevap_secondary);
1127  }
1128  } else {
1129  weights_v.push_back(weight);
1130  times_v.push_back(currentTime);
1131  secondaries_v.push_back(a_pevap_secondary);
1132  }
1133  }
1134 
1135  delete anITChannel;
1136  delete pevap_products;
1137  }
1138 }
1139