#include "Math/Functor.h"
#include "Math/Integrator.h"
#include "Math/WrappedMultiTF1.h"
#include "Math/WrappedParamFunction.h"
#include "Math/GSLMCIntegrator.h"
#include "Math/AdaptiveIntegratorMultiDim.h"
#include "Math/IOptions.h"
#include <ROOT/RDataFrame.hxx>
#include "TStopwatch.h"

using std::cout;
using std::endl;

ROOT::Math::Integrator * ig_f_1 = nullptr;
ROOT::Math::Integrator * ig_f = nullptr;
   

auto sqr(double x)
{ return x * x; }

double f_1(double x) {
    Double_t par[3]={0.504536,0.0671194,6829.58};
    Double_t val = (sqr(par[0])*par[1])/sqr(par[1]-par[0])*(exp(-par[0]*x)*((par[1]-par[0])*x-1)+exp(-par[1]*x))*par[2];
    return val;
}

double fint_1()
{
   
   
    Double_t T=60.0;

    
    Double_t fintval = ig_f_1->Integral(0.0,T);
    
    return fintval;
}

double f(double x) {
    Double_t p0=0.504536, p1=0.0671194, p2=6829.58, p3=fint_1();
    return ((sqr(p0)*p1/sqr(p1-p0))*(exp(-p0*x)*((p1-p0)*x-1.0)+exp(-p1*x))*p2)/p3;
}
//____________________________

double fint1(double t)
{
  
    Double_t tau=0.9;

    ROOT::Math::Functor1D wf_f(&f);
    ig_f->SetFunction(wf_f);
    
    Double_t fintval = ig_f->Integral(t,t+tau);
    
    return fintval;
}
//_______________________________i=6-------------------------------------------------------

double ff6(const double *x, const double *par)
{
    const int i=6;
    Double_t tau=0.9;
    Double_t T=60.0;
    Double_t t2min=x[0]+tau;
    Double_t t2max=T-(i-2.0)*tau;
    Double_t t2=(t2max-t2min)*x[1]+t2min;
    Double_t dt2=t2max-t2min;
    Double_t t3min=t2+tau;
    Double_t t3max=T-(i-3.0)*tau;
    Double_t t3=(t3max-t3min)*x[2]+t3min;
    Double_t dt3=t3max-t3min;
    Double_t t4min=t3+tau;
    Double_t t4max=T-(i-4.0)*tau;
    Double_t t4=(t4max-t4min)*x[3]+t4min;
    Double_t dt4=t4max-t4min;
    Double_t t5min=t4+tau;
    Double_t t5max=T-(i-5.0)*tau;
    Double_t t5=(t5max-t5min)*x[4]+t5min;
    Double_t dt5=t5max-t5min;
    Double_t t6min=t5+tau;
    Double_t t6max=T-(i-6.0)*tau;
    Double_t t6=(t6max-t6min)*x[5]+t6min;
    Double_t dt6=t6max-t6min;
    Double_t ffval = (TMath::Factorial(par[0])/TMath::Factorial(par[0]-i)) *f(x[0]) *f(t2) *f(t3) *f(t4) *f(t5) *f(t6) *dt2 *dt3 *dt4 *dt5 *dt6 *pow(fint1(x[0]) +fint1(t2) +fint1(t3) +fint1(t4) +fint1(t5) +fint1(t6),par[0]-i);
    return ffval;
}
double S6j(double j)
{
    const int i=6;
    Double_t tau=0.9;
    Double_t T=60.0;
    Double_t epsabs=1e-2;
    ROOT::Math::WrappedParamFunction<> wf( &ff6, i, 1);
    Double_t p[1] = {j};
    wf.SetParameters(p);
    // Create the Integrator
    ROOT::Math::GSLMCIntegrator ig(ROOT::Math::MCIntegration::kVEGAS);
    //ROOT::Math::GSLMCIntegrator ig(ROOT::Math::MCIntegration::kMISER);
    //ROOT::Math::AdaptiveIntegratorMultiDim ig;
    // Set parameters of the integration
    ig.SetFunction(wf);
    //ig.SetAbsTolerance(1.E-1);
    ig.SetRelTolerance(1.E-3);
    ROOT::Math::VegasParameters params;
    params.iterations = 2;
    params.verbose = 1;
    ig.SetParameters(params);
 
    
    ig.ExtraOptions()->Print();

    

    
    Double_t xmin[i] = {0.0,0.0,0.0,0.0,0.0,0.0};
    Double_t xmax[i] = {T-(i-1.0)*tau,1.0,1.0,1.0,1.0,1.0};

    // create internal integral
    ig_f_1 = new ROOT::Math::Integrator(ROOT::Math::IntegrationOneDim::kADAPTIVESINGULAR);
    ig_f = new ROOT::Math::Integrator(ROOT::Math::IntegrationOneDim::kADAPTIVESINGULAR);
    
    ROOT::Math::Functor1D wf_f_1(&f_1);
    ig_f_1->SetFunction(wf_f_1);
    
    ROOT::Math::Functor1D wf_f(&f);
    ig_f->SetFunction(wf_f);
    
    
    Double_t Sjvalue =ig.Integral(xmin, xmax);
    cout << "Error: " << ig.Error() << endl;
    cout << "Matrix element S66: ";

    delete ig_f_1;
    delete ig_f;
    
    return Sjvalue;
}
TStopwatch timer;
void spencerIntegrals_rootforumtest() {
    ROOT::EnableImplicitMT(0);
    timer.Start();
    cout << S6j(6) << endl;
    timer.Stop();
    printf("---RT=%7.5f m, RT=%7.5f h\n",timer.RealTime()/60.0,timer.RealTime()/60.0/60.0);
}

int main() {
   spencerIntegrals_rootforumtest();
}