tkisotope_builder.cpp

/********************************************************************************
 *   Copyright (c) : Université de Lyon 1, CNRS/IN2P3, UMR5822,                 *
 *                   IP2I, F-69622 Villeurbanne Cedex, France                   *
 *                   Normandie Université, ENSICAEN, UNICAEN, CNRS/IN2P3,       *
 *                   LPC Caen, F-14000 Caen, France                             *
 *   Contibutor(s) :                                                            *
 *      Jérémie Dudouet jeremie.dudouet@cnrs.fr [2020]                          *
 *      Diego Gruyer    diego.gruyer@cnrs.fr    [2020]                          *
 *                                                                              *
 *   Licensed under the MIT License <http://opensource.org/licenses/MIT>.       *
 *   SPDX-License-Identifier: MIT                                               *
 ********************************************************************************/
 
#include "tkisotope_builder.h"
 
#include <algorithm>
#include <iostream>
#include <fstream>
 
#include "tklog.h"
 
namespace tkn {
}
 
using namespace tkn;
using namespace std;
 
// Decay Mode           : decay modes.
// Mass excess          : Mass excess
// β2                   : quadrupole deformation parameter obtained from the B(E2) value, as explained in the 2001 work of Raman et al. Basically,
//                        β2=(4π/3ZR0^2)[B(E2)]^(1/2), where R0=1.2A^(1/3), Z and A are the number of protons and the number of nucleons respectively,
//                        and the B(E2) values are obtained from ENSDF. B(E2) values can have large uncertainties, or be defined as upper/lower limits.
//                        The cell will not be colored for the latter cases and furthermore, it will not be colored if the uncertainty of the derived β2 exceeds 50%.
// FY235U               : thermal neutron-induced fission yields for 235-Uranium. The fission yields are normalized so that their sum is equal to 2.
// FY238U               : neutron-induced fission yields for 238-Uranium. The fission yields are normalized so that their sum is equal to 2.
// FY239Pu              : thermal neutron-induced fission yields for 239-Plutonium. The fission yields are normalized so that their sum is equal to 2.
// FY241Pu              : neutron-induced fission yields for 241-Plutonium. The fission yields are normalized so that their sum is equal to 2.
// FY252Cf              : spontaneous fission yields for 252-Californium. The fission yields are normalized so that their sum is equal to 2.
// cFY235U              : cumulative fission yields
// cFY238U              : cumulative fission yields
// cFY239Pu             : cumulative fission yields
// cFY241Pu             : cumulative fission yields
// cFY252Cf             : cumulative fission yields
 
tkisotope_builder::tkisotope_builder(tkdatabase *_database, const char *_table_name) : tkdb_builder(_database, _table_name)
{
   fTable = &fDataBase->new_table(fTableName);
 
   fTable->add_column("isotope_id", "INT PRIMARY KEY NOT NULL");
   fTable->add_column("element_id", "INT NOT NULL");
   fTable->add_column("mass", "INT");
 
   // data taken from IAEA API
 
   fTable->add_column("isotope_year_discovered", "INT");
 
   add_measure("mass_excess", false);
   add_measure("electric_quadrupole", false);
   add_measure("magnetic_dipole", false);
   add_measure("radius", false);
 
   // data taken from NUDAT3 json file
 
   add_measure("abundance", false);
   add_measure("quadrupoleDeformation", false);
 
   add_measure("FY235U", false);
   add_measure("FY238U", false);
   add_measure("FY239Pu", false);
   add_measure("FY241Pu", false);
   add_measure("FY252Cf", false);
 
   add_measure("cFY235U", false);
   add_measure("cFY238U", false);
   add_measure("cFY239Pu", false);
   add_measure("cFY241Pu", false);
   add_measure("cFY252Cf", false);
 
   fTable->add_column("decay_modes", "TEXT");
 
   // columns created here but filled later by ensdf values
 
   fTable->add_column("spin", "REAL");
   fTable->add_column("parity", "INT");
   fTable->add_column("spin_parity", "TEXT");
 
   add_measure("lifetime", false);
 
   fTable->add_constraint("FOREIGN KEY (element_id) REFERENCES ELEMENT (element_id)"); // hold this link with a constraints
   fTable->write_to_database();
 
   glog << info << "Creating '" << _table_name << "' table" << do_endl;
}
 
