Introduzione al metodo Monte Carlo

Generalità e scopi

Metodo Monte Carlo: è una tecnica impiegata per calcolare probabilità e quantità correlate impiegando sequenze di numeri numeri casualicasuali. Non ha senso parlare di singoli numeri casuali, ma di sequenze

Viene utilizzato per: simulazione di fenomeni naturali simulazione di apparati sperimentali calcolo di integrali

Generazione di sequenze di numeri casuali

Procedura generale:

• Si estrae una sequenza di m numeri casuali, aventi distribuzione uniforme nell’intervallo (0,1).

• Si usa la sequenza per produrne una seconda distribuita secondo la funzione a cui si è interessati

• La sequenza costituisce il pool di dati simulati

Un esempio: la stima di π

Estrarre una coppia di numeri (a,b) uniformemente distribuiti in (0,1)

Si calcola c = √a2+b2

Al crescere del numero di eventi simulati, il rapporto

N(eventi con c<1) / N(eventi totali)

tende a

Area settore circolare/Area quadrato = π /4

Quindi π può essere il limite di

4 * N(eventi con c<1)/N(eventi totali)

Un esempio: la stima di π

Come generare sequenze di numeri casuali?

In una sequenza di numeri casuali, oltre all’uniformità della distribuzione nell’intervallo [0,1], si deve avere assenza di correlazione tra elementi successivi nella sequenza.

Esempi di generatori di numeri casuali: - Lancio di una monetina 32 volte per avere una parola di 32 bit. - Utilizzo di processi intrinsecamente casuali (arrivo di raggi cosmici, white noise di un fototubo, decadimento radioattivo,…) per associare ad essi dei numeri in un opportuno intervallo

Esistevano nel passato tabelle con sequenze di milioni di numeri casuali, a disposizione degli utenti

Questi metodi non sono ormai più usati, a favore di sequenze di numeri pseudo-casuali generati da calcolatori

Sequenze di numeri pseudo-casuali

Generatori di numeri pseudocasuali

La sequenza è riproducibile, perché l’algoritmo di generazione è deterministico

Caratteristiche a cui un generatore deve sottostare: indipendenza statistica

Periodo di ripetizione lungo ... la sequenza “sembra” casuale

Un esempio di generatori pseudo-casuali

Middle square (Von Neumann, 1946)

cifre m a numero ; di centrali cifre m 02

1 III nn

10m 0594909291333177923857721566492

InIn+1

Un esempio di generatori pseudo-casuali

Vengono generati numeri interi a m cifre. In generale gli algoritmi producono numeri interi tra 0 e MAX e restituiscono un valore in virgola mobile definito da rn=In/(MAX)

Il metodo di Von Neumann presenta dei problemi: se m è piccolo, la sequenza è breve

Usando 38 bit, si possono avere sequenze fino a 7.5105

numeri

656100008100

168100004100

44100002100

372100006100

2

2

2

2

Si riottiene 6100 !

Esistono tanti altri esempi di algoritmi per la generazione di sequenze pseudo-casuali

Distribuzione uniforme in un intervallo (a,b)

In molti problemi fisici è necessario utilizzare una sequenza di valori r’ casualmente distribuiti in un intervallo (a,b) in modo uniforme,partendo da sequenze di numeri r in (0,1).Trasformazione lineare:

r’ = a + r * (b-a)0 1

a b

r

r’

Esempi:

- Simulare l’angolo di emissione di una particella emessa isotropicamente in un piano, tra 0 e 360º

- Simulare la coordinata x di un punto sorgente in un collimatore rettangolare (avente l’altra dimensione trascurabile) θ

Esempio: valutare la distanza relativa tra 2 rivelatori

Caso 2D: rivelatori assunti come rettangoli

(dimensioni a x b)

Simulare una posizione (x1,y1) all’interno del primo rettangolo

Simulare una posizione (x2,y2) all’interno del secondo rettangolo

Valutare la distanza relativa tra i 2 punti

Ripetere per N eventi ed estrarre la distribuzione delle distanze relative, con il valore medio e RMS

