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

void Binomial()
{
  double rnd;
  int Nthrows=10000;
  int Ntotal = 300;

  float average=0.;
  float average_error=0.;

  int MeasCounts=0;
  float MeasFrequency=0.;
  float y[50],x[50],ey[50],ex[50];

  float p=0.15;
  TH1F *data = new TH1F("data","",81,-0.5,80.5);
  int nblue;
  data->Fill(36);
  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++;
      }	
    if (it == 0) {
      average=36;
    }
    else
      {
	average = ((it)*average+nblue)/(it+1);
      }
    // Log the average and its error every 100 throws
    if (it%10 == 0 && it < 500) {
      average_error = sqrt(average*(1-average/Ntotal))/sqrt(it+1);
      x[counter]=counter;
      ex[counter]=0.4;
      y[counter]=average;
      ey[counter]=average_error;
      cout << it << "  " << counter << "  " << average << "  +/-  " << average_error << endl;
      counter++;
    }

    //    cout << nblue << endl;
    data->Fill(nblue);
    if (nblue <= 36) MeasCounts++;
  }
  
  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("p");  

  c1 = new TCanvas("c1", "", 50, 50, 650, 500);
  c1->SetLeftMargin(0.12);
  c1->SetRightMargin(0.12);
  c1->SetBottomMargin(0.2);
  
  TMultiGraph *mg = new TMultiGraph();
  TGraphErrors *AverageGraph = new TGraphErrors(50,x,y,ex,ey);
  AverageGraph->SetMarkerStyle(20);
  AverageGraph->SetMarkerColor(2);
  AverageGraph->SetMarkerSize(1.);
  mg->Add(AverageGraph);
  mg->Draw("AP");



  //  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;
  MeasFrequency = (float) MeasCounts/Nthrows;
  cout << "Number of measurements bellow 36 : " << MeasCounts << "  Frequency  " << MeasFrequency << endl;
}