#include "Math/AdaptiveIntegratorMultiDim.h"
//#include "Math/GSLIntegrator.h"

double pCmVal(double m,double m1,double m2){
  double pcM = sqrt(fabs((m*m - pow(m1+m2,2))*(m*m - pow(m1-m2,2))))/(2*m);
  return pcM;
}

double bugPiImagFun(int j,double s){
  double mK = 0.493677;
  double mEta = 0.547862;
  double mPi0 = 0.1349768;
  double beta = 2;
  double m0 = 0.990;   //From Tbl. IV, Phys.Rev.D 95, 032002 (2017)
  double qj = 0;
  double gjSq = 0;
  if (j == 1){
    qj = pCmVal(sqrt(s),mEta,mPi0);
    gjSq = 0.341;       //From Tbl. IV, Phys.Rev.D 95, 032002 (2017)
  }
  if (j == 2){
    qj = pCmVal(sqrt(s),mK,mK);
    gjSq = 0.892*0.341; //From Tbl. IV, Phys.Rev.D 95, 032002 (2017)
  }
  double rhoj = 2*qj/sqrt(s);
  double Fj = exp(-beta*qj*qj);
  double bugPiImag = gjSq*rhoj*Fj;
  
  return bugPiImag;
}

double bugPiRealIntegrand(int j,double s,double sPrime){
  double mK = 0.493677;
  double mEta = 0.547862;
  double mPi0 = 0.1349768;
  double sj = 0;
  if (j == 1){
    sj = pow(mEta+mPi0,2);
  }
  if (j == 2){
    sj = pow(mK+mK,2);
  }
  double integrand = (1.0/3.14159)*bugPiImagFun(j,sPrime)/(sPrime - s);
  return integrand;
}

class  MyFitFunBugIntegrand {
public:
  double Evaluate(double *sPrimePtr, double *parPtr) {
    double sPrime = *sPrimePtr;
    int j = (int)parPtr[0];
    int s = parPtr[1];
    double fitVal = bugPiRealIntegrand(j,s,sPrime);
    return fitVal;
  }
};

void bug(int iShow){
  gROOT->Reset();
  gROOT->SetStyle("Plain");
  gStyle->SetOptStat(0);


  MyFitFunBugIntegrand * fptr1 = new MyFitFunBugIntegrand();
  int nPar = 2;
  TF1 *f1 = new TF1("f1",fptr1,&MyFitFunBugIntegrand::Evaluate,0.0,4000.0,nPar,
		    "MyFitFunBugIntegrand","Evaluate");


  TH1D *hBug = new TH1D("hBug","",300,1.2,10);

  ROOT::Math::IntegratorMultiDimOptions::SetDefaultIntegrator("AdaptiveSingular");
  //cout<<ROOT::Math::IntegratorMultiDimOptions::DefaultIntegrator()<<endl;;
  f1->FixParameter(0,1);
  f1->FixParameter(1,1.1);
  hBug->FillRandom("f1",100000);

  //double test = f1->IntegralCauchy(1.0,10,1.1);
  //cout<<"test = "<<test<<endl;
  
  if (iShow == 1) {
    hBug->Draw();
  }
}

void cTest(int iShow){
  //ROOT::Math::GSLInetagrator ig;
  ROOT::Math::IntegratorOneDim ig;
  double c = 2;
  auto f = [&](double x){return x*(x-c);};


  //ig.SetFunction(f);
  //double value = ig.IntegralCauchy(0,10,c);
  //cout<<"int = "<<value<<endl;

}