138 using namespace CLHEP;
142 #ifdef G4MULTITHREADED
148 forceDecayDirection(0.,0.,0.), forceDecayHalfAngle(0.*
deg), verboseLevel(0)
152 G4cout <<
"G4RadioactiveDecay constructor: processName = " << processName
177 #ifdef G4MULTITHREADED
214 for (DecayTableMap::iterator i =
dkmap->begin(); i !=
dkmap->end(); i++) {
224 if (((
const G4Ions*)(&aParticle))->GetExcitationEnergy() > 0.) {
return true;}
233 G4int A = ((
const G4Ions*) (&aParticle))->GetAtomicMass();
234 G4int Z = ((
const G4Ions*) (&aParticle))->GetAtomicNumber();
245 DecayTableMap::iterator table_ptr =
dkmap->find(key);
248 if (table_ptr ==
dkmap->end() ) {
250 (*dkmap)[key] = theDecayTable;
252 theDecayTable = table_ptr->second;
255 return theDecayTable;
264 for (
size_t i = 0; i < theLogicalVolumes->size(); i++) {
265 volume=(*theLogicalVolumes)[i];
266 if (volume->
GetName() == aVolume) {
274 }
else if(i == theLogicalVolumes->size()) {
275 G4cerr <<
"SelectAVolume: "<< aVolume
276 <<
" is not a valid logical volume name" <<
G4endl;
287 for (
size_t i = 0; i < theLogicalVolumes->size(); i++){
288 volume=(*theLogicalVolumes)[i];
289 if (volume->
GetName() == aVolume) {
290 std::vector<G4String>::iterator location;
297 G4cerr <<
" DeselectVolume:" << aVolume <<
" is not in the list "
302 G4cout <<
" DeselectVolume: " << aVolume <<
" is removed from list "
305 }
else if (i == theLogicalVolumes->size()) {
306 G4cerr <<
" DeselectVolume:" << aVolume
307 <<
"is not a valid logical volume name" <<
G4endl;
323 for (
size_t i = 0; i < theLogicalVolumes->size(); i++){
324 volume = (*theLogicalVolumes)[i];
374 G4cout <<
"The DecayRate Table for " << aParticleName <<
" is selected."
384 long double taotime = 0.L;
390 while (t >
SBin[nbin]) nbin++;
393 long double ltao = tao;
396 for (
G4int i = 0; i < nbin; i++) {
397 taotime += (
long double)
SProfile[i] *
398 (std::exp(-(lt-(
long double)
SBin[i+1])/ltao)-std::exp(-(lt-(
long double)
SBin[i])/ltao));
401 taotime += (
long double)
SProfile[nbin] * (1.
L-std::exp(-(lt-(
long double)
SBin[nbin])/ltao));
403 G4cout <<
" Tao time =: " <<taotime <<
" reset to zero!"<<
G4endl;
552 while ( aDecayTime >
DBin[i] ) i++;
577 G4cout <<
"G4RadioactiveDecay::GetMeanLifeTime() " <<
G4endl;
579 <<
" GeV, Mass: " << theParticle->
GetMass()/
GeV
580 <<
" GeV, Life time: " << theLife/
ns <<
" ns " <<
G4endl;
584 else if (theLife < 0.0) {meanlife =
DBL_MAX;}
585 else {meanlife = theLife;}
588 if (((
const G4Ions*)(theParticleDef))->GetExcitationEnergy() > 0. &&
589 meanlife ==
DBL_MAX) {meanlife = 0.;}
593 G4cout <<
" mean life time: " << meanlife/
s <<
" s " <<
G4endl;
615 G4cout <<
"G4RadioactiveDecay::GetMeanFreePath() " <<
G4endl;
617 <<
" GeV, Mass: " << aMass/
GeV <<
" GeV, tau: " << tau <<
" ns "
627 ed <<
"Ion has negative lifetime " << tau
628 <<
" but is not stable. Setting mean free path to DBL_MAX" <<
G4endl;
629 G4Exception(
"G4RadioactiveDecay::GetMeanFreePath()",
"HAD_RDM_011",
636 pathlength = c_light*tau*betaGamma;
642 G4cout <<
"G4Decay::GetMeanFreePath: "
644 <<
" stops, kinetic energy = "
654 G4cout <<
"mean free path: "<< pathlength/
m <<
" m" <<
G4endl;
683 #ifdef G4MULTITHREADED
693 G4int A = ((
const G4Ions*)(&theParentNucleus))->GetAtomicMass();
694 G4int Z = ((
const G4Ions*)(&theParentNucleus))->GetAtomicNumber();
695 G4int lvl = ((
const G4Ions*)(&theParentNucleus))->GetIsomerLevel();
698 #ifdef G4MULTITHREADED
699 G4AutoLock lk(&G4RadioactiveDecay::radioactiveDecayMutex);
702 DecayTableMap::iterator master_table_ptr = master_dkmap->find(key);
704 if (master_table_ptr != master_dkmap->end() ) {
705 return master_table_ptr->second;
716 if (!getenv(
"G4RADIOACTIVEDATA") ) {
717 G4cout <<
"Please setenv G4RADIOACTIVEDATA to point to the radioactive decay data files."
719 throw G4HadronicException(__FILE__, __LINE__,
" Please setenv G4RADIOACTIVEDATA to point to the radioactive decay data files.");
721 G4String dirName = getenv(
"G4RADIOACTIVEDATA");
723 std::ostringstream os;
724 os <<dirName <<
"/z" <<Z <<
".a" <<A <<
'\0';
728 std::ifstream DecaySchemeFile(file);
731 if (DecaySchemeFile) {
734 G4bool modeFirstRecord[7];
737 for (
G4int i = 0; i < nMode; i++) {
738 modeFirstRecord[i] =
true;
743 char inputChars[100]={
' '};
750 G4int levelCounter = 0;
757 while (!complete && !DecaySchemeFile.getline(inputChars, 100).eof()) {
758 inputLine = inputChars;
759 inputLine = inputLine.
strip(1);
760 if (inputChars[0] !=
'#' && inputLine.length() != 0) {
761 std::istringstream tmpStream(inputLine);
763 if (inputChars[0] ==
'P') {
766 tmpStream >> recordType >> a >> b;
771 found = (levelCounter == lvl);
778 if (inputChars[0] ==
'W') {
782 G4cout <<
" Warning in G4RadioactiveDecay::LoadDecayTable " <<
G4endl;
788 tmpStream >> theDecayMode >> a >> b >> c >> betaType;
792 if (inputLine.length() < 80) betaType =
allowed;
796 switch (theDecayMode) {
802 (
const G4Ions*)& theParentNucleus, b);
806 theDecayTable->
Insert(anITChannel);
812 if (modeFirstRecord[1]) {
813 modeFirstRecord[1] =
false;
827 for (
G4int ptn = 0; ptn < npti; ptn++) {
829 e = 1. + e0*(ptn+0.5)/100.;
830 p = std::sqrt(e*e - 1.);
831 f = p*e*(e0-e+1)*(e0-e+1);
837 f *= corrections.
ShapeFactor(betaType, p, e0-e+1.);
844 b, c*
MeV, a*MeV, 0,
FBeta, aRandomEnergy);
848 theDecayTable->
Insert(aBetaMinusChannel);
858 if (modeFirstRecord[2]) {
859 modeFirstRecord[2] =
false;
862 e0 = c*
MeV/0.510999 - 2.;
875 for (
G4int ptn = 0; ptn < npti; ptn++) {
877 e = 1. + e0*(ptn+0.5)/100.;
878 p = std::sqrt(e*e - 1.);
879 f = p*e*(e0-e+1)*(e0-e+1);
885 f *= corrections.
ShapeFactor(betaType, p, e0-e+1.);
891 &theParentNucleus, b,
892 (c-1.021998)*
MeV, a*
MeV, 0,
893 FBeta, aRandomEnergy);
897 theDecayTable->
Insert(aBetaPlusChannel);
907 if (modeFirstRecord[3]) {
908 modeFirstRecord[3] =
false;
918 theDecayTable->
Insert(aKECChannel);
925 if (modeFirstRecord[4]) {
926 modeFirstRecord[4] =
false;
935 theDecayTable->
Insert(aLECChannel);
942 if (modeFirstRecord[5]) {
943 modeFirstRecord[5] =
false;
953 theDecayTable->
Insert(aMECChannel);
959 if (modeFirstRecord[6]) {
960 modeFirstRecord[6] =
false;
970 theDecayTable->
Insert(anAlphaChannel);
981 G4Exception(
"G4RadioactiveDecay::LoadDecayTable()",
"HAD_RDM_000",
1002 if (theDecayMode !=
IT) {
1003 theBR = theChannel->
GetBR();
1004 theChannel->
SetBR(theBR*modeTotalBR[theDecayMode]/modeSumBR[theDecayMode]);
1008 DecaySchemeFile.close();
1010 if (!found && lvl > 0) {
1018 theDecayTable->
Insert(anITChannel);
1020 if (!theDecayTable) {
1024 G4cerr <<
"G4RadoactiveDecay::LoadDecayTable() : cannot find ion radioactive decay file "
1027 return theDecayTable;
1031 G4cout <<
"G4RadioactiveDecay::LoadDecayTable()" <<
G4endl;
1035 #ifdef G4MULTITHREADED
1036 (*master_dkmap)[key] = theDecayTable;
1038 return theDecayTable;
1044 if (Z < 1 || A < 2)
G4cout <<
"Z and A not valid!" <<
G4endl;
1046 std::ifstream DecaySchemeFile(filename);
1047 if (DecaySchemeFile) {
1048 G4int ID_ion = A*1000 + Z;
1052 G4cout <<
"The file " << filename <<
" does not exist!" <<
G4endl;
1059 G4int theG, std::vector<G4double> theRates,
1060 std::vector<G4double> theTaos)
1086 G4int nGeneration = 0;
1087 std::vector<G4double> rates;
1088 std::vector<G4double> taos;
1092 rates.push_back(-1.);
1095 G4int A = ((
const G4Ions*)(&theParentNucleus))->GetAtomicMass();
1096 G4int Z = ((
const G4Ions*)(&theParentNucleus))->GetAtomicNumber();
1097 G4double E = ((
const G4Ions*)(&theParentNucleus))->GetExcitationEnergy();
1099 if (tao < 0.) tao = 1e-100;
1100 taos.push_back(tao);
1128 std::vector<G4double> TP;
1129 std::vector<G4double> RP;
1146 for (j = nS; j < nT; j++) {
1153 G4cout <<
"G4RadioactiveDecay::AddDecayRateTable : daughters of ("
1154 << ZP <<
", " << AP <<
", " << EP
1155 <<
") are being calculated, generation = " << nGeneration
1159 aParentNucleus = theIonTable->
GetIon(ZP,AP,EP);
1163 for (
G4int k = 0; k < 7; k++) brs[k] = 0.0;
1166 for (i = 0; i < aTempDecayTable->
entries(); i++) {
1172 AD = ((
const G4Ions*)(theDaughterNucleus))->GetAtomicMass();
1173 ZD = ((
const G4Ions*)(theDaughterNucleus))->GetAtomicNumber();
1185 theDecayTable->
Insert(theChannel);
1187 brs[theDecayMode] += theChannel->
GetBR();
1190 brs[theDecayMode] += theChannel->
GetBR();
1193 brs[theDecayMode] += theChannel->
GetBR();
1196 brs[2] = brs[2]+brs[3]+brs[4]+brs[5];
1197 brs[3] = brs[4] =brs[5] = 0.0;
1198 for (i= 0; i<7; i++){
1204 (
const G4Ions*) aParentNucleus, brs[0]);
1205 theDecayTable->
Insert(theITChannel);
1211 brs[1], 0.*
MeV, 0.*
MeV, 1,
false, 0);
1212 theDecayTable->
Insert(theBetaMinusChannel);
1218 aParentNucleus, brs[2], 0.*
MeV, 0.*
MeV, 1,
false, 0);
1219 theDecayTable->
Insert(theBetaPlusChannel);
1227 theDecayTable->
Insert(theAlphaChannel);
1238 for (i = 0; i < theDecayTable->
entries(); i++){
1241 theBR = theChannel->
GetBR();
1245 if (theNuclearDecayChannel->
GetDecayMode() ==
IT && nGeneration == 1) {
1246 A = ((
const G4Ions*)(theDaughterNucleus))->GetAtomicMass();
1247 Z = ((
const G4Ions*)(theDaughterNucleus))->GetAtomicNumber();
1248 theDaughterNucleus=theIonTable->
GetIon(Z,A,0.);
1251 aParentNucleus != theDaughterNucleus) {
1255 if (aTempDecayTable->
entries() ) {
1256 A = ((
const G4Ions*)(theDaughterNucleus))->GetAtomicMass();
1257 Z = ((
const G4Ions*)(theDaughterNucleus))->GetAtomicNumber();
1258 E = ((
const G4Ions*)(theDaughterNucleus))->GetExcitationEnergy();
1261 if (TaoPlus <= 0.) TaoPlus = 1e-100;
1271 taos.push_back(TaoPlus);
1277 long double ta1,ta2;
1278 ta2 = (
long double)TaoPlus;
1279 for (k = 0; k < RP.size(); k++){
1280 ta1 = (
long double)TP[k];
1284 theRate = ta1/(ta1-ta2);
1286 theRate = theRate * theBR * RP[k];
1287 rates.push_back(theRate);
1294 long double aRate, aRate1;
1296 for (k = 0; k < RP.size(); k++){
1297 ta1 = (
long double)TP[k];
1301 aRate = ta2/(ta1-ta2);
1303 aRate = aRate * (
long double)(theBR * RP[k]);
1307 rates.push_back(theRate);
1319 if (nS == nT) stable =
true;
1346 std::ifstream infile ( filename, std::ios::in );
1349 ed <<
" Could not open file " << filename <<
G4endl;
1350 G4Exception(
"G4RadioactiveDecay::SetSourceTimeProfile()",
"HAD_RDM_001",
1356 while (infile >> bin >> flux ) {
1359 G4Exception(
"G4RadioactiveDecay::SetSourceTimeProfile()",
"HAD_RDM_002",
1386 std::ifstream infile(filename, std::ios::in);
1387 if (!infile)
G4Exception(
"G4RadioactiveDecay::SetDecayBias()",
"HAD_RDM_003",
1397 while (infile >> bin >> flux ) {
1400 G4Exception(
"G4RadioactiveDecay::SetDecayBias()",
"HAD_RDM_004",
1448 G4cout <<
"G4RadioactiveDecay::DecayIt : "
1450 <<
" is not selected for the RDM"<<
G4endl;
1473 G4cerr <<
"G4RadioactiveDecay::DecayIt : "
1475 <<
" is not a valid nucleus for the RDM"<<
G4endl;
1489 if (theDecayTable == 0 || theDecayTable->
entries() == 0) {
1494 G4cerr <<
"G4RadioactiveDecay::DecayIt : decay table not defined for ";
1526 if ( products->
entries() == 1) {
1555 if (temptime < 0.) temptime = 0.;
1556 finalGlobalTime += temptime;
1557 finalLocalTime += temptime;
1560 products->
Boost(ParentEnergy, ParentDirection);
1567 G4cout <<
"G4RadioactiveDecay::DecayIt : Decay vertex :";
1568 G4cout <<
" Time: " <<finalGlobalTime/
ns <<
"[ns]";
1573 G4cout <<
"G4Decay::DecayIt : decay products in Lab. Frame" <<
G4endl;
1578 for (index=0; index < numberOfSecondaries; index++) {
1592 G4cout <<
"DecayIt: Variance Reduction version " <<
G4endl;
1609 std::vector<G4double> PT;
1610 std::vector<G4double> PR;
1612 long double decayRate;
1616 G4int numberOfSecondaries;
1617 G4int totalNumberOfSecondaries = 0;
1621 std::vector<G4DynamicParticle*> secondaryparticles;
1622 std::vector<G4double> pw;
1623 std::vector<G4double> ptime;
1639 }
else if (nbin > 1) {
1664 for (j = 0; j < PT.size(); j++) {
1666 decayRate -= PR[j] * (
long double)taotime;
1688 parentNucleus = theIonTable->
GetIon(PZ,PA,PE);
1700 if (theDecayChannel == 0) {
1704 G4cerr <<
" G4RadioactiveDecay::DoIt : cannot determine decay channel ";
1705 G4cerr <<
" for this nucleus; decay as if no biasing active ";
1710 tempprods =
DoDecay(*parentNucleus);
1715 tempprods = theDecayChannel->DecayIt(tempmass);
1716 weight *= (theDecayChannel->GetBR())*(decayTable->
entries());
1719 tempprods =
DoDecay(*parentNucleus);
1723 numberOfSecondaries = tempprods->
entries();
1724 currentTime = finalGlobalTime + theDecayTime;
1725 for (index = 0; index < numberOfSecondaries; index++) {
1728 pw.push_back(weight);
1729 ptime.push_back(currentTime);
1730 secondaryparticles.push_back(asecondaryparticle);
1740 totalNumberOfSecondaries = pw.size();
1742 for (index=0; index < totalNumberOfSecondaries; index++) {
1744 ptime[index], currentPosition);
1782 if (theDecayChannel == 0) {
1784 G4cerr <<
"G4RadioactiveDecay::DoIt : can not determine decay channel";
1790 G4cerr <<
"G4RadioactiveDecay::DoIt : selected decay channel addr:";
1795 products = theDecayChannel->
DecayIt(tempmass);
1810 if (0 == products || 0 == products->
entries())
return;
1828 if (daughterType == electron || daughterType == positron ||
1829 daughterType == neutron || daughterType == gamma ||
1837 G4cout <<
"CollimateDecayProduct for daughter "
1862 dir.setPhi(dir.phi()+phi);
1863 dir.setTheta(dir.theta()+std::acos(cosTheta));
1868 G4cout <<
" ChooseCollimationDirection returns " << dir <<
G4endl;
void RegisterIsotopeTable(G4VIsotopeTable *table)
G4double halflifethreshold
G4double GetDaughterExcitation()
G4RadioactiveDecayRateVector theDecayRateTable
void SetBR(G4double value)
static G4LossTableManager * Instance()
G4double GetLocalTime() const
void InitialiseAtomicDeexcitation()
void SelectAVolume(const G4String aVolume)
std::map< G4String, G4DecayTable * > DecayTableMap
std::ostringstream G4ExceptionDescription
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
G4bool IsRateTableReady(const G4ParticleDefinition &)
const G4DynamicParticle * GetDynamicParticle() const
std::vector< G4String > ValidVolumes
G4double HalfLife() const
void SetTaos(std::vector< G4double > value)
G4String strip(G4int strip_Type=trailing, char c=' ')
void SetDecayRateC(std::vector< G4double > value)
void AddUserDecayDataFile(G4int Z, G4int A, G4String filename)
G4bool GetPDGStable() const
G4ParticleDefinition * GetIon(G4int Z, G4int A, G4int lvl=0)
static G4Alpha * Definition()
const G4ThreeVector & GetPosition() const
void AddSecondary(G4Track *aSecondary)
G4TrackStatus GetTrackStatus() const
G4RadioactiveDecayRates theDecayRateVector
static G4Electron * Definition()
void SetTouchableHandle(const G4TouchableHandle &apValue)
G4VDecayChannel * GetDecayChannel(G4int index) const
G4double GetMeanLifeTime(const G4Track &theTrack, G4ForceCondition *condition)
G4int GetVerboseLevel() const
void SetMomentumDirection(const G4ThreeVector &aDirection)
static G4NuclearLevelStore * GetInstance()
G4ParticleDefinition * GetDefinition() const
void SetSourceTimeProfile(G4String filename)
G4bool IsApplicable(const G4ParticleDefinition &)
void Boost(G4double totalEnergy, const G4ThreeVector &momentumDirection)
void SetGeneration(G4int value)
static G4Positron * Definition()
void ClearNumberOfInteractionLengthLeft()
#define G4MUTEX_INITIALIZER
const G4String & GetParticleName() const
void SetWeight(G4double aValue)
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
void SetItsRates(G4RadioactiveDecayRates arate)
G4ParticleChangeForRadDecay fParticleChangeForRadDecay
G4double GetTotalMomentum() const
virtual void Initialize(const G4Track &)
G4IonTable * GetIonTable() const
G4GLOB_DLL std::ostream G4cout
void ProposeLocalTime(G4double t)
static char * currentPosition
static const G4ThreeVector origin
void SetDecayRate(G4int, G4int, G4double, G4int, std::vector< G4double >, std::vector< G4double >)
G4DecayTable * GetDecayTable(G4ParticleDefinition *)
const G4ThreeVector & GetMomentumDirection() const
void CollimateDecayProduct(G4DynamicParticle *product)
static G4LogicalVolumeStore * GetInstance()
G4ThreeVector ChooseCollimationDirection() const
G4int GetDecayTimeBin(const G4double aDecayTime)
void SetProcessSubType(G4int)
void DeselectAllVolumes()
const G4String & GetParticleType() const
G4DecayProducts * DoDecay(G4ParticleDefinition &theParticleDef)
const G4ParticleDefinition * GetParticleDefinition() const
G4double GetGlobalTime() const
static const G4double A[nN]
G4double GetTaoTime(const G4double, const G4double)
G4NucleusLimits theNucleusLimits
const G4TouchableHandle & GetTouchableHandle() const
G4VDecayChannel * SelectADecayChannel()
G4NuclearLevelManager * GetManager(G4int Z, G4int A)
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
G4ParticleDefinition * GetDaughterNucleus()
void Insert(G4VDecayChannel *aChannel)
G4double GetMeanFreePath(const G4Track &theTrack, G4double previousStepSize, G4ForceCondition *condition)
G4RadioactiveDecayRate theDecayRate
static const G4double levelTolerance
G4RadioactiveDecayRateTable theDecayRateTableVector
std::map< G4int, G4String > theUserRadioactiveDataFiles
static const double nanosecond
G4LogicalVolume * GetLogicalVolume() const
G4double forceDecayHalfAngle
G4double GetPDGMass() const
static G4ParticleTable * GetParticleTable()
void GetDecayRateTable(const G4ParticleDefinition &)
G4RadioactiveDecaymessenger * theRadioactiveDecaymessenger
void SetNumberOfSecondaries(G4int totSecondaries)
G4DynamicParticle * PopProducts()
G4VParticleChange * pParticleChange
G4RadioactiveDecayMode GetDecayMode()
void AddDecayRateTable(const G4ParticleDefinition &)
static G4Neutron * Definition()
G4int NumberOfLevels() const
G4VPhysicalVolume * GetVolume() const
std::vector< G4RadioactivityTable * > theRadioactivityTables
G4double GetWeight() const
void SetIonName(G4String name)
G4double GetPDGLifeTime() const
G4RadioactiveDecay(const G4String &processName="RadioactiveDecay")
G4ThreeVector forceDecayDirection
G4VAtomDeexcitation * AtomDeexcitation()
G4DecayTable * LoadDecayTable(G4ParticleDefinition &theParentNucleus)
virtual G4DecayProducts * DecayIt(G4double parentMass=-1.0)=0
void ProposeTrackStatus(G4TrackStatus status)
void BuildPhysicsTable(const G4ParticleDefinition &)
G4double ShapeFactor(const G4BetaDecayType &, const G4double &p_e, const G4double &e_nu)
void SetHLThreshold(G4double hl)
void SetAnalogueMonteCarlo(G4bool r)
static const G4double alpha
const G4NuclearLevel * NearestLevel(G4double energy, G4double eDiffMax=9999.*CLHEP::GeV) const
G4double FermiFunction(const G4double &W)
void DeselectAVolume(const G4String aVolume)
void SetE(G4double value)
G4RIsotopeTable * theIsotopeTable
void AddUserDecayDataFile(G4int Z, G4int A, G4String filename)
void CollimateDecay(G4DecayProducts *products)
void SetDecayBias(G4String filename)
static G4Gamma * Definition()
void SetGoodForTrackingFlag(G4bool value=true)
G4VParticleChange * DecayIt(const G4Track &theTrack, const G4Step &theStep)
G4GLOB_DLL std::ostream G4cerr
G4int GetBaryonNumber() const