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  string line = (*comment)->getComment();
1432 
1433  // Search for this line to ensure emissions will follow the right format
1434  if (line.find("{U Energy (keV)} {U Intensity} {U Line}") != std::string::npos) {
1435  containsEmissions = true;
1436  }
1437 
1438  if (containsEmissions) {
1439  // Check for the end of multi-line records
1440  // and average them together
1441  if (line.length() < 48 ||
1442  ((xrayContinues || augerContinues) && line.at(30) != '|')) {
1443 
1444  // If we just finished going through a series of
1445  // lines that started with a "total" line at the top
1446  // then we'll check to make sure they have intensities assigned
1447  if (gotTotal) {
1448 
1449  // In the event that a zero intensity is in one of the lines
1450  // following the "total" line, ALL of those following
1451  // lines are assigned an equal fraction of the total
1452  // intensity. This is an imperfect work-around for insufficient
1453  // data. Using this method, the correct energies are used,
1454  // rather than assigning a single averaged "total" energy line.
1455  if (countNumAfterTotal > 0) {
1456  // X-rays
1457  if (lineTotalType == 0) {
1458  bool containsZeroIntensity = false;
1459  for (std::vector<double>::iterator it = xrayIntensities.end()-countNumAfterTotal; it != xrayIntensities.end(); ++it) {
1460  if (*it < epsilon) {
1461  containsZeroIntensity = true;
1462  break;
1463  }
1464  }
1465 
1466  if (containsZeroIntensity) {
1467  for (std::vector<double>::iterator it = xrayIntensities.end()-countNumAfterTotal; it != xrayIntensities.end(); ++it) {
1468  *it = lineTotalIntensity / countNumAfterTotal;
1469  }
1470  }
1471 
1472  // Auger
1473  }
1474  else if (lineTotalType == -1) {
1475  bool containsZeroIntensity = false;
1476  for (std::vector<double>::iterator it = augerIntensities.end()-countNumAfterTotal; it != augerIntensities.end(); ++it) {
1477  if (*it < epsilon) {
1478  containsZeroIntensity = true;
1479  break;
1480  }
1481  }
1482 
1483  if (containsZeroIntensity) {
1484  for (std::vector<double>::iterator it = augerIntensities.end()-countNumAfterTotal; it != augerIntensities.end(); ++it) {
1485  *it = lineTotalIntensity / countNumAfterTotal;
1486  }
1487  }
1488  }
1489  }
1490 
1491  gotTotal = false;
1492  countNumAfterTotal = 0;
1493  lineTotalIntensity = 0.;
1494  }
1495 
1496  if ((xrayContinues || augerContinues)
1497  && multilineEnergies.size() > 0) {
1498 
1499  double energySum = 0;
1500  double intensitySum = 0;
1501  unsigned int numNonzeroE = 0;
1502  unsigned int numNonzeroI = 0;
1503  for (unsigned int i=0; i < multilineEnergies.size(); ++i) {
1504  if (multilineEnergies[i] > 0) {
1505  energySum += multilineEnergies[i];
1506  numNonzeroE++;
1507  }
1508  }
1509  for (unsigned int i=0; i < multilineIntensities.size(); ++i) {
1510  if (multilineIntensities[i] > epsilon) {
1511  intensitySum += multilineIntensities[i];
1512  numNonzeroI++;
1513  }
1514  }
1515  double energy;
1516  if (numNonzeroE > 0) {
1517  energy = energySum / numNonzeroE;
1518  }
1519  double intensity;
1520  if (numNonzeroI > 0) {
1521  intensity = intensitySum / numNonzeroI;
1522  }
1523 
1524  if (numNonzeroE > 0 && numNonzeroI > 0) {
1525  if (xrayContinues) {
1526  xrayEnergies.push_back(energy);
1527  xrayIntensities.push_back(intensity);
1528  }
1529  else {
1530  augerEnergies.push_back(energy);
1531  augerIntensities.push_back(intensity);
1532  }
1533  }
1534 
1535  multilineEnergies.clear();
1536  multilineIntensities.clear();
1537  }
1538 
1539  xrayContinues = false;
1540  augerContinues = false;
1541  }
1542 
1543  // Check for records containing XRays or Auger electrons
1544  if (line.length() > 48) {
1545 
1546  string emissionLine = egsTrimString(line.substr(47));
1547 
1548  // See if the line is an XRay or Auger
1549  if (emissionLine.length() < 1 || (emissionLine.at(0) != 'X' && emissionLine.find("AUGER") == std::string::npos)) {
1550  continue;
1551  }
1552 
1553  string eStr = egsTrimString(line.substr(13, 15));
1554 
1555  // If we have a range in energy (e.g. 0.1-0.3)
1556  // Find the average
1557  size_t eDash = eStr.find('-');
1558  double energy;
1559  if (eDash!=std::string::npos) {
1560  if (eStr.length() > eDash+1) {
1561  double e1 = atof(eStr.substr(0, eDash).c_str());
1562  double e2 = atof(eStr.substr(eDash+1).c_str());
1563  energy = (e1 + e2) / 2;
1564  }
1565  else {
1566  energy = atof(eStr.substr(0, eDash).c_str());
1567  }
1568  }
1569  else {
1570  energy = atof(eStr.c_str());
1571  }
1572 
1573  // Convert the energy from keV to MeV
1574  energy /= 1000.;
1575 
1576  // Get the intensity
1577  string iStr = egsTrimString(line.substr(32, 9));
1578  double intensity = atof(iStr.c_str());
1579 
1580  // If this is a line coming after a "total" line,
1581  // increment a counter. This will be used in the
1582  // event that the lines following the "total"
1583  // have zero intensity assigned
1584  if (gotTotal && energy > epsilon) {
1585  countNumAfterTotal++;
1586  }
1587 
1588  // If this line is the total of the next lines, we will
1589  // skip this line and use the individual ones
1590  // However, record the total intensity in case we need it
1591  if (emissionLine.find("(total)") != std::string::npos) {
1592  gotTotal = true;
1593  lineTotalIntensity = intensity;
1594  if (emissionLine.find("AUGER") != std::string::npos) {
1595  lineTotalType = -1;
1596  }
1597  else {
1598  lineTotalType = 0;
1599  }
1600  continue;
1601  }
1602 
1603  // Multi-line records have a bar '|' at 30
1604  // We will store the data and average them later
1605  if (line.at(30) == '|') {
1606  if (emissionLine.at(0) == 'X') {
1607  xrayContinues = true;
1608  }
1609  else if (emissionLine.find("AUGER") != std::string::npos) {
1610  augerContinues = true;
1611  }
1612 
1613  multilineEnergies.push_back(energy);
1614  multilineIntensities.push_back(intensity);
1615 
1616  }
1617  else {
1618  if (emissionLine.at(0) == 'X') {
1619  if ((energy > epsilon && intensity > epsilon) ||
1620  (gotTotal && energy > epsilon)) {
1621  xrayEnergies.push_back(energy);
1622  xrayIntensities.push_back(intensity);
1623  }
1624  }
1625  else if (emissionLine.find("AUGER") != std::string::npos) {
1626  if ((energy > epsilon && intensity > epsilon) ||
1627  (gotTotal && energy > epsilon)) {
1628  augerEnergies.push_back(energy);
1629  augerIntensities.push_back(intensity);
1630  }
1631  }
1632  }
1633  }
1634  }
1635  }
1636 }
1637 
1638 vector<double > EGS_Ensdf::getXRayIntensities() const {
1639  return xrayIntensities;
1640 }
1641 
1642 vector<double > EGS_Ensdf::getXRayEnergies() const {
1643  return xrayEnergies;
1644 }
1645 
1646 vector<double > EGS_Ensdf::getAugerIntensities() const {
1647  return augerIntensities;
1648 }
1649 
1650 vector<double > EGS_Ensdf::getAugerEnergies() const {
1651  return augerEnergies;
1652 }
1653 
1654 vector<ParentRecord * > EGS_Ensdf::getParentRecords() const {
1655  return myParentRecords;
1656 }
1657 
1658 vector<LevelRecord * > EGS_Ensdf::getLevelRecords() const {
1659  return myLevelRecords;
1660 }
1661 
1662 vector<BetaRecordLeaf * > EGS_Ensdf::getBetaRecords() const {
1663  return myBetaRecords;
1664 }
1665 
1666 vector<GammaRecord * > EGS_Ensdf::getGammaRecords() const {
1667  return myGammaRecords;
1668 }
1669 
1670 vector<GammaRecord * > EGS_Ensdf::getMetastableGammaRecords() const {
1671  return myMetastableGammaRecords;
1672 }
1673 
1674 vector<GammaRecord * > EGS_Ensdf::getUncorrelatedGammaRecords() const {
1675  return myUncorrelatedGammaRecords;
1676 }
1677 
1678 vector<AlphaRecord * > EGS_Ensdf::getAlphaRecords() const {
1679  return myAlphaRecords;
1680 }
1681 
1682 Record::Record() {};
1683 Record::Record(vector<string> ensdf) {
1684  if (!ensdf.empty()) {
1685  lines = ensdf;
1686  }
1687 }
1688 
1689 Record::~Record() {
1690 
1691 }
1692 
1693 vector<string> Record::getRecords() const {
1694  return lines;
1695 }
1696 
1697 // Returns the double between two indices, for the first string in the ensdf
1698 // lines array. It is assumed that the characters in this range can be
1699 // converted to a double
1700 double Record::recordToDouble(int startPos, int endPos) {
1701  if (!lines.empty()) {
1702  if (lines.front().length() < startPos) {
1703  egsWarning("Record::recordToDouble: Warning: Record too short to "
1704  "contain desired quantity\n");
1705  return 0;
1706  }
1707  string record = lines.front().substr(startPos-1,
1708  endPos-startPos+1);
1709  return atof(record.c_str());
1710  }
1711  else {
1712  egsWarning("Record::recordToDouble: Error: Record is empty\n");
1713  return 0;
1714  }
1715 }
1716 
1717 // Returns the string between two indices, for the first string in the ensdf
1718 // lines array
1719 string Record::recordToString(int startPos, int endPos) {
1720  if (!lines.empty()) {
1721  if (lines.front().length() < startPos) {
1722  egsWarning("Record::recordToString: Warning: Record too short to "
1723  "contain desired quantity\n");
1724  return "";
1725  }
1726 
1727  return egsTrimString(lines.front().substr(startPos-1, endPos-startPos+1));
1728  }
1729  else {
1730  egsWarning("Record::recordToString: Error: Record is empty\n");
1731  return "";
1732  }
1733 }
1734 
1735 // Searches for a string in the ensdf array of lines for this record
1736 // If the string is found, it returns the content between the end of the
1737 // string and the next space, converting it to a double
1738 //
1739 // For example, for a line like this:
1740 // 64NI2 G KC=1.112E-4 2$LC=1.09E-5 2
1741 // If you set searchString='KC=', the return value would be 1.112E-4
1742 //
1743 // The notAfter string can be used to make sure the searchString is not
1744 // preceeded by the notAfter string. This was necessary to match "PC=" but
1745 // not "IPC="
1746 double Record::getTag(string searchString, string notAfter="") {
1747  if (lines.size() > 1) {
1748 
1749  for (int i=1; i<lines.size(); ++i) {
1750 
1751  int tagPos = lines[i].find(searchString);
1752 
1753  if (tagPos != std::string::npos) {
1754  // Make sure that the string notAfter doesn't occur before
1755  // the search string
1756  size_t notAfterPos = std::string::npos;
1757  if (notAfter.length() > 0 && tagPos-notAfter.length() > 0) {
1758  notAfterPos = lines[i].find(notAfter, tagPos-notAfter.length());
1759  }
1760 
1761  // If the "notAfter" string wasn't found, proceed and get the
1762  // data for the matched tag
1763  if (notAfterPos == std::string::npos || notAfterPos > tagPos) {
1764  tagPos += searchString.length();
1765 
1766  string record = lines[i].substr(tagPos, lines[i].find(" ",tagPos)-tagPos);
1767 
1768  return atof(record.c_str());
1769  }
1770  else {
1771  // The tag we are looking for could still exist, even though
1772  // the first match failed. Look again, this time without
1773  // worrying about the notAfter string because it was already
1774  // found and is assumed to only exist once
1775  tagPos = lines[i].find(searchString, tagPos+searchString.length());
1776 
1777  if (tagPos != std::string::npos) {
1778  tagPos += searchString.length();
1779 
1780  string record = lines[i].substr(tagPos, lines[i].find(" ",tagPos)-tagPos);
1781 
1782  return atof(record.c_str());
1783  }
1784  }
1785  }
1786  }
1787  }
1788  return 0;
1789 }
1790 
1791 // Converts uncertainties in standard format into a double
1792 // For example, 1.23E-4 (67) would be input as value='1.23E-4', and
1793 // stdUncertainty='67'. The return value would be 0.67E-4.
1794 double Record::parseStdUncertainty(string value, string stdUncertainty) {
1795  if (stdUncertainty.length() < 1) {
1796  return 0;
1797  }
1798  if (value.length() < 1) {
1799  egsInformation("Record::parseStdUncertainty: Warning: No uncertainty provided! Returning 0 uncertainty for value of %f\n", value.c_str());
1800  return 0;
1801  }
1802 
1803  if (stdUncertainty.length() > value.length()) {
1804  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());
1805  return 0;
1806  }
1807 
1808  // Determine if there is an E for scientific notation and where it is
1809  size_t sciNotLoc = value.find_last_of('E');
1810  size_t dotLoc = value.find_last_of('.');
1811 
1812  // Loop backwards through value, starting at the 'E' if there is one
1813  int startPos;
1814  if (sciNotLoc != std::string::npos) {
1815  startPos = sciNotLoc-1;
1816  }
1817  else {
1818  startPos = value.length()-1;
1819  }
1820  int j;
1821  if (stdUncertainty.length() == 2) {
1822  j = 1;
1823  }
1824  else {
1825  j = 0;
1826  }
1827  for (int i = startPos; i >= 0; --i) {
1828  if (i != dotLoc) {
1829  if (j>=0) {
1830  value[i] = stdUncertainty[j--];
1831  }
1832  else {
1833  value[i] = '0';
1834  }
1835  }
1836  }
1837 
1838  return atof(value.c_str());
1839 }
1840 
1841 string Record::getStringAfter(string searchString, size_t len) {
1842  if (lines.size() > 1) {
1843 
1844  for (int i=1; i<lines.size(); ++i) {
1845 
1846  int tagPos = lines[i].find(searchString);
1847 
1848  if (tagPos != std::string::npos) {
1849 
1850  tagPos += searchString.length();
1851 
1852  string record = lines[i].substr(tagPos, len);
1853 
1854  return record;
1855  }
1856  }
1857  }
1858  return "";
1859 }
1860 
1861 // Parse a halflife from a record
1862 // Converts the units to seconds
1863 // Returns the halflife, or a negative number upon failure
1864 double Record::parseHalfLife(int startPos, int endPos) {
1865  if (lines.empty()) {
1866  egsWarning("Record::parseHalfLife: Error: Record is empty\n");
1867  return -5;
1868  }
1869  if (lines.front().length() < startPos) {
1870  egsWarning("Record::parseHalfLife: Warning: Record too short to "
1871  "contain desired quantity\n");
1872  return -5;
1873  }
1874 
1875  string halfLifeStr = egsTrimString(lines.front().substr(startPos-1,
1876  endPos-startPos+1));
1877 
1878  // Return -1 for stable
1879  if (halfLifeStr.substr(0,5).compare("STABLE") == 0) {
1880  return -1;
1881  }
1882 
1883  // Store the length of the numeric part of the string in i
1884  unsigned int numLength;
1885  for (numLength = 0; numLength < halfLifeStr.length(); numLength++) {
1886  if (!isdigit(halfLifeStr[numLength])
1887  && halfLifeStr.at(numLength) != '.') {
1888 
1889  break;
1890  }
1891  }
1892 
1893  // If there was no numeric component return -2
1894  if (halfLifeStr.size() < numLength+2) {
1895  return -2;
1896  }
1897 
1898  // Get the numeric part
1899  double hl = atof(halfLifeStr.substr(0, numLength).c_str());
1900 
1901  // Convert to units of seconds
1902  if (halfLifeStr.size()>numLength+2) {
1903  string units = halfLifeStr.substr(numLength+1, 2);
1904  if (units.compare("Y ") == 0) {
1905  hl *= 31556925.26;
1906  }
1907  else if (units.compare("D ") == 0) {
1908  hl *= 86400;
1909  }
1910  else if (units.compare("H ") == 0) {
1911  hl *= 3600;
1912  }
1913  else if (units.compare("M ") == 0) {
1914  hl *= 60;
1915  }
1916  else if (units.compare("S ") == 0) {
1917  hl *= 1;
1918  }
1919  else if (units.compare("MS") == 0) {
1920  hl *= 1E-3;
1921  }
1922  else if (units.compare("US") == 0) {
1923  hl *= 1E-6;
1924  }
1925  else if (units.compare("NS") == 0) {
1926  hl *= 1E-9;
1927  }
1928  else if (units.compare("PS") == 0) {
1929  hl *= 1E-12;
1930  }
1931  else if (units.compare("FS") == 0) {
1932  hl *= 1E-15;
1933  }
1934  else if (units.compare("AS") == 0) {
1935  hl *= 1E-18;
1936  }
1937  else {
1938  return -3;
1939  }
1940  }
1941  else if (halfLifeStr.size()>numLength+1) {
1942  string units = halfLifeStr.substr(numLength+1, 1);
1943  if (units.compare("Y") == 0) {
1944  hl *= 31556925.26;
1945  }
1946  else if (units.compare("D") == 0) {
1947  hl *= 86400;
1948  }
1949  else if (units.compare("H") == 0) {
1950  hl *= 3600;
1951  }
1952  else if (units.compare("M") == 0) {
1953  hl *= 60;
1954  }
1955  else if (units.compare("S") == 0) {
1956  hl *= 1;
1957  }
1958  else {
1959  return -3;
1960  }
1961  }
1962  else {
1963  hl = -4;
1964  }
1965 
1966  return hl;
1967 }
1968 
1969 // Comment Record
1970 CommentRecord::CommentRecord(vector<string> ensdf):Record(ensdf) {
1971  processEnsdf();
1972 }
1973 
1974 void CommentRecord::processEnsdf() {
1975  if (!lines.empty()) {
1976 
1977  comment = lines.front();
1978  }
1979 }
1980 
1981 string CommentRecord::getComment() {
1982  return comment;
1983 }
1984 
1985 // Parent Record
1986 ParentRecord::ParentRecord(vector<string> ensdf):Record(ensdf) {
1987  processEnsdf();
1988 }
1989 
1990 void ParentRecord::processEnsdf() {
1991  halfLife = parseHalfLife(40, 49);
1992 
1993  // Ground state Q-value in keV
1994  // (total energy available for g.s. -> g.s. transition
1995  // It will always be a positive number
1996  // We convert to MeV
1997  Q = recordToDouble(65, 74) / 1000.;
1998 
1999  // If the Q was not contained in the record it returned -1
2000  if (Q == -0.001) {
2001  egsWarning("ParentRecord::processEnsdf: Warning: No Q-value given, any "
2002  "positron records will give errors\n");
2003  Q = 0.;
2004  }
2005 }
2006 
2007 double ParentRecord::getHalfLife() const {
2008  return halfLife;
2009 }
2010 
2011 double ParentRecord::getQ() const {
2012  return Q;
2013 }
2014 
2015 ParentRecord *ParentRecordLeaf::getParentRecord() const {
2016  return getBranch();
2017 }
2018 
2019 ParentRecordLeaf::ParentRecordLeaf(ParentRecord
2020  *myRecord):Leaf<ParentRecord>(myRecord) {
2021 
2022 }
2023 
2024 // Normalization Record
2025 NormalizationRecord::NormalizationRecord(vector<string> ensdf,
2026  ParentRecord *myParent):Record(ensdf), ParentRecordLeaf(myParent) {
2027  processEnsdf();
2028 }
2029 
2030 void NormalizationRecord::processEnsdf() {
2031  normalizeRelative = recordToDouble(10, 19);
2032  normalizeTransition = recordToDouble(22, 29);
2033  normalizeBranch = recordToDouble(32, 39);
2034  normalizeBeta = recordToDouble(42, 49);
2035 
2036  // If the normalization is not specified, it will get initialized to zero
2037  // Change this to 1
2038  if (normalizeRelative < epsilon) {
2039  normalizeRelative = 1;
2040  }
2041  if (normalizeTransition < epsilon) {
2042  normalizeTransition = 1;
2043  }
2044  if (normalizeBranch < epsilon) {
2045  normalizeBranch = 1;
2046  }
2047  if (normalizeBeta < epsilon) {
2048  normalizeBeta = 1;
2049  }
2050 
2051  // Get the daughter element
2052  string element = recordToString(4, 5);
2053 
2054  // Get the Z
2055  Z = setZ(element);
2056 
2057  // Load atomic relaxations
2058  relaxations = new EGS_AtomicRelaxations();
2059  relaxations->loadData(Z);
2060 
2061  // Get the number of shells
2062  nshell = relaxations->getNShell(Z);
2063 
2064  egsInformation("NormalizationRecord::processEnsdf(): Z, nshell: %d %d\n",Z,nshell);
2065 }
2066 
2067 EGS_AtomicRelaxations *NormalizationRecord::getRelaxations() const {
2068  return relaxations;
2069 }
2070 
2071 int NormalizationRecord::getNShell() const {
2072  return nshell;
2073 }
2074 
2075 double NormalizationRecord::getBindingEnergy(int shell) const {
2076  return relaxations->getBindingEnergy(Z,shell);
2077 }
2078 
2079 void NormalizationRecord::relax(int shell,
2080  EGS_Float ecut, EGS_Float pcut,
2081  EGS_RandomGenerator *rndm, double &edep,
2082  EGS_SimpleContainer<EGS_RelaxationParticle> &particles) {
2083  relaxations->relax(Z,shell,ecut,pcut,rndm,edep,particles);
2084 }
2085 
2086 // Multiplier for converting relative photon intensity to photons per 100
2087 // decays in the parent through the decay branch or to photons per 100 neutron
2088 // captures in an (n,gamma) reaction. Required if the absolute photon intensity
2089 // can be calculated
2090 double NormalizationRecord::getRelativeMultiplier() const {
2091  return normalizeRelative;
2092 }
2093 
2094 // Multiplier for convert relative transition intensity (including conversion
2095 // electrons) to transitions per 100 decays of the parent through this decay
2096 // branch or per 100 neutron captures in an (n,gamma) reaction
2097 double NormalizationRecord::getTransitionMultiplier() const {
2098  return normalizeTransition;
2099 }
2100 
2101 // Branching ratio multiplier for converting intensity per 100 decays
2102 // through this decay branch to intensity per 100 decays of the parent nuclide
2103 double NormalizationRecord::getBranchMultiplier() const {
2104  return normalizeBranch;
2105 }
2106 
2107 // Multiplier for converting relative beta- and electron capture intensities to
2108 // intensities per 100 decays through this decay branch. Required if known
2109 double NormalizationRecord::getBetaMultiplier() const {
2110  return normalizeBeta;
2111 }
2112 
2113 NormalizationRecord *NormalizationRecordLeaf::getNormalizationRecord()
2114 const {
2115  return getBranch();
2116 }
2117 
2118 NormalizationRecordLeaf::NormalizationRecordLeaf(NormalizationRecord
2119  *myRecord):Leaf<NormalizationRecord>(myRecord) {
2120 
2121 }
2122 
2123 // Level Record
2124 LevelRecord::LevelRecord() {
2125  energy = 0;
2126  halfLife = 0;
2127  disintegrationIntensity = 0;
2128 }
2129 LevelRecord::LevelRecord(vector<string> ensdf):
2130  Record(ensdf) {
2131  processEnsdf();
2132  disintegrationIntensity = 0;
2133 }
2134 
2135 void LevelRecord::processEnsdf() {
2136  energy = recordToDouble(10, 19) / 1000.; // Convert keV to MeV
2137  halfLife = parseHalfLife(40, 49);
2138 }
2139 
2140 void LevelRecord::setLevelCanDecay(bool canDecayTmp) {
2141  canDecay = canDecayTmp;
2142 }
2143 
2144 bool LevelRecord::levelCanDecay() const {
2145  return canDecay;
2146 }
2147 
2148 void LevelRecord::resetDisintegrationIntensity() {
2149  disintegrationIntensity = 0;
2150 }
2151 
2152 void LevelRecord::cumulDisintegrationIntensity(double disintIntensity) {
2153  disintegrationIntensity += disintIntensity;
2154 }
2155 
2156 double LevelRecord::getDisintegrationIntensity() const {
2157  return disintegrationIntensity;
2158 }
2159 
2160 double LevelRecord::getEnergy() const {
2161  return energy;
2162 }
2163 
2164 double LevelRecord::getHalfLife() const {
2165  return halfLife;
2166 }
2167 
2168 LevelRecord *LevelRecordLeaf::getLevelRecord() const {
2169  return getBranch();
2170 }
2171 
2172 LevelRecordLeaf::LevelRecordLeaf(LevelRecord
2173  *myRecord):Leaf<LevelRecord>(myRecord) {
2174 
2175 }
2176 
2177 // Beta Record
2178 BetaRecordLeaf::BetaRecordLeaf(vector<string> ensdf,
2179  ParentRecord *myParent,
2180  NormalizationRecord *myNormalization,
2181  LevelRecord *myLevel):
2182  Record(ensdf),
2183  ParentRecordLeaf(myParent),
2184  NormalizationRecordLeaf(myNormalization),
2185  LevelRecordLeaf(myLevel) {
2186 
2187  numSampled = 0;
2188 
2189  // Set the Z and atomic weight for the daughter of this decay
2190  string id = egsRemoveWhite(lines.front().substr(0,5));
2191  Z = setZ(id);
2192 
2193  string atomicWeight;
2194  for (unsigned int i=0; i < id.length(); ++i) {
2195  if (!isdigit(id[i])) {
2196  break;
2197  }
2198  else {
2199  atomicWeight.push_back(id[i]);
2200  }
2201  }
2202  A = atoi(atomicWeight.c_str());
2203 
2204  // Get the forbiddenness
2205  if (lines.front().length() > 77) {
2206  string lambda;
2207  lambda.push_back(lines.front().at(77));
2208  forbidden = atoi(lambda.c_str());
2209  }
2210  else {
2211  forbidden = 0;
2212  }
2213 }
2214 int BetaRecordLeaf::getCharge() const {
2215  return q;
2216 }
2217 
2218 void BetaRecordLeaf::incrNumSampled() {
2219  numSampled++;
2220 }
2221 
2222 EGS_I64 BetaRecordLeaf::getNumSampled() const {
2223  return numSampled;
2224 }
2225 
2226 unsigned short int BetaRecordLeaf::getZ() const {
2227  return Z;
2228 }
2229 
2230 unsigned short int BetaRecordLeaf::getAtomicWeight() const {
2231  return A;
2232 }
2233 
2234 unsigned short int BetaRecordLeaf::getForbidden() const {
2235  return forbidden;
2236 }
2237 
2238 void BetaRecordLeaf::setSpectrum(EGS_AliasTable *bspec) {
2239  spectrum = bspec;
2240 }
2241 
2242 EGS_AliasTable *BetaRecordLeaf::getSpectrum() const {
2243  return spectrum;
2244 }
2245 
2246 // Beta- Record
2247 BetaMinusRecord::BetaMinusRecord(vector<string> ensdf,
2248  ParentRecord *myParent,
2249  NormalizationRecord *myNormalization,
2250  LevelRecord *myLevel):
2251  BetaRecordLeaf(ensdf, myParent,
2252  myNormalization, myLevel) {
2253  processEnsdf();
2254  q = -1;
2255  myLevel->cumulDisintegrationIntensity(betaIntensity);
2256 }
2257 
2258 void BetaMinusRecord::processEnsdf() {
2259  finalEnergy = recordToDouble(10, 19) / 1000.; // Convert keV to MeV
2260  betaIntensity = recordToDouble(22, 29);
2261  string betaIntensityStr = recordToString(22, 29);
2262  string betaIntensityUncStr = recordToString(30, 31);
2263 
2264  betaIntensityUnc = parseStdUncertainty(betaIntensityStr, betaIntensityUncStr);
2265  // If the uncertainty is 0 (i.e. not specified), set it to 100%
2266  if (betaIntensityUnc == 0) {
2267  betaIntensityUnc = betaIntensity;
2268  }
2269 
2270  if (getNormalizationRecord()) {
2271  double factor = getNormalizationRecord()->getBetaMultiplier() * getNormalizationRecord()->getBranchMultiplier();
2272 
2273  betaIntensity *= factor;
2274  betaIntensityUnc *= factor;
2275  }
2276 }
2277 
2278 double BetaMinusRecord::getFinalEnergy() const {
2279  return finalEnergy;
2280 }
2281 
2282 double BetaMinusRecord::getBetaIntensity() const {
2283  return betaIntensity;
2284 }
2285 
2286 double BetaMinusRecord::getBetaIntensityUnc() const {
2287  return betaIntensityUnc;
2288 }
2289 
2290 void BetaMinusRecord::setBetaIntensity(double newIntensity) {
2291  betaIntensity = newIntensity;
2292 }
2293 
2294 // Beta+ Record (and Electron Capture)
2295 BetaPlusRecord::BetaPlusRecord(vector<string> ensdf,
2296  ParentRecord *myParent,
2297  NormalizationRecord *myNormalization,
2298  LevelRecord *myLevel):
2299  BetaRecordLeaf(ensdf, myParent,
2300  myNormalization, myLevel) {
2301  processEnsdf();
2302  q = 1;
2303  myLevel->cumulDisintegrationIntensity(betaIntensity);
2304 }
2305 
2306 void BetaPlusRecord::processEnsdf() {
2307  finalEnergy = recordToDouble(10, 19) / 1000.; // Convert keV to MeV
2308  positronIntensity = recordToDouble(22, 29);
2309  string positronIntensityStr = recordToString(22, 29);
2310  string positronIntensityUncStr = recordToString(30, 31);
2311  ecIntensity = recordToDouble(32, 39);
2312  string ecIntensityStr = recordToString(32, 39);
2313  string ecIntensityUncStr = recordToString(40, 41);
2314 
2315  positronIntensityUnc = parseStdUncertainty(positronIntensityStr, positronIntensityUncStr);
2316  // If the uncertainty is 0 (i.e. not specified), set it to 100%
2317  if (positronIntensityUnc == 0) {
2318  positronIntensityUnc = positronIntensity;
2319  }
2320 
2321  ecIntensityUnc = parseStdUncertainty(ecIntensityStr, ecIntensityUncStr);
2322  // If the uncertainty is 0 (i.e. not specified), set it to 100%
2323  if (ecIntensityUnc == 0) {
2324  ecIntensityUnc = ecIntensity;
2325  }
2326 
2327  if (getNormalizationRecord()) {
2328  double factor = getNormalizationRecord()->getBetaMultiplier() * getNormalizationRecord()->getBranchMultiplier();
2329 
2330  positronIntensity *= factor;
2331  ecIntensity *= factor;
2332  positronIntensityUnc *= factor;
2333  ecIntensityUnc *= factor;
2334  }
2335 
2336  // The total intensity for this decay branch
2337  // A decay down this branch will then be split between positron or EC
2338  betaIntensity = positronIntensity + ecIntensity;
2339 
2340  // Re-normalize the positron intensity so that
2341  // positronIntensity + ecIntensity = 1, for easy sampling
2342  // We only set positronIntensity because that's all we need
2343  positronIntensity = positronIntensity / betaIntensity;
2344 
2345  // For positrons we may need to calculate the emission energy
2346  // E = Q - level_energy - 2*mc^2
2347  if (finalEnergy == 0 && positronIntensity > epsilon) {
2348  finalEnergy = getParentRecord()->getQ()
2349  - getLevelRecord()->getEnergy() - 1.022;
2350 
2351  if (finalEnergy < 0.) {
2352  egsWarning("BetaPlusRecord::processEnsdf: Error: Final energy of "
2353  "positron could not be calculated. Setting energy to zero!\n"
2354  );
2355  finalEnergy = 0.;
2356  }
2357  }
2358 
2359  if (ecIntensity > 0) {
2360  // Get the number of shells
2361  int nshell = getNormalizationRecord()->getNShell();
2362 
2363  double icK = getTag("CK=");
2364  double icL = getTag("CL=");
2365  double icM = getTag("CM=");
2366  double icN = getTag("CN=");
2367  double icO = getTag("CO=");
2368  double icP = getTag("CP=");
2369  double icQ = getTag("CQ=");
2370 
2371  // The K shell
2372  ecShellIntensity.push_back(icK);
2373 
2374  // The L1, L2, L3 shells
2375  // TODO: Equal probability is assigned to subshells!
2376  // This assumption is an approximation
2377  int numShellsToInclude = min(4,nshell);
2378  for (unsigned int i=1; i<numShellsToInclude; ++i) {
2379  ecShellIntensity.push_back(ecShellIntensity.back() + icL/(numShellsToInclude-1));
2380  }
2381  // Count number of shells as we go
2382  // Once we hit the number of shells for this element, return
2383  if (numShellsToInclude < 4) {
2384 // for (int i=0; i<ecShellIntensity.size(); ++i) {
2385 // egsInformation("BetaPlusRecord::processEnsdf: Shell %d: P=%f\n",i,ecShellIntensity[i]);
2386 // }
2387  return;
2388  }
2389 
2390  // The M1-M5 shells
2391  numShellsToInclude = min(9,nshell);
2392  for (unsigned int i=4; i<numShellsToInclude; ++i) {
2393  ecShellIntensity.push_back(ecShellIntensity.back() + icM/(numShellsToInclude-4));
2394  }
2395  if (numShellsToInclude < 9) {
2396 // for (int i=0; i<ecShellIntensity.size(); ++i) {
2397 // egsInformation("BetaPlusRecord::processEnsdf: Shell %d: P=%f\n",i,ecShellIntensity[i]);
2398 // }
2399  return;
2400  }
2401 
2402  // The N1-7 shells
2403  numShellsToInclude = min(16,nshell);
2404  for (unsigned int i=9; i<numShellsToInclude; ++i) {
2405  ecShellIntensity.push_back(ecShellIntensity.back() + icN/(numShellsToInclude-9));
2406  }
2407  if (numShellsToInclude < 16) {
2408 // for (int i=0; i<ecShellIntensity.size(); ++i) {
2409 // egsInformation("BetaPlusRecord::processEnsdf: Shell %d: P=%f\n",i,ecShellIntensity[i]);
2410 // }
2411  return;
2412  }
2413 
2414  // The O1-7 shells
2415  numShellsToInclude = min(23,nshell);
2416  for (unsigned int i=16; i<numShellsToInclude; ++i) {
2417  ecShellIntensity.push_back(ecShellIntensity.back() + icO/(numShellsToInclude-16));
2418  }
2419 
2420  if (numShellsToInclude < 23) {
2421 // for (int i=0; i<ecShellIntensity.size(); ++i) {
2422 // egsInformation("BetaPlusRecord::processEnsdf: Shell %d: P=%f\n",i,ecShellIntensity[i]);
2423 // }
2424  return;
2425  }
2426 
2427  // The P1-3 shells
2428  numShellsToInclude = min(26,nshell);
2429  for (unsigned int i=23; i<numShellsToInclude; ++i) {
2430  ecShellIntensity.push_back(ecShellIntensity.back() + icP/(numShellsToInclude-23));
2431  }
2432 
2433  if (numShellsToInclude < 26) {
2434 // for (int i=0; i<ecShellIntensity.size(); ++i) {
2435 // egsInformation("BetaPlusRecord::processEnsdf: Shell %d: P=%f\n",i,ecShellIntensity[i]);
2436 // }
2437  return;
2438  }
2439 
2440  // The Q1 shell
2441  numShellsToInclude = 27;
2442  ecShellIntensity.push_back(ecShellIntensity.back() + icQ/(numShellsToInclude-26));
2443 
2444 // for (int i=0; i<ecShellIntensity.size(); ++i) {
2445 // egsInformation("BetaPlusRecord::processEnsdf: Shell %d: P=%f\n",i,ecShellIntensity[i]);
2446 // }
2447  return;
2448  }
2449 }
2450 
2451 void BetaPlusRecord::relax(int shell,
2452  EGS_Float ecut, EGS_Float pcut,
2453  EGS_RandomGenerator *rndm, double &edep,
2454  EGS_SimpleContainer<EGS_RelaxationParticle> &particles) {
2455  getNormalizationRecord()->relax(shell,ecut,pcut,rndm,edep,particles);
2456 }
2457 
2458 double BetaPlusRecord::getFinalEnergy() const {
2459  return finalEnergy;
2460 }
2461 
2462 double BetaPlusRecord::getBetaIntensity() const {
2463  return betaIntensity;
2464 }
2465 
2466 double BetaPlusRecord::getPositronIntensity() const {
2467  return positronIntensity;
2468 }
2469 
2470 double BetaPlusRecord::getPositronIntensityUnc() const {
2471  return positronIntensityUnc;
2472 }
2473 
2474 double BetaPlusRecord::getECIntensityUnc() const {
2475  return ecIntensityUnc;
2476 }
2477 
2478 void BetaPlusRecord::setBetaIntensity(double newIntensity) {
2479  betaIntensity = newIntensity;
2480 }
2481 
2482 void BetaPlusRecord::setPositronIntensity(double newIntensity) {
2483  positronIntensity = newIntensity;
2484 }
2485 
2486 // Gamma Record
2487 GammaRecord::GammaRecord(vector<string> ensdf,
2488  ParentRecord *myParent,
2489  NormalizationRecord *myNormalization,
2490  LevelRecord *myLevel):
2491  Record(ensdf),
2492  ParentRecordLeaf(myParent),
2493  NormalizationRecordLeaf(myNormalization),
2494  LevelRecordLeaf(myLevel) {
2495  processEnsdf();
2496  q = 0;
2497  numGammaSampled = 0;
2498  numICSampled = 0;
2499  numIPSampled = 0;
2500  multipleTransitionProb = 0;
2501 }
2502 
2503 GammaRecord::GammaRecord(GammaRecord *gamma):
2504  Record(),
2505  ParentRecordLeaf(gamma->getParentRecord()),
2506  NormalizationRecordLeaf(gamma->getNormalizationRecord()),
2507  LevelRecordLeaf(gamma->getLevelRecord()) {
2508 
2509  numGammaSampled = gamma->numGammaSampled;
2510  numICSampled = gamma->numICSampled;
2511  numIPSampled = gamma->numIPSampled;
2512  decayEnergy = gamma->decayEnergy;
2513  transitionIntensity = gamma->transitionIntensity;
2514  multipleTransitionProb = gamma->multipleTransitionProb;
2515  gammaIntensity = gamma->gammaIntensity;
2516  gammaIntensityUnc = gamma->gammaIntensityUnc;
2517  icCoeff = gamma->icCoeff;
2518  icCoeffUnc = gamma->icCoeffUnc;
2519  ipCoeff = gamma->ipCoeff;
2520  ipCoeffUnc = gamma->ipCoeffUnc;
2521  q = gamma->q;
2522  finalLevel = gamma->finalLevel;
2523 }
2524 
2525 void GammaRecord::processEnsdf() {
2526  decayEnergy = recordToDouble(10, 19) / 1000.; // Convert keV to MeV
2527  gammaIntensity = recordToDouble(22, 29);
2528  string gammaIntensityStr = recordToString(22, 29);
2529  string gammaIntensityUncStr = recordToString(30, 31);
2530 
2531  // Get the internal conversion coefficient
2532  // Note: This might not equal the sum of the individual shell intensities,
2533  // because some of the shell intensities might not be known
2534  icCoeff = recordToDouble(56, 62);
2535  string icCoeffStr = recordToString(56, 62);
2536  string icCoeffUncStr = recordToString(63, 64);
2537 
2538  // If we don't find the gamma intensity, check for the first RI=
2539  if (gammaIntensity < epsilon) {
2540  gammaIntensity = getTag("RI=");
2541  }
2542 
2543  icCoeffUnc = parseStdUncertainty(icCoeffStr, icCoeffUncStr);
2544  // If the uncertainty is 0 (i.e. not specified), set it to 100%
2545  if (icCoeffUnc == 0) {
2546  icCoeffUnc = icCoeff;
2547  }
2548 
2549  // Check for internal pair production
2550  // The value and uncertainty are stored in 11 characters after IPC=
2551  string ipCoeffStr_tmp = getStringAfter("IPC=", 11);
2552  if (ipCoeffStr_tmp.length() > 0) {
2553  string ipCoeffStr = ipCoeffStr_tmp.substr(0, 9);
2554  string ipCoeffUncStr = ipCoeffStr_tmp.substr(9, 2);
2555  ipCoeff = atof(ipCoeffStr.c_str());
2556  ipCoeffUnc = atof(ipCoeffUncStr.c_str());
2557 
2558  if (ipCoeffUnc == 0) {
2559  ipCoeffUnc = ipCoeff;
2560  }
2561  }
2562  else {
2563  ipCoeff = 0;
2564  ipCoeffUnc = 0;
2565  }
2566 
2567  // Get the transition intensity instead if gamma still zero
2568  if (gammaIntensity < epsilon) {
2569  double ti = getTag("TI ");
2570  // Calculate the gamma intensity from it
2571  gammaIntensity = ti / ((1+icCoeff) * (1+ipCoeff));
2572  }
2573 
2574  // Set the uncertainty on the gamma intensity
2575  gammaIntensityUnc = parseStdUncertainty(gammaIntensityStr, gammaIntensityUncStr);
2576  // If the uncertainty is 0 (i.e. not specified), set it to 100%
2577  if (gammaIntensityUnc == 0) {
2578  gammaIntensityUnc = gammaIntensity;
2579  }
2580 
2581  if (getNormalizationRecord()) {
2582  double factor = getNormalizationRecord()->getRelativeMultiplier() *
2583  getNormalizationRecord()->getBranchMultiplier();
2584 
2585  gammaIntensity *= factor;
2586  gammaIntensityUnc *= factor;
2587  }
2588 
2589  // Calculate the total transition intensity
2590  transitionIntensity = gammaIntensity * (1+icCoeff) * (1+ipCoeff);
2591 
2592  if (icCoeff > 0) {
2593  // Get the number of shells
2594  int nshell = getNormalizationRecord()->getNShell();
2595 
2596  double icK = getTag("KC=");
2597  double icL = getTag("LC=");
2598  double icM = getTag("MC=");
2599  double icN = getTag("NC=");
2600  double icO = getTag("OC=");
2601  double icP = getTag("PC=", "I"); // Make sure "IPC" doesn't get matched
2602  double icQ = getTag("QC=");
2603 
2604  // The K shell
2605  icIntensity.push_back(icK / icCoeff);
2606 
2607  // The L1, L2, L3 shells
2608  // TODO: Equal probability is assigned to subshells!
2609  // This assumption is an approximation
2610  int numShellsToInclude = min(4,nshell);
2611  for (unsigned int i=1; i<numShellsToInclude; ++i) {
2612  icIntensity.push_back(icIntensity.back() + (icL / icCoeff)/(numShellsToInclude-1));
2613  }
2614  // Count number of shells as we go
2615  // Once we hit the number of shells for this element, return
2616  if (numShellsToInclude < 4) {
2617 // for (int i=0; i<icIntensity.size(); ++i) {
2618 // egsInformation("GammaRecord::processEnsdf: Shell %d: P=%f\n",i,icIntensity[i]);
2619 // }
2620  return;
2621  }
2622 
2623  // The M1-M5 shells
2624  numShellsToInclude = min(9,nshell);
2625  for (unsigned int i=4; i<numShellsToInclude; ++i) {
2626  icIntensity.push_back(icIntensity.back() + (icM / icCoeff)/(numShellsToInclude-4));
2627  }
2628  if (numShellsToInclude < 9) {
2629 // for (int i=0; i<icIntensity.size(); ++i) {
2630 // egsInformation("GammaRecord::processEnsdf: Shell %d: P=%f\n",i,icIntensity[i]);
2631 // }
2632  return;
2633  }
2634 
2635  // The N1-7 shells
2636  numShellsToInclude = min(16,nshell);
2637  for (unsigned int i=9; i<numShellsToInclude; ++i) {
2638  icIntensity.push_back(icIntensity.back() + (icN / icCoeff)/(numShellsToInclude-9));
2639  }
2640  if (numShellsToInclude < 16) {
2641 // for (int i=0; i<icIntensity.size(); ++i) {
2642 // egsInformation("GammaRecord::processEnsdf: Shell %d: P=%f\n",i,icIntensity[i]);
2643 // }
2644  return;
2645  }
2646 
2647  // The O1-7 shells
2648  numShellsToInclude = min(23,nshell);
2649  for (unsigned int i=16; i<numShellsToInclude; ++i) {
2650  icIntensity.push_back(icIntensity.back() + (icO / icCoeff)/(numShellsToInclude-16));
2651  }
2652 
2653  if (numShellsToInclude < 23) {
2654 // for (int i=0; i<icIntensity.size(); ++i) {
2655 // egsInformation("GammaRecord::processEnsdf: Shell %d: P=%f\n",i,icIntensity[i]);
2656 // }
2657  return;
2658  }
2659 
2660  // The P1-3 shells
2661  numShellsToInclude = min(26,nshell);
2662  for (unsigned int i=23; i<numShellsToInclude; ++i) {
2663  icIntensity.push_back(icIntensity.back() + (icP / icCoeff)/(numShellsToInclude-23));
2664  }
2665 
2666  if (numShellsToInclude < 26) {
2667 // for (int i=0; i<icIntensity.size(); ++i) {
2668 // egsInformation("GammaRecord::processEnsdf: Shell %d: P=%f\n",i,icIntensity[i]);
2669 // }
2670  return;
2671  }
2672 
2673  // The Q1 shell
2674  numShellsToInclude = 27;
2675  icIntensity.push_back(icIntensity.back() + (icQ / icCoeff)/(numShellsToInclude-26));
2676 
2677 // for (int i=0; i<icIntensity.size(); ++i) {
2678 // egsInformation("GammaRecord::processEnsdf: Shell %d: P=%f\n",i,icIntensity[i]);
2679 // }
2680  return;
2681 
2682  }
2683 }
2684 
2685 double GammaRecord::getBindingEnergy(int shell) const {
2686  return getNormalizationRecord()->getBindingEnergy(shell);
2687 }
2688 
2689 void GammaRecord::relax(int shell,
2690  EGS_Float ecut, EGS_Float pcut,
2691  EGS_RandomGenerator *rndm, double &edep,
2692  EGS_SimpleContainer<EGS_RelaxationParticle> &particles) {
2693  getNormalizationRecord()->relax(shell,ecut,pcut,rndm,edep,particles);
2694 }
2695 
2696 double GammaRecord::getDecayEnergy() const {
2697  return decayEnergy;
2698 }
2699 
2700 double GammaRecord::getMultiTransitionProb() const {
2701  return multipleTransitionProb;
2702 }
2703 
2704 void GammaRecord::setMultiTransitionProb(double newIntensity) {
2705  multipleTransitionProb = newIntensity;
2706 }
2707 
2708 double GammaRecord::getTransitionIntensity() const {
2709  return transitionIntensity;
2710 }
2711 
2712 double GammaRecord::getGammaIntensity() const {
2713  return gammaIntensity;
2714 }
2715 
2716 double GammaRecord::getGammaIntensityUnc() const {
2717  return gammaIntensityUnc;
2718 }
2719 
2720 double GammaRecord::getICIntensity() const {
2721  return icCoeff;
2722 }
2723 
2724 double GammaRecord::getICIntensityUnc() const {
2725  return icCoeffUnc;
2726 }
2727 
2728 double GammaRecord::getIPIntensity() const {
2729  return ipCoeff;
2730 }
2731 
2732 double GammaRecord::getIPIntensityUnc() const {
2733  return ipCoeffUnc;
2734 }
2735 
2736 void GammaRecord::setTransitionIntensity(double newIntensity) {
2737  transitionIntensity = newIntensity;
2738 }
2739 
2740 void GammaRecord::setGammaIntensity(double newIntensity) {
2741  gammaIntensity = newIntensity;
2742 }
2743 
2744 void GammaRecord::setICIntensity(double newIntensity) {
2745  icCoeff = newIntensity;
2746 }
2747 
2748 int GammaRecord::getCharge() const {
2749  return q;
2750 }
2751 
2752 void GammaRecord::incrGammaSampled() {
2753  numGammaSampled++;
2754 }
2755 
2756 void GammaRecord::incrICSampled() {
2757  numICSampled++;
2758 }
2759 
2760 void GammaRecord::incrIPSampled() {
2761  numIPSampled++;
2762 }
2763 
2764 EGS_I64 GammaRecord::getGammaSampled() const {
2765  return numGammaSampled;
2766 }
2767 
2768 EGS_I64 GammaRecord::getICSampled() const {
2769  return numICSampled;
2770 }
2771 
2772 EGS_I64 GammaRecord::getIPSampled() const {
2773  return numIPSampled;
2774 }
2775 
2776 LevelRecord *GammaRecord::getFinalLevel() const {
2777  return finalLevel;
2778 }
2779 
2780 void GammaRecord::setFinalLevel(LevelRecord *newLevel) {
2781  finalLevel = newLevel;
2782 }
2783 
2784 // Alpha Record
2785 AlphaRecord::AlphaRecord(vector<string> ensdf,
2786  ParentRecord *myParent,
2787  NormalizationRecord *myNormalization,
2788  LevelRecord *myLevel):
2789  Record(ensdf),
2790  ParentRecordLeaf(myParent), NormalizationRecordLeaf(myNormalization),
2791  LevelRecordLeaf(myLevel) {
2792 
2793  processEnsdf();
2794  q = 2;
2795  numSampled = 0;
2796  myLevel->cumulDisintegrationIntensity(alphaIntensity);
2797 }
2798 
2799 void AlphaRecord::processEnsdf() {
2800  finalEnergy = recordToDouble(10, 19) / 1000.; // Convert keV to MeV
2801  alphaIntensity = recordToDouble(22, 29);
2802 
2803  string alphaIntensityStr = recordToString(22, 29);
2804  string alphaIntensityUncStr = recordToString(30, 31);
2805 
2806  alphaIntensityUnc = parseStdUncertainty(alphaIntensityStr, alphaIntensityUncStr);
2807  // If the uncertainty is 0 (i.e. not specified), set it to 100%
2808  if (alphaIntensityUnc == 0) {
2809  alphaIntensityUnc = alphaIntensity;
2810  }
2811 
2812  if (getNormalizationRecord()) {
2813  alphaIntensity *= getNormalizationRecord()->getBranchMultiplier();
2814  alphaIntensityUnc *= getNormalizationRecord()->getBranchMultiplier();
2815  }
2816 }
2817 
2818 double AlphaRecord::getFinalEnergy() const {
2819  return finalEnergy;
2820 }
2821 
2822 double AlphaRecord::getAlphaIntensity() const {
2823  return alphaIntensity;
2824 }
2825 
2826 double AlphaRecord::getAlphaIntensityUnc() const {
2827  return alphaIntensityUnc;
2828 }
2829 
2830 void AlphaRecord::setAlphaIntensity(double newIntensity) {
2831  alphaIntensity = newIntensity;
2832 }
2833 
2834 int AlphaRecord::getCharge() const {
2835  return q;
2836 }
2837 
2838 void AlphaRecord::incrNumSampled() {
2839  numSampled++;
2840 }
2841 
2842 EGS_I64 AlphaRecord::getNumSampled() const {
2843  return numSampled;
2844 }
2845 
EGS_AliasTable
A class for sampling random values from a given probability distribution using the alias table techni...
Definition: egs_alias_table.h:63
EGS_SimpleContainer< EGS_RelaxationParticle >
egs_ensdf.h
The ensdf library header file.
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::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_RandomGenerator
Base random number generator class. All random number generators should be derived from this class...
Definition: egs_rndm.h:67
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
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:134
EGS_Ensdf::~EGS_Ensdf
~EGS_Ensdf()
Destructor.
Definition: egs_ensdf.cpp:238
epsilon
const EGS_Float epsilon
The epsilon constant for floating point comparisons.
Definition: egs_functions.h:61
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
egsWarning
EGS_InfoFunction EGS_EXPORT egsWarning
Always use this function for reporting warnings.
Definition: egs_functions.cpp:135
Generated by   doxygen 1.8.5