decay_generator.C

#include "Riostream.h"
#include "TH2F.h"
#include "TROOT.h"
#include "TMath.h"
#include "TRandom3.h"
#include "TCanvas.h"
#include "TSystem.h"
#include "TLatex.h"
#include "TStyle.h"
 
#include <cctype>
#include <random>
 
#include "tkmanager.h"
#include "tknuclear_chart.h"
 
using namespace tkn;
 
struct decay_step {
   bool known = false;
   bool is_fission = false;
   int z_offset = 0;
   int n_offset = 0;
};
 
// nucleus structure containing lifetime and first decay mode
struct nucleus {
   int z, n;
   tkstring name;
   bool is_stable = false;
   double lifetime_in_sec;
   decay_step decay;
};
 
// list containing all the known nuclei
vector<nucleus> nucleus_list;
// mapping of the histogram containing the generated nuclei and associated decays
map<int, map<int, pair<int, nucleus>>> nuc_chart_mapping;
// tknuclear chart object
tknuclear_chart *tknuc_chart = nullptr;
// total elapsed time
double felapsed_time = 0;
 
bool read_emitted_particle_count(const tkstring &mode, const tkstring &particle, int &count)
{
   if (mode == particle) {
      count = 1;
      return true;
   }
 
   if (!mode.ends_with(particle)) return false;
   tkstring number = mode.substr(0, mode.length() - particle.length());
   if (!number.is_digit()) return false;
 
   count = number.atoi();
   return count > 0;
}
 
bool parse_cluster_emission(const tkstring &mode, int &z_offset, int &n_offset)
{
   if (mode.empty() || !std::isdigit(static_cast<unsigned char>(mode.front()))) return false;
 
   tkstring mass_str = mode.copy().remove_alpha();
   if (!mass_str.is_digit()) return false;
 
   tkstring symbol = mode.copy().extract_alpha();
   symbol.to_lower().capitalize();
 
   int cluster_z = 0;
   if (!gmanager->known_element(symbol, cluster_z)) return false;
 
   const int cluster_a = mass_str.atoi();
   if (cluster_a <= cluster_z) return false;
 
   z_offset -= cluster_z;
   n_offset -= cluster_a - cluster_z;
   return true;
}
 
bool apply_particle_emission(const tkstring &mode, decay_step &step)
{
   if (mode.empty()) return true;
 
   int count = 0;
   if (read_emitted_particle_count(mode, "N", count)) {
      step.n_offset -= count;
      return true;
   }
   if (read_emitted_particle_count(mode, "P", count)) {
      step.z_offset -= count;
      return true;
   }
   if (mode == "A") {
      step.z_offset -= 2;
      step.n_offset -= 2;
      return true;
   }
   if (mode == "F" || mode == "SF") {
      step.is_fission = true;
      return true;
   }
   if (parse_cluster_emission(mode, step.z_offset, step.n_offset)) return true;
 
   // NuDat sometimes reports heavy-cluster emission without a mass number
   // (for example "Si"). The animation only needs a qualitative branch.
   tkstring symbol = mode;
   symbol.to_lower().capitalize();
   int cluster_z = 0;
   if (gmanager->known_element(symbol, cluster_z) && cluster_z > 2) {
      step.is_fission = true;
      return true;
   }
 
   return false;
}
 
decay_step decode_decay_mode(tkstring mode)
{
   decay_step step;
   mode.remove_all(" ");
   mode.to_upper();
 
   if (mode.empty()) return step;
 
   if (mode == "B-") {
      step.known = true;
      step.z_offset = 1;
      step.n_offset = -1;
      return step;
   }
   if (mode == "2B-") {
      step.known = true;
      step.z_offset = 2;
      step.n_offset = -2;
      return step;
   }
   if (mode.begins_with("B-")) {
      step.known = true;
      step.z_offset = 1;
      step.n_offset = -1;
      if (apply_particle_emission(mode.substr(2), step)) return step;
   }
   if (mode.begins_with("B") && mode.length() > 1 && mode != "B+") {
      decay_step legacy_step;
      legacy_step.known = true;
      legacy_step.z_offset = 1;
      legacy_step.n_offset = -1;
      if (apply_particle_emission(mode.substr(1), legacy_step)) return legacy_step;
   }
 
   if (mode == "B+" || mode == "EC" || mode == "EC+B+") {
      step.known = true;
      step.z_offset = -1;
      step.n_offset = 1;
      return step;
   }
   if (mode == "2EC") {
      step.known = true;
      step.z_offset = -2;
      step.n_offset = 2;
      return step;
   }
   if (mode.begins_with("EC")) {
      step.known = true;
      step.z_offset = -1;
      step.n_offset = 1;
      if (apply_particle_emission(mode.substr(2), step)) return step;
   }
   if (mode.begins_with("E") && mode.length() > 1) {
      decay_step legacy_step;
      legacy_step.known = true;
      legacy_step.z_offset = -1;
      legacy_step.n_offset = 1;
      if (apply_particle_emission(mode.substr(1), legacy_step)) return legacy_step;
   }
 
   if (apply_particle_emission(mode, step)) {
      step.known = true;
      return step;
   }
 
   step.known = false;
   step.is_fission = false;
   step.z_offset = 0;
   step.n_offset = 0;
   return step;
}
 