tkisotope_builder::~tkisotope_builder() = default;
 
void tkisotope_builder::fill_isotope(const json_nucleus &_json_nuc)
{
   fIsotopeIdx++;
 
   //    glog << debug << "filling db with : " << _json_nuc.name.data() << " " << fIsotopeIdx<< do_endl;
 
   (*fTable)["isotope_id"].set_value(fIsotopeIdx);
   (*fTable)["element_id"].set_value(_json_nuc.Z + 1);
   (*fTable)["mass"].set_value(_json_nuc.A);
   (*fTable)["isotope_year_discovered"].set_value(_json_nuc.year);
 
   fill_measure("mass_excess", _json_nuc.mass_excess, false);
   fill_measure("abundance", _json_nuc.abundance, false);
   fill_measure("quadrupoleDeformation", _json_nuc.quadrupoleDeformation, false);
 
   fill_measure("FY235U", _json_nuc.FY235U, false);
   fill_measure("FY238U", _json_nuc.FY238U, false);
   fill_measure("FY239Pu", _json_nuc.FY239Pu, false);
   fill_measure("FY241Pu", _json_nuc.FY241Pu, false);
   fill_measure("FY252Cf", _json_nuc.FY252Cf, false);
 
   fill_measure("cFY235U", _json_nuc.cFY235U, false);
   fill_measure("cFY238U", _json_nuc.cFY238U, false);
   fill_measure("cFY239Pu", _json_nuc.cFY239Pu, false);
   fill_measure("cFY241Pu", _json_nuc.cFY241Pu, false);
   fill_measure("cFY252Cf", _json_nuc.cFY252Cf, false);
 
   fill_measure("electric_quadrupole", _json_nuc.electric_quadrupole, false);
   fill_measure("magnetic_dipole", _json_nuc.magnetic_dipole, false);
   fill_measure("radius", _json_nuc.radius, false);
 
   // fill the neutron lifetime
   if (_json_nuc.Z == 0 && _json_nuc.N == 1) {
      tkdb_table::measure_data_struct lifetime;
      lifetime.value = 613.9;
      lifetime.err = 0.6;
      lifetime.unit = "S";
      lifetime.filled = true;
      fill_measure("lifetime", lifetime, false);
   }
 
   tkstring decays;
   for (auto &dec : _json_nuc.decay_modes) decays.append(tkstring::form("[%s;%g]", dec.first.data(), dec.second));
   (*fTable)["decay_modes"].set_value(decays.data());
 
   fTable->push_row();
}
 