Estensione al caso 3D: Assumere i rivelatori cilindrici, estrarre 2 punti di coordinate (x1,y1,z1) e (x2,y2,z2) all’interno dei 2 cilindri e valutare la distribuzione delle distanze relative

Distribuzione uniforme entro un cerchio

R

Generare dei punti uniformemente distribuiti entro un cerchio di raggio R:

-Estrarre x tra –R e +R

-Estrarre y tra –R e +R

-Accettare la coppia (x,y) se √x2+y2 < R

altrimenti estrarre una nuova coppia (x,y)

x

y

Esempio:

Simulare il punto di emissione di una particella dall’area di una sorgente radioattiva di raggio R

Distribuzioni non uniformi

In generale non sempre è sufficiente utilizzare sequenze di numeri casuali distribuiti uniformemente

In molti problemi è necessario disporre di numeri casuali estratti con densità di probabilità diverse, quali la distribuzione normale, oppure esponenziale, o poissoniana, etc..

Esempi:- Simulare l’energia depositata da una particella in

un rivelatore, con una distribuzione gaussiana intorno ad un valore medio e con una data sigma

- Simulare la direzione di arrivo di un muone cosmico su un rivelatore, con distribuzione cos2θ

Tre possibilità:a) Distribuzione uniforme con eventi pesatib) Metodo dell’inversionec) Metodo del rigetto

Distribuzione uniforme con eventi pesati

Metodo: Usare una distribuzione uniforme della variabile in (a,b), ma assegnando un peso relativo f(x) ad ogni evento.

Vantaggi: Molto semplice da usare

Svantaggi:

Produrre eventi che non sono tutti equivalenti

Non esattamente chiaro qual è il n. totale di eventi

Nelle zone in cui f(x) ha un valore basso, questi eventi di fatto non contribuiscono, anche se sono generati egualmente (perdita di CPU time)

Metodo dell’inversione

Problema: Estrarre una sequenza di numeri casuali distribuiti secondo la funzione f(x)

Data la funzione di probabilità f(x), la funzione integrale (cumulativa) F(t) è data da:

F(t) = ∫f(t)dt

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 0 e 1) e calcolando x da:

x = F-1 (u)

Vantaggi: molto efficiente, nessun numero viene sprecato

Svantaggi: La funzione inversa può essere difficile o impossibile da valutare

Metodo dell’inversione: esempio 1

f(x) = √x

F(x) = x2

Estrarre un numero 0< u <1

Calcolare x = u2

f(x)=½x per 0<x4

Si ricava x come

x=4u2

Metodo dell’inversione: esempio 2

f(x) = exp(-x/x0)

F(x) = -(1/x0 )ln x

Estrarre un numero 0< u <1

Calcolare x = -x0 ln u

Metodo dell’inversione: esempio 3 (Breit-Wigner o Lorentziana)

f(x) = (Γ/2)/ [( x-m0)2 + (Γ/2)2 ]

Descrive la forma di una risonanza, un livello energetico,…

F(u) = m0 + (Γ/2) tan u

dove u è un numero uniformemente distribuito in –π/2 < u < + π/2

Metodo del rigetto

Obiettivo: Estrarre una sequenza di numeri casuali secondo la distribuzione f(x), nell’intervallo (x1, x2)

Metodo:

a) Generare x uniforme in (x1,x2)

b) Generare y uniforme in (0,fmax)

c) Valutare f(x)

d) Confrontare y con f(x)

e) Se y < f(x) accettare x

f) Se y > f(x) ripetere da a) in poi

x1 x2

fmax

Accettato

Rigettato

Metodo del rigetto: pregi e difetti

Pregi:

- Metodo molto semplice

- Applicabile comunque sia complicata la funzione f(x)

Difetti:

-Poco efficiente dal punto di vista della potenza di calcolo (rigetto di molti numeri)

- Poco adatto specialmente quando la funzione f(x) ha dei picchi pronunciati o delle singolarità

Metodo del rigetto: miglioramenti possibili

Si può migliorare l’efficienza del metodo effettuando il campionamento non in un “rettangolo” ma in una regione definita da una funzione g(x) maggiorante di f(x).