// nuclei list initialisation
void make_db()
{
   if (tknuc_chart) return;
   tknuc_chart = new tknuclear_chart("Decay generator", tknuclear_chart::kAll, true);
 
   // loop on all the nuclei known in tkn and fill the nucleus list
   for (auto &nuc : gmanager->get_nuclei()) {
      nucleus a_nuc;
      a_nuc.name = nuc->get_symbol();
      a_nuc.z = nuc->get_z();
      a_nuc.n = nuc->get_n();
      if (nuc->is_stable()) {
         a_nuc.is_stable = true;
      } else {
         if (!nuc->get_lifetime_measure())
            a_nuc.lifetime_in_sec = 0.;
         else
            a_nuc.lifetime_in_sec = nuc->get_lifetime();
 
         tkstring decay = nuc->get_property("decay_modes");
         auto token = decay.replace_all("]", "[").tokenize("[");
         if (token.size() == 0)
            decay = "";
         else
            decay = token.front();
         auto decay_tokens = decay.tokenize(";");
         tkstring decay_mode = decay_tokens.size() ? decay_tokens.front() : "";
         a_nuc.decay = decode_decay_mode(decay_mode);
         if (!a_nuc.decay.known && decay.length()) cout << "unknown decay: " << decay << " " << a_nuc.name << endl;
      }
      nucleus_list.push_back(a_nuc);
      nuc_chart_mapping[a_nuc.z][a_nuc.n] = make_pair(0, a_nuc);
   }
}
 
// initialisation of the nuclear chart using NEvents nuclei, taken randomly
void init(Int_t NEvents = 1e6)
{
   if (nucleus_list.empty()) make_db();
 
   // reset map
   for (auto &z : nuc_chart_mapping) {
      for (auto &n : z.second) {
         n.second.first = 0;
      }
   }
   tknuc_chart->reset();
 
   std::default_random_engine generator;
   uniform_int_distribution<Int_t> dist(0, nucleus_list.size() - 1);
 
   for (Int_t inuc = 0; inuc < NEvents; inuc++) {
      nucleus a_nuc = nucleus_list.at(dist(generator));
      nuc_chart_mapping[a_nuc.z][a_nuc.n].first += 1;
      tknuc_chart->fill(a_nuc.z, a_nuc.n);
   }
 
   tknuc_chart->set_title("T = 0");
   tknuc_chart->draw();
   felapsed_time = 0.;
 
   gSystem->Unlink("animation.gif");
   tknuc_chart->save_as("animation.gif+");
}
 
// update the map and histogram after a decay
void update_chart(int z, int n, int N)
{
   if (!(nuc_chart_mapping.count(z) && nuc_chart_mapping[z].count(n))) {
      return;
   }
   nuc_chart_mapping[z][n].first += N;
   tknuc_chart->set_value(z, n, nuc_chart_mapping[z][n].first);
}
 
// process a decay of all the nuclei of the nuclear chart considering a time dt
void process_decay(double dt, bool save = false)
{
 
   if (tknuc_chart == nullptr) {
      cout << "process first the init step to generate the T0 distribution" << endl;
      return;
   }
 
   // first, determine how many decay per nucleus needs to be done
   map<int, map<int, int>> decay_to_process;
   for (auto &z : nuc_chart_mapping) {
      for (auto &n : z.second) {
         auto &pair = n.second;
         if (pair.first == 0) continue;
         auto &a_nuc = pair.second;
         if (a_nuc.is_stable) continue;
         double Ndecay_double = pair.first - pair.first * exp(-TMath::Log(2) * dt / a_nuc.lifetime_in_sec);
         int Ndecay = TMath::Nint(Ndecay_double);
         // random decay when the number of decay to process is lower than 1
         if (Ndecay_double < 1) {
            if (gRandom->Uniform() < Ndecay_double) Ndecay = 1;
         }
         decay_to_process[z.first][n.first] = Ndecay;
      }
   }
   for (auto &z : decay_to_process) {
      for (auto &n : z.second) {
         auto &a_nuc = nuc_chart_mapping[z.first][n.first].second;
         int Ndecay = n.second;
         update_chart(a_nuc.z, a_nuc.n, -Ndecay);
         if (a_nuc.is_stable) {
            cout << "oups, a stable nucleus should not pass in this loop..." << endl;
         }
         // very approximative representation of the fission...
         if (a_nuc.decay.is_fission) {
            const int fission_z = a_nuc.z + a_nuc.decay.z_offset;
            const int fission_n = a_nuc.n + a_nuc.decay.n_offset;
            for (int i = 0; i < Ndecay; i++) {
               while (true) {
                  int newz1 = TMath::Nint(gRandom->Uniform(fission_z / 4, fission_z / 2));
                  int newn1 = TMath::Nint(gRandom->Uniform(fission_n / 4, fission_n / 2));
                  int emitted_neutrons = TMath::Nint(gRandom->Uniform(2, 7));
                  int newz2 = fission_z - newz1;
                  int newn2 = fission_n - newn1 - emitted_neutrons;
                  if (!(nuc_chart_mapping.count(newz1) && nuc_chart_mapping[newz1].count(newn1)) || !(nuc_chart_mapping.count(newz2) && nuc_chart_mapping[newz2].count(newn2))) continue;
                  update_chart(newz1, newn1, 1);
                  update_chart(newz2, newn2, 1);
                  break;
               }
            }
         } else if (a_nuc.decay.known) {
            update_chart(a_nuc.z + a_nuc.decay.z_offset, a_nuc.n + a_nuc.decay.n_offset, Ndecay);
         } else {
            // Keep the historical behaviour for nuclei without usable decay mode:
            // remove the decayed parent and do not feed a daughter nucleus.
         }
      }
   }
   felapsed_time += dt;
 
   // if save, build a gif
   if (save) {
      tkstring title = "T = ";
      if (felapsed_time < 60)
         title += tkstring::form("%.1f sec", felapsed_time);
      else if (felapsed_time < 60 * 60)
         title += tkstring::form("%.1f min", felapsed_time / 60);
      else if (felapsed_time < 60 * 60 * 24)
         title += tkstring::form("%.1f h", felapsed_time / 60 / 60);
      else if (felapsed_time < 60 * 60 * 24 * 365)
         title += tkstring::form("%.1f days", felapsed_time / 60 / 60 / 24);
      else
         title += tkstring::form("%.1g years", felapsed_time / 60 / 60 / 24 / 365);
 
      tknuc_chart->set_title(title);
      tknuc_chart->update();
      tknuc_chart->save_as("animation.gif+");
 
      gSystem->ProcessEvents();
   }
}
 
