#include "final_fitting.h"
#include "TF1.h"
#include "TH2.h"
#include "TMinuit.h"
#include "TProfile.h"

// a+bx 
double final_fitting::fit(TH2F   *spl,
                          double  low,
                          double  high, 
                          double &off) {   

 TF1 *fun = new TF1("fun","pol1",low,high);
 fun->SetLineColor(2);
 TProfile *prof = spl->ProfileX();
 prof->Fit("fun","rq");
 off = fun->GetParameter(0);

 double valor = fun->GetParameter(1);
 return valor;
}

// Fit to a first order polynomial, giving back the parameter, error in parameter and cov
double final_fitting::fitlin(TH2F   *spl,
                             double  low, 
                             double  high,
                             double &errpar,
                             double &offset, 
                             double &offset_error,
                             double &chis, 
                             double &cov00, 
                             double &cov01, 
                             double &cov11) {

 if (spl->GetEntries() > 25) {
   TF1      *fun  = new TF1("fun","pol1",low,high);
   TProfile *prof = spl->ProfileX();
   prof->Fit("fun","rq");
   offset       = fun->GetParameter(0);
   offset_error = fun->GetParError(0);
   errpar       = fun->GetParError(1);

   (fun->GetNDF() > 0) ? chis = fun->GetChisquare()/fun->GetNDF() : chis = 0.0;
    
   double emat[2][2]; // Error matrix
   gMinuit->mnemat(&emat[0][0],2);
   cov00 = emat[0][0];
   cov01 = emat[0][1];
   cov11 = emat[1][1];

   double valor = fun->GetParameter(1);
   return valor; 
 } else {
   std::cout << "Too few entries, " << spl->GetEntries() << std::endl;
   double valor = -999;
   return valor;
 }         
}

// a+bx+cx2+dx3
double final_fitting::fitlin_pol(TH2F   *spl, 
                                 double  low,
                                 double  high,
                                 double &par0, 
                                 double &par2, 
                                 double &par3, 
                                 double &epar0, 
                                 double &epar1, 
                                 double &epar2, 
                                 double &epar3) {
 
 TF1* fun = new TF1("fun","[3]*x*x*x+[2]*x*x+[1]*x+[0]",low,high);
 TProfile *prof = spl->ProfileX(); 
 prof->Fit("fun","rq"); 
 par0  = fun->GetParameter(0);
 par2  = fun->GetParameter(2);
 par3  = fun->GetParameter(3);     
 epar0 = fun->GetParError(0);
 epar1 = fun->GetParError(1);
 epar2 = fun->GetParError(2);
 epar3 = fun->GetParError(3);     

 double valor = fun->GetParameter(1);
 return valor;           
}

// Fit to a first order polynomial, giving back the parameter, error in parameter and cov
// Calculates by itself the range to fit. Used in the attenuation length
double final_fitting::fitlin2(TH2F *spl,
                              double &errpar,
                              double &offset,
                              double &offset_error, 
                              double &chis,
                              double &cov00,
                              double &cov01,
                              double &cov11) {

 double low  = -999.999;
 double high = -999.999;
 TProfile *prof = spl->ProfileX();

 for(int i=0; i<80; i++) {
   if (   prof->GetBinContent(i) < prof->GetBinContent(i+1) 
       && prof->GetBinError(i) >= prof->GetBinError(40) 
       && prof->GetBinError(40)*2 >= prof->GetBinError(i)) high = -2.0+0.05*(i+0.5);
 }
 
 for(int i=79; i>-1; i--) {
   if (   prof->GetBinContent(i) > prof->GetBinContent(i-1) 
       && prof->GetBinError(i) >= prof->GetBinError(40) 
       && prof->GetBinError(40)*2 >= prof->GetBinError(i)) low = -2.0+0.05*((i-1)+0.5);    
 }

 TF1 *fun = new TF1("fun","pol1",low,high); // Linear fit
 prof->Fit("fun","rq");
 errpar       = fun->GetParError(1);
 offset       = fun->GetParameter(0);
 offset_error = fun->GetParError(0);
 chis         = fun->GetChisquare();

 (fun->GetNDF() > 0) ? chis = fun->GetChisquare()/fun->GetNDF() : chis = 0.0;

 double emat[2][2]; // Error matrix
 gMinuit->mnemat(&emat[0][0],2);
 cov00 = emat[0][0];
 cov01 = emat[0][1];
 cov11 = emat[1][1];

 double valor = fun->GetParameter(1);
 return valor;
}

