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]                          *
 *                                                                              *
 *   Licensed under the MIT License <http://opensource.org/licenses/MIT>.       *
 *   SPDX-License-Identifier: MIT                                               *
 ********************************************************************************/
 
#include <algorithm>
#include <atomic>
#include <array>
#include <fstream>
#include <iterator>
 
#include "tkmanager.h"
#include "tklog.h"
#include "tkdatabase.h"
#include "json.hpp"
 
namespace tkn {
 
}
 
using namespace tkn;
using json = nlohmann::json;
 
tkmanager *tkmanager::the_data_manager()
{
   static tkmanager *instance = new tkmanager();
   return instance;
}
 
tkmanager::tkmanager()
{
   std::lock_guard<std::recursive_mutex> lock(get_mutex());
 
   preload_nuclei();
 
}
 
tkmanager::~tkmanager() = default;
 
bool tkmanager::is_level_scheme_loaded(const tkstring &_nuc)
{
   const auto scheme = fmap_of_level_scheme.find(_nuc);
   return scheme != fmap_of_level_scheme.end() &&
          std::atomic_load_explicit(&scheme->second, std::memory_order_acquire) != nullptr;
}
 
shared_ptr<tklevel_scheme> tkmanager::get_level_scheme(const tkn::tkstring &_nuc, int _zz, int _aa)
{
   const auto scheme = fmap_of_level_scheme.find(_nuc);
   const auto once = fmap_of_level_scheme_once.find(_nuc);
   if (scheme == fmap_of_level_scheme.end() || once == fmap_of_level_scheme_once.end())
      return make_shared<tklevel_scheme>(_nuc, _zz, _aa);
 
   std::call_once(*once->second, [&] {
      // tklevel_scheme construction still reads through the legacy global
      // gdatabase cursor.  Serialise only that construction phase; cached
      // read access below remains independent between threads.
      std::lock_guard<std::mutex> load_lock(flevel_scheme_load_mutex);
      auto loaded = make_shared<tklevel_scheme>(_nuc, _zz, _aa);
      std::atomic_store_explicit(&scheme->second, std::move(loaded), std::memory_order_release);
   });
   return std::atomic_load_explicit(&scheme->second, std::memory_order_acquire);
}
 
