Introduzione al metodo Monte Carlo
Generalità e scopi
Metodo Monte Carlo: è una tecnica impiegata per calcolare probabilità e
quantità correlate impiegando sequenze di numeri casuali.
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 pseudocasuali 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)


I n1  m cifre centrali di I n2 ;
5772156649 2
In
I 0 numero a m cifre
 33317792380594909291
In+1
m  10
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
6100 2  37210000
2100 2  4410000
4100  16810000
2
Si riottiene 6100 !
8100 2  65610000
Usando 38 bit, si possono avere sequenze fino a 7.5105 numeri
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).
r
Trasformazione lineare:
r’ = a + r * (b-a)
0
1
r’
Esempi:
a
b
- 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
Generare dei punti uniformemente distribuiti
entro un cerchio di raggio R:
y
-Estrarre x tra –R e +R
-Estrarre y tra –R e +R
-Accettare la coppia (x,y) se √x2+y2 < R
x
altrimenti estrarre una nuova coppia (x,y)
R
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 pesati
b) Metodo dell’inversione
c) 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
Si ricava x come
x=4u2
f(x)=½x per
0<x4
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)
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
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,
=110-2 s-1
t=1 s
Le fluttuazioni statistiche sono importanti
Linea continua: curva
esponenziale attesa.
Istogramma: risultato
della 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,
=310-2 s-1
t=1 s
Risultati/4: Riduzione delle fluttuazioni statistiche con il n. di eventi
Risultato per:
N0=50000,
=310-2 s-1
t=1 s
Esempio: Come generare una direzione casuale nello spazio
- Estrarre l’angolo azimutale φ
in (0, 2π)
z
θ
y
x
φ
Isotropia in
angolo solido, non
in angolo polare
- Estrarre l’angolo polare θ con
distribuzione sinθ in (-π, +π)
-Angolo solido dΩ = dcosθ dφ =
sinθ dθ dφ
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
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
Sorgente
dove S = π R2
Rivelatore
Correzione per l’accettanza:
La probabilità che una particella
vada dalla sorgente al
rivelatore (accettanza
geometrica), nel caso di
emissione isotropa, è
Se misuriamo N conteggi nel rivelatore,
dovuti a particelle emesse dalla sorgente,
quelle effettivamente emesse dalla
sorgente saranno N/ εg
εg = dΩ/4π = S/4πd2 = R2/4d2
Applicazioni ulteriori: calcoli di accettanza geometrica
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
Sorgente
Rivelatore
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
Sorgente
gamma
Pb
Un’applicazione avanzata di tecniche Monte Carlo: simulazioni GEANT
Valutare l’effetto schermante dell’edificio sui muoni
cosmici
Solai in
cemento
Tetto
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