#include "TRandom.h"
#include "TGraphErrors.h"

//---------------------------------------------------
double poisson(double nbar, double n)
//---------------------------------------------------
{
  double u = 1;
  for( int i=0 ; i<n ; i++ )
  {
    u *= nbar/(i+1);
  }
  return u*exp(-nbar);
}

void PoissonConvergence()
{
  double rnd;
  int Nthrows=50000;
  int Ntotal = 100;

  float p=0.02;
  TH1F *data = new TH1F("data","",11,-0.5,10.5);
  TH1F *Poisson = new TH1F("Poisson","",11,-0.5,10.5);
  int nblue;
  int counter=0;
  for (int it=0;it!=Nthrows;it++) {
    nblue=0;
    for (int i=0;i!=Ntotal;i++)
      {      
	rnd=rand()/2147483647.;
	if (rnd<p) nblue++;
      }	
    //    cout << nblue << endl;
    data->Fill(nblue);
  }
  
  float prob=0.;
  for (int i=0; i!=11; i++)
    {
      Poisson->SetBinContent(i+1,poisson(p*Ntotal,i));
    }    

  gROOT->LoadMacro("AtlasStyle.C"); 

  gROOT->SetStyle("Plain");
  gStyle->SetCanvasColor(0);
  gStyle->SetCanvasBorderMode(0);
  gStyle->SetPadBorderMode(0);
  gStyle->SetPadColor(0);
  gStyle->SetLineWidth(1.);
  gStyle->SetCanvasDefH(596);
  gStyle->SetCanvasDefW(596);
  gStyle->SetCanvasBorderMode(0);
  gStyle->SetPadBorderMode(0);
  gStyle->SetFrameBorderMode(0);
  gStyle->SetTitleBorderSize(0);
  gStyle->SetStatBorderSize(0);
  gStyle->SetPadTickX(1);
  gStyle->SetPadTickY(1);
  gStyle->SetPalette(1);
  gStyle->SetOptStat(0);
  
  c0 = new TCanvas("c0", "", 0, 0, 600, 450);
  c0->SetLeftMargin(0.12);
  c0->SetRightMargin(0.12);
  c0->SetBottomMargin(0.2);
  
  data->GetYaxis()->SetTitleOffset(1.5);
  data->GetXaxis()->SetTitle("Number of Blue Balls");
  data->GetYaxis()->SetTitle("Throws");
  data->SetMarkerStyle(20);
  data->SetMarkerColor(2);
  data->SetMarkerSize(1.);
  data->Draw("ep");  

  Poisson->Scale(Nthrows);
  Poisson->SetMarkerStyle(20);
  Poisson->SetMarkerColor(3);
  Poisson->SetMarkerSize(1.);
  Poisson->Draw("same");  




  //  data->Fit("gaus");
  cout << "Number of Blue balls Mean/RMS/Error " <<  data->GetMean() << "   " << data->GetRMS() << "  " << data->GetRMS()/sqrt(Nthrows) << endl;
  cout << "Probability  " <<  100*data->GetMean()/Ntotal << " +/-  " << 100*data->GetRMS()/sqrt(Nthrows)/Ntotal << " % " << endl;
}