ma richiede parecchio tempo di CPU.
Vedremo come utilizzare tali oggetti all'interno di un programma compilato.
Il generatore di numeri casuali che sta all'interno di TRandom utilizza il classico metodo della congruenza visto a lezione. Se vi ricordate questo metodo ha il problema che, se si guarda in dettaglio, esiste una correlazione tra un numero ed il suo successivo, che si dispongono in realta` su una serie di rette parallele.
Mentre continua la spiegazione, scricate 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 xN
in funzione di xN-1
Siccome le rette su cui si dispongono i numeri sono molto dense,
per vedere l'effetto bisogna ingrandire molto una regione dello spazio che
contiene i punti, ma allora solo una piccola frazione di questi punti cadra`
nella regione selezionata ed il tempo di esecuzione sara` quindi molto grande,
perche' bisogna generare molti piu` 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.
Potrete verificare che il programma compilato gira molto piu` velocemente della macro. Questa e' una caratteristica generale da tenere a mente: i linguaggi interpretati, come le macro di ROOT che vengono tradotte in istruzioni per il calcolatore da un interprete nel momento stesso in cui devono venire eseguite, sono spesso molto piu` flessibili e permissivi dei linguaggi compilati, ma molto meno efficienti nell'esecuzione delle operazioni.
Si pone comunque il problema di visualizzare i grafici prodotti, dato che il programma fatto partireda linea di comando non ha la possibilita' di farlo. Ci sono due modi per affrontare il problema:
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;
E' 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. Tali tipi di oggetti sono accessibili solo all'interno di ROOT (vedremo in seguito come utilizzarli all'interno dei nostri programmi).
Oltre alla tabella riassuntiva, potete trovare:
Una macro puo' anche contenere la definizione di una nuova funzione, come la macro medialetture.C che definisce una funzione TH1F medialetture(char* nomefile,int N) che legge dei numeri da un file e ne fa la media di N alla volta e li mette in un istogramma.
In tal caso non ha senso
eseguire
direttamente la macro (come facciamo a passare i parametri giusti alla
funzione?). Si puo' pero' caricarla in memoria con il comando:
.L nomefile
A questo punto la nuova funzione e' nota a ROOT e puo' essere
utilizzata
da linea di comando o da altre macro. Ad esempio al prompt di ROOT si
puo` dare le seguente lista di comandi:
.L medialetture.c
TRandom gen(0);
ofstream g("miofile.dat");
for (int i=0; i<10000; i++) g << gen.Rndm() << endl;
g.close();
TH1F isto=medialetture("miofile.dat",1);
isto.Draw();
Se ci si accorge che la funzione ha degli errori e si vuole modificare la macro, si puo':
N.B.: poiche' ROOT esegue automaticamente l'inclusione di molti degli include file di sistema, spesso non c'e' bisogno delle dichiarazioni di #include nel file contentente la macro.
Il teorema del limite centrale, detto anche "legge dei grandi numeri"
e' la base di tutto il trattamento statistico degli errori
e grosso modo puo' essere formulato:
Siano x1, x2, x3, ...xN numeri casuali tutti estratti da una stessa
distribuzione con valor medio Xm e varianza σ2.
Allora, per N tendente all'infinito
la quantita' (x1+x2+...+xN)/N (la media
degli N numeri) tendera' ad avere una distribuzione
gaussiana
con valor medio Xm e varianza σ2/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.
Verificare empiricamente il teorema del limite centrale:
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 e' verificata per questa distribuzione?
Una forma di integrazione possibile e' 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.
Come esercizio, provare a calcolare l'integrale tra 0 e 1 della funzione
sin2(1/x) estraendo punti a caso nel quadrato
(0,1)x(0,1). Farsi stampare 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:
Proviamo ad applicare lo stesso metodo al calcolo di volumi di oggetti bidimensionali. Provare a scrivere una macro di ROOT che calcoli il volume di una sfera n-dimensionale di raggio 1: 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 piu` complicato, la seguente macro calcola il
potenziale gravitazionale della terra, assumendo che sia una sfera di
densita` uniforme, effettuando con il metodo Monte Carlo l'integrale:
double potenziale(double x, double y, double z) {
int Npunti=1000000;
TRandom generatore(0);
double G = 6.6742E-11 ; // costante di gravitazione universale
double massa= 5.9723E+24 ; // massa della terra
double raggio=6.5781E+06 ; // raggio della terra
double potenziale=0;
double volume=0;
for ( int i=0; i<Npunti; i++ ) {
double x1 = generatore.Uniform(-raggio,raggio);
double y1 = generatore.Uniform(-raggio,raggio);
double z1 = generatore.Uniform(-raggio,raggio);
if ( (x1*x1+y1*y1+z1*z1) < (raggio*raggio) ) {
volume+=1.;
double r=sqrt( pow(x-x1,2)+pow(y-y1,2)+pow(z-z1,2) );
potenziale-= G*massa/r;
}
}
return potenziale/volume;
}
Costruire un grafico del potenziale gravitazionale dal centro della sfera,
fino ad una distanza pari a 3 volte il raggio della sfera:
Si noti che ora la normalizzazione deve essere fatta in modo diverso per tenere conto del numero di punti effettivamente utilizzati.
Provare ad adattare la funzione di cui sopra al caso in cui la massa fosse in realta` concentrata in una crosta sferica di raggio 1000 km.
In questo caso e' forse piu` utile passare direttamente in coordinate polari, generando numeri uniformemente negli intervalli:
Alla fine di tutto, compilare la relazione