Double_t fitfunc_func (Double_t *x, Double_t *par) {
    
    // function to be plotted: Polynomial for background + Gaussian signal.
    
    Double_t amp=par[0];
    Double_t mean=par[1];
    Double_t sigma=par[2];
    Double_t p1=par[3];
    Double_t x1=x[0];
    
    char string[256];
    Double_t pi=3.14159;
    
    Double_t arg = (x1-mean)*(x1-mean)/(2*sigma*sigma);
    // Double_t signal = amp*exp(-arg)/sqrt(2*pi*sigma*sigma);
    Double_t signal = amp*exp(-arg);
    
    Double_t back = p1;   // constant background
    
    sprintf(string,"p1=%g, x1=%g, amp=%g, mean=%g, sigma=%g, fitfunc=%g\n",p1,x1,amp,mean,sigma,back+signal);
    // printf ("func: %s",string);
    
    return back+signal;
}
void Gaussian_Area_Errors(void)
{
// 
// Study correllations between fitting errors when using a Gaussian fit to determine areas
//

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

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

   char string[256];
   char filename[80];
   // obtain number of events assuming uniform weights, use "W"=W, statistical weights use "W"=" ".
   // char W[1]="W";
   // char W[1]=" ";
   TString W=" ";
   Int_t j,jj;
   
   const Int_t npts1=1000;
   const Int_t npts2=8000;

   Int_t nsigma=5;           // signal region considered +/- nsigma
   Double_t mean1=0;   // Gaussian mean
   Double_t sigma1=1;  // Gaussian sigma

   TRandom1 *r1 = new TRandom1();
   TRandom1 *r2 = new TRandom1();
   TRandom1 *r3 = new TRandom1();
   TRandom1 *r4 = new TRandom1();

   // initialize random number seeds

   UInt_t  iseed = 23422777517;
   r1->SetSeed(iseed);
   iseed = 23424111595;
   r2->SetSeed(iseed);
   r3->SetSeed(iseed);
   r4->SetSeed(iseed);

   // define histograms to collect statistics

   Int_t nexpts=100;
   Int_t exp_nbins=100;
   Float_t exp_xmin=0;
   Float_t exp_xmax=npts1*2;
   Float_t exp_ymin=0;
   Float_t exp_ymax=nexpts;
   Float_t back_xmin=2000;
   Float_t back_xmax=5000;

   
   TH1D  *exp_signal = new TH1D("exp_signal","Signal",exp_nbins,exp_xmin,exp_xmax);
   TH1D  *exp_gauss = new TH1D("exp_gauss","Fitted Gauss",exp_nbins,exp_xmin,exp_xmax);
   TH1D  *exp_signal_err = new TH1D("exp_signal_err","Sqrt(Signal)",exp_nbins,0,200);
   TH1D  *exp_gauss_err = new TH1D("exp_gauss_err","Area Correct Error",exp_nbins,0,200);
   TH1D  *exp_gauss_err1 = new TH1D("exp_gauss_err1","Area (DeltaAmp/Amp)",exp_nbins,0,200);
   TH1D  *exp_gauss_err2 = new TH1D("exp_gauss_err2","Area (DeltaAmp/Amp)^2 + (DeltaSigma/Sigma)^2",exp_nbins,0,200);
 
   TH2D  *gaus_parms = new TH2D("gaus_parms","Amplitude vs Sigma",exp_nbins,0.5,1.5,exp_nbins,0,200);
   TH2D  *gaus_errors = new TH2D("gaus_errors","DeltaAmp/Amp vs DeltaSigma/Sigma",exp_nbins,0,0.1,exp_nbins,0,0.05);

   Int_t nbins=100;
   Float_t xmin=-10*sigma1;
   Float_t xmax=+10*sigma1;
   Float_t ymin=0;
   Float_t ymax=300;
    
    

   TF1 *fitfunc = new TF1("fitfunc",fitfunc_func,xmin,xmax,4);
   TF1 *fitback = new TF1("fitback","[0]",xmin,xmax);

    TCanvas *c1 = new TCanvas("c1","c1 Summary of Exptl Results",200,10,700,700);
    
    
   // loop over 'expts' to collect statistics 

   for (jj=0; jj<nexpts;jj++){

   TH1D  *hsignal = new TH1D("hsignal","Gen Signal",nbins,xmin,xmax);

    // For each experiment, generate the number of  signal and background events according to gaussian distributions.

    Double_t signal_counts = r3->Gaus(npts1,sqrt(npts1));
    Int_t npts_signal = signal_counts;
    printf ("\n\nExperiment %d, npts_signal=%d\n",jj,npts_signal);

    // generate data

   for (j=0;j<npts_signal;j++) {

       // generate Signal with Gaussian (mean1,sigma1);  npts_signal events

       Double_t z1 = r1->Gaus(mean1,sigma1);
       hsignal->Fill(z1);
       }

   Double_t amp =50;
   Double_t mean=mean1;
   Double_t sigma=sigma1;
   Double_t p1=0;
   Double_t p1_err=0;

   fitfunc->SetParameters(amp,mean,sigma,p1);
   fitfunc->FixParameter(3,p1);

   // get uncertainties using statistical weights

   // hsignal->Fit("fitfunc");
       
       
   TFitResultPtr ptr = hsignal->Fit("fitfunc","S","",xmin,xmax);
   TMatrixDSym cov = ptr->GetCovarianceMatrix(); // access to covariance matrix
   amp = ptr->Value(0);
   mean = ptr->Value(1);
   sigma = ptr->Value(2);
   Double_t row1[3], row2[3],row3[3];
   cov.ExtractRow(0,0,row1);
   cov.ExtractRow(1,0,row2);
   cov.ExtractRow(2,0,row3);
   Float_t eps11 = row1[0];
   Float_t eps22 = row2[1];
   Float_t eps33 = row3[2];
   Float_t eps12 = row1[1];
   Float_t eps13 = row1[2];
   Float_t eps23 = row2[2];
   Float_t amp_err = sqrt(eps11);
   Float_t mean_err = sqrt(eps22);
   Float_t sigma_err = sqrt(eps33);
   Double_t determinant = cov.Determinant();
   Float_t  chi2 = ptr->Chi2();
   printf ("amp=%f mean=%f sigma=%f\n",amp,mean,sigma);
   printf ("amp_err=%f mean_err=%f sigma_err=%f chi2=%f\n",amp_err,mean_err,sigma_err,chi2);
   ptr->Print("V");
   cov.Print();
   // printf ("matrix=%f %f %f %f, determinant=%f\n",row1[0],row1[1],row2[0],row2[1],determinant);


   // obtain number of events assuming uniform weights, use "W"=W, statistical weights use "W"=" ".

   // hsignal->Fit("fitfunc","W");
   hsignal->Fit("fitfunc",W);

   // get number of signal events

   amp = fitfunc->GetParameter(0);
   mean = fitfunc->GetParameter(1);
   sigma = fitfunc->GetParameter(2);
   fitfunc->SetParNames("Constant","Mean","Sigma");
   printf ("p1=%g, p1_err=%g, amp=%g, mean=%g, sigma=%g, amp_err=%f, mean_err=%f, sigma_err=%f\n",p1,p1_err, amp,mean,sigma,amp_err,mean_err,sigma_err);
   gaus_errors->Fill(sigma_err/sigma,amp_err/amp);
   gaus_parms->Fill(sigma,amp);

   Double_t bin = (xmax-xmin)/nbins;
   Double_t xlo = mean1 - nsigma*sigma1;
   Double_t xhi = mean1 + nsigma*sigma1;
   Int_t jlo = (xlo-xmin)/bin;
   Int_t jhi = (xhi-xmin)/bin;

   Double_t signal_est_err = sqrt(npts_signal);
   Double_t signal_gaus = fitfunc->Integral(xlo,xhi)/bin ;

   Double_t pi=3.14159;

   Double_t signal_gaus_err = sqrt(2*pi*amp*amp*eps33 + 2*pi*sigma*sigma*eps11 + 2*2*pi*sigma*amp*eps13)/bin;
   // Double_t signal_gaus_err = sqrt(eps11)/bin;
   printf ("amp=%g, sigma=%g, eps11=%g, eps33=%g, eps13=%g\n",amp,sigma,eps11,eps33,eps13);
   Double_t signal_gaus_err1 = signal_gaus*sqrt((amp_err/amp)*(amp_err/amp));  // ignore sigma_err as it is highly correlated with amp_err
   Double_t signal_gaus_err2 = signal_gaus*sqrt((amp_err/amp)*(amp_err/amp) + (sigma_err/sigma)*(sigma_err/sigma));  // ignore correllations
   printf ("signal_counts = %g, signal_est_err=%g, signal_gaus=%g, signal_gaus_err=%g \n",signal_counts,signal_est_err,signal_gaus,signal_gaus_err);

  // fill summary histogram with 'exptl' results

   exp_signal->Fill(signal_counts);
   exp_gauss->Fill(signal_gaus);
   exp_signal_err ->Fill(signal_est_err);
   exp_gauss_err->Fill(signal_gaus_err);
   exp_gauss_err1->Fill(signal_gaus_err1);
   exp_gauss_err2->Fill(signal_gaus_err2);


   // delete histograms and repeat to build up statistics

   /*hsignal->Delete();
   hback->Delete();
   hback_pure->Delete();
   hdata->Delete();*/

   // end of loop over number of experiments

   }

   sprintf (string,"Err Fit Option='%s'\n",W.Data());
   printf("string=%s",string);
   TLatex *t1 = new TLatex(0.16,0.92,string);
   t1->SetNDC();
   t1->SetTextSize(0.05);
   t1->Draw();

   
   TCanvas *c3 = new TCanvas("c3","c3 Summary of Exptl Results",200,10,700,700);
   c3->Divide(2,2);
   // c3->SetGridx();
   // c3->SetGridy();
   // c3->SetLogy();

   c3->cd(1);
   exp_signal->GetXaxis()->SetRangeUser(exp_xmin,exp_xmax);
   // exp_signal->GetYaxis()->SetRangeUser(exp_ymin,exp_ymax);
   exp_signal->GetXaxis()->SetTitleSize(0.07);
   exp_signal->GetYaxis()->SetTitleSize(0.07);
   exp_signal->GetXaxis()->SetTitle("Signal");
   exp_signal->GetYaxis()->SetTitle("events");
   exp_signal->GetXaxis()->SetNdivisions(505);
   exp_signal->SetMarkerColor(2);
   exp_signal->SetMarkerSize(0.5);
   exp_signal->SetLineColor(2);
   exp_signal->SetFillColor(2);
   exp_signal->SetFillStyle(11);
   exp_signal->SetMarkerStyle(20);
   exp_signal->SetTitle("");
   exp_signal->Fit("gaus");
   exp_signal->Draw();
   sprintf (string,"Err Fit Option='%s'\n",W.Data());
   printf("string=%s",string);
   t1->DrawLatex(0.16,0.92,string);
   t1->SetNDC();
   t1->SetTextSize(0.05);
   t1->Draw();

   c3->cd(2);
   exp_gauss->GetXaxis()->SetRangeUser(exp_xmin,exp_xmax);
   // exp_gauss->GetYaxis()->SetRangeUser(exp_ymin,exp_ymax);
   exp_gauss->GetXaxis()->SetTitleSize(0.07);
   exp_gauss->GetYaxis()->SetTitleSize(0.07);
   exp_gauss->GetXaxis()->SetTitle("Fitted Gaussian");
   exp_gauss->GetYaxis()->SetTitle("events");
   exp_gauss->GetXaxis()->SetNdivisions(505);
   exp_gauss->SetMarkerColor(2);
   exp_gauss->SetMarkerSize(0.5);
   exp_gauss->SetLineColor(2);
   exp_gauss->SetFillColor(2);
   exp_gauss->SetFillStyle(11);
   exp_gauss->SetMarkerStyle(20);
   exp_gauss->SetTitle("");
   exp_gauss->Fit("gaus");
   exp_gauss->Draw();
   sprintf (string,"Err Fit Option='%s'\n",W.Data());
   printf("string=%s",string);
   t1->DrawLatex(0.16,0.92,string);
   t1->SetNDC();
   t1->SetTextSize(0.05);
   t1->Draw();


   c3->cd(3);
   exp_signal_err->GetXaxis()->SetRangeUser(0,200);
   // exp_signal_err->GetYaxis()->SetRangeUser(exp_ymin,exp_ymax);
   exp_signal_err->GetXaxis()->SetTitleSize(0.07);
   exp_signal_err->GetYaxis()->SetTitleSize(0.07);
   exp_signal_err->GetXaxis()->SetTitle("Sqrt(Signal)");
   exp_signal_err->GetYaxis()->SetTitle("events");
   exp_signal_err->GetXaxis()->SetNdivisions(505);
   exp_signal_err->SetMarkerColor(2);
   exp_signal_err->SetMarkerSize(0.5);
   exp_signal_err->SetLineColor(2);
   exp_signal_err->SetFillColor(2);
   exp_signal_err->SetFillStyle(11);
   exp_signal_err->SetMarkerStyle(20);
   exp_signal_err->SetTitle("");
   exp_signal_err->Draw();
   sprintf (string,"Err Fit Option='%s'\n",W.Data());
   printf("string=%s",string);
   t1->DrawLatex(0.16,0.92,string);
   t1->SetNDC();
   t1->SetTextSize(0.05);
   t1->Draw();

   c3->cd(4);
   exp_gauss_err->GetXaxis()->SetRangeUser(0,200);
   // exp_gauss_err->GetYaxis()->SetRangeUser(exp_ymin,exp_ymax);
   exp_gauss_err->GetXaxis()->SetTitleSize(0.07);
   exp_gauss_err->GetYaxis()->SetTitleSize(0.07);
   exp_gauss_err->GetXaxis()->SetTitle("Full Fitted Gauss Err");
   exp_gauss_err->GetYaxis()->SetTitle("events");
   exp_gauss_err->GetXaxis()->SetNdivisions(505);
   exp_gauss_err->SetMarkerColor(2);
   exp_gauss_err->SetMarkerSize(0.5);
   exp_gauss_err->SetLineColor(2);
   exp_gauss_err->SetFillColor(2);
   exp_gauss_err->SetFillStyle(11);
   exp_gauss_err->SetMarkerStyle(20);
   exp_gauss_err->SetTitle("");
   exp_gauss_err->Draw();
   sprintf (string,"Err Fit Option='%s'\n",W.Data());
   printf("string=%s",string);
   t1->DrawLatex(0.16,0.92,string);
   t1->SetNDC();
   t1->SetTextSize(0.05);
   t1->Draw();

 
   TCanvas *c4 = new TCanvas("c4","c4 Summary of Exptl Results",200,10,700,700);
   c4->Divide(2,2);
   // c4->SetGridx();
   // c4->SetGridy();
   // c4->SetLogy();

   c4->cd(1);
   gaus_parms->GetXaxis()->SetRangeUser(exp_xmin,exp_xmax);
   // gaus_parms->GetYaxis()->SetRangeUser(exp_ymin,exp_ymax);
   gaus_parms->GetXaxis()->SetTitleSize(0.05);
   gaus_parms->GetYaxis()->SetTitleSize(0.05);
   gaus_parms->GetXaxis()->SetTitle("Sigma");
   gaus_parms->GetYaxis()->SetTitle("Amplitude");
   gaus_parms->GetYaxis()->SetTitleOffset(1.5);
   gaus_parms->GetXaxis()->SetNdivisions(505);
   gaus_parms->SetMarkerColor(2);
   gaus_parms->SetMarkerSize(0.5);
   gaus_parms->SetLineColor(2);
   gaus_parms->SetFillColor(2);
   gaus_parms->SetFillStyle(11);
   gaus_parms->SetMarkerStyle(20);
   gaus_parms->SetTitle("");
   gaus_parms->Draw();

   c4->cd(2);
   gaus_errors->GetXaxis()->SetRangeUser(exp_xmin,exp_xmax);
   // gaus_errors->GetYaxis()->SetRangeUser(exp_ymin,exp_ymax);
   gaus_errors->GetXaxis()->SetTitleSize(0.05);
   gaus_errors->GetYaxis()->SetTitleSize(0.05);
   gaus_errors->GetXaxis()->SetTitle("DeltaSigma/Sigma");
   gaus_errors->GetYaxis()->SetTitle("DeltaAmplitude/Amplitude");
   gaus_errors->GetYaxis()->SetTitleOffset(1.5);
   gaus_errors->GetXaxis()->SetNdivisions(505);
   gaus_errors->SetMarkerColor(2);
   gaus_errors->SetMarkerSize(0.5);
   gaus_errors->SetLineColor(2);
   gaus_errors->SetFillColor(2);
   gaus_errors->SetFillStyle(11);
   gaus_errors->SetMarkerStyle(20);
   gaus_errors->SetTitle("");
   gaus_errors->Draw();


   c4->cd(3);
   exp_gauss_err1->GetXaxis()->SetRangeUser(0,200);
   // exp_gauss_err1->GetYaxis()->SetRangeUser(exp_ymin,exp_ymax);
   exp_gauss_err1->GetXaxis()->SetTitleSize(0.07);
   exp_gauss_err1->GetYaxis()->SetTitleSize(0.07);
   exp_gauss_err1->GetXaxis()->SetTitle("Err1 DAmp/Amp");
   exp_gauss_err1->GetYaxis()->SetTitle("events");
   exp_gauss_err1->GetXaxis()->SetNdivisions(505);
   exp_gauss_err1->SetMarkerColor(2);
   exp_gauss_err1->SetMarkerSize(0.5);
   exp_gauss_err1->SetLineColor(2);
   exp_gauss_err1->SetFillColor(2);
   exp_gauss_err1->SetFillStyle(11);
   exp_gauss_err1->SetMarkerStyle(20);
   exp_gauss_err1->SetTitle("");
   exp_gauss_err1->Draw();
   sprintf (string,"Err Fit Option='%s'\n",W.Data());
   printf("string=%s",string);
   t1->DrawLatex(0.16,0.92,string);
   t1->SetNDC();
   t1->SetTextSize(0.05);
   t1->Draw();


   c4->cd(4);
   exp_gauss_err2->GetXaxis()->SetRangeUser(0,200);
   // exp_gauss_err2->GetYaxis()->SetRangeUser(exp_ymin,exp_ymax);
   exp_gauss_err2->GetXaxis()->SetTitleSize(0.07);
   exp_gauss_err2->GetYaxis()->SetTitleSize(0.07);
   exp_gauss_err2->GetXaxis()->SetTitle("Err2 DAmp/Amp + DSig/Sig");
   exp_gauss_err2->GetYaxis()->SetTitle("events");
   exp_gauss_err2->GetXaxis()->SetNdivisions(505);
   exp_gauss_err2->SetMarkerColor(2);
   exp_gauss_err2->SetMarkerSize(0.5);
   exp_gauss_err2->SetLineColor(2);
   exp_gauss_err2->SetFillColor(2);
   exp_gauss_err2->SetFillStyle(11);
   exp_gauss_err2->SetMarkerStyle(20);
   exp_gauss_err2->SetTitle("");
   exp_gauss_err2->Draw();
   sprintf (string,"Err Fit Option='%s'\n",W.Data());
   printf("string=%s",string);
   t1->DrawLatex(0.16,0.92,string);
   t1->SetNDC();
   t1->SetTextSize(0.05);
   t1->Draw();


   

   if (W == 'W') {
        sprintf(filename,"Gaussian_Area_Errors_W1_%ds.pdf(",nsigma);
        c1->SaveAs(filename);
        sprintf(filename,"Gaussian_Area_Errors_W1_%ds_c1.png",nsigma);
        c1->SaveAs(filename);
         sprintf(filename,"Gaussian_Area_Errors_W1_%ds.pdf",nsigma);
         c3->SaveAs(filename);
         sprintf(filename,"Gaussian_Area_Errors_W1_%ds_c3.png",nsigma);
         c3->SaveAs(filename);
         sprintf(filename,"Gaussian_Area_Errors_W1_%ds.pdf)",nsigma);
         c4->SaveAs(filename);
         sprintf(filename,"Gaussian_Area_Errors_W1_%ds_c4.png",nsigma);
         c4->SaveAs(filename);
         }
    else {
        sprintf(filename,"Gaussian_Area_Errors_W0_%ds.pdf(",nsigma);
        c1->SaveAs(filename);
        sprintf(filename,"Gaussian_Area_Errors_W0_%ds_c1.png",nsigma);
        c1->SaveAs(filename);
         sprintf(filename,"Gaussian_Area_Errors_W0_%ds.pdf",nsigma);
         c3->SaveAs(filename);
         sprintf(filename,"Gaussian_Area_Errors_W0_%ds_c3.png",nsigma);
         c3->SaveAs(filename);
         sprintf(filename,"Gaussian_Area_Errors_W0_%ds.pdf)",nsigma);
         c4->SaveAs(filename);
         sprintf(filename,"Gaussian_Area_Errors_W0_%ds_c4.png",nsigma);
         c4->SaveAs(filename);
         }


}