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]                          *
 *                                                                              *
 *    This software is governed by the CeCILL-B license under French law and    *
 *    abiding by the  rules of distribution of free  software.  You can use,    *
 *    modify  and/ or  redistribute  the  software under  the  terms of  the    *
 *    CeCILL-B license as circulated by CEA, CNRS and INRIA at the following    *
 *    URL \"http://www.cecill.info\".                                           *
 *                                                                              *
 *    As a counterpart to the access  to the source code and rights to copy,    *
 *    modify  and redistribute granted  by the  license, users  are provided    *
 *    only with a limited warranty  and the software's author, the holder of    *
 *    the economic  rights, and the  successive licensors have  only limited    *
 *    liability.                                                                *
 *                                                                              *
 *    In this respect, the user's attention is drawn to the risks associated    *
 *    with loading,  using, modifying  and/or developing or  reproducing the    *
 *    software by the user in light of its specific status of free software,    *
 *    that  may mean that  it is  complicated to  manipulate, and  that also    *
 *    therefore  means that it  is reserved  for developers  and experienced    *
 *    professionals having in-depth  computer knowledge. Users are therefore    *
 *    encouraged  to load  and test  the software's  suitability  as regards    *
 *    their  requirements  in  conditions  enabling the  security  of  their    *
 *    systems  and/or data to  be ensured  and, more  generally, to  use and    *
 *    operate it in the same conditions as regards security.                    *
 *                                                                              *
 *    The fact that  you are presently reading this means  that you have had    *
 *    knowledge of the CeCILL-B license and that you accept its terms.          *
 ********************************************************************************/
 
#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 239-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.FY235U,false);
    fill_measure("FY239Pu",_json_nuc.FY239Pu,false);
    fill_measure("FY241Pu",_json_nuc.FY239Pu,false);
    fill_measure("FY252Cf",_json_nuc.FY252Cf,false);
 
    fill_measure("cFY235U",_json_nuc.cFY235U,false);
    fill_measure("cFY238U",_json_nuc.cFY235U,false);
    fill_measure("cFY239Pu",_json_nuc.cFY239Pu,false);
    fill_measure("cFY241Pu",_json_nuc.cFY239Pu,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())
    {
        // 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") read_measure(_key, _val, nuc.mass_excess);
                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