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

using namespace std; 

void rha()// solo fit delle rette dai dati originali

{
	double vh[] = {0.03, 0.04, 0.04, 0.03, 0.03, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.02, 0.01, 0.01, 0.01, 0.01, 0.03, 0.03, 0.03, 0.03, 0.03, 0.03, 0.03, 0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 0.05, 0.05}; //RIFERIMENTO

	double vhB[] = {3.55,6.88,8.63,-3.44,-6.83,-8.64};

	double svhB[] = {0.02,0.02,0.02,0.02,0.02,0.02};

	double vhB2[] = {-3.47,-6.80,-8.65,3.47,6.89,8.62};

	double B[] = {0.1081,0.2157,0.2802,-0.1071,-0.2147,-0.2793};

	double sB[] = {0.0086,0.0086,0.0086,0.0086,0.0086,0.0086};

	for (int j = 0; j<6 ; ++j) 
	{
		vhB[j]=vhB[j]-0.04;
	}

	TCanvas *crh = new TCanvas("crh","Rh",10,10,600,400);
	crh->SetFillColor(0);
	crh->cd();

	TGraphErrors *VB1 =new TGraphErrors(6,B,vhB,sB,svhB);
	VB1->Draw("AP");
	TGraphErrors *VB2 =new TGraphErrors(6,B,vhB2,sB,svhB);
	VB2->SetTitle("VH(B)");
	VB2->GetXaxis()->SetTitle("B[T]");
	VB2->GetYaxis()->SetTitle("VH [mV]");
	VB2->Draw("AP");

	  VB1->Draw("SAME");
	//VB2->Draw("SAME");

	TF1 *fit1 = new TF1("fit1","[0]*x"); //B<0 (0,50A)
	fit1->SetParameter(0,+50);
    fit1->SetParameter(1,0.1);
    fit1->SetParName(0,"Coefficiente angolare");
    fit1->SetParName(1,"Intercetta");
    fit1->SetLineColor(1);
    VB1->Fit(fit1,"M");

    TF1 *fit2 = new TF1("fit2","[0]*x"); //B<0 (0,50A)
    fit2->SetParameter(0,-2);
    fit2->SetParameter(1,0.1);
    fit2->SetParName(0,"Coefficiente angolare");
    fit2->SetParName(1,"Intercetta");
    fit2->SetLineColor(2);
    VB2->Fit(fit2,"M");

	cout << " Chi quadro 5a: " << fit1->GetChisquare() << ", number of DoF: " << fit1->GetNDF() << " (Probability: " << fit1->GetProb() << ")." << endl;
    cout << " Chi quadro 5b: " << fit2->GetChisquare() << ", number of DoF: " << fit2->GetNDF() << " (Probability: " << fit2->GetProb() << ")." << endl;

	double Rh1,errRh1,Rh2,errRh2;
	Rh1=(fit1->GetParameter(0)*pow(10,-3))*pow(10,-3)/(6.00*pow(10,-3));
	errRh1=pow(pow(fit1->GetParError(0)*pow(10,-3)*pow(10,-3)/(6.00*pow(10,-3)),2)+pow(0.08*pow(10,-3)*(fit1->GetParameter(0)*pow(10,-3)*pow(10,-3)/pow(6.00*pow(10,-3),2)),2),0.5);
	Rh2=-(fit2->GetParameter(0)*pow(10,-3))*pow(10,-3)/(6.00*pow(10,-3));
	errRh2=pow(pow(fit2->GetParError(0)*pow(10,-3)*pow(10,-3)/(6.00*pow(10,-3)),2)+pow(0.08*pow(10,-3)*(fit2->GetParameter(0)*pow(10,-3)*pow(10,-3)/pow(6.00*pow(10,-3),2)),2),0.5);

	cout << "\n Media coefficiente di Hall I>0: " << Rh1 << " +- " << errRh1 << " m^3/C" << endl;

	cout << " Media coefficiente di Hall I<0: " << Rh2 << " +- " << errRh2 << " m^3/C" << endl;

	cout << " Media coefficiente di Hall : " << (Rh1/pow(errRh1,2)+Rh2/pow(errRh2,2))/(1/pow(errRh1,2)+1/pow(errRh2,2)) << " +- " << 1/pow(1/pow(errRh1,2)+1/pow(errRh2,2),0.5) << " m^3/C" << endl;
}