| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 // -------------------------------------------
28 //
29 // GEANT4 header file
30 //
31 // File name: G4NucLevel
32 //
33 // Author: V.Ivanchenko (M.Kelsey
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
58 }
59
60 G4LevelReader::G4LevelReader(G4NuclearLevelDat
61 : fData(ptr)
62 {
63 fAlphaMax = (G4float)1.e15;
64 fTimeFactor = CLHEP::second/G4Pow::GetInstan
65 fDirectory = G4String(G4FindDataDir("G4LEVEL
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::istringstr
79 {
80 stream >> x;
81 return !stream.fail();
82 }
83
84 G4bool G4LevelReader::ReadDataItem(std::istrea
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::istrea
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::istrea
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::Nor
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)
146 if(M < 5) { for(G4int i=M+4; i<=8; ++i)
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
174 G4cout << G4endl;
175 G4cout.precision(prec);
176 }
177 }
178 return vec;
179 }
180
181 const G4LevelManager*
182 G4LevelReader::CreateLevelManager(G4int Z, G4i
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() << " i
193 << Z << " A=" << A;
194 G4Exception("G4LevelReader::LevelManager
195 JustWarning, ed, "Check file path");
196 }
197 return nullptr;
198 }
199 // file is opened
200 if (fVerbose > 1) {
201 G4cout << "G4LevelReader: open file " << s
202 << Z << " A= " << A << G4endl;
203 }
204 return LevelManager(Z, A, infile);
205 }
206
207 const G4LevelManager*
208 G4LevelReader::MakeLevelManager(G4int Z, G4int
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 << "
217 << Z << " A=" << A;
218 G4Exception("G4LevelReader::LevelManager
219 FatalException, ed, "Check file path");
220 }
221 return nullptr;
222 }
223 // file is opened
224 if (fVerbose > 1) {
225 G4cout << "G4LevelReader: open external fi
226 << " for Z= " << Z << " A= " << A <
227 }
228 return LevelManager(Z, A, infile);
229 }
230
231 const G4LevelManager*
232 G4LevelReader::LevelManager(G4int Z, G4int A,
233 {
234 G4bool allLevels = fData->GetParameters()->S
235 fPol = " ";
236 G4int i = 0;
237 for (;;) {
238 infile >> i1 >> fPol; // Level number a
239 // normal end of file
240 if (infile.eof()) {
241 if (fVerbose > 1) {
242 G4cout << "### End of file Z= " << Z << " 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 << " f
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= " <
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 t
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
278 << " E(MeV)= " << fEnergy << " < " << vEn
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 }
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)
297 << " LTime(s)=" << fTime << " 2S=" <<
298 << " meta=" << vSpin[i]/100000 << " id
299 << " ntr=" << ntrans << G4endl;
300 }
301 #endif
302 vLevel[i] = nullptr;
303 if (ntrans == 0) {
304 vLevel[i] = new G4NucLevel(0, fTime, vTr
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
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 tra
344 << " from level " << i << " to " << i2
345 << " for isotope Z= " << Z << " A= "
346 << A << "; the transition probability s
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)
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
362 << " i2= " << i2 << " Etrans(MeV)= " <<
363 << " fProb= " << fProb << " MultiP= " <<
364 << " fMpRatio= " << fRatio << " 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
376 << " for transition j= " << j
377 << " Z= " << Z << " A= " << A << G
378 }
379 #endif
380 for(kk=k; kk<10; ++kk) { fICC[kk] = 0.
381 }
382 }
383 if (allLevels) {
384 vShellProbability[j] = NormalizedICCProb
385 }
386 }
387 }
388 if (ntrans > 0) {
389 G4int nt = ntrans - 1;
390 if (fVerbose > 2) {
391 G4cout << "=== New G4NucLevel: Ntran
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)ntra
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
429 if (fVerbose > 1) {
430 G4cout << "=== Reader: new manager for Z
431 << " Nlevels=" << i << " E[0]="
432 << vEnergy[0]/CLHEP::MeV << " MeV E1=
433 << vEnergy[i-1]/CLHEP::MeV << " MeV"
434 << " count1,2=" << count1 << ", " << co
435 }
436 }
437 return lman;
438 }
439