Addign emitdb resources to radixio.

radixio/CMakeLists.txt

@@ -7,6 +7,7 @@ SET(SOURCE
    cfgfile.cc
    eafstream.cc
    endian.cc
+   emitdb.cc
    decaydb.cc
    gfsfile.cc
    f71stream.cc
@@ -21,6 +22,7 @@ SET(HEADERS
    cfgfile.hh
    eafstream.hh
    endian.hh
+   emitdb.hh
    decaydb.hh
    gfsfile.hh
    f71stream.hh

radixio/create_emission_data_db.py

+import os
+
+# set globals
+XGAM_DATA = "origen.rev04.mpdkxgam.data"
+BETA_DATA = "origen.rev00.ensdf95beta.data"
+CLANG_FMT = "clang-format"
+
+def open_file(file_name=None):
+    with open(file_name, "r") as fn:
+        datalines = fn.read().splitlines()
+    if "mpdkxgam" not in os.path.split(file_name)[1]:
+        datalines = datalines[3 : ]
+    data_entries = split_entries(datalines)
+    return data_entries
+
+def isXgam(datalines=None):
+    line = datalines[0]
+    answer = bool(len(line.split()) > 5)
+    return answer
+
+def header_name(header_file_name=None):
+    if header_file_name is None:
+        header_file_name = "emitdb.hh"
+    name = os.path.splitext(os.path.split(header_file_name)[1])[0]
+    return name
+
+def struct_name(header_file_name="emitdb.hh"):
+    hname = header_name(header_file_name=header_file_name)
+    if "db" in hname:
+        name = hname.split("db")[0][0].upper() + hname.split("db")[0][1 : ].lower() +  "DBEntry"
+    else:
+        name = hname[0].upper() + hname[1 : ].lower() + "Entry"
+    return name
+
+def class_name(header_file_name="emitdb.hh"):
+    hname = header_name(header_file_name=header_file_name)
+    if "db" in hname:
+        name = hname.split("db")[0][0].upper() + hname.split("db")[0][1 : ].lower() + "Db"
+    else:
+        name = hname[0].upper() + hname[1 : ].lower()
+    return name
+
+def split_entries(datalines=None):
+    if isXgam(datalines=datalines):
+        zaid_lines = [i for i, _ in enumerate(datalines) if datalines[i][2].isdigit()]
+        entries = [datalines[zaid_lines[i] : zaid_lines[i + 1]] if i < len(zaid_lines) - 1 else datalines[zaid_lines[i] : ] for i, _ in enumerate(zaid_lines)]
+    else:
+        zaid_splits = [-1] + [i for i, _ in enumerate(datalines) if datalines[i].strip() == "0"]
+        entries = [datalines[zaid_splits[i] + 1 : zaid_splits[i + 1]] for i in range(len(zaid_splits) - 1)]
+    return entries
+
+def split_xgam_entries(datalines=None):
+    zaid_lines = [i for i, _ in enumerate(datalines) if datalines[i][2].isdigit()]
+    entries = [datalines[zaid_lines[i] : zaid_lines[i + 1]] if i < len(zaid_lines) - 1 else datalines[zaid_lines[i] : ] for i, _ in enumerate(zaid_lines)]
+    return entries
+
+def split_beta_entries(datalines=None):
+    zaid_splits = [-1] + [i for i, _ in enumerate(datalines) if datalines[i].strip() == "0"]
+    entries = [datalines[zaid_splits[i] + 1 : zaid_splits[i + 1]] for i in range(len(zaid_splits) - 1)]
+    return entries
+
+def getXgamZAID(entry=None):
+    zaid = int(entry[0][ : 7].strip())
+    return zaid
+
+def getBetaZAID(entry=None):
+    entry = [i.strip().split() for i in entry]
+    zaid = int(entry[0][0])
+    return zaid
+
+def getZ(zaid=None):
+    z = int(zaid / 10000)
+    return z
+
+def getA(zaid=None):
+    a = int((zaid - int(getZ(zaid=zaid) * 10000)) / 10)
+    return a
+
+def getM(zaid=None):
+    m = zaid - int(10000 * getZ(zaid=zaid)) - int(10 * getA(zaid=zaid))
+    return m
+
+def getXgamIsotopeInfo(entry=None):
+    zaid = getXgamZAID(entry=entry)
+    z = getZ(zaid=zaid)
+    a = getA(zaid=zaid)
+    m = getM(zaid=zaid)
+    return zaid, z, a, m
+
+def getBetaIsotopeInfo(entry=None):
+    zaid = getBetaZAID(entry=entry)
+    z = getZ(zaid=zaid)
+    a = getA(zaid=zaid)
+    m = getM(zaid=zaid)
+    return zaid, z, a, m
+
+def getXgamTotalLines(entry=None):
+    line_one = entry[0]
+    total_lines = int(line_one[10 : 19].strip())
+    return total_lines
+
+def getBetaTotalLines(entry=None):
+    entry = [i.strip().split() for i in entry]
+    beta_num = int(len(entry))
+    return beta_num
+
+def getXgamTotalEnergy(entry=None):
+    line_one = entry[0]
+    total_energy = float(line_one[60 : 70].strip())
+    return total_energy
+
+def getBetaTotalEnergy(intensities=None, energies=None, betaNum=None):
+    total_energy = sum([intensities[i] * energies[i] for i in range(betaNum)])
+    return total_energy
+
+def getXgamTotals(entry=None):
+    total_lines, total_energy = getXgamTotalLines(entry=entry), getXgamTotalEnergy(entry=entry)
+    return total_lines, total_energy
+
+def xgamZAIDS(xgam_entries=None):
+    xgam_zaids = [int(i[0][ : 7].strip()) for i in xgam_entries]
+    return xgam_zaids
+
+def betaZAIDS(beta_entries=None):
+    beta_zaids = [int(i[0].strip().split()[0]) for i in beta_entries]
+    return beta_zaids
+
+def shared_zaids(xgam_entries=None, beta_entries=None):
+    xgam_zaids, beta_zaids = xgamZAIDS(xgam_entries=xgam_entries), betaZAIDS(beta_entries=beta_entries)
+    all_zaids = list(set(xgam_zaids + beta_zaids))
+    return all_zaids
+
+def processPhotons(entry=None, revision_number=None):
+    total_lines = getXgamTotalLines(entry=entry)
+    index, position, energy, intensity = 1, 10, 0.0, 0.0
+    energies, intensities = [energy for _ in range(total_lines)], [intensity for _ in range(total_lines)]
+    lines_used = 1
+    next_line = entry[lines_used]
+    if revision_number == 3:
+        index_iterator = 2
+        index_threshold = 7
+    else:
+        index_iterator = 2.2
+        index_threshold = 7.6
+    for i in range(total_lines):
+        if int(index) == int(index_threshold):
+            lines_used += 1
+            next_line = entry[lines_used]
+            index = 1
+        position = int(10 * index)
+        energies[i] = float(next_line[position : position + 10].strip())
+        intensities[i] = float(next_line[position + 10 : position + 20].strip())
+        index += index_iterator
+    intensities_sorted = [i for _, i in sorted(zip(energies, intensities))]
+    energies_sorted = sorted(energies)
+    return energies_sorted, intensities_sorted
+
+def processBetas(entry=None):
+    beta_num = getBetaTotalLines(entry=entry)
+    entry = [i.strip().split() for i in entry]
+    intensities, energies = [[0.0 for _ in range(beta_num)] for _ in range(2)]
+    for i in range(beta_num):
+        intensities[i], energies[i] = [float(entry[i][j]) for j in range(1, 3)]
+    return energies, intensities
+
+def make_string(xgamEntry=None, betaEntry=None, xgam_revision_number=None):
+    if all([i is not None for i in [xgamEntry, betaEntry]]):
+        z, a, m = getXgamIsotopeInfo(entry=xgamEntry)[1 : ]
+        xgam_num = getXgamTotalLines(entry=xgamEntry)
+        beta_num = getBetaTotalLines(entry=betaEntry)
+        xgam_energy = getXgamTotalEnergy(entry=xgamEntry)
+        beta_energies, beta_intensities = processBetas(entry=betaEntry)
+        beta_energy = getBetaTotalEnergy(intensities=beta_intensities, energies=beta_energies, betaNum=beta_num)
+        xgam_energies, xgam_intensities = processPhotons(entry=xgamEntry, revision_number=xgam_revision_number)
+        string = "{{\n      {0},{1},{2},\n      {3},{4},{5}f,{6}f,\n      ".format(str(z), str(a), str(m), str(xgam_num), str(beta_num), str(xgam_energy), str(beta_energy))
+        string += "{"
+        for i in xgam_energies:
+            string += "{0}f,".format(i)
+        string += "},{"
+        for i in xgam_intensities:
+            string += "{0}f,".format(i)
+        string += "},{"
+        for i in beta_energies:
+            string += "{0}f,".format(i)
+        string += "},{"
+        for i in beta_intensities:
+            string += "{0}f,".format(i)
+        string += "}"
+        string += "};"
+    else:
+        if xgamEntry is None:
+            z, a, m = getBetaIsotopeInfo(entry=betaEntry)[1 : ]
+            beta_num = getBetaTotalLines(entry=betaEntry)
+            energies, intensities = processBetas(entry=betaEntry)
+            total_energy = getBetaTotalEnergy(intensities=intensities, energies=energies, betaNum=beta_num)
+            string = "{{\n      {0},{1},{2},\n      0,{3},0.0f,{4}f,\n      ".format(str(z), str(a), str(m), str(beta_num), str(total_energy))
+            string += "{},{},"
+            string += "{"
+            for i in energies:
+                string += "{0}f,".format(i)
+            string += "},{"
+            for i in intensities:
+                string += "{0}f,".format(i)
+            string += "}"
+            string += "};"
+        else:
+            z, a, m = getXgamIsotopeInfo(entry=xgamEntry)[1 : ]
+            xgam_num = getXgamTotalLines(entry=xgamEntry)
+            total_energy = getXgamTotalEnergy(entry=xgamEntry)
+            energies, intensities = processPhotons(entry=xgamEntry)
+            string = "{{\n      {0},{1},{2},\n      {3},0,{4}f,0.0f,\n      ".format(str(z), str(a), str(m), str(xgam_num), str(total_energy))
+            string += "{"
+            for i in energies:
+                string += "{0}f,".format(i)
+            string += "},{"
+            for i in intensities:
+                string += "{0}f,".format(i)
+            string += "},"
+            string += "{},{}"
+            string += "};"
+    return string
+
+def getXgamEntryByZAID(zaid=None, entries=None):
+    xgam_zaids = xgamZAIDS(xgam_entries=entries)
+    if zaid not in xgam_zaids:
+        entry = None
+    else:
+        entry = entries[xgam_zaids.index(zaid)]
+    return entry
+
+def getBetaEntryByZAID(zaid=None, entries=None):
+    beta_zaids = betaZAIDS(beta_entries=entries)
+    if zaid not in beta_zaids:
+        entry = None
+    else:
+        entry = entries[beta_zaids.index(zaid)]
+    return entry
+
+def header_file(header_file_name="emitdb.hh"):
+    hname = header_name(header_file_name=header_file_name)
+    header_def = hname.upper()
+    sname = struct_name(header_file_name=header_file_name)
+    cname = class_name(header_file_name=header_file_name)
+    lines_to_write = ["#ifndef RADIX_RADIXIO_{0}_HH_".format(header_def), "#define RADIX_RADIXIO_{0}_HH_".format(header_def), "#include <memory>", "#include <string>", "#include <unordered_map>", "#include <utility>", "#include <vector>", "\nnamespace radix\n{\n", "struct {0}\n{{".format(sname), "/**", "* symbol", "**/", "int z;", "/**", "* atomic number", "**/", "int a;", "/**", "* metastable state", "**/", "int m;", "/**", "* The variable nLines, total number of photon lines in dataset", "**/", "int nLines;", "/**", "* The variable nBetas, total number of beta lines in dataset", "**/", "int nBetas;", "/**", "* The variable xgamEnergy, total energy of photons emitted by isotope", "**/", "float xgamEnergy;", "/**", "* The variable betaEnergy, total energy of betas emitted by isotope", "**/", "float betaEnergy;", "/**", "* The energies of each photon in dataset", "**/", "std::vector<float> photonEnergies;", "/**", "* The intensities of each photon in dataset", "**/", "std::vector<float> photonIntensities;", "/**", "* The maximum endpoint energies of each beta transition in dataset", "**/", "std::vector<float> betaEnergies;", "/**", "* The intensities of each beta in dataset", "**/", "std::vector<float> betaIntensities;", "int zaid() { return z * 10000 + a * 10 + m;}", "bool operator==(const {0}& rhs) const".format(sname), "{\nreturn (z == rhs.z && a == rhs.a && m == rhs.m);\n}", "\n}}; // struct {0}\n".format(sname), "class {0}\n{{".format(cname)]
+    class_variables = [" public:", "  std::unordered_map<int,{0}> mIsotope;".format(sname), "  std::string mIsotopeTitle;"]
+    public_def = ["public:", "  {0}();".format(cname)]
+    public_methods = ["  const std::string& isotope_title() const;", "  void set_isotope_title(const std::string& title);", "  const std::unordered_map<int, EmitDBEntry>& isotopes() const;", "  std::unordered_map<int, EmitDBEntry>& isotopes();", "const EmitDBEntry& emitter(const int zaid) const;", "  const std::vector<std::vector<float>> photons(const int zaid) const;", "  std::vector<std::vector<float>> photons(const int zaid);", "  const std::vector<std::vector<float>> betas(const int zaid) const;", "  std::vector<std::vector<float>> betas(const int zaid);"]
+    end_class_namespace = ["\n}}; // class {0}".format(cname), "} // namespace radix"]
+    hash_function = ["\nnamespace std\n{", "template<>\nstruct hash<radix::{0}>".format(sname), "{{\nstd::size_t operator()(const radix::{0}& rhs) const".format(sname), "{\nusing std::size_t;", "using std::hash;", "using std::string;", "return rhs.z + rhs.a + rhs.m + rhs.nLines + rhs.nBetas;\n}", "}}; // struct hash<radix::{0}>".format(sname), "} // namespace std\n"]
+    end_header_def = ["#endif /** RADIX_RADIXIO_{0}_HH_ */".format(header_def)]
+    if class_variables:
+        lines_to_write += class_variables
+    lines_to_write += public_def
+    if public_methods:
+        lines_to_write += public_methods
+    lines_to_write += end_class_namespace + hash_function + end_header_def
+    return lines_to_write
+
+def implementation_file(header_file_name="emitdb.hh", xgam_data_file_name=XGAM_DATA, beta_data_file_name=BETA_DATA):
+    cname = class_name(header_file_name=header_file_name)
+    revision_number = int(xgam_data_file_name.split(os.path.sep)[-1].split(".")[1][-1])
+    include_lines = ["#include \"radixio/{0}\"".format(os.path.split(header_file_name)[1])]
+    namespace_line = ["\nnamespace radix\n{\n"]
+    constructor_def = ["{0}::{0}()\n{{".format(cname)]
+    all_xgam_entries, all_beta_entries = open_file(file_name=xgam_data_file_name), open_file(file_name=beta_data_file_name)
+    all_zaids = sorted(shared_zaids(xgam_entries=all_xgam_entries, beta_entries=all_beta_entries))
+    isotopes = [0 for _, _ in enumerate(all_zaids)]
+    for i, _ in enumerate(all_zaids):
+        xgam_entry = getXgamEntryByZAID(zaid=all_zaids[i], entries=all_xgam_entries)
+        beta_entry = getBetaEntryByZAID(zaid=all_zaids[i], entries=all_beta_entries)
+        new_string = make_string(xgamEntry=xgam_entry, betaEntry=beta_entry, xgam_revision_number=revision_number)
+        isotopes[i] = "  mIsotope[{0}]={1}".format(str(all_zaids[i]), new_string)
+    constructor_close = ["}"]
+    public_methods = ["const std::string& {0}::isotope_title() const\n{{".format(cname), "  return mIsotopeTitle;\n}", "void {0}::set_isotope_title(const std::string& title)\n{{".format(cname), "  mIsotopeTitle = title;\n}", "  const std::unordered_map<int, EmitDBEntry>& {0}::isotopes() const\n{{".format(cname), "  return mIsotope;\n}", "  std::unordered_map<int, EmitDBEntry>& {0}::isotopes()\n{{".format(cname), "  return mIsotope;\n}", "const EmitDBEntry& {0}::emitter(const int zaid) const\n{{".format(cname), "  const auto& it = mIsotope.find(zaid);", "  return it->second;\n}", "const std::vector<std::vector<float>> {0}::photons(const int zaid) const\n{{".format(cname), "  const auto& it = mIsotope.find(zaid);", "  std::vector<float> energies = it->second.photonEnergies;", "  std::vector<float> intensities = it->second.photonIntensities;", "  std::vector<std::vector<float>> paired = {energies, intensities};", "  return paired;\n}", "std::vector<std::vector<float>> {0}::photons(const int zaid)\n{{".format(cname), "  const auto& it = mIsotope.find(zaid);", "  std::vector<float> energies = it->second.photonEnergies;", "  std::vector<float> intensities = it->second.photonIntensities;", "  std::vector<std::vector<float>> paired = {energies, intensities};", "  return paired;\n}", "const std::vector<std::vector<float>> {0}::betas(const int zaid) const\n{{".format(cname), "  const auto& it = mIsotope.find(zaid);", "  std::vector<float> energies = it->second.betaEnergies;", "  std::vector<float> intensities = it->second.betaIntensities;", "  std::vector<std::vector<float>> paired = {energies, intensities};", "  return paired;\n}", "std::vector<std::vector<float>> {0}::betas(const int zaid)\n{{".format(cname), "  const auto& it = mIsotope.find(zaid);", "  std::vector<float> energies = it->second.betaEnergies;", "  std::vector<float> intensities = it->second.betaIntensities;", "  std::vector<std::vector<float>> paired = {energies, intensities};", "  return paired;\n}"]
+    namespace_close = ["\n} // namespace radix"]
+    implementation_lines = include_lines + namespace_line + constructor_def + isotopes + constructor_close + public_methods + namespace_close
+    return implementation_lines
+
+def write_files(implementation_file_name="emitdb.cc", implementation_lines=None, header_file_name="emitdb.hh", header_lines=None):
+    with open(header_file_name, "w") as fn:
+        for i in header_lines:
+            print >> fn, i
+    with open(implementation_file_name, "w") as fn:
+        for i in implementation_lines:
+            print >> fn, i
+    return None
+
+def format_db_files(implementation_file_name="emitdb.cc", header_file_name="emitdb.hh", clang_format=CLANG_FMT):
+    os.system("{0} {1} -i -style=file".format(clang_format, header_file_name))
+    os.system("{0} {1} -i -style=file".format(clang_format, implementation_file_name))
+    return None
+
+def create_files(implementation_file_name="emitdb.cc", header_file_name="emitdb.hh", xgam_data_file_name=XGAM_DATA, beta_data_file_name=BETA_DATA, clang_format=CLANG_FMT):
+    header_lines = header_file(header_file_name=header_file_name)
+    implementation_lines = implementation_file(header_file_name=header_file_name, xgam_data_file_name=xgam_data_file_name, beta_data_file_name=beta_data_file_name)
+    write_files(implementation_file_name=implementation_file_name, implementation_lines=implementation_lines, header_file_name=header_file_name, header_lines=header_lines)
+    format_db_files(implementation_file_name=implementation_file_name, header_file_name=header_file_name, clang_format=clang_format)
+    return None
+
+if __name__ == '__main__':
+    create_files(implementation_file_name="/home/jakeinman/projects/radix/source/radixio/emitdb.cc", header_file_name="/home/jakeinman/projects/radix/source/radixio/emitdb.hh", xgam_data_file_name=XGAM_DATA, beta_data_file_name="/home/jakeinman/software/scale/data/origen_data/origen.rev00.ensdf95beta.data", clang_format="/usr/bin/clang-format")

radixio/emitdb.hh

+#ifndef RADIX_RADIXIO_EMITDB_HH_
+#define RADIX_RADIXIO_EMITDB_HH_
+#include <memory>
+#include <string>
+#include <unordered_map>
+#include <utility>
+#include <vector>
+
+#include "radixcore/visibility.hh"
+
+namespace radix
+{
+struct EmitDBEntry
+{
+  /**
+   * symbol
+   **/
+  int z;
+  /**
+   * atomic number
+   **/
+  int a;
+  /**
+   * metastable state
+   **/
+  int m;
+  /**
+   * The variable nLines, total number of photon lines in dataset
+   **/
+  int nLines;
+  /**
+   * The variable nBetas, total number of beta lines in dataset
+   **/
+  int nBetas;
+  /**
+   * The variable xgamEnergy, total energy of photons emitted by isotope
+   **/
+  float xgamEnergy;
+  /**
+   * The variable betaEnergy, total energy of betas emitted by isotope
+   **/
+  float betaEnergy;
+  /**
+   * The energies of each photon in dataset
+   **/
+  std::vector<float> photonEnergies;
+  /**
+   * The intensities of each photon in dataset
+   **/
+  std::vector<float> photonIntensities;
+  /**
+   * The maximum endpoint energies of each beta transition in dataset
+   **/
+  std::vector<float> betaEnergies;
+  /**
+   * The intensities of each beta in dataset
+   **/
+  std::vector<float> betaIntensities;
+  int zaid() { return z * 10000 + a * 10 + m; }
+  bool operator==(const EmitDBEntry& rhs) const
+  {
+    return (z == rhs.z && a == rhs.a && m == rhs.m);
+  }
+
+};  // struct EmitDBEntry
+
+class RADIX_PUBLIC EmitDb
+{
+ public:
+  std::unordered_map<int, EmitDBEntry> mIsotope;
+  std::string mIsotopeTitle;
+
+ public:
+  EmitDb();
+  const std::string& isotope_title() const;
+  void set_isotope_title(const std::string& title);
+  const std::unordered_map<int, EmitDBEntry>& isotopes() const;
+  std::unordered_map<int, EmitDBEntry>& isotopes();
+  const EmitDBEntry& emitter(const int zaid) const;
+  const std::vector<std::vector<float>> photons(const int zaid) const;
+  std::vector<std::vector<float>> photons(const int zaid);
+  const std::vector<std::vector<float>> betas(const int zaid) const;
+  std::vector<std::vector<float>> betas(const int zaid);
+
+};  // class EmitDb
+}  // namespace radix
+
+namespace std
+{
+template <>
+struct hash<radix::EmitDBEntry>
+{
+  std::size_t operator()(const radix::EmitDBEntry& rhs) const
+  {
+    using std::hash;
+    using std::size_t;
+    using std::string;
+    return rhs.z + rhs.a + rhs.m + rhs.nLines + rhs.nBetas;
+  }
+};  // struct hash<radix::EmitDBEntry>
+}  // namespace std
+
+#endif /** RADIX_RADIXIO_EMITDB_HH_ */