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