#include <iostream>
#include <TGraphErrors.h>     
#include <TAxis.h>            
#include <TCanvas.h>          
#include <TF1.h>              
#include <iomanip>
#include <math.h>
#include <TBox.h>
#include <TPaveText.h>

using namespace std; 
void ist(){
 //gStyle->SetOptFit(1);
  float Bd1[] =  {0.364, 0.346, 0.322, 0.295,0.266, 0.234, 0.204, 0.171, 0.139, 0.106, 0.074, 0.04, 0.008, -0.024, -0.057,-0.09, -0.121, -0.154, -0.185, -0.218, -0.249, -0.28, -0.31, -0.339, -0.363};

  float sBd1[] =  {0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001};

  float id1[] = {1.80, 1.65, 1.50, 1.35, 1.20, 1.05, 0.90, 0.75, 0.60, 0.45, 0.30, 0.15, 0.00, -0.15, -0.30, -0.45, -0.60, -0.75, -0.90, -1.05, -1.20,-1.35,-1.50,-1.65,-1.80};

  float sid1[] = {0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01};

  float Bs1[] =  {-0.363, -0.345, -0.321, -0.294, -0.264, -0.233, -0.201, -0.17, -0.138, -0.105, -0.071, -0.035, -0.006, 0.025,	0.057,	0.09, 0.124, 0.155, 0.187, 0.22, 0.252,	0.281, 0.311, 0.339, 0.364};

  float sBs1[] =  {0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001};

  float is1[] = {-1.80,	-1.65, -1.50, -1.35, -1.20, -1.05, -0.90, -0.75, -0.60, -0.45, -0.30, -0.15, 0.00, 0.15, 0.30, 0.45, 0.60, 0.75, 0.90, 1.05, 1.20, 1.35, 1.50, 1.65, 1.80};

  float sis1[] = {0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01,0.01, 0.01};

  float Bd2[] =  {0.364, 0.337,	0.304, 0.264, 0.225, 0.182, 0.138, 0.089, 0.05,	0.007, -0.034, -0.078, -0.123, -0.167, -0.207, -0.252, -0.29, -0.328, -0.364};

  float sBd2[] =  {0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001};

  float id2[] = {1.80, 1.60, 1.40, 1.20, 1.00, 0.80, 0.60, 0.40, 0.20, 0.00, -0.20, -0.40, -0.60, -0.80, -1.00, -1.20, -1.40, -1.60, -1.80};

  float sid2[] = {0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01};
 
  float Bs2[] =  {-0.364, -0.337, -0.303, -0.263, -0.223, -0.181, -0.138, -0.095, -0.049, -0.007, 0.035, 0.08, 0.123, 0.167, 0.208, 0.251, 0.291, 0.329, 0.364};

  float sBs2[] =  {0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001, 0.001};

  float is2[] = {-1.80,	-1.60, -1.40, -1.20, -1.00, -0.80, -0.60, -0.40, -0.20,	0.00, 0.20, 0.40, 0.60,	0.80, 1.00, 1.20, 1.40,1.60, 1.80};

  float sis2[] = {0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01, 0.01};
 
  const int n1misure = 25;
  
  const int n2misure = 25; 

  const int n3misure = 19;
  
  const int n4misure = 19;

  
  TCanvas *cBid1 = new TCanvas("cBid1","prima discesa",10,10,600,400);
  TPaveText *r=new TPaveText(0.8,0.8,1,1,"brNDC");
  r->SetBorderSize(1);
  cBid1->SetFillColor(0);
  cBid1->cd();
  TGraphErrors *gBid1 =new TGraphErrors(n1misure,id1,Bd1,sid1,sBd1);
  gBid1->SetMarkerSize(0.6);
  gBid1->SetMarkerStyle(21);
  gBid1->SetTitle("prima discesa");
  gBid1->GetXaxis()->SetTitle("Intensità di corrente [A]");
  gBid1->GetYaxis()->SetTitle("Campo magnetico [T] ");
  gBid1->Draw("AP");

  
  
  TCanvas *cBis1 = new TCanvas("cBis1","prima salita",610,10,600,400);
  cBis1->SetFillColor(0);
  cBis1->cd();
  TPaveText *f=new TPaveText(0.8,0.8,1,1,"brNDC");
  f->SetBorderSize(1);
  TGraphErrors *gBis1 =new TGraphErrors(n2misure,is1,Bs1,sis1,sBs1);
  gBis1->SetMarkerSize(0.6);
  gBis1->SetMarkerStyle(21);
  gBis1->SetTitle("prima salita");
  gBis1->GetXaxis()->SetTitle("Intensità di corrente [A]");
  gBis1->GetYaxis()->SetTitle("Campo magnetico [T] ");
  gBis1->Draw("AP");

  
  TCanvas *cBid2 = new TCanvas("cBid2","seconda discesa",10,610,600,400);
  cBid2->SetFillColor(0);
  cBid2->cd();
  TPaveText *s=new TPaveText(0.8,0.8,1,1,"brNDC");
  s->SetBorderSize(1);
  TGraphErrors *gBid2 =new TGraphErrors(n3misure,id2,Bd2,sid2,sBd2);
  gBid2->SetMarkerSize(0.6);
  gBid2->SetMarkerStyle(21);
  gBid2->SetTitle("seconda discesa");
  gBid2->GetXaxis()->SetTitle("Intensità di corrente [A]");
  gBid2->GetYaxis()->SetTitle("Campo magnetico [T] ");
  gBid2->Draw("AP");
  
  
  TCanvas *cBis2 = new TCanvas("cBi2","seconda salita",610,610,600,400);
  cBis2->SetFillColor(0);
  cBis2->cd();
  TPaveText *d=new TPaveText(0.8,0.8,1,1,"brNDC");
  d->SetBorderSize(1);
  TGraphErrors *gBis2 =new TGraphErrors(n4misure,is2,Bs2,sis2,sBs2);
  gBis2->SetMarkerSize(0.6);
  gBis2->SetMarkerStyle(21);
  gBis2->SetTitle("seconda salita");
  gBis2->GetXaxis()->SetTitle("Intensità di corrente [A]");
  gBis2->GetYaxis()->SetTitle("Campo magnetico [T] ");
  gBis2->Draw("AP");
 
  cout << "\n\n Fit  \n " <<endl;

  TF1 *fit1 = new TF1("fit1","[0]*x+[1]",-1.4,1.4);
  fit1->SetParameter(0,0.2);
  fit1->SetParameter(1,1);
  fit1->SetParName(0,"Coefficiente angolare");
  fit1->SetParName(1,"Intercetta");
  fit1->SetLineColor(3);
 

  TF1 *fit2 = new TF1("fit2","[0]*x+[1]",-1.4,1.4);
  fit2->SetParameter(0,0.2);
  fit2->SetParameter(1,1);
  fit2->SetParName(0,"Coefficiente angolare");
  fit2->SetParName(1,"Intercetta");
  fit2->SetLineColor(4);
  

  TF1 *fit3 = new TF1("fit3","[0]*x+[1]",-1.3,1.3);
  fit3->SetParameter(0,0.2);
  fit3->SetParameter(1,1);
  fit3->SetParName(0,"Coefficiente angolare");
  fit3->SetParName(1,"Intercetta");
  fit3->SetLineColor(5);
  

  TF1 *fit4 = new TF1("fit4","[0]*x+[1]",-1.3,1.3);
  fit4->SetParameter(0,0.2);
  fit4->SetParameter(1,1);
  fit4->SetParName(0,"Coefficiente angolare");
  fit4->SetParName(1,"Intercetta");
  fit4->SetLineColor(6);
  
  
  gBis1->Fit(fit1,"R");
  double chi2,prob,df,a,erra,Int,errInt;
  chi2=fit1->GetChisquare();
  df=fit1->GetNDF() ;
  prob=fit1->GetProb();
  a=fit1->GetParameter(0);
  erra=fit1->GetParError(0);
  Int=fit1->GetParameter(1);
  errInt=fit1->GetParError(1);
  r->AddText(Form("#chi^{2} / DF = %3.1lf / %1.0lf ",chi2,df));
  r->AddText(Form("Prob = %3.4lf  ",prob));
  r->AddText(Form("m = %3.2lf #pm %3.2lf  ",a,erra));
  r->AddText(Form("b = %3.2lf #pm %3.2lf  ",Int,errInt));
 r->Draw();

  gBid1->Fit(fit2,"R");
  
  chi2=fit2->GetChisquare();
  df=fit2->GetNDF() ;
  prob=fit2->GetProb();
  f->AddText(Form("#chi^{2} / DF = %3.1lf / %1.0lf ",chi2,df));
  f->AddText(Form("Prob = %3.4lf  ",prob));
  f->Draw();

  gBis2->Fit(fit3,"R");
  
  chi2=fit3->GetChisquare();
  df=fit3->GetNDF() ;
  prob=fit3->GetProb();
  s->AddText(Form("#chi^{2} / DF = %3.1lf / %1.0lf ",chi2,df));
  s->AddText(Form("Prob = %3.4lf  ",prob));
  s->Draw();

  gBid2->Fit(fit4,"R");
 
  chi2=fit4->GetChisquare();
  df=fit4->GetProb();
  d->AddText(Form("#chi^{2} / DF = %3.1lf / %1.0lf ",chi2,df));
  d->AddText(Form("Prob = %3.4lf  ",prob));
  d->Draw();
 
  cout <<"\n\n Conti \n" << endl;

 /* double a = fit1->GetParameter(0);
  double b = fit1->GetParameter(1);
  double c = fit2->GetParameter(0);
  double d = fit2->GetParameter(1);
  double e = fit3->GetParameter(0);
  double f = fit3->GetParameter(1);
  double g = fit4->GetParameter(0);
  double h = fit4->GetParameter(1);*/
  double m1,errm1,m2,errm2,b1,errb1,b2,errb2;

  m1=((fit1->GetParameter(0))/pow((fit1->GetParError(0)),2)+(fit3->GetParameter(0))/pow(fit3->GetParError(0),2))/pow((1/pow((fit1->GetParError(0)),2)+(1/pow((fit3->GetParError(0)),2))),1);
  errm1=1/pow((1/pow((fit1->GetParError(0)),2)+(1/pow((fit3->GetParError(0)),2))),0.5);
  m2=((fit2->GetParameter(0))/pow((fit2->GetParError(0)),2)+(fit4->GetParameter(0))/pow(fit4->GetParError(0),2))/pow((1/pow((fit2->GetParError(0)),2)+(1/pow((fit4->GetParError(0)),2))),1);
  errm2=1/pow((1/pow((fit2->GetParError(0)),2)+(1/pow((fit4->GetParError(0)),2))),0.5);
  b1=(fit1->GetParameter(1)+fit3->GetParameter(1))/2;
  errb1=(fit1->GetParameter(1)-fit3->GetParameter(1)+fit1->GetParError(1)+fit3->GetParError(1))/2; //1/pow((1/pow((fit1->GetParError(1)),2)+(1/pow((fit3->GetParError(1)),2))),0.5);
  b2=(fit2->GetParameter(1)+fit4->GetParameter(1))/2;
  errb2=(fit2->GetParameter(1)-fit4->GetParameter(1)+fit2->GetParError(1)+fit4->GetParError(1))/2;

  cout << "\n Media coeff. angolare discese: " << m1 << "+-"<< errm1 <<endl; 
  cout << "\n Media coeff. angolare salite: " << m2 << "+-"<< errm2 <<endl;
  cout << "\n Media intercetta discese: " <<b1<< "+-"<< errb1 <<endl; 
  cout << "\n Media intercetta salite: " << b2 << "+-"<< errb2 << "\n" <<endl;

  double m,b,errm,errb;
  m=(m1/pow(errm1,2)+m2/pow(errm2,2))/(1/pow(errm1,2)+1/pow(errm2,2));
  errm=1/pow(1/pow(errm1,2)+1/pow(errm2,2),0.5);
  b=(b1+b2)/2;
  errb=(-b1+errb1+b2+errb2)/2;
  cout << "\n Media medie coeff. angolare: " << m << "+-" << errm <<endl;
  cout << "\n Media intercetta : " << b << "+-" << errb <<endl; 
  
}