// Esempio di analisi dei residui
//
#define NMISURE 10
void Residui1(bool verbose=false) {

    gStyle->SetOptFit(111111);
  
    // Ipotiziamo che tra due variabili esiste una relazione
    // del tipo y=5x+2 che rappresenta un certo processo
    TRandom R(14567);
    double xTrue[NMISURE],x[NMISURE];
    double yTrue[NMISURE],y[NMISURE];
    double ex[NMISURE],ey[NMISURE];
    double stepX=1.5;
    double errY=0.1;
    double xm1[NMISURE];
    double ym1[NMISURE];
    double exm1[NMISURE],eym1[NMISURE];


    TH1D* hpulls = new TH1D("hPulls","Pulls",100.,-5.,5.);
    TH1D* hpullst = new TH1D("hPullst","Pulls (truth)",100.,-5.,5.);
    TH1D* hpullsC = new TH1D("hPullsC","Pulls (corrected)",100.,-5.,5.);
    TH1D* hpullsCC = new TH1D("hPullsCC","Pulls (corrected with err)",100.,-5.,5.);

    int nIter= 10000;
    for (int iter=0; iter<nIter; ++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 gaussiani
	  // con sigma =0.1
	  double eps=R.Gaus(0,errY);
	  y[i]=yTrue[i]+eps;
	  x[i]=xTrue[i];
	  // Stimo correttamente l'errore per la regressione
	  ex[i]=0.;
	  ey[i]=errY*1.;

	  //
	  xm1[i]=x[i];
	  ym1[i]=y[i];
	}
	
	// costruisco il grafico
	TGraphErrors *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");
	// definisco una funzione parametrica
	TF1 * fun= new TF1("fun","[0]+[1]*x",0,20);
	// 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);
 	// estraggo e stampo la matrice di covarianza
	double * covMatrix = fitter->GetCovarianceMatrix();
	if (verbose){
	  printf("\n p1=%f p0=%f\n",p1,p0);
	  printf("\n c11=%f   c12=%f ",covMatrix[0],covMatrix[1]);
	  printf("\n c21=%f   c22=%f\n ",covMatrix[2],covMatrix[3]);	
	  printf("\n p0+/-errp0 = %f +/- %f\n",p0,sqrt(covMatrix[0]));
	  printf("\n p1+/-errp1 = %f +/- %f\n",p1,sqrt(covMatrix[3]));
	}
	// Calcolo i residui con relativo errore tenendo conto della
	// covarianza
	//  residuo = y-f(x)
	double residuoC=0.;
	double ErrResiduoC=0.;
	double residuo[NMISURE];
	double residuot[NMISURE];
	double ErrResiduo[NMISURE]={0};
	double tmpx=0; double tmpy=0.;
	for (int i=0; i<NMISURE; i++){
	  residuo[i]=y[i]-p1*x[i]-p0;
	  residuot[i]=y[i]-5.*x[i]-2.;
	  ErrResiduo[i]=ey[i];
	  hpulls->Fill(residuo[i]/ErrResiduo[i]);
	  hpullst->Fill(residuot[i]/ErrResiduo[i]);
	  //residuo[i] = residuo[i]/ErrResiduo[i];

	  tmpx=xm1[i];
	  tmpy=ym1[i];
	  xm1[i]=xm1[NMISURE-1];
	  ym1[i]=ym1[NMISURE-1];
	  for (int j=0;j<NMISURE; ++j) std::cout<<" ("<<  x[j]<<", "<<  y[j]<<") ";
	  std::cout<<std::endl;
	  for (int j=0;j<NMISURE; ++j) std::cout<<" ("<<xm1[j]<<", "<<ym1[j]<<") ";
	  std::cout<<std::endl;
	  TGraphErrors *grm1= new TGraphErrors(NMISURE-1,xm1,ym1,ex,ey);
	  grm1->Fit("fun","q");
	  // estraggo i risultati della regressione
	  fitter = TVirtualFitter::GetFitter();
	  // Estraggo e stampo il valore dei parametri.
	  double p1=fitter->GetParameter(1);
	  double p0=fitter->GetParameter(0);
	  covMatrix = fitter->GetCovarianceMatrix();
	  
	  residuoC=tmpy-p1*tmpx-p0;
	  ErrResiduoC = sqrt( ey[i]*ey[i] + covMatrix[3]*tmpx*tmpx + covMatrix[0] + 2.*covMatrix[1]*tmpx);
	  hpullsC->Fill(residuoC/ErrResiduo[i]);
	  hpullsCC->Fill(residuoC/ErrResiduoC);
	  xm1[i]=tmpx;
	  ym1[i]=tmpy;
	  
	}


	if (iter==nIter-1)
	  {
	    // Grafico residui
	    TGraphErrors *gr2= new TGraphErrors(NMISURE,x,residuo,ex,ErrResiduo);
	    TGraphErrors *gr2t= new TGraphErrors(NMISURE,x,residuot,ex,ErrResiduo);
	    gr2->SetTitle("Residui vs X");
	    gr2->SetMarkerColor(4);
	    gr2->SetMarkerStyle(20);
	    gr2->SetMarkerSize(0.8);
	    gr2->Draw("AP");
	    gr2t->SetTitle("Residui vs X");
	    gr2t->SetMarkerColor(2);
	    gr2t->SetMarkerStyle(21);
	    gr2t->SetMarkerSize(0.8);
	    gr2t->Draw("P");
	  }
      }
    
    TCanvas * c = new TCanvas("c","c",800,600);
    hpulls->Fit("gaus");
    
    TVirtualFitter * fitter = TVirtualFitter::GetFitter();
    double p2=fitter->GetParameter(2);
    double p1=fitter->GetParameter(1);
    double p0=fitter->GetParameter(0);
    
    double ep2=fitter->GetParError(2);
    double ep1=fitter->GetParError(1);
    double ep0=fitter->GetParError(0);

    printf("\n Fit of the pulls: p0 (normalization) = %f +/- %f ",p0,ep0);
    printf("\n Fit of the pulls: p1 (mean)          = %f +/- %f ",p1,ep1);
    printf("\n Fit of the pulls: p2 (sigma)         = %f +/- %f\n\n ",p2,ep2);


    hpullst->Fit("gaus");

    c->Clear();
    c->Divide(2,2);
    c->cd(1);
    hpulls->Draw();
    c->cd(2);
    hpullst->Draw();
    c->cd(3);
    

    fitter = TVirtualFitter::GetFitter();
    p2=fitter->GetParameter(2);
    p1=fitter->GetParameter(1);
    p0=fitter->GetParameter(0);
    
    ep2=fitter->GetParError(2);
    ep1=fitter->GetParError(1);
    ep0=fitter->GetParError(0);
    
    printf("\n Fit of the pulls (teo): p0 (normalization) = %f +/- %f ",p0,ep0);
    printf("\n Fit of the pulls (teo): p1 (mean)          = %f +/- %f ",p1,ep1);
    printf("\n Fit of the pulls (teo): p2 (sigma)         = %f +/- %f\n\n ",p2,ep2);


    hpullsC->Fit("gaus");
    fitter = TVirtualFitter::GetFitter();
    p2=fitter->GetParameter(2);
    p1=fitter->GetParameter(1);
    p0=fitter->GetParameter(0);
    
    ep2=fitter->GetParError(2);
    ep1=fitter->GetParError(1);
    ep0=fitter->GetParError(0);
    
    printf("\n Fit of the pulls (teo): p0 (normalization) = %f +/- %f ",p0,ep0);
    printf("\n Fit of the pulls (teo): p1 (mean)          = %f +/- %f ",p1,ep1);
    printf("\n Fit of the pulls (teo): p2 (sigma)         = %f +/- %f\n\n ",p2,ep2);
    hpullsC->Draw();


    c->cd(4);
    hpullsCC->Fit("gaus");
    fitter = TVirtualFitter::GetFitter();
    p2=fitter->GetParameter(2);
    p1=fitter->GetParameter(1);
    p0=fitter->GetParameter(0);
    
    ep2=fitter->GetParError(2);
    ep1=fitter->GetParError(1);
    ep0=fitter->GetParError(0);
    
    printf("\n Fit of the pulls (teo): p0 (normalization) = %f +/- %f ",p0,ep0);
    printf("\n Fit of the pulls (teo): p1 (mean)          = %f +/- %f ",p1,ep1);
    printf("\n Fit of the pulls (teo): p2 (sigma)         = %f +/- %f\n\n ",p2,ep2);
    hpullsCC->Draw();

    
}