#include <iostream>
#include <fstream>
#include <stdio.h>
#include <sstream>
#include <stdarg.h>
#include <cmath>
#include <iomanip>
#include <cstdlib>
#include <omp.h>
#include <TROOT.h>
#include <TF1.h>

double ff_lt(double* x, double* par)
{
 double  lt2 = x[0]*x[0]; // l_t^2
 double& mg = par[0];    // m_g,min
 double& mu = par[1];   // \mu
 return x[0] * log(1. + lt2/(exp(1.0) * mg*mg)) / pow(lt2 + mu*mu, 2);
}

int main(int argc, char **argv)
{
 int nthreads = atoi(argv[1]);
 std::cout << "running with " << nthreads << " threads\n";
 int rseed = 438468301;
 double mg = 0.3, mu = 0.4;
 const int Ncasc = 10000;
 omp_set_num_threads(nthreads);
 ROOT::EnableThreadSafety();
 // initializing the TF1 objects for all threads
 TF1* flt [nthreads];
 for(int i=0; i<nthreads; i++) {
  char flt_ROOT_name [20];
  sprintf(flt_ROOT_name, "dsigmadlt_%i", i);
  flt[i] = new TF1(flt_ROOT_name, ff_lt, 0., 10.0, 2);
 }
 std::cout << "init done\n";
 // parallel loop
 #pragma omp parallel for schedule(static, 1000)
 for(int icasc=0; icasc<Ncasc; icasc++) {
   int thread_id = omp_get_thread_num();
   flt[thread_id]->SetParameters(mg, mu);
   double lt = flt[thread_id]->GetRandom();
 } // end parallel loop
 return 0;
}