// calculate the decay for _time seconds, with a delta t dt. Process nsave pictures in the gif
void elapse_time_in_seconds(double _time, double dt = 0., int nsave = 5)
{
 
   if (tknuc_chart == nullptr) {
      cout << "process first the init step to generate the T0 distribution" << endl;
      return;
   }
 
   if (dt == 0.) dt = _time / 1000.;
   int ndt = 1;
   int save = (_time / dt) / nsave;
   for (int i = 1; i * dt < _time; i++) {
      if (ndt % save == 0)
         process_decay(dt, true);
      else
         process_decay(dt);
      ndt++;
   }
   process_decay(dt, true);
}
 
void elapse_time_in_hours(double time, double dt = 0.)
{
 
   time = time * 3600.;
   if (dt != 0.) dt = dt * 3600.;
   int nsave = 4;
   elapse_time_in_seconds(time, dt, nsave);
}
 
void elapse_time_in_days(double time, double dt = 0.)
{
 
   time = time * 3600. * 24;
   if (dt != 0.) dt = dt * 3600. * 24;
   int nsave = 2;
   elapse_time_in_seconds(time, dt, nsave);
}
 
void elapse_time_in_years(double time, double dt = 0.)
{
 
   time = time * 3600. * 24 * 365;
   if (dt != 0.) dt = dt * 3600. * 24 * 365;
   int nsave = 1;
   elapse_time_in_seconds(time, dt, nsave);
}
 
// decay generation during 5e9 years
void decay_generator_animation()
{
   init(1e9);
   elapse_time_in_seconds(1);       // 1sec
   elapse_time_in_seconds(9);       // 10sec
   elapse_time_in_seconds(20);      // 30sec
   elapse_time_in_seconds(30);      // 1min
   elapse_time_in_seconds(60);      // 2min
   elapse_time_in_seconds(60 * 8);  // 10min
   elapse_time_in_seconds(60 * 20); // 30min
   elapse_time_in_seconds(60 * 30); // 1h
   elapse_time_in_hours(1);         // 2h
   elapse_time_in_hours(8);         // 10h
   elapse_time_in_hours(14);        // 1d
   elapse_time_in_days(1);          // 2d
   elapse_time_in_days(8);          // 10d
   elapse_time_in_days(40);         // 50d
   elapse_time_in_days(50);         // 100d
   elapse_time_in_days(100);        // 200d
   elapse_time_in_days(165);        // 1y
   elapse_time_in_years(1);         // 2y
   elapse_time_in_years(4);         // 5y
   elapse_time_in_years(5);         // 10y
   elapse_time_in_years(40);        // 50y
   elapse_time_in_years(50);        // 100y
   elapse_time_in_years(400);       // 500y
   elapse_time_in_years(500);       // 1000y
   elapse_time_in_years(9000);      // 10000y
   elapse_time_in_years(40000);     // 50000y
   elapse_time_in_years(50000);     // 100000y
   elapse_time_in_years(400000);    // 500000y
   elapse_time_in_years(500000);    // 1e6y
   elapse_time_in_years(4e6);       // 5e6y
   elapse_time_in_years(5e6);       // 1e7y
   elapse_time_in_years(4e7);       // 5e7y
   elapse_time_in_years(5e7);       // 1e8y
   elapse_time_in_years(4e8);       // 5e8y
   elapse_time_in_years(5e8);       // 1e9y
   elapse_time_in_years(4e9);       // 1e9y
   tknuc_chart->save_as("animation.gif++");
}