Lezione 4

Introduzione

Scopo di questa sessione di laboratorio è prendere confidenza con il trattamento di numeri di numeri pseudocasuali. Questi saranno usati estensivamente nella lezione successiva per lo sviluppo di una simulazione Monte Carlo.

Per fare questi esercizi utilizzeremo principalmente gli oggetti TRandom di ROOT.

Classe TRandom

La classe TRandom permette di costruire oggetti in grado di generare numeri pseudocasuali.

Il generatore di numeri casuali che sta al suo interno di TRandom utilizza un metodo di congruenza lineare, non particolarmente efficace. Per verificarlo, scaricate e fate eseguire la macro
correlazione.C:
questa macro usa un oggetto TRandom per generare numeri uniformemente nell'intervallo (0,1) e crea un istogramma di x_N in funzione di x_N-1

Si può osservare come i punti generati non riempiano il piano, ma si dispongono su un reticolo, a causa della bassa periodicità del generatore. Per vedere l'effetto bisogna ingrandire molto una regione dello spazio che contiene i punti, ma allora solo una piccola frazione di questi punti cadrà nella regione selezionata ed il tempo di esecuzione sarà quindi molto grande, perché bisogna generare molti più punti di quelli utili.

L'esecuzione della macro dovrebbe dare un grafico simile al seguente, anche se con meno statistica:

ma richiede parecchio tempo di CPU.

Il generatore TRandom è adeguato per questi esercizi, ma per ogni applicazione più seria, è meglio utilizzare una delle altre classi messe a disposizione da RooT con generatori di qualità maggiore (TRandom2 e TRandom3), che hanno la stessa interfaccia utente.

Metodi degli oggetti TRandom

TRandom(unsigned int seed); (costruttore)

costruisce un oggetto TRandom, inizializzando il seme del generatore di numeri casuali a seed.

Se seed=0, allora il seme viene generato a partire dall'orologio del sistema

