// ----------------------------------------------------------------------------------- //
/*
  ROOT macro for generating Binomial and Poisson distributions and testing their
  differences and likenesses.

  For more information, see:
    Barlow: 3.2 + 3.3
    http://en.wikipedia.org/wiki/Binomial_distribution
    http://en.wikipedia.org/wiki/Poisson_distribution

  Author: Troels C. Petersen (NBI/CERN)
  Email:  Troels.Petersen@cern.ch
  Date:   12th of September 2010
*/
// ----------------------------------------------------------------------------------- //

#include <TFile.h>
#include <TH1F.h>
#include <TRandom3.h>


//---------------------------------------------------------------------------------- 
double sqr(double a) {
  return a*a;
}


//---------------------------------------------------------------------------------- 
void BinomialPoisson() 
//---------------------------------------------------------------------------------- 
{
  // Setting of general plotting style:
  gStyle->SetCanvasColor(0);
  gStyle->SetFillColor(0);
  // Setting what to be shown in statistics box:
  gStyle->SetOptStat("e");
  gStyle->SetOptFit(1111);

  // Set parameters:
  int Nexperiments = 1000;
  int Ntrials =  10;
  double prob = 0.3;

  TRandom3 r;

  // Make histograms:
  TH1F *Hist_Bino = new TH1F( "Hist_Bino", "Binomial Distribution", 10, -0.5, 9.5);
  TH1F *Hist_Pois = new TH1F( "Hist_Pois", "Poisson Distribution",  10, -0.5, 9.5);


//----------------------------------------------------------------------------------
// Loop over process:
//----------------------------------------------------------------------------------

  for (int i=0; i < Nexperiments; i++) {
    // Get random number from Binomial distributed:
    double x1 = r.Binomial(Ntrials,prob);
    Hist_Bino->Fill(x1);

    // Get random number from Poisson distributed:
    double lambda = double(Ntrials) * prob;
    double x2 = r.Poisson(lambda);
    Hist_Pois->Fill(x2);

    // Print a few entries (just to check output and see that the program runs!):
    if (i < 5) printf("  Random Binomial and Poisson: %3d  %3d \n", x1, x2);
  }


//---------------------------------------------------------------------------------- 
// Plot histograms on screen:
//---------------------------------------------------------------------------------- 

  // Making a new window (canvas):
  // --------------------------------------------------------------------
  TCanvas* c0 = new TCanvas("c0", "", 100, 20, 600, 400);
  // c0->SetLogy();

  // Setting line width and color and axis titles of histograms:
  Hist_Bino->GetXaxis()->SetTitle("Random number (Binomial/Poisson)");
  Hist_Bino->GetYaxis()->SetTitle("Frequency");
  Hist_Bino->SetLineWidth(2);
  Hist_Bino->SetLineColor(2);
  Hist_Pois->SetLineWidth(2);
  Hist_Pois->SetLineColor(4);

  // Fitting histogram with constant:
  TF1 *fitPois = new TF1("fitPois","[0] * TMath::Exp(-[1]) * TMath::Power([1],x) / TMath::Gamma(x+1.0)", -0.5, 9.5);
  fitPois->SetParameters(double(Nexperiments), lambda);
  fitPois->SetLineStyle(3);
  // Hist_Bino->Fit("fitPois", "RL");

  // Drawing histograms with errors in the same plot:
  Hist_Bino->Draw("e");
  Hist_Pois->Draw("e same");


  // Legend:
  TLegend* legend = new TLegend(0.20, 0.18, 0.45, 0.33);
  legend->SetLineColor(1);
  legend->SetFillColor(0);
  legend->AddEntry(Hist_Bino, "Binomial", "LE");
  legend->AddEntry(Hist_Pois, "Poisson", "LE");
  legend->Draw();


  c0->Update();
  // c0->SaveAs("Dist_BinoPois.png");




//---------------------------------------------------------------------------------- 
// Write histograms to file:
//---------------------------------------------------------------------------------- 
  /*
  TFile *file = new TFile("results.root","RECREATE","Histograms");

  Hist_Bino->Write();
  Hist_Pois->Write();
  
  file->Write();
  file->Close();
  */
}



//---------------------------------------------------------------------------------- 
/*

Yes, once again, read through the program, and see if you understand what is going on!!!

Questions:
----------
 1) With the initial settings, does the Binomial distribution approach the Poisson?
    Think about how you quantify this. Do you include the errors?
    To decrease the errors, try increasing the number of experiments to 10000. How much
    smaller should they then be in average?

 2) Try to increase the number of trials by a factor 10, but decrease the probability
    of succes with a factor 10, keeping lambda constant. How good is the correspondance
    now? Try to repeat this for other values. When do you get a reasonably good
    correspondance (i.e. when can you use a Poisson instead of a Binomial)?
    Note that you might want to change the histograms, and perhaps also compare them on
    a logarithmic scale (uncomment line 77).

 3) Try to fit the Poisson distribution with a Gaussian distribution (copy this from
    e.g. "TableMeasurements"). How large does lambda need to be for the fit to become
    "reasonable"?


Advanced questions:
-------------------
 1) Look up the "Chi2Test" method for 1D histograms, and apply it to the two histograms
    in the problem (using 1000 experiments). For what values does the probability for
    the two to be the same become "reasonable" (say above 1%)?

*/
//----------------------------------------------------------------------------------