tkmanager.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 <algorithm>
#include <fstream>
#include <iterator>
 
#include "tkmanager.h"
#include "tklog.h"
#include "tkdatabase.h"
#include "json.hpp"
 
namespace tkn {
 
}
 
using namespace tkn;
using json = nlohmann::json;
 
tkmanager* tkmanager::g_data_manager = nullptr;
 
tkmanager* tkmanager::the_data_manager()
{
    if(g_data_manager) return g_data_manager;
 
    std::lock_guard<std::recursive_mutex> lock(get_mutex());
 
    if(g_data_manager == nullptr) {
        g_data_manager = new tkmanager();
    }
 
    return g_data_manager;
}
 
tkmanager::tkmanager()
{
    std::lock_guard<std::recursive_mutex> lock(get_mutex());
 
    preload_nuclei();
 
    g_data_manager = this;
}
 
tkmanager::~tkmanager()
{
    g_data_manager = nullptr;
}
 
shared_ptr<tklevel_scheme> tkmanager::get_level_scheme(const tkn::tkstring &_nuc, int _zz, int _aa)
{
    if(fmap_of_level_scheme.count(_nuc)) return fmap_of_level_scheme[_nuc];
 
    std::lock_guard<std::recursive_mutex> lock(get_mutex());
 
    auto [it, inserted] = fmap_of_level_scheme.try_emplace(
        _nuc, make_shared<tklevel_scheme>(_nuc,_zz,_aa));
    return it->second;
}
 