// Fit to a first order polynomial, giving back the parameter, error in parameter and cov
// Calculates by itself the range to fit. Used in the time walk
double final_fitting::fitlin3(TH2F   *spl,
                              double &errpar,
                              double &offset,
                              double &offset_error,  
                              double &chis,
                              double &cov00,
                              double &cov01, 
                              double &cov11) { 

 double low  = -999.999;
 double high = -999.999;
 TProfile *prof = spl->ProfileX();

 for(int i=2; i<25; i++) {
   //if(   prof->GetBinContent(i) < prof->GetBinContent(i+1) 
   //   && prof->GetBinError(i) >= prof->GetBinError(13)    
   //   && prof->GetBinError(13)*3 >= prof->GetBinError(i)) high = -0.05+0.004*(i+0.5);
   if (   prof->GetBinEntries(i) > 10
       && prof->GetBinEntries(i+1) > 10) {low = -0.05+0.004*(i+0.5);break;}

 }

 for(int i=23; i>0; i--) {
   //if(   prof->GetBinContent(i) > prof->GetBinContent(i-1) 
   //   && prof->GetBinError(i)    >= prof->GetBinError(13)
   //   && prof->GetBinError(13)*3 >= prof->GetBinError(i)) low = -0.05+0.004*((i-1)+0.5);

   if (   prof->GetBinEntries(i) > 10
       && prof->GetBinEntries(i-1) > 10) {high = -0.05+0.004*(i+0.5);break;}
 }

 TF1 *fun = new TF1("fun","[0]+[1]*x",low,high);
 prof->Fit("fun","rq");
 offset       = fun->GetParameter(0);
 offset_error = fun->GetParError(0);
 errpar       = fun->GetParError(1);

 (fun->GetNDF() > 0) ? chis = fun->GetChisquare()/fun->GetNDF() : chis = 0.0;

 double emat[2][2]; // Error matrix
 gMinuit->mnemat(&emat[0][0],2);
 cov00 = emat[0][0];
 cov01 = emat[0][1];
 cov11 = emat[1][1];
  
 double valor = fun->GetParameter(1);
 return valor;           
}

// Fit to a gaussian function, giving back the parameter, error in parameter
double final_fitting::gaus1(TH1F   *spl,
                            double  low,
                            double  high,
                            double &errvalor, 
                            double &sigma,
                            double &errsig,
                            double &chis) {
 
 TF1 *fun2 = new TF1("fun2","gaus",low,high);
 spl->Fit("fun2","rq");
 sigma    = fun2->GetParameter(2);
 errsig   = fun2->GetParError(2);
 errvalor = fun2->GetParError(1);

 (fun2->GetNDF() > 0) ? chis = fun2->GetChisquare()/fun2->GetNDF() : chis = 0.0;
 
 double valor = fun2->GetParameter(1);
 return valor;
}

// Fit to a+bx+c/(d-x)
double final_fitting::fit_p(TH2F   *spl,
                            double  low,
                            double  high, 
                            double &par0, 
                            double &par2, 
                            double &par3,
                            double &epar0,
                            double &epar1,
                            double &epar2,
                            double &epar3,
                            double &chis){   

 TF1 *fun1 = new TF1("fun1","[0]+[1]*x+[2]/([3]-x)",low,high);
 fun1->SetLineColor(2);
 TProfile *prof = spl->ProfileX();
 prof->Fit("fun1","rq");

 (fun1->GetNDF() > 0) ? chis = fun1->GetChisquare()/fun1->GetNDF() : chis = 0.0;
 
 par0  = fun1->GetParameter(0);
 par2  = fun1->GetParameter(2);
 par3  = fun1->GetParameter(3);
 epar0 = fun1->GetParError(0);
 epar1 = fun1->GetParError(1);
 epar2 = fun1->GetParError(2);
 epar3 = fun1->GetParError(3);

 double valor = fun1->GetParameter(1);
 return valor;
}

// Fit to a+bx+cx2+dx3
double final_fitting::mainfit(TGraph *graph, 
                              int     evt, 
                              double &errpar,
                              double &offset,
                              double &offset_error,
                              double &chis,
                              double &cov00,
                              double &cov01, 
                              double &cov11,
                              int     i,
                              int     j) {

 char title[50];

 TF1 *fun1 = new TF1("fun1","[0]+[1]*x+[2]*x*x+[3]*x*x*x",-0.04,0.04);   
 graph->Set(evt);  
 graph->Fit("fun1","rq");

 if(j == 1) sprintf(title,"graph_e_%d_%d",i,i+1);
 else sprintf(title,"graph_w_%d_%d",i,i+1);
 graph->Write(title);

 offset       = fun1->GetParameter(0);
 offset_error = fun1->GetParError(0);
 errpar       = fun1->GetParError(1);

 (fun1->GetNDF() > 0) ? chis = fun1->GetChisquare()/fun1->GetNDF() : chis = 0.0;

 double emat[4][4];    
 gMinuit->mnemat(&emat[0][0],4);
 cov00 = emat[0][0];
 cov01 = emat[0][1];
 cov11 = emat[1][1];

 double valor = fun1->GetParameter(1);
 return valor;
}

// Fit a linear graph
double final_fitting::linearGraph(TGraph *graph,
                                  double &errpar,
                                  double &offset,
                                  double &offset_error, 
	                          double &chis,
                                  double &cov00,
                                  double &cov01,
                                  double &cov11,
                                  int     i,
                                  int     side)
{
 char title[100];
 sprintf(title,"fun%d_%d_%d",side,i,i+1);
 TF1 *fun = new TF1("fun","pol1");
 graph->Fit(fun,"nqw");
 fun->Write(title);
 offset       = fun->GetParameter(0);
 offset_error = fun->GetParError(0);
 errpar       = fun->GetParError(1);

 (fun->GetNDF() > 0) ? chis = fun->GetChisquare()/fun->GetNDF() : chis = 0.0;

 double emat[2][2]; // Error matrix
 gMinuit->mnemat(&emat[0][0],2);
 cov00 = emat[0][0];
 cov01 = emat[0][1];
 cov11 = emat[1][1];
  
 double valor = fun->GetParameter(1);
 fun->Clear();
 return valor;         
}