#include <TF2.h>
#include <TMath.h>
TF2 *resolutionF; //function for the X and Y resolution
TF2 *efficiencyF; //function for the efficiency
TF2 *convolutionintegrandF; //function for the convolution integrand
TF2 *amplitudeF; //function for the amplitude squared
TF2 *conv; //function for the amplitude squared
Double_t amp(Double_t *xv,Double_t *par);
Double_t convolutionintegrand(Double_t *xv,Double_t *par);
Double_t efficiency(Double_t *xv,Double_t *par);
Double_t convolution(Double_t *xv,Double_t *par);


void fit() 
{
   //resolution function
  resolutionF = new TF2("resolutionF","(xgaus(0)+xgaus(3))*(ygaus(6)+ygaus(9))",-0.11,+0.11,-0.16,+0.16); 
  Double_t parres[12];
  parres[0]=6.41771216145510878e+06/3.92099798187156848e+05;
  parres[1]=-8.67549998305522569e-05;
  parres[2]=2.11706881979254226e-02;
  parres[3]=3.62357508006337332e+05/3.92099798187156848e+05;
  parres[4]=-5.61828174605189860e-05;
  parres[5]=6.11327681796687769e-02;
  parres[6]=4.19619853696381208e+06/3.89888112109451613e+05;
  parres[7]=-1.16415240256723986e-03;
  parres[8]=3.17138377573694605e-02;
  parres[9]=2.61556065862704825e+05/3.89888112109451613e+05;
  parres[10]=3.16121669021575576e-03;
  parres[11]=9.44218001135947449e-02;
  resolutionF->SetParameters(parres);

  //make the amplitude function
  amplitudeF=new TF2("amplitudeF",amp,-1,1,-1,1,7); 
//make the convolution integrand function
  convolutionintegrandF=new TF2("convolutionintegrandF",convolutionintegrand,-1.11,1.11,-1.2,1.2,9); 
  efficiencyF=new TF2("efficiencyF",efficiency,-1,1,-1,1); 

  conv = new TF2("conv",convolution,-1.,1.,-1.,1.,7); 
  Double_t start1[7] = {-1.01,0.14,0.08,0.14,-0.426823,0,0};
   conv->SetParameters(start1);
   cout<<Form("Integral = %f",conv->Integral(-1,1,-1,1))<<endl;
   //conv->Draw("colz");
}


Double_t amp(Double_t *xv,Double_t *par)
{
   Double_t x = xv[0];
  Double_t y = xv[1];
  
    return ((1+par[0]*y+par[1]*y*y+par[2]*x+par[3]*x*x+par[4]*x*y+par[5]*y*y*y+par[6]*x*x*y)); 
 
 }

Double_t convolutionintegrand(Double_t *xv,Double_t *par)
{
  Double_t x = xv[0];
  Double_t y = xv[1];
  
  amplitudeF->SetParameters(par);

  return amplitudeF->Eval(par[7]-x,par[8]-y)*resolutionF->Eval(x,y);
  
 }

Double_t efficiency(Double_t *xv,Double_t *par)
{
  Double_t x = xv[0];
  Double_t y = xv[1];
 
 Double_t m0=0.547862;   
Double_t m3=0.1349766;    
 Double_t m1=0.13957018, m2=0.13957018;
  Double_t Q=(m0-m1-m2-m3);
  Double_t T3=Q/3.0*(y+1);
  Double_t T1=Q/6.0*(TMath::Sqrt(3)*x+2-y);
  Double_t T2=Q/6.0*(-TMath::Sqrt(3)*x+2-y);
  Double_t p1=TMath::Sqrt(T1*(T1+2*m1));
  Double_t p2=TMath::Sqrt(T2*(T2+2*m2));
  Double_t p3=TMath::Sqrt(T3*(T3+2*m3));
  Double_t L=TMath::Power(Q+m3-T1-T2,2)-m3*m3-p1*p1-p2*p2;
  Double_t R=2*p1*p2;
  Int_t    tpm0=0; 
  if(TMath::Abs(L)<=R)tpm0=1;
  return tpm0;
}

Double_t convolution(Double_t *xv,Double_t *par)
{
  Double_t x = xv[0];
  Double_t y = xv[1];
  
  convolutionintegrandF->SetParameters(par);
  convolutionintegrandF->SetParameter(7,x);
  convolutionintegrandF->SetParameter(8,y);
  return efficiencyF->Eval(x,y)*convolutionintegrandF->Integral(-0.10585,+0.10585,-0.15857,+0.15857);
}