void plot_box_smear(void)
{
// 
// plot smeared box distribution on [-a/2,a/2] with Gaussian width sigma
//

// #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(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 100000;
   # define nbins 100;

   // define parameters

   Double_t a=20;
   Double_t sigma=2;
   Double_t norm=1;

   // define limits

   Double_t xmin=-2*a;
   Double_t xmax=2*a;

   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 *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 *box = new TH1F("box","random uniform smear",nbins,xmin,xmax);

   for (j=0;j<npts;j++) {
      Double_t rand = boxsmear->GetRandom();
      box->Fill(rand);
      }

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

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

   sprintf(string,"Entries = %.0f  nbins = %.0f  #Delta x = %.1f\n",box->GetEntries(),(Float_t )nbins,(xmax-xmin)/nbins);
   t1 = new TLatex(0.15,0.92,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->Draw();

   c1->SaveAs("plot_box_smear_c1.png");
   c3->SaveAs("plot_box_smear_c3.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);

   // distr = 1;
    // if (x1 < -a/2 || x1 > a/2) distr = 0;

    return distr;
}