#if !defined(__CINT__) || defined(__MAKECINT__)
#include <Riostream.h>
#include <TF1.h>
#include <TFile.h>
#include <TH1F.h>
#include <TMath.h>

#include <TRandom3.h>
#include <TCanvas.h>
#endif

const double knorm = 1000.;  // number of experiments

// Function declaration 
double Decay(int n0,double alpha, double Delt,double timetot);  // Decay simulation 

double poissonian(double xx, double norm, double p); // POISSONIAN density 
double binomial(double xx, double norm, double ntrials, double prob);   // BINOMIAL density
void poisson(int n0 = 1000,double alpha = 2.5e-5, double Delt = 1.,double time=100.); //main function 

///////////////////////////////////////////////////////////////////////////
void poisson(int n0,double alpha, double Delt,double time){

  unsigned int seed;
  cout<<"Insert the seed for random numbers\n";
  cin>>seed;
  gRandom->SetSeed(seed);
  
  TString title = Form("Number of deacys in %4.1f s",time);
  int nbins = static_cast<int>(n0*alpha*time*4);
  cout<<"numer of bins"<<nbins<<endl;
  double high = nbins-0.5;
  TH1F *ndecays = new TH1F("ndecays",title,nbins,-0.5,high);

  //theoretical function (poissonian)
  double mu = alpha*n0*time;
  title =Form("Poissonian distribution with parameter %10.4g ",mu);
  
  TH1F *fteo = new TH1F("fteo",title,nbins,-0.5,high);
      
  //theoretical function 2(binomial)
  double ntrials = time/Delt;
  double prob = alpha*n0*Delt;
  title = Form("Binomial with parametris p= %10.4g e m = %8.2g",prob,ntrials);

  TH1F *fteo2 = new TH1F("fteo2",title,nbins,-0.5,high);
    
  cout<<"Numer of nuclei: "<<n0<<endl;
  cout<<"Parameter of poissonian distribution "<<mu<<endl;
  cout<<"Parameters of binomial distribution p= "<<prob<<" m ="<< ntrials<<endl;


  // Fill histo with theoretical functions
  for (double x=0.;x<high;x+=1.)
    fteo->Fill(x,poissonian(x,knorm,mu));
  for (double x=0.;x<high;x+=1.)
    //    fteo2->Fill(x,binomial(x,knorm,ntrials,prob)); //ramona
    fteo2->Fill(x,binomial(x,knorm,ntrials,prob)); //ramona
  // Simulation of the decay
  for(int exper = 0; exper<knorm; exper++){
    if((exper%100)==0)cout<<"Processinf event n. "<<exper<<endl;
    double nodecay=Decay(n0,alpha,Delt,time);
    ndecays->Fill(nodecay);
  }

  // SAVE histograms 
  
  TFile *file = new TFile("poisson.root","recreate");
  ndecays->Write();
  fteo->Write();
  fteo2->Write();
  file->Close(); 
  cout<<endl<<endl;

  fteo->SetLineColor(2);
  fteo2->SetLineColor(4);

  fteo->Draw("histc");
  new TCanvas();
  fteo2->Draw("histc");
  new TCanvas();
  ndecays->Draw("hist");
  fteo->Draw("histcsame");
  fteo2->Draw("histcsame");
  cout<<"Done - Histohrams are accessibile"<<endl;
  cout<<"in file poisson.root"<<endl;
}
///////////////////////////////////////////////////////////////////////////
double Decay(int n0,double alpha, double Delt,double timetot){
  double prob = alpha*Delt; //probabilita
  int n0init = n0; 
  for(double time=0.; time<timetot; time+=Delt){
    int ndec = 0;
    for(int nuclei=0; nuclei<n0;nuclei++)if(gRandom->Rndm()<prob)ndec++;
    n0-=ndec;
  }
  return static_cast<double>(n0init-n0);
}
///////////////////////////////////////////////////////////////////////////
double poissonian(double xx, double norm, double p){
  return norm*(TMath::PoissonI(xx,p));
  //  return norm*(TMath::Power(param,xx)*TMath::Exp(-param-Gamma(xx+1)));
}
///////////////////////////////////////////////////////////////////////////
double binomial(double xx, double norm, double ntrials, double prob)  {

  Double_t n = TMath::Binomial(ntrials,xx);
  n *= TMath::Power(prob,xx)*TMath::Power((1-prob),(ntrials-xx));
  n *= norm;
  return n;
}