int tkisotope_builder::fill_database(const char *_json_filename, const char *_gs_filename, int _only_charge, int _only_mass)
{
   glog.set_class("db_isotope_builder");
   glog.set_method(tkstring::form("fill_database('%s','%s')", _json_filename, _gs_filename));
 
   tkdb_table &fTable = (*fDataBase)[fTableName.data()];
 
   tkstring message = "Building '" + fTable.get_name() + "' table";
 
   std::ifstream ifs(_json_filename);
   if (!ifs.good()) {
      glog << error_v << "data file cannot be opened" << do_endl;
      glog.clear();
      return 1;
   }
 
   json jf = json::parse(ifs);
 
   auto &nuclides = jf.at("nuclides");
   int NNucs = nuclides.size();
   int inuc = 0;
 
   map<int, map<int, json_nucleus>> nuclei;
   tkstring nread;
 
   for (auto &nuc : nuclides) {
      json_nucleus json_nuc = read_nucleus(nuc);
 
      bool do_push = ((_only_charge == 0) && (_only_mass == 0));
      if ((_only_charge == json_nuc.Z) && (_only_mass == 0)) do_push = true;
      if ((_only_charge == json_nuc.Z) && (_only_mass == json_nuc.A)) do_push = true;
      if ((_only_charge == 0) && (_only_mass == json_nuc.A)) do_push = true;
 
      if (do_push) {
         nuclei[json_nuc.Z][json_nuc.A] = json_nuc;
         nread += tkstring::form("(%d,%d)", json_nuc.Z, json_nuc.A);
      }
      glog.progress_bar(NNucs, inuc, message.data());
      inuc++;
   }
 
   tkn::tkstring line = "";
   fstream file(_gs_filename, ios::in);
   if (!file.is_open()) {
      glog << error_v << "IAEA ground state file cannot be opened: '" << _gs_filename << "'" << do_endl;
      glog << error_v << "mass_excess, radius, electric_quadrupole and magnetic_dipole will be missing from the database." << do_endl;
      glog << error_v << "Download it manually from: https://nds.iaea.org/relnsd/v1/data?fields=ground_states&nuclides=all" << do_endl;
      glog << error_v << "and place it at: " << _gs_filename << do_endl;
      glog.clear();
      return 1;
   } else {
      // read the first line to fill the <number,name> column map
      if (!getline(file, line)) return 1;
      int ic = 0;
 
      std::map<int, tkn::tkstring> gs_columns;
      for (const auto &cname : line.tokenize(",")) gs_columns[ic++] = cname;
 
      // now read all lines and fill the <name, value> map
      while (getline(file, line)) {
         std::map<tkn::tkstring, tkn::tkstring> values;
         line.replace_all(",,", ", ,");
         line.replace_all(",,", ", ,");
         ic = 0;
         for (auto val : line.tokenize(",")) {
            val.remove_all(" ");
            values[gs_columns[ic]] = val;
            ic++;
         }
         tkstring dum = tkstring::form("(%d,%d)", values["z"].atoi(), values["n"].atoi() + values["z"].atoi());
 
         if (!nread.contains(tkstring::form("(%d,%d)", values["z"].atoi(), values["n"].atoi() + values["z"].atoi()))) continue;
         auto &nuc = nuclei[values["z"].atoi()][values["n"].atoi() + values["z"].atoi()];
         if (!values["massexcess"].is_empty()) nuc.mass_excess = {values["massexcess"].atof(), values["unc_me"].atof(), -1., -1., "keV", values["me_systematics"].equal_to("Y") ? ";ERR=SY" : "", true};
         if (!values["radius"].is_empty()) nuc.radius = {values["radius"].atof(), values["unc_r"].atof(), -1., -1., "fm", "", true};
         if (!values["electric_quadrupole"].is_empty()) nuc.electric_quadrupole = {values["electric_quadrupole"].atof(), values["unc_eq"].atof(), -1., -1., "barn", "", true};
         if (!values["magnetic_dipole"].is_empty()) nuc.magnetic_dipole = {values["magnetic_dipole"].atof(), values["unc_md"].atof(), -1., -1., "mun", "", true};
         if (!values["discovery"].is_empty()) nuc.year = values["discovery"].atoi();
      }
   }
 
   fDataBase->exec_sql("BEGIN TRANSACTION");
 
   inuc = 0;
   for (auto &z : nuclei) {
      for (auto &a : z.second) {
         //            glog << debug << "filling db with : " << a.second.name.data() << do_endl;
 
         fill_isotope(a.second);
         //            json_nuc.print();
         fDataBase->exec_sql(tkstring::form("CREATE VIEW nuc%s as select * from isotope INNER JOIN element on isotope.element_id=element.element_id where element.charge=%d AND isotope.mass=%d", a.second.name.data(), a.second.Z, a.second.A));
         glog.progress_bar(NNucs, inuc, message.data());
         inuc++;
      }
   }
 
   fDataBase->exec_sql("END TRANSACTION");
 
   fDataBase->exec_sql("CREATE INDEX element_index ON isotope(element_id)");
   fDataBase->exec_sql("CREATE INDEX mass_index ON isotope(mass)");
 
   glog.clear();
 
   return 0;
}
 
