Geant4-11
G4Radioactivation.cc
Go to the documentation of this file.
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
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"
66#include "G4NuclearLevelData.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
81#include "G4LossTableManager.hh"
85
86#include <vector>
87#include <sstream>
88#include <algorithm>
89#include <fstream>
90
91using namespace CLHEP;
92
94 : G4RadioactiveDecay(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
127void 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
159G4bool
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
172void
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
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 = G4RadioactiveDecay::GetMeanLifeTime(theTrack, 0);
321 return meanlife;
322}
323
324
325void
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
342CalculateChainsFromParent(const G4ParticleDefinition& theParentNucleus)
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 = G4RadioactiveDecayModeSize;
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[BetaPlus] = brs[BetaPlus]+brs[KshellEC]+brs[LshellEC]+brs[MshellEC]+brs[NshellEC]; // Combine beta+ and EC
499 brs[KshellEC] = brs[LshellEC] = brs[MshellEC] = brs[NshellEC] = 0.0;
500 for (G4int i = 0; i < nMode; i++) { // loop over decay modes
501 if (brs[i] > 0.) {
502 switch (i) {
503 case IT:
504 // Decay mode is isomeric transition
505 theITChannel = new G4ITDecay(aParentNucleus, brs[IT], 0.0, 0.0,
507
508 summedDecayTable->Insert(theITChannel);
509 break;
510
511 case BetaMinus:
512 // Decay mode is beta-
513 theBetaMinusChannel = new G4BetaMinusDecay(aParentNucleus, brs[BetaMinus],
514 0.*MeV, 0.*MeV,
516 summedDecayTable->Insert(theBetaMinusChannel);
517 break;
518
519 case BetaPlus:
520 // Decay mode is beta+ + EC.
521 theBetaPlusChannel = new G4BetaPlusDecay(aParentNucleus, brs[BetaPlus],
522 0.*MeV, 0.*MeV,
524 summedDecayTable->Insert(theBetaPlusChannel);
525 break;
526
527 case Alpha:
528 // Decay mode is alpha.
529 theAlphaChannel = new G4AlphaDecay(aParentNucleus, brs[Alpha], 0.*MeV,
530 0.*MeV, noFloat);
531 summedDecayTable->Insert(theAlphaChannel);
532 break;
533
534 case Proton:
535 // Decay mode is proton.
536 theProtonChannel = new G4ProtonDecay(aParentNucleus, brs[Proton], 0.*MeV,
537 0.*MeV, noFloat);
538 summedDecayTable->Insert(theProtonChannel);
539 break;
540
541 case Neutron:
542 // Decay mode is neutron.
543 theNeutronChannel = new G4NeutronDecay(aParentNucleus, brs[Neutron], 0.*MeV,
544 0.*MeV, noFloat);
545 summedDecayTable->Insert(theNeutronChannel);
546 break;
547
548 case SpFission:
549 // Decay mode is spontaneous fission
550 theFissionChannel = new G4SFDecay(aParentNucleus, brs[SpFission], 0.*MeV,
551 0.*MeV, noFloat);
552 summedDecayTable->Insert(theFissionChannel);
553 break;
554
555 case BDProton:
556 // Not yet implemented
557 break;
558
559 case BDNeutron:
560 // Not yet implemented
561 break;
562
563 case Beta2Minus:
564 // Not yet implemented
565 break;
566
567 case Beta2Plus:
568 // Not yet implemented
569 break;
570
571 case Proton2:
572 // Not yet implemented
573 break;
574
575 case Neutron2:
576 // Not yet implemented
577 break;
578
579 case Triton:
580 // Decay mode is Triton.
581 theTritonChannel = new G4TritonDecay(aParentNucleus, brs[Triton], 0.*MeV,
582 0.*MeV, noFloat);
583 summedDecayTable->Insert(theTritonChannel);
584 break;
585
586 default:
587 break;
588 }
589 }
590 }
591
592 // loop over all branches in summedDecayTable
593 //
594 for (G4int i = 0; i < summedDecayTable->entries(); i++){
595 theChannel = summedDecayTable->GetDecayChannel(i);
596 theNuclearDecayChannel = static_cast<G4NuclearDecay*>(theChannel);
597 theBR = theChannel->GetBR();
598 theDaughterNucleus = theNuclearDecayChannel->GetDaughterNucleus();
599
600 // First check if the decay of the original nucleus is an IT channel,
601 // if true create a new ground-state nucleus
602 if (theNuclearDecayChannel->GetDecayMode() == IT && nGeneration == 1) {
603
604 A = ((const G4Ions*)(theDaughterNucleus))->GetAtomicMass();
605 Z = ((const G4Ions*)(theDaughterNucleus))->GetAtomicNumber();
606 theDaughterNucleus=theIonTable->GetIon(Z,A,0.);
607 }
608 if (IsApplicable(*theDaughterNucleus) && theBR > 0.0 &&
609 aParentNucleus != theDaughterNucleus) {
610 // need to make sure daughter has decay table
611 parentDecayTable = GetDecayTable1(theDaughterNucleus);
612 if (parentDecayTable->entries() ) {
613 A = ((const G4Ions*)(theDaughterNucleus))->GetAtomicMass();
614 Z = ((const G4Ions*)(theDaughterNucleus))->GetAtomicNumber();
615 E = ((const G4Ions*)(theDaughterNucleus))->GetExcitationEnergy();
616
617 TaoPlus = theDaughterNucleus->GetPDGLifeTime();
618 if (TaoPlus <= 0.) TaoPlus = 1e-100;
619
620 // first set the taos, one simply need to add to the parent ones
621 taos.clear();
622 taos = TP; // load lifetimes of all previous generations
623 size_t k;
624 //check that TaoPlus differs from other taos from at least 1.e5 relative difference
625 //for (k = 0; k < TP.size(); k++){
626 //if (std::abs((TaoPlus-TP[k])/TP[k])<1.e-5 ) TaoPlus=1.00001*TP[k];
627 //}
628 taos.push_back(TaoPlus); // add daughter lifetime to list
629 // now calculate the coefficiencies
630 //
631 // they are in two parts, first the less than n ones
632 // Eq 4.24 of the TN
633 Acoeffs.clear();
634 long double ta1,ta2;
635 ta2 = (long double)TaoPlus;
636 for (k = 0; k < RP.size(); k++){
637 ta1 = (long double)TP[k]; // loop over lifetimes of all previous generations
638 if (ta1 == ta2) {
639 theRate = 1.e100;
640 } else {
641 theRate = ta1/(ta1-ta2);
642 }
643 theRate = theRate * theBR * RP[k];
644 Acoeffs.push_back(theRate);
645 }
646
647 // the second part: the n:n coefficiency
648 // Eq 4.25 of the TN. Note Yn+1 is zero apart from Y1 which is -1
649 // as treated at line 1013
650 theRate = 0.;
651 long double aRate, aRate1;
652 aRate1 = 0.L;
653 for (k = 0; k < RP.size(); k++){
654 ta1 = (long double)TP[k];
655 if (ta1 == ta2 ) {
656 aRate = 1.e100;
657 } else {
658 aRate = ta2/(ta1-ta2);
659 }
660 aRate = aRate * (long double)(theBR * RP[k]);
661 aRate1 += aRate;
662 }
663 theRate = -aRate1;
664 Acoeffs.push_back(theRate);
665 SetDecayRate (Z,A,E,nGeneration,Acoeffs,taos);
667 nEntry++;
668 } // there are entries in the table
669 } // nuclide is OK to decay
670 } // end of loop (i) over decay table branches
671
672 delete summedDecayTable;
673
674 } // Getting contents of decay rate vector (end loop on j)
675 nS = nT;
676 nT = nEntry;
677 if (nS == nT) stable = true;
678 } // while nuclide is not stable
679
680 // end of while loop
681 // the calculation completed here
682
683
684 // fill the first part of the decay rate table
685 // which is the name of the original particle (isotope)
686 chainsFromParent.SetIonName(theParentNucleus.GetParticleName());
687
688 // now fill the decay table with the newly completed decay rate vector
690
691 // finally add the decayratetable to the tablevector
693}
694
696// //
697// SetSourceTimeProfile //
698// read in the source time profile function (histogram) //
699// //
701
703{
704 std::ifstream infile ( filename, std::ios::in );
705 if (!infile) {
707 ed << " Could not open file " << filename << G4endl;
708 G4Exception("G4Radioactivation::SetSourceTimeProfile()", "HAD_RDM_001",
709 FatalException, ed);
710 }
711
712 G4double bin, flux;
713 NSourceBin = -1;
714
715 G4int loop = 0;
716 while (infile >> bin >> flux) { /* Loop checking, 01.09.2015, D.Wright */
717 loop++;
718 if (loop > 10000) {
719 G4Exception("G4Radioactivation::SetSourceTimeProfile()", "HAD_RDM_100",
720 JustWarning, "While loop count exceeded");
721 break;
722 }
723
724 NSourceBin++;
725 if (NSourceBin > 99) {
726 G4Exception("G4Radioactivation::SetSourceTimeProfile()", "HAD_RDM_002",
727 FatalException, "Input source time file too big (>100 rows)");
728
729 } else {
730 SBin[NSourceBin] = bin * s; // Convert read-in time to ns
731 SProfile[NSourceBin] = flux; // Dimensionless
732 }
733 }
734
735 AnalogueMC = false;
736 infile.close();
737
738#ifdef G4VERBOSE
739 if (GetVerboseLevel() > 2)
740 G4cout <<" Source Timeprofile Nbin = " << NSourceBin <<G4endl;
741#endif
742}
743
745// //
746// SetDecayBiasProfile //
747// read in the decay bias scheme function (histogram) //
748// //
750
752{
753 std::ifstream infile(filename, std::ios::in);
754 if (!infile) G4Exception("G4Radioactivation::SetDecayBias()", "HAD_RDM_001",
755 FatalException, "Unable to open bias data file" );
756
757 G4double bin, flux;
758 G4int dWindows = 0;
759 G4int i ;
760
762
763 NDecayBin = -1;
764
765 G4int loop = 0;
766 while (infile >> bin >> flux ) { /* Loop checking, 01.09.2015, D.Wright */
767 NDecayBin++;
768 loop++;
769 if (loop > 10000) {
770 G4Exception("G4Radioactivation::SetDecayBias()", "HAD_RDM_100",
771 JustWarning, "While loop count exceeded");
772 break;
773 }
774
775 if (NDecayBin > 99) {
776 G4Exception("G4Radioactivation::SetDecayBias()", "HAD_RDM_002",
777 FatalException, "Input bias file too big (>100 rows)" );
778 } else {
779 DBin[NDecayBin] = bin * s; // Convert read-in time to ns
780 DProfile[NDecayBin] = flux; // Dimensionless
781 if (flux > 0.) {
782 decayWindows[NDecayBin] = dWindows;
783 dWindows++;
785 theRadioactivityTables.push_back(rTable);
786 }
787 }
788 }
789 for ( i = 1; i<= NDecayBin; i++) DProfile[i] += DProfile[i-1]; // Cumulative flux vs i
790 for ( i = 0; i<= NDecayBin; i++) DProfile[i] /= DProfile[NDecayBin];
791 // Normalize so entries increase from 0 to 1
792 // converted to accumulated probabilities
793
794 AnalogueMC = false;
795 infile.close();
796
797#ifdef G4VERBOSE
798 if (GetVerboseLevel() > 2)
799 G4cout <<" Decay Bias Profile Nbin = " << NDecayBin <<G4endl;
800#endif
801}
802
804// //
805// DecayIt //
806// //
808
811{
812 // Initialize G4ParticleChange object, get particle details and decay table
815 const G4DynamicParticle* theParticle = theTrack.GetDynamicParticle();
816 const G4ParticleDefinition* theParticleDef = theParticle->GetDefinition();
817
818 // First check whether RDM applies to the current logical volume
819 if (!isAllVolumesMode) {
820 if (!std::binary_search(ValidVolumes.begin(), ValidVolumes.end(),
821 theTrack.GetVolume()->GetLogicalVolume()->GetName())) {
822#ifdef G4VERBOSE
823 if (GetVerboseLevel()>0) {
824 G4cout <<"G4RadioactiveDecay::DecayIt : "
825 << theTrack.GetVolume()->GetLogicalVolume()->GetName()
826 << " is not selected for the RDM"<< G4endl;
827 G4cout << " There are " << ValidVolumes.size() << " volumes" << G4endl;
828 G4cout << " The Valid volumes are " << G4endl;
829 for (size_t i = 0; i< ValidVolumes.size(); i++)
830 G4cout << ValidVolumes[i] << G4endl;
831 }
832#endif
834
835 // Kill the parent particle.
840 }
841 }
842
843 // Now check if particle is valid for RDM
844 if (!(IsApplicable(*theParticleDef) ) ) {
845 // Particle is not an ion or is outside the nucleuslimits for decay
846#ifdef G4VERBOSE
847 if (GetVerboseLevel() > 1) {
848 G4cout << "G4RadioactiveDecay::DecayIt : "
849 << theParticleDef->GetParticleName()
850 << " is not an ion or is outside (Z,A) limits set for the decay. "
851 << " Set particle change accordingly. "
852 << G4endl;
853 }
854#endif
856
857 // Kill the parent particle
862 }
863
864 G4DecayTable* theDecayTable = GetDecayTable1(theParticleDef);
865
866 if (theDecayTable == 0 || theDecayTable->entries() == 0) {
867 // No data in the decay table. Set particle change parameters
868 // to indicate this.
869#ifdef G4VERBOSE
870 if (GetVerboseLevel() > 1) {
871 G4cout << "G4RadioactiveDecay::DecayIt : "
872 << "decay table not defined for "
873 << theParticleDef->GetParticleName()
874 << ". Set particle change accordingly. "
875 << G4endl;
876 }
877#endif
879
880 // Kill the parent particle.
885
886 } else {
887 // Data found. Try to decay nucleus
888 if (AnalogueMC) {
890
891 } else {
892 // Proceed with decay using variance reduction
893 G4double energyDeposit = 0.0;
894 G4double finalGlobalTime = theTrack.GetGlobalTime();
895 G4double finalLocalTime = theTrack.GetLocalTime();
896 G4int index;
897 G4ThreeVector currentPosition;
898 currentPosition = theTrack.GetPosition();
899
900 G4IonTable* theIonTable;
901 G4ParticleDefinition* parentNucleus;
902
903 // Get decay chains for the given nuclide
904 if (!IsRateTableReady(*theParticleDef)) CalculateChainsFromParent(*theParticleDef);
905 GetChainsFromParent(*theParticleDef);
906
907 // Declare some of the variables required in the implementation
908 G4int PZ;
909 G4int PA;
910 G4double PE;
911 G4String keyName;
912 std::vector<G4double> PT;
913 std::vector<G4double> PR;
914 G4double taotime;
915 long double decayRate;
916
917 size_t i;
918 G4int numberOfSecondaries;
919 G4int totalNumberOfSecondaries = 0;
920 G4double currentTime = 0.;
921 G4int ndecaych;
922 G4DynamicParticle* asecondaryparticle;
923 std::vector<G4DynamicParticle*> secondaryparticles;
924 std::vector<G4double> pw;
925 std::vector<G4double> ptime;
926 pw.clear();
927 ptime.clear();
928
929 // Now apply the nucleus splitting
930 for (G4int n = 0; n < NSplit; n++) {
931 // Get the decay time following the decay probability function
932 // supplied by user
933 G4double theDecayTime = GetDecayTime();
934 G4int nbin = GetDecayTimeBin(theDecayTime);
935
936 // calculate the first part of the weight function
937 G4double weight1 = 1.;
938 if (nbin == 1) {
939 weight1 = 1./DProfile[nbin-1]
940 *(DBin[nbin]-DBin[nbin-1])/NSplit; // width of window in ns
941 } else if (nbin > 1) {
942 // Go from nbin to nbin-2 because flux entries in file alternate between 0 and 1
943 weight1 = 1./(DProfile[nbin]-DProfile[nbin-2])
944 *(DBin[nbin]-DBin[nbin-1])/NSplit;
945 // weight1 = (probability of choosing one of the bins)*(time width of bin)/NSplit
946 }
947 // it should be calculated in seconds
948 weight1 /= s ;
949
950 // loop over all the possible secondaries of the nucleus
951 // the first one is itself.
952 for (i = 0; i < theDecayRateVector.size(); i++) {
953 PZ = theDecayRateVector[i].GetZ();
954 PA = theDecayRateVector[i].GetA();
955 PE = theDecayRateVector[i].GetE();
956 PT = theDecayRateVector[i].GetTaos();
957 PR = theDecayRateVector[i].GetDecayRateC();
958
959 // The array of arrays theDecayRateVector contains all possible decay
960 // chains of a given parent nucleus (ZP,AP,EP) to a given descendant
961 // nuclide (Z,A,E).
962 //
963 // theDecayRateVector[0] contains the decay parameters of the parent
964 // nucleus
965 // PZ = ZP
966 // PA = AP
967 // PE = EP
968 // PT[] = {TP}
969 // PR[] = {RP}
970 //
971 // theDecayRateVector[1] contains the decay of the parent to the first
972 // generation daughter (Z1,A1,E1).
973 // PZ = Z1
974 // PA = A1
975 // PE = E1
976 // PT[] = {TP, T1}
977 // PR[] = {RP, R1}
978 //
979 // theDecayRateVector[2] contains the decay of the parent to the first
980 // generation daughter (Z1,A1,E1) and the decay of the first
981 // generation daughter to the second generation daughter (Z2,A2,E2).
982 // PZ = Z2
983 // PA = A2
984 // PE = E2
985 // PT[] = {TP, T1, T2}
986 // PR[] = {RP, R1, R2}
987 //
988 // theDecayRateVector[3] may contain a branch chain
989 // PZ = Z2a
990 // PA = A2a
991 // PE = E2a
992 // PT[] = {TP, T1, T2a}
993 // PR[] = {RP, R1, R2a}
994 //
995 // and so on.
996
997 // Calculate the decay rate of the isotope. decayRate is the
998 // radioactivity of isotope (PZ,PA,PE) at 'theDecayTime'
999 // it will be used to calculate the statistical weight of the
1000 // decay products of this isotope
1001
1002 // For each nuclide, calculate all the decay chains which can reach
1003 // the parent nuclide
1004 decayRate = 0.L;
1005 for (G4int j = 0; j < G4int(PT.size() ); j++) {
1006 taotime = ConvolveSourceTimeProfile(theDecayTime,PT[j]);
1007 decayRate -= PR[j] * (long double)taotime; // PRs are Acoeffs, taotime is inverse time
1008 // Eq.4.23 of of the TN
1009 // note the negative here is required as the rate in the
1010 // equation is defined to be negative,
1011 // i.e. decay away, but we need positive value here.
1012
1013 // G4cout << j << "\t"<< PT[j]/s << "\t" << PR[j] << "\t" << decayRate << G4endl;
1014 }
1015
1016 // At this point any negative decay rates are probably small enough
1017 // (order 10**-30) that negative values are likely due to cancellation
1018 // errors. Set them to zero.
1019 if (decayRate < 0.0) decayRate = 0.0;
1020
1021 // G4cout <<theDecayTime/s <<"\t"<<nbin<<G4endl;
1022 // G4cout << theTrack.GetWeight() <<"\t"<<weight1<<"\t"<<decayRate<< G4endl;
1023
1024 // Add isotope to the radioactivity tables
1025 // One table for each observation time window specified in
1026 // SetDecayBias(G4String filename)
1027
1029 ->AddIsotope(PZ,PA,PE,weight1*decayRate,theTrack.GetWeight());
1030
1031 // Now calculate the statistical weight
1032 // One needs to fold the source bias function with the decaytime
1033 // also need to include the track weight! (F.Lei, 28/10/10)
1034 G4double weight = weight1*decayRate*theTrack.GetWeight();
1035
1036 // decay the isotope
1038 parentNucleus = theIonTable->GetIon(PZ,PA,PE);
1039
1040 // Create a temprary products buffer.
1041 // Its contents to be transfered to the products at the end of the loop
1042 G4DecayProducts* tempprods = 0;
1043
1044 // Decide whether to apply branching ratio bias or not
1045 if (BRBias) {
1046 G4DecayTable* decayTable = GetDecayTable1(parentNucleus);
1047 ndecaych = G4int(decayTable->entries()*G4UniformRand());
1048 G4VDecayChannel* theDecayChannel = decayTable->GetDecayChannel(ndecaych);
1049
1050 if (theDecayChannel == 0) {
1051 // Decay channel not found.
1052
1053 if (GetVerboseLevel() > 0) {
1054 G4cout << " G4RadioactiveDecay::DoIt : cannot determine decay channel ";
1055 G4cout << " for this nucleus; decay as if no biasing active. ";
1056 G4cout << G4endl;
1057 decayTable ->DumpInfo();
1058 }
1059
1060 tempprods = DoDecay(*parentNucleus); // DHW 6 Dec 2010 - do decay as if no biasing
1061 // to avoid deref of temppprods = 0
1062 } else {
1063 // A decay channel has been identified, so execute the DecayIt.
1064 G4double tempmass = parentNucleus->GetPDGMass();
1065 tempprods = theDecayChannel->DecayIt(tempmass);
1066 weight *= (theDecayChannel->GetBR())*(decayTable->entries());
1067 }
1068 } else {
1069 tempprods = DoDecay(*parentNucleus);
1070 }
1071
1072 // save the secondaries for buffers
1073 numberOfSecondaries = tempprods->entries();
1074 currentTime = finalGlobalTime + theDecayTime;
1075 for (index = 0; index < numberOfSecondaries; index++) {
1076 asecondaryparticle = tempprods->PopProducts();
1077 if (asecondaryparticle->GetDefinition()->GetPDGStable() ) {
1078 pw.push_back(weight);
1079 ptime.push_back(currentTime);
1080 secondaryparticles.push_back(asecondaryparticle);
1081 }
1082 // Generate gammas and Xrays from excited nucleus, added by L.Desorgher
1083 else if (((const G4Ions*)(asecondaryparticle->GetDefinition()))
1084 ->GetExcitationEnergy() > 0. && weight > 0.) { //Compute the gamma
1085 G4ParticleDefinition* apartDef = asecondaryparticle->GetDefinition();
1086 AddDeexcitationSpectrumForBiasMode(apartDef,weight,currentTime,pw,
1087 ptime,secondaryparticles);
1088 }
1089 }
1090
1091 delete tempprods;
1092 } // end of i loop
1093 } // end of n loop
1094
1095 // now deal with the secondaries in the two stl containers
1096 // and submmit them back to the tracking manager
1097 totalNumberOfSecondaries = pw.size();
1098 fParticleChangeForRadDecay.SetNumberOfSecondaries(totalNumberOfSecondaries);
1099 for (index=0; index < totalNumberOfSecondaries; index++) {
1100 G4Track* secondary = new G4Track(secondaryparticles[index],
1101 ptime[index], currentPosition);
1102 secondary->SetGoodForTrackingFlag();
1103 secondary->SetTouchableHandle(theTrack.GetTouchableHandle());
1104 secondary->SetWeight(pw[index]);
1106 }
1107
1108 // Kill the parent particle
1112 // Reset NumberOfInteractionLengthLeft.
1114 } // end VR decay
1115
1117 } // end of data found branch
1118}
1119
1120
1121// Add gamma, X-ray, conversion and auger electrons for bias mode
1122void
1124 G4double weight,G4double currentTime,
1125 std::vector<double>& weights_v,
1126 std::vector<double>& times_v,
1127 std::vector<G4DynamicParticle*>& secondaries_v)
1128{
1129 G4double elevel=((const G4Ions*)(apartDef))->GetExcitationEnergy();
1130 G4double life_time=apartDef->GetPDGLifeTime();
1131 G4ITDecay* anITChannel = 0;
1132
1133 while (life_time < halflifethreshold && elevel > 0.) {
1134 anITChannel = new G4ITDecay(apartDef, 100., elevel, elevel, photonEvaporation);
1135 G4DecayProducts* pevap_products = anITChannel->DecayIt(0.);
1136 G4int nb_pevapSecondaries = pevap_products->entries();
1137
1138 G4DynamicParticle* a_pevap_secondary = 0;
1139 G4ParticleDefinition* secDef = 0;
1140 for (G4int ind = 0; ind < nb_pevapSecondaries; ind++) {
1141 a_pevap_secondary= pevap_products->PopProducts();
1142 secDef = a_pevap_secondary->GetDefinition();
1143
1144 if (secDef->GetBaryonNumber() > 4) {
1145 elevel = ((const G4Ions*)(secDef))->GetExcitationEnergy();
1146 life_time = secDef->GetPDGLifeTime();
1147 apartDef = secDef;
1148 if (secDef->GetPDGStable() ) {
1149 weights_v.push_back(weight);
1150 times_v.push_back(currentTime);
1151 secondaries_v.push_back(a_pevap_secondary);
1152 }
1153 } else {
1154 weights_v.push_back(weight);
1155 times_v.push_back(currentTime);
1156 secondaries_v.push_back(a_pevap_secondary);
1157 }
1158 }
1159
1160 delete anITChannel;
1161 delete pevap_products;
1162 }
1163}
1164
@ allowed
@ JustWarning
@ FatalException
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.cc:35
std::ostringstream G4ExceptionDescription
Definition: G4Exception.hh:40
G4ForceCondition
#define noFloat
Definition: G4Ions.hh:112
G4RadioactiveDecayMode
@ G4RadioactiveDecayModeSize
static constexpr double s
Definition: G4SIunits.hh:154
static constexpr double nanosecond
Definition: G4SIunits.hh:138
static constexpr double MeV
Definition: G4SIunits.hh:200
@ fStopAndKill
double G4double
Definition: G4Types.hh:83
bool G4bool
Definition: G4Types.hh:86
int G4int
Definition: G4Types.hh:85
const G4int Z[17]
const G4double A[17]
#define G4endl
Definition: G4ios.hh:57
G4GLOB_DLL std::ostream G4cout
#define G4UniformRand()
Definition: Randomize.hh:52
G4int entries() const
G4DynamicParticle * PopProducts()
G4VDecayChannel * GetDecayChannel(G4int index) const
void Insert(G4VDecayChannel *aChannel)
Definition: G4DecayTable.cc:53
void DumpInfo() const
G4int entries() const
G4ParticleDefinition * GetDefinition() const
virtual G4DecayProducts * DecayIt(G4double)
Definition: G4ITDecay.cc:73
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
Definition: G4IonTable.cc:522
Definition: G4Ions.hh:52
G4double LifeTime(size_t i) const
size_t NearestLevelIndex(G4double energy, size_t index=0) const
G4double NearestLevelEnergy(G4double energy, size_t index=0) const
size_t NumberOfTransitions() const
const G4String & GetName() const
G4double GetDaughterExcitation()
G4ParticleDefinition * GetDaughterNucleus()
G4RadioactiveDecayMode GetDecayMode()
const G4LevelManager * GetLevelManager(G4int Z, G4int A)
static G4NuclearLevelData * GetInstance()
virtual void Initialize(const G4Track &)
void AddSecondary(G4Track *aSecondary)
G4bool GetPDGStable() const
G4double GetPDGLifeTime() const
const G4String & GetParticleName() const
G4IonTable * GetIonTable() const
static G4ParticleTable * GetParticleTable()
G4RadioactiveDecayRates theDecayRateVector
G4int GetDecayTimeBin(const G4double aDecayTime)
void AddDeexcitationSpectrumForBiasMode(G4ParticleDefinition *apartDef, G4double weight, G4double currenTime, std::vector< double > &weights_v, std::vector< double > &times_v, std::vector< G4DynamicParticle * > &secondaries_v)
void SetDecayBias(G4String filename)
G4bool IsRateTableReady(const G4ParticleDefinition &)
virtual void ProcessDescription(std::ostream &outFile) const
std::vector< G4RadioactivityTable * > theRadioactivityTables
G4Radioactivation(const G4String &processName="Radioactivation")
G4VParticleChange * DecayIt(const G4Track &theTrack, const G4Step &theStep)
G4RadioactiveDecayRatesToDaughter ratesToDaughter
G4double ConvolveSourceTimeProfile(const G4double, const G4double)
G4RadioactiveDecayChainsFromParent chainsFromParent
void CalculateChainsFromParent(const G4ParticleDefinition &)
G4double DProfile[100]
G4RadioactiveDecayParentChainTable theParentChainTable
void SetDecayRate(G4int, G4int, G4double, G4int, std::vector< G4double >, std::vector< G4double >)
G4DecayTable * GetDecayTable1(const G4ParticleDefinition *)
void SetSourceTimeProfile(G4String filename)
G4double SProfile[100]
void GetChainsFromParent(const G4ParticleDefinition &)
G4double GetMeanLifeTime(const G4Track &theTrack, G4ForceCondition *condition)
G4RadioactivationMessenger * theRadioactivationMessenger
void SetItsRates(G4RadioactiveDecayRates arate)
void SetTaos(std::vector< G4double > value)
void SetDecayRateC(std::vector< G4double > value)
void DecayAnalog(const G4Track &theTrack)
std::vector< G4String > ValidVolumes
G4DecayTable * LoadDecayTable(const G4ParticleDefinition &theParentNucleus)
G4double GetMeanLifeTime(const G4Track &theTrack, G4ForceCondition *condition)
G4DecayProducts * DoDecay(const G4ParticleDefinition &theParticleDef)
static const G4double levelTolerance
G4ParticleChangeForRadDecay fParticleChangeForRadDecay
G4int GetVerboseLevel() const
G4bool IsApplicable(const G4ParticleDefinition &)
G4PhotonEvaporation * photonEvaporation
Definition: G4Step.hh:62
G4VPhysicalVolume * GetVolume() const
G4double GetWeight() const
void SetWeight(G4double aValue)
const G4ThreeVector & GetPosition() const
void SetTouchableHandle(const G4TouchableHandle &apValue)
G4double GetGlobalTime() const
G4double GetLocalTime() const
const G4DynamicParticle * GetDynamicParticle() const
const G4TouchableHandle & GetTouchableHandle() const
void SetGoodForTrackingFlag(G4bool value=true)
G4double GetBR() const
virtual G4DecayProducts * DecayIt(G4double parentMass=-1.0)=0
void ProposeTrackStatus(G4TrackStatus status)
void ProposeWeight(G4double finalWeight)
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
void SetNumberOfSecondaries(G4int totSecondaries)
G4LogicalVolume * GetLogicalVolume() const
void ClearNumberOfInteractionLengthLeft()
Definition: G4VProcess.hh:424
Definition: DoubConv.h:17
static const G4double AP[5]
Definition: paraMaker.cc:42
#define ns
Definition: xmlparse.cc:614