/////////////////////////
// binomiale
/////////////////////////

#include "TStyle.h"
#include "TTimer.h"
#include "Riostream.h"
#include "TCanvas.h"
#include "TGraph.h"
#include "TROOT.h"
#include "TGraphErrors.h"
#include "TFrame.h"
#include "TAxis.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TRandom.h"
#include "TF1.h"
#include "TFormula.h"
#include "TFile.h"
#include "TTree.h"
#include "TMath.h"
#include "TMultiGraph.h"
#include "TPaveLabel.h"
#include "TLatex.h"
#include "TLegend.h"

Double_t constant (Double_t* x, Double_t* par) {

    return par[0];

}


Double_t BinomialCoeff( double N, double k){
  return ROOT::Math::tgamma(N+1) / (ROOT::Math::tgamma(k+1)*ROOT::Math::tgamma(N-k+1));
}

double BinomialFunc(Double_t *x, Double_t *par){
  double Ntot = par[0];
  double N = par[1];
  double p = par[2];

  return Ntot*BinomialCoeff(N,x[0]) * TMath::Power(p,x[0]) * TMath::Power(1-p,N-x[0]);
}

double PoissonFunc(Double_t *x, Double_t *par){
  double Ntot = par[0];
  double nu = par[1];

  return Ntot * TMath::Exp(-nu) * TMath::Power(nu,x[0]) / ROOT::Math::tgamma(x[0]+1);
}