void tkmanager::preload_nuclei()
{
    // try with isotope.*, element.*
    gdatabase->begin("isotope.*, element.*","isotope INNER JOIN element on isotope.element_id=element.element_id");
 
    while(gdatabase->next()) {
        tkdb_table& element = (*gdatabase)["ELEMENT"];
        tkdb_table& iso_table = (*gdatabase)["ISOTOPE"];
 
        shared_ptr<tknucleus> anucleus = make_shared<tknucleus>();
 
        tkstring name  = element["symbol"].get_value();
        anucleus->felement_symbol = name;
        anucleus->felement_name = element["name"].get_value();
 
        int zz = element["charge"].get_value().atoi();
        int aa = iso_table["mass"].get_value().atoi();
 
        gdebug << aa << name << " loaded" << do_endl;
 
        name.prepend(iso_table["mass"].get_value());
 
        anucleus->set_name(name.data());
        anucleus->fA = aa;
        anucleus->fZ = zz;
        anucleus->fN = aa - zz;
        anucleus->fis_known = true;
 
        // load element properties
        for (auto &col : element.get_columns()) {
            tkstring colname = col.first;
            if(colname=="element_id") continue;
            if(colname.ends_with("_unit")) continue;
 
            tkstring colname_unit = colname + "_unit";
            tkstring unit("");
            tkstring value = col.second.get_value();
            if(element.get_columns().count(colname_unit)) {
                unit = element.get_columns().at(colname_unit).get_value();
            }
            if( colname == "atomic_mass" ||
                colname == "atomic_radius_van_der_Waals" ||
                colname == "boiling_point" ||
                colname == "density" ||
                colname == "ionization_energy" ||
                colname == "melting_point"
                ) {
                shared_ptr<tkmeasure> data = make_shared<tkmeasure>(value.atof(),unit);
                data->set_info(tkproperty::kKnown);
                data->set_type(colname);
                anucleus->add_property(colname,data);
            }
            else if(colname.begins_with("XRay")) {
                unit = "keV";
                shared_ptr<tkmeasure> data = make_shared<tkmeasure>(value.atof(),unit);
                data->set_info(tkproperty::kKnown);
                data->set_type(colname);
                if(value.atof()>0) anucleus->add_property(colname,data);
            }
            else {
                anucleus->add_property_str(colname,value,unit);
            }
        }
 
        // load isotope properties
        vector<tkstring> fNames = {
            "lifetime",                                                         // From ENSDF
            "abundance", "quadrupoleDeformation", "FY235U", "FY238U", "FY239Pu", "FY241Pu", "FY252Cf", "cFY235U", "cFY238U", "cFY239Pu", "cFY241Pu", "cFY252Cf", // From NUDAT
            "mass_excess", "radius", "magnetic_dipole", "electric_quadrupole"   // From IAEA-NDS API
        };
 
        for(auto &property_name: fNames) {
            tkdb_table::measure_data_struct data_struct;
            iso_table.read_measure(data_struct,property_name);
            if(data_struct.filled) {
                shared_ptr<tkmeasure> data = make_shared<tkmeasure>(data_struct.value,data_struct.unit,data_struct.err);
                if(data_struct.info.contains("ERR=SY")) data->set_info(tkproperty::kSystematic);
                else data->set_info(tkproperty::kKnown);
                data->set_type(property_name);
 
                if(property_name == "lifetime") {
                    if(!data_struct.info.is_empty()) data->set_info_tag(data_struct.info);
                    if(data->get_info_tag().contains("STABLE")) {
                        anucleus->fis_stable = true;
                        data->set_error(0);
                        anucleus->add_property(property_name,data,"STABLE","");
                    }
                    else anucleus->add_property(property_name,data, tkstring::form("%g",data->get_value()));
                }
                else {
                    anucleus->add_property(property_name,data);
                }
            }
        }
 
        anucleus->add_property_str("isotope_year_discovered",iso_table["isotope_year_discovered"].get_value(),"");
 
        tkstring decay_modes = iso_table["decay_modes"].get_value();
        anucleus->add_property_str("decay_modes",decay_modes,"");
 
        tkstring spin_parity_str = iso_table["spin_parity"].get_value();
        double spin = iso_table["spin"].get_value().atof();
        int parity  = iso_table["parity"].get_value().atoi();
 
        anucleus->fspin_parity.set_spin(spin);
        anucleus->fspin_parity.set_parity(parity);
        anucleus->fspin_parity.set_from_str(spin_parity_str);
 
        anucleus->add_property_str("spin_parity",spin_parity_str,"");
 
        fmap_of_nuclei[name] = anucleus;
        fmap_of_nuclei_per_z[anucleus->fZ].push_back(anucleus);
        fmap_of_nuclei_per_a[anucleus->fA].push_back(anucleus);
        fmap_of_nuclei_per_n[anucleus->fN].push_back(anucleus);
        fmap_of_nuclei_per_z_and_a[anucleus->fZ][anucleus->fA] = anucleus;
        fvector_of_nuclei.push_back(anucleus);
 
        if(!fmap_of_symbols.count(anucleus->get_z()))
            fmap_of_symbols[anucleus->get_z()] = anucleus->get_element_symbol();
    }
 
    // separation energy calculation
    if(!get_nucleus(0,1)||!get_nucleus(1,1)||!get_nucleus(2,4)) return;
    auto neutron_ME = fmap_of_nuclei_per_z_and_a.at(0).at(1)->get("mass_excess");
    auto proton_ME = fmap_of_nuclei_per_z_and_a.at(1).at(1)->get("mass_excess");
    auto alpha_ME = fmap_of_nuclei_per_z_and_a.at(2).at(4)->get("mass_excess");
    tkmeasure electron_mass(0.51099895000,"MeV",0.00000000015);
    for(auto &nuc: fmap_of_nuclei) {
        auto ME = nuc.second->get("mass_excess");
        if(!ME) continue;
        // Sn calculation: Sn(A,Z) = − ME(A,Z) + ME(A−1,Z) + ME(n)
        auto daughter = get_nucleus(nuc.second->get_z(),nuc.second->get_a()-1);
        if(daughter && daughter->get("mass_excess")) {
            shared_ptr<tkmeasure> Sn = make_shared<tkmeasure>(-*ME + *daughter->get("mass_excess") + *neutron_ME);
            Sn->set_type("neutronSeparationEnergy");
            nuc.second->add_property("neutronSeparationEnergy",Sn);
        }
        // S2n calculation: S2n(A,Z) = − ME(A,Z) + ME(A−2,Z) + 2*ME(n)
        daughter = get_nucleus(nuc.second->get_z(),nuc.second->get_a()-2);
        if(daughter && daughter->get("mass_excess")) {
            shared_ptr<tkmeasure> S2n = make_shared<tkmeasure>(-*ME + *daughter->get("mass_excess") + 2.*(*neutron_ME));
            S2n->set_type("twoNeutronSeparationEnergy");
            nuc.second->add_property("twoNeutronSeparationEnergy",S2n);
        }
        // Sp calculation: Sp(A,Z) = − ME(A,Z) + ME(A−1,Z-1) + ME(p)
        daughter = get_nucleus(nuc.second->get_z()-1,nuc.second->get_a()-1);
        if(daughter && daughter->get("mass_excess")) {
            shared_ptr<tkmeasure> Sp = make_shared<tkmeasure>(-*ME + *daughter->get("mass_excess") + *proton_ME);
            Sp->set_type("protonSeparationEnergy");
            nuc.second->add_property("protonSeparationEnergy",Sp);
        }
        // S2p calculation: S2p(A,Z) = − ME(A,Z) + ME(A−2,Z-2) + 2*ME(p)
        daughter = get_nucleus(nuc.second->get_z()-2,nuc.second->get_a()-2);
        if(daughter && daughter->get("mass_excess")) {
            shared_ptr<tkmeasure> S2p = make_shared<tkmeasure>(-*ME + *daughter->get("mass_excess") + 2.*(*proton_ME));
            S2p->set_type("twoProtonSeparationEnergy");
            nuc.second->add_property("twoProtonSeparationEnergy",S2p);
        }
 
        // Qalpha
        daughter = get_nucleus(nuc.second->get_z()-2,nuc.second->get_a()-4);
        if(daughter && daughter->get("mass_excess")) {
            shared_ptr<tkmeasure> Q = make_shared<tkmeasure>( *ME - *daughter->get("mass_excess") - *alpha_ME);
            Q->set_type("Qalpha");
            nuc.second->add_property("Qalpha",Q);
        }
        // QbetaMinus
        daughter = get_nucleus(nuc.second->get_z()+1,nuc.second->get_a());
        if(daughter && daughter->get("mass_excess")) {
            shared_ptr<tkmeasure> Q = make_shared<tkmeasure>( *ME - *daughter->get("mass_excess"));
            Q->set_type("QbetaMinus");
            nuc.second->add_property("QbetaMinus",Q);
        }
        // QbetaMinusOneNeutronEmission
        daughter = get_nucleus(nuc.second->get_z()+1,nuc.second->get_a()-1);
        if(daughter && daughter->get("mass_excess")) {
            shared_ptr<tkmeasure> Q = make_shared<tkmeasure>( *ME - *daughter->get("mass_excess") - *neutron_ME);
            Q->set_type("QbetaMinusOneNeutronEmission");
            nuc.second->add_property("QbetaMinusOneNeutronEmission",Q);
        }
        // QbetaMinusTwoNeutronEmission
        daughter = get_nucleus(nuc.second->get_z()+1,nuc.second->get_a()-2);
        if(daughter && daughter->get("mass_excess")) {
            shared_ptr<tkmeasure> Q = make_shared<tkmeasure>( *ME - *daughter->get("mass_excess") - 2.*(*neutron_ME) );
            Q->set_type("QbetaMinusTwoNeutronEmission");
            nuc.second->add_property("QbetaMinusTwoNeutronEmission",Q);
        }
        // QdeltaAlpha: 0.5*(Qa(Z+2,N+2)-Qa(Z,N))
        daughter = get_nucleus(nuc.second->get_z()+2,nuc.second->get_a()+4);
        if(daughter && daughter->get("mass_excess") && nuc.second->get("Qalpha")) {
            auto Qa = nuc.second->get("Qalpha");
            auto Qa2 = *daughter->get("mass_excess") - *ME - *alpha_ME;
            shared_ptr<tkmeasure> Q = make_shared<tkmeasure>( 0.5*(Qa2 - *Qa ));
            Q->set_type("QdeltaAlpha");
            nuc.second->add_property("QdeltaAlpha",Q);
        }
        // QdoubleBetaMinus
        daughter = get_nucleus(nuc.second->get_z()+2,nuc.second->get_a());
        if(daughter && daughter->get("mass_excess")) {
            shared_ptr<tkmeasure> Q = make_shared<tkmeasure>( *ME - *daughter->get("mass_excess"));
            Q->set_type("QdoubleBetaMinus");
            nuc.second->add_property("QdoubleBetaMinus",Q);
        }
        // QelectronCapture
        daughter = get_nucleus(nuc.second->get_z()-1,nuc.second->get_a());
        if(daughter && daughter->get("mass_excess")) {
            shared_ptr<tkmeasure> Q = make_shared<tkmeasure>( *ME - *daughter->get("mass_excess"));
            Q->set_type("QelectronCapture");
            nuc.second->add_property("QelectronCapture",Q);
        }
        // QdoubleElectronCapture
        daughter = get_nucleus(nuc.second->get_z()-2,nuc.second->get_a());
        if(daughter && daughter->get("mass_excess")) {
            shared_ptr<tkmeasure> Q = make_shared<tkmeasure>( *ME - *daughter->get("mass_excess"));
            Q->set_type("QdoubleElectronCapture");
            nuc.second->add_property("QdoubleElectronCapture",Q);
        }
        // QelectronCaptureOneProtonEmission
        daughter = get_nucleus(nuc.second->get_z()-2,nuc.second->get_a()-1);
        if(daughter && daughter->get("mass_excess")) {
            shared_ptr<tkmeasure> Q = make_shared<tkmeasure>( *ME - *daughter->get("mass_excess") - *proton_ME);
            Q->set_type("QelectronCaptureOneProtonEmission");
            nuc.second->add_property("QelectronCaptureOneProtonEmission",Q);
        }
        // QpositronEmission
        daughter = get_nucleus(nuc.second->get_z()-1,nuc.second->get_a());
        if(daughter && daughter->get("mass_excess")) {
            auto q_positron = make_shared<tkmeasure>( *ME - *daughter->get("mass_excess") - 2.*electron_mass);
            q_positron->set_type("QpositronEmission");
            nuc.second->add_property("QpositronEmission",q_positron);
        }
 
        // binding energy
        auto q_binding = make_shared<tkmeasure>( 1./nuc.second->get_a()*(nuc.second->get_z()*(*proton_ME) + nuc.second->get_n()*(*neutron_ME) - *ME));
        q_binding->set_type("binding_energy_overA");
        nuc.second->add_property("binding_energy_overA",q_binding);
 
        // pairingGap PG=0.5 x (-1)N x ( 2BE(Z,N) - BE(Z,N-1) -BE(Z,N+1) )
        auto nuc_minus = get_nucleus(nuc.second->get_z(),nuc.second->get_a()-1);
        auto nuc_plus = get_nucleus(nuc.second->get_z(),nuc.second->get_a()+1);
        if(nuc_minus && nuc_plus && nuc_minus->get("mass_excess") && nuc_plus->get("mass_excess")) {
            auto q_pairing = make_shared<tkmeasure>( 0.5 * pow(-1,nuc.second->get_n()) * ( -2* (*ME) + *nuc_minus->get("mass_excess") + *nuc_plus->get("mass_excess")));
            q_pairing->set_type("pairingGap");
            nuc.second->add_property("pairingGap",q_pairing);
        }
    }
}
 
