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