// Simulazione fit lineare
#define NMISURE 10
void Regressione() {
    // Ipotiziamo che tra due variabili esiste una relazione
    // del tipo y=5x+3 che rappresenta un certo processo
    // fisico
    // assumiamo che solo le y siano affette da errori
    // di misura distributiti gaussianamente
    // Consideriamo il caso di 10 misure sperimentali
    // Vari casi affrontati.
    // Serve in questo caso a definire il package di minimizzazione
   TVirtualFitter::SetDefaultFitter("Minuit");
// generatore Random
   TRandom3 R(32525);
    double xTrue[NMISURE],x[NMISURE];
    double yTrue[NMISURE],y[NMISURE];
    double ex[NMISURE],  ey[NMISURE];
    double stepX=1.5;

    TH2D* hpar  = new TH2D("hpar","p1 vs p0",100, -3., 7., 100, 4.,6.);
    TH1D* hchi2 = new TH1D("hchi2","hchi2",100,0.,30.);
    TH1D* hpchi = new TH1D("hpchi","hpchi2",100,-0.1,1.1);
    
    // definisco una funzione parametrica
    TF1 * fun= new TF1("fun","[0]+[1]*x",0,20);
    TGraphErrors *gr;
    for (int iter=0; iter<5000; ++iter)
      {
	for (int i=0; i<NMISURE; i++) {
	  // questi sono i valori "veri" del processo fisico
	  xTrue[i]=i*stepX;
	  yTrue[i]=5*xTrue[i]+2;
	  // assume che gli errori di misura sulla y
	  // siano distribuiti gaussianamente con sigma=0.1
	  double eps=R.Gaus(0,2.);
	  // questo e' il valore misurato di Y.
	  y[i]=yTrue[i]+eps;
	  eps=R.Gaus(0,0.1);
	  // questo e' il valore misurato di X.
	  x[i]=xTrue[i]+eps;
	  // Errori sulle x nulli e sulle y sempre uguali a 1, 0.1, 0.01
	  ex[i]=0.1;
	  ey[i]=2.;
	}
	// costruisco il grafico
	gr= new TGraphErrors(NMISURE,x,y,ex,ey);
	// aspetto estetico.....
	gr->SetTitle("Y vs X");
	gr->SetMarkerColor(4);
	gr->SetMarkerStyle(21);
	gr->SetMarkerSize(0.2);
	// Disegno i punti "A" sta per disegna gli assi, "P"
	gr->Draw("AP");
	// La uso per fare un fit
	gr->Fit("fun");
	// estraggo i risultati della regressione
	TVirtualFitter * fitter = TVirtualFitter::GetFitter();
	// Estraggo e stampo il valore dei parametri.
	double p1=fitter->GetParameter(1);
	double p0=fitter->GetParameter(0);
	printf("\n********************** p1=%f p0=%f\n",p1,p0);
	// estraggo e stampo la matrice di covarianza
	double * covMatrix = fitter->GetCovarianceMatrix();
	printf("\n c11=%f   c12=%f\n",covMatrix[0],covMatrix[1]);
	printf(" c21=%f   c22=%f\n ",covMatrix[2],covMatrix[3]);
	printf("\n p0=%f +/- %f\n",p0,sqrt(covMatrix[0]));
	printf(" p1=%f +/- %f\n\n",p1,sqrt(covMatrix[3]));

	hpar->Fill(p0,p1);
	double chi2 = fun->GetChisquare();
	double pchi2 = fun->GetProb();
	hchi2->Fill( chi2 );
	hpchi->Fill( pchi2 );
      }
    TCanvas* c = new TCanvas("c","c",800,600);
    c->Divide(2,2);
    c->cd(1);
    //hchi2->Draw();
    gr->Draw("AP");
    c->cd(2);
    hpar->Draw("zcol");
    c->cd(3);
    TH1D* hpx = hpar->ProjectionX(); hpx->Draw();
    c->cd(4);
    TH1D* hpy = hpar->ProjectionY(); hpy->Draw();
    TCanvas* c1 = new TCanvas("c1","c1",1000,400);
    c1->Divide(2,1);
    c1->cd(1);
    hchi2->Draw();
    //gr->Draw("AP");
    c1->cd(2);
    hpchi->Draw();

}