void tkmanager::preload_level_schemes(bool _verbose)
{
    size_t nnuc=fvector_of_nuclei.size(),inuc=0;
    for(const auto &nuc : fvector_of_nuclei){
        if(_verbose) glog.progress_bar(nnuc,inuc,"... pre-loading");
        get_level_scheme(nuc->get_symbol(),nuc->get_z(),nuc->get_a());
        inuc++;
    }
}

vector<shared_ptr<tknucleus> > tkmanager::get_nuclei(const std::function<bool (shared_ptr<tknucleus>)> &_selection)
{
    vector<shared_ptr<tknucleus>> res;
    std::copy_if(fvector_of_nuclei.begin(), fvector_of_nuclei.end(), std::back_inserter(res),
                 [&_selection](const auto &nuc) { return _selection(nuc); });
    return res;
}
 
bool tkmanager::known_element(tkstring _nuc, int &_z) {
    tkstring symbol=_nuc.extract_alpha();
    for(auto &nuc_test: fmap_of_nuclei) {
        if(nuc_test.first.copy().extract_alpha() == symbol) {
            _z = nuc_test.second->get_z();
            return true;
        }
    }
    return false;
}
 
int tkmanager::get_new_level_id()
{
    std::lock_guard<std::recursive_mutex> lock(get_mutex());
    fmax_level_id++;
    return fmax_level_id;
}
 
