EGSnrc C++ class library  Report PIRS-898 (2021)
Iwan Kawrakow, Ernesto Mainegra-Hing, Frederic Tessier, Reid Townson and Blake Walters
egs_ensdf.cpp
Go to the documentation of this file.
1 /*
2 ###############################################################################
3 #
4 # EGSnrc egs++ ensdf
5 # Copyright (C) 2015 National Research Council Canada
6 #
7 # This file is part of EGSnrc.
8 #
9 # EGSnrc is free software: you can redistribute it and/or modify it under
10 # the terms of the GNU Affero General Public License as published by the
11 # Free Software Foundation, either version 3 of the License, or (at your
12 # option) any later version.
13 #
14 # EGSnrc is distributed in the hope that it will be useful, but WITHOUT ANY
15 # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
16 # FOR A PARTICULAR PURPOSE. See the GNU Affero General Public License for
17 # more details.
18 #
19 # You should have received a copy of the GNU Affero General Public License
20 # along with EGSnrc. If not, see <http://www.gnu.org/licenses/>.
21 #
22 ###############################################################################
23 #
24 # Author: Reid Townson, 2016
25 #
26 # Contributors:
27 #
28 ###############################################################################
29 */
30 
31 
38 #include "egs_ensdf.h"
39 
40 map<string, unsigned short int> getElementMap() {
41  map<string, unsigned short int> elementTable;
42  elementTable["H"] = 1;
43  elementTable["HE"] = 2;
44  elementTable["LI"] = 3;
45  elementTable["BE"] = 4;
46  elementTable["B"] = 5;
47  elementTable["C"] = 6;
48  elementTable["N"] = 7;
49  elementTable["O"] = 8;
50  elementTable["F"] = 9;
51  elementTable["NE"] = 10;
52  elementTable["NA"] = 11;
53  elementTable["MG"] = 12;
54  elementTable["AL"] = 13;
55  elementTable["SI"] = 14;
56  elementTable["P"] = 15;
57  elementTable["S"] = 16;
58  elementTable["CL"] = 17;
59  elementTable["AR"] = 18;
60  elementTable["K"] = 19;
61  elementTable["CA"] = 20;
62  elementTable["SC"] = 21;
63  elementTable["TI"] = 22;
64  elementTable["V"] = 23;
65  elementTable["CR"] = 24;
66  elementTable["MN"] = 25;
67  elementTable["FE"] = 26;
68  elementTable["CO"] = 27;
69  elementTable["NI"] = 28;
70  elementTable["CU"] = 29;
71  elementTable["ZN"] = 30;
72  elementTable["GA"] = 31;
73  elementTable["GE"] = 32;
74  elementTable["AS"] = 33;
75  elementTable["SE"] = 34;
76  elementTable["BR"] = 35;
77  elementTable["KR"] = 36;
78  elementTable["RB"] = 37;
79  elementTable["SR"] = 38;
80  elementTable["Y"] = 39;
81  elementTable["ZR"] = 40;
82  elementTable["NB"] = 41;
83  elementTable["MO"] = 42;
84  elementTable["TC"] = 43;
85  elementTable["RU"] = 44;
86  elementTable["RH"] = 45;
87  elementTable["PD"] = 46;
88  elementTable["AG"] = 47;
89  elementTable["CD"] = 48;
90  elementTable["IN"] = 49;
91  elementTable["SN"] = 50;
92  elementTable["SB"] = 51;
93  elementTable["TE"] = 52;
94  elementTable["I"] = 53;
95  elementTable["XE"] = 54;
96  elementTable["CS"] = 55;
97  elementTable["BA"] = 56;
98  elementTable["LA"] = 57;
99  elementTable["CE"] = 58;
100  elementTable["PR"] = 59;
101  elementTable["ND"] = 60;
102  elementTable["PM"] = 61;
103  elementTable["SM"] = 62;
104  elementTable["EU"] = 63;
105  elementTable["GD"] = 64;
106  elementTable["TB"] = 65;
107  elementTable["DY"] = 66;
108  elementTable["HO"] = 67;
109  elementTable["ER"] = 68;
110  elementTable["TM"] = 69;
111  elementTable["YB"] = 70;
112  elementTable["LU"] = 71;
113  elementTable["HF"] = 72;
114  elementTable["TA"] = 73;
115  elementTable["W"] = 74;
116  elementTable["RE"] = 75;
117  elementTable["OS"] = 76;
118  elementTable["IR"] = 77;
119  elementTable["PT"] = 78;
120  elementTable["AU"] = 79;
121  elementTable["HG"] = 80;
122  elementTable["TL"] = 81;
123  elementTable["PB"] = 82;
124  elementTable["BI"] = 83;
125  elementTable["PO"] = 84;
126  elementTable["AT"] = 85;
127  elementTable["RN"] = 86;
128  elementTable["FR"] = 87;
129  elementTable["RA"] = 88;
130  elementTable["AC"] = 89;
131  elementTable["TH"] = 90;
132  elementTable["PA"] = 91;
133  elementTable["U"] = 92;
134  elementTable["NP"] = 93;
135  elementTable["PU"] = 94;
136  elementTable["AM"] = 95;
137  elementTable["CM"] = 96;
138  elementTable["BK"] = 97;
139  elementTable["CF"] = 98;
140  elementTable["ES"] = 99;
141  elementTable["FM"] = 100;
142  elementTable["MD"] = 101;
143  elementTable["NO"] = 102;
144  elementTable["LR"] = 103;
145  elementTable["RF"] = 104;
146  elementTable["DB"] = 105;
147  elementTable["SG"] = 106;
148  elementTable["BH"] = 107;
149  elementTable["HS"] = 108;
150  elementTable["MT"] = 109;
151  elementTable["DS"] = 110;
152  elementTable["RG"] = 111;
153  elementTable["CN"] = 112;
154  elementTable["UUT"] = 113;
155  elementTable["UUQ"] = 114;
156  elementTable["UUP"] = 115;
157  elementTable["UUH"] = 116;
158  elementTable["UUS"] = 117;
159  elementTable["UUO"] = 118;
160 
161  return elementTable;
162 }
163 
164 unsigned short int findZ(string element) {
165 
166  transform(element.begin(), element.end(), element.begin(), ::toupper);
167 
168  map<string, unsigned short int> elementMap = getElementMap();
169 
170  if (elementMap.find(element) != elementMap.end()) {
171  return elementMap[element];
172  }
173  else {
174  return 0;
175  }
176 }
177 
178 unsigned short int setZ(string id) {
179 
180  string element;
181  for (unsigned int i=0; i < id.length(); ++i) {
182  if (!isdigit(id[i])) {
183  element.push_back(id[i]);
184  }
185  }
186 
187  unsigned short int Z = findZ(element);
188  if (Z == 0) {
189  egsWarning("setZ: Warning: Element does not exist "
190  "in our data (%s)\n", element.c_str());
191  }
192 
193  return Z;
194 }
195 
196 EGS_Ensdf::EGS_Ensdf(const string nuclide, const string ensdf_filename, const string relaxType, const bool allowMultiTrans, int verbosity) {
197 
198  verbose = verbosity;
199  relaxationType = relaxType;
200  allowMultiTransition = allowMultiTrans;
201 
202  if (ensdf_file.is_open()) {
203  ensdf_file.close();
204  }
205 
206  radionuclide = nuclide.substr(0, nuclide.find_last_of("."));
207 
208  // The parent element
209  string element = radionuclide.substr(0, radionuclide.find("-"));
210  Z = setZ(element);
211 
212  egsInformation("EGS_Ensdf::EGS_Ensdf: Nuclide: "
213  "%s\n",nuclide.c_str());
214  egsInformation("EGS_Ensdf::EGS_Ensdf: Now loading ensdf file: "
215  "\"%s\"\n",ensdf_filename.c_str());
216 
217  ensdf_file.open(ensdf_filename.c_str(),ios::in);
218  if (!ensdf_file.is_open()) {
219  egsWarning("\nEGS_Ensdf::EGS_Ensdf: failed to open ensdf file %s"
220  " for reading\n\n",ensdf_filename.c_str());
221  return;
222  }
223 
224  string line;
225  vector<string> ensdf;
226  while (getline(ensdf_file, line)) {
227  ensdf.push_back(line);
228  }
229 
230  if (ensdf_file.is_open()) {
231  ensdf_file.close();
232  }
233 
234  // Parse the ensdf data
235  parseEnsdf(ensdf);
236 }
237 
238 EGS_Ensdf::~EGS_Ensdf() {
239  if (ensdf_file.is_open()) {
240  ensdf_file.close();
241  }
242 
243  for (vector<ParentRecord * >::iterator it = myParentRecords.begin();
244  it!=myParentRecords.end(); it++) {
245  delete *it;
246  *it=0;
247  }
248  myParentRecords.clear();
249  for (vector<NormalizationRecord * >::iterator it =
250  myNormalizationRecords.begin();
251  it!=myNormalizationRecords.end(); it++) {
252  delete *it;
253  *it=0;
254  }
255  myNormalizationRecords.clear();
256  for (vector<LevelRecord * >::iterator it =
257  myLevelRecords.begin();
258  it!=myLevelRecords.end(); it++) {
259  delete *it;
260  *it=0;
261  }
262  myLevelRecords.clear();
263  for (vector<BetaMinusRecord * >::iterator it =
264  myBetaMinusRecords.begin();
265  it!=myBetaMinusRecords.end(); it++) {
266  delete *it;
267  *it=0;
268  }
269  myBetaMinusRecords.clear();
270  for (vector<BetaPlusRecord * >::iterator it =
271  myBetaPlusRecords.begin();
272  it!=myBetaPlusRecords.end(); it++) {
273  delete *it;
274  *it=0;
275  }
276  myBetaPlusRecords.clear();
277  for (vector<GammaRecord * >::iterator it =
278  myGammaRecords.begin();
279  it!=myGammaRecords.end(); it++) {
280  delete *it;
281  *it=0;
282  }
283  myGammaRecords.clear();
284  for (vector<AlphaRecord * >::iterator it =
285  myAlphaRecords.begin();
286  it!=myAlphaRecords.end(); it++) {
287  delete *it;
288  *it=0;
289  }
290  myAlphaRecords.clear();
291  for (vector<GammaRecord * >::iterator it =
292  myMetastableGammaRecords.begin();
293  it!=myMetastableGammaRecords.end(); it++) {
294  delete *it;
295  *it=0;
296  }
297  myMetastableGammaRecords.clear();
298  for (vector<GammaRecord * >::iterator it =
299  myUncorrelatedGammaRecords.begin();
300  it!=myUncorrelatedGammaRecords.end(); it++) {
301  delete *it;
302  *it=0;
303  }
304  myUncorrelatedGammaRecords.clear();
305 }
306 
307 string egsRemoveWhite(string myString) {
308  string result = "";
309 
310  for (unsigned int i = 0; i<myString.size(); i++) {
311  if (!(myString[i]==' ' || myString[i]=='\n' || myString[i]=='\t')) {
312  result += myString[i];
313  }
314  }
315 
316  return result;
317 }
318 
319 string egsTrimString(string myString) {
320  int start = -1;
321  int end = myString.size();
322  while (myString[++start]==' ');
323  while (myString[--end]==' ');
324  return myString.substr(start,end-start+1);
325 }
326 
327 // Parse an ensdf file to create a decay structure
328 void EGS_Ensdf::parseEnsdf(vector<string> ensdf) {
329  /* IDs of recordStack
330  * 0 Identification (not used)
331  * 1 History (not used)
332  * 2 Q-value(not used)
333  * 3 Cross-Reference (not used)
334  * 4 Comment
335  * 5 Parent
336  * 6 Normalization
337  * 7 Level
338  * 8 Beta-
339  * 9 EC + Beta+
340  * 10 Alpha
341  * 11 Delayed Particle (not used)
342  * 12 Gamma
343  * 13 Reference (not used)
344  * */
345  for (int i = 0; i < 14; i++) {
346  recordStack.push_back(vector<string>());
347  }
348 
349  // Loop over each line
350  // When we recognize a line as containing an important record,
351  // add it to the recordStack
352  // Any time we get to a new record (line[5]==' '), call buildRecords()
353  for (vector<string>::iterator it = ensdf.begin(); it!=ensdf.end(); it++) {
354 
355  string line = *it;
356 
357  if (verbose) {
358  egsInformation("EGS_Ensdf::parseEnsdf: %s\n", line.c_str());
359  }
360 
361  if (line[6]==' ' && line[7]==' ' && line[8]==' ') {
362  // Identification
363 
364  }
365  else if (line[6]==' ' && line[7]=='H' && line[8]==' ') {
366  // History
367 
368  }
369  else if (line[6]== ' ' && line[7]=='Q' && line[8]==' ') {
370  // Q-value
371 
372  }
373  else if (line[6]==' ' && line[7]=='X') {
374  // Cross-Reference
375 
376  }
377  else if ((line[6]=='C' || line[6]=='D' || line[6]=='T' ||
378  line[6]=='c' || line[6]=='d' || line[6]=='t')) {
379 
380  if (line[7]=='G') {
381  // If this is related to a gamma record, keep it
382  recordStack[12].push_back(line);
383  }
384  else {
385  // General comment
386  recordStack[4].push_back(line);
387  }
388 
389  }
390  else if (line[6]==' ' && line[7]=='P') {
391  //Parent
392  if (line[5]==' ') {
393  buildRecords();
394  }
395  recordStack[5].push_back(line);
396 
397  }
398  else if (line[6]==' ' && line[7]=='N') {
399  // Normalization
400  if (line[5]==' ') {
401  buildRecords();
402  }
403  recordStack[6].push_back(line);
404 
405  }
406  if (line[6]==' ' && line[7]=='L' && line[8]== ' ') {
407  // Level
408  if (line[5]==' ') {
409  buildRecords();
410  }
411  recordStack[7].push_back(line);
412 
413  }
414  else if (line[6]==' ' && line[7]=='B' && line[8]==' ') {
415  // Beta-
416  if (line[5]==' ') {
417  buildRecords();
418  }
419  recordStack[8].push_back(line);
420 
421  }
422  else if (line[6]==' ' && line[7]=='E' && line[8]==' ') {
423  // Beta+ and Electron Capture
424  if (line[5]==' ') {
425  buildRecords();
426  }
427  recordStack[9].push_back(line);
428 
429  }
430  else if (line[6]==' ' && line[7]=='A' && line[8]==' ') {
431  // Alpha
432  if (line[5]==' ') {
433  buildRecords();
434  }
435  recordStack[10].push_back(line);
436 
437  }
438  else if (line[6]==' ' && (line[7]=='D' || line[7]==' ') &&
439  (line[8]=='N' || line[8]=='P' || line[8]=='A')) {
440  // Delayed Particle
441  if (line[5]==' ') {
442  buildRecords();
443  }
444  recordStack[11].push_back(line);
445 
446  }
447  else if (line[6]==' ' && line[7]=='G' && line[8]==' ') {
448  // Gamma
449  if (line[5]==' ') {
450  buildRecords();
451  }
452  recordStack[12].push_back(line);
453  }
454  }
455 
456  // Build the records into objects
457  if (!recordStack.empty()) {
458  buildRecords();
459  }
460 
461  // Get X-ray and auger emissions from comments
462  if (relaxationType == "ensdf") {
463  if (verbose) {
464  egsInformation("EGS_Ensdf::parseEnsdf: Checking for x-rays and Auger...\n");
465  }
466 
467  getEmissionsFromComments();
468 
469  if (verbose > 1) {
470  egsInformation("EGS_Ensdf::parseEnsdf: Done checking for x-rays and Auger.\n");
471  }
472  }
473 
474  // Get rid of very low emission probability particles
475  double minimumIntensity = 1e-10;
476  for (vector<BetaMinusRecord * >::iterator it = myBetaMinusRecords.begin();
477  it!=myBetaMinusRecords.end();) {
478  if ((*it)->getBetaIntensity() <= minimumIntensity) {
479  if (verbose) {
480  egsInformation("EGS_Ensdf::parseEnsdf: Removing beta- due to small intensity (%.1e < %.1e)\n",(*it)->getBetaIntensity(),minimumIntensity);
481  }
482  myBetaMinusRecords.erase(it);
483  }
484  else {
485  it++;
486  }
487  }
488  for (vector<BetaPlusRecord * >::iterator it = myBetaPlusRecords.begin();
489  it!=myBetaPlusRecords.end();) {
490  if ((*it)->getBetaIntensity() <= minimumIntensity) {
491  if (verbose) {
492  egsInformation("EGS_Ensdf::parseEnsdf: Removing beta+ due to small intensity (%.1e < %.1e)\n",(*it)->getBetaIntensity(),minimumIntensity);
493  }
494  myBetaPlusRecords.erase(it);
495  }
496  else {
497  it++;
498  }
499  }
500  for (vector<AlphaRecord *>::iterator it = myAlphaRecords.begin();
501  it != myAlphaRecords.end();) {
502  if ((*it)->getAlphaIntensity() <= minimumIntensity) {
503  if (verbose) {
504  egsInformation("EGS_Ensdf::parseEnsdf: Removing alpha due to small intensity (%.1e < %.1e)\n",(*it)->getAlphaIntensity(),minimumIntensity);
505  }
506  myAlphaRecords.erase(it);
507  }
508  else {
509  it++;
510  }
511  }
512 
513  // Search through the gamma records for any with unknown levels
514  // or with low emission probability
515  bool printedWarning = false;
516  for (vector<GammaRecord * >::iterator it = myGammaRecords.begin();
517  it!=myGammaRecords.end();) {
518 
519  if ((*it)->getTransitionIntensity() <= minimumIntensity) {
520  if (verbose) {
521  egsInformation("EGS_Ensdf::parseEnsdf: Removing gamma due to small intensity (%.1e < %.1e)\n",(*it)->getTransitionIntensity(),minimumIntensity);
522  }
523  // Throw away gammas with low probability
524  // Erase the gamma record object
525  myGammaRecords.erase(it);
526 
527  }
528  else if ((*it)->getLevelRecord()->getEnergy() < epsilon) {
529  // Some gamma may be emitted but the energy level is not known
530  // This is reported in the lnhb data as decays from the -1 level
531  // Since we cannot correlate the emission with a change of energy
532  // states of the daughter, we will treat this transition independently
533 
534  if (!printedWarning) {
535  egsWarning("EGS_Ensdf::parseEnsdf: Warning: Switching internal transition with unknown decay level to uncorrelated event (the emissions will still occur, but uncorrelated with disintegrations).\n");
536  printedWarning = true;
537  }
538 
539  myUncorrelatedGammaRecords.push_back(new GammaRecord(*it));
540 
541  egsInformation("EGS_Ensdf::parseEnsdf: Uncorrelated gamma (E,I): %f %f\n", myUncorrelatedGammaRecords.back()->getDecayEnergy(), myUncorrelatedGammaRecords.back()->getTransitionIntensity());
542 
543  // Erase the gamma record object
544  myGammaRecords.erase(it);
545  }
546  else {
547  ++it;
548  }
549  }
550 
551  // Combine the beta- and beta+ records together
552  for (vector<BetaMinusRecord * >::iterator it = myBetaMinusRecords.begin();
553  it!=myBetaMinusRecords.end(); it++) {
554 
555  myBetaRecords.push_back(*it);
556  }
557 
558  for (vector<BetaPlusRecord * >::iterator it = myBetaPlusRecords.begin();
559  it!=myBetaPlusRecords.end(); it++) {
560 
561  myBetaRecords.push_back(*it);
562  }
563 
564  // Print out a summary of the decays
565  egsInformation("\nEGS_Ensdf::parseEnsdf: Summary of %s emissions:\n", radionuclide.c_str());
566  egsInformation("========================\n");
567  egsInformation("Energy | Intensity per 100 decays\n");
568  if (myBetaRecords.size()) {
569  egsInformation("Beta records:\n");
570  for (vector<BetaRecordLeaf *>::iterator beta = myBetaRecords.begin();
571  beta != myBetaRecords.end(); beta++) {
572  egsInformation("%f %f\n", (*beta)->getFinalEnergy(), (*beta)->getBetaIntensity());
573  }
574  }
575  if (myAlphaRecords.size()) {
576  egsInformation("Alpha records:\n");
577  for (vector<AlphaRecord *>::iterator alpha = myAlphaRecords.begin();
578  alpha != myAlphaRecords.end(); alpha++) {
579  egsInformation("%f %f\n", (*alpha)->getFinalEnergy(), (*alpha)->getAlphaIntensity());
580  }
581  }
582  if (myGammaRecords.size()) {
583  egsInformation("Gamma records (E,Igamma,Ice,Ipp):\n");
584  for (vector<GammaRecord *>::iterator gamma = myGammaRecords.begin();
585  gamma != myGammaRecords.end(); gamma++) {
586  double icI = 0;
587  double ipI = 0;
588  if ((*gamma)->getICIntensity() > 0) {
589  icI = (*gamma)->getGammaIntensity()*(1+(*gamma)->getICIntensity()) - (*gamma)->getGammaIntensity();
590  }
591  if ((*gamma)->getIPIntensity() > 0) {
592  ipI = (*gamma)->getTransitionIntensity() - (*gamma)->getGammaIntensity() - icI;
593  }
594  egsInformation("%f %f %.4e %.4e\n", (*gamma)->getDecayEnergy(), (*gamma)->getGammaIntensity(), icI, ipI);
595  }
596  }
597  if (myUncorrelatedGammaRecords.size()) {
598  egsInformation("Uncorrelated gamma records (E,Igamma,Ice,Ipp):\n");
599  for (vector<GammaRecord *>::iterator gamma = myUncorrelatedGammaRecords.begin();
600  gamma != myUncorrelatedGammaRecords.end(); gamma++) {
601  double icI = 0;
602  double ipI = 0;
603  if ((*gamma)->getICIntensity() > 0) {
604  icI = (*gamma)->getGammaIntensity()*(1+(*gamma)->getICIntensity()) - (*gamma)->getGammaIntensity();
605  }
606  if ((*gamma)->getIPIntensity() > 0) {
607  ipI = (*gamma)->getTransitionIntensity() - (*gamma)->getGammaIntensity() - icI;
608  }
609  egsInformation("%f %f %.4e %.4e\n", (*gamma)->getDecayEnergy(), (*gamma)->getGammaIntensity(), icI, ipI);
610  }
611  }
612  if (xrayEnergies.size() > 0) {
613  egsInformation("X-Ray records:\n");
614  for (unsigned int i=0; i < xrayEnergies.size(); ++i) {
615  egsInformation("%f %f\n", xrayEnergies[i], xrayIntensities[i]);
616  }
617  }
618  if (augerEnergies.size() > 0) {
619  egsInformation("Auger records:\n");
620  for (unsigned int i=0; i < augerEnergies.size(); ++i) {
621  egsInformation("%f %f\n", augerEnergies[i], augerIntensities[i]);
622  }
623  }
624  egsInformation("=== End of summary ===\n");
625 
626  // Determine the final level that the gammas decay towards
627  // We have to use the gamma decay energy to guess at the resulting
628  // energy state of the radionuclide
629  for (vector<GammaRecord *>::iterator gamma = myGammaRecords.begin();
630  gamma != myGammaRecords.end(); gamma++) {
631 
632  double energy = (*gamma)->getDecayEnergy();
633  double guessedLevelEnergy =
634  ((*gamma)->getLevelRecord()->getEnergy() - energy);
635 
636  if (verbose) {
637  egsInformation("EGS_Ensdf::parseEnsdf: Gamma "
638  "(LevelE,E,GuessedE): "
639  "%f %f %f\n",(*gamma)->getLevelRecord()->getEnergy(),
640  energy, guessedLevelEnergy);
641  }
642 
643  double bestMatch = 1E10;
644  LevelRecord *level = 0;
645  for (vector<LevelRecord * >::iterator it = myLevelRecords.begin();
646  it!=myLevelRecords.end(); it++) {
647 
648  double testMatch = fabs((*it)->getEnergy()-guessedLevelEnergy);
649 
650  if (testMatch < bestMatch &&
651  (testMatch < guessedLevelEnergy*0.3 || testMatch < 20)) {
652 
653  bestMatch = testMatch;
654  level = (*it);
655  }
656  }
657  if (bestMatch == 1E10) {
658  egsWarning("EGS_Ensdf::parseEnsdf: Warning: Could "
659  "not find a level with energy matching decay "
660  "of gamma with energy E=%f; "
661  "assuming final level is ground state\n",energy);
662  (*gamma)->setFinalLevel(myLevelRecords.front());
663  }
664  else {
665  (*gamma)->setFinalLevel(level);
666  }
667 
668  if (verbose) {
669  egsInformation("EGS_Ensdf::parseEnsdf: Gamma (final level E, I, Igamma): "
670  "%f %f %f\n",level->getEnergy(), (*gamma)->getTransitionIntensity(), (*gamma)->getGammaIntensity());
671  }
672 
673  (*gamma)->getFinalLevel()->cumulDisintegrationIntensity((*gamma)->getTransitionIntensity());
674  }
675 
676  // Get the total gamma transition intensity of each level
677  unsigned int j = 0;
678  vector<double> totalLevelIntensity;
679  totalLevelIntensity.resize(myLevelRecords.size());
680  for (vector<LevelRecord * >::iterator it = myLevelRecords.begin();
681  it!=myLevelRecords.end(); ++it) {
682 
683  totalLevelIntensity[j] = 0;
684  for (vector<GammaRecord *>::iterator gamma = myGammaRecords.begin();
685  gamma != myGammaRecords.end(); gamma++) {
686 
687  if ((*gamma)->getLevelRecord() == (*it)) {
688  totalLevelIntensity[j] += (*gamma)->getTransitionIntensity();
689  }
690  }
691  ++j;
692  }
693 
694  // For each parent record, search through all disintegration types
695  // (betas, alphas) to see if any exist
696  // If no disintegrations exist for the parent, but we do have internal
697  // transition (IT) gammas, then we must have a metastable radionuclide
698  // In this case, add the gammas to the myMetastableGammaRecords vector
699  if (verbose) {
700  egsInformation("EGS_Ensdf::parseEnsdf: Checking for metastable radionuclides...\n");
701  }
702  for (vector<ParentRecord * >::iterator parent = myParentRecords.begin();
703  parent!=myParentRecords.end(); parent++) {
704 
705  bool gotDisint = false;
706  for (vector<BetaRecordLeaf *>::iterator beta = myBetaRecords.begin();
707  beta != myBetaRecords.end(); beta++) {
708 
709  if ((*beta)->getParentRecord() == *parent) {
710  gotDisint = true;
711  break;
712  }
713  }
714  if (!gotDisint) {
715  for (vector<AlphaRecord *>::iterator alpha = myAlphaRecords.begin();
716  alpha != myAlphaRecords.end(); alpha++) {
717 
718  if ((*alpha)->getParentRecord() == *parent) {
719  gotDisint = true;
720  break;
721  }
722  }
723  }
724 
725  // No disintegrations, so this must be metastable
726  if (!gotDisint) {
727  // We're going to need to "fake" disintegrations toward each level
728  j=0;
729  for (vector<LevelRecord * >::iterator it = myLevelRecords.begin();
730  it!=myLevelRecords.end(); ++it) {
731 
732  double disintIntensity = (*it)->getDisintegrationIntensity();
733  if (disintIntensity < epsilon) {
734  bool gotDecayToLevel = false;
735  for (vector<GammaRecord *>::iterator gamma = myGammaRecords.begin();
736  gamma < myGammaRecords.end(); ++gamma) {
737 
738  // For each gamma matching the current level and parent
739  if ((*gamma)->getParentRecord() == *parent && (*gamma)->getLevelRecord() == (*it)) {
740 
741  // Once per level, we create a zero energy disintegration
742  // with intensity equal to the gamma transition intensity
743  // leaving the level
744  if (!gotDecayToLevel) {
745  gotDecayToLevel = true;
746 
747  // Push a copy of the gamma record
748  myMetastableGammaRecords.push_back(new GammaRecord(*gamma));
749 
750  // Set the transition intensity to be the disintegration
751  // intensity leading towards this level
752  myMetastableGammaRecords.back()->setTransitionIntensity(totalLevelIntensity[j]);
753  }
754  }
755  }
756 
757  if (verbose && myMetastableGammaRecords.size() > 0) {
758  egsInformation("EGS_Ensdf::parseEnsdf: Metastable nuclide "
759  "detected.\n");
760  }
761  }
762 
763  ++j;
764  }
765  }
766  }
767  if (verbose && myMetastableGammaRecords.size() < 1) {
768  egsInformation("EGS_Ensdf::parseEnsdf: No metastable nuclides "
769  "detected.\n");
770  }
771 }
772 
773 // Create record objects from the arrays
774 void EGS_Ensdf::buildRecords() {
775  ParentRecord *lastParent = 0;
776  if (!myParentRecords.empty()) {
777  lastParent = myParentRecords.back();
778  }
779  NormalizationRecord *lastNormalization = 0;
780  if (!myNormalizationRecords.empty()) {
781  lastNormalization = myNormalizationRecords.back();
782  }
783  LevelRecord *lastLevel;
784  if (!myLevelRecords.empty()) {
785  if (!previousParent || previousParent == lastParent) {
786  lastLevel = myLevelRecords.back();
787  }
788  else {
789  lastLevel = new LevelRecord();
790  }
791  }
792  else {
793  lastLevel = new LevelRecord();
794  }
795 
796  for (int i = 0; i < recordStack.size(); i++) {
797  if (!recordStack[i].empty() && recordStack[i].front().length() > 5) {
798  if (i==0) {
799 
800  }
801  else if (i==1) {
802 
803  }
804  else if (i==2) {
805 
806  }
807  else if (i==3) {
808 
809  }
810  else if (i==4) {
811  myCommentRecords.push_back(new CommentRecord(recordStack[i]));
812  }
813  else if (i==5) {
814  myParentRecords.push_back(new ParentRecord(recordStack[i]));
815  }
816  else if (i==6) {
817  myNormalizationRecords.push_back(new
818  NormalizationRecord(recordStack[i], lastParent));
819  }
820  else if (i==7) {
821  myLevelRecords.push_back(new LevelRecord(recordStack[i]));
822  previousParent = lastParent;
823  }
824  else if (i==8) {
825  myBetaMinusRecords.push_back(new
826  BetaMinusRecord(recordStack[i], lastParent,
827  lastNormalization, lastLevel));
828  }
829  else if (i==9) {
830  myBetaPlusRecords.push_back(new
831  BetaPlusRecord(recordStack[i], lastParent,
832  lastNormalization, lastLevel));
833  }
834  else if (i==10) {
835  myAlphaRecords.push_back(new
836  AlphaRecord(recordStack[i], lastParent,
837  lastNormalization, lastLevel));
838  }
839  else if (i==11) {
840  egsWarning("EGS_Ensdf::buildRecords: Warning: Delayed particle not "
841  "supported! Further development required.\n");
842  }
843  else if (i==12) {
844  myGammaRecords.push_back(new
845  GammaRecord(recordStack[i], lastParent,
846  lastNormalization, lastLevel));
847  }
848 
849  recordStack[i].clear();
850  }
851  }
852 }
853 
854 // Normalize intensities for alpha, beta, gamma objects
855 void EGS_Ensdf::normalizeIntensities() {
856  if (verbose) {
857  egsInformation("EGS_Ensdf::normalizeIntensities: Normalizing the "
858  "emission intensities to allow for spectrum sampling "
859  "routines...\n");
860  }
861 
862  // Add up the beta, alpha, xray and auger decay intensities
863  double totalDecayIntensity = 0;
864  double totalDecayIntensityUnc = 0;
865  double lastIntensity = 0;
866 
867  for (vector<BetaMinusRecord * >::iterator it = myBetaMinusRecords.begin();
868  it!=myBetaMinusRecords.end(); it++) {
869 
870  if (verbose > 1) {
871  egsInformation("EGS_Ensdf::normalizeIntensities: Beta- (E,I): %f %f\n",
872  (*it)->getFinalEnergy(), (*it)->getBetaIntensity());
873  }
874 
875  totalDecayIntensity += (*it)->getBetaIntensity();
876  totalDecayIntensityUnc += (*it)->getBetaIntensityUnc();
877  }
878  for (vector<BetaPlusRecord * >::iterator it = myBetaPlusRecords.begin();
879  it!=myBetaPlusRecords.end(); it++) {
880 
881  if (verbose > 1) {
882  egsInformation("EGS_Ensdf::normalizeIntensities: Beta+/EC (E,I): %f %f\n",
883  (*it)->getFinalEnergy(), (*it)->getBetaIntensity());
884  }
885 
886  totalDecayIntensity += (*it)->getBetaIntensity();
887  totalDecayIntensityUnc += (*it)->getPositronIntensityUnc() + (*it)->getECIntensityUnc();
888  }
889  for (vector<AlphaRecord *>::iterator alpha = myAlphaRecords.begin();
890  alpha != myAlphaRecords.end(); alpha++) {
891 
892  if (verbose > 1) {
893  egsInformation("EGS_Ensdf::normalizeIntensities: Alpha (E,I): %f %f\n",
894  (*alpha)->getFinalEnergy(), (*alpha)->getAlphaIntensity());
895  }
896 
897  totalDecayIntensity += (*alpha)->getAlphaIntensity();
898  totalDecayIntensityUnc += (*alpha)->getAlphaIntensityUnc();
899  }
900  for (vector<GammaRecord *>::iterator gamma = myMetastableGammaRecords.begin();
901  gamma != myMetastableGammaRecords.end(); gamma++) {
902 
903  if (verbose > 1) {
904  egsInformation("EGS_Ensdf::normalizeIntensities: MetastableGamma (I): %f\n",
905  (*gamma)->getTransitionIntensity());
906  }
907 
908  totalDecayIntensity += (*gamma)->getTransitionIntensity();
909  totalDecayIntensityUnc += (*gamma)->getGammaIntensityUnc() + (*gamma)->getICIntensityUnc() + (*gamma)->getIPIntensityUnc();
910  }
911 
912  // Check that the branch probabilities add up to one
913  double branchSum = 0;
914  for (vector<NormalizationRecord * >::iterator norm =
915  myNormalizationRecords.begin();
916  norm!=myNormalizationRecords.end(); norm++) {
917  branchSum += (*norm)->getBranchMultiplier();
918  }
919  // Currently there is only 1 case in the LNHB ensdf data where this is true
920  // It is for Cf-252 fission events
921  if (branchSum < 1-epsilon) {
922  egsWarning("\nEGS_Ensdf::normalizeIntensities: Warning: The branching ratios of this nuclide add to less than 1 (%f). The leftover probability will be assigned to fission events. These events will return a zero energy particle and be counted as disintegrations. This is expected for Cf-252 in the LNHB collection.\n\n",branchSum);
923 
924  // Add the fission probability to the total decay intensity
925  totalDecayIntensity /= branchSum;
926  }
927  else if (branchSum > 1+epsilon) {
928  egsWarning("\nEGS_Ensdf::normalizeIntensities: Warning: The branching ratios of this nuclide add to greater than 1 (%f). This will result in overall emission rates being incorrect (e.g. number of emissions per 100 decays) when compared against the input.\n\n",branchSum);
929  }
930 
931  // At this stage, totalDecayIntensity would ideally be 100, to represent
932  // that all the decay probabilities added to 100%. However, this is not usually
933  // achieved, due to uncertainties in the values. Thus we must adjust the numbers
934  // somehow in order to model decays.
935  // Option 1: Do nothing special, and the decays will get normalized to to force
936  // them to add to 100%. This evenly spreads the discrepancy over all decays,
937  // despite the fact that some of them might have higher uncertainties.
938  // Option 2: Spread the discrepancy over all decays, scaled by the uncertainty
939  // of each decay.
940 
941  // Here we do option 2:
942  // Check that totalDecayIntensity != 100
943  if (totalDecayIntensity > 100 + epsilon || totalDecayIntensity < 100 - epsilon) {
944 
945  // This is the discrepancy of the total decay intensity we found in the file,
946  // that we will have to account for. This difference will be spread over
947  // all of the decays, by adjusting the decay intensities proportional
948  // to their uncertainty
949  decayDiscrepancy = 100 - totalDecayIntensity;
950 
951  egsInformation("EGS_Ensdf::normalizeIntensities: Warning: The sum of the decay probabilities for this nuclide does not equal 100\%! In order for modeling to proceed, this must be accounted for. The discrepancy of %f has been distributed over all decays, proportional to the corresponding uncertainties. Note that this will also change the internal transition intensities, since they depend on the decays.\n", decayDiscrepancy);
952 
953  // Reset the disintegration intensities for each level, we will
954  // need to recalculate this in order to adjust the transition
955  // intensities after adjusting the decay intensities
956  for (vector<LevelRecord * >::iterator it = myLevelRecords.begin();
957  it!=myLevelRecords.end(); it++) {
958  (*it)->resetDisintegrationIntensity();
959  }
960 
961  for (vector<BetaMinusRecord * >::iterator it = myBetaMinusRecords.begin();
962  it!=myBetaMinusRecords.end(); it++) {
963 
964  (*it)->setBetaIntensity((*it)->getBetaIntensity() + (*it)->getBetaIntensityUnc() / totalDecayIntensityUnc * decayDiscrepancy);
965 
966  // Add to the intensity toward this level
967  (*it)->getLevelRecord()->cumulDisintegrationIntensity((*it)->getBetaIntensity());
968  }
969 
970  for (vector<BetaPlusRecord * >::iterator it = myBetaPlusRecords.begin();
971  it!=myBetaPlusRecords.end(); it++) {
972 
973  // The positron intensity is already normalized with the electron
974  // capture intensity so they add to 1. So we multiply by the beta
975  // intensity to get the correct units on each, and calculate
976  // the new values
977  double newPositronIntensity = (*it)->getBetaIntensity() * (*it)->getPositronIntensity() + (*it)->getPositronIntensityUnc() / totalDecayIntensityUnc * decayDiscrepancy;
978  double newECIntensity = (*it)->getBetaIntensity() * (1-(*it)->getPositronIntensity()) + (*it)->getECIntensityUnc() / totalDecayIntensityUnc * decayDiscrepancy;
979 
980  // Set the total intensity for this decay based on the new values
981  (*it)->setBetaIntensity(newPositronIntensity + newECIntensity);
982 
983  // Normalize the positron intensity for sampling (either positron
984  // decay or EC occurs, so they add to 1)
985  (*it)->setPositronIntensity(newPositronIntensity / (*it)->getBetaIntensity());
986 
987  // Add to the intensity toward this level
988  (*it)->getLevelRecord()->cumulDisintegrationIntensity((*it)->getBetaIntensity());
989  }
990 
991  for (vector<AlphaRecord * >::iterator it = myAlphaRecords.begin();
992  it!=myAlphaRecords.end(); it++) {
993 
994  (*it)->setAlphaIntensity((*it)->getAlphaIntensity() + (*it)->getAlphaIntensityUnc() / totalDecayIntensityUnc * decayDiscrepancy);
995 
996  // Add to the intensity toward this level
997  (*it)->getLevelRecord()->cumulDisintegrationIntensity((*it)->getAlphaIntensity());
998  }
999 
1000  for (vector<GammaRecord *>::iterator it = myMetastableGammaRecords.begin();
1001  it != myMetastableGammaRecords.end(); it++) {
1002 
1003  double icI = 0;
1004  double ipI = 0;
1005  if ((*it)->getICIntensity() > 0) {
1006  icI = (*it)->getGammaIntensity()*(1+(*it)->getICIntensity()) - (*it)->getGammaIntensity();
1007  }
1008  if ((*it)->getIPIntensity() > 0) {
1009  ipI = (*it)->getTransitionIntensity() - (*it)->getGammaIntensity() - ((*it)->getGammaIntensity()*(1+(*it)->getICIntensity()) - (*it)->getGammaIntensity());
1010  }
1011 
1012  // We need to calculate the original intensities, adjust
1013  // based on the uncertainty scaling, and normalize again
1014  double newGammaIntensity = (*it)->getGammaIntensity() + (*it)->getGammaIntensityUnc() / totalDecayIntensityUnc * decayDiscrepancy;
1015  double newICIntensity = icI + (*it)->getICIntensityUnc() / totalDecayIntensityUnc * decayDiscrepancy;
1016  double newIPIntensity = ipI + (*it)->getIPIntensityUnc() / totalDecayIntensityUnc * decayDiscrepancy;
1017 
1018  // Set the intensities for this decay based on the new values
1019  (*it)->setTransitionIntensity(newGammaIntensity + newICIntensity + newIPIntensity);
1020  (*it)->setGammaIntensity(newGammaIntensity);
1021  (*it)->setICIntensity((newICIntensity+(*it)->getGammaIntensity()) / (*it)->getGammaIntensity() - 1);
1022 
1023  // Add to the intensity toward this level
1024  (*it)->getLevelRecord()->cumulDisintegrationIntensity((*it)->getTransitionIntensity());
1025  }
1026 
1027  // For regular internal transitions there is no normalization needed
1028  for (vector<GammaRecord *>::iterator gamma = myGammaRecords.begin();
1029  gamma != myGammaRecords.end(); gamma++) {
1030 
1031  // Add the contribution of this transition *toward* a different level
1032  (*gamma)->getFinalLevel()->cumulDisintegrationIntensity((*gamma)->getTransitionIntensity());
1033  }
1034 
1035  // Print out a summary of the decays
1036  egsInformation("\nEGS_Ensdf::normalizeIntensities: Summary of %s decays (adjusted by %f):\n", radionuclide.c_str(), decayDiscrepancy);
1037  egsInformation("========================\n");
1038  egsInformation("Energy | Intensity per 100 decays\n");
1039  if (myBetaRecords.size()) {
1040  egsInformation("Beta records:\n");
1041  for (vector<BetaRecordLeaf *>::iterator beta = myBetaRecords.begin();
1042  beta != myBetaRecords.end(); beta++) {
1043  egsInformation("%f %f\n", (*beta)->getFinalEnergy(), (*beta)->getBetaIntensity());
1044  }
1045  }
1046  if (myAlphaRecords.size()) {
1047  egsInformation("Alpha records:\n");
1048  for (vector<AlphaRecord *>::iterator alpha = myAlphaRecords.begin();
1049  alpha != myAlphaRecords.end(); alpha++) {
1050  egsInformation("%f %f\n", (*alpha)->getFinalEnergy(), (*alpha)->getAlphaIntensity());
1051  }
1052  }
1053  if (myGammaRecords.size()) {
1054  egsInformation("Gamma records (E,Igamma,Ice,Ipp):\n");
1055  for (vector<GammaRecord *>::iterator gamma = myGammaRecords.begin();
1056  gamma != myGammaRecords.end(); gamma++) {
1057  double icI = 0;
1058  double ipI = 0;
1059  if ((*gamma)->getICIntensity() > 0) {
1060  icI = (*gamma)->getGammaIntensity()*(1+(*gamma)->getICIntensity()) - (*gamma)->getGammaIntensity();
1061  }
1062  if ((*gamma)->getIPIntensity() > 0) {
1063  ipI = (*gamma)->getTransitionIntensity() - (*gamma)->getGammaIntensity() - icI;
1064  }
1065  egsInformation("%f %f %.4e %.4e\n", (*gamma)->getDecayEnergy(), (*gamma)->getGammaIntensity(), icI, ipI);
1066  }
1067  }
1068  if (myUncorrelatedGammaRecords.size()) {
1069  egsInformation("Uncorrelated gamma records (E,Igamma,Ice,Ipp):\n");
1070  for (vector<GammaRecord *>::iterator gamma = myUncorrelatedGammaRecords.begin();
1071  gamma != myUncorrelatedGammaRecords.end(); gamma++) {
1072  double icI = 0;
1073  double ipI = 0;
1074  if ((*gamma)->getICIntensity() > 0) {
1075  icI = (*gamma)->getGammaIntensity()*(1+(*gamma)->getICIntensity()) - (*gamma)->getGammaIntensity();
1076  }
1077  if ((*gamma)->getIPIntensity() > 0) {
1078  ipI = (*gamma)->getTransitionIntensity() - (*gamma)->getGammaIntensity() - icI;
1079  }
1080  egsInformation("%f %f %.4e %.4e\n", (*gamma)->getDecayEnergy(), (*gamma)->getGammaIntensity(), icI, ipI);
1081  }
1082  }
1083  if (xrayEnergies.size() > 0) {
1084  egsInformation("X-Ray records:\n");
1085  for (unsigned int i=0; i < xrayEnergies.size(); ++i) {
1086  egsInformation("%f %f\n", xrayEnergies[i], xrayIntensities[i]);
1087  }
1088  }
1089  if (augerEnergies.size() > 0) {
1090  egsInformation("Auger records:\n");
1091  for (unsigned int i=0; i < augerEnergies.size(); ++i) {
1092  egsInformation("%f %f\n", augerEnergies[i], augerIntensities[i]);
1093  }
1094  }
1095  egsInformation("=== End of summary ===\n");
1096 
1097  // Now we have guaranteed that the decay probabilities add to 100%
1098  totalDecayIntensity = 100;
1099  }
1100 
1101  // Check for instances where the intensity of disintegrations that
1102  // lead towards a particular excited daughter level is less than
1103  // the intensities of gamma transitions from it.
1104  if (allowMultiTransition) {
1105  if (verbose) {
1106  egsInformation("EGS_Ensdf::normalizeIntensities: Comparing the cumulative disintegration intensity of each level with the gamma transition intensities... \n");
1107  }
1108 
1109  // Get the total transition intensity *away* from each level
1110  unsigned int j = 0;
1111  vector<double> totalLevelIntensity;
1112  totalLevelIntensity.resize(myLevelRecords.size());
1113  for (vector<LevelRecord * >::iterator it = myLevelRecords.begin();
1114  it!=myLevelRecords.end(); ++it) {
1115 
1116  totalLevelIntensity[j] = 0;
1117  for (vector<GammaRecord *>::iterator gamma = myGammaRecords.begin();
1118  gamma != myGammaRecords.end(); gamma++) {
1119 
1120  if ((*gamma)->getLevelRecord() == (*it)) {
1121  totalLevelIntensity[j] += (*gamma)->getTransitionIntensity();
1122  }
1123  }
1124  ++j;
1125  }
1126 
1127  // Iterate in reverse through the levels (highest level first)
1128  // This ensures we modify the intensities feeding the lower levels first
1129  j = myLevelRecords.size()-1;
1130  for (vector<LevelRecord * >::reverse_iterator it = myLevelRecords.rbegin();
1131  it!=myLevelRecords.rend(); ++it) {
1132 
1133  double disintIntensity = (*it)->getDisintegrationIntensity();
1134 
1135  if (verbose) {
1136  egsInformation("EGS_Ensdf::normalizeIntensities: (Level, ItoLevel, IfromLevel): %d %f %f\n", j, disintIntensity, totalLevelIntensity[j]);
1137  }
1138 
1139  // Notice that we don't do this if disintIntensity==0
1140  if (disintIntensity > epsilon && totalLevelIntensity[j] > disintIntensity + epsilon) {
1141  for (vector<GammaRecord *>::iterator gamma = myGammaRecords.begin();
1142  gamma != myGammaRecords.end(); gamma++) {
1143 
1144  if ((*gamma)->getLevelRecord() == (*it)) {
1145  double multipleTransitionProb = (1.-disintIntensity/totalLevelIntensity[j]);
1146 
1147  (*gamma)->setMultiTransitionProb(multipleTransitionProb);
1148 
1149  if (verbose) {
1150  egsInformation("EGS_Ensdf::normalizeIntensities: Multiple gamma transition probability (E,I): %f %f\n",(*gamma)->getDecayEnergy(), multipleTransitionProb);
1151  }
1152 
1153  // Reduce the transition intensities
1154  (*gamma)->setTransitionIntensity(
1155  (*gamma)->getTransitionIntensity() * disintIntensity/totalLevelIntensity[j]
1156  );
1157  (*gamma)->setGammaIntensity(
1158  (*gamma)->getGammaIntensity() * disintIntensity/totalLevelIntensity[j]
1159  );
1160  }
1161  }
1162  }
1163 
1164  --j;
1165  }
1166  }
1167 
1168  for (vector<GammaRecord *>::iterator gamma = myUncorrelatedGammaRecords.begin();
1169  gamma != myUncorrelatedGammaRecords.end(); gamma++) {
1170  totalDecayIntensity += (*gamma)->getTransitionIntensity();
1171  }
1172 
1173  // Add X-rays and Auger electrons that are read from the ENSDF file
1174  // and going to be used instead of modelling correlated atomic relaxations.
1175  // These are modeled like independent decays, though they don't count
1176  // as disintegration events in the fluence.
1177  for (unsigned int i=0; i < xrayIntensities.size(); ++i) {
1178  if (verbose > 1) {
1179  egsInformation("EGS_Ensdf::normalizeIntensities: XRay (E,I): %f %f\n",
1180  xrayEnergies[i], xrayIntensities[i]);
1181  }
1182 
1183  totalDecayIntensity += xrayIntensities[i];
1184  }
1185  for (unsigned int i=0; i < augerIntensities.size(); ++i) {
1186  if (verbose > 1) {
1187  egsInformation("EGS_Ensdf::normalizeIntensities: Auger (E,I): %f %f\n",
1188  augerEnergies[i], augerIntensities[i]);
1189  }
1190 
1191  totalDecayIntensity += augerIntensities[i];
1192  }
1193 
1194  if (verbose) {
1195  egsInformation("EGS_Ensdf::normalizeIntensities: totalDecayIntensity: "
1196  "%f\n\n",totalDecayIntensity);
1197  egsInformation("EGS_Ensdf::normalizeIntensities: "
1198  "Calculating renormalized intensities...\n");
1199  }
1200 
1201  // Normalize beta emission intensities
1202  for (vector<BetaRecordLeaf *>::iterator beta = myBetaRecords.begin();
1203  beta != myBetaRecords.end(); beta++) {
1204 
1205  (*beta)->setBetaIntensity(
1206  (*beta)->getBetaIntensity() / totalDecayIntensity);
1207 
1208  if ((beta - myBetaRecords.begin()) > 0) {
1209  (*beta)->setBetaIntensity(
1210  (*beta)->getBetaIntensity() + (*(beta-1))->getBetaIntensity());
1211  }
1212  lastIntensity = (*beta)->getBetaIntensity();
1213 
1214  if (verbose) {
1215  egsInformation("EGS_Ensdf::normalizeIntensities: Beta (E,I): %f %f\n",
1216  (*beta)->getFinalEnergy(), (*beta)->getBetaIntensity());
1217  }
1218  }
1219 
1220  // Normalize alpha emission intensities
1221  for (vector<AlphaRecord *>::iterator alpha = myAlphaRecords.begin();
1222  alpha != myAlphaRecords.end(); alpha++) {
1223 
1224  (*alpha)->setAlphaIntensity(
1225  (*alpha)->getAlphaIntensity() / totalDecayIntensity);
1226 
1227  if ((alpha - myAlphaRecords.begin()) == 0 && lastIntensity > epsilon) {
1228  (*alpha)->setAlphaIntensity(
1229  (*alpha)->getAlphaIntensity() + lastIntensity);
1230  }
1231  else if ((alpha - myAlphaRecords.begin()) > 0) {
1232  (*alpha)->setAlphaIntensity(
1233  (*alpha)->getAlphaIntensity() +
1234  (*(alpha-1))->getAlphaIntensity());
1235  }
1236  lastIntensity = (*alpha)->getAlphaIntensity();
1237 
1238  if (verbose) {
1239  egsInformation("EGS_Ensdf::normalizeIntensities: Alpha (E,I): %f %f\n",
1240  (*alpha)->getFinalEnergy(), (*alpha)->getAlphaIntensity());
1241  }
1242  }
1243 
1244  // Normalize metastable gamma transition intensities
1245  for (vector<GammaRecord *>::iterator gamma = myMetastableGammaRecords.begin();
1246  gamma != myMetastableGammaRecords.end(); gamma++) {
1247 
1248  (*gamma)->setTransitionIntensity(
1249  (*gamma)->getTransitionIntensity() / totalDecayIntensity);
1250 
1251  if ((gamma - myMetastableGammaRecords.begin()) == 0 && lastIntensity > epsilon) {
1252  (*gamma)->setTransitionIntensity(
1253  (*gamma)->getTransitionIntensity() + lastIntensity);
1254  }
1255  else if ((gamma - myMetastableGammaRecords.begin()) > 0) {
1256  (*gamma)->setTransitionIntensity(
1257  (*gamma)->getTransitionIntensity() +
1258  (*(gamma-1))->getTransitionIntensity());
1259  }
1260  lastIntensity = (*gamma)->getTransitionIntensity();
1261 
1262  if (verbose) {
1263  egsInformation("EGS_Ensdf::normalizeIntensities: MetastableGamma (I): %f\n",
1264  (*gamma)->getTransitionIntensity());
1265  }
1266  }
1267 
1268  // Normalize uncorrelated internal transitions
1269  for (vector<GammaRecord *>::iterator gamma = myUncorrelatedGammaRecords.begin();
1270  gamma != myUncorrelatedGammaRecords.end(); gamma++) {
1271  (*gamma)->setTransitionIntensity(
1272  (*gamma)->getTransitionIntensity() / totalDecayIntensity);
1273 
1274  if ((gamma - myUncorrelatedGammaRecords.begin()) == 0 && lastIntensity > epsilon) {
1275  (*gamma)->setTransitionIntensity(
1276  (*gamma)->getTransitionIntensity() + lastIntensity);
1277  }
1278  else if ((gamma - myUncorrelatedGammaRecords.begin()) > 0) {
1279  (*gamma)->setTransitionIntensity(
1280  (*gamma)->getTransitionIntensity() +
1281  (*(gamma-1))->getTransitionIntensity());
1282  }
1283  lastIntensity = (*gamma)->getTransitionIntensity();
1284 
1285  if (verbose) {
1286  egsInformation("EGS_Ensdf::normalizeIntensities: UncorrelatedGamma (I): %f\n",
1287  (*gamma)->getTransitionIntensity());
1288  }
1289  }
1290 
1291  // Normalize XRay emission intensities
1292  for (unsigned int i=0; i < xrayIntensities.size(); ++i) {
1293 
1294  xrayIntensities[i] /= totalDecayIntensity;
1295 
1296  if (i==0 && lastIntensity > epsilon) {
1297  xrayIntensities[i] += lastIntensity;
1298  }
1299  else if (i > 0) {
1300  xrayIntensities[i] += xrayIntensities[i-1];
1301  }
1302  lastIntensity = xrayIntensities[i];
1303 
1304  if (verbose) {
1305  egsInformation("EGS_Ensdf::normalizeIntensities: XRay (E,I): %f %f\n",
1306  xrayEnergies[i], xrayIntensities[i]);
1307  }
1308  }
1309 
1310  // Normalize auger emission intensities
1311  for (unsigned int i=0; i < augerIntensities.size(); ++i) {
1312 
1313  augerIntensities[i] /= totalDecayIntensity;
1314 
1315  if (i==0 && lastIntensity > epsilon) {
1316  augerIntensities[i] += lastIntensity;
1317  }
1318  else if (i > 0) {
1319  augerIntensities[i] += augerIntensities[i-1];
1320  }
1321  lastIntensity = augerIntensities[i];
1322 
1323  if (verbose) {
1324  egsInformation("EGS_Ensdf::normalizeIntensities: Auger (E,I): %f %f\n",
1325  augerEnergies[i], augerIntensities[i]);
1326  }
1327  }
1328 
1329  // Get the gamma transition intensities
1330  // and the total intensity for each level
1331  unsigned int j = 0;
1332  vector<double> totalLevelIntensity;
1333  totalLevelIntensity.resize(myLevelRecords.size());
1334  for (vector<LevelRecord * >::iterator it = myLevelRecords.begin();
1335  it!=myLevelRecords.end(); it++) {
1336 
1337  double disintIntensity = (*it)->getDisintegrationIntensity();
1338 
1339  totalLevelIntensity[j] = 0;
1340  for (vector<GammaRecord *>::iterator gamma = myGammaRecords.begin();
1341  gamma != myGammaRecords.end(); gamma++) {
1342 
1343  if ((*gamma)->getLevelRecord() == (*it)) {
1344  totalLevelIntensity[j] += (*gamma)->getTransitionIntensity();
1345  }
1346  }
1347 
1348  if (verbose > 1) {
1349  egsInformation("EGS_Ensdf::normalizeIntensities: "
1350  "totalLevelIntensity: %f\n", totalLevelIntensity[j]);
1351  }
1352 
1353  if (disintIntensity > epsilon && totalLevelIntensity[j] < disintIntensity + epsilon) {
1354  totalLevelIntensity[j] = disintIntensity;
1355  }
1356  ++j;
1357  }
1358 
1359  // Normalize transition intensities over each level
1360  j = 0;
1361  for (vector<LevelRecord * >::iterator it = myLevelRecords.begin();
1362  it!=myLevelRecords.end(); it++) {
1363 
1364  unsigned int i = 0;
1365  bool levelCanDecay = false;
1366  for (vector<GammaRecord *>::iterator gamma = myGammaRecords.begin();
1367  gamma != myGammaRecords.end(); gamma++) {
1368 
1369  if ((*gamma)->getLevelRecord() == (*it)) {
1370  levelCanDecay = true;
1371 
1372  (*gamma)->setGammaIntensity(
1373  (*gamma)->getGammaIntensity() /
1374  (*gamma)->getTransitionIntensity());
1375 
1376  (*gamma)->setICIntensity(
1377  (*gamma)->getGammaIntensity() * (1+(*gamma)->getICIntensity()));
1378 
1379  if (totalLevelIntensity[j] > epsilon) {
1380  (*gamma)->setTransitionIntensity(
1381  (*gamma)->getTransitionIntensity() /
1382  totalLevelIntensity[j]);
1383  }
1384 
1385  if (i > 0) {
1386  (*gamma)->setTransitionIntensity(
1387  (*gamma)->getTransitionIntensity() +
1388  (*(gamma-1))->getTransitionIntensity());
1389  }
1390  ++i;
1391 
1392  if (verbose > 1) {
1393  egsInformation("EGS_Ensdf::normalizeIntensities: "
1394  "Gamma (level,E,I,Igamma,Ice): "
1395  "%d %f %f %f %f\n",
1396  j,(*gamma)->getDecayEnergy(), (*gamma)->getTransitionIntensity(),
1397  (*gamma)->getGammaIntensity(), (*gamma)->getICIntensity());
1398  }
1399  }
1400  }
1401 
1402  // Set whether or not the level can decay
1403  // If there are no gammas that decay from this excited energy state,
1404  // then it will effectively stay in this state forever
1405  // In practice this means we don't have to sample for the emission of
1406  // transition photons
1407  (*it)->setLevelCanDecay(levelCanDecay);
1408 
1409  ++j;
1410  }
1411 }
1412 
1413 void EGS_Ensdf::getEmissionsFromComments() {
1414  if (verbose) {
1415  egsInformation("EGS_Ensdf::getEmissionsFromComments: Attempting to obtain x-ray and Auger emissions from the ENSDF comments. This assumes a particular comment format...\n");
1416  }
1417 
1418  bool containsEmissions = false;
1419  bool xrayContinues = false;
1420  bool augerContinues = false;
1421  bool gotTotal = false;
1422  vector<double> multilineEnergies,
1423  multilineIntensities;
1424  double lineTotalIntensity = 0;
1425  unsigned int countNumAfterTotal = 0;
1426  int lineTotalType;
1427 
1428  for (vector<CommentRecord *>::iterator comment = myCommentRecords.begin();
1429  comment != myCommentRecords.end(); comment++) {
1430 
1431  auto commentSet = (*comment)->getComments();
1432  for (auto c = commentSet.begin(); c != commentSet.end(); ++c) {
1433  string line = *c;
1434 
1435  // Search for this line to ensure emissions will follow the right format
1436  if (line.find("{U Energy (keV)} {U Intensity} {U Line}") != std::string::npos) {
1437  containsEmissions = true;
1438  }
1439 
1440  if (containsEmissions) {
1441  // Check for the end of multi-line records
1442  // and average them together
1443  if (line.length() < 48 ||
1444  ((xrayContinues || augerContinues) && line.at(30) != '|')) {
1445 
1446  // If we just finished going through a series of
1447  // lines that started with a "total" line at the top
1448  // then we'll check to make sure they have intensities assigned
1449  if (gotTotal) {
1450 
1451  // In the event that a zero intensity is in one of the lines
1452  // following the "total" line, ALL of those following
1453  // lines are assigned an equal fraction of the total
1454  // intensity. This is an imperfect work-around for insufficient
1455  // data. Using this method, the correct energies are used,
1456  // rather than assigning a single averaged "total" energy line.
1457  if (countNumAfterTotal > 0) {
1458  // X-rays
1459  if (lineTotalType == 0) {
1460  bool containsZeroIntensity = false;
1461  for (std::vector<double>::iterator it = xrayIntensities.end()-countNumAfterTotal; it != xrayIntensities.end(); ++it) {
1462  if (*it < epsilon) {
1463  containsZeroIntensity = true;
1464  break;
1465  }
1466  }
1467 
1468  if (containsZeroIntensity) {
1469  for (std::vector<double>::iterator it = xrayIntensities.end()-countNumAfterTotal; it != xrayIntensities.end(); ++it) {
1470  *it = lineTotalIntensity / countNumAfterTotal;
1471  }
1472  }
1473 
1474  // Auger
1475  }
1476  else if (lineTotalType == -1) {
1477  bool containsZeroIntensity = false;
1478  for (std::vector<double>::iterator it = augerIntensities.end()-countNumAfterTotal; it != augerIntensities.end(); ++it) {
1479  if (*it < epsilon) {
1480  containsZeroIntensity = true;
1481  break;
1482  }
1483  }
1484 
1485  if (containsZeroIntensity) {
1486  for (std::vector<double>::iterator it = augerIntensities.end()-countNumAfterTotal; it != augerIntensities.end(); ++it) {
1487  *it = lineTotalIntensity / countNumAfterTotal;
1488  }
1489  }
1490  }
1491  }
1492 
1493  gotTotal = false;
1494  countNumAfterTotal = 0;
1495  lineTotalIntensity = 0.;
1496  }
1497 
1498  if ((xrayContinues || augerContinues)
1499  && multilineEnergies.size() > 0) {
1500 
1501  double energySum = 0;
1502  double intensitySum = 0;
1503  unsigned int numNonzeroE = 0;
1504  unsigned int numNonzeroI = 0;
1505  for (unsigned int i=0; i < multilineEnergies.size(); ++i) {
1506  if (multilineEnergies[i] > 0) {
1507  energySum += multilineEnergies[i];
1508  numNonzeroE++;
1509  }
1510  }
1511  for (unsigned int i=0; i < multilineIntensities.size(); ++i) {
1512  if (multilineIntensities[i] > epsilon) {
1513  intensitySum += multilineIntensities[i];
1514  numNonzeroI++;
1515  }
1516  }
1517  double energy;
1518  if (numNonzeroE > 0) {
1519  energy = energySum / numNonzeroE;
1520  }
1521  double intensity;
1522  if (numNonzeroI > 0) {
1523  intensity = intensitySum / numNonzeroI;
1524  }
1525 
1526  if (numNonzeroE > 0 && numNonzeroI > 0) {
1527  if (xrayContinues) {
1528  xrayEnergies.push_back(energy);
1529  xrayIntensities.push_back(intensity);
1530  }
1531  else {
1532  augerEnergies.push_back(energy);
1533  augerIntensities.push_back(intensity);
1534  }
1535  }
1536 
1537  multilineEnergies.clear();
1538  multilineIntensities.clear();
1539  }
1540 
1541  xrayContinues = false;
1542  augerContinues = false;
1543  }
1544 
1545  // Check for records containing XRays or Auger electrons
1546  if (line.length() > 48) {
1547 
1548  string emissionLine = egsTrimString(line.substr(47));
1549 
1550  // See if the line is an XRay or Auger
1551  if (emissionLine.length() < 1 || (emissionLine.at(0) != 'X' && emissionLine.find("AUGER") == std::string::npos)) {
1552  continue;
1553  }
1554 
1555  string eStr = egsTrimString(line.substr(13, 15));
1556 
1557  // If we have a range in energy (e.g. 0.1-0.3)
1558  // Find the average
1559  size_t eDash = eStr.find('-');
1560  double energy;
1561  if (eDash!=std::string::npos) {
1562  if (eStr.length() > eDash+1) {
1563  double e1 = atof(eStr.substr(0, eDash).c_str());
1564  double e2 = atof(eStr.substr(eDash+1).c_str());
1565  energy = (e1 + e2) / 2;
1566  }
1567  else {
1568  energy = atof(eStr.substr(0, eDash).c_str());
1569  }
1570  }
1571  else {
1572  energy = atof(eStr.c_str());
1573  }
1574 
1575  // Convert the energy from keV to MeV
1576  energy /= 1000.;
1577 
1578  // Get the intensity
1579  string iStr = egsTrimString(line.substr(32, 9));
1580  double intensity = atof(iStr.c_str());
1581 
1582  // If this is a line coming after a "total" line,
1583  // increment a counter. This will be used in the
1584  // event that the lines following the "total"
1585  // have zero intensity assigned
1586  if (gotTotal && energy > epsilon) {
1587  countNumAfterTotal++;
1588  }
1589 
1590  // If this line is the total of the next lines, we will
1591  // skip this line and use the individual ones
1592  // However, record the total intensity in case we need it
1593  if (emissionLine.find("(total)") != std::string::npos) {
1594  gotTotal = true;
1595  lineTotalIntensity = intensity;
1596  if (emissionLine.find("AUGER") != std::string::npos) {
1597  lineTotalType = -1;
1598  }
1599  else {
1600  lineTotalType = 0;
1601  }
1602  continue;
1603  }
1604 
1605  // Multi-line records have a bar '|' at 30
1606  // We will store the data and average them later
1607  if (line.at(30) == '|') {
1608  if (emissionLine.at(0) == 'X') {
1609  xrayContinues = true;
1610  }
1611  else if (emissionLine.find("AUGER") != std::string::npos) {
1612  augerContinues = true;
1613  }
1614 
1615  multilineEnergies.push_back(energy);
1616  multilineIntensities.push_back(intensity);
1617 
1618  }
1619  else {
1620  if (emissionLine.at(0) == 'X') {
1621  if ((energy > epsilon && intensity > epsilon) ||
1622  (gotTotal && energy > epsilon)) {
1623  xrayEnergies.push_back(energy);
1624  xrayIntensities.push_back(intensity);
1625  }
1626  }
1627  else if (emissionLine.find("AUGER") != std::string::npos) {
1628  if ((energy > epsilon && intensity > epsilon) ||
1629  (gotTotal && energy > epsilon)) {
1630  augerEnergies.push_back(energy);
1631  augerIntensities.push_back(intensity);
1632  }
1633  }
1634  }
1635  }
1636  }
1637  }
1638  }
1639 }
1640 
1641 vector<double > EGS_Ensdf::getXRayIntensities() const {
1642  return xrayIntensities;
1643 }
1644 
1645 vector<double > EGS_Ensdf::getXRayEnergies() const {
1646  return xrayEnergies;
1647 }
1648 
1649 vector<double > EGS_Ensdf::getAugerIntensities() const {
1650  return augerIntensities;
1651 }
1652 
1653 vector<double > EGS_Ensdf::getAugerEnergies() const {
1654  return augerEnergies;
1655 }
1656 
1657 vector<ParentRecord * > EGS_Ensdf::getParentRecords() const {
1658  return myParentRecords;
1659 }
1660 
1661 vector<LevelRecord * > EGS_Ensdf::getLevelRecords() const {
1662  return myLevelRecords;
1663 }
1664 
1665 vector<BetaRecordLeaf * > EGS_Ensdf::getBetaRecords() const {
1666  return myBetaRecords;
1667 }
1668 
1669 vector<GammaRecord * > EGS_Ensdf::getGammaRecords() const {
1670  return myGammaRecords;
1671 }
1672 
1673 vector<GammaRecord * > EGS_Ensdf::getMetastableGammaRecords() const {
1674  return myMetastableGammaRecords;
1675 }
1676 
1677 vector<GammaRecord * > EGS_Ensdf::getUncorrelatedGammaRecords() const {
1678  return myUncorrelatedGammaRecords;
1679 }
1680 
1681 vector<AlphaRecord * > EGS_Ensdf::getAlphaRecords() const {
1682  return myAlphaRecords;
1683 }
1684 
1685 Record::Record() {};
1686 Record::Record(vector<string> ensdf) {
1687  if (!ensdf.empty()) {
1688  lines = ensdf;
1689  }
1690 }
1691 
1692 Record::~Record() {
1693 
1694 }
1695 
1696 vector<string> Record::getRecords() const {
1697  return lines;
1698 }
1699 
1700 // Returns the double between two indices, for the first string in the ensdf
1701 // lines array. It is assumed that the characters in this range can be
1702 // converted to a double
1703 double Record::recordToDouble(int startPos, int endPos) {
1704  if (!lines.empty()) {
1705  if (lines.front().length() < startPos) {
1706  egsWarning("Record::recordToDouble: Warning: Record too short to "
1707  "contain desired quantity\n");
1708  return 0;
1709  }
1710  string record = lines.front().substr(startPos-1,
1711  endPos-startPos+1);
1712  return atof(record.c_str());
1713  }
1714  else {
1715  egsWarning("Record::recordToDouble: Error: Record is empty\n");
1716  return 0;
1717  }
1718 }
1719 
1720 // Returns the string between two indices, for the first string in the ensdf
1721 // lines array
1722 string Record::recordToString(int startPos, int endPos) {
1723  if (!lines.empty()) {
1724  if (lines.front().length() < startPos) {
1725  egsWarning("Record::recordToString: Warning: Record too short to "
1726  "contain desired quantity\n");
1727  return "";
1728  }
1729 
1730  return egsTrimString(lines.front().substr(startPos-1, endPos-startPos+1));
1731  }
1732  else {
1733  egsWarning("Record::recordToString: Error: Record is empty\n");
1734  return "";
1735  }
1736 }
1737 
1738 // Searches for a string in the ensdf array of lines for this record
1739 // If the string is found, it returns the content between the end of the
1740 // string and the next space, converting it to a double
1741 //
1742 // For example, for a line like this:
1743 // 64NI2 G KC=1.112E-4 2$LC=1.09E-5 2
1744 // If you set searchString='KC=', the return value would be 1.112E-4
1745 //
1746 // The notAfter string can be used to make sure the searchString is not
1747 // preceeded by the notAfter string. This was necessary to match "PC=" but
1748 // not "IPC="
1749 double Record::getTag(string searchString, string notAfter="") {
1750  if (lines.size() > 1) {
1751 
1752  for (int i=1; i<lines.size(); ++i) {
1753 
1754  int tagPos = lines[i].find(searchString);
1755 
1756  if (tagPos != std::string::npos) {
1757  // Make sure that the string notAfter doesn't occur before
1758  // the search string
1759  size_t notAfterPos = std::string::npos;
1760  if (notAfter.length() > 0 && tagPos-notAfter.length() > 0) {
1761  notAfterPos = lines[i].find(notAfter, tagPos-notAfter.length());
1762  }
1763 
1764  // If the "notAfter" string wasn't found, proceed and get the
1765  // data for the matched tag
1766  if (notAfterPos == std::string::npos || notAfterPos > tagPos) {
1767  tagPos += searchString.length();
1768 
1769  string record = lines[i].substr(tagPos, lines[i].find(" ",tagPos)-tagPos);
1770 
1771  return atof(record.c_str());
1772  }
1773  else {
1774  // The tag we are looking for could still exist, even though
1775  // the first match failed. Look again, this time without
1776  // worrying about the notAfter string because it was already
1777  // found and is assumed to only exist once
1778  tagPos = lines[i].find(searchString, tagPos+searchString.length());
1779 
1780  if (tagPos != std::string::npos) {
1781  tagPos += searchString.length();
1782 
1783  string record = lines[i].substr(tagPos, lines[i].find(" ",tagPos)-tagPos);
1784 
1785  return atof(record.c_str());
1786  }
1787  }
1788  }
1789  }
1790  }
1791  return 0;
1792 }
1793 
1794 // Converts uncertainties in standard format into a double
1795 // For example, 1.23E-4 (67) would be input as value='1.23E-4', and
1796 // stdUncertainty='67'. The return value would be 0.67E-4.
1797 double Record::parseStdUncertainty(string value, string stdUncertainty) {
1798  if (stdUncertainty.length() < 1) {
1799  return 0;
1800  }
1801  if (value.length() < 1) {
1802  egsInformation("Record::parseStdUncertainty: Warning: No uncertainty provided! Returning 0 uncertainty for value of %f\n", value.c_str());
1803  return 0;
1804  }
1805 
1806  if (stdUncertainty.length() > value.length()) {
1807  egsInformation("Record::parseStdUncertainty: Warning: Number of digits in uncertainty greater than number of digits in value. Returning 0 uncertainty for value of %f\n", value.c_str());
1808  return 0;
1809  }
1810 
1811  // Determine if there is an E for scientific notation and where it is
1812  size_t sciNotLoc = value.find_last_of('E');
1813  size_t dotLoc = value.find_last_of('.');
1814 
1815  // Loop backwards through value, starting at the 'E' if there is one
1816  int startPos;
1817  if (sciNotLoc != std::string::npos) {
1818  startPos = sciNotLoc-1;
1819  }
1820  else {
1821  startPos = value.length()-1;
1822  }
1823  int j;
1824  if (stdUncertainty.length() == 2) {
1825  j = 1;
1826  }
1827  else {
1828  j = 0;
1829  }
1830  for (int i = startPos; i >= 0; --i) {
1831  if (i != dotLoc) {
1832  if (j>=0) {
1833  value[i] = stdUncertainty[j--];
1834  }
1835  else {
1836  value[i] = '0';
1837  }
1838  }
1839  }
1840 
1841  return atof(value.c_str());
1842 }
1843 
1844 string Record::getStringAfter(string searchString, size_t len) {
1845  if (lines.size() > 1) {
1846 
1847  for (int i=1; i<lines.size(); ++i) {
1848 
1849  int tagPos = lines[i].find(searchString);
1850 
1851  if (tagPos != std::string::npos) {
1852 
1853  tagPos += searchString.length();
1854 
1855  string record = lines[i].substr(tagPos, len);
1856 
1857  return record;
1858  }
1859  }
1860  }
1861  return "";
1862 }
1863 
1864 // Parse a halflife from a record
1865 // Converts the units to seconds
1866 // Returns the halflife, or a negative number upon failure
1867 double Record::parseHalfLife(int startPos, int endPos) {
1868  if (lines.empty()) {
1869  egsWarning("Record::parseHalfLife: Error: Record is empty\n");
1870  return -5;
1871  }
1872  if (lines.front().length() < startPos) {
1873  egsWarning("Record::parseHalfLife: Warning: Record too short to "
1874  "contain desired quantity\n");
1875  return -5;
1876  }
1877 
1878  string halfLifeStr = egsTrimString(lines.front().substr(startPos-1,
1879  endPos-startPos+1));
1880 
1881  // Return -1 for stable
1882  if (halfLifeStr.substr(0,5).compare("STABLE") == 0) {
1883  return -1;
1884  }
1885 
1886  // Store the length of the numeric part of the string in i
1887  unsigned int numLength;
1888  for (numLength = 0; numLength < halfLifeStr.length(); numLength++) {
1889  if (!isdigit(halfLifeStr[numLength])
1890  && halfLifeStr.at(numLength) != '.') {
1891 
1892  break;
1893  }
1894  }
1895 
1896  // If there was no numeric component return -2
1897  if (halfLifeStr.size() < numLength+2) {
1898  return -2;
1899  }
1900 
1901  // Get the numeric part
1902  double hl = atof(halfLifeStr.substr(0, numLength).c_str());
1903 
1904  // Convert to units of seconds
1905  if (halfLifeStr.size()>numLength+2) {
1906  string units = halfLifeStr.substr(numLength+1, 2);
1907  if (units.compare("Y ") == 0) {
1908  hl *= 31556925.26;
1909  }
1910  else if (units.compare("D ") == 0) {
1911  hl *= 86400;
1912  }
1913  else if (units.compare("H ") == 0) {
1914  hl *= 3600;
1915  }
1916  else if (units.compare("M ") == 0) {
1917  hl *= 60;
1918  }
1919  else if (units.compare("S ") == 0) {
1920  hl *= 1;
1921  }
1922  else if (units.compare("MS") == 0) {
1923  hl *= 1E-3;
1924  }
1925  else if (units.compare("US") == 0) {
1926  hl *= 1E-6;
1927  }
1928  else if (units.compare("NS") == 0) {
1929  hl *= 1E-9;
1930  }
1931  else if (units.compare("PS") == 0) {
1932  hl *= 1E-12;
1933  }
1934  else if (units.compare("FS") == 0) {
1935  hl *= 1E-15;
1936  }
1937  else if (units.compare("AS") == 0) {
1938  hl *= 1E-18;
1939  }
1940  else {
1941  return -3;
1942  }
1943  }
1944  else if (halfLifeStr.size()>numLength+1) {
1945  string units = halfLifeStr.substr(numLength+1, 1);
1946  if (units.compare("Y") == 0) {
1947  hl *= 31556925.26;
1948  }
1949  else if (units.compare("D") == 0) {
1950  hl *= 86400;
1951  }
1952  else if (units.compare("H") == 0) {
1953  hl *= 3600;
1954  }
1955  else if (units.compare("M") == 0) {
1956  hl *= 60;
1957  }
1958  else if (units.compare("S") == 0) {
1959  hl *= 1;
1960  }
1961  else {
1962  return -3;
1963  }
1964  }
1965  else {
1966  hl = -4;
1967  }
1968 
1969  return hl;
1970 }
1971 
1972 // Comment Record
1973 CommentRecord::CommentRecord(vector<string> ensdf):Record(ensdf) {
1974  processEnsdf();
1975 }
1976 
1977 void CommentRecord::processEnsdf() {
1978  if (!lines.empty()) {
1979 
1980  comments = lines;
1981  }
1982 }
1983 
1984 vector<string> CommentRecord::getComments() {
1985  return comments;
1986 }
1987 
1988 // Parent Record
1989 ParentRecord::ParentRecord(vector<string> ensdf):Record(ensdf) {
1990  processEnsdf();
1991 }
1992 
1993 void ParentRecord::processEnsdf() {
1994  halfLife = parseHalfLife(40, 49);
1995 
1996  // Ground state Q-value in keV
1997  // (total energy available for g.s. -> g.s. transition
1998  // It will always be a positive number
1999  // We convert to MeV
2000  Q = recordToDouble(65, 74) / 1000.;
2001 
2002  // If the Q was not contained in the record it returned -1
2003  if (Q == -0.001) {
2004  egsWarning("ParentRecord::processEnsdf: Warning: No Q-value given, any "
2005  "positron records will give errors\n");
2006  Q = 0.;
2007  }
2008 }
2009 
2010 double ParentRecord::getHalfLife() const {
2011  return halfLife;
2012 }
2013 
2014 double ParentRecord::getQ() const {
2015  return Q;
2016 }
2017 
2018 ParentRecord *ParentRecordLeaf::getParentRecord() const {
2019  return getBranch();
2020 }
2021 
2022 ParentRecordLeaf::ParentRecordLeaf(ParentRecord
2023  *myRecord):Leaf<ParentRecord>(myRecord) {
2024 
2025 }
2026 
2027 // Normalization Record
2028 NormalizationRecord::NormalizationRecord(vector<string> ensdf,
2029  ParentRecord *myParent):Record(ensdf), ParentRecordLeaf(myParent) {
2030  processEnsdf();
2031 }
2032 
2033 void NormalizationRecord::processEnsdf() {
2034  normalizeRelative = recordToDouble(10, 19);
2035  normalizeTransition = recordToDouble(22, 29);
2036  normalizeBranch = recordToDouble(32, 39);
2037  normalizeBeta = recordToDouble(42, 49);
2038 
2039  // If the normalization is not specified, it will get initialized to zero
2040  // Change this to 1
2041  if (normalizeRelative < epsilon) {
2042  normalizeRelative = 1;
2043  }
2044  if (normalizeTransition < epsilon) {
2045  normalizeTransition = 1;
2046  }
2047  if (normalizeBranch < epsilon) {
2048  normalizeBranch = 1;
2049  }
2050  if (normalizeBeta < epsilon) {
2051  normalizeBeta = 1;
2052  }
2053 
2054  // Get the daughter element
2055  string element = recordToString(4, 5);
2056 
2057  // Get the Z
2058  Z = setZ(element);
2059 
2060  // Load atomic relaxations
2061  relaxations = new EGS_AtomicRelaxations();
2062  relaxations->loadData(Z);
2063 
2064  // Get the number of shells
2065  nshell = relaxations->getNShell(Z);
2066 
2067  egsInformation("NormalizationRecord::processEnsdf(): Z, nshell: %d %d\n",Z,nshell);
2068 }
2069 
2070 EGS_AtomicRelaxations *NormalizationRecord::getRelaxations() const {
2071  return relaxations;
2072 }
2073 
2074 int NormalizationRecord::getNShell() const {
2075  return nshell;
2076 }
2077 
2078 double NormalizationRecord::getBindingEnergy(int shell) const {
2079  return relaxations->getBindingEnergy(Z,shell);
2080 }
2081 
2082 void NormalizationRecord::relax(int shell,
2083  EGS_Float ecut, EGS_Float pcut,
2084  EGS_RandomGenerator *rndm, double &edep,
2085  EGS_SimpleContainer<EGS_RelaxationParticle> &particles) {
2086  relaxations->relax(Z,shell,ecut,pcut,rndm,edep,particles);
2087 }
2088 
2089 // Multiplier for converting relative photon intensity to photons per 100
2090 // decays in the parent through the decay branch or to photons per 100 neutron
2091 // captures in an (n,gamma) reaction. Required if the absolute photon intensity
2092 // can be calculated
2093 double NormalizationRecord::getRelativeMultiplier() const {
2094  return normalizeRelative;
2095 }
2096 
2097 // Multiplier for convert relative transition intensity (including conversion
2098 // electrons) to transitions per 100 decays of the parent through this decay
2099 // branch or per 100 neutron captures in an (n,gamma) reaction
2100 double NormalizationRecord::getTransitionMultiplier() const {
2101  return normalizeTransition;
2102 }
2103 
2104 // Branching ratio multiplier for converting intensity per 100 decays
2105 // through this decay branch to intensity per 100 decays of the parent nuclide
2106 double NormalizationRecord::getBranchMultiplier() const {
2107  return normalizeBranch;
2108 }
2109 
2110 // Multiplier for converting relative beta- and electron capture intensities to
2111 // intensities per 100 decays through this decay branch. Required if known
2112 double NormalizationRecord::getBetaMultiplier() const {
2113  return normalizeBeta;
2114 }
2115 
2116 NormalizationRecord *NormalizationRecordLeaf::getNormalizationRecord()
2117 const {
2118  return getBranch();
2119 }
2120 
2121 NormalizationRecordLeaf::NormalizationRecordLeaf(NormalizationRecord
2122  *myRecord):Leaf<NormalizationRecord>(myRecord) {
2123 
2124 }
2125 
2126 // Level Record
2127 LevelRecord::LevelRecord() {
2128  energy = 0;
2129  halfLife = 0;
2130  disintegrationIntensity = 0;
2131 }
2132 LevelRecord::LevelRecord(vector<string> ensdf):
2133  Record(ensdf) {
2134  processEnsdf();
2135  disintegrationIntensity = 0;
2136 }
2137 
2138 void LevelRecord::processEnsdf() {
2139  energy = recordToDouble(10, 19) / 1000.; // Convert keV to MeV
2140  halfLife = parseHalfLife(40, 49);
2141 }
2142 
2143 void LevelRecord::setLevelCanDecay(bool canDecayTmp) {
2144  canDecay = canDecayTmp;
2145 }
2146 
2147 bool LevelRecord::levelCanDecay() const {
2148  return canDecay;
2149 }
2150 
2151 void LevelRecord::resetDisintegrationIntensity() {
2152  disintegrationIntensity = 0;
2153 }
2154 
2155 void LevelRecord::cumulDisintegrationIntensity(double disintIntensity) {
2156  disintegrationIntensity += disintIntensity;
2157 }
2158 
2159 double LevelRecord::getDisintegrationIntensity() const {
2160  return disintegrationIntensity;
2161 }
2162 
2163 double LevelRecord::getEnergy() const {
2164  return energy;
2165 }
2166 
2167 double LevelRecord::getHalfLife() const {
2168  return halfLife;
2169 }
2170 
2171 LevelRecord *LevelRecordLeaf::getLevelRecord() const {
2172  return getBranch();
2173 }
2174 
2175 LevelRecordLeaf::LevelRecordLeaf(LevelRecord
2176  *myRecord):Leaf<LevelRecord>(myRecord) {
2177 
2178 }
2179 
2180 // Beta Record
2181 BetaRecordLeaf::BetaRecordLeaf(vector<string> ensdf,
2182  ParentRecord *myParent,
2183  NormalizationRecord *myNormalization,
2184  LevelRecord *myLevel):
2185  Record(ensdf),
2186  ParentRecordLeaf(myParent),
2187  NormalizationRecordLeaf(myNormalization),
2188  LevelRecordLeaf(myLevel) {
2189 
2190  numSampled = 0;
2191 
2192  // Set the Z and atomic weight for the daughter of this decay
2193  string id = egsRemoveWhite(lines.front().substr(0,5));
2194  Z = setZ(id);
2195 
2196  string atomicWeight;
2197  for (unsigned int i=0; i < id.length(); ++i) {
2198  if (!isdigit(id[i])) {
2199  break;
2200  }
2201  else {
2202  atomicWeight.push_back(id[i]);
2203  }
2204  }
2205  A = atoi(atomicWeight.c_str());
2206 
2207  // Get the forbiddenness
2208  if (lines.front().length() > 77) {
2209  string lambda;
2210  lambda.push_back(lines.front().at(77));
2211  forbidden = atoi(lambda.c_str());
2212  }
2213  else {
2214  forbidden = 0;
2215  }
2216 }
2217 int BetaRecordLeaf::getCharge() const {
2218  return q;
2219 }
2220 
2221 void BetaRecordLeaf::incrNumSampled() {
2222  numSampled++;
2223 }
2224 
2225 EGS_I64 BetaRecordLeaf::getNumSampled() const {
2226  return numSampled;
2227 }
2228 
2229 unsigned short int BetaRecordLeaf::getZ() const {
2230  return Z;
2231 }
2232 
2233 unsigned short int BetaRecordLeaf::getAtomicWeight() const {
2234  return A;
2235 }
2236 
2237 unsigned short int BetaRecordLeaf::getForbidden() const {
2238  return forbidden;
2239 }
2240 
2241 void BetaRecordLeaf::setSpectrum(EGS_AliasTable *bspec) {
2242  spectrum = bspec;
2243 }
2244 
2245 EGS_AliasTable *BetaRecordLeaf::getSpectrum() const {
2246  return spectrum;
2247 }
2248 
2249 // Beta- Record
2250 BetaMinusRecord::BetaMinusRecord(vector<string> ensdf,
2251  ParentRecord *myParent,
2252  NormalizationRecord *myNormalization,
2253  LevelRecord *myLevel):
2254  BetaRecordLeaf(ensdf, myParent,
2255  myNormalization, myLevel) {
2256  processEnsdf();
2257  q = -1;
2258  myLevel->cumulDisintegrationIntensity(betaIntensity);
2259 }
2260 
2261 void BetaMinusRecord::processEnsdf() {
2262  finalEnergy = recordToDouble(10, 19) / 1000.; // Convert keV to MeV
2263  betaIntensity = recordToDouble(22, 29);
2264  string betaIntensityStr = recordToString(22, 29);
2265  string betaIntensityUncStr = recordToString(30, 31);
2266 
2267  betaIntensityUnc = parseStdUncertainty(betaIntensityStr, betaIntensityUncStr);
2268  // If the uncertainty is 0 (i.e. not specified), set it to 100%
2269  if (betaIntensityUnc == 0) {
2270  betaIntensityUnc = betaIntensity;
2271  }
2272 
2273  if (getNormalizationRecord()) {
2274  double factor = getNormalizationRecord()->getBetaMultiplier() * getNormalizationRecord()->getBranchMultiplier();
2275 
2276  betaIntensity *= factor;
2277  betaIntensityUnc *= factor;
2278  }
2279 }
2280 
2281 double BetaMinusRecord::getFinalEnergy() const {
2282  return finalEnergy;
2283 }
2284 
2285 double BetaMinusRecord::getBetaIntensity() const {
2286  return betaIntensity;
2287 }
2288 
2289 double BetaMinusRecord::getBetaIntensityUnc() const {
2290  return betaIntensityUnc;
2291 }
2292 
2293 void BetaMinusRecord::setBetaIntensity(double newIntensity) {
2294  betaIntensity = newIntensity;
2295 }
2296 
2297 // Beta+ Record (and Electron Capture)
2298 BetaPlusRecord::BetaPlusRecord(vector<string> ensdf,
2299  ParentRecord *myParent,
2300  NormalizationRecord *myNormalization,
2301  LevelRecord *myLevel):
2302  BetaRecordLeaf(ensdf, myParent,
2303  myNormalization, myLevel) {
2304  processEnsdf();
2305  q = 1;
2306  myLevel->cumulDisintegrationIntensity(betaIntensity);
2307 }
2308 
2309 void BetaPlusRecord::processEnsdf() {
2310  finalEnergy = recordToDouble(10, 19) / 1000.; // Convert keV to MeV
2311  positronIntensity = recordToDouble(22, 29);
2312  string positronIntensityStr = recordToString(22, 29);
2313  string positronIntensityUncStr = recordToString(30, 31);
2314  ecIntensity = recordToDouble(32, 39);
2315  string ecIntensityStr = recordToString(32, 39);
2316  string ecIntensityUncStr = recordToString(40, 41);
2317 
2318  positronIntensityUnc = parseStdUncertainty(positronIntensityStr, positronIntensityUncStr);
2319  // If the uncertainty is 0 (i.e. not specified), set it to 100%
2320  if (positronIntensityUnc == 0) {
2321  positronIntensityUnc = positronIntensity;
2322  }
2323 
2324  ecIntensityUnc = parseStdUncertainty(ecIntensityStr, ecIntensityUncStr);
2325  // If the uncertainty is 0 (i.e. not specified), set it to 100%
2326  if (ecIntensityUnc == 0) {
2327  ecIntensityUnc = ecIntensity;
2328  }
2329 
2330  if (getNormalizationRecord()) {
2331  double factor = getNormalizationRecord()->getBetaMultiplier() * getNormalizationRecord()->getBranchMultiplier();
2332 
2333  positronIntensity *= factor;
2334  ecIntensity *= factor;
2335  positronIntensityUnc *= factor;
2336  ecIntensityUnc *= factor;
2337  }
2338 
2339  // The total intensity for this decay branch
2340  // A decay down this branch will then be split between positron or EC
2341  betaIntensity = positronIntensity + ecIntensity;
2342 
2343  // Re-normalize the positron intensity so that
2344  // positronIntensity + ecIntensity = 1, for easy sampling
2345  // We only set positronIntensity because that's all we need
2346  positronIntensity = positronIntensity / betaIntensity;
2347 
2348  // For positrons we may need to calculate the emission energy
2349  // E = Q - level_energy - 2*mc^2
2350  if (finalEnergy == 0 && positronIntensity > epsilon) {
2351  finalEnergy = getParentRecord()->getQ()
2352  - getLevelRecord()->getEnergy() - 1.022;
2353 
2354  if (finalEnergy < 0.) {
2355  egsWarning("BetaPlusRecord::processEnsdf: Error: Final energy of "
2356  "positron could not be calculated. Setting energy to zero!\n"
2357  );
2358  finalEnergy = 0.;
2359  }
2360  }
2361 
2362  if (ecIntensity > 0 && getNormalizationRecord()) {
2363  // Get the number of shells
2364  int nshell = getNormalizationRecord()->getNShell();
2365 
2366  double icK = getTag("CK=");
2367  double icL = getTag("CL=");
2368  double icM = getTag("CM=");
2369  double icN = getTag("CN=");
2370  double icO = getTag("CO=");
2371  double icP = getTag("CP=");
2372  double icQ = getTag("CQ=");
2373 
2374  // The K shell
2375  ecShellIntensity.push_back(icK);
2376 
2377  // The L1, L2, L3 shells
2378  // TODO: Equal probability is assigned to subshells!
2379  // This assumption is an approximation
2380  int numShellsToInclude = min(4,nshell);
2381  for (unsigned int i=1; i<numShellsToInclude; ++i) {
2382  ecShellIntensity.push_back(ecShellIntensity.back() + icL/(numShellsToInclude-1));
2383  }
2384  // Count number of shells as we go
2385  // Once we hit the number of shells for this element, return
2386  if (numShellsToInclude < 4) {
2387 // for (int i=0; i<ecShellIntensity.size(); ++i) {
2388 // egsInformation("BetaPlusRecord::processEnsdf: Shell %d: P=%f\n",i,ecShellIntensity[i]);
2389 // }
2390  return;
2391  }
2392 
2393  // The M1-M5 shells
2394  numShellsToInclude = min(9,nshell);
2395  for (unsigned int i=4; i<numShellsToInclude; ++i) {
2396  ecShellIntensity.push_back(ecShellIntensity.back() + icM/(numShellsToInclude-4));
2397  }
2398  if (numShellsToInclude < 9) {
2399 // for (int i=0; i<ecShellIntensity.size(); ++i) {
2400 // egsInformation("BetaPlusRecord::processEnsdf: Shell %d: P=%f\n",i,ecShellIntensity[i]);
2401 // }
2402  return;
2403  }
2404 
2405  // The N1-7 shells
2406  numShellsToInclude = min(16,nshell);
2407  for (unsigned int i=9; i<numShellsToInclude; ++i) {
2408  ecShellIntensity.push_back(ecShellIntensity.back() + icN/(numShellsToInclude-9));
2409  }
2410  if (numShellsToInclude < 16) {
2411 // for (int i=0; i<ecShellIntensity.size(); ++i) {
2412 // egsInformation("BetaPlusRecord::processEnsdf: Shell %d: P=%f\n",i,ecShellIntensity[i]);
2413 // }
2414  return;
2415  }
2416 
2417  // The O1-7 shells
2418  numShellsToInclude = min(23,nshell);
2419  for (unsigned int i=16; i<numShellsToInclude; ++i) {
2420  ecShellIntensity.push_back(ecShellIntensity.back() + icO/(numShellsToInclude-16));
2421  }
2422 
2423  if (numShellsToInclude < 23) {
2424 // for (int i=0; i<ecShellIntensity.size(); ++i) {
2425 // egsInformation("BetaPlusRecord::processEnsdf: Shell %d: P=%f\n",i,ecShellIntensity[i]);
2426 // }
2427  return;
2428  }
2429 
2430  // The P1-3 shells
2431  numShellsToInclude = min(26,nshell);
2432  for (unsigned int i=23; i<numShellsToInclude; ++i) {
2433  ecShellIntensity.push_back(ecShellIntensity.back() + icP/(numShellsToInclude-23));
2434  }
2435 
2436  if (numShellsToInclude < 26) {
2437 // for (int i=0; i<ecShellIntensity.size(); ++i) {
2438 // egsInformation("BetaPlusRecord::processEnsdf: Shell %d: P=%f\n",i,ecShellIntensity[i]);
2439 // }
2440  return;
2441  }
2442 
2443  // The Q1 shell
2444  numShellsToInclude = 27;
2445  ecShellIntensity.push_back(ecShellIntensity.back() + icQ/(numShellsToInclude-26));
2446 
2447 // for (int i=0; i<ecShellIntensity.size(); ++i) {
2448 // egsInformation("BetaPlusRecord::processEnsdf: Shell %d: P=%f\n",i,ecShellIntensity[i]);
2449 // }
2450  return;
2451  }
2452 }
2453 
2454 void BetaPlusRecord::relax(int shell,
2455  EGS_Float ecut, EGS_Float pcut,
2456  EGS_RandomGenerator *rndm, double &edep,
2457  EGS_SimpleContainer<EGS_RelaxationParticle> &particles) {
2458  getNormalizationRecord()->relax(shell,ecut,pcut,rndm,edep,particles);
2459 }
2460 
2461 double BetaPlusRecord::getFinalEnergy() const {
2462  return finalEnergy;
2463 }
2464 
2465 double BetaPlusRecord::getBetaIntensity() const {
2466  return betaIntensity;
2467 }
2468 
2469 double BetaPlusRecord::getPositronIntensity() const {
2470  return positronIntensity;
2471 }
2472 
2473 double BetaPlusRecord::getPositronIntensityUnc() const {
2474  return positronIntensityUnc;
2475 }
2476 
2477 double BetaPlusRecord::getECIntensityUnc() const {
2478  return ecIntensityUnc;
2479 }
2480 
2481 void BetaPlusRecord::setBetaIntensity(double newIntensity) {
2482  betaIntensity = newIntensity;
2483 }
2484 
2485 void BetaPlusRecord::setPositronIntensity(double newIntensity) {
2486  positronIntensity = newIntensity;
2487 }
2488 
2489 // Gamma Record
2490 GammaRecord::GammaRecord(vector<string> ensdf,
2491  ParentRecord *myParent,
2492  NormalizationRecord *myNormalization,
2493  LevelRecord *myLevel):
2494  Record(ensdf),
2495  ParentRecordLeaf(myParent),
2496  NormalizationRecordLeaf(myNormalization),
2497  LevelRecordLeaf(myLevel) {
2498  processEnsdf();
2499  q = 0;
2500  numGammaSampled = 0;
2501  numICSampled = 0;
2502  numIPSampled = 0;
2503  multipleTransitionProb = 0;
2504 }
2505 
2506 GammaRecord::GammaRecord(GammaRecord *gamma):
2507  Record(),
2508  ParentRecordLeaf(gamma->getParentRecord()),
2509  NormalizationRecordLeaf(gamma->getNormalizationRecord()),
2510  LevelRecordLeaf(gamma->getLevelRecord()) {
2511 
2512  numGammaSampled = gamma->numGammaSampled;
2513  numICSampled = gamma->numICSampled;
2514  numIPSampled = gamma->numIPSampled;
2515  decayEnergy = gamma->decayEnergy;
2516  transitionIntensity = gamma->transitionIntensity;
2517  multipleTransitionProb = gamma->multipleTransitionProb;
2518  gammaIntensity = gamma->gammaIntensity;
2519  gammaIntensityUnc = gamma->gammaIntensityUnc;
2520  icCoeff = gamma->icCoeff;
2521  icCoeffUnc = gamma->icCoeffUnc;
2522  ipCoeff = gamma->ipCoeff;
2523  ipCoeffUnc = gamma->ipCoeffUnc;
2524  q = gamma->q;
2525  finalLevel = gamma->finalLevel;
2526 }
2527 
2528 void GammaRecord::processEnsdf() {
2529  decayEnergy = recordToDouble(10, 19) / 1000.; // Convert keV to MeV
2530  gammaIntensity = recordToDouble(22, 29);
2531  string gammaIntensityStr = recordToString(22, 29);
2532  string gammaIntensityUncStr = recordToString(30, 31);
2533 
2534  // Get the internal conversion coefficient
2535  // Note: This might not equal the sum of the individual shell intensities,
2536  // because some of the shell intensities might not be known
2537  icCoeff = recordToDouble(56, 62);
2538  string icCoeffStr = recordToString(56, 62);
2539  string icCoeffUncStr = recordToString(63, 64);
2540 
2541  // If we don't find the gamma intensity, check for the first RI=
2542  if (gammaIntensity < epsilon) {
2543  gammaIntensity = getTag("RI=");
2544  }
2545 
2546  icCoeffUnc = parseStdUncertainty(icCoeffStr, icCoeffUncStr);
2547  // If the uncertainty is 0 (i.e. not specified), set it to 100%
2548  if (icCoeffUnc == 0) {
2549  icCoeffUnc = icCoeff;
2550  }
2551 
2552  // Check for internal pair production
2553  // The value and uncertainty are stored in 11 characters after IPC=
2554  string ipCoeffStr_tmp = getStringAfter("IPC=", 11);
2555  if (ipCoeffStr_tmp.length() > 0) {
2556  string ipCoeffStr = ipCoeffStr_tmp.substr(0, 9);
2557  string ipCoeffUncStr = ipCoeffStr_tmp.substr(9, 2);
2558  ipCoeff = atof(ipCoeffStr.c_str());
2559  ipCoeffUnc = atof(ipCoeffUncStr.c_str());
2560 
2561  if (ipCoeffUnc == 0) {
2562  ipCoeffUnc = ipCoeff;
2563  }
2564  }
2565  else {
2566  ipCoeff = 0;
2567  ipCoeffUnc = 0;
2568  }
2569 
2570  // Get the transition intensity instead if gamma still zero
2571  if (gammaIntensity < epsilon) {
2572  double ti = getTag("TI ");
2573  // Calculate the gamma intensity from it
2574  gammaIntensity = ti / ((1+icCoeff) * (1+ipCoeff));
2575  }
2576 
2577  // Set the uncertainty on the gamma intensity
2578  gammaIntensityUnc = parseStdUncertainty(gammaIntensityStr, gammaIntensityUncStr);
2579  // If the uncertainty is 0 (i.e. not specified), set it to 100%
2580  if (gammaIntensityUnc == 0) {
2581  gammaIntensityUnc = gammaIntensity;
2582  }
2583 
2584  if (getNormalizationRecord()) {
2585  double factor = getNormalizationRecord()->getRelativeMultiplier() *
2586  getNormalizationRecord()->getBranchMultiplier();
2587 
2588  gammaIntensity *= factor;
2589  gammaIntensityUnc *= factor;
2590  }
2591 
2592  // Calculate the total transition intensity
2593  transitionIntensity = gammaIntensity * (1+icCoeff) * (1+ipCoeff);
2594 
2595  if (icCoeff > 0 && getNormalizationRecord()) {
2596  // Get the number of shells
2597  int nshell = getNormalizationRecord()->getNShell();
2598 
2599  double icK = getTag("KC=");
2600  double icL = getTag("LC=");
2601  double icM = getTag("MC=");
2602  double icN = getTag("NC=");
2603  double icO = getTag("OC=");
2604  double icP = getTag("PC=", "I"); // Make sure "IPC" doesn't get matched
2605  double icQ = getTag("QC=");
2606 
2607  // The K shell
2608  icIntensity.push_back(icK / icCoeff);
2609 
2610  // The L1, L2, L3 shells
2611  // TODO: Equal probability is assigned to subshells!
2612  // This assumption is an approximation
2613  int numShellsToInclude = min(4,nshell);
2614  for (unsigned int i=1; i<numShellsToInclude; ++i) {
2615  icIntensity.push_back(icIntensity.back() + (icL / icCoeff)/(numShellsToInclude-1));
2616  }
2617  // Count number of shells as we go
2618  // Once we hit the number of shells for this element, return
2619  if (numShellsToInclude < 4) {
2620 // for (int i=0; i<icIntensity.size(); ++i) {
2621 // egsInformation("GammaRecord::processEnsdf: Shell %d: P=%f\n",i,icIntensity[i]);
2622 // }
2623  return;
2624  }
2625 
2626  // The M1-M5 shells
2627  numShellsToInclude = min(9,nshell);
2628  for (unsigned int i=4; i<numShellsToInclude; ++i) {
2629  icIntensity.push_back(icIntensity.back() + (icM / icCoeff)/(numShellsToInclude-4));
2630  }
2631  if (numShellsToInclude < 9) {
2632 // for (int i=0; i<icIntensity.size(); ++i) {
2633 // egsInformation("GammaRecord::processEnsdf: Shell %d: P=%f\n",i,icIntensity[i]);
2634 // }
2635  return;
2636  }
2637 
2638  // The N1-7 shells
2639  numShellsToInclude = min(16,nshell);
2640  for (unsigned int i=9; i<numShellsToInclude; ++i) {
2641  icIntensity.push_back(icIntensity.back() + (icN / icCoeff)/(numShellsToInclude-9));
2642  }
2643  if (numShellsToInclude < 16) {
2644 // for (int i=0; i<icIntensity.size(); ++i) {
2645 // egsInformation("GammaRecord::processEnsdf: Shell %d: P=%f\n",i,icIntensity[i]);
2646 // }
2647  return;
2648  }
2649 
2650  // The O1-7 shells
2651  numShellsToInclude = min(23,nshell);
2652  for (unsigned int i=16; i<numShellsToInclude; ++i) {
2653  icIntensity.push_back(icIntensity.back() + (icO / icCoeff)/(numShellsToInclude-16));
2654  }
2655 
2656  if (numShellsToInclude < 23) {
2657 // for (int i=0; i<icIntensity.size(); ++i) {
2658 // egsInformation("GammaRecord::processEnsdf: Shell %d: P=%f\n",i,icIntensity[i]);
2659 // }
2660  return;
2661  }
2662 
2663  // The P1-3 shells
2664  numShellsToInclude = min(26,nshell);
2665  for (unsigned int i=23; i<numShellsToInclude; ++i) {
2666  icIntensity.push_back(icIntensity.back() + (icP / icCoeff)/(numShellsToInclude-23));
2667  }
2668 
2669  if (numShellsToInclude < 26) {
2670 // for (int i=0; i<icIntensity.size(); ++i) {
2671 // egsInformation("GammaRecord::processEnsdf: Shell %d: P=%f\n",i,icIntensity[i]);
2672 // }
2673  return;
2674  }
2675 
2676  // The Q1 shell
2677  numShellsToInclude = 27;
2678  icIntensity.push_back(icIntensity.back() + (icQ / icCoeff)/(numShellsToInclude-26));
2679 
2680 // for (int i=0; i<icIntensity.size(); ++i) {
2681 // egsInformation("GammaRecord::processEnsdf: Shell %d: P=%f\n",i,icIntensity[i]);
2682 // }
2683  return;
2684 
2685  }
2686 }
2687 
2688 double GammaRecord::getBindingEnergy(int shell) const {
2689  return getNormalizationRecord()->getBindingEnergy(shell);
2690 }
2691 
2692 void GammaRecord::relax(int shell,
2693  EGS_Float ecut, EGS_Float pcut,
2694  EGS_RandomGenerator *rndm, double &edep,
2695  EGS_SimpleContainer<EGS_RelaxationParticle> &particles) {
2696  getNormalizationRecord()->relax(shell,ecut,pcut,rndm,edep,particles);
2697 }
2698 
2699 double GammaRecord::getDecayEnergy() const {
2700  return decayEnergy;
2701 }
2702 
2703 double GammaRecord::getMultiTransitionProb() const {
2704  return multipleTransitionProb;
2705 }
2706 
2707 void GammaRecord::setMultiTransitionProb(double newIntensity) {
2708  multipleTransitionProb = newIntensity;
2709 }
2710 
2711 double GammaRecord::getTransitionIntensity() const {
2712  return transitionIntensity;
2713 }
2714 
2715 double GammaRecord::getGammaIntensity() const {
2716  return gammaIntensity;
2717 }
2718 
2719 double GammaRecord::getGammaIntensityUnc() const {
2720  return gammaIntensityUnc;
2721 }
2722 
2723 double GammaRecord::getICIntensity() const {
2724  return icCoeff;
2725 }
2726 
2727 double GammaRecord::getICIntensityUnc() const {
2728  return icCoeffUnc;
2729 }
2730 
2731 double GammaRecord::getIPIntensity() const {
2732  return ipCoeff;
2733 }
2734 
2735 double GammaRecord::getIPIntensityUnc() const {
2736  return ipCoeffUnc;
2737 }
2738 
2739 void GammaRecord::setTransitionIntensity(double newIntensity) {
2740  transitionIntensity = newIntensity;
2741 }
2742 
2743 void GammaRecord::setGammaIntensity(double newIntensity) {
2744  gammaIntensity = newIntensity;
2745 }
2746 
2747 void GammaRecord::setICIntensity(double newIntensity) {
2748  icCoeff = newIntensity;
2749 }
2750 
2751 int GammaRecord::getCharge() const {
2752  return q;
2753 }
2754 
2755 void GammaRecord::incrGammaSampled() {
2756  numGammaSampled++;
2757 }
2758 
2759 void GammaRecord::incrICSampled() {
2760  numICSampled++;
2761 }
2762 
2763 void GammaRecord::incrIPSampled() {
2764  numIPSampled++;
2765 }
2766 
2767 EGS_I64 GammaRecord::getGammaSampled() const {
2768  return numGammaSampled;
2769 }
2770 
2771 EGS_I64 GammaRecord::getICSampled() const {
2772  return numICSampled;
2773 }
2774 
2775 EGS_I64 GammaRecord::getIPSampled() const {
2776  return numIPSampled;
2777 }
2778 
2779 LevelRecord *GammaRecord::getFinalLevel() const {
2780  return finalLevel;
2781 }
2782 
2783 void GammaRecord::setFinalLevel(LevelRecord *newLevel) {
2784  finalLevel = newLevel;
2785 }
2786 
2787 // Alpha Record
2788 AlphaRecord::AlphaRecord(vector<string> ensdf,
2789  ParentRecord *myParent,
2790  NormalizationRecord *myNormalization,
2791  LevelRecord *myLevel):
2792  Record(ensdf),
2793  ParentRecordLeaf(myParent), NormalizationRecordLeaf(myNormalization),
2794  LevelRecordLeaf(myLevel) {
2795 
2796  processEnsdf();
2797  q = 2;
2798  numSampled = 0;
2799  myLevel->cumulDisintegrationIntensity(alphaIntensity);
2800 }
2801 
2802 void AlphaRecord::processEnsdf() {
2803  finalEnergy = recordToDouble(10, 19) / 1000.; // Convert keV to MeV
2804  alphaIntensity = recordToDouble(22, 29);
2805 
2806  string alphaIntensityStr = recordToString(22, 29);
2807  string alphaIntensityUncStr = recordToString(30, 31);
2808 
2809  alphaIntensityUnc = parseStdUncertainty(alphaIntensityStr, alphaIntensityUncStr);
2810  // If the uncertainty is 0 (i.e. not specified), set it to 100%
2811  if (alphaIntensityUnc == 0) {
2812  alphaIntensityUnc = alphaIntensity;
2813  }
2814 
2815  if (getNormalizationRecord()) {
2816  alphaIntensity *= getNormalizationRecord()->getBranchMultiplier();
2817  alphaIntensityUnc *= getNormalizationRecord()->getBranchMultiplier();
2818  }
2819 }
2820 
2821 double AlphaRecord::getFinalEnergy() const {
2822  return finalEnergy;
2823 }
2824 
2825 double AlphaRecord::getAlphaIntensity() const {
2826  return alphaIntensity;
2827 }
2828 
2829 double AlphaRecord::getAlphaIntensityUnc() const {
2830  return alphaIntensityUnc;
2831 }
2832 
2833 void AlphaRecord::setAlphaIntensity(double newIntensity) {
2834  alphaIntensity = newIntensity;
2835 }
2836 
2837 int AlphaRecord::getCharge() const {
2838  return q;
2839 }
2840 
2841 void AlphaRecord::incrNumSampled() {
2842  numSampled++;
2843 }
2844 
2845 EGS_I64 AlphaRecord::getNumSampled() const {
2846  return numSampled;
2847 }
2848 
EGS_AliasTable
A class for sampling random values from a given probability distribution using the alias table techni...
Definition: egs_alias_table.h:64
EGS_AtomicRelaxations
Definition: egs_atomic_relaxations.h:61
EGS_AtomicRelaxations::getBindingEnergy
EGS_Float getBindingEnergy(int Z, int shell)
Definition: egs_atomic_relaxations.cpp:529
EGS_AtomicRelaxations::getNShell
int getNShell(int Z)
Definition: egs_atomic_relaxations.cpp:533
EGS_AtomicRelaxations::loadData
int loadData(int Z)
Definition: egs_atomic_relaxations.cpp:511
EGS_AtomicRelaxations::relax
void relax(int Z, int sh, EGS_Float ecut, EGS_Float pcut, EGS_RandomGenerator *rndm, double &edep, EGS_SimpleContainer< EGS_RelaxationParticle > &particles)
Definition: egs_atomic_relaxations.cpp:523
EGS_Ensdf::~EGS_Ensdf
~EGS_Ensdf()
Destructor.
Definition: egs_ensdf.cpp:238
EGS_Ensdf::EGS_Ensdf
EGS_Ensdf(const string nuclide, const string ensdf_filename="", const string relaxType="eadl", const bool allowMultiTrans=false, int verbosity=1)
Construct an ensdf object.
Definition: egs_ensdf.cpp:196
EGS_RandomGenerator
Base random number generator class. All random number generators should be derived from this class.
Definition: egs_rndm.h:67
EGS_SimpleContainer< EGS_RelaxationParticle >
egs_ensdf.h
The ensdf library header file.
egsInformation
EGS_InfoFunction EGS_EXPORT egsInformation
Always use this function for reporting the progress of a simulation and any other type of information...
Definition: egs_functions.cpp:135
epsilon
const EGS_Float epsilon
The epsilon constant for floating point comparisons.
Definition: egs_functions.h:62
egsWarning
EGS_InfoFunction EGS_EXPORT egsWarning
Always use this function for reporting warnings.
Definition: egs_functions.cpp:136
Generated by   doxygen 1.9.1