tkn-thread-example.cpp

#include <thread>
#include <getopt.h>
#include <random>
 
#include "tkmanager.h"
 
#ifdef HAS_ROOT
#include "TH2F.h"
#include "TFile.h"
#include "TROOT.h"
#endif
 
#define no_argument 0
#define required_argument 1
#define optional_argument 2
 
using namespace tkn;
using namespace std;
 
#ifdef HAS_ROOT
TH2F *hNucChart = nullptr;
TH2F *hNucChartBetaminus = nullptr;
TH2F *hNucChartBetaplus = nullptr;
TH2F *hNucChartStable = nullptr;
TH1F *hRandExcitedState = nullptr;
#endif
 
void a_thread(int ithread, long NEvents);
void print_help();
 
// To compile this code, use: g++ tkn-thread-example.cpp -o tkn-thread `tkn-config --cflags --linklibs`
 
void print_help() {
    glog << blue << high_intensity << " *******************************************" << do_endl;
    glog << blue << high_intensity << " ***    tkn-thread-example user-guide    ***" << do_endl;
    glog << blue << high_intensity << " *******************************************" << do_endl;
    glog << do_endl;
    glog << "this utility allows to test the multi thread compatibility" << do_endl;
    glog << do_endl;
    glog  << "supported options:" << do_endl;
    glog << " --evts N     number of random nuclei to process (default 1e6)"  << do_endl;
    glog << " --workers N  number of workers to be used (default 1)"  << do_endl;
    glog << do_endl;
}
 
int main(int argc, char **argv) {
 
    int option_index = 0;
    static struct option long_options[] = {
        {"evts",       required_argument,       nullptr,    'e'},
        {"workers",    required_argument,       nullptr,    'w'},
        {"help",       no_argument,             nullptr,    'h'},
        {nullptr,      0,                       nullptr,    0}};
 
    long NEvents=1e6;
    int NWorkers = 1;
 
    while (true) {
        int c = getopt_long (argc, argv, "ewh", long_options, &option_index);
        if (c == -1) break;
        switch (c) {
        case 'e':
            NEvents = atof(optarg);
            break;
        case 'w':
            NWorkers = atoi(optarg);
            break;
        case 'h':
            print_help();
            exit(EXIT_SUCCESS);
        }
    }
    cout << "***********************************************************************************" << endl;
    cout << "* Randomly produce a nucleus and extract a random excited state in multi thread   *" << endl;
    cout << "***********************************************************************************" << endl;
    cout<<endl;
 
#ifdef HAS_ROOT
    ROOT::EnableImplicitMT();
    hNucChart = new TH2F("NucChart","NucChart",90,0,180,60,0,120);
    hNucChartBetaminus = new TH2F("NucChartBetaminus","NucChartBetaminus",90,0,180,60,0,120);
    hNucChartBetaplus = new TH2F("NucChartBetaplus","NucChartBetaplus",90,0,180,60,0,120);
    hNucChartStable = new TH2F("NucChartStable","NucChartStable",90,0,180,60,0,120);
    hRandExcitedState = new TH1F("RandExcitedState","RandExcitedState",20000,0,20000);
#endif
 
    cout << "Nevents  : " << NEvents << endl;
    cout << "N workers: " << NWorkers << endl;
 
    // to remove printouts due to empty level schemes
    glog.set_warnings(false);
 
    auto tsload = std::chrono::high_resolution_clock::now();
 
    // preload one all level schemes for fair multi thread comparison
    gmanager->preload_level_schemes(true);
    cout<<endl;
 
    auto tstart = std::chrono::high_resolution_clock::now();
 
    vector<thread> threads;
    threads.reserve(NWorkers);
    for(int i=0 ; i<NWorkers ; i++) {
        threads.emplace_back(thread(a_thread, i, NEvents/NWorkers));
    }
    for(int i=0 ; i<NWorkers ; i++) {
        threads.at(i).join();
    }
 
#ifdef HAS_ROOT
    auto *fout = new TFile("tkn-thread.root","recreate");
    hNucChart->Write();
    hNucChartBetaminus->Write();
    hNucChartBetaplus->Write();
    hNucChartStable->Write();
    hRandExcitedState->Write();
    glog << info << "tkn-thread.root created with content: " << do_endl;
    fout->ls();
    fout->Close();
#endif
 
    auto tstop = std::chrono::high_resolution_clock::now();
    auto dtload = std::chrono::duration<double, std::milli>(tstart-tsload).count() * 0.001;
    auto dt_process = std::chrono::duration<double, std::milli>(tstop-tstart).count() * 0.001;
    auto dt = std::chrono::duration<double, std::milli>(tstop-tsload).count() * 0.001;
 
    cout<<endl;
    cout << fixed << setprecision(4) << "CPU time to load the full db  : " << dtload << " s\n";
    cout << fixed << setprecision(4) << "CPU time to process the events: " << dt_process << " s\n";
    cout << fixed << setprecision(4) << "Total CPU time                : " << dt << " s\n";
}
 
void a_thread(int ithread, long NEvents) {
 
    // Create a random number generator engine
    std::random_device rd;
    std::mt19937 gen(rd());
 
    // Create a uniform integer distribution on the number of nuclei
    std::uniform_int_distribution<> dist_nuc(0, gmanager->get_nuclei().size()-1);
 
    int Fact = NEvents/100.;
    for(long entry=0 ; entry<NEvents ; entry++) {
        if(ithread==0 && entry%Fact==0) {
            cout<<"\r"<<"Analysis progress : "<<setw(3)<<setprecision(3)<<(double)entry/(double)NEvents*100<<" %"<<flush;
        }
        long inuc = dist_nuc(gen);
        auto nuc = gmanager->get_nuclei().at(inuc);
        if(nuc->get_level_scheme()->get_levels().size()>1) {
 
            // Create a uniform integer distribution on the number of levels
            std::uniform_int_distribution<> dist_lev(0, nuc->get_level_scheme()->get_levels().size()-1);
 
            long ilvl = dist_lev(gen);
            double elev = nuc->get_level_scheme()->get_levels().at(ilvl)->get_energy(tkn::tkunit_manager::keV,true);
#ifdef HAS_ROOT
            hNucChart->Fill(nuc->get_n(),nuc->get_z());
            if(nuc->has_property("QbetaMinus") && nuc->get("QbetaMinus")->get_value()>0) hNucChartBetaminus->Fill(nuc->get_n(),nuc->get_z());
            if(nuc->has_property("QpositronEmission") && nuc->get("QpositronEmission")->get_value()>0) hNucChartBetaplus->Fill(nuc->get_n(),nuc->get_z());
            if(nuc->is_stable()) hNucChartStable->Fill(nuc->get_n(),nuc->get_z());
            if(elev>0) hRandExcitedState->Fill(elev);
#endif
        }
    }
}