int tkmanager::get_new_decay_id()
{
    std::lock_guard<std::recursive_mutex> lock(get_mutex());
    fmax_decay_id++;
    return fmax_decay_id;
}

 *   - `S1p`, `S1p_sigma`
 *   - `S2p`, `S2p_sigma`
 *
 * @param json_path Path to the JSON file (default: "%s/databases/driplines.json").
 *
 * @note After loading, all points are sorted by (Z,N) for deterministic order.
 *
 * @note Data extracted from the supplementary material of:
 *       L. Neufcourt, Y. Cao, S. A. Giuliani, W. Nazarewicz, E. Olsen, O. B. Tarasov,
 *       "Quantified limits of the nuclear landscape", Phys. Rev. C 101, 044307 (2020).
 *
 * @code
 * // Example usage:
 * gmanager->load_drip_lines("path/to/nuclei_quantities.json");
 * auto s2n = gmanager->get_drip_line("S2n");
 * @endcode
 */
void tkmanager::load_drip_lines(const tkstring& json_path)
{
    std::ifstream f(json_path);
    if (!f.is_open()) {
        glog << error << "tkmanager::load_drip_lines: cannot open " << json_path << do_endl;
        return;
    }
 
    json j;
    try {
        f >> j;
    } catch (const std::exception& e) {
        glog << error << "tkmanager::load_drip_lines: JSON parse error in "
             << json_path << " : " << e.what() << do_endl;
        return;
    }
 
    if (!j.is_array()) {
        glog << error << "tkmanager::load_drip_lines: root is not an array in "
             << json_path << do_endl;
        return;
    }
 
    fdrip_lines.clear();
 
    auto push_if_present = [&](const json& o, const char* qname, const char* sname){
        if (!o.contains(qname)) return;
        if (o[qname].is_null()) return;
 
        tkn_drip_point p;
        p.Z = o.value("Z", 0);
        p.N = o.value("N", 0);
        p.value_mev = o.value(qname, std::numeric_limits<double>::quiet_NaN());
        // sigma may be absent → keep NaN
        if (o.contains(sname) && !o[sname].is_null()) {
            p.sigma_mev = o.value(sname, std::numeric_limits<double>::quiet_NaN());
        }
        p.type = qname;
 
        fdrip_lines[p.type].push_back(std::move(p));
    };
 
    try {
        for (const auto& o : j) {
            if (!o.is_object()) continue;
            // Z and N must exist; skip otherwise
            if (!o.contains("Z") || !o.contains("N")) continue;
 
            push_if_present(o, "S1n", "S1n_sigma");
            push_if_present(o, "S2n", "S2n_sigma");
            push_if_present(o, "S1p", "S1p_sigma");
            push_if_present(o, "S2p", "S2p_sigma");
        }
    } catch (const std::exception& e) {
        glog << error << "tkmanager::load_drip_lines: error while iterating JSON: "
             << e.what() << do_endl;
        fdrip_lines.clear();
        return;
    }
 
    // Sort deterministically by Z, then N
    for (auto& kv : fdrip_lines) {
        auto& vec = kv.second;
        std::sort(vec.begin(), vec.end(), [](const tkn_drip_point& a, const tkn_drip_point& b){
            if (a.Z != b.Z) return a.Z < b.Z;
            return a.N < b.N;
        });
    }
 
    glog << info << "load_drip_lines: loaded " << j.size()
         << " nuclei objects from " << json_path << do_endl;
}
 