void tkmanager::preload_nuclei()
{
   // open() already reported the failure; do not query a null connection.
   if (!gdatabase->get_sql_db()) return;
 
   fvector_of_nuclei.reserve(4000);
 
   tkdb_table &element = (*gdatabase)["ELEMENT"];
   tkdb_table &iso_table = (*gdatabase)["ISOTOPE"];
 
   static const std::array<tkstring, 17> property_names = {
      "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 primary AME and IAEA tables
   };
   std::array<tkdb_table::measure_columns, property_names.size()> property_columns;
   for (std::size_t index = 0; index < property_names.size(); ++index)
      property_columns[index] = iso_table.get_measure_columns(property_names[index]);
 
   // try with isotope.*, element.*
   gdatabase->begin("isotope.*, element.*", "isotope INNER JOIN element on isotope.element_id=element.element_id");
 
   while (gdatabase->next()) {
      shared_ptr<tknucleus> anucleus = make_shared<tknucleus>();
 
      tkstring name = element["symbol"].get_value();
      anucleus->felement_symbol = name;
      anucleus->felement_name = element["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;
      anucleus->fis_known = true;
 
      // 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();
         // SQL NULL values are represented by an empty tkdb_column.  Do not
         // expose them as real properties with a numeric value of zero.
         if (value.is_empty()) continue;
         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
      for (std::size_t index = 0; index < property_names.size(); ++index) {
         const auto &property_name = property_names[index];
         tkdb_table::measure_data_struct data_struct;
         iso_table.read_measure(data_struct, property_columns[index]);
         if (data_struct.filled) {
            shared_ptr<tkmeasure> data;
            if (data_struct.err_low > 0. && data_struct.err_high > 0.)
               data = make_shared<tkmeasure>(data_struct.value, data_struct.unit,
                                             data_struct.err_low, data_struct.err_high);
            else
               data = make_shared<tkmeasure>(data_struct.value, data_struct.unit, data_struct.err);
            data->set_info(tkproperty::kKnown);
            if (!data_struct.info.is_empty()) data->set_info_tag(data_struct.info);
            data->set_type(property_name);
            if (property_name == "lifetime")
               data->set_type(gunits->get_type(data->get_unit_key()) == gunits->kEnergy ? "decay width" : "half-life");
 
            if (property_name == "lifetime") {
               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);
            }
         }
      }
 
      tkstring discovery_year = iso_table["isotope_year_discovered"].get_value();
      if (!discovery_year.is_empty())
         anucleus->add_property_str("isotope_year_discovered", discovery_year, "");
 
      tkstring decay_modes = iso_table["decay_modes"].get_value();
      if (!decay_modes.is_empty())
         anucleus->add_property_str("decay_modes", decay_modes, "");
 
      tkstring spin_parity_str = iso_table["spin_parity"].get_value();
      if (!spin_parity_str.is_empty()) {
         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);
      fmap_of_level_scheme.emplace(name, nullptr);
      fmap_of_level_scheme_once.emplace(name, make_shared<std::once_flag>());
 
      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);
      auto alpha_daughter = get_nucleus(nuc.second->get_z() - 2, nuc.second->get_a() - 4);
      if (daughter && daughter->get("mass_excess") && alpha_daughter && alpha_daughter->get("mass_excess")) {
         // Simplify before propagating errors: the alpha mass cancels and the
         // central mass occurs twice, not as two independent measurements.
         shared_ptr<tkmeasure> Q = make_shared<tkmeasure>(
            0.5 * (*daughter->get("mass_excess") + *alpha_daughter->get("mass_excess") - 2. * (*ME)));
         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")) {
         auto q_positron = make_shared<tkmeasure>(*ME - *daughter->get("mass_excess") - 2. * electron_mass);
         q_positron->set_type("QpositronEmission");
         nuc.second->add_property("QpositronEmission", q_positron);
      }
 
      // binding energy
      auto q_binding = make_shared<tkmeasure>(1. / nuc.second->get_a() * (nuc.second->get_z() * (*proton_ME) + nuc.second->get_n() * (*neutron_ME) - *ME));
      q_binding->set_type("binding_energy_overA");
      nuc.second->add_property("binding_energy_overA", q_binding);
 
      // 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")) {
         auto q_pairing = make_shared<tkmeasure>(0.5 * pow(-1, nuc.second->get_n()) * (-2 * (*ME) + *nuc_minus->get("mass_excess") + *nuc_plus->get("mass_excess")));
         q_pairing->set_type("pairingGap");
         nuc.second->add_property("pairingGap", q_pairing);
      }
   }
}
 
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;
   std::copy_if(fvector_of_nuclei.begin(), fvector_of_nuclei.end(), std::back_inserter(res),
                [&_selection](const auto &nuc) { return _selection(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;
}
 
void tkmanager::set_max_level_id(int _id)
{
   std::lock_guard<std::recursive_mutex> lock(get_mutex());
   if (_id > fmax_level_id) fmax_level_id = _id;
}
 
void tkmanager::set_max_decay_id(int _id)
{
   std::lock_guard<std::recursive_mutex> lock(get_mutex());
   if (_id > fmax_decay_id) fmax_decay_id = _id;
}

void tkmanager::load_drip_lines(const tkstring &json_path)
{
   std::ifstream f(json_path);
   if (!f.is_open()) {
      glog << error << "tkmanager::load_drip_lines: cannot open " << json_path << do_endl;
      return;
   }
 
   json j;
   try {
      f >> j;
   } catch (const std::exception &e) {
      glog << error << "tkmanager::load_drip_lines: JSON parse error in "
           << json_path << " : " << e.what() << do_endl;
      return;
   }
 
   if (!j.is_array()) {
      glog << error << "tkmanager::load_drip_lines: root is not an array in "
           << json_path << do_endl;
      return;
   }
 
   fdrip_lines.clear();
 
   auto push_if_present = [&](const json &o, const char *qname, const char *sname) {
      if (!o.contains(qname)) return;
      if (o[qname].is_null()) return;
 
      tkn_drip_point p;
      p.Z = o.value("Z", 0);
      p.N = o.value("N", 0);
      p.value_mev = o.value(qname, std::numeric_limits<double>::quiet_NaN());
      // sigma may be absent → keep NaN
      if (o.contains(sname) && !o[sname].is_null()) {
         p.sigma_mev = o.value(sname, std::numeric_limits<double>::quiet_NaN());
      }
      p.type = qname;
 
      fdrip_lines[p.type].push_back(std::move(p));
   };
 
   try {
      for (const auto &o : j) {
         if (!o.is_object()) continue;
         // Z and N must exist; skip otherwise
         if (!o.contains("Z") || !o.contains("N")) continue;
 
         push_if_present(o, "S1n", "S1n_sigma");
         push_if_present(o, "S2n", "S2n_sigma");
         push_if_present(o, "S1p", "S1p_sigma");
         push_if_present(o, "S2p", "S2p_sigma");
      }
   } catch (const std::exception &e) {
      glog << error << "tkmanager::load_drip_lines: error while iterating JSON: "
           << e.what() << do_endl;
      fdrip_lines.clear();
      return;
   }
 
   // Sort deterministically by Z, then N
   for (auto &kv : fdrip_lines) {
      auto &vec = kv.second;
      std::sort(vec.begin(), vec.end(), [](const tkn_drip_point &a, const tkn_drip_point &b) {
         if (a.Z != b.Z) return a.Z < b.Z;
         return a.N < b.N;
      });
   }
 
   glog << info << "load_drip_lines: loaded " << j.size()
        << " nuclei objects from " << json_path << do_endl;
}

std::vector<tkn_drip_point> tkmanager::get_drip_line(const std::string &_type, bool reduce_per_Z)
{
   std::lock_guard<std::recursive_mutex> lock(get_mutex());
   if (fdrip_lines.empty()) {
      // charge par défaut si pas déjà fait
      load_drip_lines(tkstring::Form("%s/driplines.json", TKN_DB_DIR));
   }
 
   auto it = fdrip_lines.find(_type);
   if (it == fdrip_lines.end()) {
      glog << warning << "get_drip_line: unknown type '" << _type << "'" << do_endl;
      return {};
   }
 
   const auto &vec = it->second;
 
   if (!reduce_per_Z) {
      return vec;
   }
 
   // --- Réduction : garder seulement le point avec min |value| par Z
   std::unordered_map<int, tkn_drip_point> best_by_Z;
   for (const auto &p : vec) {
      auto itZ = best_by_Z.find(p.Z);
      if (itZ == best_by_Z.end()) {
         best_by_Z[p.Z] = p;
      } else {
         const auto &cur = itZ->second;
         if (std::isnan(cur.value_mev) ||
             (!std::isnan(p.value_mev) && std::abs(p.value_mev) < std::abs(cur.value_mev))) {
            best_by_Z[p.Z] = p;
         }
      }
   }
 
   std::vector<tkn_drip_point> reduced;
   reduced.reserve(best_by_Z.size());
   std::transform(best_by_Z.begin(), best_by_Z.end(), std::back_inserter(reduced),
                  [](const auto &kv) { return kv.second; });
 
   std::sort(reduced.begin(), reduced.end(), [](const tkn_drip_point &a, const tkn_drip_point &b) {
      if (a.Z != b.Z) return a.Z < b.Z;
      return a.N < b.N;
   });
 
   return reduced;
}
 
#ifdef HAS_ROOT
ClassImp(tkmanager)
#endif