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