tklevel_scheme.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.          *
 ********************************************************************************/
 
#ifdef HAS_DEBUG
#include <chrono>
#endif
 
#include "tkdb_table.h"
#include "tklevel_scheme.h"
#include "tkstring.h"
#include "tklog.h"
#include "tkdatabase.h"
#include "tkmanager.h"
 
namespace tkn {
}
 
using namespace tkn;
 
tklevel_scheme::tklevel_scheme(const tkstring &_nuc, int _zz, int _aa):fnucleus(_nuc), fA(_aa), fZ(_zz)
{
    if(fA<0||fZ<0) return;
    init();
}
 
std::vector<shared_ptr<tklevel>> tklevel_scheme::get_levels(const std::function<bool (shared_ptr<tklevel>)> &_selection)
{
    vector<shared_ptr<tklevel>> res;
    for(const auto &lvl : get_levels()) if(_selection(lvl)) res.push_back(lvl);
    return res;
}
 
std::vector<shared_ptr<tkdecay>> tklevel_scheme::get_decays(const std::function<bool (shared_ptr<tkdecay>)> &_selection)
{
    vector<shared_ptr<tkdecay>> res;
    for(const auto &dec : get_decays()) if(_selection(dec)) res.push_back(dec);
    return res;
}
 
 
void tklevel_scheme::print(const tkstring &_data, const tkstring &_option)
{
    if(_data.is_empty()) {
        glog << info << " Level scheme of " << fnucleus << do_endl;
        glog << info << " Options : 'dataset', 'level', 'decay'" << do_endl;
    }
    if(_data.contains("dataset")){
        glog << info << "Available datasets are : " << do_endl;
        for(auto &dts : fmap_of_dataset) glog << info << dts.second->get_name() << " (" << dts.first << ")" << do_endl;
        glog << comment << "Current dataset is '" << get_dataset()->get_name() << "'" << " (" << get_dataset()->get_id() << ")" << do_endl;
    }
    if(_data.contains("level")||_data.contains("decay")){
        fmap_of_dataset[fdatasetid]->print(_data, _option);
    }
    if(_data == "*") {
        int _current_data_set = fdatasetid;
        for(auto &dataset: fmap_of_dataset) {
            glog << info << "dataset : " << dataset.second->get_name() << do_endl;
            select_dataset(dataset.first);
            dataset.second->print("leveldecay");
        }
        // to select back the current dataset
        select_dataset(_current_data_set);
    }
}
 
bool tklevel_scheme::select_dataset(const tkstring &_dataset_name)
{
    if(fmap_of_dataset_name_id.count(_dataset_name)==0) {
        glog << warning << "dataset '" << _dataset_name << "' do not exist for " << fnucleus << do_endl;
        print("dataset");
        return false;
    }
    fdatasetid = fmap_of_dataset_name_id[_dataset_name];
    return select_dataset(fdatasetid);
}
 
bool tklevel_scheme::select_dataset(int _dataset_id)
{
    if(fmap_of_dataset.count(_dataset_id)==0) {
        glog << warning << "dataset with id '" << _dataset_id << "' do not exist for " << fnucleus << do_endl;
        print("dataset");
        return false;
    }
 
    fdatasetid = _dataset_id;
    fmap_of_dataset[fdatasetid]->load_dataset();
    return true;
}
 
shared_ptr<tklevel> tklevel_scheme::get_level(const tkstring &_name, bool _exact)
{
    // get a specific level [spin][parity][rank] (ex: 0+1 for the first 0+ state)
    string spin;
    int rank;
 
    const char* sep = "";
    if(_name.contains("+"))      {sep = "+";}
    else if(_name.contains("-")) {sep = "-";}
 
    auto array = _name.tokenize(sep);
    if(array.size()!=2) {
        glog << error << " <nuclear_level_scheme::get_level_energy()>  invalid argument '" << _name << "'. Syntax: '[spin][parity][rank]'." << do_endl;
        return nullptr;
    }
 
    spin = array.at(0);
    rank = array.at(1).atoi();
 
    int idx = 0;
    for(const auto &lvlp : get_levels([](auto lvl) {return (lvl->get_spin_parity()!=nullptr);}))  {
        auto jpistr = lvlp->get_spin_parity()->get_jpi_str();
        bool ok=false;
        if(!_exact) {
            if(jpistr.equal_to(tkstring::form("%s%s",spin.data(),sep)) ||
                    jpistr.equal_to(tkstring::form("(%s%s)",spin.data(),sep)) ||
                    jpistr.equal_to(tkstring::form("(%s)%s",spin.data(),sep)) ||
                    jpistr.equal_to(tkstring::form("%s(%s)",spin.data(),sep))) {
                ok = true;
            }
        }
        else if(jpistr.equal_to(tkstring::form("%s%s",spin.data(),sep))) {
            ok = true;
        }
        if(ok) {
            idx++;
            if(idx==rank) return lvlp;
        }
    }
    return nullptr;
}
 
