//
// esempio di utilizzo del generatore di numeri random in root
//
// Includo le librerire degli "oggetti" che  servono
#include "TRandom.h"
#include "TH1D.h"
#include "TFile.h"
//
//  Calcolo integrale con il metodo di reiezione
//
// calcoliamo l'integrale di questa funzione
double myf(double x) {
     return 2.*x*x+1.;
}
void Reiezione(){
// creo un "oggetto" generatore di numeri random
   TRandom R(56335);

   // con N=10^6 RMS ~ 0.88 
   // con N=10^4 RMS ~ 8.8 
   // con N=10^2 RMS ~ 88. 
   
   //int N=1000000/100;
   int N=1000000;
   int M = 500;
// creo un oggetto file lo collego al file di nome Reiezione.roots
// e gli dico di crearlo (ricrearlo) quindi è un output
   TFile f("Integrale_Reiezione.root","RECREATE");
   TH2D * h1 = new TH2D ("h1","Tutti i Punti Estratti",80,0.,11.,80,0.,210);
   TH2D * h2 = new TH2D ("h2","Punti Contenuti",80,0.,11.,80,0.,210);
   //double dx = 5.*10;
   double dx = 5.;
   TH1D * harea = new TH1D ("harea","Stima dell'area -  metodo della Reiezione",100,675.-dx,675.+dx);




   for (int j=0; j<M; ++j) {
     // inizializzo sequenza di conteggio
     int In=0;
     for (int i=0; i<N ; i++) {
       double x,y;
       // Uso l'oggetto generatore di numeri random per generare numeri random in maniera uniforme
       x=R.Uniform()*9+1;
       y=R.Uniform()*201;
       
       // Riempio l'oggetto histogramma con tutti numeri random generati
       if (j==1) h1->Fill(x,y);
      // Riempio l'oggetto histogramma con solo i punti contenuti
      if (y<myf(x)) {
	if (j==1) h2->Fill(x,y);
	In++;
      }
     }
// Calcolo area
     double area=9.*201.*In/N;
     double succRate = double(In)/double(N);
     double errSuccRate = sqrt(succRate*(1.-succRate)/(double)N);
     double errore = 9.*201.*errSuccRate;
     printf("\nSequenza n %i di %i  ---- elementi nella sequenza %i \n", j, M, N);
     printf("Ho estratto n=%i valori utili su un totale di %i  --- efficienza = %g +/- %g \n",In, N, succRate, errSuccRate);
     printf("Integrale = %lf\n",area);
     printf("Errore su Integrale = %g  \n",errore);
     harea->Fill(area);
   }
// scrivo gli oggetti histogramma nel file
   h1->Write();
   h2->Write();
   harea->Write();
// chiudo il file attraverso il metodo dell'oggetto di tipo file.
   f.Close();
   printf("Istogrammi salvati nel file %s \n",f.GetName());
}