91 using namespace CLHEP;
98 G4cout <<
"G4Radioactivation constructor: processName = " << processName
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";
147 DecayTableMap::iterator table_ptr =
dkmap->find(key);
150 if (table_ptr ==
dkmap->end() ) {
152 if(theDecayTable) (*dkmap)[key] = theDecayTable;
154 theDecayTable = table_ptr->second;
156 return theDecayTable;
185 G4cout <<
"The DecayRate Table for " << aParticleName <<
" is selected."
208 while (t >
SBin[nbin]) {
211 G4Exception(
"G4Radioactivation::ConvolveSourceTimeProfile()",
212 "HAD_RDM_100",
JustWarning,
"While loop count exceeded");
224 for (
G4int i = 0; i < nbin; i++) {
227 convolvedTime +=
SProfile[i] * std::exp((
SBin[i] - t)/tau) *
231 (std::exp(-(t-
SBin[i+1])/tau)-std::exp(-(t-
SBin[i])/tau));
235 convolvedTime -=
SProfile[nbin] * std::expm1((
SBin[nbin] - t)/tau);
238 if (convolvedTime < 0.) {
239 G4cout <<
" Convolved time =: " << convolvedTime <<
" reset to zero! " <<
G4endl;
246 G4cout <<
" Convolved time: " << convolvedTime <<
G4endl;
248 return convolvedTime;
273 G4Exception(
"G4Radioactivation::GetDecayTime()",
"HAD_RDM_100",
294 while (aDecayTime >
DBin[i] ) {
298 G4Exception(
"G4Radioactivation::GetDecayTimeBin()",
"HAD_RDM_100",
327 G4int theG, std::vector<G4double> theCoefficients,
328 std::vector<G4double> theTaos)
353 G4int nGeneration = 0;
355 std::vector<G4double> taos;
358 std::vector<G4double> Acoeffs;
361 Acoeffs.push_back(-1.);
363 G4int A = ((
const G4Ions*)(&theParentNucleus))->GetAtomicMass();
364 G4int Z = ((
const G4Ions*)(&theParentNucleus))->GetAtomicNumber();
365 G4double E = ((
const G4Ions*)(&theParentNucleus))->GetExcitationEnergy();
367 if (tao < 0.) tao = 1
e-100;
402 std::vector<G4double> TP;
403 std::vector<G4double> RP;
407 G4int nearestLevelIndex = 0;
415 const G4int nMode = 12;
425 G4Exception(
"G4Radioactivation::CalculateChainsFromParent()",
"HAD_RDM_100",
430 for (j = nS; j < nT; j++) {
438 G4cout <<
"G4RadioactiveDecay::CalculateChainsFromParent: daughters of ("
439 << ZP <<
", " << AP <<
", " << EP
440 <<
") are being calculated, generation = " << nGeneration
447 aParentNucleus = theIonTable->
GetIon(ZP,AP,EP);
460 for (
G4int k = 0;
k < nMode;
k++) brs[
k] = 0.0;
463 for (
G4int i = 0; i < parentDecayTable->
entries(); i++) {
465 theNuclearDecayChannel =
static_cast<G4NuclearDecay*
>(theChannel);
469 AD = ((
const G4Ions*)(theDaughterNucleus))->GetAtomicMass();
470 ZD = ((
const G4Ions*)(theDaughterNucleus))->GetAtomicNumber();
485 summedDecayTable->
Insert(theChannel);
487 brs[theDecayMode] += theChannel->
GetBR();
490 brs[theDecayMode] += theChannel->
GetBR();
493 brs[theDecayMode] += theChannel->
GetBR();
498 brs[2] = brs[2]+brs[3]+brs[4]+brs[5]+brs[6];
499 brs[3] = brs[4] = brs[5] = brs[6] = 0.0;
500 for (
G4int i = 0; i < nMode; i++) {
505 theITChannel =
new G4ITDecay(aParentNucleus, brs[0], 0.0, 0.0,
508 summedDecayTable->
Insert(theITChannel);
516 summedDecayTable->
Insert(theBetaMinusChannel);
524 summedDecayTable->
Insert(theBetaPlusChannel);
531 summedDecayTable->
Insert(theAlphaChannel);
538 summedDecayTable->
Insert(theProtonChannel);
545 summedDecayTable->
Insert(theNeutronChannel);
550 theFissionChannel =
new G4SFDecay(aParentNucleus, brs[10], 0.*
MeV,
552 summedDecayTable->
Insert(theFissionChannel);
558 summedDecayTable->
Insert(theTritonChannel);
569 for (
G4int i = 0; i < summedDecayTable->
entries(); i++){
571 theNuclearDecayChannel =
static_cast<G4NuclearDecay*
>(theChannel);
572 theBR = theChannel->
GetBR();
577 if (theNuclearDecayChannel->
GetDecayMode() ==
IT && nGeneration == 1) {
579 A = ((
const G4Ions*)(theDaughterNucleus))->GetAtomicMass();
580 Z = ((
const G4Ions*)(theDaughterNucleus))->GetAtomicNumber();
581 theDaughterNucleus=theIonTable->
GetIon(Z,A,0.);
584 aParentNucleus != theDaughterNucleus) {
587 if (parentDecayTable->
entries() ) {
588 A = ((
const G4Ions*)(theDaughterNucleus))->GetAtomicMass();
589 Z = ((
const G4Ions*)(theDaughterNucleus))->GetAtomicNumber();
590 E = ((
const G4Ions*)(theDaughterNucleus))->GetExcitationEnergy();
593 if (TaoPlus <= 0.) TaoPlus = 1
e-100;
603 taos.push_back(TaoPlus);
610 ta2 = (
long double)TaoPlus;
611 for (k = 0; k < RP.size(); k++){
612 ta1 = (
long double)TP[k];
616 theRate = ta1/(ta1-ta2);
618 theRate = theRate * theBR * RP[
k];
619 Acoeffs.push_back(theRate);
626 long double aRate, aRate1;
628 for (k = 0; k < RP.size(); k++){
629 ta1 = (
long double)TP[k];
633 aRate = ta2/(ta1-ta2);
635 aRate = aRate * (
long double)(theBR * RP[k]);
639 Acoeffs.push_back(theRate);
647 delete summedDecayTable;
652 if (nS == nT) stable =
true;
682 ed <<
" Could not open file " << filename <<
G4endl;
683 G4Exception(
"G4Radioactivation::SetSourceTimeProfile()",
"HAD_RDM_001",
691 while (infile >> bin >> flux) {
694 G4Exception(
"G4Radioactivation::SetSourceTimeProfile()",
"HAD_RDM_100",
701 G4Exception(
"G4Radioactivation::SetSourceTimeProfile()",
"HAD_RDM_002",
729 if (!infile)
G4Exception(
"G4Radioactivation::SetDecayBias()",
"HAD_RDM_001",
741 while (infile >> bin >> flux ) {
745 G4Exception(
"G4Radioactivation::SetDecayBias()",
"HAD_RDM_100",
751 G4Exception(
"G4Radioactivation::SetDecayBias()",
"HAD_RDM_002",
799 G4cout <<
"G4RadioactiveDecay::DecayIt : "
801 <<
" is not selected for the RDM"<<
G4endl;
823 G4cout <<
"G4RadioactiveDecay::DecayIt : "
825 <<
" is not an ion or is outside (Z,A) limits set for the decay. "
826 <<
" Set particle change accordingly. "
841 if (theDecayTable == 0 || theDecayTable->
entries() == 0) {
846 G4cout <<
"G4RadioactiveDecay::DecayIt : "
847 <<
"decay table not defined for "
849 <<
". Set particle change accordingly. "
887 std::vector<G4double> PT;
888 std::vector<G4double> PR;
890 long double decayRate;
893 G4int numberOfSecondaries;
894 G4int totalNumberOfSecondaries = 0;
898 std::vector<G4DynamicParticle*> secondaryparticles;
899 std::vector<G4double> pw;
900 std::vector<G4double> ptime;
916 }
else if (nbin > 1) {
980 for (
G4int j = 0; j <
G4int(PT.size() ); j++) {
982 decayRate -= PR[j] * (
long double)taotime;
994 if (decayRate < 0.0) decayRate = 0.0;
1004 ->AddIsotope(PZ,PA,PE,weight1*decayRate,theTrack.
GetWeight());
1013 parentNucleus = theIonTable->
GetIon(PZ,PA,PE);
1025 if (theDecayChannel == 0) {
1029 G4cout <<
" G4RadioactiveDecay::DoIt : cannot determine decay channel ";
1030 G4cout <<
" for this nucleus; decay as if no biasing active. ";
1035 tempprods =
DoDecay(*parentNucleus);
1040 tempprods = theDecayChannel->DecayIt(tempmass);
1041 weight *= (theDecayChannel->GetBR())*(decayTable->
entries());
1044 tempprods =
DoDecay(*parentNucleus);
1048 numberOfSecondaries = tempprods->
entries();
1049 currentTime = finalGlobalTime + theDecayTime;
1050 for (index = 0; index < numberOfSecondaries; index++) {
1053 pw.push_back(weight);
1054 ptime.push_back(currentTime);
1055 secondaryparticles.push_back(asecondaryparticle);
1059 ->GetExcitationEnergy() > 0. && weight > 0.) {
1062 ptime,secondaryparticles);
1072 totalNumberOfSecondaries = pw.size();
1074 for (index=0; index < totalNumberOfSecondaries; index++) {
1076 ptime[index], currentPosition);
1100 std::vector<double>& weights_v,
1101 std::vector<double>& times_v,
1102 std::vector<G4DynamicParticle*>& secondaries_v)
1104 G4double elevel=((
const G4Ions*)(apartDef))->GetExcitationEnergy();
1108 while (life_time < halflifethreshold && elevel > 0.) {
1115 for (
G4int ind = 0; ind < nb_pevapSecondaries; ind++) {
1120 elevel = ((
const G4Ions*)(secDef))->GetExcitationEnergy();
1124 weights_v.push_back(weight);
1125 times_v.push_back(currentTime);
1126 secondaries_v.push_back(a_pevap_secondary);
1129 weights_v.push_back(weight);
1130 times_v.push_back(currentTime);
1131 secondaries_v.push_back(a_pevap_secondary);
1136 delete pevap_products;