shared_ptr<tklevel> tklevel_scheme::get_level(double _energy, tkstring _offset)
{
    shared_ptr<tklevel> clostest_lev;
    double best_ediff=1e6;
    for(const auto &lev: get_levels()) {
        if(!_offset.is_empty() && !_offset.contains("*")) {
            if(_offset != lev->get_offset_bandhead()) continue;
        }
        if(_offset.is_empty() && !lev->get_offset_bandhead().is_empty()) continue;
        if(_offset.contains("dec") && lev->get_decays_down().size()==0) continue;
        double ediff = abs(lev->get_energy(tkunit_manager::units_keys::keV,true)-_energy);
        if(ediff<best_ediff) {
            best_ediff=ediff;
            clostest_lev = lev;
        }
    }
    return clostest_lev;
}
 
void tklevel_scheme::init()
{
    gdatabase->begin("DISTINCT dataset.dataset_id, dataset.dataset_name","level INNER JOIN dataset on level.dataset_id=dataset.dataset_id INNER JOIN isotope on level.isotope_id=isotope.isotope_id INNER JOIN element on isotope.element_id=element.element_id",tkstring::form("element.charge=%d AND isotope.mass=%d",fZ,fA), tkstring::form("read_dataset_%s",fnucleus.data()));
 
    while(gdatabase->next(tkstring::form("read_dataset_%s",fnucleus.data()))) {
        tkdb_table& datasets = (*gdatabase)["DATASET"];
        fmap_of_dataset_name_id[datasets["dataset_name"].get_value()] = datasets["dataset_id"].get_value().atoi();
        fmap_of_dataset[datasets["dataset_id"].get_value().atoi()] = make_shared<tkdataset>(tkdataset(fnucleus,fZ,fA,datasets["dataset_name"].get_value(),datasets["dataset_id"].get_value().atoi()));
    }
    tkstring default_dataset = fnucleus.data();
    default_dataset += " : ADOPTED LEVELS, GAMMAS";
    gdebug << "looking for dataset : " << default_dataset << do_endl;
    bool found = (fmap_of_dataset_name_id.count(default_dataset)>0);
    gdebug << (found?"":"not ") << "found " << (found?fmap_of_dataset_name_id[default_dataset.data()]:-1) << do_endl;
 
    if(found) select_dataset(fmap_of_dataset_name_id[default_dataset.data()]);
    else{
        for(const auto &dts : fmap_of_dataset) {
            glog << warning << tkstring::form("%5s",fnucleus.data()) << " -- no 'ADOPTED LEVELS' dataset : '" << dts.second->get_name() << "' used instead."  << do_endl;
            select_dataset(dts.first);
            break;
        }
    }
    if(fmap_of_dataset.empty()) {
        tkstring datasetname = "EMPTY";
        fdatasetid = 0;
        fmap_of_dataset_name_id[datasetname] = fdatasetid;
        fmap_of_dataset[fdatasetid] =  make_shared<tkdataset>(tkdataset(fnucleus,fZ,fA,datasetname,fdatasetid));
    }
 
    floaded = true;
}
 