void tkisotope_builder::read_measure(const tkstring &key, const json &entry, tkdb_table::measure_data_struct &_data, bool _fromstd)
{
   auto to_number = [](const json &value) {
      if (value.is_string()) return value.get<tkstring>().atof();
      if (value.is_number()) return value.get<double>();
      return 0.0;
   };
   bool value_set = false;
   for (auto &_item : entry.items()) {
      const auto &_key = _item.key();
      const auto &_val = _item.value();
      if (_key == "formats") continue;
      if (_key == "uncertainty") {
         if (_val.is_object()) {
            const auto type = _val.value("type", "");
            if (type == "symmetric")
               _data.err = to_number(_val.at("value"));
            else if (type == "range") {
               if (_val.contains("upperLimit") && _val.contains("lowerLimit")) {
                  const double upper = to_number(_val.at("upperLimit"));
                  const double lower = to_number(_val.at("lowerLimit"));
                  if (!value_set) _data.value = 0.5 * (upper + lower);
                  _data.err = 0.0;
                  _data.err_low = _data.value - lower;
                  _data.err_high = upper - _data.value;
               }
            }
         } else {
            _data.err = to_number(_val);
         }
      } else if (_key == "value") {
         _data.value = to_number(_val);
         value_set = true;
      } else if (_key == "unit")
         _data.unit = _val.get<tkstring>();
      else if (_key == "evaluatorInput")
         continue;
      else if (_key == "isSystematics" && (_val == "true" || (_val.is_boolean() && _val.get<bool>())))
         _data.info.append(";ERR=SY");
      else {
         glog << error << key << ": unkown key: " << _key << " type: " << _item.value().type_name() << do_endl;
         cout << _val << endl;
      }
   }
   if (_fromstd) {
      // extract value and err from nudat formated values
      if (entry.contains("formats") && entry.at("formats").contains("STD")) {
         tkstring format = entry.at("formats").at("STD").get<tkstring>();
         auto tokens = format.tokenize();
         _data.value = tokens.front().atof();
         if (tokens.size() == 2) {
            _data.err = tokens.back().atof();
         }
      }
   }
   _data.filled = true;
}
 
