| 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 // $Id: G4LevelReader.cc 88407 2015-02-18 09:18:44Z vnivanch $
26 // 27 //
27 // ------------------------------------------- 28 // -------------------------------------------------------------------
28 // 29 //
29 // GEANT4 header file 30 // GEANT4 header file
30 // 31 //
31 // File name: G4NucLevel 32 // File name: G4NucLevel
32 // 33 //
33 // Author: V.Ivanchenko (M.Kelsey 34 // Author: V.Ivanchenko (M.Kelsey reading method is used)
34 // 35 //
35 // Creation date: 4 January 2012 36 // Creation date: 4 January 2012
36 // 37 //
37 // Modifications: 38 // Modifications:
38 // 39 //
39 // ------------------------------------------- 40 // -------------------------------------------------------------------
40 41
41 #include "G4LevelReader.hh" 42 #include "G4LevelReader.hh"
42 #include "G4NucleiProperties.hh" 43 #include "G4NucleiProperties.hh"
43 #include "G4NucLevel.hh" 44 #include "G4NucLevel.hh"
44 #include "G4NuclearLevelData.hh" <<
45 #include "G4DeexPrecoParameters.hh" <<
46 #include "G4SystemOfUnits.hh" 45 #include "G4SystemOfUnits.hh"
47 #include "G4Pow.hh" 46 #include "G4Pow.hh"
>> 47 #include <vector>
48 #include <fstream> 48 #include <fstream>
49 #include <sstream> 49 #include <sstream>
50 50
51 namespace << 51 static const G4float x_energy = 0.000001;
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 52
60 G4LevelReader::G4LevelReader(G4NuclearLevelDat << 53 G4String G4LevelReader::fTrans[] = {
61 : fData(ptr) << 54 "1-", "1+", "2-", "2+", "3-", "3+", "4-", "4+", "5-", "5+"};
>> 55
>> 56 G4LevelReader::G4LevelReader()
>> 57 : fMinProbability(1.e-8),fVerbose(0)
62 { 58 {
63 fAlphaMax = (G4float)1.e15; <<
64 fTimeFactor = CLHEP::second/G4Pow::GetInstan 59 fTimeFactor = CLHEP::second/G4Pow::GetInstance()->logZ(2);
65 fDirectory = G4String(G4FindDataDir("G4LEVEL << 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
>> 74 fEnergy = fCurrEnergy = fTrEnergy = fProb = fTime = fSpin = fAlpha = 0.0;
>> 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
>> 93 G4LevelReader::~G4LevelReader()
>> 94 {}
66 95
67 vTrans.resize(fTransMax,0); << 96 const G4LevelManager*
68 vRatio.resize(fTransMax,0.0f); << 97 G4LevelReader::CreateLevelManager(G4int Z, G4int A)
69 vGammaCumProbability.resize(fTransMax,0.0f); << 98 {
70 vGammaProbability.resize(fTransMax,0.0f); << 99 std::ostringstream ss;
71 vShellProbability.resize(fTransMax,nullptr); << 100 ss << "/z" << Z << ".a" << A;
72 << 101 G4String st = G4String(ss.str());
73 vEnergy.resize(fLevelMax,0.0); << 102 fFile = fDirectory + st;
74 vSpin.resize(fLevelMax,0); << 103 return MakeLevelManager(Z, A, fFile);
75 vLevel.resize(fLevelMax,nullptr); <<
76 } 104 }
77 105
78 G4bool G4LevelReader::ReadData(std::istringstr << 106 const G4LevelManager*
>> 107 G4LevelReader::MakeLevelManager(G4int Z, G4int A, const G4String& filename)
79 { 108 {
80 stream >> x; << 109 vEnergy.clear();
81 return !stream.fail(); << 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;
>> 142 fCurrEnergy = DBL_MAX;
>> 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;
>> 155 fTrEnergy *= 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) {
>> 189 vGammaCumProbability[i] *= fNorm1;
>> 190 vGammaECumProbability[i] *= fNorm2;
>> 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,
>> 214 vGammaCumProbability,
>> 215 vGammaECumProbability,
>> 216 vGammaProbability,
>> 217 vTrans,
>> 218 vShellProbability));
>> 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 }
>> 228 fCurrEnergy = fEnergy;
>> 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; }
>> 237 fProb = std::max(fProb, fMinProbability);
>> 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) {
>> 279 vShellProbability.push_back(NormalizedICCProbability(Z));
>> 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) {
>> 297 man = new G4LevelManager(vEnergy,vTime,vTimeg,vSpin,vLevel);
>> 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;
82 } 309 }
83 310
84 G4bool G4LevelReader::ReadDataItem(std::istrea 311 G4bool G4LevelReader::ReadDataItem(std::istream& dataFile, G4double& x)
85 { 312 {
86 x = 0.0; 313 x = 0.0;
87 for(G4int i=0; i<nbufmax; ++i) { buffer[i] = << 314 for(G4int i=0; i<20; ++i) { buffer[i] = ' '; }
88 G4bool okay = true; 315 G4bool okay = true;
89 dataFile >> buffer; 316 dataFile >> buffer;
90 if(dataFile.fail()) { okay = false; } 317 if(dataFile.fail()) { okay = false; }
91 else { x = std::strtod(buffer, 0); } << 318 else { x = 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 319
104 return okay; 320 return okay;
105 } 321 }
106 322
107 G4bool G4LevelReader::ReadDataItem(std::istrea << 323 G4bool G4LevelReader::ReadDataItem(std::istream& dataFile, G4String& x)
108 { 324 {
109 ix = 0; <<
110 for(G4int i=0; i<nbuf2; ++i) { buff2[i] = ' <<
111 G4bool okay = true; 325 G4bool okay = true;
112 dataFile >> buff2; << 326 dataFile >> bufp;
113 if(dataFile.fail()) { okay = false; } 327 if(dataFile.fail()) { okay = false; }
114 else { ix = std::atoi(buff2); } << 328 else { x = G4String(bufp, 2); }
115 <<
116 return okay; 329 return okay;
117 } 330 }
118 331
119 const std::vector<G4float>* G4LevelReader::Nor 332 const std::vector<G4float>* G4LevelReader::NormalizedICCProbability(G4int Z)
120 { 333 {
121 std::vector<G4float>* vec = nullptr; 334 std::vector<G4float>* vec = nullptr;
122 G4int LL = 3; 335 G4int LL = 3;
123 G4int M = 5; 336 G4int M = 5;
124 G4int N = 1; 337 G4int N = 1;
125 if(Z <= 27) { << 338 if(Z <= 4) {
>> 339 LL = 1;
126 M = N = 0; 340 M = N = 0;
127 if(Z <= 4) { << 341 } else if(Z <= 6) {
128 LL = 1; << 342 LL = 2;
129 } else if(Z <= 6) { << 343 M = N = 0;
130 LL = 2; << 344 } else if(Z <= 10) {
131 } else if(Z <= 10) { << 345 M = N = 0;
132 } else if(Z <= 12) { << 346 } else if(Z <= 12) {
133 M = 1; << 347 M = 1;
134 } else if(Z <= 17) { << 348 N = 0;
135 M = 2; << 349 } else if(Z <= 17) {
136 } else if(Z == 18) { << 350 M = 2;
137 M = 3; << 351 N = 0;
138 } else if(Z <= 20) { << 352 } else if(Z == 18) {
139 M = 3; << 353 M = 3;
140 N = 1; << 354 N = 0;
141 } else { << 355 } else if(Z <= 20) {
142 M = 4; << 356 M = 3;
143 N = 1; << 357 } else if(Z <= 27) {
144 } << 358 M = 4;
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 } 359 }
149 G4float norm = 0.0f; << 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; }
150 for(G4int i=0; i<10; ++i) { 364 for(G4int i=0; i<10; ++i) {
151 norm += fICC[i]; 365 norm += fICC[i];
152 fICC[i] = norm; 366 fICC[i] = norm;
153 } 367 }
154 if(norm == 0.0f && fAlpha > 0.0f) { << 368 if(norm > 0.0) {
155 fICC[0] = norm = 1.0f; << 369 norm = 1.0/norm;
156 } <<
157 if(norm > 0.0f) { <<
158 norm = 1.0f/norm; <<
159 vec = new std::vector<G4float>; 370 vec = new std::vector<G4float>;
160 G4float x; << 371 vec->reserve(10);
161 for(G4int i=0; i<10; ++i) { << 372 for(G4int i=0; i<9; ++i) {
162 x = fICC[i]*norm; << 373 vec->push_back((G4float)(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 } 374 }
>> 375 vec->push_back(1.0f);
177 } 376 }
178 return vec; 377 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 } 378 }
439 379