void tkdataset::check_yrast(shared_ptr<tklevel> _lvl)
{
    // check Yrast
    if(_lvl->get_spin_parity()->get_spin().is_known()) {
        int TwoJ = (int)(_lvl->get_spin_parity()->get_spin().get_value()*2+0.5);
        if(!fyrastmap_exact.count(TwoJ)) {
            fyrastmap_exact[TwoJ] = _lvl;
            _lvl->fisYrast_exact = true;
            //                cout<<"exact yrast: 2J:"<<TwoJ<<" "<<lvl->get_spin_parity_str()<<" "<<lvl->get_energy()<<endl;
        }
        if(!fyrastmap_uncertain.count(TwoJ)) {
            fyrastmap_uncertain[TwoJ] = _lvl;
            _lvl->fisYrast_uncertain = true;
            //                cout<<"uncertain yrast: 2J:"<<TwoJ<<" "<<lvl->get_spin_parity_str()<<" "<<lvl->get_energy()<<endl;
        }
    }
    else if(!_lvl->get_spin_parity_str().is_empty() && !_lvl->get_spin_parity_str().contains(",")) {
        tkstring tmp = _lvl->get_spin_parity_str().copy().remove_all("(").remove_all(")").remove_all("+").remove_all("-");
        int TwoJ;
        if(tmp.contains("/2")) TwoJ = tmp.remove_all("/2").atoi();
        else TwoJ = tmp.atoi()*2;
        if(!fyrastmap_uncertain.count(TwoJ)) {
            fyrastmap_uncertain[TwoJ] = _lvl;
            _lvl->fisYrast_uncertain = true;
            //                cout<<"uncertain yrast: 2J:"<<TwoJ<<" "<<lvl->get_spin_parity_str()<<" "<<lvl->get_energy()<<endl;
        }
    }
}
 
