void plot_binomial(void)
{
// 
// plot various probability density distributions
//

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

   // define parameters

   Double_t mean=5;
   Double_t sigma=1;
   Double_t norm=1;
   Int_t ntot = npts;
   Double_t prob=0.5;

   // define limits

   Double_t xmin=0;
   Double_t xmax=ntot;

   sprintf (string,"Gaussian: mean=%f, sigma=%f, norm=%f\n",mean,sigma,norm);
   printf ("plot_binomial: %s",string);

   TF1 *func = new TF1("func",func,xmin,xmax,3);
   func->SetParameters(mean,sigma,norm);

   TF1 *binfunc = new TF1("binfunc",binfunc,xmin,xmax,2);
   binfunc->SetParameters(ntot,prob);

   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 = 10;
   func->GetXaxis()->SetRangeUser(xmin,xmax);
   func->GetYaxis()->SetRangeUser(ymin,ymax);
   func->GetXaxis()->SetTitleSize(0.05);
   func->GetYaxis()->SetTitleSize(0.05);
   func->GetYaxis()->SetTitleOffset(1.5);
   func->GetXaxis()->SetTitle("variable x");
   func->GetYaxis()->SetTitle("probability");
   func->GetXaxis()->SetNdivisions(5);
   // func->SetMarkerColor(1);
   // func->SetMarkerStyle(21);
   func->SetTitle("");
   func->Draw("");

   sprintf(string,"#mu = %.1f\n",mean);
   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();

   TLine *lmean = new TLine(mean,ymin,mean,ymax);
   // lmean->SetNDC();
   lmean->Draw();
   TLine *lsig1 = new TLine(mean-sigma,ymin,mean-sigma,ymax);
   // lsig1->Draw();
   TLine *lsig2 = new TLine(mean+sigma,ymin,mean+sigma,ymax);
   // lsig2->Draw();

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

   TH1F *pdf = new TH1F("pdf","random generation",nbins,xmin,xmax);
   TH1F *pdfhist = new TH1F("pdfhist","Fill hist with pdf",nbins,xmin,xmax);

   TRandom1 *r = new TRandom1();
   for (j=0;j<npts;j++) {
      Int_t rbin = r->Binomial(ntot,prob);
      pdf->Fill((Double_t) rbin);
      }

    for (j=0;j<ntot;j++) {
         Double_t bin = (xmax-xmin)/npts;
         Double_t x = xmin + j*bin;
         pdfhist->Fill(x,binfunc(x)*(Double_t) ntot);
        }

   // pdf->Fit("func","","",xmin,xmax);

   TLegend *leg = new TLegend(0.2,0.75,0.85,0.9);
   leg->AddEntry(pdf,"Random generation","l");
   leg->AddEntry(pdfhist,"PDF normalized to N","p");

   pdf->SetTitle("");
   // pdf->GetXaxis()->SetRangeUser(1.0,2.5);
   pdf->GetYaxis()->SetRangeUser(ymin,ymax);
   pdf->GetXaxis()->SetTitleSize(0.05);
   pdf->GetYaxis()->SetTitleSize(0.05);
   pdf->GetYaxis()->SetTitleOffset(1.5);
   pdf->GetYaxis()->SetTitle("events per bin");
   pdf->GetXaxis()->SetTitle("random variable n");
   pdf->GetXaxis()->SetNdivisions(505);
   pdf->DrawCopy();
   // binfunc->Draw("same");
   pdfhist->SetMarkerColor(2);
   pdfhist->SetMarkerStyle(21);
   pdfhist->Draw("samep");
   leg->Draw();

   sprintf(string,"N = %d  Probability= %.2f\n",ntot,prob);
   t1 = new TLatex(0.15,0.92,string);
   t1->SetTextColor(1);
   t1->SetNDC();
   t1->Draw();

   c1->SaveAs("plot_binomial_c1.png");
   c3->SaveAs("plot_binomial_c3.png");
     
}
Double_t func(Double_t *x, Double_t *par) {

   // function to be plotted

   Double_t p1=par[0];
   Double_t p2=par[1];
   Double_t p3=par[2];
   Double_t x1=x[0];

   char string[256];
   Double_t pi=3.14159;

   Double_t distr = exp(-(x1-p1)*(x1-p1)/(2*p2*p2))/(sqrt(2*pi)*p2);

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

   return p3*distr;
}
Double_t binfunc(Double_t *x, Double_t *par) {

   // return real value of binomial distribution

   Double_t ntot=par[0];
   Double_t p=par[1];
   Double_t x1=x[0];

   char string[256];
   Double_t pi=3.14159;

   Double_t q = 1-p;
   Double_t x2 = ntot-x1;

   Double_t distr = log(TMath::Gamma(ntot+1)) - log(TMath::Gamma(x1+1)) - log(TMath::Gamma(x2+1)) + x1*log(p) + x2*log(q);

   sprintf(string,"ntot=%g, p=%g, x1=%g, x2=%g, distr=%g\n",ntot,p,x1,x2,exp(distr));
   // printf ("func: %s",string);

   Double_t factor=1;
   return factor*exp(distr);
}