#if !defined(__CINT__) || defined(__MAKECINT__)
#include <stdio.h>
#include <iostream>
#include <TF1.h>
#include <TFile.h>
#include <TH1F.h>
#include <TMath.h>
#include <TRandom3.h>
#include <TROOT.h>
#endif

// DICHIARAZIONE DELLE FUNZIONI
Float_t gamma(Float_t xx);   // FUNZIONE GAMMA (FATTORIALE)
Float_t decay(Int_t n0,Float_t alpha, Float_t Delt,Float_t timetot);  // SIM. DECADIMENTI
Float_t poissonian(Double_t xx, Double_t norm, Double_t param); // DENSITA POISSONIANA
Float_t binomial(Double_t xx, Double_t norm, Double_t prob, Double_t ntrials); //DENSITA BINOMIALE
void poisson(Int_t n0 = 500,Float_t alpha = 4.e-5, Float_t Delt = 10.,Float_t time=100.);

///////////////////////////////////////////////////////////////////////////
void poisson(Int_t n0,Float_t alpha, Float_t Delt,Float_t time){
  // gRandom punta ad un oggetto della classe TRandom, instanziato per
  // default all'avvio di root.
  // Questo oggetto viene qui eliminato e sostituito con un oggetto
  // della classe TRandom3 (generatore  Mersenne-Twister)
  delete gRandom;  
  gRandom = new TRandom3();

  // booking degli istogrammi
  char title[50];
  sprintf(title,"Numero di decadimenti in %4.1f s",time);
  Int_t nbins = static_cast<Int_t>(n0*0.05);
  Float_t high = nbins-0.5;
  TH1F *ndecays = new TH1F("ndecays",title,nbins,-0.5,high);
  //funzione teorica (poissoniana)
  Double_t nu = static_cast<Double_t>(alpha*n0*time);
  sprintf(title,"Poissoniana con parametro %10.4g ",nu);
  TH1F *fteo = new TH1F("fteo",title,nbins,-0.5,high);
  //funzione teorica2 (binomiale)
  Double_t ntrials = static_cast<Double_t>(time/Delt);
  Double_t prob = static_cast<Double_t>(alpha*n0*Delt);
  sprintf(title,"Binomiale con parametri p= %10.4g e m = %8.2g",prob,ntrials);
  cout<<title<<endl;
  TH1F *fteo2 = new TH1F("fteo2",title,nbins,-0.5,high);

  cout<<"Numero nuclei: "<<n0<<endl;
  cout<<"Parametro della poissoniana "<<nu<<endl;
  cout<<"Limite superiore per gli istogrammi "<<high<<endl;

  Double_t norm = 1000.;  // numero di esperimenti (--> normalizzazione curva teorica)
  for (Float_t x=0.;x<high;x+=1.)fteo->Fill(x,poissonian(x,norm,nu));
  for (Float_t x=0.;x<high;x+=1.)fteo2->Fill(x,binomial(x,norm,prob,ntrials));

  for(Int_t exper = 0; exper<norm; exper++){
    Float_t nodecay=decay(n0,alpha,Delt,time);
    ndecays->Fill(nodecay);
  }

  // SALVATAGGIO DEGLI ISTOGRAMMI
  TFile *file = new TFile("poisson.root","recreate");
  ndecays->Write();
  fteo->Write();
  fteo2->Write();
  file->Close(); 
}
///////////////////////////////////////////////////////////////////////////
Float_t decay(Int_t n0,Float_t alpha, Float_t Delt,Float_t timetot){
  Float_t prob = alpha*Delt; //probabilita
  Int_t n0init = n0; 
  for(Float_t time=0.; time<timetot; time+=Delt){
    for(Int_t nuclei=0; nuclei<n0;nuclei++)if(gRandom->Rndm()<prob)n0--;
  }
  return static_cast<Float_t>(n0init-n0);
}
///////////////////////////////////////////////////////////////////////////
Float_t poissonian(Double_t xx, Double_t norm, Double_t param){
  Double_t fact = TMath::Exp(gamma(xx+1));
  return norm*(TMath::Power(param,xx)*TMath::Exp(-param)/fact);
}
///////////////////////////////////////////////////////////////////////////
Float_t binomial(Double_t xx, Double_t norm, Double_t prob, Double_t ntrials){
  Double_t log1=gamma(ntrials+1);
  Double_t log2=gamma(xx+1);
  Double_t log3=gamma(ntrials-xx+1);
  Double_t bcoeff=TMath::Exp(log1-log2-log3);
  return norm*TMath::Power(prob,xx)*TMath::Power(1-prob,ntrials-xx)*bcoeff;
}
///////////////////////////////////////////////////////////////////////////
Float_t gamma(Float_t xx){
  // returns log(gamma(xx)) - tratta da  "Numerical Recipes"

  if(xx<1){
    cout<<"function gamma: argomento non valido "<<xx<<endl;
    return -1.;
  }
  double x,y,tmp,ser;
  static double cof[6]={76.18009172947146,-86.50532032941677,
			24.01409824083091,-1.231739572450155,
			0.1208650973866179e-2,-0.5395239384953e-5};
  int j;

  y=x=xx;
  tmp=x+5.5;
  tmp -= (x+0.5)*TMath::Log(tmp);
  ser=1.000000000190015;
  for (j=0;j<=5;j++) ser += cof[j]/++y;
  return -tmp+log(2.5066282746310005*ser/x);

}