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