Se si è in grado di generare numeri casuali distribuiti secondo g(x) per esempio con il metodo dell’inversione, la procedura dell’inversione diventa:

1. Si estrae x secondo g(x)2. Si estrae u con

distribuzione uniforme tra 0 e g(x)

3. Se u<f(x) si accetta x

Un esercizio: simulare il fenomeno del decadimento radioattivo

Un esercizio: simulare il fenomeno del decadimento radioattivo

Risultati/1

Risultato per: N0=100, =110-2 s-1

t=1 sLe fluttuazioni statistiche sono importanti

Linea continua: curva Linea continua: curva esponenziale attesa. esponenziale attesa. Istogramma: risultato Istogramma: risultato della simulazionedella simulazione

Risultati/2: Confronto tra 2 diverse sequenze di numeri casuali

Due diverse simulazioni, con differenti sequenze di numeri casuali, danno risultati leggermente diversi. Sono visibili le fluttuazioni statistiche.

Risultati/3: Riduzione delle fluttuazioni statistiche con il n. di eventi

Risultato per: N0=5000, =310-2 s-1

t=1 s

Risultati/4: Riduzione delle fluttuazioni statistiche con il n. di eventi

Risultato per: N0=50000, =310-2 s-1

t=1 s

Esempio: Come generare una direzione casuale nello spazio

x

z

y

θ

φ

- Estrarre l’angolo azimutale φ in (0, 2π)

- Estrarre l’angolo polare θ con distribuzione sinθ in (-π, +π)

-Angolo solido dΩ = dcosθ dφ = sinθ dθ dφ

Isotropia in angolo solido, non in angolo polare

Esempio: Emissione di particelle da una sorgente a forma di disco, con emissione isotropa nello spazio

Esempio: Il moto molecolare in 2 dimensioni

D = diametro molecolare (Angstrom)

P = Pressione (atmosfere)

T = Temperatura assoluta (K)

Libero cammino medio di una molecola in un gas:

P( l) = exp (-l/ L)

Θ uniforme in (0,2π)

Uso del teorema del limite centrale per una distribuzione gaussiana

Per generare una sequenza di numeri con distribuzione gaussiana, oltre che i metodi precedenti, si può sfruttare il teorema del limite centrale:

Per N molto grande, la somma di N numeri casuali xi è una variabile g distribuita secondo una distribuzione normale:

g = Σ (xi - ½ )

Assumendo N=12 (già sufficiente) la distribuzione gaussiana ha media 0 e deviazione standard 1.

Applicazioni ulteriori: calcoli di accettanza geometrica

Sorgente Rivelatore

Per una sorgente puntiforme e un rivelatore con superficie S a forma di calotta sferica di raggio R, posto a distanza d, l’angolo solido è

dΩ = S /d2

dove S = π R2

La probabilità che una particella vada dalla sorgente al rivelatore (accettanza geometrica), nel caso di emissione isotropa, è

εg = dΩ/4π = S/4πd2 = R2/4d2Correzione per l’accettanza:

Se misuriamo N conteggi nel rivelatore, dovuti a particelle emesse dalla sorgente, quelle effettivamente emesse dalla sorgente saranno N/ εg

Applicazioni ulteriori: calcoli di accettanza geometrica

Sorgente Rivelatore

Per una sorgente estesa e un rivelatore di forma qualunque, la probabilità che una particella vada dalla sorgente al rivelatore (accettanza geometrica), può essere valutata con tecniche Monte Carlo

Correzione per l’accettanza:

Estrarre un punto entro le dimensioni della sorgente

Estrarre una direzione di volo della particella

Verificare se essa intercetta il rivelatore

Ripetere N volte e valutare il rapporto εg = Nacc/Ngen

Utilizzo del metodo Monte Carlo per il calcolo di integrali multidimensionali

Un metodo numerico semplice per valutare l’ integrale definito della funzione f(x) in (a,b) è calcolare:

(b-a)/n Σ f(xi)

dove le xi sono equidistanziate in (a,b).

Il corrispondente metodo Monte Carlo consiste nel valutare la funzione f(x) in corrispondenza ad una sequenza di numeri casuali xi uniformemente distribuiti in (a,b).

