void plot_background(void)
{
// 
// plot a Poisson-generated signal over a polynomial background and fit data to extract signal.
//

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

gROOT->Reset();
//TTree *Bfield = (TTree *) gROOT->FindObject("Bfield");
gStyle->SetPalette(1,0);
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];
   Int_t j,jj;
   
   #define npts1 10000;
   #define npts2 100;

   Double_t Pmean1=10;   // Poisson mean
   Double_t Pmean2=13;   // Poisson mean
   Double_t sigma1=15;
   Double_t sigma2=1;
   TRandom1 *r1 = new TRandom1();
   TRandom1 *r2 = new TRandom1();

   // define histograms

   Int_t nbins=100;
   Float_t xmin=0;
   Float_t xmax=30;
   Float_t ymin=0;
   Float_t ymax=300;
   TH1D  *back= new TH1D("back","Back",nbins,xmin,xmax);
   TH1D  *data= new TH1D("data","Signal and background",nbins,xmin,xmax);

   //
   TCanvas *c1 = new TCanvas("c1","c1 Plot Poisson",200,10,700,700);
   c1->SetBorderMode(0);
   c1->SetFillColor(0);
   // c1->SetGridx();
   // c1->SetGridy();
   // c1->SetLogy();

    // generate data

   UInt_t  iseed = 23424111515;
   r1->SetSeed(iseed);
   // r1->SetSeed();
   for (j=0;j<npts1;j++) {

       // generate background data assuming Gaussian form with wide width

       Double_t z1 = r1->Gaus(Pmean1,sigma1);
       // Int_t n1 = r1->Poisson(Pmean1);
       // back->Fill((Double_t) n1);
       back->Fill(z1);
       data->Fill(z1);
       // data->Fill((Double_t) n1);
       }

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

       // generate Gaussian signal and add to data

       // Int_t n2 = r2->Poisson(Pmean2);
       Double_t  z2 = r2->Gaus(Pmean2,sigma2);
       // data->Fill((Double_t) n2);
      data->Fill(z2);
       }

   // define fit function, assumed to be polynomial
  
   Double_t p1 = 50;
   Double_t p2 = -1;
   Double_t p3 = 0;

   /*Double_t p1 = 100;
   Double_t p2 = Pmean2;
   Double_t p3 = sigma2;*/

   Double_t amp =600;
   Double_t mean=Pmean2;
   Double_t sigma=sigma2;

   TF1 *fitfunc = new TF1("fitfunc",fitfunc,xmin,xmax,6);
   fitfunc->SetParameters(p1,p2,p3,amp,mean,sigma);

   TLegend *leg = new TLegend(0.2,0.75,0.85,0.9);
   leg->AddEntry(data,"Data (Signal+Background)","p");
   leg->AddEntry(back,"Generated Background","l");

   data->GetXaxis()->SetRangeUser(xmin,xmax);
   data->GetYaxis()->SetRangeUser(ymin,ymax);
   data->GetXaxis()->SetTitleSize(0.07);
   data->GetYaxis()->SetTitleSize(0.07);
   data->GetXaxis()->SetTitle("variable x");
   data->GetYaxis()->SetTitle("events");
   data->SetMarkerColor(2);
   data->SetMarkerSize(0.5);
   data->SetLineColor(2);
   data->SetFillColor(2);
   data->SetFillStyle(11);
   data->SetMarkerStyle(20);
   data->SetTitle("");
   // data->Fit("fitfunc","","",5,20);
   // data->Fit("gaus");
   // data->Draw("Ape");

   //
   // data->Fit("fitfunc","0");
   // data->Fit("fitfunc");
   p1 = fitfunc->GetParameter(0);
   p2 = fitfunc->GetParameter(1);
   p3 = fitfunc->GetParameter(2);

   // define background

   TF1 *fitback = new TF1("fitback","[0]+[1]*x+[2]*x*x",xmin,xmax);
   // TF1 *fitback = new TF1("fitback","[0]*exp(-(x-[1])*(x-[1])/(2*[2]*[2]))/sqrt(2*3.14159*[2]*[2])",xmin,xmax);
   fitback->SetParameters(p1,p2,p3);
   back->SetLineColor(4);
   back->Fit("fitback");
   p1 = fitback->GetParameter(0);
   p2 = fitback->GetParameter(1);
   p3 = fitback->GetParameter(2);
   fitfunc->FixParameter(0,p1);
   fitfunc->FixParameter(1,p2);
   fitfunc->FixParameter(2,p3);
   // data->Fit("fitfunc","0");
   data->Fit("fitfunc");
   data->Draw("pe");
   back->Draw("same");
   fitback->Draw("same");
   leg->AddEntry(fitback,"Fitted Background","l");
   leg->SetTextSize(0.04);
   leg->Draw();

   // get number of signal events

   amp = fitfunc->GetParameter(3);
   mean = fitfunc->GetParameter(4);
   sigma = fitfunc->GetParameter(5);
   Double_t amp_err = fitfunc->GetParError(3);
   Double_t mean_err = fitfunc->GetParError(4);
   Double_t sigma_err = fitfunc->GetParError(5);
   printf ("p1=%g, p2=%g, p3=%g, amp=%g, mean=%g, sigma=%g, amp_err=%f, mean_err=%f, sigma_err=%f\n",p1,p2,p3,amp,mean,sigma,amp_err,mean_err,sigma_err);

   Double_t bin = (xmax-xmin)/nbins;
   Double_t xlo = mean - 3*sigma;
   Double_t xhi = mean + 3*sigma;
   Int_t jlo = (xlo-xmin)/bin;
   Int_t jhi = (xhi-xmin)/bin;

   printf ("xlo = %g, xhi = %g, jlo=%d, jhi=%d \n",xlo, xhi, jlo, jhi);

   // recompute xlo and xhi to bin centers

   xlo = back->GetBinLowEdge(jlo+1);
   xhi = back->GetBinCenter(jhi+1)+back->GetBinWidth(jhi+1);

   Double_t data_gaus = fitfunc->Integral(xlo,xhi)/bin;
   printf ("bin=%g, data_gaus=%g\n",bin,data_gaus);
   Double_t back_est= fitback->Integral(xlo,xhi)/bin;
   Double_t back_counts, data_counts, signal_counts;

   printf ("xlo = %g, xhi = %g, jlo=%d, jhi=%d, back_est = %g\n",xlo, xhi, jlo, jhi, back_est);
  
   for(j=jlo;j<jhi+1;j++) {
       Double_t x = back->GetBinCenter(j+1);
       Double_t entryb = back->GetBinContent(j+1);
       back_counts = back_counts + entryb;
       Double_t x = data->GetBinCenter(j+1);
       Double_t entry = data->GetBinContent(j+1);
       data_counts = data_counts + entry;
       signal_counts = data_counts - back_counts;
       Double_t back_fit = fitback->Eval(x);
;
       printf ("j+1=%d, x=%g, data=%g, back=%g, back_fit=%g, signal_counts=%g\n",j+1,x,entry,entryb,back_fit,signal_counts);
       }

   Double_t signal_est = data_counts - back_est;
   Double_t signal_est_err = sqrt(data_counts);
   Double_t signal_gaus = data_gaus- back_est;
   Double_t signal_gaus_err = signal_gaus*sqrt((amp_err/amp)*(amp_err/amp));
   printf ("data_counts = %g, back_counts = %g, signal_counts = %g, back_est=%g, signal_est=%g, signal_gaus=%g\n",data_counts,back_counts,signal_counts,back_est,signal_est,signal_gaus);
   //printf ("data_counts = %g, back_counts = %g, signal_counts = %g, back_est=%g, signal_est=%g, signal_gaus=%g\n",data_counts,back_counts,signal_counts,back_est,signal_est);

   // ouput information

   sprintf (string,"Interval of interest #pm 3#sigma(%.1f-%.1f)\n",xlo,xhi);
   printf("string=%s",string);
   t1 = new TLatex(0.5,0.68,string);
   t1->SetNDC();
   t1->SetTextSize(0.03);
   t1->Draw();
   sprintf (string,"Data events= %.1f\n",data_counts);
   printf("string=%s",string);
   t1 = new TLatex(0.5,0.64,string);
   t1->SetNDC();
   t1->SetTextSize(0.03);
   t1->Draw();
   sprintf (string,"Gen back events= %.1f\n",back_counts);
   printf("string=%s",string);
   t1 = new TLatex(0.5,0.60,string);
   t1->SetNDC();
   t1->SetTextSize(0.03);
   t1->Draw();
   sprintf (string,"Gen signal events= %.1f\n",signal_counts);
   printf("string=%s",string);
   t1 = new TLatex(0.5,0.56,string);
   t1->SetNDC();
   t1->SetTextSize(0.03);
   t1->Draw();
   sprintf (string,"Background estimate= %.1f#pm%.1f\n",back_est,sqrt(back_est));
   printf("string=%s",string);
   t1 = new TLatex(0.5,0.52,string);
   t1->SetNDC();
   t1->SetTextSize(0.03);
   t1->Draw();
   sprintf (string,"Signal estimate= %.1f#pm%.1f\n",signal_est,signal_est_err);
   printf("string=%s",string);
   t1 = new TLatex(0.5,0.48,string);
   t1->SetNDC();
   t1->SetTextSize(0.03);
   t1->Draw();
   sprintf (string,"Gauss estimate= %.1f#pm%.1f \n",signal_gaus,signal_gaus_err);
   printf("string=%s",string);
   t1 = new TLatex(0.5,0.44,string);
   t1->SetNDC();
   t1->SetTextSize(0.03);
   t1->Draw();

   sprintf (string,"Gen signal mean=%.1f\n",Pmean2);
   printf("string=%s",string);
   t1 = new TLatex(0.16,0.68,string);
   t1->SetNDC();
   t1->SetTextSize(0.03);
   t1->Draw();
   sprintf (string,"Gen signal sigma=%.1f\n",sigma2);
   printf("string=%s",string);
   t1 = new TLatex(0.16,0.64,string);
   t1->SetNDC();
   t1->SetTextSize(0.03);
   t1->Draw();
   sprintf (string,"Fit signal mean=%.1f#pm%.1f\n",mean,mean_err);
   printf("string=%s",string);
   t1 = new TLatex(0.16,0.60,string);
   t1->SetNDC();
   t1->SetTextSize(0.03);
   t1->Draw();
   sprintf (string,"Fit signal sigma=%.1f#pm%.1f\n",sigma,sigma_err);
   printf("string=%s",string);
   t1 = new TLatex(0.16,0.56,string);
   t1->SetNDC();
   t1->SetTextSize(0.03);
   t1->Draw();
   
   c1->SaveAs("plot_background_c1.pdf");
   c1->SaveAs("plot_background_c1.png");

}
Double_t fitfunc(Double_t *x, Double_t *par) {

   // function to be plotted: Polynomial for background + Gaussian signal.

   Double_t p1=par[0];
   Double_t p2=par[1];
   Double_t p3=par[2];
   Double_t amp=par[3];
   Double_t mean=par[4];
   Double_t sigma=par[5];
   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 back = p1 + p2*x1 + p3*x1*x1;   // polynomial background
   // Double_t arg2 = (x1-p2)*(x1-p2)/(2*p3*p3);
   // Double_t back =p1*exp(-arg2)/sqrt(2*pi*p3*p3);
   Double_t signal = amp*exp(-arg)/sqrt(2*pi*sigma*sigma);

   sprintf(string,"p1=%g, p2=%g, p3=%g, x1=%g, amp=%g, mean=%g, sigma=%g, fitfunc=%g\n",p1,p2,p3,x1,amp,mean,sigma,back+signal);
   // printf ("func: %s",string);

   return back+signal;
}