tkisotope_builder::json_nucleus tkisotope_builder::read_nucleus(json &entry)
{
   json_nucleus nuc;
   bool pause = false;
   for (auto &item : entry.items()) {
      auto &key = item.key();
      auto &val = item.value();
      if (key == "name")
         nuc.name = val;
      else if (key == "z")
         nuc.Z = val;
      else if (key == "n")
         nuc.N = val;
      else if (key == "a")
         nuc.A = val;
      else if (key == "bindingEnergy")
         continue;
      else if (key == "bindingEnergyLDMFit")
         continue;
      else if (key == "levels") {
         // pick the ground state level when isomers are present
         auto is_ground_state = [](const json &level) {
            if (!level.contains("energy")) return false;
            const auto &energy = level.at("energy");
            if (!energy.is_object() || !energy.contains("value")) return false;
            const auto &energy_val = energy.at("value");
            if (energy_val.is_string()) return energy_val.get<tkstring>().atof() == 0.0;
            if (energy_val.is_number()) return energy_val.get<double>() == 0.0;
            return false;
         };
         const json *level_ptr = nullptr;
         const auto it = std::find_if(val.begin(), val.end(), is_ground_state);
         bool picked_ground = (it != val.end());
         if (picked_ground)
            level_ptr = &(*it);
         else if (!val.empty())
            level_ptr = &val.front();
         if (!level_ptr) continue;
         for (auto &_item : level_ptr->items()) {
            const auto &_key = _item.key();
            const auto &_val = _item.value();
            if (_key == "decayModes") {
               if (_val.is_object() && (_val.contains("observed") || _val.contains("predicted"))) {
                  if (_val.contains("observed")) {
                     for (const auto &mode : _val.at("observed")) {
                        tkstring type = mode.value("mode", "");
                        tkstring ratio = "-1";
                        if (mode.contains("value")) {
                           if (mode.at("value").is_string())
                              ratio = mode.at("value").get<tkstring>();
                           else if (mode.at("value").is_number())
                              ratio = tkstring::form("%g", mode.at("value").get<double>());
                        }
                        nuc.decay_modes.emplace_back(type, ratio.atof());
                     }
                  }
                  if (_val.contains("predicted")) {
                     for (const auto &mode : _val.at("predicted")) {
                        tkstring type = mode.value("mode", "");
                        tkstring ratio = "-1";
                        nuc.decay_modes.emplace_back(type, ratio.atof());
                     }
                  }
               } else {
                  for (auto &mode : _val.items()) {
                     if (mode.value().size() == 3) {
                        tkstring type = mode.value().at("name");
                        tkstring symbol = mode.value().at("symbol");
                        tkstring ratio;
                        // unkown ratios
                        if (symbol == "?") {
                           symbol = "=";
                           ratio = "-1";
                        } else
                           ratio = mode.value().at("value");
                        nuc.decay_modes.emplace_back(type, ratio.atof());
                     } else {
                        tkstring type = mode.value().at("name");
                        type.remove_all("?");
                        tkstring symbol = "=";
                        tkstring ratio = "-1";
                        nuc.decay_modes.emplace_back(type, ratio.atof());
                     }
                  }
               }
            } else if (_key == "energy") {
               if (_val.empty()) continue;
               if (_val.size() < 2) continue;
               double ener = 0.0;
               const auto &energy_val = _val.at("value");
               if (energy_val.is_string())
                  ener = energy_val.get<tkstring>().atof();
               else if (energy_val.is_number())
                  ener = energy_val.get<double>();
               if (ener != 0. && picked_ground) {
                  glog << error << nuc.name << ": level is not ground state ! " << _val << do_endl;
               }
            } else if (_key == "massExcess")
               continue; // mass_excess is sourced exclusively from the IAEA CSV
            else if (_key == "abundance") {
               read_measure(_key, _val, nuc.abundance, true);
               nuc.abundance.unit = "%";
            }
            // include in the ensdf parsing
            else if (_key == "halflife" || _key == "spinAndParity" || _key == "spinParity")
               continue;
            else {
               glog << error << "levels: unkown key: " << _key << " type: " << _item.value().type_name() << do_endl;
               cout << _val << endl;
            }
         }
      } else if (key == "pairingGap")
         continue;
      else if (key == "qValues")
         continue;
      else if (key == "separationEnergies")
         continue;
      else if (key == "webLinks")
         continue;
      else if (key == "BE4DBE2")
         continue;
      else if (key == "excitedStateEnergies") {
         for (auto &_item : val.items()) {
            const auto &_key = _item.key();
            const auto &_val = _item.value();
            if (_key == "firstExcitedStateEnergy") continue;
            if (_key == "firstTwoPlusEnergy") continue;
            if (_key == "firstFourPlusEnergy") continue;
            if (_key == "firstFourPlusOverFirstTwoPlusEnergy") continue;
            if (_key == "firstThreeMinusEnergy") continue;
            glog << error << key << ": unkown key: " << _key << " type: " << _item.value().type_name() << do_endl;
            cout << _val << endl;
         }
      } else if (key == "fissionYields") {
         for (auto &_item : val.items()) {
            const auto &_key = _item.key();
            const auto &_val = _item.value();
            if (_key == "FY235U")
               read_measure(_key, _val, nuc.FY235U);
            else if (_key == "FY238U")
               read_measure(_key, _val, nuc.FY238U);
            else if (_key == "FY239Pu")
               read_measure(_key, _val, nuc.FY239Pu);
            else if (_key == "FY241Pu")
               read_measure(_key, _val, nuc.FY241Pu);
            else if (_key == "FY252Cf")
               read_measure(_key, _val, nuc.FY252Cf);
            else if (_key == "cFY235U")
               read_measure(_key, _val, nuc.cFY235U);
            else if (_key == "cFY238U")
               read_measure(_key, _val, nuc.cFY238U);
            else if (_key == "cFY239Pu")
               read_measure(_key, _val, nuc.cFY239Pu);
            else if (_key == "cFY241Pu")
               read_measure(_key, _val, nuc.cFY241Pu);
            else if (_key == "cFY252Cf")
               read_measure(_key, _val, nuc.cFY252Cf);
            else {
               glog << error << key << ": unkown key: " << _key << " type: " << _item.value().type_name() << do_endl;
               cout << _val << endl;
            }
         }
      } else if (key == "quadrupoleDeformation")
         read_measure(key, val, nuc.quadrupoleDeformation);
      else if (key == "thermalNeutronCapture")
         continue;
      else if (key == "thermalNeutronFission")
         continue;
      else if (key == "discovery")
         continue;
      else {
         glog << warning << "unknown key: " << key << do_endl;
         cout << "value: " << endl;
         cout << val << endl;
         pause = true;
      }
   }
   if (pause) {
      nuc.print();
      cin.get();
   }
   return nuc;
}
 
#ifdef HAS_ROOT
ClassImp(tkisotope_builder);
#endif