#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 <TRandom3.h>
#include <TF1.h>
#include <Math/DistSampler.h>
#include <Math/Factory.h>

using namespace std;
// Use the Math namespace
using namespace ROOT::Math;

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]);
 cout << "running with " << nthreads << " threads\n";
 int rseed = 438468301;
 double mg = 0.3, mu = 0.4;
 const double dt = 0.2;
 const int Ncasc = 10000;
 omp_set_num_threads(nthreads);
 ROOT::EnableThreadSafety();
 // initializing the TF1 objects for all threads
 TRandom3* rnd [nthreads];
 TF1* flt [nthreads];
 DistSampler * sampler [nthreads];

 for(int i=0; i<nthreads; i++) {
  rnd[i] = new TRandom3(1231753 + 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);
  sampler[i] = Factory::CreateDistSampler();
  sampler[i]->SetFunction(*flt[i],1);  // last argument: 1-dim
  sampler[i]->SetRandom(rnd[i]);
  sampler[i]->Init();
 }
 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();
   sampler[thread_id]->SetRange(0.,10.);
   bool ret = sampler[thread_id]->Init();
   double lt[3];
   sampler[thread_id]->Sample(lt);
 } // end parallel loop
 return 0;
}