Geant4 Cross Reference |
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 // 26 // 27 // ------------------------------------------------------------------- 28 // 29 // GEANT4 header file 30 // 31 // File name: G4NucLevel 32 // 33 // Author: V.Ivanchenko (M.Kelsey reading method is used) 34 // 35 // Creation date: 4 January 2012 36 // 37 // Modifications: 38 // 39 // ------------------------------------------------------------------- 40 41 #include "G4LevelReader.hh" 42 #include "G4NucleiProperties.hh" 43 #include "G4NucLevel.hh" 44 #include "G4NuclearLevelData.hh" 45 #include "G4DeexPrecoParameters.hh" 46 #include "G4SystemOfUnits.hh" 47 #include "G4Pow.hh" 48 #include <fstream> 49 #include <sstream> 50 51 namespace 52 { 53 const G4int countmax = 4; 54 const G4int nfloting = 13; 55 const G4double eTolarence = 2*CLHEP::eV; 56 const G4String fFloatingLevels[13] = { 57 "-", "+X", "+Y", "+Z", "+U", "+V", "+W", "+R", "+S", "+T", "+A", "+B", "+C"}; 58 } 59 60 G4LevelReader::G4LevelReader(G4NuclearLevelData* ptr) 61 : fData(ptr) 62 { 63 fAlphaMax = (G4float)1.e15; 64 fTimeFactor = CLHEP::second/G4Pow::GetInstance()->logZ(2); 65 fDirectory = G4String(G4FindDataDir("G4LEVELGAMMADATA")); 66 67 vTrans.resize(fTransMax,0); 68 vRatio.resize(fTransMax,0.0f); 69 vGammaCumProbability.resize(fTransMax,0.0f); 70 vGammaProbability.resize(fTransMax,0.0f); 71 vShellProbability.resize(fTransMax,nullptr); 72 73 vEnergy.resize(fLevelMax,0.0); 74 vSpin.resize(fLevelMax,0); 75 vLevel.resize(fLevelMax,nullptr); 76 } 77 78 G4bool G4LevelReader::ReadData(std::istringstream& stream, G4double& x) 79 { 80 stream >> x; 81 return !stream.fail(); 82 } 83 84 G4bool G4LevelReader::ReadDataItem(std::istream& dataFile, G4double& x) 85 { 86 x = 0.0; 87 for(G4int i=0; i<nbufmax; ++i) { buffer[i] = ' '; } 88 G4bool okay = true; 89 dataFile >> buffer; 90 if(dataFile.fail()) { okay = false; } 91 else { x = std::strtod(buffer, 0); } 92 return okay; 93 } 94 95 G4bool G4LevelReader::ReadDataItem(std::istream& dataFile, G4float& x) 96 { 97 x = 0.0f; 98 for(G4int i=0; i<nbuf1; ++i) { buff1[i] = ' '; } 99 G4bool okay = true; 100 dataFile >> buff1; 101 if(dataFile.fail()) { okay = false; } 102 else { x = std::atof(buff1); } 103 104 return okay; 105 } 106 107 G4bool G4LevelReader::ReadDataItem(std::istream& dataFile, G4int& ix) 108 { 109 ix = 0; 110 for(G4int i=0; i<nbuf2; ++i) { buff2[i] = ' '; } 111 G4bool okay = true; 112 dataFile >> buff2; 113 if(dataFile.fail()) { okay = false; } 114 else { ix = std::atoi(buff2); } 115 116 return okay; 117 } 118 119 const std::vector<G4float>* G4LevelReader::NormalizedICCProbability(G4int Z) 120 { 121 std::vector<G4float>* vec = nullptr; 122 G4int LL = 3; 123 G4int M = 5; 124 G4int N = 1; 125 if(Z <= 27) { 126 M = N = 0; 127 if(Z <= 4) { 128 LL = 1; 129 } else if(Z <= 6) { 130 LL = 2; 131 } else if(Z <= 10) { 132 } else if(Z <= 12) { 133 M = 1; 134 } else if(Z <= 17) { 135 M = 2; 136 } else if(Z == 18) { 137 M = 3; 138 } else if(Z <= 20) { 139 M = 3; 140 N = 1; 141 } else { 142 M = 4; 143 N = 1; 144 } 145 if(LL < 3) { for(G4int i=LL+1; i<=4; ++i) { fICC[i] = 0.0f; } } 146 if(M < 5) { for(G4int i=M+4; i<=8; ++i) { fICC[i] = 0.0f; } } 147 if(N < 1) { fICC[9] = 0.0f; } 148 } 149 G4float norm = 0.0f; 150 for(G4int i=0; i<10; ++i) { 151 norm += fICC[i]; 152 fICC[i] = norm; 153 } 154 if(norm == 0.0f && fAlpha > 0.0f) { 155 fICC[0] = norm = 1.0f; 156 } 157 if(norm > 0.0f) { 158 norm = 1.0f/norm; 159 vec = new std::vector<G4float>; 160 G4float x; 161 for(G4int i=0; i<10; ++i) { 162 x = fICC[i]*norm; 163 if(x > 0.995f || 9 == i) { 164 vec->push_back(1.0f); 165 break; 166 } 167 vec->push_back(x); 168 } 169 if (fVerbose > 3) { 170 G4long prec = G4cout.precision(3); 171 G4cout << "# InternalConv: "; 172 std::size_t nn = vec->size(); 173 for(std::size_t i=0; i<nn; ++i) { G4cout << " " << (*vec)[i]; } 174 G4cout << G4endl; 175 G4cout.precision(prec); 176 } 177 } 178 return vec; 179 } 180 181 const G4LevelManager* 182 G4LevelReader::CreateLevelManager(G4int Z, G4int A) 183 { 184 std::ostringstream ss; 185 ss << fDirectory << "/z" << Z << ".a" << A; 186 std::ifstream infile(ss.str(), std::ios::in); 187 188 // file is not opened 189 if (!infile.is_open()) { 190 if(fVerbose > 1) { 191 G4ExceptionDescription ed; 192 ed << "Regular file " << ss.str() << " is not opened! Z=" 193 << Z << " A=" << A; 194 G4Exception("G4LevelReader::LevelManager(..)","had014", 195 JustWarning, ed, "Check file path"); 196 } 197 return nullptr; 198 } 199 // file is opened 200 if (fVerbose > 1) { 201 G4cout << "G4LevelReader: open file " << ss.str() << " for Z= " 202 << Z << " A= " << A << G4endl; 203 } 204 return LevelManager(Z, A, infile); 205 } 206 207 const G4LevelManager* 208 G4LevelReader::MakeLevelManager(G4int Z, G4int A, const G4String& filename) 209 { 210 std::ifstream infile(filename, std::ios::in); 211 212 // file is not opened 213 if (!infile.is_open()) { 214 if(fVerbose > 1) { 215 G4ExceptionDescription ed; 216 ed << "External file " << filename << " is not opened! Z=" 217 << Z << " A=" << A; 218 G4Exception("G4LevelReader::LevelManager(..)","had014", 219 FatalException, ed, "Check file path"); 220 } 221 return nullptr; 222 } 223 // file is opened 224 if (fVerbose > 1) { 225 G4cout << "G4LevelReader: open external file " << filename 226 << " for Z= " << Z << " A= " << A << G4endl; 227 } 228 return LevelManager(Z, A, infile); 229 } 230 231 const G4LevelManager* 232 G4LevelReader::LevelManager(G4int Z, G4int A, std::ifstream& infile) 233 { 234 G4bool allLevels = fData->GetParameters()->StoreICLevelData(); 235 fPol = " "; 236 G4int i = 0; 237 for (;;) { 238 infile >> i1 >> fPol; // Level number and floating level 239 // normal end of file 240 if (infile.eof()) { 241 if (fVerbose > 1) { 242 G4cout << "### End of file Z= " << Z << " A= " << A 243 << " Nlevels= " << i << G4endl; 244 } 245 break; 246 } 247 // start reading new level data 248 #ifdef G4VERBOSE 249 if(fVerbose > 2) { 250 G4cout << "New line: i1= " << i1 << " fPol= <" << fPol << "> " << G4endl; 251 } 252 #endif 253 // read new level data 254 if (!(ReadDataItem(infile, fEnergy) && 255 ReadDataItem(infile, fTime) && 256 ReadDataItem(infile, fSpin) && 257 ReadDataItem(infile, ntrans))) { 258 if (fVerbose > 1) { 259 G4cout << "### Incomplete end of file Z= " << Z << " A= " << A 260 << " Nlevels= " << i << G4endl; 261 } 262 break; 263 } 264 fEnergy *= CLHEP::keV; 265 for (k=0; k<nfloting; ++k) { 266 if (fPol == fFloatingLevels[k]) { 267 break; 268 } 269 } 270 // if a previous level has higher energy the current should be ignored 271 // data with wrong level should red anyway 272 if (0 < i) { 273 if (fEnergy < vEnergy[i-1]) { 274 #ifdef G4VERBOSE 275 ++count1; 276 if (count1 < countmax && fVerbose > 0) { 277 G4cout << "### G4LevelReader: broken level " << i 278 << " E(MeV)= " << fEnergy << " < " << vEnergy[i-1] 279 << " for isotope Z= " << Z << " A= " 280 << A << " level energy increased" << G4endl; 281 } 282 #endif 283 // for any case 284 fEnergy = vEnergy[i-1] + eTolarence; 285 } 286 } 287 vEnergy[i] = fEnergy; 288 if (fTime > 0.0) { fTime *= fTimeFactor; } 289 else if (fTime < 0.0) { fTime = DBL_MAX; } 290 291 G4int twos = G4lrint(fSpin + fSpin); 292 twos = std::max(twos, -100); 293 vSpin[i] = 100 + twos + k*100000; 294 #ifdef G4VERBOSE 295 if (fVerbose > 2) { 296 G4cout << " Level #" << i1 << " E(MeV)=" << fEnergy/CLHEP::MeV 297 << " LTime(s)=" << fTime << " 2S=" << vSpin[i] 298 << " meta=" << vSpin[i]/100000 << " idx=" << i 299 << " ntr=" << ntrans << G4endl; 300 } 301 #endif 302 vLevel[i] = nullptr; 303 if (ntrans == 0) { 304 vLevel[i] = new G4NucLevel(0, fTime, vTrans, 305 vGammaCumProbability, 306 vGammaProbability, 307 vRatio, 308 vShellProbability); 309 } else if (ntrans > 0) { 310 311 G4bool isTransOK = true; 312 if (ntrans > fTransMax) { 313 fTransMax = ntrans; 314 vTrans.resize(fTransMax); 315 vRatio.resize(fTransMax); 316 vGammaCumProbability.resize(fTransMax); 317 vGammaProbability.resize(fTransMax); 318 vShellProbability.resize(fTransMax); 319 } 320 fNorm1 = 0.0f; 321 for (G4int j=0; j<ntrans; ++j) { 322 323 if (!(ReadDataItem(infile, i2) && 324 ReadDataItem(infile, fTransEnergy) && 325 ReadDataItem(infile, fProb) && 326 ReadDataItem(infile, tnum) && 327 ReadDataItem(infile, vRatio[j]) && 328 ReadDataItem(infile, fAlpha))) { 329 #ifdef G4VERBOSE 330 ++count2; 331 if (count2 < countmax && fVerbose > 0) { 332 G4cout << "### Fail to read transition j=" << j 333 << " j=" << j << " i2=" << i2 334 << " Z=" << Z << " A=" << A << G4endl; 335 } 336 #endif 337 isTransOK = false; 338 } 339 if (i2 >= i) { 340 #ifdef G4VERBOSE 341 ++count2; 342 if (count2 < countmax) { 343 G4cout << "### G4LevelReader: broken transition " << j 344 << " from level " << i << " to " << i2 345 << " for isotope Z= " << Z << " A= " 346 << A << "; the transition probability set to zero" << G4endl; 347 } 348 #endif 349 isTransOK = false; 350 fProb = 0.0f; 351 } 352 vTrans[j] = i2*10000 + tnum; 353 fAlpha = std::min(std::max(fAlpha,0.f), fAlphaMax); 354 G4float x = 1.0f + fAlpha; 355 fNorm1 += x*fProb; 356 vGammaCumProbability[j] = fNorm1; 357 vGammaProbability[j] = 1.0f/x; 358 vShellProbability[j] = nullptr; 359 if (fVerbose > 2) { 360 G4long prec = G4cout.precision(4); 361 G4cout << "### Transition #" << j << " to level " << i2 362 << " i2= " << i2 << " Etrans(MeV)= " << fTransEnergy*CLHEP::keV 363 << " fProb= " << fProb << " MultiP= " << tnum 364 << " fMpRatio= " << fRatio << " fAlpha= " << fAlpha 365 << G4endl; 366 G4cout.precision(prec); 367 } 368 if (fAlpha > 0.0f) { 369 for (k=0; k<10; ++k) { 370 if (!ReadDataItem(infile,fICC[k])) { 371 isTransOK = false; 372 #ifdef G4VERBOSE 373 ++count2; 374 if (count2 < countmax) { 375 G4cout << "### G4LevelReader: fail to read conversion coeff k= " << k 376 << " for transition j= " << j 377 << " Z= " << Z << " A= " << A << G4endl; 378 } 379 #endif 380 for(kk=k; kk<10; ++kk) { fICC[kk] = 0.f; } 381 } 382 } 383 if (allLevels) { 384 vShellProbability[j] = NormalizedICCProbability(Z); 385 } 386 } 387 } 388 if (ntrans > 0) { 389 G4int nt = ntrans - 1; 390 if (fVerbose > 2) { 391 G4cout << "=== New G4NucLevel: Ntrans=" << ntrans 392 << " Time(ns)=" << fTime 393 << " IdxTrans=" << vTrans[nt]/10000 394 << " isOK=" << isTransOK 395 << G4endl; 396 } 397 if (0.0f < fNorm1) { fNorm1 = 1.0f/fNorm1; } 398 for (k=0; k<nt; ++k) { 399 vGammaCumProbability[k] *= fNorm1; 400 #ifdef G4VERBOSE 401 if (fVerbose > 3) { 402 G4cout << "Probabilities[" << k 403 << "]= " << vGammaCumProbability[k] 404 << " " << vGammaProbability[k] 405 << " idxTrans= " << vTrans[k]/10000 406 << G4endl; 407 } 408 #endif 409 } 410 vGammaCumProbability[nt] = 1.0f; 411 vLevel[i] = new G4NucLevel((std::size_t)ntrans, fTime, vTrans, 412 vGammaCumProbability, 413 vGammaProbability, 414 vRatio, 415 vShellProbability); 416 } 417 } 418 ++i; 419 if (i == fLevelMax) { 420 fLevelMax += 10; 421 vEnergy.resize(fLevelMax, 0.0); 422 vSpin.resize(fLevelMax, 0); 423 vLevel.resize(fLevelMax, nullptr); 424 } 425 } 426 G4LevelManager* lman = nullptr; 427 if (1 <= i) { 428 lman = new G4LevelManager(Z, A, (std::size_t)i, vEnergy, vSpin, vLevel); 429 if (fVerbose > 1) { 430 G4cout << "=== Reader: new manager for Z=" << Z << " A=" << A 431 << " Nlevels=" << i << " E[0]=" 432 << vEnergy[0]/CLHEP::MeV << " MeV E1=" 433 << vEnergy[i-1]/CLHEP::MeV << " MeV" 434 << " count1,2=" << count1 << ", " << count2 << G4endl; 435 } 436 } 437 return lman; 438 } 439