void tkdataset::load_dataset()
{
    if(floaded) return;
 
#ifdef HAS_DEBUG
    auto t_0 = std::chrono::high_resolution_clock::now();
    gdebug  << "loading levels for dataset '" << fname << "' of " << fnucleus << do_endl;
#endif
 
    gdatabase->begin("DISTINCT level.*","level",tkstring::form("dataset_id=%d",fid));
 
    while(gdatabase->next()) {
        tkdb_table *levels = &(*gdatabase)["LEVEL"];
        tkdb_column *col = &(*levels)["level_id"];
        int id = col->get_value().atoi();
 
        gmanager->set_max_level_id(id);
 
        tkdb_table::measure_data_struct energy_struc;
        levels->read_measure(energy_struc, "level_energy");
 
        shared_ptr<tklevel> lvl = make_shared<tklevel>(id,energy_struc);
        lvl->fbelongs_to_nucleus = fnucleus;
 
        double spin = -1.;
        int parity = -1.;
        tkstring spin_parity_str = "";
 
        col = &(*levels)["level_spin_parity"];
        if(!col->is_empty()) spin_parity_str = col->get_value();
        col = &(*levels)["level_spin"];
        if(!col->is_empty()) spin = col->get_value().atof();
        col = &(*levels)["level_parity"];
        if(!col->is_empty()) parity = col->get_value().atoi();
 
        if(!spin_parity_str.is_empty()) lvl->set_jpi(spin,parity,spin_parity_str);
 
        tkdb_table::measure_data_struct lifetime_struc;
        levels->read_measure(lifetime_struc, "level_lifetime");
        if(lifetime_struc.filled) lvl->set_lifetime(lifetime_struc);
 
        col = &(*levels)["level_comment"];
        if(!col->is_empty()) lvl->set_comment(col->get_value());
        col = &(*levels)["level_uncertain"];
        if(!col->is_empty()) lvl->set_uncertain_level(col->get_value());
 
        flevels.push_back(lvl);
 
        check_yrast(lvl);
    }
    for(const auto &lvl: flevels) {
        fmapoflevels[lvl->get_id()] = lvl;
    }
 
#ifdef HAS_DEBUG
    auto t_1 = std::chrono::high_resolution_clock::now();
    gdebug  << "loading decays for dataset '" << fname << "' of " << fnucleus << do_endl;
#endif
 
    gdatabase->begin("DISTINCT decay.*","decay INNER JOIN level on decay.level_from_id=level.level_id ",tkstring::form("dataset_id=%d",fid));
 
    while(gdatabase->next()) {
        int dectype=-1;
        int level_from = -1.;
        int level_to = -1.;
 
        tkdb_table *decay = &(*gdatabase)["DECAY"];
        tkdb_column *col = &(*decay)["decay_id"];
        int id = col->get_value().atoi();
 
        gmanager->set_max_decay_id(id);
 
        col = &(*decay)["decay_type"];
        if(!col->is_empty()) dectype = col->get_value().atoi();
        col = &(*decay)["level_from_id"];
        if(!col->is_empty()) level_from = col->get_value().atoi();
        col = &(*decay)["level_to_id"];
        if(!col->is_empty()) level_to = col->get_value().atoi();
 
        shared_ptr<tkdecay> dec = nullptr;
 
        tkdb_table::measure_data_struct energy_struc;
        decay->read_measure(energy_struc, "decay_energy");
 
        if(dectype == tkdecay::kgamma) {
            dec = make_shared<tkgammadecay>(id, level_from, level_to, energy_struc);
            shared_ptr<tkgammadecay> gamma = dynamic_pointer_cast<tkgammadecay>(dec);
 
            tkdb_table::measure_data_struct relative_intensity_struc;
            decay->read_measure(relative_intensity_struc, "decay_relative_intensity");
            if(relative_intensity_struc.filled) gamma->set_relative_intensity(relative_intensity_struc);
 
            col = &(*decay)["decay_multipolarity"];
            if(!col->is_empty()) gamma->set_multipolarity(col->get_value());
 
            tkdb_table::measure_data_struct mixing_ratio_struc;
            decay->read_measure(mixing_ratio_struc, "decay_mixing_ratio");
            if(mixing_ratio_struc.filled) gamma->set_mixing_ratio(mixing_ratio_struc);
 
            tkdb_table::measure_data_struct conv_coeff_struc;
            decay->read_measure(conv_coeff_struc, "decay_conv_coeff");
            if(conv_coeff_struc.filled) gamma->set_conv_coeff(conv_coeff_struc);
 
            // transition probabilities
            tkdb_table::measure_data_struct transprob_struc;
            decay->read_measure(transprob_struc, "decay_BE");
            if(transprob_struc.filled) gamma->set_transition_probability(transprob_struc);
            transprob_struc.clear();
            decay->read_measure(transprob_struc, "decay_BEW");
            if(transprob_struc.filled) gamma->set_transition_probability(transprob_struc);
            transprob_struc.clear();
            decay->read_measure(transprob_struc, "decay_BM");
            if(transprob_struc.filled) gamma->set_transition_probability(transprob_struc);
            transprob_struc.clear();
            decay->read_measure(transprob_struc, "decay_BMW");
            if(transprob_struc.filled) gamma->set_transition_probability(transprob_struc);
        }
        else continue;
 
        col = &(*decay)["decay_comment"];
        if(!col->is_empty()) dec->set_comment(col->get_value());
        col = &(*decay)["decay_uncertain"];
        if(!col->is_empty()) dec->set_uncertain_decay(col->get_value());
 
        fdecays.push_back(dec);
        fdecays.back()->set_levels(fmapoflevels.at(level_from),fmapoflevels.at(level_to));
    }
    for(auto &dec: fdecays) {
        fmapoflevels.at(dec->get_level_from_id())->add_decay_down(dec);
        fmapoflevels.at(dec->get_level_to_id())->add_decay_up(dec);
    }
 
 
#ifdef HAS_DEBUG
    auto t_2 = std::chrono::high_resolution_clock::now();
 
    auto dt01 = std::chrono::duration<double, std::milli>(t_1-t_0).count() * 0.001; // open db time
    auto dt12 = std::chrono::duration<double, std::milli>(t_2-t_1).count() * 0.001; // first read loop with views
 
    gdebug << std::fixed << std::setprecision(4) << "time level: "
           << dt01
           << " s" << do_endl;
 
    gdebug << std::fixed << std::setprecision(4) << "time decays: "
           <<dt12
          << " s" << do_endl;
#endif
 
    floaded = true;
}
 