std::vector<tkn_drip_point> tkmanager::get_drip_line(const std::string& _type, bool reduce_per_Z)
{
    if (fdrip_lines.empty()) {
        // charge par défaut si pas déjà fait
        load_drip_lines(tkstring::Form("%s/databases/driplines.json", TKN_SYS));
    }
 
    auto it = fdrip_lines.find(_type);
    if (it == fdrip_lines.end()) {
        glog << warning << "get_drip_line: unknown type '" << _type << "'" << do_endl;
        return {};
    }
 
    const auto& vec = it->second;
 
    if (!reduce_per_Z) {
        return vec;
    }
 
    // --- Réduction : garder seulement le point avec min |value| par Z
    std::unordered_map<int, tkn_drip_point> best_by_Z;
    for (const auto& p : vec) {
        auto itZ = best_by_Z.find(p.Z);
        if (itZ == best_by_Z.end()) {
            best_by_Z[p.Z] = p;
        } else {
            const auto& cur = itZ->second;
            if (std::isnan(cur.value_mev) ||
                (!std::isnan(p.value_mev) && std::abs(p.value_mev) < std::abs(cur.value_mev))) {
                best_by_Z[p.Z] = p;
            }
        }
    }
 
    std::vector<tkn_drip_point> reduced;
    reduced.reserve(best_by_Z.size());
    std::transform(best_by_Z.begin(), best_by_Z.end(), std::back_inserter(reduced),
                   [](const auto &kv) { return kv.second; });
 
    std::sort(reduced.begin(), reduced.end(), [](const tkn_drip_point& a, const tkn_drip_point& b){
        if (a.Z != b.Z) return a.Z < b.Z;
        return a.N < b.N;
    });
 
    return reduced;
}
 
#ifdef HAS_ROOT
ClassImp(tkmanager)
#endif