(con granularita` di un secondo, per cui fare attenzione: se si dichiarano insieme diversi oggetti TRandom,

molto probabilmente questi genereranno la stessa sequenza di numeri).

void SetSeed(unsigned int seed);

definisce il seme del generatore di numeri casuali, se invocato con argomento 0,

lo costruisce a partire dall'orologio di sistema.

double Rndm();

restituisce un numero distribuito uniformemente nell'intervallo [0,1]

int Poisson(double m);

restituisce un numero intero secondo la distribuzione di Poisson con valor medio m

double Gaus(double media, double sigma);

restituisce un numero secondo una distribuzione gaussiana con i valori dati di valor medio e sigma.

double Uniform(double min, double max);

restituisce un numero secondo una distribuzione uniforme nell'intervallo corrispondente ai valori dati.

double Exp(double tau);

restituisce un numero secondo una distribuzione esponenziale (1/τ)exp(-x/τ) con il valore dato del parametro τ.

double BreitWigner(double media, double gamma);

restituisce un numero secondo una distribuzione di Breit-Wigner (tipica dei fenomeni di risonanza: (2/π) Γ/((x-x_medio)²+Γ²)

con i valori dati dei parametri x_medio e Γ.

Utilizzo degli oggetti TRandom

Per definire un oggetto di tipo TRandom, basta dichiararlo:
TRandom nomeOggetto;
eventualmente utilizzando un costruttore per inizializzare il seme del generatore:
TRandom nomeOggetto(0);.

Per usare i metodi ad esso associati, si usa la solita sintassi degli oggetti:
nomeOggetto.nomeMetodo(argomenti);
Per fare un esempio, per stampare 100 numeri interi estratti da una distribuzione di Poisson con valor medio 3, potremmo scrivere il seguente frammento di codice in RooT:


TRandom gen(0);
for (int i=0; i<100; i++) cout << gen.Poisson(3.) << endl;

È possibile anche costruire puntatori ad oggetti, utilizzanto l'operatore new abbinato ad un costruttore:
TRandom *gen = new TRandom(0);
ma a questo punto per accedere ai metodi, bisogna usare l'operatore freccia:
gen->Poisson(3.)

Potete trovare una tabella riassuntiva con la selezione dei tipi di oggetti e dei relativi metodi che saranno utili per questa lezione e per la parte di simulazione.

Oltre alla tabella riassuntiva, potete trovare:

referenza completa delle classi di ROOT
Guida dell'utente
Lezioni di Alfio Lazzaro per il corso di analisi dati in fisica (prof. Palombo)

Teorema del limite centrale

Il teorema del limite centrale, detto anche legge dei grandi numeri è la base di tutto il trattamento statistico degli errori e grosso modo puo' essere formulato:
Siano x₁, x₂, x₃, ...x_N numeri casuali tutti estratti da una stessa distribuzione con valor medio X_m e varianza σ². Allora, per N tendente all'infinito la quantita' (x₁+x₂+...+x_N)/N (la media degli N numeri) tendera' ad avere una distribuzione gaussiana con valor medio X_m e varianza σ²/N.

Il teorema di fatto dice che, se si fanno tante misure, anche con errori non gaussiani, la media alla fine avra' una distribuzione gaussiana ed un'incertezza (r.m.s. ovvero la radice quadrata della varianza) che scala con la radice quadrata del numero di misure mediate.

Esercizio

Verificare empiricamente il teorema del limite centrale:

costruire un programma che usa le librerie di ROOT che scriva su di un file dei numeri distribuiti uniformemente tra 0 e 1.
costruire un programma per leggere il file e mettere i valori ottenuti in un istogramma, verificare che la distribuzione sia corretta, ben contenuta nell'istogramma e segnare il suo valor medio ed r.m.s.
leggere i numeri del file facendone le medie a due a due, fare un istogramma delle medie.
ripetere il punto 3 facendo le media a gruppi di 4 e 8: verificare che la distribuzione diventi sempre più gaussiana come forma e che l'rms scali nel modo indicato sopra.

N.B.: potrebbe essere utile strutturare il programma di lettura per i punti da 2 a 3 in modo che prenda come input il numero di misure da mediare (da linea di comando o da terminale) ed effettui il ciclo di lettura di conseguenza.

Integrazione con metodi Monte Carlo

Una forma di integrazione possibile è quella di estrarre delle numeri a caso in un area che contiene il grafico della funzione da integrale. La frazione dei punti che cadono all'interno del grafico della funzione moltiplicata per l'area totale su cui si estraggono i punti, da` il valore dell'integrale.

Questo metodo puo` essere utile nel caso di funzioni rapidamente oscillanti o nel caso di integrali multidimensionali (trasparenze).

Una funzione oscillante

Come esercizio, calcolare l'integrale tra 0 e 2 della funzione sin²(1/x) con i seguenti metodi:

Metodo hit or miss
estraendo punti a caso nel quadrato (0,2)x(0,1), calcolare la frazione di punti che cadono all'interno del grafico della funzione (hits) e stimare l'integrale usando la relazione
I=A_rettangoloN_hits/N_punti.
Metodo del valor medio
scegliendo a caso punti nell'intervallo di validità della funzione (0,2), e calcolando il valor medio della funzione integranda f(x) sui punti estratti, approssimare l'integrale con
I=A_intervallo<f(x)>

In entrambi di casi, farsi dare dal programma il valore dell'integrale aggiornato ogni 100 estrazioni: facendo un grafico di questo numero dovreste vedere un fenomeno del tipo di quello presentato nella figura di seguito, in cui, dopo alcune oscillazioni iniziali, l'integrale converge:

Integrali multidimensionali

Proviamo ad applicare gli stessi principi al calcolo di volumi di oggetti multi-dimensionali. Provare a scrivere una macro di ROOT che calcoli il volume di una sfera n-dimensionale di raggio 1 con il metodo hit or miss: fate un grafico del valore ottenuto in funzione del numero di punti estratti, come nel caso dell'integrale precedente.

Per comodita` sarebbe meglio fare il grafico della differenza tra il valore calcolato del volume della sfera e quello vero:

Si noti che in ROOT si puo` accedere ai valori di π ed alla funzione Γ attraverso le chiamate alle funzioni:
double TMath::Pi(); double TMath::Gamma(double);

Andando su qualcosa di più complicato, la macro Potenziale.C calcola il potenziale gravitazionale della terra, assumendo che sia una sfera di densità uniforme, effettuando con il metodo del valor medio l'integrale:

Costruire un grafico del potenziale gravitazionale dal centro della sfera, fino ad una distanza pari a 3 volte il raggio della sfera:

Esercizi opzionali

Test di `TRandom2` e `TRandom3`

Come detto in precedenza, RooT fornisce delle classi TRandom2 e TRandom3, che hanno le stesse proprieta` degli oggetti TRandom, salvo dei generatori che adottano particolari accorgimenti per ridurre le correlazioni: se siete riusciti a costruire il programma compilato, provate a farlo rigirare usando TRandom2 o TRandom3 al posto del semplice TRandom.

Ulteriori veriche del limite centrale

Provare a ripetere l'esercizio del limite centrale, con altre distribuzioni di partenza:

un'esponenziale con τ=0.5
una Breit-Wigner con x_medio=0.5 e Γ=0.2.

Nel caso della Breit-Wigner si vede che il teorema del limite centrale non funziona e la r.m.s. non dimisuisce facendo la media. Quale delle ipotesi del teorema non è verificata per questa distribuzione?

Mediana

Per distribuzioni come la Breit-Wigner, abbiamo visto che le condizione del teorema del limite centrale non sono soddisfatte. Per avere una stima del punto centrale della distribuzione, aziché il valor medio, conviene prendere la mediana, definita con la seguente procedura:

si leggono N punti;
si mettono in ordine (sfruttare la funzione ordina della prima lezione)
se N e' dispari, la mediana e' il punto centrale, di indice (N-1)/2;
se N e' pari, la mediana e' il punto di mezzo tra i punti di indice N/2-1 e N/2.

Si dovrebbe osservare che, mentre la varianza della media non cambia, quella della mediana scende rapidamente con il numero di misure.

Potenziale generato da una crosta sferica

Provare ad adattare la funzione che effettua l'integrale multidimensionale del potenziale gravitazionale generato da una massa sferica al caso in cui la massa fosse in realtà concentrata in una crosta sferica di raggio 1000 km.

In questo caso è forse più utile passare direttamente in coordinate polari, generando numeri uniformemente negli intervalli:

R_Terra-1000 km < r < R_Terra
0 < θ < π
0 < φ < 2π

e ricordandosi che la funzione integranda in coordinate polari cambia forma e diventa:

Laboratorio di Calcolo 2