88 using namespace CLHEP;
95 G4cout <<
"G4Radioactivation constructor: processName = " << processName
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";
167 G4cout <<
"The DecayRate Table for " << aParticleName <<
" is selected."
191 ed <<
" While count exceeded " <<
G4endl;
192 while (t >
SBin[nbin]) {
195 G4Exception(
"G4RadioactiveDecay::ConvolveSourceTimeProfile()",
208 for (
G4int i = 0; i < nbin; i++) {
211 convolvedTime +=
SProfile[i] * std::exp((
SBin[i] - t)/tau) *
215 (std::exp(-(t-
SBin[i+1])/tau)-std::exp(-(t-
SBin[i])/tau));
219 convolvedTime -=
SProfile[nbin] * std::expm1((
SBin[nbin] - t)/tau);
222 if (convolvedTime < 0.) {
223 G4cout <<
" Convolved time =: " << convolvedTime <<
" reset to zero! " <<
G4endl;
230 G4cout <<
" Convolved time: " << convolvedTime <<
G4endl;
232 return convolvedTime;
252 ed <<
" While count exceeded " <<
G4endl;
268 G4cout <<
" Decay time: " <<decaytime/
s <<
"[s]" <<G4endl;
280 ed <<
" While count exceeded " <<
G4endl;
281 while ( aDecayTime >
DBin[i] ) {
296 G4int theG, std::vector<G4double> theCoefficients,
297 std::vector<G4double> theTaos)
322 G4int nGeneration = 0;
324 std::vector<G4double> taos;
327 std::vector<G4double> Acoeffs;
330 Acoeffs.push_back(-1.);
332 G4int A = ((
const G4Ions*)(&theParentNucleus))->GetAtomicMass();
333 G4int Z = ((
const G4Ions*)(&theParentNucleus))->GetAtomicNumber();
334 G4double E = ((
const G4Ions*)(&theParentNucleus))->GetExcitationEnergy();
336 if (tao < 0.) tao = 1
e-100;
368 std::vector<G4double> TP;
369 std::vector<G4double> RP;
373 G4int nearestLevelIndex = 0;
381 const G4int nMode = 9;
389 ed <<
" While count exceeded " <<
G4endl;
398 for (j = nS; j < nT; j++) {
406 G4cout <<
"G4RadioactiveDecay::CalculateChainsFromParent: daughters of ("
407 << ZP <<
", " << AP <<
", " << EP
408 <<
") are being calculated, generation = " << nGeneration
415 aParentNucleus = theIonTable->
GetIon(ZP,AP,EP);
428 for (
G4int k = 0; k < nMode; k++) brs[k] = 0.0;
431 for (i = 0; i < parentDecayTable->
entries(); i++) {
433 theNuclearDecayChannel =
static_cast<G4NuclearDecay*
>(theChannel);
438 AD = ((
const G4Ions*)(theDaughterNucleus))->GetAtomicMass();
439 ZD = ((
const G4Ions*)(theDaughterNucleus))->GetAtomicNumber();
445 if (std::abs(daughterExcitation - nearestEnergy) <
levelTolerance) {
451 summedDecayTable->
Insert(theChannel);
453 brs[theDecayMode] += theChannel->
GetBR();
456 brs[theDecayMode] += theChannel->
GetBR();
459 brs[theDecayMode] += theChannel->
GetBR();
463 brs[2] = brs[2]+brs[3]+brs[4]+brs[5];
464 brs[3] = brs[4] =brs[5] = 0.0;
465 for (i= 0; i<nMode; i++){
470 theITChannel =
new G4ITDecay(aParentNucleus, brs[0], 0.0, 0.0,
473 summedDecayTable->
Insert(theITChannel);
481 summedDecayTable->
Insert(theBetaMinusChannel);
489 summedDecayTable->
Insert(theBetaPlusChannel);
496 summedDecayTable->
Insert(theAlphaChannel);
503 summedDecayTable->
Insert(theProtonChannel);
509 summedDecayTable->
Insert(theNeutronChannel);
519 for (i = 0; i < summedDecayTable->
entries(); i++){
521 theNuclearDecayChannel =
static_cast<G4NuclearDecay*
>(theChannel);
522 theBR = theChannel->
GetBR();
527 if (theNuclearDecayChannel->
GetDecayMode() ==
IT && nGeneration == 1) {
529 A = ((
const G4Ions*)(theDaughterNucleus))->GetAtomicMass();
530 Z = ((
const G4Ions*)(theDaughterNucleus))->GetAtomicNumber();
531 theDaughterNucleus=theIonTable->
GetIon(Z,A,0.);
534 aParentNucleus != theDaughterNucleus) {
538 if (parentDecayTable->
entries() ) {
539 A = ((
const G4Ions*)(theDaughterNucleus))->GetAtomicMass();
540 Z = ((
const G4Ions*)(theDaughterNucleus))->GetAtomicNumber();
541 E = ((
const G4Ions*)(theDaughterNucleus))->GetExcitationEnergy();
544 if (TaoPlus <= 0.) TaoPlus = 1
e-100;
554 taos.push_back(TaoPlus);
561 ta2 = (
long double)TaoPlus;
562 for (k = 0; k < RP.size(); k++){
563 ta1 = (
long double)TP[k];
567 theRate = ta1/(ta1-ta2);
569 theRate = theRate * theBR * RP[k];
570 Acoeffs.push_back(theRate);
577 long double aRate, aRate1;
579 for (k = 0; k < RP.size(); k++){
580 ta1 = (
long double)TP[k];
584 aRate = ta2/(ta1-ta2);
586 aRate = aRate * (
long double)(theBR * RP[k]);
590 Acoeffs.push_back(theRate);
602 if (nS == nT) stable =
true;
632 ed <<
" Could not open file " << filename <<
G4endl;
633 G4Exception(
"G4RadioactiveDecay::SetSourceTimeProfile()",
"HAD_RDM_001",
642 ed <<
" While count exceeded " <<
G4endl;
644 while (infile >> bin >> flux ) {
653 G4Exception(
"G4RadioactiveDecay::SetSourceTimeProfile()",
"HAD_RDM_002",
681 if (!infile)
G4Exception(
"G4RadioactiveDecay::SetDecayBias()",
"HAD_RDM_003",
694 ed <<
" While count exceeded " <<
G4endl;
696 while (infile >> bin >> flux ) {
705 G4Exception(
"G4RadioactiveDecay::SetDecayBias()",
"HAD_RDM_004",
752 G4cout <<
"G4RadioactiveDecay::DecayIt : "
754 <<
" is not selected for the RDM"<<
G4endl;
777 G4cerr <<
"G4RadioactiveDecay::DecayIt : "
779 <<
" is not a valid nucleus for the RDM"<<
G4endl;
792 if (theDecayTable == 0 || theDecayTable->
entries() == 0) {
797 G4cerr <<
"G4RadioactiveDecay::DecayIt : decay table not defined for ";
829 std::vector<G4double> PT;
830 std::vector<G4double> PR;
832 long double decayRate;
836 G4int numberOfSecondaries;
837 G4int totalNumberOfSecondaries = 0;
841 std::vector<G4DynamicParticle*> secondaryparticles;
842 std::vector<G4double> pw;
843 std::vector<G4double> ptime;
859 }
else if (nbin > 1) {
924 for (j = 0; j < PT.size(); j++) {
928 decayRate -= PR[j] * (
long double)taotime;
940 if (decayRate < 0.0) decayRate = 0.0;
965 ->AddIsotope(PZ,PA,PE,weight1*decayRate,theTrack.
GetWeight());
974 parentNucleus = theIonTable->
GetIon(PZ,PA,PE);
986 if (theDecayChannel == 0) {
990 G4cerr <<
" G4RadioactiveDecay::DoIt : cannot determine decay channel ";
991 G4cerr <<
" for this nucleus; decay as if no biasing active ";
996 tempprods =
DoDecay(*parentNucleus);
1001 tempprods = theDecayChannel->DecayIt(tempmass);
1002 weight *= (theDecayChannel->GetBR())*(decayTable->
entries());
1005 tempprods =
DoDecay(*parentNucleus);
1009 numberOfSecondaries = tempprods->
entries();
1010 currentTime = finalGlobalTime + theDecayTime;
1011 for (index = 0; index < numberOfSecondaries; index++) {
1014 pw.push_back(weight);
1015 ptime.push_back(currentTime);
1016 secondaryparticles.push_back(asecondaryparticle);
1019 else if (((
const G4Ions*)(asecondaryparticle->
GetDefinition()))->GetExcitationEnergy()>0. && weight>0.){
1031 totalNumberOfSecondaries = pw.size();
1033 for (index=0; index < totalNumberOfSecondaries; index++) {
1035 ptime[index], currentPosition);
1058 std::vector<double>& weights_v,
1059 std::vector<double>& times_v,
1060 std::vector<G4DynamicParticle*>& secondaries_v)
1062 G4double elevel=((
const G4Ions*)(apartDef))->GetExcitationEnergy();
1066 while (life_time <halflifethreshold && elevel>0.) {
1073 for (
G4int ind = 0; ind < nb_pevapSecondaries; ind++) {
1078 elevel = ((
const G4Ions*)(secDef))->GetExcitationEnergy();
1082 weights_v.push_back(weight);
1083 times_v.push_back(currentTime);
1084 secondaries_v.push_back(a_pevap_secondary);
1087 weights_v.push_back(weight);
1088 times_v.push_back(currentTime);
1089 secondaries_v.push_back(a_pevap_secondary);
const G4LevelManager * GetLevelManager(G4int Z, G4int A)
G4double GetLocalTime() const
G4int GetDecayTimeBin(const G4double aDecayTime)
G4DecayTable * GetDecayTable(const G4ParticleDefinition *)
std::ostringstream G4ExceptionDescription
G4LogicalVolume * GetLogicalVolume() const
static G4ParticleTable * GetParticleTable()
static constexpr double MeV
void CalculateChainsFromParent(const G4ParticleDefinition &)
G4double GetPDGLifeTime() const
void SetWeight(G4double aValue)
void GetChainsFromParent(const G4ParticleDefinition &)
virtual G4DecayProducts * DecayIt(G4double)
void ProposeWeight(G4double finalWeight)
void AddSecondary(G4Track *aSecondary)
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)
void SetItsRates(G4RadioactiveDecayRates arate)
G4DynamicParticle * PopProducts()
G4double GetPDGMass() const
G4double GetWeight() const
G4RadioactiveDecayParentChainTable theParentChainTable
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
G4int GetBaryonNumber() const
void SetTaos(std::vector< G4double > value)
G4Radioactivation(const G4String &processName="Radioactivation")
G4VDecayChannel * GetDecayChannel(G4int index) const
std::vector< G4String > ValidVolumes
G4DecayProducts * DoDecay(const G4ParticleDefinition &theParticleDef)
void SetE(G4double value)
G4double GetGlobalTime() const
double A(double temperature)
void SetTouchableHandle(const G4TouchableHandle &apValue)
G4double ConvolveSourceTimeProfile(const G4double, const G4double)
const G4ThreeVector & GetPosition() const
G4double LifeTime(size_t i) const
G4ParticleChangeForRadDecay fParticleChangeForRadDecay
void ProposeLocalTime(G4double t)
G4RadioactiveDecayChainsFromParent chainsFromParent
G4int GetVerboseLevel() const
static constexpr double nanosecond
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)
G4RadioactiveDecayRates theDecayRateVector
std::vector< G4RadioactivityTable * > theRadioactivityTables
void SetDecayBias(G4String filename)
size_t NumberOfTransitions() const
G4bool IsRateTableReady(const G4ParticleDefinition &)
void SetProcessSubType(G4int)
G4double GetDaughterExcitation()
G4GLOB_DLL std::ostream G4cout
void ClearNumberOfInteractionLengthLeft()
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
void AddDeexcitationSpectrumForBiasMode(G4ParticleDefinition *apartDef, G4double weight, G4double currenTime, std::vector< double > &weights_v, std::vector< double > ×_v, std::vector< G4DynamicParticle * > &secondaries_v)
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
void SetGeneration(G4int value)
G4RadioactiveDecayMode GetDecayMode()
virtual void Initialize(const G4Track &)
void ProposeTrackStatus(G4TrackStatus status)
const G4DynamicParticle * GetDynamicParticle() const
G4RadioactivationMessenger * theRadioactivationMessenger
G4RadioactiveDecayRatesToDaughter ratesToDaughter
G4ParticleDefinition * GetDaughterNucleus()
size_t NearestLevelIndex(G4double energy, size_t index=0) const
G4double halflifethreshold
void SetIonName(G4String name)
const G4String & GetName() const
G4VPhysicalVolume * GetVolume() const
void SetNumberOfSecondaries(G4int totSecondaries)
virtual void ProcessDescription(std::ostream &outFile) const