tkensdf_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 "tkensdf_builder.h"
 
#include <iostream>
#include <unistd.h>
#include <cmath>
 
#include "tklog.h"
#include "tkensdf_level_rec.h"
#include "tkensdf_gamma_rec.h"
 
namespace tkn {
}
 
using namespace tkn;
using namespace std;
 
tkensdf_builder::tkensdf_builder(tkdatabase *_database, const tkstring &_input_folder):
    fDataBase(_database),
    fensdf_reader(_input_folder)
{
}
 
tkensdf_builder::~tkensdf_builder() = default;
 
int tkensdf_builder::fill_database()
{
    //--------------------------------------------------------
    // read files and fill database
    //--------------------------------------------------------
 
    fDataBase->begin("charge,mass,symbol,isotope_id","isotope INNER JOIN element ON isotope.element_id=element.element_id","charge>0");
 
    tkstring message;
 
    if(flevel_builder) {
        message = "Building '" + flevel_builder->get_table().get_name() + "' table";
    }
 
    int oldZ=0;
    while(fDataBase->next())
    {
        tkstring nuc = tkstring::form("%s%s",(*fDataBase)["ISOTOPE"]["mass"].get_value().data(),(*fDataBase)["ELEMENT"]["symbol"].get_value().data());
        tkstring zz = (*fDataBase)["ELEMENT"]["charge"].get_value();
        tkstring aa = (*fDataBase)["ISOTOPE"]["mass"].get_value();
        fIsotopeIndex = ((tkstring)(*fDataBase)["ISOTOPE"]["isotope_id"].get_value()).atoi();
 
        int curZ = zz.atoi();
        if(curZ<1) continue;
 
        // open the ensdf file
        bool has_levels = fensdf_reader.open_nuc(nuc,kensdf);
        //        fensdf_reader.print_datasets();
 
        // loop over the datasets
        for(auto &dataset: *fensdf_reader.get_datasets()) {
            if(dataset.fdsid.begins_with("COMMENTS")) continue;
            int datasetIdx = add_dataset(dataset,nuc.data());
            read_dataset(&dataset, datasetIdx);
        }
 
        // open the xundl file
        bool has_levelsxundl = fensdf_reader.open_nuc(nuc,kxundl);
        // loop over the datasets
        for(auto &dataset: *fensdf_reader.get_datasets()) {
            if(dataset.fdsid.begins_with("COMMENTS")) continue;
            int datasetIdx = add_dataset(dataset,nuc.data());
            read_dataset(&dataset, datasetIdx);
        }
 
        if(curZ!=oldZ) {
            glog.progress_bar(117,curZ,message.data()); oldZ=curZ;
        }
 
        if(has_levels||has_levelsxundl){
            fDataBase->exec_sql(tkstring::form("CREATE VIEW lvl%s as select * from 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 where element.charge=%s AND isotope.mass=%s",nuc.data(),zz.data(),aa.data()));
            if(fdecay_builder) fDataBase->exec_sql(tkstring::form("CREATE VIEW dec%s as select * from decay INNER JOIN level on decay.level_from_id=level.level_id 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 where element.charge=%s AND isotope.mass=%s",nuc.data(),zz.data(),aa.data()));
        }
    }
 
    fDataBase->exec_sql("CREATE INDEX isotope_index ON level(isotope_id)");
    fDataBase->exec_sql("CREATE INDEX dataset_index ON level(dataset_id)");
 
    fDataBase->exec_sql("CREATE INDEX levelfrom_index ON decay(level_from_id)");
    fDataBase->exec_sql("CREATE INDEX levelto_index ON decay(level_to_id)");
 
    fensdf_reader.print_record_counters();
 
    return 0;
}
 
int tkensdf_builder::add_dataset(tkensdf_ident_rec &_dataset, const tkstring &_comment)
{
    tkstring newname = _dataset.fdsid.copy().remove_all_extra_white_space();
    newname.replace_all("'","`");
 
    // use begins_with to handle at the same time:
    // - ADOPTED LEVELS  (no known gammas)
    // - ADOPTED LEVELS, GAMMAS (known gammas)
    if(_dataset.fdsid.begins_with("ADOPTED LEVELS") && !_dataset.fdsid.contains("TENTATIVE")) {
        newname = tkstring::form("%s : ADOPTED LEVELS, GAMMAS",_comment.data());//.prepend(tkstring::form("%s : ",_comment.data()));
        _dataset.fis_adopted = true;
    }
    else _dataset.fis_adopted = false;
 
    (*fDataBase)["DATASET"]["dataset_name"].set_value(newname);
    (*fDataBase)["DATASET"]["dataset_comment"].set_value(_comment);
    (*fDataBase)["DATASET"]["dataset_id"].set_value(fDatasetIndex++);
    (*fDataBase)["DATASET"].push_row();
 
    return (fDatasetIndex-1);
}
 
