00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00055
00056 #include "G4DecayTable.hh"
00057 #include "G4SystemOfUnits.hh"
00058 #include "G4ParticleTable.hh"
00059 #include "G4IsotopeProperty.hh"
00060 #include "G4RIsotopeTable.hh"
00061
00062 #include "G4HadronicException.hh"
00063 #include "G4NuclearLevelStore.hh"
00064
00065
00066
00067
00068
00069
00070
00071
00072
00073
00074 #include "G4ios.hh"
00075 #include "globals.hh"
00076 #include <iomanip>
00077 #include <fstream>
00078 #include <sstream>
00079
00080 const G4double G4RIsotopeTable::levelTolerance = 2.0*keV;
00081
00083
00084 G4RIsotopeTable::G4RIsotopeTable()
00085 {
00086
00087
00088 theUserRadioactiveDataFiles.clear();
00089 }
00090
00092
00093 G4RIsotopeTable::~G4RIsotopeTable()
00094 {
00095 fIsotopeList.clear();
00096 fIsotopeNameList.clear();
00097 }
00099
00100 G4int G4RIsotopeTable::GetVerboseLevel() const
00101 {
00102 return G4ParticleTable::GetParticleTable()->GetVerboseLevel();
00103 }
00105
00106 G4bool G4RIsotopeTable::FindIsotope(G4IsotopeProperty* )
00107 {
00108
00109
00110 return true;
00111 }
00113
00114 G4IsotopeProperty* G4RIsotopeTable::GetIsotope(G4int Z, G4int A, G4double E)
00115 {
00116 G4String fname = GetIsotopeName(Z, A, E);
00117 G4int j = -1;
00118 for (G4int i = 0 ; i< Entries(); i++) {
00119 if(fIsotopeNameList[i] == fname) j = i;}
00120 if (j >=0) {
00121 if (GetVerboseLevel()>1) {
00122 G4cout <<"G4RIsotopeTable::GetIsotope No. : ";
00123 G4cout <<j<<G4endl;
00124 }
00125 return GetIsotope(j);}
00126
00127 else{
00128 G4double meanlife = GetMeanLifeTime(Z, A, E);
00129
00130
00131
00132
00133
00134
00135 G4IsotopeProperty* fProperty = new G4IsotopeProperty();
00136
00137 fProperty->SetLifeTime(meanlife);
00138 fProperty->SetAtomicNumber(Z);
00139 fProperty->SetAtomicMass(A);
00140
00141 fProperty->SetEnergy(E);
00142
00143 fProperty->SetiSpin(0);
00144
00145 fProperty->SetDecayTable(0);
00146
00147 fIsotopeList.push_back(fProperty);
00148 fname = GetIsotopeName(Z, A, E);
00149 fIsotopeNameList.push_back(fname);
00150 if (GetVerboseLevel()>1) {
00151 G4cout <<"G4RIsotopeTable::GetIsotope create: ";
00152 G4cout <<fname <<G4endl;
00153 }
00154 return fProperty;
00155
00156 }
00157 }
00159
00160 G4String G4RIsotopeTable::GetIsotopeName(G4int Z, G4int A, G4double E)
00161 {
00162 std::ostringstream os;
00163 os.setf(std::ios::fixed);
00164 os <<"A"<< A << "Z" << Z <<'[' << std::setprecision(1) << E/keV << ']';
00165 G4String name = os.str();
00166 if (GetVerboseLevel()>1) {
00167 G4cout <<"G4RIsotopeTable::GetIsotope Name: ";
00168 G4cout <<name <<G4endl;
00169 }
00170 return name;
00171 }
00172
00173
00174 G4double G4RIsotopeTable::GetMeanLifeTime(G4int Z, G4int A, G4double& aE)
00175 {
00176
00177 G4double lifetime = -1.0;
00178
00179
00180
00181 G4String file= theUserRadioactiveDataFiles[1000*A+Z];
00182 if (file ==""){
00183 if (!getenv("G4RADIOACTIVEDATA")) {
00184 G4cout << "Please setenv G4RADIOACTIVEDATA to point to the radioactive decay data files." << G4endl;
00185 throw G4HadronicException(__FILE__, __LINE__,
00186 "Please setenv G4RADIOACTIVEDATA to point to the radioactive decay data files.");
00187 }
00188 G4String dirName = getenv("G4RADIOACTIVEDATA");
00189
00190 std::ostringstream os;
00191 os <<dirName <<"/z" <<Z <<".a" <<A ;
00192 file = os.str();
00193 }
00194 std::ifstream DecaySchemeFile(file);
00195
00196 G4bool found_in_raddecay_data(false);
00197 if (!DecaySchemeFile) {
00198 if (GetVerboseLevel()>1) {
00199 G4cout <<"G4RIsotopeTable::GetMeanLife() : "
00200 <<"cannot find ion radioactive decay file: "
00201 <<file <<G4endl;
00202 }
00203 } else {
00204 char inputChars[100]={' '};
00205 G4String inputLine;
00206 G4String recordType("");
00207 G4double a(0.0);
00208 G4double b(0.0);
00209
00210 while (!found_in_raddecay_data && !DecaySchemeFile.getline(inputChars, 100).eof()) {
00211 inputLine = inputChars;
00212 inputLine = inputLine.strip(1);
00213
00214 if (inputChars[0] != '#' && inputLine.length() != 0) {
00215 std::istringstream tmpstream(inputLine);
00216 tmpstream >> recordType >> a >> b;
00217 if (recordType == "P") {
00218 if (std::abs(a*keV-aE) < levelTolerance) {
00219 found_in_raddecay_data = true;
00220 lifetime = b/0.693147*s ;
00221 }
00222 }
00223 }
00224 }
00225 DecaySchemeFile.close();
00226 }
00227
00228 if (!found_in_raddecay_data && aE) {
00229 G4double half_life=-1.;
00230 lifetime = 1.0E-20*s;
00231
00232
00233
00234
00235 const G4NuclearLevel* aLevel =
00236 G4NuclearLevelStore::GetInstance()->GetManager(Z, A)
00237 ->NearestLevel(aE,levelTolerance);
00238 if (aLevel) {
00239 half_life = aLevel->HalfLife();
00240 lifetime = half_life/0.693147;
00241 }
00242
00243 if (GetVerboseLevel()>1 && half_life<0) {
00244 G4cout << "G4RIsotopeTable::GetMeanLife() : ";
00245 G4cout << "cannot find ion of required excitation E = " << aE << G4endl;
00246 G4cout << "state in radioactive or photoevaporation data file " << G4endl;
00247 G4cout <<"The nucleus is assumed to be IT decayed with life = 1E-20 s" << G4endl;
00248 G4cout <<" -----------* THIS MAY CAUSE PROBLEM IN ITS DECAY-----------" << G4endl;
00249 }
00250 }
00251
00252 if (!found_in_raddecay_data && !aE) {
00253 if (GetVerboseLevel()>1) {
00254 G4cout <<"G4RIsotopeTable::GetMeanLife() : ";
00255 G4cout <<"cannot find ion of required excitation E = " << aE << G4endl;
00256 G4cout <<"state in radioactive or photoevaporation data file" <<G4endl;
00257 G4cout <<"The nucleus is assumed to be stable" <<G4endl;
00258 lifetime = -1.0;
00259 }
00260 }
00261
00262 if (GetVerboseLevel()>1) {
00263 G4cout <<"G4RIsotopeTable::GetMeanLifeTime: ";
00264 G4cout <<lifetime << " for " << GetIsotopeName(Z, A, aE) <<G4endl;
00265 }
00266 return lifetime;
00267 }
00269
00270 void G4RIsotopeTable::AddUserDecayDataFile(G4int Z, G4int A,G4String filename)
00271 { if (Z<1 || A<2) {
00272 G4cout<<"Z and A not valid!"<<G4endl;
00273 }
00274
00275 std::ifstream DecaySchemeFile(filename);
00276 if (DecaySchemeFile){
00277 G4int ID_ion=A*1000+Z;
00278 theUserRadioactiveDataFiles[ID_ion]=filename;
00279 }
00280 else {
00281 G4cout<<"The file "<<filename<<" does not exist!"<<G4endl;
00282 }
00283 }
00284