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]                          *
 *                                                                              *
 *   Licensed under the MIT License <http://opensource.org/licenses/MIT>.       *
 *   SPDX-License-Identifier: MIT                                               *
 ********************************************************************************/
 
#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;
         tkstring date = dataset.fdate.copy().remove_all(" ");
         if (date.length() == 6) date = date.substr(0, 4) + "-" + date.substr(4, 2);
         int datasetIdx = add_dataset(dataset, nuc.data(), get_dataset_comment(&dataset), date, "ENSDF");
         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;
         tkstring date = dataset.fdate.copy().remove_all(" ");
         if (date.length() == 6) date = date.substr(0, 4) + "-" + date.substr(4, 2);
         int datasetIdx = add_dataset(dataset, nuc.data(), get_dataset_comment(&dataset), date, "XUNDL");
         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 dataset_isotope_index ON dataset(isotope_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 &_nucleus, const tkstring &_dataset_comment, const tkstring &_dataset_date, const tkstring &_dataset_source)
{
   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", _nucleus.data()); //.prepend(tkstring::form("%s : ",_nucleus.data()));
      _dataset.fis_adopted = true;
   } else
      _dataset.fis_adopted = false;
 
   (*fDataBase)["DATASET"]["dataset_name"].set_value(newname);
   (*fDataBase)["DATASET"]["isotope_id"].set_value(fIsotopeIndex);
   if (!_dataset_comment.is_empty()) (*fDataBase)["DATASET"]["dataset_comment"].set_value(_dataset_comment);
   if (!_dataset_date.is_empty()) (*fDataBase)["DATASET"]["dataset_date"].set_value(_dataset_date);
   if (!_dataset_source.is_empty()) (*fDataBase)["DATASET"]["dataset_source"].set_value(_dataset_source);
   (*fDataBase)["DATASET"]["dataset_id"].set_value(fDatasetIndex++);
   (*fDataBase)["DATASET"].push_row();
 
   return (fDatasetIndex - 1);
}
 
tkstring tkensdf_builder::get_dataset_comment(tkensdf_ident_rec *_data_set)
{
   tkstring record;
   bool is_record = fensdf_reader.first_record(_data_set, record);
 
   tkensdf_record the_record;
   tkensdf_record dataset_comment;
 
   while (is_record) {
      if (!the_record.set_record(record)) {
         is_record = fensdf_reader.next_record(_data_set, record);
         continue;
      }
 
      if (the_record.get_record_type() == tkensdf_record::klevel) break;
      if (the_record.is_comment()) dataset_comment.add_comment_record(record, the_record.is_continuation_record());
 
      is_record = fensdf_reader.next_record(_data_set, record);
   }
 
   return dataset_comment.get_comment_record();
}
 
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
               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