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