#include <string.h>
#include <iostream>
#include <vector>
#include <math.h>
#include "TFile.h"
#include "TH1.h"
#include "TF1.h"
#include <stdio.h>
#include <gsl/gsl_integration.h>

double pCmVal(double m,double m1,double m2){
  double pcMSq = (m*m - pow(m1+m2,2))*(m*m - pow(m1-m2,2))/(4*m*m);
  double pcM = 0;
  if (pcMSq > 0) pcM = sqrt(pcMSq);
  return pcM;
}
double buggPiImagFun(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 buggPiImag = gjSq*rhoj*Fj;

  //buggPiImag = 79-0.0450401*s;

  return buggPiImag;
}

double buggPiRealIntegrand(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)*buggPiImagFun(j,sPrime)/(sPrime - s);
  //integrand = 1.0/(sPrime - s);
  return integrand;
}

struct buggStruct_t {
  int j;
  double s;
};

double f (double sPrime, void * params) {
  buggStruct_t buggStruct = *(buggStruct_t *) params;
  double f = buggPiRealIntegrand(buggStruct.j,buggStruct.s,sPrime);
  return f;
}

void printUsage();

int main (int argc, char **argv){
  char  inFileName[120];
  char  outFileName[120];
  char *argptr;

  sprintf(outFileName,"./myOutFile.root");
  for (int i=1; i<argc; i++) {
    argptr = argv[i];
    if (*argptr == '-') {
      argptr++;
      switch (*argptr) {
      case 'h':
	printUsage();
  	break;
      case 'o':
	strcpy(outFileName,++argptr);
	break;
      default:
	printUsage();
	break;
      }
    } else {
      strcpy(inFileName,argptr);
    }
  }
  TFile *outFile = new TFile(outFileName,"RECREATE");

  int nBinsX = 100;
  double xLow = 0.5;
  double xHigh = 4.0;
  TH1D *hImPi1     = new TH1D("hImPi1","",nBinsX,xLow,xHigh);
  TH1D *hImPi2     = new TH1D("hImPi2","",nBinsX,xLow,xHigh);
  TH1D *hRePi1     = new TH1D("hRePi1","",nBinsX,xLow,xHigh);
  TH1D *hRePi2     = new TH1D("hRePi2","",nBinsX,xLow,xHigh);
  

  ///////////////////////////////  
  buggStruct_t buggStruct;

  int nWin = pow(10,9);

  
  double result, error;
  gsl_function F;
  F.function = &f;
  F.params = &buggStruct;

  double mK = 0.493677;
  double mEta = 0.547862;
  double mPi0 = 0.1349768;
  double sCh1 = pow(mEta+mPi0,2);
  double sCh2 = pow(mK+mK,2);

  std::cout<<"sCh1, sCh2 = "<<sCh1<<" , "<<sCh2<<std::endl;

  //int nBinsX = hImPi1->GetNbinsX();
  for(Int_t i=1;i<=nBinsX;i++){

    gsl_integration_workspace * w = gsl_integration_workspace_alloc (nWin);

    double sVal = hImPi1->GetBinCenter(i);
    double PiImag1 = buggPiImagFun(1,sVal);
    double PiImag2 = buggPiImagFun(2,sVal);
    if (sVal > sCh1) hImPi1->SetBinContent(i,PiImag1);
    if (sVal > sCh2) hImPi2->SetBinContent(i,PiImag2);
    

    buggStruct.s = sVal;
    double sDelta = 0.2;
    double sMax = 10.0;
    double sLowTmp  = sVal-sDelta;
    double sHighTmp = sVal+sDelta;
    

    double PiReal1 = 0;
    double PiReal2 = 0;
    double PiRealError1 = 0;
    double PiRealError2 = 0;
    double sChj = 0;
    //j=1    
    buggStruct.j = 1;
    sChj = sCh1;
    int bType = 2;
    if (sVal < sChj) bType = 0;
    if (sVal > sChj+2*sDelta) bType = 3;
    if (bType == 0){
      //gsl_integration_qags (&F, sCh1, sMax, 0, 1e-3, 100000, w, &result, &error);
      //PiReal1 = result;
      //PiRealError1 = error;
    } 
    if (bType == 3){
      double pts[3] = {sVal-sDelta,sVal,sVal+sDelta};

      std::cout<<"j,sLowA1,sVal,sHigh1 = "<<buggStruct.j<<" , "<<pts[0]<<" , "
               <<pts[1]<<" , "<<pts[2]<<std::endl;

      gsl_integration_qagp (&F, pts, 3, 0.5, 1, nWin,w, &result, &error);
      //PiReal1 = result;
      //PiRealError1 = error;
      //gsl_integration_qags (&F, sVal+sChj, sMax, 1, 1, nWin,w, &result, &error);
      //gsl_integration_qagiu (&F, sCh1+sChj, 0.5, 1, nWin,
      //		     w, &result, &error);
      //PiReal1 += result;
      //PiRealError1 = sqrt(pow(PiRealError1,2)+pow(error,2));
    } 


    hRePi1->SetBinContent(i,PiReal1);
    hRePi1->SetBinError(i,PiRealError1);

    //j=2
    //buggStruct.j = 2;
    //sChj = sCh2;
    //bType = 2;
    //if (sVal < sChj) bType = 0;
    //if (bType == 0){

    //gsl_integration_qags (&F, sCh2, sMax, 0, 1e-3, 100000,w, &result, &error);
    //PiReal2 = result;
    //PiRealError2 = error;
    //} 
    //hRePi2->SetBinContent(i,PiReal2);
    //hRePi2->SetBinError(i,PiRealError2);

    gsl_integration_workspace_free (w);
    
    //std::cout<<"i, PiReal2 = "<<i<<" , "<<PiReal2<<std::endl;
    //double eDiv2 = PiRealError2/PiReal2;

    //bool plotMe = true;
    //if (fabs(eDiv2) > 0.1) plotMe = false;

    //hRePi2->SetBinContent(i,0);
    //hRePi2->SetBinError(i,0);
    //if (plotMe == true){
    //hRePi2->SetBinContent(i,PiReal2);
    //hRePi2->SetBinError(i,PiRealError2);
    //}

  }


  outFile->cd();
  hImPi1->Write();
  hImPi2->Write();
  hRePi1->Write();
  hRePi2->Write();
  outFile->Write();
  outFile->Close();
  return 0;
}

void printUsage(){};