| 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 "G4NuclearLevelData.hh"
45 #include "G4DeexPrecoParameters.hh" 46 #include "G4DeexPrecoParameters.hh"
46 #include "G4SystemOfUnits.hh" 47 #include "G4SystemOfUnits.hh"
47 #include "G4Pow.hh" 48 #include "G4Pow.hh"
>> 49 #include <vector>
48 #include <fstream> 50 #include <fstream>
49 #include <sstream> 51 #include <sstream>
50 52
51 namespace << 53 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 54 "-", "+X", "+Y", "+Z", "+U", "+V", "+W", "+R", "+S", "+T", "+A", "+B", "+C"};
58 } <<
59 55
60 G4LevelReader::G4LevelReader(G4NuclearLevelDat 56 G4LevelReader::G4LevelReader(G4NuclearLevelData* ptr)
61 : fData(ptr) << 57 : fData(ptr),fMinProbability(1.e-8),fVerbose(0),fLevelMax(632),fTransMax(30)
62 { 58 {
63 fAlphaMax = (G4float)1.e15; << 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.0; }
>> 72 for(G4int i=0; i<nbufmax; ++i) { buffer[i] = ' '; }
>> 73 for(G4int i=0; i<nbuf2; ++i) { buff2[i] = ' '; }
>> 74 bufp[0] = bufp[1] = ' ';
66 75
>> 76 fEnergy = fCurrEnergy = fTrEnergy = fProb = fTime =
>> 77 fSpin = fAlpha = fRatio = 0.0;
>> 78 fNorm1 = fNorm2 = 0.0f;
>> 79
>> 80 vIndex.resize(fTransMax,0);
67 vTrans.resize(fTransMax,0); 81 vTrans.resize(fTransMax,0);
68 vRatio.resize(fTransMax,0.0f); 82 vRatio.resize(fTransMax,0.0f);
69 vGammaCumProbability.resize(fTransMax,0.0f); 83 vGammaCumProbability.resize(fTransMax,0.0f);
>> 84 vGammaECumProbability.resize(fTransMax,0.0f);
70 vGammaProbability.resize(fTransMax,0.0f); 85 vGammaProbability.resize(fTransMax,0.0f);
71 vShellProbability.resize(fTransMax,nullptr); 86 vShellProbability.resize(fTransMax,nullptr);
>> 87 vMpRatio.resize(fTransMax,0.0f);
72 88
73 vEnergy.resize(fLevelMax,0.0); << 89 vEnergy.resize(fLevelMax,0.0f);
>> 90 vTime.resize(fLevelMax,0.0f);
>> 91 vTimeg.resize(fLevelMax,0.0f);
74 vSpin.resize(fLevelMax,0); 92 vSpin.resize(fLevelMax,0);
75 vLevel.resize(fLevelMax,nullptr); 93 vLevel.resize(fLevelMax,nullptr);
>> 94 vMeta.resize(fLevelMax,0);
>> 95 vIndexDB.resize(fLevelMax,-1);
>> 96 }
>> 97
>> 98 const G4LevelManager*
>> 99 G4LevelReader::CreateLevelManager(G4int Z, G4int A)
>> 100 {
>> 101 std::ostringstream ss;
>> 102 ss << "/z" << Z << ".a" << A;
>> 103 G4String st = G4String(ss.str());
>> 104 fFile = fDirectory + st;
>> 105 return MakeLevelManager(Z, A, fFile);
>> 106 }
>> 107
>> 108 const G4LevelManager*
>> 109 G4LevelReader::MakeLevelManager(G4int Z, G4int A, const G4String& filename)
>> 110 {
>> 111 vEnergy.resize(1,0.0f);
>> 112 vTime.resize(1,FLT_MAX);
>> 113 vTimeg.resize(1,FLT_MAX);
>> 114 vSpin.resize(1,0);
>> 115 vMeta.resize(1,0);
>> 116 vLevel.resize(1,nullptr);
>> 117
>> 118 std::ifstream infile(filename, std::ios::in);
>> 119
>> 120 // file is not opened
>> 121 if (!infile.is_open()) {
>> 122 if (fVerbose > 0) {
>> 123 G4cout << " G4LevelReader: fail open file for Z= "
>> 124 << Z << " A= " << A
>> 125 << " <" << filename << ">" << G4endl;
>> 126 }
>> 127
>> 128 } else {
>> 129
>> 130 if (fVerbose > 0) {
>> 131 G4cout << "G4LevelReader: open file for Z= "
>> 132 << Z << " A= " << A
>> 133 << " <" << filename << ">" << G4endl;
>> 134 }
>> 135 // read line by line
>> 136 G4bool end = false;
>> 137 G4bool next = true;
>> 138 G4int nline = -1;
>> 139 G4String xl = "- ";
>> 140 fCurrEnergy = DBL_MAX;
>> 141 do {
>> 142 fPol = " ";
>> 143 ++nline;
>> 144 G4bool isOK = (ReadDataItem(infile,fEnergy) &&
>> 145 ReadDataItem(infile,fTrEnergy) &&
>> 146 ReadDataItem(infile,fProb) &&
>> 147 ReadDataItem(infile,fPol) &&
>> 148 ReadDataItem(infile,fTime) &&
>> 149 ReadDataItem(infile,fSpin) &&
>> 150 ReadDataItem(infile,fAlpha));
>> 151
>> 152 fEnergy *= CLHEP::keV;
>> 153 fTrEnergy *= CLHEP::keV;
>> 154
>> 155 if(isOK) {
>> 156 for(G4int i=0; i<10; ++i) {
>> 157 isOK = (isOK && (ReadDataItem(infile,fICC[i])));
>> 158 }
>> 159 }
>> 160 if(!isOK) {
>> 161 end = true;
>> 162 next = false;
>> 163 }
>> 164 if(fVerbose > 1) {
>> 165 G4cout << "#Line " << nline << " " << fEnergy << " " << fTrEnergy
>> 166 << " " << fProb << " " << fPol << " " << fTime << " "
>> 167 << fSpin << " " << fAlpha << G4endl;
>> 168 G4cout << " ";
>> 169 for(G4int i=0; i<10; ++i) { G4cout << fICC[i] << " "; }
>> 170 G4cout << G4endl;
>> 171 }
>> 172 // end of nuclear level data
>> 173 if(end || fEnergy > fCurrEnergy) {
>> 174 size_t nn = vTrans.size();
>> 175 if(nn > 0) {
>> 176 --nn;
>> 177 if(fVerbose > 1) {
>> 178 G4cout << "Reader: new level E= " << fEnergy
>> 179 << " Ntransitions= " << nn+1 << " fNorm1= " << fNorm1
>> 180 << " fNorm2= " << fNorm2 << G4endl;
>> 181 }
>> 182 if(fNorm1 > 0.0f) {
>> 183 fNorm1 = 1.0f/fNorm1;
>> 184 vTimeg.push_back(((G4float)(fTime*fTimeFactor))*fNorm2*fNorm1);
>> 185 fNorm2 = 1.0f/fNorm2;
>> 186 for(size_t i=0; i<nn; ++i) {
>> 187 vGammaCumProbability[i] *= fNorm1;
>> 188 vGammaECumProbability[i] *= fNorm2;
>> 189 if(fVerbose > 2) {
>> 190 G4cout << "Probabilities[" << i
>> 191 << "]= " << vGammaCumProbability[i]
>> 192 << " " << vGammaECumProbability[i]
>> 193 << " idxTrans= " << vIndex[i]
>> 194 << G4endl;
>> 195 }
>> 196 }
>> 197 vGammaCumProbability[nn] = 1.0f;
>> 198 vGammaECumProbability[nn] = 1.0f;
>> 199 if(fVerbose > 2) {
>> 200 G4cout << "Probabilities[" << nn << "]= " << vGammaCumProbability[nn]
>> 201 << " " << vGammaECumProbability[nn]
>> 202 << " IdxTrans= " << vIndex[nn]
>> 203 << G4endl;
>> 204 }
>> 205 vMeta.push_back(0);
>> 206 //case of X-level
>> 207 } else {
>> 208 vMeta.push_back(1);
>> 209 vTimeg.push_back(0.0f);
>> 210 vGammaCumProbability[0] = 0.0f;
>> 211 vGammaECumProbability[0] = 0.0f;
>> 212 }
>> 213
>> 214 vLevel.push_back(new G4NucLevel(vIndex.size(),
>> 215 vIndex,
>> 216 vTrans,
>> 217 vGammaCumProbability,
>> 218 vGammaECumProbability,
>> 219 vGammaProbability,
>> 220 vMpRatio,
>> 221 vShellProbability));
>> 222 vIndex.clear();
>> 223 vTrans.clear();
>> 224 vGammaCumProbability.clear();
>> 225 vGammaECumProbability.clear();
>> 226 vGammaProbability.clear();
>> 227 vShellProbability.clear();
>> 228 vMpRatio.clear();
>> 229 }
>> 230 if(!end) { next = true; }
>> 231 }
>> 232 fCurrEnergy = fEnergy;
>> 233 // begin nuclear level data
>> 234 if(next) {
>> 235 if(fVerbose > 2) {
>> 236 G4cout << "== Reader: begin of new level E= " << fEnergy << G4endl;
>> 237 }
>> 238 // protection for bad level energy
>> 239 size_t nn = vEnergy.size();
>> 240 G4float ener = (G4float)fEnergy;
>> 241 if(0 < nn && vEnergy[nn-1] > ener) { ener = vEnergy[nn-1]; }
>> 242 vEnergy.push_back(ener);
>> 243 vTime.push_back((G4float)(fTime*fTimeFactor));
>> 244 if(fSpin > 20.0) { fSpin = 0.0; }
>> 245 fProb = std::max(fProb, fMinProbability);
>> 246 vSpin.push_back((G4int)(fSpin+fSpin));
>> 247 fNorm1 = 0.0f;
>> 248 fNorm2 = 0.0f;
>> 249 next = false;
>> 250 }
>> 251 // continue filling level data
>> 252 if(!end) {
>> 253 if(fProb > 0.0) {
>> 254 // by default transition to a ground state
>> 255 G4float efinal = std::max((G4float)(fEnergy - fTrEnergy),0.0f);
>> 256 G4float elevel = 0.0f;
>> 257 size_t idxLevel = 0;
>> 258 G4int tnum = 0;
>> 259 // do not check initial energy
>> 260 size_t nn = vEnergy.size();
>> 261 static const G4float x_energy = (G4float)(0.1*CLHEP::eV);
>> 262 if(1 < nn) {
>> 263 G4float ediffMin = fEnergy;
>> 264 for(size_t i=0; i<nn-1; ++i) {
>> 265 G4float ediff = std::abs(efinal - vEnergy[i]);
>> 266 /*
>> 267 G4cout << "Elevel[" << i << "]= " << vEnergy[i]
>> 268 << " Efinal= " << efinal
>> 269 << " Ediff= " << ediff
>> 270 << " EdiffMin= " << ediffMin << G4endl;
>> 271 */
>> 272 if(ediff < ediffMin) {
>> 273 ediffMin = ediff;
>> 274 elevel = vEnergy[i];
>> 275 idxLevel = i;
>> 276 if(ediff <= x_energy) { break; }
>> 277 }
>> 278 }
>> 279 if(std::abs(vEnergy[nn-1] - elevel) < x_energy) { tnum = 1; }
>> 280 }
>> 281 G4double x = 1.0 + fAlpha;
>> 282 fNorm1 += (G4float)fProb;
>> 283 fNorm2 += (G4float)(fProb*x);
>> 284 vIndex.push_back(idxLevel);
>> 285 vGammaCumProbability.push_back(fNorm1);
>> 286 vGammaECumProbability.push_back(fNorm2);
>> 287 vGammaProbability.push_back((G4float)(1.0/x));
>> 288 vMpRatio.push_back(0.0f);
>> 289 vTrans.push_back(tnum);
>> 290 if(fAlpha > 0.0) {
>> 291 vShellProbability.push_back(NormalizedICCProbability(Z));
>> 292 } else {
>> 293 vShellProbability.push_back(nullptr);
>> 294 }
>> 295 }
>> 296 }
>> 297 if(nline > 10000) {
>> 298 G4cout << "G4LevelReader: Line #" << nline << " Z= " << Z << " A= "
>> 299 << " this file is too long - stop loading" << G4endl;
>> 300 end = true;
>> 301 }
>> 302 // Loop checking, 05-Aug-2015, Vladimir Ivanchenko
>> 303 } while (!end);
>> 304 infile.close();
>> 305 }
>> 306
>> 307 G4LevelManager* man = nullptr;
>> 308 if(vEnergy.size() >= 2) {
>> 309 man = new G4LevelManager(vEnergy.size(),vEnergy,vTime,vTimeg,vSpin,vMeta,vLevel);
>> 310 if(fVerbose > 0) {
>> 311 G4cout << "=== Reader: new manager for Z= " << Z << " A= " << A
>> 312 << " Nlevels= " << vEnergy.size() << " E[0]= "
>> 313 << vEnergy[0]/CLHEP::MeV << " MeV Emax= "
>> 314 << man->MaxLevelEnergy()/CLHEP::MeV << " MeV "
>> 315 << " S: " << vEnergy.size() << " " << vTime.size()
>> 316 << " " << vTimeg.size() << " " << vSpin.size() << " " << vLevel.size()
>> 317 << G4endl;
>> 318 }
>> 319 }
>> 320 return man;
76 } 321 }
77 322
78 G4bool G4LevelReader::ReadData(std::istringstr 323 G4bool G4LevelReader::ReadData(std::istringstream& stream, G4double& x)
79 { 324 {
80 stream >> x; 325 stream >> x;
81 return !stream.fail(); << 326 return stream.fail() ? false : true;
82 } 327 }
83 328
84 G4bool G4LevelReader::ReadDataItem(std::istrea 329 G4bool G4LevelReader::ReadDataItem(std::istream& dataFile, G4double& x)
85 { 330 {
86 x = 0.0; 331 x = 0.0;
87 for(G4int i=0; i<nbufmax; ++i) { buffer[i] = 332 for(G4int i=0; i<nbufmax; ++i) { buffer[i] = ' '; }
88 G4bool okay = true; 333 G4bool okay = true;
89 dataFile >> buffer; 334 dataFile >> buffer;
90 if(dataFile.fail()) { okay = false; } 335 if(dataFile.fail()) { okay = false; }
91 else { x = std::strtod(buffer, 0); } << 336 else { x = strtod(buffer, 0); }
>> 337
92 return okay; 338 return okay;
93 } 339 }
94 340
95 G4bool G4LevelReader::ReadDataItem(std::istrea << 341 G4bool G4LevelReader::ReadDataItem(std::istream& dataFile, G4int& ix)
96 { 342 {
97 x = 0.0f; << 343 ix = 0;
98 for(G4int i=0; i<nbuf1; ++i) { buff1[i] = ' << 344 for(G4int i=0; i<nbuf2; ++i) { buff2[i] = ' '; }
99 G4bool okay = true; 345 G4bool okay = true;
100 dataFile >> buff1; << 346 dataFile >> buff2;
101 if(dataFile.fail()) { okay = false; } 347 if(dataFile.fail()) { okay = false; }
102 else { x = std::atof(buff1); } << 348 else { ix = atoi(buff2); }
103 349
104 return okay; 350 return okay;
105 } 351 }
106 352
107 G4bool G4LevelReader::ReadDataItem(std::istrea << 353 G4bool G4LevelReader::ReadDataItem(std::istream& dataFile, G4String& x)
108 { 354 {
109 ix = 0; <<
110 for(G4int i=0; i<nbuf2; ++i) { buff2[i] = ' <<
111 G4bool okay = true; 355 G4bool okay = true;
112 dataFile >> buff2; << 356 dataFile >> bufp;
113 if(dataFile.fail()) { okay = false; } 357 if(dataFile.fail()) { okay = false; }
114 else { ix = std::atoi(buff2); } << 358 else { x = G4String(bufp, 2); }
115 <<
116 return okay; 359 return okay;
117 } 360 }
118 361
119 const std::vector<G4float>* G4LevelReader::Nor 362 const std::vector<G4float>* G4LevelReader::NormalizedICCProbability(G4int Z)
120 { 363 {
121 std::vector<G4float>* vec = nullptr; 364 std::vector<G4float>* vec = nullptr;
122 G4int LL = 3; 365 G4int LL = 3;
123 G4int M = 5; 366 G4int M = 5;
124 G4int N = 1; 367 G4int N = 1;
125 if(Z <= 27) { 368 if(Z <= 27) {
126 M = N = 0; 369 M = N = 0;
127 if(Z <= 4) { 370 if(Z <= 4) {
128 LL = 1; 371 LL = 1;
129 } else if(Z <= 6) { 372 } else if(Z <= 6) {
130 LL = 2; 373 LL = 2;
131 } else if(Z <= 10) { 374 } else if(Z <= 10) {
132 } else if(Z <= 12) { 375 } else if(Z <= 12) {
133 M = 1; 376 M = 1;
134 } else if(Z <= 17) { 377 } else if(Z <= 17) {
135 M = 2; 378 M = 2;
136 } else if(Z == 18) { 379 } else if(Z == 18) {
137 M = 3; 380 M = 3;
138 } else if(Z <= 20) { 381 } else if(Z <= 20) {
139 M = 3; 382 M = 3;
140 N = 1; 383 N = 1;
141 } else { 384 } else {
142 M = 4; 385 M = 4;
143 N = 1; 386 N = 1;
144 } 387 }
145 if(LL < 3) { for(G4int i=LL+1; i<=4; ++i) << 388 if(LL < 3) { for(G4int i=LL+1; i<=4; ++i) { fICC[i] = 0.0; } }
146 if(M < 5) { for(G4int i=M+4; i<=8; ++i) << 389 if(M < 5) { for(G4int i=M+4; i<=8; ++i) { fICC[i] = 0.0; } }
147 if(N < 1) { fICC[9] = 0.0f; } << 390 if(N < 1) { fICC[9] = 0.0; }
148 } 391 }
149 G4float norm = 0.0f; << 392 G4float norm = 0.0;
150 for(G4int i=0; i<10; ++i) { 393 for(G4int i=0; i<10; ++i) {
151 norm += fICC[i]; 394 norm += fICC[i];
152 fICC[i] = norm; 395 fICC[i] = norm;
153 } 396 }
154 if(norm == 0.0f && fAlpha > 0.0f) { <<
155 fICC[0] = norm = 1.0f; <<
156 } <<
157 if(norm > 0.0f) { 397 if(norm > 0.0f) {
158 norm = 1.0f/norm; 398 norm = 1.0f/norm;
159 vec = new std::vector<G4float>; 399 vec = new std::vector<G4float>;
160 G4float x; 400 G4float x;
161 for(G4int i=0; i<10; ++i) { 401 for(G4int i=0; i<10; ++i) {
162 x = fICC[i]*norm; << 402 x = (G4float)(fICC[i]*norm);
163 if(x > 0.995f || 9 == i) { << 403 if(x > 0.995f) {
164 vec->push_back(1.0f); 404 vec->push_back(1.0f);
165 break; 405 break;
166 } 406 }
167 vec->push_back(x); 407 vec->push_back(x);
168 } 408 }
169 if (fVerbose > 3) { << 409 if (fVerbose > 2) {
170 G4long prec = G4cout.precision(3); << 410 G4int prec = G4cout.precision(3);
171 G4cout << "# InternalConv: "; 411 G4cout << "# InternalConv: ";
172 std::size_t nn = vec->size(); << 412 G4int nn = vec->size();
173 for(std::size_t i=0; i<nn; ++i) { G4cout << 413 for(G4int i=0; i<nn; ++i) { G4cout << " " << (*vec)[i]; }
174 G4cout << G4endl; 414 G4cout << G4endl;
175 G4cout.precision(prec); 415 G4cout.precision(prec);
176 } 416 }
177 } 417 }
178 return vec; 418 return vec;
179 } 419 }
180 420
181 const G4LevelManager* 421 const G4LevelManager*
182 G4LevelReader::CreateLevelManager(G4int Z, G4i << 422 G4LevelReader::CreateLevelManagerNEW(G4int Z, G4int A)
183 { 423 {
184 std::ostringstream ss; 424 std::ostringstream ss;
185 ss << fDirectory << "/z" << Z << ".a" << A; << 425 ss << "/correlated_gamma/z" << Z << ".a" << A;
186 std::ifstream infile(ss.str(), std::ios::in) << 426 G4String st = G4String(ss.str());
>> 427 fFile = fDirectory + st;
>> 428 std::ifstream infile(fFile, std::ios::in);
187 429
188 // file is not opened 430 // file is not opened
189 if (!infile.is_open()) { 431 if (!infile.is_open()) {
190 if(fVerbose > 1) { << 432 if (fVerbose > 0) {
191 G4ExceptionDescription ed; << 433 G4cout << " G4LevelReader: fail open file for Z= "
192 ed << "Regular file " << ss.str() << " i << 434 << Z << " A= " << A
193 << Z << " A=" << A; << 435 << " <" << fFile << ">" << G4endl;
194 G4Exception("G4LevelReader::LevelManager <<
195 JustWarning, ed, "Check file path"); <<
196 } 436 }
197 return nullptr; 437 return nullptr;
198 } 438 }
199 // file is opened << 439 if (fVerbose > 0) {
200 if (fVerbose > 1) { << 440 G4cout << "G4LevelReader: open file for Z= "
201 G4cout << "G4LevelReader: open file " << s << 441 << Z << " A= " << A
202 << Z << " A= " << A << G4endl; << 442 << " <" << fFile << ">" << G4endl;
203 } 443 }
204 return LevelManager(Z, A, infile); << 444 return LevelManager(Z, A, 0, infile);
205 } 445 }
206 446
207 const G4LevelManager* 447 const G4LevelManager*
208 G4LevelReader::MakeLevelManager(G4int Z, G4int << 448 G4LevelReader::MakeLevelManagerNEW(G4int Z, G4int A,
>> 449 const G4String& filename)
209 { 450 {
210 std::ifstream infile(filename, std::ios::in) 451 std::ifstream infile(filename, std::ios::in);
211 << 452
212 // file is not opened 453 // file is not opened
213 if (!infile.is_open()) { 454 if (!infile.is_open()) {
214 if(fVerbose > 1) { << 455 if (fVerbose > 0) {
215 G4ExceptionDescription ed; << 456 G4cout << " G4LevelReader: fail open file for Z= "
216 ed << "External file " << filename << " << 457 << Z << " A= " << A
217 << Z << " A=" << A; << 458 << " <" << filename << ">" << G4endl;
218 G4Exception("G4LevelReader::LevelManager <<
219 FatalException, ed, "Check file path"); <<
220 } 459 }
221 return nullptr; 460 return nullptr;
222 } 461 }
223 // file is opened << 462 if (fVerbose > 0) {
224 if (fVerbose > 1) { << 463 G4cout << "G4LevelReader: open file for Z= "
225 G4cout << "G4LevelReader: open external fi << 464 << Z << " A= " << A
226 << " for Z= " << Z << " A= " << A < << 465 << " <" << filename << ">" << G4endl;
227 } 466 }
228 return LevelManager(Z, A, infile); << 467 return LevelManager(Z, A, 0, infile);
229 } 468 }
230 469
231 const G4LevelManager* 470 const G4LevelManager*
232 G4LevelReader::LevelManager(G4int Z, G4int A, << 471 G4LevelReader::LevelManager(G4int Z, G4int A, G4int nlev,
>> 472 std::ifstream& infile)
233 { 473 {
234 G4bool allLevels = fData->GetParameters()->S << 474 G4bool allLevels = fParam->StoreAllLevels();
235 fPol = " "; << 475 G4float emax = fData->GetMaxLevelEnergy(Z, A);
236 G4int i = 0; << 476
237 for (;;) { << 477 static const G4double fkev = CLHEP::keV;
238 infile >> i1 >> fPol; // Level number a << 478 G4int nlevels = (0 == nlev) ? fLevelMax : nlev;
239 // normal end of file << 479 if(fVerbose > 0) {
240 if (infile.eof()) { << 480 G4cout << "## New isotope Z= " << Z << " A= " << A;
241 if (fVerbose > 1) { << 481 if(nlevels < fLevelMax) { G4cout << " Nlevels= " << nlevels; }
>> 482 G4cout << G4endl;
>> 483 }
>> 484 if(nlevels > fLevelMax) {
>> 485 fLevelMax = nlevels;
>> 486 vEnergy.resize(fLevelMax,0.0f);
>> 487 vTime.resize(fLevelMax,0.0f);
>> 488 vTimeg.resize(fLevelMax,0.0f);
>> 489 vSpin.resize(fLevelMax,0);
>> 490 vLevel.resize(fLevelMax,0);
>> 491 vMeta.resize(fLevelMax,0);
>> 492 vIndexDB.resize(fLevelMax,-1);
>> 493 }
>> 494 G4int ntrans(0), i(0), i1, i2, i3, j, k;
>> 495 G4String xf(" ");
>> 496 G4float x, x1;
>> 497
>> 498 for(G4int ii=0; ii<nlevels; ++ii) {
>> 499
>> 500 infile >> i1 >> xf;
>> 501 if(infile.eof()) {
>> 502 if(fVerbose > 1) {
242 G4cout << "### End of file Z= " << Z << " A= 503 G4cout << "### End of file Z= " << Z << " A= " << A
243 << " Nlevels= " << i << G4endl; << 504 << " Nlevels= " << ii << G4endl;
244 } 505 }
245 break; 506 break;
246 } 507 }
247 // start reading new level data << 508 if(!(ReadDataItem(infile,fEnergy) &&
248 #ifdef G4VERBOSE << 509 ReadDataItem(infile,fTime) &&
249 if(fVerbose > 2) { << 510 ReadDataItem(infile,fSpin) &&
250 G4cout << "New line: i1= " << i1 << " f << 511 ReadDataItem(infile,ntrans))) {
251 } << 512 if(fVerbose > 1) {
252 #endif << 513 G4cout << "### End of file Z= " << Z << " A= " << A
253 // read new level data << 514 << " Nlevels= " << ii << G4endl;
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 } 515 }
262 break; 516 break;
263 } 517 }
264 fEnergy *= CLHEP::keV; << 518 fTime = std::max(fTime, 0.0);
265 for (k=0; k<nfloting; ++k) { << 519 fEnergy *= fkev;
266 if (fPol == fFloatingLevels[k]) { << 520 for(k=0; k<nfloting; ++k) {
>> 521 if(xf == fFloatingLevels[k]) {
267 break; 522 break;
268 } 523 }
269 } 524 }
270 // if a previous level has higher energy t << 525
271 // data with wrong level should red anyway << 526 // if a previous level has not transitions it may be ignored
272 if (0 < i) { << 527 if(0 < ii) {
273 if (fEnergy < vEnergy[i-1]) { << 528 // do not store level without transitions
274 #ifdef G4VERBOSE << 529 if(!allLevels && 0 == k && 0 == ntrans) { continue; }
275 ++count1; << 530
276 if (count1 < countmax && fVerbose > 0) { << 531 // protection
277 G4cout << "### G4LevelReader: broken level << 532 if(fEnergy < vEnergy[i-1]) {
278 << " E(MeV)= " << fEnergy << " < " << vEn << 533 G4cout << "### G4LevelReader: broken level " << ii
279 << " for isotope Z= " << Z << " A= " << 534 << " E(MeV)= " << fEnergy << " < " << vEnergy[i-1]
280 << A << " level energy increased" << G4en << 535 << " for isotope Z= " << Z << " A= "
281 } << 536 << A << " level energy increased" << G4endl;
282 #endif << 537 fEnergy = vEnergy[i-1];
283 // for any case <<
284 fEnergy = vEnergy[i-1] + eTolarence; <<
285 } 538 }
>> 539 // upper limit
>> 540 if(fEnergy > emax) { break; }
286 } 541 }
287 vEnergy[i] = fEnergy; << 542 vEnergy[i] = (G4float)fEnergy;
288 if (fTime > 0.0) { fTime *= fTimeFactor; << 543 vTime[i] = (G4float)(fTime*fTimeFactor);
289 else if (fTime < 0.0) { fTime = DBL_MAX; } << 544 vTimeg[i] = vTime[i];
290 << 545 if(fSpin > 20.0) { fSpin = 0.0; }
291 G4int twos = G4lrint(fSpin + fSpin); << 546 vSpin[i] = (G4int)(fSpin + fSpin);
292 twos = std::max(twos, -100); << 547 vMeta[i] = k;
293 vSpin[i] = 100 + twos + k*100000; << 548 vIndexDB[ii] = i;
294 #ifdef G4VERBOSE << 549 if(fVerbose > 1) {
295 if (fVerbose > 2) { << 550 G4cout << " Level #" << i1 << " E(MeV)= " << fEnergy/CLHEP::MeV
296 G4cout << " Level #" << i1 << " E(MeV) << 551 << " LTime(s)= " << fTime << " 2S= " << vSpin[i]
297 << " LTime(s)=" << fTime << " 2S=" << << 552 << " meta= " << vMeta[i] << " idx= " << i << " ii= " << ii
298 << " meta=" << vSpin[i]/100000 << " id << 553 << " ntr= " << ntrans << G4endl;
299 << " ntr=" << ntrans << G4endl; <<
300 } 554 }
301 #endif <<
302 vLevel[i] = nullptr; 555 vLevel[i] = nullptr;
303 if (ntrans == 0) { << 556 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 557
311 G4bool isTransOK = true; << 558 // there are transitions
312 if (ntrans > fTransMax) { << 559 if(ntrans > fTransMax) {
313 fTransMax = ntrans; 560 fTransMax = ntrans;
>> 561 vIndex.resize(fTransMax);
314 vTrans.resize(fTransMax); 562 vTrans.resize(fTransMax);
315 vRatio.resize(fTransMax); 563 vRatio.resize(fTransMax);
316 vGammaCumProbability.resize(fTransMax); 564 vGammaCumProbability.resize(fTransMax);
>> 565 vGammaECumProbability.resize(fTransMax);
317 vGammaProbability.resize(fTransMax); 566 vGammaProbability.resize(fTransMax);
318 vShellProbability.resize(fTransMax); 567 vShellProbability.resize(fTransMax);
>> 568 vMpRatio.resize(fTransMax);
319 } 569 }
320 fNorm1 = 0.0f; << 570 fNorm1 = fNorm2 = 0.0f;
321 for (G4int j=0; j<ntrans; ++j) { << 571 j = 0;
>> 572 for(G4int jj=0; jj<ntrans; ++jj) {
322 573
323 if (!(ReadDataItem(infile, i2) && << 574 if(!(ReadDataItem(infile,i2) &&
324 ReadDataItem(infile, fTransEnergy) && << 575 ReadDataItem(infile,fTrEnergy) &&
325 ReadDataItem(infile, fProb) && << 576 ReadDataItem(infile,fProb) &&
326 ReadDataItem(infile, tnum) && << 577 ReadDataItem(infile,vTrans[j]) &&
327 ReadDataItem(infile, vRatio[j]) && << 578 ReadDataItem(infile,fRatio) &&
328 ReadDataItem(infile, fAlpha))) { << 579 ReadDataItem(infile,fAlpha))) {
329 #ifdef G4VERBOSE << 580 //infile >>i2 >> fTrEnergy >> fProb >> vTrans[j] >> fRatio >> fAlpha;
330 ++count2; << 581 //if(infile.fail()) {
331 if (count2 < countmax && fVerbose > 0) { << 582 if(fVerbose > 0) {
332 G4cout << "### Fail to read transition j << 583 G4cout << "### Fail to read transition j= " << j
333 << " j=" << j << " i2=" << i2 << 584 << " Z= " << Z << " A= " << A << G4endl;
334 << " Z=" << Z << " A=" << A << G4endl; <<
335 } 585 }
336 #endif << 586 break;
337 isTransOK = false; <<
338 } 587 }
339 if (i2 >= i) { << 588 if(i2 >= ii) {
340 #ifdef G4VERBOSE << 589 G4cout << "### G4LevelReader: broken transition " << j
341 ++count2; << 590 << " from level " << ii << " to " << i2
342 if (count2 < countmax) { << 591 << " for isotope Z= " << Z << " A= "
343 G4cout << "### G4LevelReader: broken tra << 592 << A << " - use ground level" << G4endl;
344 << " from level " << i << " to " << i2 << 593 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 } 594 }
352 vTrans[j] = i2*10000 + tnum; << 595 i3 = vIndexDB[std::abs(i2)];
353 fAlpha = std::min(std::max(fAlpha,0.f) << 596
354 G4float x = 1.0f + fAlpha; << 597 if(i3 >= 0) {
355 fNorm1 += x*fProb; << 598 vIndex[j] = i3;
356 vGammaCumProbability[j] = fNorm1; << 599 x = 1.0f + std::max((G4float)fAlpha,0.0f);
357 vGammaProbability[j] = 1.0f/x; << 600 x1= (G4float)fProb;
358 vShellProbability[j] = nullptr; << 601 fNorm1 += x1;
359 if (fVerbose > 2) { << 602 fNorm2 += x*x1;
360 G4long prec = G4cout.precision(4); << 603 vGammaCumProbability[j] = fNorm1;
361 G4cout << "### Transition #" << j << " to << 604 vGammaECumProbability[j]= fNorm2;
362 << " i2= " << i2 << " Etrans(MeV)= " << << 605 vGammaProbability[j] = 1.0f/x;
363 << " fProb= " << fProb << " MultiP= " << << 606 vRatio[j] = (G4float)fRatio;
364 << " fMpRatio= " << fRatio << " fAlpha= << 607 vShellProbability[j] = nullptr;
365 << G4endl; << 608 if(fVerbose > 1) {
366 G4cout.precision(prec); << 609 fTrEnergy *= fkev;
>> 610 G4int prec = G4cout.precision(4);
>> 611 G4cout << "### Transition #" << j << " to level " << vIndex[j]
>> 612 << " i2= " << i2 << " Etrans(MeV)= " << fTrEnergy
>> 613 << " fProb= " << fProb << " MultiP= " << vTrans[j]
>> 614 << " fMpRatio= " << fRatio << " fAlpha= " << fAlpha
>> 615 << G4endl;
>> 616 G4cout.precision(prec);
>> 617 }
367 } 618 }
368 if (fAlpha > 0.0f) { << 619 if(fAlpha > 0.0) {
369 for (k=0; k<10; ++k) { << 620 for(k=0; k<10; ++k) {
370 if (!ReadDataItem(infile,fICC[k])) { << 621 //infile >> fICC[k];
371 isTransOK = false; << 622 if(!ReadDataItem(infile,fICC[k])) {
372 #ifdef G4VERBOSE << 623 //if(infile.fail()) {
373 ++count2; << 624 if(fVerbose > 0) {
374 if (count2 < countmax) { << 625 G4cout << "### Fail to read convertion coeff k= " << k
375 G4cout << "### G4LevelReader: fail to read <<
376 << " for transition j= " << j 626 << " for transition j= " << j
377 << " Z= " << Z << " A= " << A << G 627 << " Z= " << Z << " A= " << A << G4endl;
378 } 628 }
379 #endif << 629 break;
380 for(kk=k; kk<10; ++kk) { fICC[kk] = 0. <<
381 } 630 }
382 } 631 }
383 if (allLevels) { << 632 if(i3 >= 0) {
384 vShellProbability[j] = NormalizedICCProb 633 vShellProbability[j] = NormalizedICCProbability(Z);
>> 634 if(!vShellProbability[j]) { vGammaProbability[j] = 1.0f; }
385 } 635 }
386 } 636 }
>> 637 if(i3 >= 0) { ++j; }
387 } 638 }
388 if (ntrans > 0) { << 639 if(j > 0) {
389 G4int nt = ntrans - 1; << 640 if(0.0f < fNorm1) {
390 if (fVerbose > 2) { << 641 fNorm1 = 1.0f/fNorm1;
391 G4cout << "=== New G4NucLevel: Ntran << 642 vTimeg[i] *= fNorm2*fNorm1;
392 << " Time(ns)=" << fTime << 643 fNorm2 = 1.0f/fNorm2;
393 << " IdxTrans=" << vTrans[nt]/10000 << 644 }
394 << " isOK=" << isTransOK << 645 G4int nt = j - 1;
395 << G4endl; << 646 for(k=0; k<nt; ++k) {
396 } << 647 vGammaCumProbability[k] *= fNorm1;
397 if (0.0f < fNorm1) { fNorm1 = 1.0f/fNorm1; } << 648 vGammaECumProbability[k] *= fNorm2;
398 for (k=0; k<nt; ++k) { << 649 if(fVerbose > 2) {
399 vGammaCumProbability[k] *= fNorm1; <<
400 #ifdef G4VERBOSE <<
401 if (fVerbose > 3) { <<
402 G4cout << "Probabilities[" << k 650 G4cout << "Probabilities[" << k
403 << "]= " << vGammaCumProbability[k] 651 << "]= " << vGammaCumProbability[k]
404 << " " << vGammaProbability[k] << 652 << " " << vGammaECumProbability[k]
405 << " idxTrans= " << vTrans[k]/10000 << 653 << " idxTrans= " << vIndex[k]
406 << G4endl; 654 << G4endl;
407 } 655 }
408 #endif <<
409 } 656 }
410 vGammaCumProbability[nt] = 1.0f; << 657 vGammaCumProbability[nt] = 1.0f;
411 vLevel[i] = new G4NucLevel((std::size_t)ntra << 658 vGammaECumProbability[nt] = 1.0f;
>> 659 if(fVerbose > 2) {
>> 660 G4cout << "Probabilities[" << nt << "]= "
>> 661 << vGammaCumProbability[nt]
>> 662 << " " << vGammaECumProbability[nt]
>> 663 << " IdxTrans= " << vIndex[nt]
>> 664 << G4endl;
>> 665 }
>> 666 vLevel[i] = new G4NucLevel((size_t)j,
>> 667 vIndex,
>> 668 vTrans,
412 vGammaCumProbability, 669 vGammaCumProbability,
>> 670 vGammaECumProbability,
413 vGammaProbability, 671 vGammaProbability,
414 vRatio, << 672 vMpRatio,
415 vShellProbability); 673 vShellProbability);
416 } 674 }
417 } 675 }
418 ++i; 676 ++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 } 677 }
426 G4LevelManager* lman = nullptr; 678 G4LevelManager* lman = nullptr;
427 if (1 <= i) { << 679 if(1 < i) {
428 lman = new G4LevelManager(Z, A, (std::size << 680 lman = new G4LevelManager((size_t)i,vEnergy,vTime,vTimeg,
429 if (fVerbose > 1) { << 681 vSpin,vMeta,vLevel);
430 G4cout << "=== Reader: new manager for Z << 682 if(fVerbose > 0) {
431 << " Nlevels=" << i << " E[0]=" << 683 G4cout << "=== Reader: new manager for Z= " << Z << " A= " << A
432 << vEnergy[0]/CLHEP::MeV << " MeV E1= << 684 << " Nlevels= " << i << " E[0]= "
433 << vEnergy[i-1]/CLHEP::MeV << " MeV" << 685 << vEnergy[0]/CLHEP::MeV << " MeV Emax= "
434 << " count1,2=" << count1 << ", " << co << 686 << vEnergy[i-1]/CLHEP::MeV << " MeV "
>> 687 << G4endl;
435 } 688 }
>> 689 }
>> 690 for(G4int ii=0; ii<nlevels; ++ii) {
>> 691 vIndexDB[ii] = -1;
436 } 692 }
437 return lman; 693 return lman;
438 } 694 }
439 695