#include "globalesI.h"
#include <gsl/gsl_errno.h>

Int_t i=0;
Double_t par[8] = {42.5, 25.5, 0.45, 0.565, 2.25, 30.0, 0.10, -130.27};
extern Double_t rhopp(Double_t s , /*const*/ Double_t *par);
double func_rhopp(double *x , double * p) { 
return rhopp(x[0],p);
}

//READ DATA

void leerdata(){
  
  ifstream data;
  data.open("pp_elastic.reim");
  
  Int_t i = 0;
  
  if(!data) {
    cout<<"no data file found "<<endl;
  } else {
    while (!data.eof()&&i<117) {
      
      data >> i;
      data >> plab[i];
      data >> plab_min[i];
      data >> plab_max[i];
      data >> sigma[i];
      data >> stat_err1[i];
      data >> stat_err2[i];
      data >> sys_err1[i];
      data >> sys_err2[i];     
      i++;
      
    }
  }
  
  for (i=1; i<117; i++){
    
    elab[i] = sqrt(pow(plab[i],2)+pow(mp,2));
    s[i]=2*elab[i]*mp+pow(mp,2)+pow(mp,2); 
    statistical[i]=0.5*(stat_err1[i]+stat_err2[i]);
    systematic[i]= 0.5*(sys_err1[i]+sys_err2[i])*sigma[i]/100;      
    error_sigma[i]=sqrt(pow(statistical[i],2)+pow(systematic[i],2));
    
  }
  
data.close();
}

Double_t f(Double_t u, void * params){
  Double_t *par =  (Double_t *) params;
  Double_t A = 1-par[7]*u;
  Double_t B = par[7]*par[7]*u*u;
  Double_t C = A*A;
  Double_t D = 1/(2*par[8])* pow((C-B),0.5);
  Double_t f = D*(par[0]*pow(u,par[2]))*(2-(2*par[7]/par[9]))/((par[7]*u-1)*(1-2*par[7]*u+par[9]*u));  
  return f;
}


Double_t f1(Double_t u, void * params){
 Double_t *par =  (Double_t *) params;
 Double_t A = 1-par[7]*u;
 Double_t B = par[7]*par[7]*u*u;
 Double_t C = A*A;
 Double_t D = 1/(2*par[8])*pow((C-B),0.5);
 Double_t f1 = D*(par[5]*(2-(2*par[7])/par[9]))/((par[7]*u-1)*(1-2*par[7]*u+par[9]*u));  
 return f1;
}

Double_t f2(Double_t u, void * params){
Double_t *par =  (Double_t *) params;
Double_t f2; 
Double_t A = 1-par[7]*u*u;
Double_t B = par[7]*par[7]*u*u*u*u;
Double_t C = A*A;
Double_t D = pow ((C-B),0.5);//el momentum
Double_t E = pow(u,2*par[3]+1);//el primer termino con A2
Double_t F = -2*log(u);//logaritmo natural 
Double_t G = pow(F,par[4]);//lo elevo a la gamm
Double_t H = pow(u,2*par[3]-1);//el segundo termino con A2
if (u!=0.0){
       f2 = D*(4.0*par[1]*0.5*par[7]*E+par[6]*u*G*(2*par[9]-2*par[7])-2*par[1]*H)/
         ((par[8]*par[9])*(par[7]*u*u-1)*(u+pow(par[9],-0.5))*(1-2*par[7]*u*u+par[9]*u*u));
       //             cout<<f<<endl;  
}
else f2=0.0;
return f2;
}

Double_t integral(Double_t s, const Double_t * parm){ 

  Double_t par[11];

 par[0] = parm[0];
 par[1] = parm[1];
 par[2] = parm[2];
 par[3] = parm[3];
 par[4] = parm[4];
 par[5] = parm[5];
 par[6] = parm[6];
 par[7] = 2*mp*mp;
 par[8] = mp;
 par[9] = s;
 par[10]= parm[7];
 
gsl_integration_workspace * w = gsl_integration_workspace_alloc (1800);
gsl_integration_workspace * w1 = gsl_integration_workspace_alloc (1800);
gsl_integration_workspace_free(w1);
gsl_integration_workspace * w2 = gsl_integration_workspace_alloc (1800);
gsl_integration_workspace_free(w2);

  int niter = 1000;
  Double_t liminf = 0.0;
  Double_t limsup = 0.25/(mp*mp);
  Double_t limsup1 = 0.532894;
  Double_t errabs = 1e-1;
  Double_t errrel = 1e-1;
  Double_t result, error;
  Double_t result1, error1;
  Double_t result2, error2;
  
  gsl_function F;
  F.function = &f;
  F.params = &par;

  gsl_function F1;
  F1.function = &f1;
  F1.params = &par;
 
  gsl_function F2;
  F2.function = &f2;
  F2.params = &par;

  par[9]=s;
  gsl_integration_qawc (&F, liminf,limsup, 1.0/s, errabs, errrel , niter, w, &result, &error);
  gsl_integration_qawc (&F1, liminf,limsup, 1.0/s, errabs, errrel , niter, w1, &result1, &error1);
  gsl_integration_qawc (&F2, liminf,limsup1, sqrt(1/s), errabs, errrel , niter, w2, &result2, &error2);
  Double_t suma = result+result1+result2; 
  return suma;
}


Double_t sigmateopp (Double_t s, const void* params){
Double_t *par = (Double_t*) params; 
Double_t sigmateopp = par[0]*pow(s,-par[2])-par[1]*pow(s,-par[3])+par[5]+par[6]*pow(log(s),par[4]);
return sigmateopp;
}

//DEFINICION DE LA FUNCION

Double_t rhopp(Double_t s, /*const*/ Double_t *parm){
Double_t B = parm[7]/*-52.27*/;
Double_t p = (1/(2*mp)*sqrt((s-2*ma*ma)*(s-2*ma*ma)-4*mp*mp*ma*ma));
Double_t rhopp = integral(s,parm);
rhopp = (rhopp*((s-ma*ma-mp*mp)/(M_PI*p)))+B/p;
return rhopp/sigmateopp(s,parm);  

}

//______________________________________________________________________________________________________



void chisq()

{

leerdata();

Double_t par[8] = {42.5, 25.5, 0.45, 0.565, 2.25, 30.0, 0.10, -130.27};

Int_t nrhopp = 116;
Double_t chisq = 0;
Double_t delta;
Double_t f;
Double_t s_min = 5;
for (Int_t i=0;i<nrhopp;i++){
      if (sqrt(s[i])>s_min){
         delta  = (sigma[i]-rhopp(s[i],par) )/error_sigma[i];
         chisq += delta*delta;
      }
   }
   f = chisq/116;
cout<<f<<endl;
}  

int main (){
  chisq();
}