bool tkensdf_builder::read_dataset(tkensdf_ident_rec *_data_set, int _dataset_idx)
{
    glog.set_class("db_ensdf_builder");
    glog.set_method(tkstring::form("read_levels(%s,%d)", _data_set->fdsid.data(),_dataset_idx));
 
    tkensdf_record the_record;
    tkensdf_level_rec the_level_record;
    tkensdf_gamma_rec the_gamma_record;
 
    tkstring record;
    bool is_record = fensdf_reader.first_record(_data_set,record);
 
    std::vector<tkensdf_level_rec> flist_of_levels;
 
    char * sErrMsg = nullptr;
    sqlite3_exec(fDataBase->get_sql_db(), "BEGIN TRANSACTION", nullptr, nullptr, &sErrMsg);
 
    // loop over the dataset and fill temporary collections.
    while( is_record ) {
 
        if(!the_record.set_record(record)) {
            is_record = fensdf_reader.next_record(_data_set,record);
            glog << error_v << "error in : " << _data_set->fnuclide.data()<<": "<<_data_set->fdsid.data()<<" line: "<<_data_set->fcurrentPosition<< do_endl;
            continue;
        }
 
        // get the record type
        tkensdf_record::record_type the_type = the_record.get_record_type();
 
        // check and analyse the global comments records before the first level
        if(flist_of_levels.empty() && the_record.is_comment()) {
            if(record.contains("T$LABEL",tkstring::ECaseCompare::kIgnoreCase) || record.contains("T         LABEL",tkstring::ECaseCompare::kIgnoreCase)) {
                // skip cases where time info is replaced by cross section data
                if(record.contains("mb/sr",tkstring::ECaseCompare::kIgnoreCase)) the_level_record.fglobal_time_unit = "skip";
                else if(record.contains("d|s/d|W",tkstring::ECaseCompare::kIgnoreCase)) the_level_record.fglobal_time_unit = "skip";
                else if(record.contains("DS/DW",tkstring::ECaseCompare::kIgnoreCase)) the_level_record.fglobal_time_unit = "skip";
                else if(record.contains("d(SIGMA)/d(OMEGA)",tkstring::ECaseCompare::kIgnoreCase)) the_level_record.fglobal_time_unit = "skip";
                // skip cases where time info is replaced by Pygmy Dipole Resonance (PDR) strength (SIS)
                else if(record.contains("S{-IS}",tkstring::ECaseCompare::kIgnoreCase)) the_level_record.fglobal_time_unit = "skip";
                // no real idea of what it is..
                else if(record.contains("G{-|g0",tkstring::ECaseCompare::kIgnoreCase)) the_level_record.fglobal_time_unit = "skip";
                // here define the global witdth units
                else if(record.contains("mev",tkstring::ECaseCompare::kIgnoreCase)) the_level_record.fglobal_time_unit = "MEV";
                else if(record.contains("kev",tkstring::ECaseCompare::kIgnoreCase)) the_level_record.fglobal_time_unit = "KEV";
                else if(record.contains("ev",tkstring::ECaseCompare::kIgnoreCase)) the_level_record.fglobal_time_unit = "EV";
            }
        }
 
        if(fdecay_builder && the_type == tkensdf_record::kgamma && the_level_record.get_energy().value>0.) {
 
            the_gamma_record.clear();
            the_gamma_record.set_record(record);
 
            tkensdf_level_rec *level_to = nullptr;
 
            // check if the following records are comments or continuation records on the level
            is_record = fensdf_reader.next_record(_data_set,record);
            the_record.set_record(record);
            while( is_record && (the_record.is_continuation_record() || the_record.is_comment()) ) {
                if(the_record.is_comment()) the_gamma_record.add_comment_record(record,the_record.is_continuation_record());
                else if(the_record.is_continuation_record()) the_gamma_record.add_continuation_record(record);
                is_record = fensdf_reader.next_record(_data_set,record);
                the_record.set_record(record);
            }
 
            the_gamma_record.analyse_record();
 
            // in specific cases that are not treated, the gamma energy is set to 0
            if(the_gamma_record.get_energy().value==0.) continue;
 
            // if the level-to is manually specified:
            if(the_gamma_record.is_final_level_set()) {
 
                // if this is an unplaces gamma-ray, we skip it
                if(the_gamma_record.get_final_level_energy()<0)
                    continue;
 
                // else we get the exact energy in the list
                for(auto it = flist_of_levels.rbegin()+1 ; it != flist_of_levels.rend() ; ++it) {
                    auto level =  &(*it);
                    if(the_gamma_record.get_final_level_energy() == level->get_energy().value) {
                        level_to = level;
                    }
                }
                if(level_to == nullptr) the_gamma_record.set_final_level_set(false);
            }
            if(level_to == nullptr) {
                // if the level-to is not specified, or has not been found, we take the closest one
                double lev_to_approx_energy = the_level_record.get_energy().value - the_gamma_record.get_energy().value;
 
                // looking for the id of the level to
                double last_diff = numeric_limits<double>::max();
                for(auto it = flist_of_levels.rbegin()+1 ; it != flist_of_levels.rend() ; ++it) {
                    auto level =  &(*it);
                    if(the_level_record.get_energy_offset() != level->get_energy_offset()) continue;
                    double diff = abs(lev_to_approx_energy-level->get_energy().value);
                    if(diff>last_diff) break;
                    level_to = level;
                    last_diff = diff;
                }
 
                if(level_to == nullptr) {
                    gdebug <<"Level_to not found !"<< do_endl;
#ifdef HAS_DEBUG
                    cout<<_data_set->fnuclide.data()<<": "<<_data_set->fdsid.data()<<" line: "<<_data_set->fcurrentPosition<<endl;
                    cout<< "Level from record: ";
                    if(the_level_record.is_energy_offset()) cout << the_level_record.get_energy_offset() << " + ";
                    cout << the_level_record.get_energy().value << " " << the_level_record.get_energy().err <<endl;
                    cout<<" record -> " << the_level_record.get_record() << endl;
                    cout<< "gamma record     : " << the_gamma_record.get_energy().value << " " << the_gamma_record.get_energy().err <<endl;
                    cout<<" record -> " << the_gamma_record.get_record() << endl;
                    if(the_gamma_record.is_final_level_set()) cout<<"gamma levels manually set !" << endl;
                    cout<< " Looking for a level at: ";
                    if(the_level_record.is_energy_offset()) cout << the_level_record.get_energy_offset() << " + ";
                    cout << lev_to_approx_energy<< " keV in:"<<endl;                    for(auto it = flist_of_levels.rbegin()+1 ; it != flist_of_levels.rend() ; ++it) {
                        if(it->get_energy_offset()!=the_level_record.get_energy_offset()) continue;
                        if(it->is_energy_offset()) cout << it->get_energy_offset() << " + ";
                        cout << it->get_energy().value << " keV" << endl;
                    }
#endif
                }
                if(!level_to) continue;
                // if no level-to found closest to 100keV above error bars, the gamma-ray is skipped
                bool too_high_diff = ((last_diff-level_to->get_energy().err-the_level_record.get_energy().err-the_gamma_record.get_energy().err) > 100);
                if(too_high_diff) {
                    gdebug <<"High energy difference !"<< do_endl;
#ifdef HAS_DEBUG
                    cout<<_data_set->fnuclide.data()<<": "<<_data_set->fdsid.data()<<" line: "<<_data_set->fcurrentPosition<<endl;
                    cout<< "Level from record: ";
                    if(the_level_record.is_energy_offset()) cout << the_level_record.get_energy_offset() << " + ";
                    cout << the_level_record.get_energy().value << " " << the_level_record.get_energy().err <<endl;
                    cout<<" record -> " << the_level_record.get_record() << endl;
                    cout<< "gamma record     : " << the_gamma_record.get_energy().value << " " << the_gamma_record.get_energy().err <<endl;
                    cout<<" record -> " << the_gamma_record.get_record() << endl;
                    if(the_gamma_record.is_final_level_set()) cout<<"gamma levels manually set !" << endl;
                    cout<< " Looking for a level at: ";
                    if(the_level_record.is_energy_offset()) cout << the_level_record.get_energy_offset() << " + ";
                    cout << lev_to_approx_energy<< " keV in:"<<endl;
                    for(auto it = flist_of_levels.rbegin()+1 ; it != flist_of_levels.rend() ; ++it) {
                        if(it->get_energy_offset()!=the_level_record.get_energy_offset()) continue;
                        if(it->is_energy_offset()) cout << it->get_energy_offset() << " + ";
                        cout << it->get_energy().value << " keV" << endl;
                    }
                    cout<< " Found a level at: " << level_to->get_energy().value<< " keV ==> Diff= "<<last_diff<< " keV" << endl;
#endif
                    continue;
                }
            }
 
            // fill level info in db
            fdecay_builder->fill_gamma(the_level_record.flevelid, level_to->flevelid, the_gamma_record);
        }
        else if(the_type == tkensdf_record::klevel) {
            the_level_record.clear();
            the_level_record.set_record(record);
 
            // check if the following records are comments or continuation records on the level
            is_record = fensdf_reader.next_record(_data_set,record);
            the_record.set_record(record);
            while( is_record && (the_record.is_continuation_record() || the_record.is_comment()) ) {
                if(the_record.is_comment()) the_level_record.add_comment_record(record,the_record.is_continuation_record());
                else if(the_record.is_continuation_record()) the_level_record.add_continuation_record(record);
                is_record = fensdf_reader.next_record(_data_set,record);
                the_record.set_record(record);
            }
 
            the_level_record.analyse_record();
 
            // fill level info in db
            if(flevel_builder && (the_level_record.get_energy().value>=0.)) {
 
                // if no ground state, add a dummy level a 0 energy
                if(flist_of_levels.size()==0 && (the_level_record.get_energy().value>0. || the_level_record.is_energy_offset())) {
                    tkensdf_level_rec dummy_rec;
                    dummy_rec.set_record("Dummy  L    0.0                                                                 ");
                    dummy_rec.analyse_record();
                    flevel_builder->fill_level(_dataset_idx, fIsotopeIndex, dummy_rec, _data_set->fis_adopted);
                    dummy_rec.flevelid = flevel_builder->get_level_id();
                    flist_of_levels.emplace_back(dummy_rec);
                }
 
                flevel_builder->fill_level(_dataset_idx, fIsotopeIndex, the_level_record, _data_set->fis_adopted);
                the_level_record.flevelid = flevel_builder->get_level_id();
                flist_of_levels.emplace_back(the_level_record);
 
                // check that the energy is not higher thant the previous one
                if(the_level_record.is_energy_offset()) {
                    std::vector<pair<int,tkensdf_level_rec*>> vec;
                    int idx=0;
                    for(auto it = flist_of_levels.begin() ; it != flist_of_levels.end() ; it++, idx++) {
                        if(it->get_energy_offset() == the_level_record.get_energy_offset()) {
                            vec.push_back({idx,(&(*it))});
                        }
                    }
                    if(vec.size()>1 && (vec.at(vec.size()-1).second->get_energy().value < vec.at(vec.size()-2).second->get_energy().value)) {
                        gdebug <<"Inverted levels in: "<<_data_set->fnuclide.data()<<": "<<_data_set->fdsid.data() << " [ "<< flist_of_levels.at(vec.at(vec.size()-1).first).get_energy_offset() << "+" << flist_of_levels.at(vec.at(vec.size()-1).first).get_energy().value << " ; " << flist_of_levels.at(vec.at(vec.size()-2).first).get_energy_offset() << "+" << flist_of_levels.at(vec.at(vec.size()-2).first).get_energy().value<<"]"<<do_endl;
                        swap(flist_of_levels[vec.at(vec.size()-1).first],flist_of_levels[vec.at(vec.size()-2).first]);
                        cin.get();
                    }
                }
                else if(flist_of_levels.size()>1 && (flist_of_levels.at(flist_of_levels.size()-1).get_energy().value < flist_of_levels.at(flist_of_levels.size()-2).get_energy().value)) {
                    gdebug <<"Inverted levels in: "<<_data_set->fnuclide.data()<<": "<<_data_set->fdsid.data() << " ["<< flist_of_levels.at(flist_of_levels.size()-1).get_energy().value << " ; " << flist_of_levels.at(flist_of_levels.size()-2).get_energy().value<<"]"<<do_endl;
                    swap(flist_of_levels[flist_of_levels.size()-1],flist_of_levels[flist_of_levels.size()-2]);
                }
            }
        }
        else {
            is_record = fensdf_reader.next_record(_data_set,record);
        }
    } // ok
    sqlite3_exec(fDataBase->get_sql_db(), "END TRANSACTION", nullptr, nullptr, &sErrMsg);
 
 
    glog.clear();
 
    return true;
}
 
#ifdef HAS_ROOT
ClassImp(tkensdf_builder);
#endif