void linear_fit_errors(void)
{
// Explore uncertainties coming from using correlated parameters from a linear fit

// #include <TMath.h>
#include <TRandom.h>

gROOT->Reset();
//TTree *Bfield = (TTree *) gROOT->FindObject("Bfield");
gStyle->SetPalette(1,0);
gStyle->SetOptStat(kTRUE);
gStyle->SetOptStat(111111);
// gStyle->SetOptFit(kFALSE);
gStyle->SetOptFit(kTRUE);
gStyle->SetOptFit(1111);
gStyle->SetPadRightMargin(0.15);
gStyle->SetPadLeftMargin(0.15);
gStyle->SetPadBottomMargin(0.15);
gStyle->SetFillColor(0);
//

   char string[256];
   Int_t j,jj;

   # define nevnts 1000;

   // set basis for "data"
 
   Int_t npts=10;
   // add dummy point on end to be able to plot beyond last data
   Float_t x[10]={1,2,3,4,5,6,7,8,9,10};
   Float_t y[10]; 
   Float_t ex[10]={0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01};
   Float_t ey[10]={0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5,0.5};
   Float_t yoffset=5;
   Float_t yslope=1;

   Float_t xmin = 0;
   Float_t xmax = 10;
   Float_t ymin = 0;
   Float_t ymax = 20;

  // define histograms to keep statistics

  Int_t nbins=100;
  TH1F *fita = new TH1F("fita","Fitted intercept",nbins,0,10);
  TH1F *fitb = new TH1F("fitb","Fitted slope",nbins,0,2);
  TH2D *deltay = new TH2D("deltay","Fitted function-true",nbins,0,10,nbins,-2,2);
  TH1F *estsiga = new TH1F("estsiga","Estimated siga",nbins,0,1);
  TH1F *estsigb = new TH1F("estsigb","Estimated sigb",nbins,0,0.1);

  // loop over events, generating data with random uncertainties

  for (jj=0;jj<nevnts;jj++) {
 
   // generate "data"

  TRandom1 *r = new TRandom1();
   for (j=0;j<npts;j++) {
        y[j] = yoffset + x[j]*yslope +  r->Gaus(0,ey[j]);
       }

   TGraphErrors *gr = new TGraphErrors (npts,x,y,ex,ey);

   // fit linear function to data

   TFitResultPtr ptr = gr->Fit("pol1","S","",xmin,xmax);
   TMatrixDSym cov = ptr->GetCovarianceMatrix(); // access to covariance matrix
   Float_t p0 = ptr->Value(0);
   Float_t p1 = ptr->Value(1);
   Double_t row1[2], row2[2];
   cov.ExtractRow(0,0,row1);
   cov.ExtractRow(1,0,row2);
   Float_t eps11 = row1[0];
   Float_t eps22 = row2[1];
   Float_t eps12 = row1[1];
   Float_t p0err = sqrt(eps11);
   Float_t p1err = sqrt(eps22);
   Double_t determinant = cov.Determinant();
   Float_t  chi2 = ptr->Chi2();
   // printf ("2 p0 %f p1 %f p0err %f p1err %f chi2 %f\n",p0,p1,p0err,p1err,chi2);
   // ptr->Print("V");
   // cov.Print();
   // printf ("matrix=%f %f %f %f, determinant=%f\n",row1[0],row1[1],row2[0],row2[1],determinant);
   
   fita->Fill(p0);
   fitb->Fill(p1);
   estsiga->Fill(p0err);
   estsigb->Fill(p1err);
   Int_t kmax=10;
   for (Int_t k=0;k<kmax+1;k++) {
          xval = k*(xmax-xmin)/kmax;
          Float_t yval = p0+p1*xval;
          Float_t ytrue = yoffset+yslope*xval;
          Float_t ydiff = yval - ytrue;
          deltay->Fill(xval,ydiff);
          // printf ("p0 %f p1 %f xval %f yoffset %f yslope %f yval %f ytrue %f ydiff %f\n",p0,p1,xval,yoffset, yslope,yval,ytrue,ydiff);
          }

   }

  // use last values of covariance matrix to plot estimate errors on y position

  TF1 *uncert = new TF1("uncert",uncert_func,0,10,3);
  uncert->SetLineColor(2);
  uncert->SetParameters(eps11,eps22,eps12);
  
   //
   TCanvas *c2 = new TCanvas("c2 linear fit errors","Plot linear function",200,10,700,700);
   // c2->SetGridx();
   // c2->SetGridy();
   c2->SetBorderMode(0);
   c2->SetFillColor(0);
   c2->Divide(2,2);
  
   c2->cd(1);
   /*c2_1->SetBorderMode(0);
   c2_1->SetFillColor(0);
   c2_1->SetGridx();
   c2_1->SetGridy();*/
   fita->Fit("gaus");
   fita->Draw();
  
   c2->cd(2);
   /*c2_2->SetBorderMode(0);
   c2_2->SetFillColor(0);
   c2_2->SetGridx();
   c2_2->SetGridy();*/
   fitb->Fit("gaus");
   fitb->Draw();
  
   c2->cd(3);
   /*c2_3->SetBorderMode(0);
   c2_3->SetFillColor(0);
   c2_3->SetGridx();
   c2_3->SetGridy();*/
   estsiga->Draw();
  
   c2->cd(4);
   /*c2_4->SetBorderMode(0);
   c2_4->SetFillColor(0);
   c2_4->SetGridx();
   c2_4->SetGridy();*/
   estsigb->Draw();


   TCanvas *c3 = new TCanvas("c3 linear fit errors","Plot linear function",200,10,700,700);
   // c3->SetGridx();
   // c3->SetGridy();
   c3->SetBorderMode(0);
   c3->SetFillColor(0);
   deltay->Draw("box");


   // Fit slices and plot fitted 
   deltay->FitSlicesY();
    
   TCanvas *c4 = new TCanvas("c4 linear fit errors","Plot linear function",200,10,700,700);
   // c4->SetGridx();
   // c4->SetGridy();
   c4->SetBorderMode(0);
   c4->SetFillColor(0);
   c4->Divide(2,2);
  
   c4->cd(1);
   /*c4_1->SetBorderMode(0);
   c4_1->SetFillColor(0);
   c4_1->SetGridx();
   c4_1->SetGridy();*/
   deltay->Draw("box");
  
   c4->cd(2);
   /*c4_2->SetBorderMode(0);
   c4_2->SetFillColor(0);
   c4_2->SetGridx();
   c4_2->SetGridy();*/
   deltay_0->Draw();
  
   c4->cd(3);
   /*c4_3->SetBorderMode(0);
   c4_3->SetFillColor(0);
   c4_3->SetGridx();
   c4_3->SetGridy();*/
   deltay_1->Draw();
  
   c4->cd(4);
   /*c4_4->SetBorderMode(0);
   c4_4->SetFillColor(0);
   c4_4->SetGridx();
   c4_4->SetGridy();*/
   deltay_2->Draw();
   uncert->Draw("Samep");

   //
   TCanvas *c1 = new TCanvas("c1 linear fit","Plot linear function",200,10,700,500);
   c1->SetGridx();
   c1->SetGridy();
   c1->SetBorderMode(0);
   c1->SetFillColor(0);

   gr->GetXaxis()->SetRangeUser(xmin,xmax);
   gr->GetYaxis()->SetRangeUser(ymin,ymax);
   gr->GetXaxis()->SetTitleSize(0.07);
   gr->GetYaxis()->SetTitleSize(0.07);
   gr->GetXaxis()->SetTitle("x_{i}");
   gr->GetYaxis()->SetTitle("y_{i}");
   gr->SetMarkerColor(2);
   gr->SetMarkerStyle(20);
   gr->SetTitle("");

   gr->Draw("AP");
 
   sprintf (string,"Fit function y(x) = a + bx\n");
   printf("string=%s",string);
   t1 = new TLatex(0.2,0.85,string);
   t1->SetNDC();
   t1->SetTextSize(0.05);
   t1->Draw();

   sprintf (string,"a = %.2f #pm%.2f\n",p0,p0err);
   printf("string=%s",string);
   t1 = new TLatex(0.2,0.78,string);
   t1->SetNDC();
   t1->SetTextSize(0.05);
   t1->Draw();
   sprintf (string,"b = %.2f #pm%.2f\n",p1,p1err);
   printf("string=%s",string);
   t1 = new TLatex(0.2,0.72,string);
   t1->SetNDC();
   t1->SetTextSize(0.05);
   t1->Draw();

   // c1->SaveAs("linear_fit_errors_c1.pdf");
   c1->SaveAs("linear_fit_errors_c1.png");
}
Double_t uncert_func (Double_t *x, Double_t *par) {

     // return uncertainty in y0 = x0*a + b,  value as computed using covariance matrix elements

     Double_t eps11=par[0];
     Double_t eps22=par[1];
     Double_t eps12=par[2];
     Double_t x1=x[0]; 

     char string[256];

     Double_t sigmay2=0;
    
     sigmay2 = 1*eps11 + x1*x1*eps22 + 2*x1*eps12;

     return sqrt(sigmay2);
}