void tkdataset::print(const tkstring &_data, const tkstring &_option)
{
    if(_data.contains("level") && _data.contains("decay")) {
        glog << info << "dataset '" << fname << "' contains " << flevels.size() << " levels" << do_endl;
        glog << info << "dataset '" << fname << "' contains " << fdecays.size() << " decays" << do_endl;
        for(auto &lev : flevels) {
            if(_option.contains("yrastt") && !lev->is_yrast(true)) continue;
            else if(_option.contains("yrast") && !_option.contains("yrastt") && !lev->is_yrast()) continue;
            lev->print(_option);
            for(auto &dec : lev->get_decays_down()) {
                cout<<setw(13)<<""<<"-> ";dec->print("quiet;levto");
            }
        }
    }
    else if(_data.contains("level")) {
        if(_option.contains("tab")){
            std::cout << std::setw(15) << left << "symbol" << std::setw(15) << "levelid"
                      << std::setw(15) << left << "energy (keV)"
                      << std::setw(15) << left << "energy error"
                      << std::setw(15) << left << "spin parity"
                      << std::setw(15) << left << "t1/2"
                      << std::setw(15) << left << "t1/2 error"
                      << std::setw(15) << left << "t1/2 unit" << std::endl;}
        else  glog << info << "dataset '" << fname << "' contains " << flevels.size() << " levels:" << do_endl;
        for(auto &lev : flevels) {
            if(_option.contains("yrastt") && !lev->is_yrast(true)) continue;
            else if(_option.contains("yrast") && !_option.contains("yrastt") && !lev->is_yrast()) continue;
            if(_option.contains("tab")) std::cout << std::setw(15)<< left << fnucleus;
            lev->print(_option);
        }
    }
    else if(_data.contains("decay")) {
        if(_option.contains("tab")) std::cout << std::setw(15) << left << "symbol" << std::setw(15) << left << "decayid" << std::setw(15) << left << "energy (keV)" << std::setw(15) << left << "energy error (keV)"   << std::endl;
        else glog << info << "dataset '" << fname << "' contains " << fdecays.size() << " decays:" << do_endl;
        for(auto &dec : fdecays) {
            if(_option.contains("tab")) std::cout << std::setw(15) << left << fnucleus;
            dec->print(_option);
        }
    }
}
 
shared_ptr<tklevel> tkdataset::add_level(double _ener, double _unc, tkstring _unit, tkstring _jpistr)
{
    tkspin spin;
    tkparity parity;
    double spin_val = -1.;
    int parity_val = -1.;
 
    tkdb_table::measure_data_struct lvl_properties{_ener,_unc,0,0,_unit};
    shared_ptr<tklevel> lvl = make_shared<tklevel>(gmanager->get_new_level_id(),lvl_properties);
 
    spin.set(_jpistr);parity.set(_jpistr);
    if(spin.is_known()) spin_val = spin.get_value();
    if(parity.is_known()) parity_val = (parity.is_parity(tkparity::eparity::kParityMinus));
 
    lvl->set_jpi(spin_val,parity_val,_jpistr);
 
    // needs to generate a unique id for the level
    lvl->fbelongs_to_nucleus = fnucleus;
 
    // find the good position in the vector to be sorted by energy
    auto it = std::lower_bound(flevels.begin(), flevels.end(), lvl,
                               [](const std::shared_ptr<tklevel>& a, const std::shared_ptr<tklevel>& b) {
                                   return a->get_energy() < b->get_energy();
                               });
    flevels.insert(it, lvl);
 
    check_yrast(lvl);
    fmapoflevels[lvl->get_id()] = lvl;
 
    return lvl;
}
 
shared_ptr<tkgammadecay> tkdataset::add_gamma_decay(shared_ptr<tklevel> _lvlfrom, shared_ptr<tklevel> _lvlto, double _ener, double _unc)
{
    if(_ener==0.) _ener = _lvlfrom->get_energy() - _lvlto->get_energy();
 
    tkdb_table::measure_data_struct energy_struc{_ener,_unc,0.,0.,"keV"};
 
    shared_ptr<tkgammadecay> decay = make_shared<tkgammadecay>(gmanager->get_new_decay_id(), _lvlfrom->get_id(), _lvlto->get_id(), energy_struc);
    decay->set_levels(_lvlfrom,_lvlto);
    fdecays.push_back(decay);
    _lvlfrom->add_decay_down(decay);
    _lvlto->add_decay_up(decay);
 
    return decay;
}
 
 
#ifdef HAS_ROOT
ClassImp(tklevel_scheme);
#endif