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 <random>
 
#include "tkmanager.h"
#include "tknuclear_chart.h"
 
using namespace tkn;
 
// enumeration of the available decay modes
enum decay_types{kEF,kEP,kBN,kB2N,kB3N,kstable,kproton, kdiproton, kneutron, kdineutron, kalpha, kbeta_plus, kdouble_beta_plus, kdouble_beta_minus ,kbeta_minus, kfission, kunknown};
 
// nucleus structure containing lifetime and first decay mode
struct nucleus {
    int z,n;
    tkstring name;
    bool is_stable = false;
    double lifetime_in_sec;
    decay_types decay_type;
};
 
// 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;
 
// 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;
            a_nuc.decay_type = kstable;
        }
        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();
            if(decay.begins_with("B-;")) a_nuc.decay_type = kbeta_minus;
            else if(decay.begins_with("2B-")) a_nuc.decay_type = kdouble_beta_minus;
            else if(decay.begins_with("B+")) a_nuc.decay_type = kbeta_plus;
            else if(decay.begins_with("EC")) a_nuc.decay_type = kbeta_plus;
            else if(decay.begins_with("2EC")) a_nuc.decay_type = kdouble_beta_plus;
            else if(decay.begins_with("P")) a_nuc.decay_type = kproton;
            else if(decay.begins_with("2P")) a_nuc.decay_type = kdiproton;
            else if(decay.begins_with("N")) a_nuc.decay_type = kneutron;
            else if(decay.begins_with("2N")) a_nuc.decay_type = kdineutron;
            else if(decay.begins_with("A")) a_nuc.decay_type = kalpha;
            else if(decay.begins_with("SF")) a_nuc.decay_type = kfission;
            else if(decay.begins_with("EF")) a_nuc.decay_type = kEF;
            else if(decay.begins_with("EP")) a_nuc.decay_type = kEP;
            else if(decay.begins_with("BN")) a_nuc.decay_type = kBN;
            else if(decay.begins_with("B2N")) a_nuc.decay_type = kB2N;
            else if(decay.begins_with("B3N")) a_nuc.decay_type = kB3N;
            else {
                a_nuc.decay_type = kunknown;
                if(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.decay_type==kstable) {
                cout<<"oups, a stable nucleus should not pass in this loop..."<<endl;
            }
            if(a_nuc.decay_type==kproton) {
                update_chart(a_nuc.z-1,a_nuc.n,Ndecay);
            }
            else if(a_nuc.decay_type==kdiproton) {
                update_chart(a_nuc.z-2,a_nuc.n,Ndecay);
            }
            else if(a_nuc.decay_type==kneutron) {
                update_chart(a_nuc.z,a_nuc.n-1,Ndecay);
            }
            else if(a_nuc.decay_type==kdineutron) {
                update_chart(a_nuc.z,a_nuc.n-2,Ndecay);
            }
            else if(a_nuc.decay_type==kalpha) {
                update_chart(a_nuc.z-2,a_nuc.n-2,Ndecay);
            }
            else if(a_nuc.decay_type==kbeta_plus) {
                update_chart(a_nuc.z-1,a_nuc.n+1,Ndecay);
            }
            else if(a_nuc.decay_type==kEP) {
                update_chart(a_nuc.z-2,a_nuc.n+1,Ndecay);
            }
            else if(a_nuc.decay_type==kdouble_beta_plus) {
                update_chart(a_nuc.z-1,a_nuc.n+1,Ndecay);
            }
            else if(a_nuc.decay_type==kbeta_minus) {
                update_chart(a_nuc.z+1,a_nuc.n-1,Ndecay);
            }
            else if(a_nuc.decay_type==kdouble_beta_minus) {
                update_chart(a_nuc.z+2,a_nuc.n-2,Ndecay);
            }
            else if(a_nuc.decay_type==kBN) {
                update_chart(a_nuc.z+1,a_nuc.n-2,Ndecay);
            }
            else if(a_nuc.decay_type==kB2N) {
                update_chart(a_nuc.z+1,a_nuc.n-3,Ndecay);
            }
            else if(a_nuc.decay_type==kB3N) {
                update_chart(a_nuc.z+1,a_nuc.n-4,Ndecay);
            }
            // very approximative representation of the fission...
            else if(a_nuc.decay_type==kfission || a_nuc.decay_type==kEF) {
                for(int i=0 ; i<Ndecay ; i++) {
                    while(true) {
                        int newz1 = TMath::Nint(gRandom->Uniform(a_nuc.z/4,a_nuc.z/2));
                        int newn1 = TMath::Nint(gRandom->Uniform(a_nuc.n/4,a_nuc.n/2));
                        int emitted_neutrons = TMath::Nint(gRandom->Uniform(2,7));
                        int newz2 = a_nuc.z-newz1;
                        if(a_nuc.decay_type==kEF) newz2 -= 1;
                        int newn2 = a_nuc.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_type==kunknown) {
            }
            else {
                cout<<"unknown decay for nucleus " << a_nuc.name <<endl;
            }
        }
    }
    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(double i=dt ; i<_time ; i+=dt) {
        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++");
}