int main(){

	using namespace std;

  // Parametri per la simulazione.
  int     Nrep  = 400;  // Numero di variabli casuali da generare.
  int     N     = 30;   // Il valore di N nella funzione di distribuzione binomiale.
  double  p     = 0.1;  // La probabilita' di osservare un esito favorevole della variabile casuale.
  int     k;
  double  numb;

  // Parametri per l'estetica dei grafici.
  int lf=133;  // fonts scales
  int tf=13;   // fonts label
  int cf=23;  // legend label
  int ls=20;   // size  scales
  int ts=25;   // size  label
  int cs=25;   // size  label

  double xOffset = 1.2;
  double yOffset = 1.0;
  double ms=1.2;
  int mStylepos=20, mStyleneg=22;

  TLatex *t2 = new TLatex();
  t2 = new TLatex();
  t2->SetNDC();
  t2->SetTextAlign(12);
  t2->SetTextFont(lf);
  t2->SetTextSize(ls);
  t2->SetTextColor(kBlack);

	// Istogramma da riempire con le variabili casuali generate.
  TH1F *hbino = new TH1F("hbino","",11,0,10);
  hbino->GetXaxis()->SetNdivisions(10);
  hbino->GetYaxis()->SetNdivisions(10);
  hbino->GetXaxis()->SetLabelFont(lf);
  hbino->GetYaxis()->SetLabelFont(lf);
  hbino->GetXaxis()->SetLabelSize(ls);
  hbino->GetYaxis()->SetLabelSize(ls);
  hbino->GetXaxis()->SetTitleFont(tf);
  hbino->GetYaxis()->SetTitleFont(tf);
  hbino->GetXaxis()->SetTitleSize(ts);
  hbino->GetYaxis()->SetTitleSize(ts);
  hbino->GetXaxis()->SetTitleOffset(1.5);
  hbino->GetYaxis()->SetTitleOffset(1.5);
  hbino->GetYaxis()->SetTitle( "counts" );
  hbino->GetXaxis()->SetTitle( "k" );
  hbino->SetFillColorAlpha(kOrange+7,0.4);
  hbino->SetLineColor(kBlack);
  hbino->SetLineWidth(2);

  // Generatore di numeri casuali.
  TRandom *rnd = new TRandom();

  // Istogramma per i numeri casuali > p.
  TH1F *hout = new TH1F("hout","",100,0,1.4);
  // Istogramma per i numeri casuali < p.
  TH1F *hin = new TH1F("hin","",100,0,1.4);
  hout->GetXaxis()->SetNdivisions(10);
  hout->GetYaxis()->SetNdivisions(10);
  hout->GetXaxis()->SetLabelFont(lf);
  hout->GetYaxis()->SetLabelFont(lf);
  hout->GetXaxis()->SetLabelSize(ls);
  hout->GetYaxis()->SetLabelSize(ls);
  hout->GetXaxis()->SetTitleFont(tf);
  hout->GetYaxis()->SetTitleFont(tf);
  hout->GetXaxis()->SetTitleSize(ts);
  hout->GetYaxis()->SetTitleSize(ts);
  hout->GetXaxis()->SetTitleOffset(1.5);
  hout->GetYaxis()->SetTitleOffset(1.5);
  hout->GetYaxis()->SetTitle( "counts" );
  hout->GetXaxis()->SetTitle( "rnd" );
  hout->SetStats(kFALSE);
  hout->SetFillColorAlpha(kRed,0.4);
  hin->SetFillColorAlpha(kGreen,0.4);

  // La funzione di distribuzione binomiale.
  TF1* fbino = new TF1("fbino",BinomialFunc,0.001,10,3);
  fbino -> SetParameters(Nrep,N,p);
  fbino->SetLineColor(kBlack);
  fbino->SetLineStyle(2);
  fbino->SetLineWidth(4);

  // La funzione di distribuzione di Poisson.
  TF1* fpois = new TF1("fpois",PoissonFunc,0.001,10,2);
  fpois -> SetParameters(Nrep,N*p);
  fpois->SetLineColor(kRed);
  fpois->SetLineStyle(4);
  fpois->SetLineWidth(4);

  // Funzioni costanti per il calcolo dei conteggi attesi nelle
  // distribuzioni dei numeri casuali con distribuzione uniforme
  // che sono minori e maggiori di p.
  TF1* fconst1 = new TF1("fconst1",constant,0.001,p,1);
  fconst1 -> SetParameter(0,Nrep);
  fconst1->SetLineColor(kGreen);
  fconst1->SetLineStyle(1);
  fconst1->SetLineWidth(2);

  TF1* fconst2 = new TF1("fconst2",constant,p,1,1);
  fconst2 -> SetParameter(0,Nrep);
  fconst2->SetLineColor(kRed);
  fconst2->SetLineStyle(1);
  fconst2->SetLineWidth(2);

  // Leggende associate agli istogrammi.
  TLegend *leg = new TLegend(0.6,0.50,0.9,0.7);
  leg->AddEntry(hbino,"num. events < p","f");
  leg->AddEntry(fbino,"expected binomial","l");
  leg->AddEntry(fpois,"expected poisson","l");
  leg->SetTextFont(12);
  leg->SetTextSize(0.03);
  leg->SetLineColor(1);
  leg->SetLineStyle(1);
  leg->SetLineWidth(0);
  leg->SetFillColor(0);
  leg->SetFillStyle(0);

  TLegend *leg2 = new TLegend(0.679,0.45,0.9,0.7);
  leg2->AddEntry(hin,"rnd < p","f");
  leg2->AddEntry(hout,"rnd > p","f");
  leg2->AddEntry(fconst1,"exp. unif. [0,p]","l");
  leg2->AddEntry(fconst2,"exp. unif. (p,1]","l");
  leg2->SetTextFont(12);
  leg2->SetTextSize(0.05);
  leg2->SetLineColor(1);
  leg2->SetLineStyle(1);
  leg2->SetLineWidth(0);
  leg2->SetFillColor(0);
  leg2->SetFillStyle(0);

  TLegend *leg2in = new TLegend(0.15,0.75,0.4,0.85);
  leg2in->SetTextFont(12);
  leg2in->SetTextSize(0.06);
  leg2in->SetLineColor(1);
  leg2in->SetLineStyle(1);
  leg2in->SetLineWidth(0);
  leg2in->SetFillColor(0);
  leg2in->SetFillStyle(0);

  TLegend *leg2out = new TLegend(0.5,0.75,0.8,0.85);
  leg2out->SetTextFont(12);
  leg2out->SetTextSize(0.06);
  leg2out->SetLineColor(1);
  leg2out->SetLineStyle(1);
  leg2out->SetLineWidth(0);
  leg2out->SetFillColor(0);
  leg2out->SetFillStyle(0);

  TLegend *leg2tot = new TLegend(0.32,0.60,0.45,0.70);
  leg2tot->SetTextFont(12);
  leg2tot->SetTextSize(0.06);
  leg2tot->SetLineColor(1);
  leg2tot->SetLineStyle(1);
  leg2tot->SetLineWidth(0);
  leg2tot->SetFillColor(0);
  leg2tot->SetFillStyle(0);


  TH1F *hout_run = (TH1F*)hout->Clone("hout_run");
  TH1F *hin_run = (TH1F*)hin->Clone("hin_run");
  hout_run->SetAxisRange(0,3,"Y");


  // Finestra per disegnare gli istogrammi dei numeri casuali
  // con distribuzione uniforme generati, e per confrontare
  // con i conteggi attesi.
  TCanvas * cside = new TCanvas("cside","uniform",600,800);
  cside->cd();
  gPad->SetBottomMargin(0.2);
  gPad->SetLeftMargin(0.2);
  cside->Divide(1,2);
  cside->cd(1);
  hout_run->Draw();
  hin_run->Draw("same");
  leg2->Draw("same");
  cside->cd(2);
  hout->Draw();
  hin->Draw("same");
  leg2->Draw("same");
  fconst1->Draw("same");
  fconst2->Draw("same");

  // Finestra per il grafico dell'istogramma delle variabili
  // casuali con funzione di distribuzione binomiale, e per il confronto
  // tra istogramma e conteggi attesi assumendo una funzione di distribuzione
  // binomiale e di Poisson.
  TCanvas *c1 = new TCanvas("c1","binomial",600,600);
  c1->cd();
  c1->Update();
  gPad->SetBottomMargin(0.2);
  gPad->SetLeftMargin(0.2);
  hbino->Draw();
  leg->Draw("same");
  double dxin = hin->GetBinWidth(10);
  double dxout = hin->GetBinWidth(10);
  for(int irep=0;irep<Nrep;irep++){ // Ciclo esterno per la generazione di Nrep variabili casuali binomiali.
    k = 0;
    hin_run->Reset();
    hout_run->Reset();
    for(int i=0; i<N; i++){   // Ciclo interno per generare N variabili con distribuzione uniforme in [0:1]. 
      numb = rnd -> Uniform(1.0); // Genero una variabile unifiorme in [0:1].
      if( numb < p) { // Conto quante variabili casuali sono minori di p.
          k += 1; // Il valore di k alla fine del conteggio, e' il valore della variabile casuale binomiale.
          hin->Fill(numb);
          hin_run->Fill(numb);
      } else {
        hout->Fill(numb);
        hout_run->Fill(numb);
      }
    }
    // Commandi per disegnare tutti i grafici dopo la generazione di ogni variabile casuale k.
    gSystem->ProcessEvents();
    c1->cd();
    c1->Update();
    hbino -> Fill(k);
    fbino -> SetParameters(hbino->GetEntries(),N,p);
    fpois -> SetParameters(hbino->GetEntries(),N*p);
    hbino->Draw("HIST,E");
    fbino->Draw("same");
    fpois->Draw("same");
    leg->Draw("same");
    cside->cd(1);
    cside->Update();
    leg2in->Clear();
    leg2in->AddEntry(hin_run,Form("rnd < p: %d",int(hin_run->GetEntries())),"f");
    leg2out->Clear();
    leg2out->AddEntry(hout_run,Form("rnd > p: %d",int(hout_run->GetEntries())),"f");
    leg2tot->Clear();
    leg2tot->AddEntry(hin_run,Form("total ( = N ): %d",int(hin_run->GetEntries())+int(hout_run->GetEntries())),"");
    leg2tot->AddEntry(hin_run,Form("p = %f",p),"");
    hout_run->Draw();
    hin_run->Draw("same");
    leg2in->Draw("same");
    leg2out->Draw("same");
    leg2tot->Draw("same");
    cside->cd(2);
    cside->Update();
    hout->Draw();
    hin->Draw("same");
    fconst1 -> SetParameter(0,hin->GetEntries()*dxin/p);
    fconst2 -> SetParameter(0,hout->GetEntries()*dxout/(1-p));
    fconst1->Draw("same");
    fconst2->Draw("same");
    leg2->Draw("same");
  }

  return 0;
}