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]                          *
 *                                                                              *
 *   Licensed under the MIT License <http://opensource.org/licenses/MIT>.       *
 *   SPDX-License-Identifier: MIT                                               *
 ********************************************************************************/
 
#ifdef HAS_DEBUG
#include <algorithm>
#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, bool _load_default_dataset) : fnucleus(_nuc), fA(_aa), fZ(_zz)
{
   if (fA < 0 || fZ < 0) return;
   init(_load_default_dataset);
}

std::vector<shared_ptr<tklevel>> tklevel_scheme::get_levels(const std::function<bool(shared_ptr<tklevel>)> &_selection)
{
   vector<shared_ptr<tklevel>> res;
   const auto &levels = get_levels();
   std::copy_if(levels.begin(), levels.end(), std::back_inserter(res),
                [&_selection](const auto &lvl) { return _selection(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;
   const auto &decays = get_decays();
   std::copy_if(decays.begin(), decays.end(), std::back_inserter(res),
                [&_selection](const auto &dec) { return _selection(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")) {
      auto print_group = [&](const tkstring &source_label) {
         bool header_printed = false;
         for (auto &dts : fmap_of_dataset) {
            if (dts.second->get_source() != source_label) continue;
            if (!header_printed) {
               glog << info << "--- " << source_label << " ---" << do_endl;
               header_printed = true;
            }
            glog << info << dts.second->get_name() << " (" << dts.first << ")";
            if (dts.second->has_date()) glog << " [" << dts.second->get_date() << "]";
            glog << do_endl;
            if (_option.contains("com") && dts.second->has_comment()) {
               const size_t comment_prefix_size = 13;
               const size_t content_col = comment_prefix_size + 2;
               glog << comment << "  " << wrap_text(dts.second->get_comment(), content_col, content_col) << do_endl;
            }
         }
      };
      print_group("ENSDF");
      print_group("XUNDL");
      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)
{
   const bool require_decay = _offset.contains("dec");
   _offset.remove_all("dec").strip_all_extra_white_space();
   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 (require_decay && lev->get_decays_down().empty()) 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;
}
 
shared_ptr<tklevel> tklevel_scheme::get_ground_state()
{
   const auto name = fnucleus + " : ADOPTED LEVELS, GAMMAS";
   const auto id = fmap_of_dataset_name_id.find(name);
   if (id == fmap_of_dataset_name_id.end()) return nullptr;
   const auto dataset = fmap_of_dataset.at(id->second);
   dataset->load_dataset();
   for (const auto &level : dataset->get_levels()) {
      if (!level->is_energy_offset() && level->get_energy() == 0.) return level;
   }
   return nullptr;
}
 
void tklevel_scheme::init(bool _load_default_dataset)
{
   gdatabase->begin("dataset.dataset_id, dataset.dataset_name, dataset.dataset_comment, dataset.dataset_date, dataset.dataset_source",
                    "dataset INNER JOIN isotope on dataset.isotope_id=isotope.isotope_id INNER JOIN element on isotope.element_id=element.element_id",
                    tkstring::form("element.charge=%d AND isotope.mass=%d AND EXISTS (SELECT 1 FROM level WHERE level.dataset_id=dataset.dataset_id)", fZ, fA),
                    tkstring::form("read_dataset_%s", fnucleus.data()));
 
   while (gdatabase->next(tkstring::form("read_dataset_%s", fnucleus.data()))) {
      tkdb_table &datasets = (*gdatabase)["DATASET"];
      int ds_id = datasets["dataset_id"].get_value().atoi();
      fmap_of_dataset_name_id[datasets["dataset_name"].get_value()] = ds_id;
      auto ds = make_shared<tkdataset>(tkdataset(fnucleus, fZ, fA, datasets["dataset_name"].get_value(), ds_id));
      if (!datasets["dataset_comment"].is_empty()) ds->fdataset_comment = datasets["dataset_comment"].get_value();
      if (!datasets["dataset_date"].is_empty()) ds->fdataset_date = datasets["dataset_date"].get_value();
      if (!datasets["dataset_source"].is_empty()) ds->fdataset_source = datasets["dataset_source"].get_value();
      fmap_of_dataset[ds_id] = ds;
   }
   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) {
      fdatasetid = fmap_of_dataset_name_id[default_dataset.data()];
      if (_load_default_dataset) select_dataset(fdatasetid);
   } 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;
         fdatasetid = dts.first;
         if (_load_default_dataset) select_dataset(fdatasetid);
         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("level.*", "level", tkstring::form("dataset_id=%d", fid));
 
   tkdb_table *levels = &(*gdatabase)["LEVEL"];
   const auto energy_columns = levels->get_measure_columns("level_energy");
   const auto lifetime_columns = levels->get_measure_columns("level_lifetime");
   const auto magnetic_dipole_columns =
      levels->get_measure_columns("level_magnetic_dipole");
   const auto electric_quadrupole_columns =
      levels->get_measure_columns("level_electric_quadrupole");
 
   while (gdatabase->next()) {
      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, energy_columns);
 
      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, lifetime_columns);
      if (lifetime_struc.filled) lvl->set_lifetime(lifetime_struc);
 
      tkdb_table::measure_data_struct magnetic_dipole_struc;
      levels->read_measure(magnetic_dipole_struc, magnetic_dipole_columns);
      if (magnetic_dipole_struc.filled) lvl->set_magnetic_dipole(magnetic_dipole_struc);
 
      tkdb_table::measure_data_struct electric_quadrupole_struc;
      levels->read_measure(electric_quadrupole_struc, electric_quadrupole_columns);
      if (electric_quadrupole_struc.filled) lvl->set_electric_quadrupole(electric_quadrupole_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("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;
      bool first_level_printed = false;
      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 (first_level_printed) cout << std::endl;
         first_level_printed = true;
         lev->print(_option);
         for (auto &dec : lev->get_decays_down()) {
            tkstring decay_option = "quiet;levto";
            if (_option.contains("com")) decay_option += ";com";
            cout << setw(13) << "" << "-> ";
            dec->print(decay_option);
         }
      }
   } 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