//---------this program is to draw  the 2D graph of en ,angle,cs ,cserror
//-------and fit them with legendre 

#include <TMath.h>
#include <TGraph2DErrors.h>
#include <TStyle.h>
#include <TSystem.h>
#include "Math/SpecFunc.h" 
#include <TCanvas.h>
#include <TF2.h>
#include <TMultiGraph.h>
#include <iostream>
#include <TF1.h>
#include <fstream>

using std::cout; 
using std::endl; 

//----------first ,construct the function using to fit


//=====================================================================
//   the sub-function 0
Double_t fit_0(  Double_t *x, Double_t *par)
{
      Double_t fitval=par[0]*TMath::Exp(par[1]*x[0])+par[2]*TMath::Exp(par[3]*x[0])+(par[4]+par[5]*x[0])/(0.9+par[6]*(x[0]-par[7])*(x[0]-par[7]));
      return fitval;
 }


//====================================================================== 
//   the sub-function 1
Double_t  fit_1(  Double_t *x, Double_t *par) 
{
      Double_t fitval=(par[0]*TMath::Exp(par[1]*x[0])+par[2]*TMath::Exp(par[3]*x[0]))*ROOT::Math::assoc_legendre(1,0,cos(x[1]*TMath::Pi()/180));  
      return fitval;
}

//======================================================================
//   the sub-function 2
Double_t  fit_2(  Double_t *x, Double_t *par) 
{
      Double_t fitval=(par[0]*TMath::Exp(par[1]*x[0])+par[2]*TMath::Exp(par[3]*x[0]))*ROOT::Math::assoc_legendre(2,0,cos(x[1]*TMath::Pi()/180));  
      return fitval;
}

//======================================================================
//   the sub-function 3
Double_t  fit_3(  Double_t *x, Double_t *par) 
{
      Double_t fitval=(par[0]*TMath::Exp(par[1]*x[0])+par[2]*TMath::Exp(par[3]*x[0]))*ROOT::Math::assoc_legendre(3,0,cos(x[1]*TMath::Pi()/180));  
      return fitval;
}

//======================================================================
//   the sub-function 4
Double_t  fit_4(  Double_t *x, Double_t *par) 
{
      Double_t fitval=(par[0]*TMath::Exp(par[1]*x[0])+par[2]*TMath::Exp(par[3]*x[0]))*ROOT::Math::assoc_legendre(4,0,cos(x[1]*TMath::Pi()/180));  
                    
      return fitval;
}


//======================================================================
//sum of the fit function

Double_t fit_f( Double_t *x, Double_t *par)
{
    Double_t fitval=fit_3(x,&par[16])+fit_2(x,&par[12])+ fit_1(x,&par[8])+fit_0(x,par)+fit_4(x,&par[20]);
                  ;
               
    return fitval;
}



//==============================================================
//--------------------------------------------------------------

//            construct the main function 


//--------------------------------------------------------------
//==============================================================


void fit2Dg_of_diff()

{
   const Int_t counter=80;     // number of data point
	 Double_t en[counter],angle[counter],cs[counter],cserror[counter];
	 Double_t ex[counter]={0.0};
	 Double_t ey[counter]={0.0};
	 double zero;
	 int i;
	 
	 
	
	 gSystem->Load("libMathMore");  
	 gStyle->SetOptFit(1111);
	 
	  
	  
	 
	 
	 // open the data file and read them
	 ifstream alldata("c:\\root\\3he\\all_normalize_data\\all_diff_data.txt");
	 if(!alldata)
	 	  cout<<"No file"<<endl;
	 
	 for(i=0;i<counter;i++){
	 
	    alldata>>en[i]>>angle[i]>>cs[i]>>zero>>cserror[i];
	    
	    cout<<en[i]<<"  "<<angle[i]<<"  "<<cs[i]<<"  "<<ex[i]<<" "<<cserror[i]<<endl;
	  }


//---------------------------------------
      
    TCanvas *c1 = new TCanvas("c1", "Fitting a 2DG ");  
    
      
    
	  TGraph2DErrors *f=new TGraph2DErrors(counter,en,angle,cs,ex,ey,cserror);
	    
	  f->Draw("surf1");
	  
	  
	  
	  
	  
	  
	  
	  
	  
	  //fit
	  //TFile ff("2d.root","new");
	  
	 TF2 *f2=new TF2("fit",fit_f,4,150,0,180,24);
	 
	 
	 for(i=0;i<24;i++) {
	    f2->SetParameter(i,1.);
	    f2->SetParLimits(i,-10000,10000);
         }
	   
	  
	 f->Fit("fit","EV");
	 //f->Fit("fit","E");
	 
	
	 
	 Double_t par[24];
	 Double_t err[24];
	 ofstream outpar("c:\\root\\3he\\final par.dat");
	  if(!outpar)
	  	cout<<"no file"<<endl;
	 
	 f2->GetParameters(par);
	 Double_t *x=f2->GetParErrors();
	 for(i=0;i<24;i++)
	 	  outpar<<i<<"    "<< par[i]<<"      "<<*(x+i)<<endl;
	 	  
	    
	  
	  
	}