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 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 background_subtraction(void)
{
// 
// Generate a Gaussian signal over a constant background.  Initial region of interest used is +/- 3 sigma
//                       Estimate background from fitted interval outside region of interest.
//                       Estimate signal and uncertainties and compare with generated distributions
//

// #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];
   // 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=3;           // signal region considered +/- nsigma
   Double_t mean1=0;   // Gaussian mean
   Double_t sigma1=1;
   TRandom1 *r1 = new TRandom1();
   TRandom1 *r2 = new TRandom1();
   TRandom1 *r3 = new TRandom1();
   TRandom1 *r4 = new TRandom1();

   // initialize random number seeds

   UInt_t  iseed = 23424111513;
   r1->SetSeed(iseed);
   iseed = 23424111591;
   r2->SetSeed(iseed);
   // r1->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::SetDefaultSumw2(kTRUE);
   TH1D  *exp_signal = new TH1D("exp_signal","Data - Fitted Back",exp_nbins,exp_xmin,exp_xmax);
   TH1D  *exp_gauss = new TH1D("exp_signal","Fitted Gauss",exp_nbins,exp_xmin,exp_xmax);
   TH1D  *exp_signal_err = new TH1D("exp_signal_err","Data - Fitted Back",exp_nbins,0,200);
   TH1D  *exp_gauss_err = new TH1D("exp_signal_err","Fitted Gauss",exp_nbins,0,200);

   TH1D  *gen_back = new TH1D("gen_back","Generated Back",exp_nbins,back_xmin,back_xmax);
   TH1D  *exp_back = new TH1D("exp_back","Fitted Back",exp_nbins,back_xmin,back_xmax);
   TH1D  *exp_back_err = new TH1D("exp_back_err","Fitted Back Err",exp_nbins,0,200);
 
   TH2D  *gaus_errors = new TH2D("gaus_erros","DeltaAmp/Amp vs DeltaSigma/Sigma",exp_nbins,0,0.1,exp_nbins,0,0.1);

   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);
    TH1D  *hsignal;
    TH1D  *hback;
    TH1D  *hback_pure;
    TH1D  *hdata;
    Double_t bin;
    Double_t xlo;
    Double_t xhi;
    Double_t back_counts=0;
    Double_t  data_counts=0;
    Double_t signal_counts=0;
    Double_t data_gaus;
    Double_t back_est;
    Double_t signal_est;
    Double_t signal_est_err;
    Double_t signal_gaus;
    Double_t signal_gaus_err;
    Double_t back_est_err;
    Double_t p1 = 50;
    
    Double_t amp;
    Double_t mean;
    Double_t sigma;
    Double_t amp_err;
    Double_t mean_err;
    Double_t sigma_err;
    
    
    
   // loop over 'expts' to collect statistics 

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

   hsignal = new TH1D("hsignal","Gen Signal",nbins,xmin,xmax);
   hback = new TH1D("hback","Gen Back",nbins,xmin,xmax);
   hback_pure = new TH1D("hback_pure","Gen Back Excluding Signal Region",nbins,xmin,xmax);
   hdata = new TH1D("hdata","Gen Data: Signal + background",nbins,xmin,xmax);

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

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

    // 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);
       hdata->Fill(z1);
       }

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

       // generate uniform background; same interval; npts_back events

       Double_t  z2 = xmin+(xmax-xmin)*r2->Uniform();
       hback->Fill(z2);
       hdata->Fill(z2);
       }

   // define fit function, assumed to be constant + Gaussian
  
   p1 = 50;

   amp =50;
   mean=mean1;
   sigma=sigma1;

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

   bin = (xmax-xmin)/nbins;
   xlo = mean1 - nsigma*sigma1;
   xhi = mean1 + nsigma*sigma1;
   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 = hback->GetBinLowEdge(jlo+1);
   xhi = hback->GetBinLowEdge(jhi)+hback->GetBinWidth(jhi);

   // first obtain background fit excluding signal region


   for(j=1;j<jlo+1;j++) {
       Int_t j1 = j;
       Double_t x = hdata->GetBinCenter(j1);
       Double_t entryb = hdata->GetBinContent(j1);
       hback_pure->Fill(x,entryb);
       // printf ("j1=%d, entryb=%g\n",j1,entryb);
       }

   for(j=jhi+1;j<nbins+1;j++) {
       Int_t j1 = j;
       Double_t x = hdata->GetBinCenter(j1);
       Double_t entryb = hdata->GetBinContent(j1);
       hback_pure->Fill(x,entryb);
       // printf ("j1=%d, entryb=%g\n",j1,entryb);
       }

   // Fit background excluding signal region

   fitback->SetParameter(0,p1);
   hback->SetLineColor(4);

   // get uncertainties using statistical weights

   hback_pure->Fit("fitback");
   Double_t p1_err = fitback->GetParError(0);

   // obtain number of events assuming uniform weights

   // hback_pure->Fit("fitback","W");
   hback_pure->Fit("fitback",W);
   p1 = fitback->GetParameter(0);
   fitfunc->FixParameter(3,p1);

   // get uncertainties using statistical weights

   hdata->Fit("fitfunc");
   amp_err = fitfunc->GetParError(0);
   mean_err = fitfunc->GetParError(1);
   sigma_err = fitfunc->GetParError(2);

   // obtain number of events assuming uniform weights

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

   // fit full background with constant weight

   /*hback->Fit("fitback","W");
   Double_t p1_full = fitback->GetParameter(0);
   Double_t p1_full_err = fitback->GetParError(0);
   printf ("p1_full=%g, p1_full_err=%g\n",p1_full,p1_full_err);*/

   // get number of signal events

   amp = fitfunc->GetParameter(0);
   mean = fitfunc->GetParameter(1);
   sigma = fitfunc->GetParameter(2);
   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(amp_err/amp,sigma_err/sigma);

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

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

       // printf ("j1=%d, x=%g, data=%g, back_counts=%g, back_fit=%g, signal_counts=%g\n",j1,x,entry,back_counts,back_fit,signal_counts);
       }

   signal_est = data_counts - back_est;
   signal_est_err = sqrt(data_counts);
   signal_gaus = data_gaus- back_est;
   // signal_gaus_err = signal_gaus*sqrt((amp_err/amp)*(amp_err/amp));  // Equivalent to line below
   signal_gaus_err = amp_err/bin;  // Note that fit definition includes normalization so sigma error plays no role.
   printf ("data_counts = %g, back_counts = %g, signal_counts = %g, back_fit=%g, back_est=%g, signal_est=%g, signal_gaus=%g\n",data_counts,back_counts,signal_counts,back_fit, back_est, signal_est,signal_gaus);
   back_est_err=sqrt(back_fit);

  // fill summary histogram with 'exptl' results

   exp_signal->Fill(signal_est);
   exp_gauss->Fill(signal_gaus);
   exp_signal_err->Fill(signal_est_err);
   exp_gauss_err->Fill(signal_gaus_err);
   gen_back->Fill(back_counts);
   exp_back->Fill(back_fit);
   exp_back_err->Fill(back_est_err);


   // delete histograms and repeat to build up statistics

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

   // end of loop over number of experiments

   }


   //
   TCanvas *c1 = new TCanvas("c1","c1 Gen Signal and Back",200,10,700,700);
   // c1->SetGridx();
   // c1->SetGridy();
   // c1->SetLogy();


   TLegend *leg = new TLegend(0.5,0.80,0.85,0.9);
   leg->AddEntry(hdata,"Data (Signal+Background)","p");
   leg->AddEntry(hback,"Generated Background","l");

   hdata->GetXaxis()->SetRangeUser(xmin,xmax);
   hdata->GetYaxis()->SetRangeUser(ymin,ymax);
   hdata->GetXaxis()->SetTitleSize(0.07);
   hdata->GetYaxis()->SetTitleSize(0.07);
   hdata->GetXaxis()->SetTitle("variable x");
   hdata->GetYaxis()->SetTitle("events");
   hdata->SetMarkerColor(2);
   hdata->SetMarkerSize(0.5);
   hdata->SetLineColor(2);
   hdata->SetFillColor(2);
   hdata->SetFillStyle(11);
   hdata->SetMarkerStyle(20);
   hdata->SetTitle("");
   hdata->SetStats(kFALSE);

   hdata->Draw("pe");
   hback->Draw("same");
   fitback->Draw("same");
   leg->AddEntry(fitback,"Fitted Background","l");
   leg->SetTextSize(0.02);
   leg->Draw();
   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.025);
   t1->Draw();


   // ouput information

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

   sprintf (string,"Gen signal mean=%.2f\n",mean1);
   printf("string=%s",string);
   t1->DrawLatex(0.16,0.76,string);
   t1->SetNDC();
   t1->SetTextSize(0.02);
   t1->Draw();
   sprintf (string,"Gen signal sigma=%.2f\n",sigma1);
   printf("string=%s",string);
   t1->DrawLatex(0.16,0.72,string);
   t1->SetNDC();
   t1->SetTextSize(0.02);
   t1->Draw();
   sprintf (string,"Fit signal mean=%.2f#pm%.2f\n",mean,mean_err);
   printf("string=%s",string);
   t1->DrawLatex(0.16,0.68,string);
   t1->SetNDC();
   t1->SetTextSize(0.02);
   t1->Draw();
   sprintf (string,"Fit signal sigma=%.2f#pm%.2f\n",sigma,sigma_err);
   printf("string=%s",string);
   t1->DrawLatex(0.16,0.64,string);
   t1->SetNDC();
   t1->SetTextSize(0.02);
   t1->Draw();


   TCanvas *c2 = new TCanvas("c2","c2 Data excluding signal region",200,10,700,700);
   c2->Divide(2,2);
   // c2->SetGridx();
   // c2->SetGridy();
   // c2->SetLogy();

   c2->cd(1);
   hback_pure->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();


   c2->cd(2);
   gaus_errors->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("Data - Fitted Back");
   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("Data - Fitted Back Err");
   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("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);
   gen_back->GetXaxis()->SetRangeUser(back_xmin,back_xmax);
   // gen_back->GetYaxis()->SetRangeUser(back_ymin,back_ymax);
   gen_back->GetXaxis()->SetTitleSize(0.07);
   gen_back->GetYaxis()->SetTitleSize(0.07);
   gen_back->GetXaxis()->SetTitle("Gen Back");
   gen_back->GetYaxis()->SetTitle("events");
   gen_back->GetXaxis()->SetNdivisions(505);
   gen_back->SetMarkerColor(2);
   gen_back->SetMarkerSize(0.5);
   gen_back->SetLineColor(2);
   gen_back->SetFillColor(2);
   gen_back->SetFillStyle(11);
   gen_back->SetMarkerStyle(20);
   gen_back->SetTitle("");
   gen_back->Fit("gaus");
   gen_back->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(2);
   exp_back->GetXaxis()->SetRangeUser(back_xmin,back_xmax);
   // exp_back->GetYaxis()->SetRangeUser(back_ymin,back_ymax);
   exp_back->GetXaxis()->SetTitleSize(0.07);
   exp_back->GetYaxis()->SetTitleSize(0.07);
   exp_back->GetXaxis()->SetTitle("Fitted Back");
   exp_back->GetYaxis()->SetTitle("events");
   exp_back->GetXaxis()->SetNdivisions(505);
   exp_back->SetMarkerColor(2);
   exp_back->SetMarkerSize(0.5);
   exp_back->SetLineColor(2);
   exp_back->SetFillColor(2);
   exp_back->SetFillStyle(11);
   exp_back->SetMarkerStyle(20);
   exp_back->SetTitle("");
   exp_back->Fit("gaus");
   exp_back->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_back_err->GetXaxis()->SetRangeUser(0,200);
   // exp_back_err->GetYaxis()->SetRangeUser(0,200);
   exp_back_err->GetXaxis()->SetTitleSize(0.07);
   exp_back_err->GetYaxis()->SetTitleSize(0.07);
   exp_back_err->GetXaxis()->SetTitle("Fitted Back Err");
   exp_back_err->GetYaxis()->SetTitle("events");
   exp_back_err->GetXaxis()->SetNdivisions(505);
   exp_back_err->SetMarkerColor(2);
   exp_back_err->SetMarkerSize(0.5);
   exp_back_err->SetLineColor(2);
   exp_back_err->SetFillColor(2);
   exp_back_err->SetFillStyle(11);
   exp_back_err->SetMarkerStyle(20);
   exp_back_err->SetTitle("");
   exp_back_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();
   
   if (W == 'W') {
         sprintf(filename,"background_subtraction_W1_%ds.pdf(",nsigma);
         c1->SaveAs(filename);
         sprintf(filename,"background_subtraction_W1_%ds_c1.png",nsigma);
         c1->SaveAs(filename);
         sprintf(filename,"background_subtraction_W1_%ds.pdf",nsigma);
         c3->SaveAs(filename);
         sprintf(filename,"background_subtraction_W1_%ds_c3.png",nsigma);
         c3->SaveAs(filename);
         sprintf(filename,"background_subtraction_W1_%ds.pdf)",nsigma);
         c4->SaveAs(filename);
         sprintf(filename,"background_subtraction_W1_%ds_c4.png",nsigma);
         c4->SaveAs(filename);
         }
    else {
         sprintf(filename,"background_subtraction_W0C_%ds.pdf(",nsigma);
         c1->SaveAs(filename);
         sprintf(filename,"background_subtraction_W0C_%ds_c1.png",nsigma);
         c1->SaveAs(filename);
         sprintf(filename,"background_subtraction_W0C_%ds.pdf",nsigma);
         c3->SaveAs(filename);
         sprintf(filename,"background_subtraction_W0C_%ds_c3.png",nsigma);
         c3->SaveAs(filename);
         sprintf(filename,"background_subtraction_W0C_%ds.pdf)",nsigma);
         c4->SaveAs(filename);
         sprintf(filename,"background_subtraction_W0C_%ds_c4.png",nsigma);
         c4->SaveAs(filename);
         }


}