void gen_stat_test(void)
{
// 
// generate test statistics assuming a uniform background in two dimensions and a signal which is the sum of the two variables.
//

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

gROOT->Reset();
//TTree *Bfield = (TTree *) gROOT->FindObject("Bfield");
gStyle->SetPalette(1,0);
gStyle->SetOptStat(kTRUE);
// gStyle->SetOptFit(kTRUE);
// gStyle->SetOptFit(1111);
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 npts 10000;
   # define nbins 100;

   // define parameters

   Double_t a=10;
   Double_t sigma=0.25;
   Double_t norm=1;
   Double_t tau=a/3;

   // define limits

   Double_t xmin=-5*sigma;
   Double_t xmax=a+5*sigma;

   sprintf (string,"Box length =%f, sigma=%f, norm=%f\n",a,sigma,norm);
   printf ("plot_uniform_smear: %s",string);

   TF1 *boxsmear = new TF1("boxsmear",boxsmear,xmin,xmax,2);
   boxsmear->SetParameters(a,sigma);


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

   Double_t ymin = 0;
   Double_t ymax = 2;
   boxsmear->GetXaxis()->SetRangeUser(xmin,xmax);
   boxsmear->GetYaxis()->SetRangeUser(ymin,ymax);
   boxsmear->GetXaxis()->SetTitleSize(0.05);
   boxsmear->GetYaxis()->SetTitleSize(0.05);
   boxsmear->GetYaxis()->SetTitleOffset(1.5);
   boxsmear->GetXaxis()->SetTitle("variable x");
   boxsmear->GetYaxis()->SetTitle("probability");
   boxsmear->GetXaxis()->SetNdivisions(5);
   // boxsmear->SetMarkerColor(1);
   // boxsmear->SetMarkerStyle(21);
   boxsmear->SetTitle("");
   boxsmear->Draw("");
   boxsmear->SetLineColor(2);

   sprintf(string,"box size= %.1f\n",a);
   t1 = new TLatex(0.2,0.85,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->Draw();
   sprintf(string,"#sigma = %.1f\n",sigma);
   t1 = new TLatex(0.2,0.80,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->Draw();
   

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

   TH1F *back1 = new TH1F("back1","Background 1",nbins,xmin,xmax);


   // back1->Fit("boxsmear","","",xmin,xmax);

   back1->SetTitle("");
   // back1->GetXaxis()->SetRangeUser(1.0,2.5);
   // back1->GetYaxis()->SetRangeUser(ymin,ymax);
   back1->GetXaxis()->SetTitleSize(0.05);
   back1->GetYaxis()->SetTitleSize(0.05);
   back1->GetYaxis()->SetTitleOffset(1.5);
   back1->GetYaxis()->SetTitle("events per bin");
   back1->GetXaxis()->SetTitle("Background 1");
   back1->GetXaxis()->SetNdivisions(505);
   back1->SetMarkerColor(2);
   back1->SetMarkerStyle(21);
   back1->Draw();

   sprintf(string,"Entries = %.0f\n",back1->GetEntries());
   t1 = new TLatex(0.15,0.92,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->Draw();

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

   TH1F *back2 = new TH1F("back2","Background 2",nbins,xmin,xmax);
   TH1F *backtot = new TH1F("backtot","Background 1+Background 2",nbins,xmin,2*xmax);
   TH2F *back = new TH2F("back","background correlation",nbins,xmin,xmax,nbins,xmin,xmax);

   TRandom1 *r = new TRandom1();

   for (j=0;j<npts;j++) {
      // Double_t b1 = boxsmear->GetRandom();
      Double_t b1 = r->Exp(tau);
      back1->Fill(b1);
      // Double_t b2 = boxsmear->GetRandom();
      Double_t b2 = r->Exp(tau);
      back2->Fill(b2);
      back->Fill(b1,b2);
      backtot->Fill(b1+b2);
      }

   // back2->Fit("boxsmear","","",xmin,xmax);

   back2->SetTitle("");
   // back2->GetXaxis()->SetRangeUser(1.0,2.5);
   // back2->GetYaxis()->SetRangeUser(ymin,ymax);
   back2->GetXaxis()->SetTitleSize(0.05);
   back2->GetYaxis()->SetTitleSize(0.05);
   back2->GetYaxis()->SetTitleOffset(1.5);
   back2->GetYaxis()->SetTitle("events per bin");
   back2->GetXaxis()->SetTitle("Background 2");
   back2->GetXaxis()->SetNdivisions(505);
   back2->SetMarkerColor(2);
   back2->SetMarkerStyle(21);
   back2->Draw();

   sprintf(string,"Entries = %.0f\n",back2->GetEntries());
   t1 = new TLatex(0.15,0.92,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->Draw();

   Double_t Etot = a;
  

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

   TH1F *sig1 = new TH1F("sig1","signal1",nbins,xmin,xmax);

   // sig1->Fit("boxsmear","","",xmin,xmax);

   sig1->SetTitle("");
   // sig1->GetXaxis()->SetRangeUser(1.0,2.5);
   // sig1->GetYaxis()->SetRangeUser(ymin,ymax);
   sig1->GetXaxis()->SetTitleSize(0.05);
   sig1->GetYaxis()->SetTitleSize(0.05);
   sig1->GetYaxis()->SetTitleOffset(1.5);
   sig1->GetYaxis()->SetTitle("events per bin");
   sig1->GetXaxis()->SetTitle("Signal 1");
   sig1->GetXaxis()->SetNdivisions(505);
   sig1->SetMarkerColor(2);
   sig1->SetMarkerStyle(21);
   sig1->Draw();

   sprintf(string,"Entries = %.0f\n",sig1->GetEntries());
   t1 = new TLatex(0.15,0.92,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->Draw();

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

   TH1F *sig2 = new TH1F("sig2","signal2",nbins,xmin,xmax);
   TH1F *sigtot = new TH1F("sigtot","signal1+signal2",nbins,xmin,2*xmax);
   TH2F *sig = new TH2F("sig","signal correlation",nbins,xmin,xmax,nbins,xmin,xmax);

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

      // pick energies

      Double_t x = r->Uniform();
      Double_t E10 = Etot*x;
      Double_t E20 = Etot*(1-x);
      Double_t E1, E2;
      E1= E10 + r->Gaus(0,sigma);
      E2 = E20 + r->Gaus(0,sigma);
      sig1->Fill(E1);
      sig2->Fill(E2);
      sig->Fill(E1,E2);
      sigtot->Fill(E1+E2);
      }

   // sig2->Fit("boxsmear","","",xmin,xmax);

   sig2->SetTitle("");
   // sig2->GetXaxis()->SetRangeUser(1.0,2.5);
   // sig2->GetYaxis()->SetRangeUser(ymin,ymax);
   sig2->GetXaxis()->SetTitleSize(0.05);
   sig2->GetYaxis()->SetTitleSize(0.05);
   sig2->GetYaxis()->SetTitleOffset(1.5);
   sig2->GetYaxis()->SetTitle("events per bin");
   sig2->GetXaxis()->SetTitle("Signal 2");
   sig2->GetXaxis()->SetNdivisions(505);
   sig2->SetMarkerColor(2);
   sig2->SetMarkerStyle(21);
   sig2->Draw();

   sprintf(string,"Entries = %.0f\n",sig2->GetEntries());
   t1 = new TLatex(0.15,0.92,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->Draw();

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


   // sig->Fit("boxsmear","","",xmin,xmax);

   sig->SetTitle("");
   // sig->GetXaxis()->SetRangeUser(1.0,2.5);
   sig->GetXaxis()->SetTitleSize(0.05);
   sig->GetYaxis()->SetTitleSize(0.05);
   sig->GetYaxis()->SetTitleOffset(1.5);
   sig->GetYaxis()->SetTitle("Signal 2 or Background 2");
   sig->GetXaxis()->SetTitle("Signal1 or Background 2");
   sig->GetXaxis()->SetNdivisions(505);
   sig->SetMarkerColor(2);
   sig->SetMarkerStyle(20);
   sig->SetMarkerSize(0.3);
   sig->SetLineColor(2);
   sig->Draw();


   back->SetTitle("");
   // back->GetXaxis()->SetRangeUser(1.0,2.5);
   back->GetXaxis()->SetTitleSize(0.05);
   back->GetYaxis()->SetTitleSize(0.05);
   back->GetYaxis()->SetTitleOffset(1.5);
   back->GetYaxis()->SetTitle("events per bin");
   back->GetXaxis()->SetTitle("random variable x");
   back->GetXaxis()->SetNdivisions(505);
   back->SetMarkerColor(4);
   back->SetMarkerStyle(20);
   back->SetMarkerSize(0.3);
   back->SetLineColor(4);
   back->Draw("same");

   sprintf(string,"Entries = %.0f\n",back->GetEntries());
   t1 = new TLatex(0.15,0.92,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->Draw();


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

   sigtot->SetTitle("");
   // sigtot->GetXaxis()->SetRangeUser(1.0,2.5);
   sigtot->GetXaxis()->SetTitleSize(0.05);
   sigtot->GetYaxis()->SetTitleSize(0.05);
   sigtot->GetYaxis()->SetTitleOffset(1.5);
   sigtot->GetYaxis()->SetTitle("events per bin");
   sigtot->GetXaxis()->SetTitle("Sum sig1 + sig2");
   sigtot->GetXaxis()->SetNdivisions(505);
   sigtot->SetMarkerColor(2);
   sigtot->SetMarkerStyle(20);
   sigtot->SetMarkerSize(0.3);
   sigtot->SetLineColor(2);
   sigtot->Draw();

   backtot->SetTitle("");
   // backtot->GetXaxis()->SetRangeUser(1.0,2.5);
   backtot->GetXaxis()->SetTitleSize(0.05);
   backtot->GetYaxis()->SetTitleSize(0.05);
   backtot->GetYaxis()->SetTitleOffset(1.5);
   backtot->GetYaxis()->SetTitle("events per bin");
   backtot->GetXaxis()->SetTitle("Sum sig1 + sig2");
   backtot->GetXaxis()->SetNdivisions(505);
   backtot->SetMarkerColor(4);
   backtot->SetMarkerStyle(20);
   backtot->SetMarkerSize(0.3);
   backtot->SetLineColor(4);
   backtot->Draw("same");

   c1->SaveAs("gen_stat_test_c1.png");
   c2->SaveAs("gen_stat_test_c2.png");
   c3->SaveAs("gen_stat_test_c3.png");
   c4->SaveAs("gen_stat_test_c4.png");
   c5->SaveAs("gen_stat_test_c5.png");
   c6->SaveAs("gen_stat_test_c6.png");
   c7->SaveAs("gen_stat_test_c7.png");
     
}
Double_t boxsmear(Double_t *x, Double_t *par) {

   // box [-a/2,a/2] smeared with Gaussian with width sigma

   Double_t a=par[0];
   Double_t sigma=par[1];
   Double_t x1=x[0]-a/2;

   char string[256];
   Double_t pi=3.14159;

   Double_t arg;
   Double_t distr;

   if (x1 > 0) {
         arg = (x1-a/2)/sigma;
         distr = 0.5-0.5*TMath::Erf(arg/sqrt(2));
         }
        else {
         arg = (x1+a/2)/sigma;
         distr = 0.5 +0.5*TMath::Erf(arg/sqrt(2));
        }


   sprintf(string,"a=%g, sigma=%g, x=%g, distr=%g\n",a,sigma,x1,distr);
   // printf ("func: %s",string);


    return distr;
}