#include <TTree.h>
#include <TBranch.h>
#include <TROOT.h>
#include <TChain.h>
#include <TSystem.h>
#include <TFile.h>
#include <TH2F.h>
#include <TF1.h>
#include <TH2.h>
#include <TH2I.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TMath.h>
#include <TRandom.h>
#include <TNamed.h>
#include <TGraph.h>
#include <TGraphErrors.h>
#include <TGraph2DErrors.h>
#include <TCut.h>
#include <TGraph2D.h>
#include <TMultiGraph.h>
#include <TF2.h>
#include <TLegend.h>

#include <TMinuit.h>

#include <string>
#include <iostream>
#include <iomanip>
#include <fstream>
#include <cmath>

using namespace std;

Double_t fitf(Double_t *x, Double_t *par) 
{ 
  //x[0] contains the x value 
  //x[1] contains the y value 
  //par[6] 
  // [0]+[1]*x+[2]*x*x+[3]*y+[4]*y*y+[5]*x*y 
  Double_t fitval(0); 
  
  if ( sqrt(x[0]*x[0]+x[1]*x[1])<90. ) 
  { 
    fitval=par[0]+par[1]*x[0]+par[2]*x[0]*x[0]+par[3]*x[1]+par[4]*x[1]*x[1]+par[5]*x[0]*x[1];    
  } else {fitval=0;} 
  

  return fitval; 
} 

void macro()
{
  gStyle->SetPalette(1,0);
  
  TFile *a = new TFile("graph.root");
  TGraph2D *g = (TGraph2D*)a->Get("Graph2D");

  //TFile *a = new TFile("grapherr.root");  
  //TGraph2DErrors *g = (TGraph2DErrors*)a->Get("Graph2D");
  
  TCanvas *c1 = new TCanvas("c1","c1",500,500);
  c1->cd(1);
    
  g->Draw("surf1");
  g->SetTitle("Data");
  g->GetXaxis()->SetTitle("X(deg)");
  g->GetXaxis()->SetTitleOffset(2.0);
  g->GetYaxis()->SetTitle("Y(deg)");
  g->GetYaxis()->SetTitleOffset(2.0);
  g->GetZaxis()->SetTitle("mu100");
  g->GetZaxis()->SetTitleOffset(1.2);  
  gPad->SetTheta(40.2113);
  gPad->SetPhi(-54.79899);

 
  //TF2 *f2 = new TF2("f2", "[0]+[1]*x+[2]*x*x+[3]*y+[4]*y*y+[5]*x*y",0.,90,0.,90.);

  TF2 *f2 = new TF2("f2",fitf,0.,90,0.,90.,6);
  f2->SetParameters(1.,1.,1.,1.,1.,1.);
  g->Fit("f2");
  
  double par2D[6];
  f2->GetParameters(par2D);
  cout << "Fit parameters = " << par2D[0] << ", " << par2D[1] << ", " << par2D[2] << "," 
	<< par2D[3] << ", " << par2D[4] << ", " << par2D[5] << endl;


  TF2 *fit2 = (TF2*)g->FindObject("f2");
  fit2->Draw("surf1");
  fit2->SetTitle("Minuit fit result on the Graph2D points");
  fit2->GetXaxis()->SetTitle("X");
  fit2->GetXaxis()->SetTitleOffset(2.0);
  fit2->GetYaxis()->SetTitle("Y");
  fit2->GetYaxis()->SetTitleOffset(2.0);
  fit2->GetZaxis()->SetTitle("Z");
  fit2->GetZaxis()->SetTitleOffset(1.2);
  g->Draw("same p0");
  
  gPad->SetTheta(40.2113);
  gPad->SetPhi(-54.79899);


}