ROOT Tutorial
Parte 3
Generatori di numeri casuali
In ROOT numeri pseudo - casuali sono generati usando le classi
TRandom. Esistono 4 differenti classi (TRandom, TRandom1,
TRandom2 e TRandom3), ciascuna delle quali implementa un differente
tipo di generatore di numeri casuali.
Il metodo SetSeed permette di scegliere il seme.
gRandom->SetSeed(Int_t seed);
Se non specificato, il seme è stabilito dal clock del PC
Esempio: distribuzione uniforme
void randomUni(){
Double_t x, y;
TH2D *h2 = new TH2D(“h2”,”Uniform Distribution”, 100,0,10,100,0,10);
for(Int_t i=0; i < 100000; i++){
x = gRandom->Uniform(0,10);
y = gRandom->Uniform(0,10);
h2->Fill(x,y);
}
gStyle->SetPalette(1);
h2->Draw(“colz”);
}
Distribuzioni non uniformi
Tre possibilità (vedi lezione MC prof. Riggi):
a)
b)
c)
Distribuzione uniforme con eventi pesati
Metodo dell’inversione
Metodo del rigetto
Esempio: distribuzione eventi secondo sin
Proviamo:

Metodo dell’inversione
 Metodo del rigetto
Metodo dell’inversione
Problema: Estrarre una sequenza di numeri casuali distribuiti secondo
la funzione f(x) = sin(x)
Data la funzione di probabilità f(x), la funzione integrale (cumulativa)
F(t) è data da:
F(t) = ∫f(t)dt = ∫sin(t)dt = -cos(t)
Se la funzione integrale è invertibile, si può generare la sequenza di
numeri casuali x, distribuiti secondo f(x), partendo da un numero
casuale u (distribuito tra -1 e 1) e calcolando x da:
x = F-1 (u) = - arccos(u)
(se  < x < 2)
x = F-1 (u) = arccos(u)
(se 0 < x < )
Oppure
Metodo dell’inversione
void inversion()
{
Int_t nev = 1000000;
//sin(theta) distribution
TH1D *h1 = new TH1D("h1","sin(theta) distribution", 100,0,TMath::Pi());
TStopwatch t;
t.Start();
Double_t ux;
for(Int_t n = 0; n < nev; n++) {
ux = TMath::ACos(gRandom->Uniform(-1,1));
h1->Fill(ux);
}
t.Stop();
t.Print("m");
h1->Draw();
}
Metodo del rigetto
Problema: Estrarre una sequenza di numeri casuali distribuiti secondo
la funzione f(x) = sin(x) nell’intervallo (0, )
Metodo:
a.
b.
c.
d.
Generare ux uniforme in (0, )
Generare uy uniforme in (0,fmax) dove fmax = 1
Valutare f(ux)
Confrontare uy con f(ux)
Se uy < f(ux) accettare ux
Se uy > f(ux) ripetere da a. in poi
Metodo del rigetto
void rejection(){
Int_t nev = 1000000;
Int_t n=0;
//sin(theta) distribution
TH1D *h1 = new TH1D("h1","sin(theta) distribution", 100,0,TMath::Pi());
TStopwatch t;
t.Start();
Double_t ux; // between 0 and Pi/2
Double_t f_ux;
Double_t uy; // between -1 and 1
do {
ux = gRandom->Uniform(0, TMath::Pi());
f_ux = TMath::Sin(ux);
uy = gRandom->Uniform(0, 1);
if(uy > f_ux) continue; //rejection
n++;
h1->Fill(ux);
}
while (n < nev);
t.Stop();
t.Print("m");
h1->Draw();
}
Distribuzione non uniformi disponibili in ROOT
Esercizio: simulazione MC per valutare la distanza relativa tra due
contatori geiger
Ipotesi:
1) Contatori Geiger bidimensionali
r = 0.5 Geiger radius in cm
l = 6 Geiger length in cm
2) Distanza nominale d = 5 cm
3) Estrazione uniforme dei punti
Possibili miglioramenti:
1) Contatori Geiger in 3D
DA VERIFICARE!
void geigerMC(){
Double_t x1, y1, x2, y2;
//Define the experimental setup
const int d = 5; // Distance between geigers
Double_t r = 0.5; // Geiger radius in cm
Double_t l = 6; // Geiger length in cm
Double_t distance;
TH1D *h1 = new TH1D("h1","Distances distribution", 100,0,10);
TH2D *h2 = new TH2D("h2","Points distribution", 150, -5, 10, 100, -5, 5);
for(Int_t i=0; i < 100000; i++){
x1 = gRandom->Uniform(-r,r);
y1 = gRandom->Uniform(-l/2,l/2);
x2 = gRandom->Uniform(d-r,d+r);
y2 = gRandom->Uniform(-l/2,l/2);
h2->Fill(x1,y1);
h2->Fill(x2,y2);
distance = TMath::Sqrt((x1-x2)**2+(y1-y2)**2);
h1->Fill(distance);
}
TCanvas *c1 = new TCanvas();
h2->Draw("colz");
TCanvas *c2 = new TCanvas();
h1->Draw();
}
5 cm
3 cm
15 cm
Esempio: distribuzione uniforme in x, gaussiana in y
void randomGaus(){
Double_t x, y;
TH2D *h2 = new TH2D(“h2”,”Uniform Distribution”, 100,0,10,100,0,10);
for(Int_t i=0; i < 100000; i++){
x = gRandom->Uniform(0,10);
y = gRandom->Gaus(5,1);
h2->Fill(x,y);
}
gStyle->SetPalette(1);
h2->Draw();
}
h2->Draw(“colz”);
h2->Draw(“lego”);
Esempio: distribuzione uniforme nello spazio
Estrazione uniforme (,)?
void sfera(){
Int_t raggio=1;
Double_t teta, phi;
TH3D *hsemisfera = new TH3D("hsemisfera","hsemisfera",200,1,1,200,-1,1,100,0,1);
for(Int_t i=0; i<10000; i++){
teta = gRandom->Uniform(0,TMath::Pi()/2); //semisfera
// teta = gRandom->Uniform(0,TMath::Pi()); // sfera
phi = gRandom->Uniform(0,2*TMath::Pi());
Double_t x = raggio * TMath::Cos(phi)*TMath::Sin(teta);
Double_t y = raggio * TMath::Sin(phi)*TMath::Sin(teta);
Double_t z = raggio * TMath::Cos(teta);
hsemisfera->Fill(x,y,z);
}
hsemisfera->Draw();
}
Estrazione uniforme in  e in sin()
void sfera2(){
Int_t raggio=1;
Double_t teta, phi;
TH3D *hsemisfera = new TH3D("hsemisfera","hsemisfera",200,1,1,200,-1,1,100,0,1);
for(Int_t i=0; i<10000; i++){
teta = TMath::ACos(gRandom->Uniform(0,1)); //Semisfera
//teta = TMath::ACos(-1 + 2*gRandom->Uniform(0,1)); //Sfera
phi = gRandom->Uniform(0,2*TMath::Pi());
Double_t x = raggio * TMath::Cos(phi)*TMath::Sin(teta);
Double_t y = raggio * TMath::Sin(phi)*TMath::Sin(teta);
Double_t z = raggio * TMath::Cos(teta);
hsemisfera->Fill(x,y,z);
}
hsemisfera->Draw();
}
Estrarre numeri casuali secondo una distribuzione f(x) nota
TF1 *f = new TF1(“f”,”TMath::Gaus(x,0,1)”,-5,5);
TH1D *h1 = new TH1D("h1","",100,-5,5);
for(int i=0; i<1000; i++) h1->Fill(f->GetRandom());
h1->Draw();
ROOT Files
Si tratta di file salvati in un formato leggibile da ROOT. Può contenere
vari tipi di oggetti.
Esempio:
TFile f(“esempio.root","recreate");
h0->Write();
f.Close();
Per creare un file ROOT
e salvare un istogramma
Opzioni: NEW, CREATE, RECREATE, UPDATE, READ
TFile f("demo.root");
TBrowser browser;
Per aprire un file ROOT
e vederne il contenuto