Geant4  10.02.p02
G4LevelReader.cc
Go to the documentation of this file.
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 // $Id: G4LevelReader.cc 88407 2015-02-18 09:18:44Z vnivanch $
27 //
28 // -------------------------------------------------------------------
29 //
30 // GEANT4 header file
31 //
32 // File name: G4NucLevel
33 //
34 // Author: V.Ivanchenko (M.Kelsey reading method is used)
35 //
36 // Creation date: 4 January 2012
37 //
38 // Modifications:
39 //
40 // -------------------------------------------------------------------
41 
42 #include "G4LevelReader.hh"
43 #include "G4NucleiProperties.hh"
44 #include "G4NucLevel.hh"
45 #include "G4SystemOfUnits.hh"
46 #include "G4Pow.hh"
47 #include <vector>
48 #include <fstream>
49 #include <sstream>
50 
51 static const G4float x_energy = 0.000001;
52 
54  "1-", "1+", "2-", "2+", "3-", "3+", "4-", "4+", "5-", "5+"};
55 
57  : fMinProbability(1.e-8),fVerbose(0)
58 {
60  char* directory = getenv("G4LEVELGAMMADATA");
61  if(directory) {
62  fDirectory = directory;
63  } else {
64  G4Exception("G4LevelReader()","had0707",FatalException,
65  "Environment variable G4LEVELGAMMADATA is not defined");
66  fDirectory = "";
67  }
68  fFile = fDirectory + "/z100.a200";
69  fPol = " ";
70  for(G4int i=0; i<10; ++i) { fICC[i] = 0.0; }
71  for(G4int i=0; i<20; ++i) { buffer[i] = ' '; }
72  bufp[0] = bufp[1] = ' ';
73 
75  fNorm1 = fNorm2 = 0.0f;
76 
77  size_t nn = 10;
78  vTransEnergy.reserve(nn);
79  vGammaCumProbability.reserve(nn);
80  vGammaECumProbability.reserve(nn);
81  vGammaProbability.reserve(nn);
82  vShellProbability.reserve(nn);
83  vTrans.reserve(nn);
84 
85  nn = 100;
86  vEnergy.reserve(nn);
87  vTime.reserve(nn);
88  vTimeg.reserve(nn);
89  vSpin.reserve(nn);
90  vLevel.reserve(nn);
91 }
92 
94 {}
95 
96 const G4LevelManager*
98 {
99  std::ostringstream ss;
100  ss << "/z" << Z << ".a" << A;
101  G4String st = G4String(ss.str());
102  fFile = fDirectory + st;
103  return MakeLevelManager(Z, A, fFile);
104 }
105 
106 const G4LevelManager*
108 {
109  vEnergy.clear();
110  vTime.clear();
111  vTimeg.clear();
112  vSpin.clear();
113  vLevel.clear();
114 
115  vEnergy.push_back(0.0f);
116  vTime.push_back(FLT_MAX);
117  vTimeg.push_back(FLT_MAX);
118  vSpin.push_back(0);
119  vLevel.push_back(nullptr);
120 
121  std::ifstream infile(filename, std::ios::in);
122 
123  // file is not opened
124  if (!infile.is_open()) {
125  if (fVerbose > 0) {
126  G4cout << " G4LevelReader: fail open file for Z= "
127  << Z << " A= " << A
128  << " <" << filename << ">" << G4endl;
129  }
130 
131  } else {
132 
133  if (fVerbose > 0) {
134  G4cout << "G4LevelReader: open file for Z= "
135  << Z << " A= " << A
136  << " <" << filename << ">" << G4endl;
137  }
138  // read line by line
139  G4bool end = false;
140  G4bool next = true;
141  G4int nline = -1;
143  do {
144  fPol = " ";
145  ++nline;
146  G4bool isOK = (ReadDataItem(infile,fEnergy) &&
147  ReadDataItem(infile,fTrEnergy) &&
148  ReadDataItem(infile,fProb) &&
149  ReadDataItem(infile,fPol) &&
150  ReadDataItem(infile,fTime) &&
151  ReadDataItem(infile,fSpin) &&
152  ReadDataItem(infile,fAlpha));
153 
154  fEnergy *= CLHEP::keV;
156 
157  if(isOK) {
158  for(G4int i=0; i<10; ++i) {
159  isOK = (isOK && (ReadDataItem(infile,fICC[i])));
160  }
161  }
162  if(!isOK) {
163  end = true;
164  next = false;
165  }
166  if(fVerbose > 1) {
167  G4cout << "#Line " << nline << " " << fEnergy << " " << fTrEnergy
168  << " " << fProb << " " << fPol << " " << fTime << " "
169  << fSpin << " " << fAlpha << G4endl;
170  G4cout << " ";
171  for(G4int i=0; i<10; ++i) { G4cout << fICC[i] << " "; }
172  G4cout << G4endl;
173  }
174  // end of nuclear level data
175  if(end || fEnergy > fCurrEnergy) {
176  size_t nn = vTransEnergy.size();
177  if(nn > 0) {
178  if(fVerbose > 1) {
179  G4cout << "Reader: new level E= " << fCurrEnergy
180  << " Ntransitions= " << nn << " fNorm1= " << fNorm1
181  << " fNorm2= " << fNorm2 << G4endl;
182  }
183  if(fNorm1 > 0.0f) {
184  --nn;
185  fNorm1 = 1.0f/fNorm1;
186  vTimeg.push_back(((G4float)(fTime*fTimeFactor))*fNorm2*fNorm1);
187  fNorm2 = 1.0f/fNorm2;
188  for(size_t i=0; i<nn; ++i) {
191  if(fVerbose > 2) {
192  G4cout << "Probabilities[" << i << "]= " << vGammaCumProbability[i]
193  << " " << vGammaECumProbability[i]
194  << " Etran= " << vTransEnergy[i]
195  << G4endl;
196  }
197  }
198  vGammaCumProbability[nn] = 1.0f;
199  vGammaECumProbability[nn] = 1.0f;
200  if(fVerbose > 2) {
201  G4cout << "Probabilities[" << nn << "]= " << vGammaCumProbability[nn]
202  << " " << vGammaECumProbability[nn]
203  << " Etran= " << vTransEnergy[nn]
204  << G4endl;
205  }
206  //case of X-level
207  } else {
208  vTimeg.push_back(0.0f);
209  vGammaCumProbability[0] = 0.0f;
210  vGammaECumProbability[0] = 0.0f;
211  }
212 
213  vLevel.push_back(new G4NucLevel(vTransEnergy,
217  vTrans,
219  vTransEnergy.clear();
220  vGammaCumProbability.clear();
221  vGammaECumProbability.clear();
222  vGammaProbability.clear();
223  vShellProbability.clear();
224  vTrans.clear();
225  }
226  if(!end) { next = true; }
227  }
229  // begin nuclear level data
230  if(next) {
231  if(fVerbose > 2) {
232  G4cout << "== Reader: begin of new level E= " << fEnergy << G4endl;
233  }
234  vEnergy.push_back((G4float)fEnergy);
235  vTime.push_back((G4float)(fTime*fTimeFactor));
236  if(fSpin > 20.0) { fSpin = 0.0; }
238  vSpin.push_back(G4lrint(2*fSpin));
239  fNorm1 = 0.0f;
240  fNorm2 = 0.0f;
241  next = false;
242  }
243  // continue filling level data
244  if(!end) {
245  if(fProb > 0.0) {
246  // by default transition to a ground state
247  G4float efinal = (G4float)(fEnergy - fTrEnergy);
248  G4float elevel = 0.0f;
249  // do not check initial energy
250  G4int nn = vEnergy.size() - 1;
251  if(0 < nn) {
252  G4float ediffMin = std::abs(efinal);
253  for(G4int i=0; i<nn; ++i) {
254  G4float ediff = std::abs(efinal - vEnergy[i]);
255  /*
256  G4cout << "Elevel[" << i << "]= " << vEnergy[i]
257  << " Efinal= " << efinal
258  << " Ediff= " << ediff
259  << " EdiffMin= " << ediffMin << G4endl;
260  */
261  if(ediff < ediffMin) {
262  ediffMin = ediff;
263  elevel = vEnergy[i];
264  }
265  }
266  }
267  G4double x = 1.0 + fAlpha;
268  fNorm1 += (G4float)fProb;
269  fNorm2 += (G4float)(fProb*x);
270  vTransEnergy.push_back(elevel);
271  vGammaCumProbability.push_back(fNorm1);
272  vGammaECumProbability.push_back(fNorm2);
273  vGammaProbability.push_back((G4float)(1.0/x));
274  G4int tnum = 0;
275  if(std::abs(vEnergy[nn] - elevel) < x_energy) { tnum = 1; }
276  //for(; tnum<10; ++tnum) { if(fTrans[tnum] == fPol) { break; } }
277  vTrans.push_back(tnum);
278  if(fAlpha > 0.0) {
280  } else {
281  vShellProbability.push_back(nullptr);
282  }
283  }
284  }
285  if(nline > 10000) {
286  G4cout << "G4LevelReader: Line #" << nline << " Z= " << Z << " A= "
287  << " this file is too long - stop loading" << G4endl;
288  end = true;
289  }
290  // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
291  } while (!end);
292  infile.close();
293  }
294 
295  G4LevelManager* man = nullptr;
296  if(vEnergy.size() >= 2) {
298  if(fVerbose > 0) {
299  G4cout << "=== Reader: new manager for Z= " << Z << " A= " << A
300  << " Nlevels= " << vEnergy.size() << " E[0]= "
301  << vEnergy[0]/CLHEP::MeV << " MeV Emax= "
302  << man->MaxLevelEnergy()/CLHEP::MeV << " MeV "
303  << " S: " << vEnergy.size() << " " << vTime.size()
304  << " " << vTimeg.size() << " " << vSpin.size() << " " << vLevel.size()
305  << G4endl;
306  }
307  }
308  return man;
309 }
310 
311 G4bool G4LevelReader::ReadDataItem(std::istream& dataFile, G4double& x)
312 {
313  x = 0.0;
314  for(G4int i=0; i<20; ++i) { buffer[i] = ' '; }
315  G4bool okay = true;
316  dataFile >> buffer;
317  if(dataFile.fail()) { okay = false; }
318  else { x = strtod(buffer, 0); }
319 
320  return okay;
321 }
322 
323 G4bool G4LevelReader::ReadDataItem(std::istream& dataFile, G4String& x)
324 {
325  G4bool okay = true;
326  dataFile >> bufp;
327  if(dataFile.fail()) { okay = false; }
328  else { x = G4String(bufp, 2); }
329  return okay;
330 }
331 
332 const std::vector<G4float>* G4LevelReader::NormalizedICCProbability(G4int Z)
333 {
334  std::vector<G4float>* vec = nullptr;
335  G4int LL = 3;
336  G4int M = 5;
337  G4int N = 1;
338  if(Z <= 4) {
339  LL = 1;
340  M = N = 0;
341  } else if(Z <= 6) {
342  LL = 2;
343  M = N = 0;
344  } else if(Z <= 10) {
345  M = N = 0;
346  } else if(Z <= 12) {
347  M = 1;
348  N = 0;
349  } else if(Z <= 17) {
350  M = 2;
351  N = 0;
352  } else if(Z == 18) {
353  M = 3;
354  N = 0;
355  } else if(Z <= 20) {
356  M = 3;
357  } else if(Z <= 27) {
358  M = 4;
359  }
360  G4double norm = 0.0;
361  if(LL < 3) { for(G4int i=LL+1; i<=4; ++i) { fICC[i] = 0.0; } }
362  if(M < 5) { for(G4int i=M+4; i<=8; ++i) { fICC[i] = 0.0; } }
363  if(N < 1) { fICC[9] = 0.0; }
364  for(G4int i=0; i<10; ++i) {
365  norm += fICC[i];
366  fICC[i] = norm;
367  }
368  if(norm > 0.0) {
369  norm = 1.0/norm;
370  vec = new std::vector<G4float>;
371  vec->reserve(10);
372  for(G4int i=0; i<9; ++i) {
373  vec->push_back((G4float)(fICC[i]*norm));
374  }
375  vec->push_back(1.0f);
376  }
377  return vec;
378 }
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
static const double MeV
Definition: G4SIunits.hh:211
std::vector< G4int > vSpin
G4double fCurrEnergy
G4double fICC[10]
std::vector< const G4NucLevel * > vLevel
float G4float
Definition: G4Types.hh:77
const std::vector< G4float > * NormalizedICCProbability(G4int Z)
std::vector< G4float > vTransEnergy
const G4LevelManager * MakeLevelManager(G4int Z, G4int A, const G4String &filename)
G4double fEnergy
int G4int
Definition: G4Types.hh:78
G4String fFile
G4double logZ(G4int Z) const
Definition: G4Pow.hh:166
std::vector< G4float > vTime
G4GLOB_DLL std::ostream G4cout
double A(double temperature)
bool G4bool
Definition: G4Types.hh:79
G4bool ReadDataItem(std::istream &dataFile, G4double &x)
std::vector< G4float > vTimeg
static const double second
Definition: G4SIunits.hh:156
G4double fAlpha
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
static G4String fTrans[10]
std::vector< const std::vector< G4float > * > vShellProbability
G4double fTime
int G4lrint(double ad)
Definition: templates.hh:163
T max(const T t1, const T t2)
brief Return the largest of the two arguments
std::vector< G4int > vTrans
const G4double x[NPOINTSGL]
static const G4int LL[nN]
std::vector< G4float > vGammaCumProbability
G4double fSpin
#define G4endl
Definition: G4ios.hh:61
G4double fTrEnergy
G4String fDirectory
static const double keV
Definition: G4SIunits.hh:213
const G4LevelManager * CreateLevelManager(G4int Z, G4int A)
std::vector< G4float > vGammaECumProbability
double G4double
Definition: G4Types.hh:76
static const G4float x_energy
G4double fTimeFactor
#define FLT_MAX
Definition: templates.hh:99
#define DBL_MAX
Definition: templates.hh:83
G4float MaxLevelEnergy() const
std::vector< G4float > vEnergy
std::vector< G4float > vGammaProbability
G4double fMinProbability
G4double fProb