#if !defined(__CINT__) || defined(__MAKECINT__)
#include <TF1.h>
#include <TFile.h>
#include <TH1F.h>
#include <TMath.h>
#include <TRandom3.h>
#include <TROOT.h>
#endif
void decay(Int_t n0 = 100,Float_t alpha = 0.01, Float_t Delt = 1.);
Double_t exponential(Double_t *x, Double_t *par);

void decay(Int_t n0,Float_t alpha, Float_t Delt){
  delete gRandom;
  gRandom = new TRandom3();
  Int_t Nbins = 300; // numero di intervalli di tempo
  Float_t timetot = Delt*Nbins;  // tempo totale
  // histogram booking
  TH1F *h1 = new TH1F("h1","Remaining nuclei",Nbins,0.,timetot);
  //funzione teorica
  TF1 *fteo = new TF1("fteo",exponential,0.,timetot,2);
  //  fteo->SetParameters(static_cast<Double_t>(n0),static_cast<Double_t>(alpha));
  Double_t N0 = static_cast<Double_t>(n0);
  Double_t ALPHA = static_cast<Double_t>(alpha);
  fteo->SetParameters(N0,ALPHA);
  fteo->SetParNames("normalizzazione","coefficiente");
  Float_t prob = alpha*Delt; //probabilita
 
  for(Float_t time=0.; time<timetot; time+=Delt){
    for(Int_t nuclei=0; nuclei<n0;nuclei++)if(gRandom->Rndm()<prob)n0--;
    h1->Fill(time,static_cast<Float_t>(n0));
  }
  TFile *file = new TFile("decay.root","recreate");
  h1->Write();
  fteo->Write();
  file->Close(); 
}

Double_t exponential(Double_t *x, Double_t *par){
  Double_t xx = x[0];
  return par[0]*exp(-par[1]*xx);
}