Simulazione Montecarlo di decadimento nucleare.
vorrei un consiglio a grandi linee su un programma che devo scrivere.
abbiamo un nucleo che emette dei fotoni a diverse energie, ognuno con una distribuzione angolare $f_E (theta)$ e un rate d'emissione costante nel tempo (R= n decadimenti/secondo). da questo poi devo ricavare delle cose che non ci importano al momento.
come organizzereste un programma di simulazione Montecarlo per questa situazione? (facciamo finta che sia in C++)
io pensavo di creare una classe evento, che per ogni evento mi dice l'angolo di emissione e il tempo. quando creo un evento, l'angolo viene attribuito con un generatore di numeri casuali adattato alla funzione $f_E (theta)$ e allo stesso modo gli viene dato casualmente il tempo di emissione.
qui c'è il problema, come faccio a creare n eventi al secondo partendo un generatore di numeri casuali uniforme?
forse è un dubbio scemo e c'è una soluzione banalissima, ma proprio non mi viene in mente nulla se non una cosa un po' macchinosa, che comporta creare molti più eventi e poi filtraril in seguito.
grazie dei consigli
abbiamo un nucleo che emette dei fotoni a diverse energie, ognuno con una distribuzione angolare $f_E (theta)$ e un rate d'emissione costante nel tempo (R= n decadimenti/secondo). da questo poi devo ricavare delle cose che non ci importano al momento.
come organizzereste un programma di simulazione Montecarlo per questa situazione? (facciamo finta che sia in C++)
io pensavo di creare una classe evento, che per ogni evento mi dice l'angolo di emissione e il tempo. quando creo un evento, l'angolo viene attribuito con un generatore di numeri casuali adattato alla funzione $f_E (theta)$ e allo stesso modo gli viene dato casualmente il tempo di emissione.
qui c'è il problema, come faccio a creare n eventi al secondo partendo un generatore di numeri casuali uniforme?
forse è un dubbio scemo e c'è una soluzione banalissima, ma proprio non mi viene in mente nulla se non una cosa un po' macchinosa, che comporta creare molti più eventi e poi filtraril in seguito.
grazie dei consigli
Risposte
Se il problema è ottenere una distribuzione particolare da una distribuzione uniforme dovresti utilizzare il metodo di Metropolis-Hastings. Sinceramente ora non ricordo molto bene i dettagli, ma serve proprio a questo scopo.
"Eredir":
Se il problema è ottenere una distribuzione particolare da una distribuzione uniforme dovresti utilizzare il metodo di Metropolis-Hastings. Sinceramente ora non ricordo molto bene i dettagli, ma serve proprio a questo scopo.
no no quello lo so fare (più che altro lo sa fare il nostro amico ROOT
)rispiego il problema, avevo divagato un po' e centrato poco la questione.
ho una classe evento(angolo,tempo) . quando creo un evento() vengono assegnati casualmente l'angolo (secondo la nostra distribuzione) e il tempo. qui è il problema: so che ho un rate R=n eventi/secondi. Voglio creare in successione gli eventi in un tempo di T secondi (saranno circa R*T, ma con una piccola distribuzione probabilistica). Come faccio operativamente a distribure gli eventi con il rate voluto? So che è un processo poissoniano, mi chiedo proprio come fare questa cosa operativamente. Forse mi perdo in un bicchier d'acqua perchè oggi sono un po' stanco, ma... boh!
ti dico come farei io, sperando vada bene per te.
Intanto il rate di decadimento è determinato da due grandezze: il numero di nuclei, la costante di decadimento $lambda$. Ti serve quest'ultima, che è data (sono sicuro che lo sai) da $lambda=-((dN(t))/(dt))/(N(t))$ in cui $N(t)$ è il numero di nuclei radioattivi non decaduti all'istante $t$.
Ovviamente, avrai deciso, prima di cominciare, la discretizzazione sia spaziale sia temporale, individuando i passi di integrazione $d theta$ e $dt$.
Fatto ciò, tu hai che la probabilità che un nucleo decada nel tempo $dt$ è $dp=lambdadt$. Per cui, richiamato un numero random x, lo confronti con $lambdadt$, con la condizione che se $x<=lambdadt$ allora il decadimento è avvenuto, altrimenti no, e reiteri avanzando il tempo di $dt$. Se hai N nuclei, dovrai ripetere ciò per N volte, prima di avanzare il tempo di $dt$.
Se il decadimento è avvenuto, hai da individuare l'angolo, e lo farai in maniera non troppo diversa (con l'accortezza di esserti prima costruito l'integrale di $f_E( theta)$, la funzione di distribuzione da confrontare con un altro numero casuale $x$).
Intanto il rate di decadimento è determinato da due grandezze: il numero di nuclei, la costante di decadimento $lambda$. Ti serve quest'ultima, che è data (sono sicuro che lo sai) da $lambda=-((dN(t))/(dt))/(N(t))$ in cui $N(t)$ è il numero di nuclei radioattivi non decaduti all'istante $t$.
Ovviamente, avrai deciso, prima di cominciare, la discretizzazione sia spaziale sia temporale, individuando i passi di integrazione $d theta$ e $dt$.
Fatto ciò, tu hai che la probabilità che un nucleo decada nel tempo $dt$ è $dp=lambdadt$. Per cui, richiamato un numero random x, lo confronti con $lambdadt$, con la condizione che se $x<=lambdadt$ allora il decadimento è avvenuto, altrimenti no, e reiteri avanzando il tempo di $dt$. Se hai N nuclei, dovrai ripetere ciò per N volte, prima di avanzare il tempo di $dt$.
Se il decadimento è avvenuto, hai da individuare l'angolo, e lo farai in maniera non troppo diversa (con l'accortezza di esserti prima costruito l'integrale di $f_E( theta)$, la funzione di distribuzione da confrontare con un altro numero casuale $x$).
"kinder":
Fatto ciò, tu hai che la probabilità che un nucleo decada nel tempo $dt$ è $dp=lambdadt$. Per cui, richiamato un numero random x, lo confronti con $lambdadt$, con la condizione che se $x<=lambdadt$ allora il decadimento è avvenuto, altrimenti no, e reiteri avanzando il tempo di $dt$. Se hai N nuclei, dovrai ripetere ciò per N volte, prima di avanzare il tempo di $dt$.
bravo, ottima idea.
mi ero instabilmente ostinato all'idea di generare gli eventi con un costruttore "puro", del tipo evento(), invece il tempo può essere assegnato alla classe esternamente da una serie di cicli come dici tu.
la discretizzazione spaziale non mi è necessaria, quando creo l'evento la classe gli assegna l'angolo in base alla distribuzione data, però per il tempo la tua procedura è sicuramente migliore. poi dopo devo vedere quando gli eventi colpiscono i rilevatori (che hanno una data distribuzione spaziale) in coincidenza con la loro probabilità di interazione eccetera eccetera, speriamo mi funzioni tutto
per adesso grazie mille kinder!
il metodo Montecarlo è concettualmente semplice, e scarica l'onere del calcolo sul computer; sono sicuro che ti funzionerà. In questi casi il nemico maggiore è il baco software.
Visto che devi simulare i conteggi che ottieni con un rivelatore, la rogna maggiore sarà la simulazione della sua geometria, molto più della sua posizione spaziale, a meno che non fai opportune ipotesi semplificative.
Visto che devi simulare i conteggi che ottieni con un rivelatore, la rogna maggiore sarà la simulazione della sua geometria, molto più della sua posizione spaziale, a meno che non fai opportune ipotesi semplificative.
Premetto che io di decadimento radioattivo non conosco quasi nulla, comunque ho dovuto risolvere un problema molto simile per il campionamento delle collisioni in un DSMC.
Sinceramente io il campionamento lo farei in un modo un po' diverso. Più che altro perché, magari mi sbaglio, mi sembra che così non si riesca a rendere la frequenza indipendente dal numero di atomi: cosa abbastanza importante a meno che non simuli tutti gli atomi, ovvero il tasso di emissione dipende dal numero di atomi che vengono usati nella simulazione, cosa che non dovrebbe essere visto che in un Monte Carlo il numero di "atomi" o "molecole" non ha alcun significato fisico, ma ha soltanto un ruolo statistico.
Io all'epoca avevo risolto (nel senso che avevo trovato che era il metodo standard nel Montecarlo e ho copiato
) così: fissato Dt, il passo temporale, ad ogni step creo una nuova variabile Tc che mi tenga conto del momento in cui avviene la collisione o, in questo caso, l'emissione radioattiva, quindi campiono i tempi da una variabile di Poisson (nel senso di processo di Poisson), avendo fissato preventivamente il numero di collisioni (la frequenza delle emissioni). Il codice (FORTRAN77) è questo:
dove tau è la funzione (scusate il pessimo Inglese nei commenti, ma tanto l'ho usato/scritto solo io questo codice):
in oltre facevo un controllo su ncoll verificando a posteriori che non ci fosse troppo discostamento dalla frequenza voluta (può accadere per Dt troppo piccoli). La chiave del metodo era proprio il garantire che le proprietà statistiche fossero perfettamente verificate, altrimenti si rischia che il tutto non converga.
Comunque, ripeto, magari tutto questo non centra nulla con il vostro caso, che non conosco, anche se mi sembrava simile...
Sinceramente io il campionamento lo farei in un modo un po' diverso. Più che altro perché, magari mi sbaglio, mi sembra che così non si riesca a rendere la frequenza indipendente dal numero di atomi: cosa abbastanza importante a meno che non simuli tutti gli atomi, ovvero il tasso di emissione dipende dal numero di atomi che vengono usati nella simulazione, cosa che non dovrebbe essere visto che in un Monte Carlo il numero di "atomi" o "molecole" non ha alcun significato fisico, ma ha soltanto un ruolo statistico.
Io all'epoca avevo risolto (nel senso che avevo trovato che era il metodo standard nel Montecarlo e ho copiato
) così: fissato Dt, il passo temporale, ad ogni step creo una nuova variabile Tc che mi tenga conto del momento in cui avviene la collisione o, in questo caso, l'emissione radioattiva, quindi campiono i tempi da una variabile di Poisson (nel senso di processo di Poisson), avendo fissato preventivamente il numero di collisioni (la frequenza delle emissioni). Il codice (FORTRAN77) è questo:
Tc=0.d0
15 continue
Tc=Tc+tau(avncoll2)
ncoll=ncoll+1
if(Tc.le.Dt) goto 15
dove tau è la funzione (scusate il pessimo Inglese nei commenti, ma tanto l'ho usato/scritto solo io questo codice):
c ====================================================================
double precision function tau(nu)
c ====================================================================
c This functions generates random numbers from an exponential
c distribution with paramether "nu"
c
include "param.f"
c
double precision nu,tmp
c
10 continue
tmp=1.d0-rand(idum)
if(tmp.le.1.d-16) goto 10
tau=-dlog(tmp)/nu
c
return
end
in oltre facevo un controllo su ncoll verificando a posteriori che non ci fosse troppo discostamento dalla frequenza voluta (può accadere per Dt troppo piccoli). La chiave del metodo era proprio il garantire che le proprietà statistiche fossero perfettamente verificate, altrimenti si rischia che il tutto non converga.
Comunque, ripeto, magari tutto questo non centra nulla con il vostro caso, che non conosco, anche se mi sembrava simile...
"kinder":
Visto che devi simulare i conteggi che ottieni con un rivelatore, la rogna maggiore sarà la simulazione della sua geometria, molto più della sua posizione spaziale, a meno che non fai opportune ipotesi semplificative.
ah si, quello è un altro casino. oggi siamo stati in 3 a sbatterci su problemi di geometria euclidea da vergognarsi
il problema finale è simulare le coincidenze tra due rilevatori (uno dei quali fa da gate energetico), non è banalissimo. a titolo di informazione al momento sto ottenendo ciò che m'aspettavo, ti ringrazio nuovamente per il consiglio che mi ha dato l'incipit. oggi mi han consigliato di abbandonare questo sistema di C++/ROOT e passare al GEANT4, che è un pacchetto specifico per l'interazione particelle-materia, lo conosci o hai mai provato? per david_e: due piccole note sul ragionamento iniziale. spesso sul campione c'è già segnata l'attività complessiva di tale campione, quindi quanti atomi ci siano poco importa. invece con una semplificazione molto grezza si può dire che osservare un decadimento su 100 atomi in un minuto ha la stessa probabilità di osservare il decadimento di un atomo in 100 minuti. (altrimenti non si potrebbe cercare il decadimento del protone che ha un tempo di dimezzamento molto più grande della vita dell'universo).
sul programma invece ammetto che fatico molto a leggere il FORTRAN, cosa significano i parametri dentro gli if?
ncoll sarebbe il numero di eventi (va poi decrementato di 1 a fine ciclo) che poi verranno generati.
Tc e' l'istante in cui avviene il singolo evento (si esclude la possibilita' di eventi contemporanei).
Come ho detto prima, magari mi sbaglio: non conosco bene il problema e non sono certo un esperto di Montecarlo. Comunque penso che agli atomi della simulazione Montecarlo vada assegnata una probabilita di decadimento che non sia quella reale del singolo atomo, ma che deve essere riscalata in modo che la frequenza sia la stessa del limite per $N \to \infty$, quella poi che si osserva sperimentalmente se lavori nel limite del numero di atomi molto grande. Altrimenti rischi di simulare un campione con $N$ atomi, con $N$ il numero di particelle della simulazione, che ovviamente ha un comportamento diverso dal campione sperimentale che, penso, abbia un numero di atomi paragonabile col numero di Avogadro.
Se devi simulare un numero conosciuto, piccolo di atomi allora il discorso e' diverso...
Tc e' l'istante in cui avviene il singolo evento (si esclude la possibilita' di eventi contemporanei).
Come ho detto prima, magari mi sbaglio: non conosco bene il problema e non sono certo un esperto di Montecarlo. Comunque penso che agli atomi della simulazione Montecarlo vada assegnata una probabilita di decadimento che non sia quella reale del singolo atomo, ma che deve essere riscalata in modo che la frequenza sia la stessa del limite per $N \to \infty$, quella poi che si osserva sperimentalmente se lavori nel limite del numero di atomi molto grande. Altrimenti rischi di simulare un campione con $N$ atomi, con $N$ il numero di particelle della simulazione, che ovviamente ha un comportamento diverso dal campione sperimentale che, penso, abbia un numero di atomi paragonabile col numero di Avogadro.
Se devi simulare un numero conosciuto, piccolo di atomi allora il discorso e' diverso...
"wedge":
oggi mi han consigliato di abbandonare questo sistema di C++/ROOT e passare al GEANT4, che è un pacchetto specifico per l'interazione particelle-materia, lo conosci o hai mai provato?
Non lo conosco. Non so neppure se esisteva quando facevo il ricercatore (più di 12 anni fa). Comunque, in genere mi facevo le cose in casa, in FORTRAN, incluse le routine specificatamente matematiche..
Accedi a tutti gli appunti
Tutor AI: studia meglio e in meno tempo