La precisione del metodo Monte Carlo è in questo caso proporzionale a 1/√n, mentre per molti altri metodi numerici è dell’ordine di 1/n2, cioè molto migliore.

Perché allora il metodo Monte Carlo?

Soprattutto negli integrali multidimensionali (dimensione d)

In questo caso, la precisione di altri metodi numerici è n(-2/d) , mentre per il metodo Monte Carlo è sempre 1/√n. Quindi da d=4 in poi conviene.

Inoltre scegliere i valori da usare con altri metodi numerici può essere molto complicato, mentre è semplice con il metodo Monte Carlo

Un’applicazione avanzata di tecniche Monte Carlo: simulazioni GEANT

Nella progettazione di un apparato di rivelazione, dalla fisica nucleare alla fisica medica, alle applicazioni nello spazio, occorre simulare l’interazione di particelle con i materiali.

Un codice di simulazione generale: il software GEANT, esistente in versione GEANT3 (derivante dal Fortran) e in versione GEANT4 (C++)

Sviluppato al CERN nell’ambito della fisica particellare

Liberamente utilizzabile, ma molto complesso!

Un’applicazione avanzata di tecniche Monte Carlo: simulazioni GEANT

Valutare l’influenza che un blocco di Pb ha sulla risposta di uno scintillatore usato per la spettrometria gamma

Scintillatore

Pb

Sorgente gamma

Un’applicazione avanzata di tecniche Monte Carlo: simulazioni GEANT

Valutare l’effetto schermante dell’edificio sui muoni cosmici

TettoSolai in cemento

Pareti

Un’applicazione avanzata di tecniche Monte Carlo: simulazioni GEANT

Valutare la dose assorbita da un tessuto del corpo umano a seguito di radiazioni indotte da un fascio di particelle a scopo terapeutico

Fascio di particelle ionizzanti

Un’applicazione avanzata di tecniche Monte Carlo: simulazioni GEANT

Valutare la dose assorbita da astronauti nello spazio in viaggi interplanetari

Radiazioni cosmiche

Struttura di GEANT

Definire la geometria e i materiali del sistema da simulare, mediante delle forme predefinite

Struttura di GEANT

Definire la geometria e i materiali del sistema da simulare, mediante delle forme predefinite

Struttura di GEANT

Posizionare e orientare nello spazio i volumi che costituiscono il proprio sistema, all’interno di un volume “madre”

Anche decine di migliaia di singoli volumi… come per i rivelatori di LHC

Struttura di GEANT

Scegliere la/le particelle da generare (“vertice”)

Ad esempio una singola particella con una certa distribuzione in energia e angolo, generata in una data regione dello spazio

(una sorgente radioattiva)

Oppure, tutte le particelle generate da un modello teorico che descrive una collisione nucleare

Struttura di GEANT

Valutare tutte le possibili interazioni della particella primaria e di quelle secondarie prodotte nei rivelatori, in base alla conoscenza dei processi fisici incorporate in GEANT

Struttura di GEANT

Estrarre le informazioni prodotte nei rivelatori (energia depositata, tempo di arrivo delle particelle, posizione,…) e costruirne opportuni istogrammi

Un esempio di analisi su dati simulati

Diverse particelle (p, π, K in una certa proporzione, con distribuzione in energia di tipo esponenziale)

Rivelatore costituito da 3 rivelatori sottili e da un rivelatore spesso in silicio

Output: Un file di Root con 50000 eventi simulati che contengono:

N.evento, Tipo di particella, Impulso, E1, E2, E3, Etot (θ)

Esercitazione: Costruire gli spettri dell’energia depositata in ciascuno dei rivelatori

Un esempio di analisi su dati simulati

Esercitazioni:

1 Costruire gli spettri dell’energia depositata in ciascuno dei rivelatori

2 Costruire gli spettri selezionando le particelle incidenti

3 Costruire dei plot bidimensionali Ej-Etot per j=1, 2, 3 (matrici DeltaE-E)

4 Studiare l’effetto dell’impulso delle particelle incidenti

5 Studiare l’effetto dell’inclinazione della traccia