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 "tkmanager.h"
#include "tklog.h"
#include "tkdatabase.h"
 
namespace tkn {
}
 
using namespace tkn;
 
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());
 
    if(fmap_of_level_scheme.count(_nuc)==0) {
        fmap_of_level_scheme[_nuc] = make_shared<tklevel_scheme>(_nuc,_zz,_aa);
    }
 
    return fmap_of_level_scheme[_nuc];
}
 
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;
 
        // 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")) {
            shared_ptr<tkmeasure> Q = make_shared<tkmeasure>( *ME - *daughter->get("mass_excess") - 2.*electron_mass);
            Q->set_type("QpositronEmission");
            nuc.second->add_property("QpositronEmission",Q);
        }
 
        // binding energy
        shared_ptr<tkmeasure> Q = make_shared<tkmeasure>( 1./nuc.second->get_a()*(nuc.second->get_z()*(*proton_ME) + nuc.second->get_n()*(*neutron_ME) - *ME));
        Q->set_type("binding_energy_overA");
        nuc.second->add_property("binding_energy_overA",Q);
 
        // 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")) {
            shared_ptr<tkmeasure> Q = make_shared<tkmeasure>( 0.5 * pow(-1,nuc.second->get_n()) * ( -2* (*ME) + *nuc_minus->get("mass_excess") + *nuc_plus->get("mass_excess")));
            Q->set_type("pairingGap");
            nuc.second->add_property("pairingGap",Q);
        }
    }
}
 
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;
    for(const auto &nuc : fvector_of_nuclei)
        if(_selection(nuc)) res.push_back(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;
}
 
#ifdef HAS_ROOT
ClassImp(tkmanager)
#endif