Sommario si ha Q = m nc 2 m pc 2 m e c 2 (1.6) 1.3.2 Decadimento + Nel seguente processo avviene un...

Sommario

Il problema che si vuole analizzare in questo elaborato finale è relativo alla produ-zione di radiofarmaci. Queste sostanze sono fondamentali per la cura di determinatemalattie e rappresentano una delle più importanti conquiste della medicina nucleare.La loro produzione però non è semplice e comporta delle notevoli difficoltà, una dellequali verrà qui analizzata.Quando un campione viene creato infatti, all’interno dello stesso non è presente so-lo l’elemento desiderato per la produzione del farmaco, ma possono essere presentianche altre sostanze, create per diversi motivi, che hanno la possibilità di risultaresia innocue per l’organismo umano, sia dannose, se non addirittura velenose.

Si rende quindi necessario sviluppare un sistema di controllo su questi campionidi farmaci prodotti per determinarne la composizione specifica e le quantità relativedi ogni elemento al fine di indicarne l’eventuale tossicità.I radiofarmaci hanno come principale caratteristica quella di essere composti daelementi con nuclei radioattivi, in particolare i decadimenti interessati sono di tipoβ. Ogni specie di radionuclide ha uno spettro particolare di decadimento ricavatosperimentalmente e attraverso quest’ultimo è possibile identificarlo.Verrà qui presentata una simulazione via software della produzione di un campionedi radiofarmaco composto da due elementi, lo scopo sarà quello di risalire alle suecomponenti attraverso dei fit matematici dello spettro elaborato dal software.

Indice

1 Introduzione al decadimento β 11.1 Panoramica generale . . . . . . . . . . . . . . . . . . . . . . . . . . . 1

1.1.1 Tempi di vita del decadimento β . . . . . . . . . . . . . . . . 21.2 Spettro β e introduzione al neutrino . . . . . . . . . . . . . . . . . . 31.3 Tipi di decadimento β ed energia rilasciata . . . . . . . . . . . . . . 3

1.3.1 Decadimento β− . . . . . . . . . . . . . . . . . . . . . . . . . 41.3.2 Decadimento β+ . . . . . . . . . . . . . . . . . . . . . . . . . 41.3.3 Cattura elettronica (EC) . . . . . . . . . . . . . . . . . . . . 5

1.4 Spettro β statistico . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51.5 Massa del neutrino . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71.6 Distorsione coulombiana dello spettro . . . . . . . . . . . . . . . . . 81.7 Il rapporto di densità elettronica F (Z,W ) . . . . . . . . . . . . . . . 91.8 Test χ2 per l’attendibilità dei risultati ottenuti . . . . . . . . . . . . 11

1.8.1 Introduzione statistica . . . . . . . . . . . . . . . . . . . . . . 11

2 Distribuzione di probabilità 142.1 Distribuzione di probabilità uniforme . . . . . . . . . . . . . . . . . . 142.2 Distribuzione di probabilità non uniforme . . . . . . . . . . . . . . . 17

2.2.1 Mappatura a canali . . . . . . . . . . . . . . . . . . . . . . . 19

3 Simulazione dello spettro β− 213.1 Elementi radioattivi considerati: 9039Y51 e

9038Sr52 . . . . . . . . . . . . 21

3.2 Distribuzione di probabilità . . . . . . . . . . . . . . . . . . . . . . . 223.3 Campionamento dei valori generati . . . . . . . . . . . . . . . . . . . 253.4 Fit con dati sperimentali . . . . . . . . . . . . . . . . . . . . . . . . . 27

Capitolo 1

Introduzione al decadimento β

1.1 Panoramica generale

Come abbiamo già detto, il processo di decadimento radioattivo che coinvolge iradionuclidi in questione è di tipo β. Questo processo consiste nella trasformazionedi particolari nuclei in altri, con emissione di particelle β dotate di carica elettrica.Queste possono essere osservate attraverso l’applicazione di un campo magneticoche ne devia la traiettoria, come si può vedere in fig.1.1.

Figura 1.1: Curvatura di particelle e onde prodotte dalle emissioni radioattive.

Le particelle β hanno carica negativa e sono state identificate negli elettroni, leparticelle α (nuclei di elio) e i raggi γ (radiazione elettromagnetica) sono prodottidi altri tipi di decadimenti nucleari e hanno carica rispettivamente positiva e nulla.Gli elettroni in questione non provengono dagli orbitali elettronici dell’atomo, mavengono prodotti all’interno del nucleo nell’atto dell’emissione. Da analisi chimicherisulta che la carica totale del nucleo padre differisce di una unità rispetto a quelladel nucleo figlio e tale cambiamento di carica consiste nella trasformazione di unnucleone, in particolare del neutrone, in un protone e un elettrone che si allontana:n → p+e-.

Un procedimento simile accade anche ai protoni che si trasformano in un neutrone e

1

Capitolo 1. Introduzione al decadimento β

un antielettrone (positrone): p→ n+e+, anche se questo non avviene per un protonelibero.Ciò è dovuto al fatto che il protone ha meno massa del neutrone e quindi si haMpc

2 < Mnc2, perciò l’energia non è disponibile. Tuttavia se il protone che si deve

trasformare, e il neutrone che si deve creare, si trovano in uno stato legato in unastruttura atomica, allora l’energia potrebbe essere disponibile.

In generale il decadimento β consiste nella trasformazione nucleare seguente: Z → Z ± 1,con Z numero atomico, e N → N ∓ 1, con N numero di neutroni, cos̀ı che il numerodi massa atomica A = Z +N rimanga invariato.

1.1.1 Tempi di vita del decadimento β

L’osservazione del processo di decadimento β può essere messa in pratica attraversola selezione di isotopi per mezzo di procedure chimico-fisiche, oppure, come accadepiù spesso, attraverso l’osservazione di elettroni β emessi ad un tasso specifico perogni specie di un radioisotopo.Il tasso di emissione risulta essere proporzionale al numero N di nuclei padri presentiin ogni momento che ancora non sono decaduti secondo la relazione

−dNdt

= λN (1.1)

dove λ è la costante di decadimento caratteristica dell’isotopo, che è tale perchél’attività di decadimento è un esponenziale decrescente con il tempo

−dNdt

= λN0e−λt (1.2)

Tale risultato segue proprio dal fatto che λ = − dNNdt sia una costante.Il numero dei nuclei di un isotopo β-radioattivo sono ridotti di una frazione 1e inogni intervallo di tempo τ = 1λ , dove τ è il tempo di vita medio dell’isotopo.

Il tempo che occorre affinché metà dei radioisotopi del campione decada è dettotempo di dimezzamento

t 12

= τ ln2 = 0.693 · τ (1.3)

Il range di variazione per le emivite degli isotopi esistenti è molto vasto, ad esempioinfatti il 5B

127 , un isotopo del boro con 5 protoni e 7 neutroni, ha un’emivita di 0.027

s, mentre il 19K4021 ha un’emivita di 1.5 · 109anni.

2


1.2 Spettro β e introduzione al neutrino

Un nucleo β-radioattivo ha una massa iniziale definita Mi e il nucleo figlio prodottodal decadimento ha una massa Mf . Come si è detto detto il nucleo padre differiscedal nucleo figlio per la sostituzione di un neutrone con un protone n↔p.Il processo di decadimento emette quindi una quantità di energia pari aW0 = (Mi − Mf )c2, parte di questa energia viene utilizzata per creare la massam della particella emessa e l’altra parte viene trasformata in energia cinetica dellaparticella stessa. Tuttavia solo una parte infinitesima delle particelle emesse dal cam-pione segue questo schema di ripartizione dell’energia, cioè solo pochissime hannoun’energia cinetica T = W0 −mc2.

Figura 1.2: Spettro sperimentale delle energie.

Nella fig.1.2 viene mostrato l’andamento della frazione di decadimenti che pro-ducono gli elettroni con energia cinetica T, su una scala per unità di energia. Ilproblema è determinare cosa accade all’energia che rimane libera W0 −mc2 − T intutti i casi in cui l’elettrone ha a disposizione l’energia cinetica T < W0 −mc2.Per risolvere questa questione Pauli propose l’idea dell’esistenza di una particellafino ad allora sconosciuta. Questa doveva essere di carica neutra e priva di qualsiasiinterazione con la materia, altrimenti sarebbe stata scoperta con gli esperimenti fattiin precedenza. Fermi chiamò questa particella neutrino.

1.3 Tipi di decadimento β ed energia rilasciata

Il processo fino ad ora descritto non è unico ma ne esistono di tre generi, che sonoβ−, β+ e cattura elettronica EC.

3


1.3.1 Decadimento β−

In questo processo si ha la trasformazione del neutrone all’interno del nucleo in unprotone, un elettrone libero e un antineutrino ν̄e, ovvero l’antiparticella del neutrino.

n→ p+ e− + ν̄e (1.4)

Il decadimento può essere quindi schematizzato come segue

AZXN →AZ+1 X

′N−1 + e

− + ν̄e

Q-valore

Per fare delle considerazioni energetiche sul decadimento si può indicare la massa diun atomo come

Matomo(X) · c2 = MN (X) · c2 + Z ·me · c2 −Z∑i=1

Bi

dove con Bi si è indicata l’energia di legame dei singoli elettroni con il nucleo e conZ il numero atomico, mentre MN (X) indica la massa del nucleo.La differenza di massa che si viene a creare da atomo padre a atomo figlio, moltipli-cata per c2, è detta Q-valore del decadimento e rappresenta l’energia rilasciata.Eseguendo tale differenza gli elementi uguali si elidono a vicenda e rimane

Q = {MN (Z,A)− [MN (Z + 1, A) +me +mν̄e ]} · c2 (1.5)

dove con MN (Z,A) e MN (Z + 1, A) sono state indicate rispettivamente la massadel nucleo padre e quella del nucleo figlio. La massa dell’antineutrino in questo casoviene trascurata con buon grado di approssimazione. La (1.5) può essere ridotta aisoli elementi in gioco, ovvero il neutrone, il protone e l’antineutrino visto che tuttoil resto rimane invariato.In questo primo caso di β− si ha

Qβ− = mnc2 −mpc2 −me−c2 (1.6)

1.3.2 Decadimento β+

Nel seguente processo avviene un altro tipo di trasformazione nucleare, ovvero unprotone si trasforma in un neutrone, un positrone (antiparticella dell’elettrone) e unneutrino.

p→ n+ e+ + νe (1.7)

Anche in questo caso si può schematizzare il processo nel seguente modo

AZXN →AZ−1 X

′N+1 + e

+ + νe

4


il numero dei protoni diminuisce di una unità, quello dei neutroni aumenta di unae il numero di massa A rimane invariato. Come detto precedentemente questodecadimento non può avvenire per un protone libero, ma si può verificare in unnucleo dove l’energia in eccesso per creare la differenza di massa tra neutrone eprotone e per la creazione del positrone può essere disponibile.In questo caso il Q− valore è

Qβ+ = mpc2 −mnc2 −me+c2 (1.8)

1.3.3 Cattura elettronica (EC)

Questo tipo di trasformazione avviene per i nuclei instabili ricchi di protoni. Infattiun elettrone che si trova negli orbitali più interni, come quello s, si lega a un protonedel nucleo per creare un neutrone legato e un neutrino. L’elettrone che si trasformalascia un orbitale libero che viene occupato da un elettrone che si trovava in unorbitale più esterno. Questo processo porta all’emissione di raggi X.

p+ e− → n+ νe (1.9)

Il numero di protoni diminuisce di una unità, mentre quello dei neutroni aumentadi una.

AZXN + e

− →AZ−1 X′N+1 + νe

Mentre per il Q-valore si ottiene

QEC = mpc2 +me−c

2 −mnc2 (1.10)

1.4 Spettro β statistico

Si consideri un nucleo che si divide in tre parti, con associate quantità di moto p, qe P che rispettano la relazione

P + p + q = 0 (1.11)

Inoltre la conservazione dell’energia impone che

Mic2 ≡Mfc2 + E +W +K (1.12)

Dove Mi e Mf sono rispettivamente la massa iniziale e quella finale. E rappresental’energia di rinculo del nucleo, che può essere trascurata come da evidenze speri-mentali, W e K sono le energie rispettivamente dell’elettrone e del neutrino checorrispondono alle quantità di moto p e q. Esse includono anche l’energia a riposodi queste particelle

W =√

(m2ec2 + p2)c2 (1.13)

K =√

(µc2 + q2)c2 (1.14)

5


dove µ è la massa del neutrino che consideriamo per completezza. Si assume chel’elettrone e il neutrino si dividano l’energia ’statisticamente’. Questo si traduce nelfatto che la probabilità di manifestazione di un sottogruppo di stati sia proporzionaleproprio al numero degli stati stessi.Il numero di stati per un elettrone liberato con quantità di moto p in un intervallodp = dpxdpydpz è proporzionale alla grandezza dell’intervallo dp e, analogamente,il numero di stati del neutrino in un intervallo dq è proporzionale a dq.Se la quantità dλλ indica la frazione dei decadimenti nei quali le quantità di motoconsiderate sono ricavate, allora

dλ

λ=

dpdq∫...∫

(dp)(dq)(1.15)

Al denominatore si ha l’intero spazio disponibile delle quantità di moto che rispet-tano la condizione energetica

K +W = W0 (1.16)

dato che il nucleo pesante assorbe assorbe la quantità di moto di cui dispongono leparticelle più leggere senza avere un guadagno di energia apprezzabile.Si può fare un cambio di coordinate passando a quelle sferiche

dq = q2dqdΩν

dp = p2dpdΩe

Ων , Ωe rappresentano gli angoli solidi verso i quali sono dirette le quantità di motonegli intervalli dp e dq. Si può sostituire qdq con KdK

c2.

La condizione sull’energia (1.16) impone che l’intervallo di integrazione dK dell’e-nergia del neutrino sia limitata a un elemento dK → 0. Si può quindi concludereche la (1.15) assume la forma

dλ

λ=

qKp2dp∫ p00 qKp

2dp(1.17)

Questa rappresenta la frazione dei decadimenti che producono elettroni con quan-tità di moto comprese nell’intervallo dp. Il limite superiore di integrazione p0 =

Qc

è la quantità di moto massima dell’elettrone ottenuta quando quella del neutrino siannulla q → 0.

La frazione dei decadimenti che produce elettroni con energia nell’intervallo dWsi ottiene sostituendo pdp = WdW

c2nella (1.17), da cui si ha che

dλ

dW=

λ

f0W (W0 −W )

√W 2 −m2c4

√(W0 −W )2 − µ2c4

1

(mc2)5(1.18)

dove f0(p0) =∫ p0

0

qKp2dp

(mc)5c. Il limite superiore di integrazione

p0 =1

c·√

(W0 − µc2)2 −m2c4

6


è la massima quantità di moto dell’elettrone ottenuta quando il neutrino ha unaquantità di moto che si annulla, ovvero K → µc2.

Lo spettro (1.18) si annulla quando l’elettrone non riceve energia cinetica, ovveroquando W = mec

2, ma anche quando non ne riceve il neutrino W = Q+mec2, con

Q = Te +K e Te = W −mec2 è l’energia cinetica dell’elettrone.

Lo spettro (1.18) è mostrato in fig.1.3

Figura 1.3: Spettro derivante dalla ripartizione statistica dell’energia tra elettronee neutrino.

1.5 Massa del neutrino

Il fatto che la maggior parte degli elettroni osservati abbia meno della metà dell’e-nergia di end-point, ovvero l’energia massima data da Emax = (W + K) − mec2,evidenzia il fatto che il neutrino debba avere una massa a riposo minore di quelladell’elettrone.Il più leggero dei tre frammenti del decadimento che partecipano all’equilibrio ca-suale delle quantità di moto deve avere in media più energia per avere una quantitàdi moto comparabile.In fig.1.4 i grafici di (1.18) per vari valori di massa 0 ≤ µ ≤ m e l’energia cineticaspecifica di end-point W0 − (m+ µ)c2 = 0.036mc2 è caratteristica del trizio H3.Un confronto dettagliato tra lo spettro sperimentale di H3, vicino all’end-point, condiversi valori attesi per la massa a riposo del neutrino, portarono Langer e Moffata concludere che µ < 0.0005m. Si può dedurre quindi che la massa a riposo delneutrino sia trascurabile e considerata nulla per i nostri scopi.

7


Figura 1.4: Paragone con varie masse associate al neutrino.

1.6 Distorsione coulombiana dello spettro

Osservando gli spettri sperimentali si può notare che essi non coincidono perfetta-mente con (1.18), ciò accade soprattutto con i nuclei pesanti. Quando una particellaatomica viene emessa nel processo di decadimento, un elettrone o un positrone, essasubisce l’effetto della carica positiva del nucleo che può rispettivamente rallentare oaccelerare la particella.In fig.1.4 sono riportati gli spettri di decadimento per una atomo di Cu64. Le lineetratteggiate rappresentano i dati sperimentali, mentre le linee continue quelli teori-ci. Da qui si può apprezzare l’effetto di attrazione (o repulsione) coulombiana sullaparticella emessa.

Figura 1.5: Confronto tra spettro teorico (tratteggiato) e dati statistici (continuo)

Si osserva anche che le particelle più veloci sono le meno influenzate dalla cari-ca elettrostatica del nucleo. Risulta possibile però incorporare questo effetto nella(1.18) prendendo in considerazione la densità spaziale |Ψc|2 dell’elettrone, dove Ψc

8


rappresenta la soluzione dell’equazione di Schrödinger. La densità spaziale descrivel’elettrone in uno stato di moto da quando compare in r = 0, cioè nel nucleo +Zefino a quando raggiunge la quantità di moto finale all’esterno del raggio di azionedella forza elettrostatica nucleare (r →∞).Questa densità può essere comparata con quella di una particella libera non soggettaa nessuna forza, |Ψ0|2, in uno stato in cui la quantità di moto non si modifica. Infig.1.6 viene mostrato l’andamento del rapporto tra queste due quantità.

Figura 1.6: Confronto tra particella libera e soggetta alla froza elettrostatica

Si nota che la densità nel nucleo varia con la carica nucleare presente. Il rapportoΨcΨ0

tende a 1 per la distanza che va ad infinito perché |Ψc| → |Ψ0|, dove l’intensitàdella forza elettrostatica tende a svanire. Le diverse funzioni d’onda degli elettronisono quindi eccitate proporzionalmente alla loro densità nel nucleo.

Il rapporto di densità elettronica è definito come F (Z,W ) =

∣∣∣∣Ψc(0)Ψ0∣∣∣∣2 e si può

assumere che dλ ∼ F (Z,W ) cos̀ı che la (1.18) può essere corretta nel seguente modo

dλ

dW=F (Z,W )

f(Z,Q)W (Q+mec

2 −W )2√W 2 −m2ec4(mec2)5

(1.19)

dove f(Z,Q) =∫ Q+mec2mec2

dWcpW (Q+mec2 −W )2F (Z,W ) · 1

(mec2)5

1.7 Il rapporto di densità elettronica F (Z,W )

Una stima di F (Z,W ) =

∣∣∣∣Ψc(0)Ψ0∣∣∣∣2 può essere ricavata dalla soluzione dell’equazione

di Schrödinger o Dirac (caso relativistico) per un elettrone immerso in un campoelettrostatico.

9


L’elettrone β emesso nel processo di decadimento solitamente ha una lunghezza d’on-da associata che è molto maggiore delle dimensioni nucleari, quindi si può trattarela carica del nucleo come puntiforme. Allo stesso tempo le lunghezze d’onda per lagran parte dello spettro sono minori delle dimensioni atomiche, cos̀ı che lo screeningcon elettroni atomici può essere trascurato.Si può considerare l’approssimazione non relativistica per atomi con carica non trop-po elevata e utilizzare la soluzione dell’equazione di Schrödinger per un elettroneimmerso in un campo elettrostatico generato da una carica puntiforme, ovvero Ze

2

r .

Segue che

F (Z,W ) =2πν

1− e−2πν(1.20)

dove ν = ±Ze2

~v per e∓, v è la velocità dell’elettrone (o positrone) e Z è il numero di

carica del nucleo figlio. Si nota che F = 1 per Z = 0 e per v → ∞, come dovrebbeessere.

10


1.8 Test χ2 per l’attendibilità dei risultati ottenuti

I radionuclidi che sono necessari alla produzione di radiofarmaci vengono prodottiattraverso il bombardamento di campioni di interesse tramite protoni o nuclei dideuterio 2H. Da questo procedimento si generano non soltanto gli elementi sperati,ma anche altre sostanze che possono essere dannose per l’organismo umano.Risulta di particolare interesse determinare le varie concentrazioni degli elementiprodotti perché alcuni di essi possono essere tossici anche in quantità minime, quin-di bisogna fornire i dati con una incertezza che deve essere inferiore ad una soglia disicurezza. Ad esempio, se in un campione viene indicato che la presenza di un certoelemento è pari a 0 mg con un’incertezza di 1 mg, ma è noto che quell’elementostesso risulta tossico per quantità anche leggermene inferiori a 1 mg, allora è chiaroche questo risultato non può essere utilizzato.

Un altro aspetto da considerare per i risultati che si otterranno nel prossimo ca-pitolo riguarda l’attendibilità delle percentuali di concentrazione degli elementi delcampione. Quando si fornisce un dato è necessario avere uno strumento che indichi,come in questo caso, la possibilità che la percentuale ricada con una probabilità piùo meno elevata nel range indicato dalle previsioni teoriche, per fornire un grado diattendibilità della misura. Viene cos̀ı in aiuto uno strumento statistico, il χ2, che haproprio lo scopo di verificare che i valori osservati si adattino alle previsioni teoriche.

1.8.1 Introduzione statistica

Si consideri un campione in cui vengono osservati una serie di eventi E1, E2, E3, . . . , ENle cui frequenze sperimentali di osservazione sono f1, f2, f3, . . . , fN .Si considerino poi le frequenze F1, F2, F3, . . . , FN , che sono le frequenze di manife-stazione degli stessi eventi che ci si aspetterebbe dalle leggi della statistica.

Tabella 1.1: Tabella riassuntiva.

Evento E1 E2 E3 ... EN

Frequenze osservate f1 f2 f3 ... fNFrequenze teoriche F1 F2 F3 ... FN

La variabile χ2 può essere definita come una somma sulle differenze al quadratotra le frequenze osservate e quelle teoriche per ogni evento Ei, pesate sulle frequenzeteoriche

χ2k =k∑i=1

(fi − Fi)2

Fi(1.21)

Nella fig.1.7 viene mostrata la distribuzione di Pearson della funzione χ2. Il p-valuerappresenta la probabilità di ottenere un valore maggiore o uguale a quello calcolato.

11


Figura 1.7: Distribuzione cumulativa della variabile χ2.

Dalla teoria dei cumulanti si può ottenere l’andamento della funzione densità diprobabilità fk(x) per la distribuzione χ

2 mostrata in fig.1.8 e si nota che tale distri-buzione tende ad una Gaussiana all’aumentare dei gradi di libertà k.

Figura 1.8: Densità di probabilità della distribuzione χ2.

Considerando un numero fissato di eventi N , la variabile (1.21) si distribuisce conk − 1 gradi di libertà, mentre se N rappresenta una variabile aleatoria (ad esempiodi Poisson) allora i gradi di libertà sono k. Una proprietà fondamentale della distri-

12


buzione χ2 è che la somma di due distribuzioni con gradi di libertà diversi è ugualea una distribuzione χ2 con grado di libertà dato dalla somme degli altri due

χ2n + χ2m = χ

2n+m (1.22)

Si può ora utilizzare questo strumento per determinare l’attendibilità delle frequenzeottenute di una serie di eventi imponendo un valore di tolleranza, partendo con ilfissare un errore α che può essere tollerato per la misura. Si confronta poi con latavola delle distribuzioni χ2 in fig.1.9. Solitamente viene considerato 0.05 come ilvalore critico che distingue risultati significativi da quelli non significativi.

Figura 1.9: Tabella della probabilità di andare oltre il valore critico.

A titolo di esempio, nel caso più semplice del lancio di una moneta si ha unsolo grado di libertà. Calcolando la variabile χ2 con i valori delle frequenze attesee osservate si può avere un’idea se la moneta sia truccata o meno fissando comevalore critico 0.05 e considerando la relativa colonna. Se χ2 < 3.84, allora si puòaccettare l’ipotesi nulla, ovvero che per i vari eventi non ci sia una differenza statisticasignificativa tra i valori osservati e quelli attesi, mentre se χ2 > 3.84 allora non sipuò accettare l’ipotesi nulla e quindi c’è una differenza significativa.

13

Capitolo 2

Distribuzione di probabilità

2.1 Distribuzione di probabilità uniforme

Si consideri una semplice distribuzione di probabilità uniforme di valori che vannoda a a b come mostrato in fig.2.1 dove nell’ascissa è indicato il valore della densitàdi probabilità, che è appunto una costante.

a b

1

b− a

x0

ρ(x)

Figura 2.1: Distribuzione uniforme

La relativa funzione di ripartizione F (x) è mostrata in fig.2.1. La funzione

14

Capitolo 2. Distribuzione di probabilità

x0

F (x)

a b

1

Figura 2.2: Distribuzione uniforme

densità di probabilità deve rispettare la relazione∫ baρ(x)dx = 1 (2.1)

ovvero deve essere normalizzata sull’intervallo in questione.Si consideri ora la variabile

x = f(y)

come funzione di un’altra variabile y che si può ricavare invertendo la precedenterelazione ottenendo

y = f−1(x)

La (2.1) diventa quindi ∫ f−1(b)f−1(a)

ρ(f(y))df

dydy = 1 (2.2)

Imponendo le seguenti relazioni sui nuovi estremi di integrazione

f−1(a) = 0 ⇒ f(0) = a

f−1(b) = 1 ⇒ f(1) = b

dalla (2.1) segue che

ρ(f(y))df(y)

dy= 1 → ρ(f(y))f ′(y) = 1 (2.3)

Una forma per la x = f(y) che rispetti le condizioni sopra indicate può essere

x ≡ f(y) = a+ (b− a)y (2.4)

15


Poiché è necessario ricavare un’espressione per la ρ(f(y)), seguendo la (2.3) si puòderivare la (2.4) ottenendo

df(y)

dy=d[a+ (b− a)y]

dy= (b− a) (2.5)

Invertendo la (2.3) si trova che

ρ(f)df(y)

dy= 1 (2.6a)

f ′(y) = b− a (2.6b)

ρ(f(y)) =1

b− a(2.6c)

La (2.6a) è un’equazione differenziale che può essere risolta per separazione dellevariabili

df

dy=

1

ρ(f)

ρ(f)df = dy

(2.7)

da cui si ricava che ∫ρ(f)df =

∫dy (2.8)

Indicando con P(x) una primitiva della funzione integranda a sinistra della (2.8)risulta

P(x) = y + c (2.9)

con c costante di integrazione. Imponendo le condizioni

P(a) = 0

P(b) = 1

si può dire che una funzione che rispetti queste relazioni sia

P(x) = x− ab− a

= y (2.10)

come detto precedentemente. Inserendo questo tipo di distribuzione in un pro-gramma software che genera un numero arbitrario di valori random (in questo caso10.000) compresi tra 0 e 2 si ottiene il grafico in fig.2.3. Si nota che la disposizionedei dati non è perfettamente omogenea perché comunque è presente una fluttuazionestatistica.

16


0 0.2 0.4 0.6 0.8

1 1.2 1.4 1.6 1.8

2

Figura 2.3: Grafico di valori random con distribuzione uniforme.

2.2 Distribuzione di probabilità non uniforme

Consideriamo ora una distribuzione di probabilità non uniforme del tipo in fig.2.4

x0

ρ(x)

2

1

Figura 2.4: Distribuzione non uniforme

che segue la relazione

ρ(x) = 1− x2

17


Una primitiva di questa funzione P̃(x) risulta essere

P̃(x) = x− x2

4= y (2.11)

Si pone per semplicità P̃(x) = P(x).

Risolvendo la (2.11) per x si ha

x2

4− x+ y = 0 (2.12a)

2± 2√

1− y (2.12b)

Delle due di (2.12b) si prende quella con il meno perché la P(x) deve rispettare leseguenti condizioni

P(a) = 0 → P−1(0) = aP(b) = 1 → P−1(1) = b

(2.13)

in questo caso si ha a = 0 e b = 2. Lanciando nuovamente il programma generatoredi numeri random e seguendo ora questo tipo di distribuzione si ottiene il grafico infig.2.5. In questo caso si può apprezzare il fatto che la maggiore concentrazione di

0 0.2 0.4 0.6 0.8

1 1.2 1.4 1.6 1.8

2

Figura 2.5: Grafico di valori random con distribuzione non uniforme.

punti si abbia in prossimità del valore 0, come era stato previsto.

18


2.2.1 Mappatura a canali

Seguendo la densità di probabilità appena descritta per la distribuzione di valo-ri random, si organizzano ora questi risultati in vari sottoinsiemi dell’intervallo dipartenza, detti appunto canali

20jj-1

Figura 2.6: Divisione dell’intervallo in canali

Il numero di canali nc è arbitrario e in questo caso è stato preso nc = 100 (infig.2.6 ne è mostrato un numero indicativo).Il canale j − esimo lo si può indicare come

ncj = [2

nc(j − 1), 2

ncj] j = 1, ..., nc

quindi il valore random ξi può disporsi in uno specifico canale a seconda della suaparte intera e decimale (come tutti gli altri che ricadono nello stesso range). Si haquindi un’operazione di conteggio dei casi che ricadono in ciascun intervallo che puòessere rappresentata da degli istogrammi.Eseguendo il programma per 10.000 dati si osserva il grafico in fig.2.7 nel quale sipuò apprezzare una certa variazione statistica rispetto alla linea di distribuzioneteorica rappresentata in verde.Questo problema si può in parte risolvere aumentando il numero di valori generatidal software, ad esempio 1.000.000, e in questo caso (fig.2.8) si può notare il fattoche i valori ora si raggruppano molto meglio attorno alla distribuzione prevista.

19


0

0.2

0.4

0.6

0.8

1

1.2

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

Figura 2.7: Grafico dello schema a canali per 10.000 valori con distribuzione nonuniforme.

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9

1

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

Figura 2.8: Grafico dello schema a canali per 1.000.000 valori con distribuzione nonuniforme.

20

Capitolo 3

Simulazione dello spettro β−

3.1 Elementi radioattivi considerati: 9039Y51 e9038Sr52

In questo capitolo verrà analizzato il problema esposto inizialmente, ovvero dellapossibile presenza di isotopi contaminanti all’interno di campioni di sostanze uti-lizzate per cure oncologiche. Il campione che si prende in esame è composto inparticolare da 9039Y51 (un isotopo dell’ittrio con 39 protoni e 51 neutroni) e da

9038Sr52

(un isotopo dello stronzio con 38 protoni e 52 neutroni). Il primo elemento ha unaapplicazione medica notevole, infatti è utilizzato sia per cure radioterapiche sia peril trattamento del carcinoma epatocellulare, il tipo più comune dei tumori al fegatoche è una delle principali cause di morte nelle persone con cirrosi. 90Y viene utiliz-zato nella radioterapia interna selettiva che consiste nell’iniettare nelle arterieche conducono all’organo infetto delle microsfere contenenti proprio questo elementoche si va poi a depositare nelle cavità del fegato.

Per quanto riguarda il nostro caso, consideriamo che nel campione sia presente lo9038Sr52 come contaminante. Lo

9038Sr52 è il più importante isotopo radioattivo dello

stronzio, viene creato per scopi industriali attraverso la fissione nucleare e si trovaanche in natura a causa dell’uomo essendo uno dei componenti del fallout nuclearedovuto ai test atomici atmosferici compiuti tra gli anni ’50 e ’60. Poiché è fruttodella fissione nucleare è contenuto anche negli scarti del combustibile dei reattorinucleari e grandi quantità di questo elemento sono state rilasciate nell’ambiente du-rante gli incidenti alle centrali nucleari di Chernobyl nel 1986 e di Fukushima nel2011.Questo elemento può essere inalato, ma le maggiori preoccupazioni per la salutederivano dall’ingestione. Una volta nell’organismo esso va a sostituirsi al calcio eviene incorporato nei denti e nelle ossa dove può causare forme cancerose sia neltessuto osseo che in quello circostante. Questa sostituzione avviene perché la confi-gurazione elettronica esterna, cioè quella che crea i legami tra gli elementi, è la stessa.

Entrambi radioisotopi in questione subiscono un procedimento di decadimento β−

21

Capitolo 3. Simulazione dello spettro β−

ma si differiscono in particolare per i tempi di emivita e per le energie rilasciate.Infatti 9039Y51 ha un tempo di dimezzamento di 64.053 ore mentre

9038Sr52 dimezza i

suoi radionuclidi in 28.9 anni. 9038Sr52 decade nel9039Y51 con

Q− valore(90Sr→90Y ) = 545.96 keV

mentre 9039Y51 decade nello stato fondamentale di9040Zr50 (un isotopo dello zirconio

considerato stabile) con

Q− valore(90Y→90Zr) = 2278.5 keV

Questi decadimenti possono essere visualizzati in fig.3.1

Figura 3.1: Decadimenti di 9038Sr52 e di9039Y51

3.2 Distribuzione di probabilità

Come si è visto nel Capitolo 2, data una certa distribuzione di probabilità, generan-do dei dati casuali questi andranno a disporsi secondo tale relazione.Attraverso un programma software in linguaggio Fortran 90, messo a punto con ilrelatore di questo elaborato finale, è possibile simulare la distribuzione delle proba-bilità con variabile l’energia emessa dal decadimento. Il programma per calcolare ladistribuzione degli eventi secondo le energie ha bisogno di alcuni parametri che sono

• momento angolare totale sia del nucleo padre che del nucleo figlio;

• numeri relativistici delle shell di occupazione dei protoni e degli elettroni;

• Q-valore dei processi;

22


• numero di massa degli elementi.

Il programma considera anche l’interazione dell’elettrone emesso con quelli che occu-pano gli orbitali atomici attraverso il metodo di Hartee-Fock e per un numero fissatodi eventi si ottengono i grafici mostrati in fig.3.2 e in fig.3.3. Il codice utilizzato è ilseguente

program distrib

implicit none

integer(kind=8) i,ne

integer j,np,nc

double precision, allocatable :: ep(:),sp(:),pp(:),ec(:)

double precision x,y,em,f

external f

read(5,*) np,nc,em

read(5,’(i15)’) ne

allocate(ep(np),sp(np),pp(np),ec(nc))

do i=1,np

read(5,*) ep(i),sp(i)

end do

! Genera la funzione integrale

pp(1)=0

do i=2,np

pp(i)=pp(i-1)+(ep(i)-ep(i-1))*(sp(i)+sp(i-1))/2

end do

! Rinormalizza la funzione

sp=sp/pp(np)

pp=pp/pp(np)

open(10,file=’sp’)

do i=1,np

write(10,*) ep(i),sp(i),pp(i)

end do

close(10)

ec=0

do i=1,ne

call random_number(y)

x=f(y,ep,pp,np)

write(8,*) x

j=int(x/em*nc)+1

ec(j)=ec(j)+1

end do

open(9,file=’out’)

do j=1,nc

write(9,*)

(2*j-1.d0)*em/(2*nc),ec(j)*nc/(em*ne),nint(ec(j)),

(f((j-0.d0)*em/nc,pp,ep,np)-f((j-1.d0)*em/nc,pp,ep,np))*ne

end do

close(9)

deallocate(ep,sp,pp,ec)

23


end program

double precision function f(y,ep,pp,np)

implicit none

integer i,ia,ib,ii,np

double precision y,ep(np),pp(np)

ia=1

ib=np

if (y.lt.pp(ia)) then

write(6,*) pp(ia),’>’,y

f=ep(ia)

return

end if

if (y.gt.pp(ib)) then

write(6,*) y,’>’,pp(ib)

f=ep(ib)

return

end if

do while (ib-ia.gt.1)

ii=(ib+ia)/2

if (pp(ii).gt.y) then

ib=ii

else

ia=ii

endif

end do

f=(y-pp(ia))/(pp(ib)-pp(ia))*(ep(ib)-ep(ia))+ep(ia)

end function

In questo caso si è utilizzato lo schema a canali visto nel Capitolo 2.Il programma che si utilizza genera una serie di punti ricavati in modo non analitico cheformano una rappresentazione discreta della funzione primitiva della densità di probabilità.Attraverso il software Octave è possibile interpolare linearmente questi valori per avere lafunzione primitiva. Dalla relazione (2.13) si può cos̀ı generare la funzione densità di proba-bilità teorica che è mostrata nei grafici in fig.3.2 e in fig.3.3. Si nota come nel caso di 90Y lafluttuazione intorno alla distribuzione teorica sia molto meno apprezzabile rispetto al casodi 90Sr dovuto alla grande differenza del numero di valori generati (come si vedrà tra poco).Come nell’esempio mostrato nel Capitolo 2 è ora possibile generare un numero arbitrariodi eventi che andranno a disporsi secondo questo schema probabilistico. Tale procedimentogenera comunque un errore dovuto all’approssimazione ma può essere trascurato dato chei valori su cui si va ad interpolare sono 1000 e la funzione risulta ben riprodotta. La di-stribuzione di probabilità ora quindi è ricavata e non data per nota come nel caso visto nelcapitolo precedente.I due grafici non rappresentano lo spettro di decadimento reale degli elementi analizzatiperché l’integrale su tutto l’intervallo è 1 e non rappresenta il rate di decadimento comenegli spettri sperimentali.

24


0

0.0001

0.0002

0.0003

0.0004

0.0005

0.0006

0.0007

0 500 1000 1500 2000 2500

Freq

uenc

y

Electron energy (keV)

90Y

FrequencyDensity probability

Figura 3.2: Distribuzione di probabilità per 90Y.

0

0.0005

0.001

0.0015

0.002

0.0025

0.003

0.0035

0 500 1000 1500 2000 2500

Freq

uenc

y


90Sr

FrequencyDensity probability

Figura 3.3: Distribuzione di probabilità per 90Sr.

3.3 Campionamento dei valori generati

È noto che un certo evento ha la probabilità di manifestarsi pari a

prob =#casifavorevoli

#eventi

25


Conoscendo quindi la probabilità (dai grafici) e il numero di eventi (impostato dall’utente delprogramma) si possono ottenere i casi favorevoli per ogni singolo canale, ovvero i decadimentiche avvengono per un certo range di energia, che in questo caso sarebbe 10 keV dato che ilfondo scala del grafico è 2290 keV ed è diviso in 229 canali. Conoscendo il numero di eventitotali si può ottenere il numero degli eventi favorevoli, cioè quelli compresi in un singolocanale moltiplicando l’area dell’istogramma per l’ampiezza dell’intervallo e per il numero dicasi totali generati

N = A ·∆E ·Ntot (3.1)

Ad esempio in fig.3.2 si consideri il canale che ha densità di probabilità di 6·10−4 e unnumero totale di eventi pari a 399 milioni, si ottiene che

N = 0.0006 · 10 · 399.000.000 w 2.394.000

Per 9038Sr52 sono stati scelti 100.000 eventi, mentre per9039Y51 399.000.000 in quanto il numero

di eventi è inversamente proporzionale alla vita media e il rapporto tra le due vite medieè circa 3995, cioè la vita media di 9038Sr52 è circa 3995 volte quella di

9039Y51. Quindi il rate

dell’ittrio è 3995 volte più alto di quello dello stronzio e se avessimo una quantità di atomiuguali, per ogni atomo di 90Sr che decade, ne decadrebbero mediamente 3995 di 90Y. Su100.000 eventi del primo ne vengono circa 399.500.000 per il secondo.Da qui si possono fare delle stime: se ad esempio si avesse 1 ng di 90Y, dividendo per lamassa del protone 1.67 · 10−27g e per il numero di nucleoni 90, si otterrebbe che il numerodi atomi nel campione è approssimativamente 6.7 · 1015.

Il tempo di vita di 9039Y51 si può ricavare dalla (1.3) ed è circa 332.396 s. Facendo l’in-verso del tempo di vita si ottiene il rate di decadimento. Moltiplicando questo valore peril numero di atomi si ottiene il numero di decadimenti in 1 s, che per l’ittrio è circa 20miliardi al secondo. Ora il numero di eventi che è stato generato dal simulatore è moltominore (circa 2100 ) di quelli che avvengono in un secondo. Per ovviare a questo problema siconsidera dunque un campione inferiore a 1 ng, ad esempio 1 pg e considerando un tempodi osservazione di 20 s, per 90Y si avranno circa 400 milioni di eventi. Per 9038Sr52 si fa unacosa analoga, ottenendo circa 2 milioni di eventi.

Una volta ottenuto il numero di eventi per ogni canale, si possono sommare per ottene-re lo spettro completo dei conteggi per i singoli elementi considerati separatamente, comemostrato in fig.3.4 e in fig.3.5.

26


0

500000

1x106

1.5x106

2x106

2.5x106

0 500 1000 1500 2000 2500

Coun

ts


90Sr

Figura 3.4: Conteggi per 90Y.

0

500

1000

1500

2000

2500

3000

0 500 1000 1500 2000 2500

Coun

ts


90Sr

Figura 3.5: Conteggi per 90Sr.

3.4 Fit con dati sperimentali

Considerando un campione da analizzare composto da i due elementi secondo la relazionedi massa 1 : 1 (1 pg di 90Y e 1 pg di 90Sr), si può generare uno spettro di conteggi di questocampione.La quantità S corrisponde a χ2 e si ha l’obbiettivo ora di minimizzarla per ottenere i valori

27


migliori di a e b.nc∑j=1

[(a · fYj + b · fSrj )− Fj

]2Fj

= S (3.2)

con fYj frequenza teorica di conteggi per il canale j−esimo relativamente a 90Y consideratosingolarmente, fSrj frequenza teorica per il canale j−esimo relativamente a 90Sr consideratosingolarmente, Fj frequenza sperimentale per il canale j − esimo del campione compostodalla somma di entrambi gli elementi. La differenza tra la frequenza dei dati generati(sperimentali) e quella dei dati teorici può essere apprezzata nei grafici in fig.3.6 e 3.7,mentre il conteggio degli eventi generati per il campione composto è mostrato in figura 3.8.La sommatoria viene pesata sulle frequenze teoriche invece che su quelle sperimentali persemplificare il conteggio dato che l’errore di si commette non è rilevante. In questi casi siè imposto un numero di conteggi di circa 400 milioni per l’ittrio e di circa 100.000 per lostronzio.

0

500000

1x106

1.5x106

2x106

2.5x106

0 500 1000 1500 2000

Coun

ts


90Sr

Exp.Theor.

Figura 3.6: Differenza tra frequenze teoriche e quelle sperimentali per 90Y.

28


0

500

1000

1500

2000

2500

3000

0 100 200 300 400 500

Coun

ts


90Sr

Exp.Theor.

Figura 3.7: Differenza tra frequenze teoriche e quelle sperimentali per 90Sr.

0

500000

1x106

1.5x106

2x106

2.5x106

0 500 1000 1500 2000

Coun

ts


90Sr

Exp.

Figura 3.8: Conteggi totali sperimentali del campione composto da entrambi glielementi.

29


Al fine di ottenere i valori di a e b che minimizzino quantità S sia rispetto all’ittrio cheallo stronzio si esegue la derivata prima e si pone uguale a 0 come in (3.3).

�2∑

j

(a · fYj + b · fSrj − Fj

F 2j

)· fYj = 0

�2∑

j

(a · fYj + b · fSrj − Fj

F 2j

)· fSrj = 0

(3.3)

Le frequenze teoriche in questione sono state ottenute moltiplicando le aree con base lalarghezza di ogni intervallo moltiplicato per il numero di eventi generati. Il sistema (3.3) siriduce a risolvere l’equazione

Cx = B

con C e B matrici di coefficienti sommate su j.

C =∑j

[fYj f

Yj f

Srj f

Yj

fYj fSrj f

Srj f

Srj

](3.4)

B =∑j

[fYj FjfSrj Fj

](3.5)

quindi x = C−1B, da cui si ottiene che i coefficienti a e b della (3.2) risultano

a = 1.00002

b = 0.92791

Si osserva che la quantità di 90Y che ritorna dal fit si discosta da quella teorica di 0.002%,mentre la quantità di 90Sr si discosta da quella prevista del 7.2% circa.Considerano ora un rapporto di quantità 1 : 100 tra ittrio e stronzio si ottengono i valori

a = 1.0000

b = −6.2003

ovvero un’incertezza su 90Y di circa lo 0% e su 90Sr di circa il 6.2%.

Un metodo per ridurre questa incertezza è aumentare il campionamento dei valori di 10volte, passando dai 400 milioni per 90Y a 4 miliardi e da 100 mila per 90Sr a 1 milione.I valori ottenuti in questo caso sono

a = 0.9999

b = 1.01979

da cui risulta un errore di circa lo 0% per 90Y e del 2% per 90Sr. In questo caso si ottengonorisultati simili considerando un rapporto di quantità tra ittrio e stronzio di 1 : 100 ottenendoun errore di circa 0% per l’ittrio e del 3% per lo stronzio.

30

Conclusioni

Il procedimento discusso in questo lavoro è utile quindi a stabilire con un buon grado diapprossimazione il contenuto di un campione composto da due elementi che sono soggettial decadimento β−, ma con le opportune modifiche può essere utilizzato per composti piùcomplessi. Come si è visto i risultati sono migliori all’aumentare del numero di eventi ge-nerati ottenendo una maggiore precisione sul risultato finale, ma questo comporta una piùelevata potenza di calcolo richiesta al simulatore.

Il discorso fatto fino ad ora è valido quindi per una misura simulata attraverso un software,mentre a livello sperimentale possono sorgere diversi problemi. Alcuni di questi possonoessere il rumore dell’elettronica dello strumento con cui si esegue il conteggio degli eventi,la possibilità che un evento si disponga anche in un altro canale attiguo a quello reale conuna certa frequenza, rumore statistico e altri eventi, come ad esempio un conteggio falsoa causa di un segnale elettrico che parte dall’interno della macchina misuratrice causandoquindi un falso conteggio. Questi elementi andrebbero presi in considerazione nell’applicareil procedimento e bisognerebbe considerarli come fattori di correzione. Inoltre, per avere unmaggior numero di dati occorre far durare l’esperimento per un tempo prefissato che puòessere anche molto lungo.

31

Bibliografia

• E.J. Konopinski, The Theory of Beta Radioactivity, Oxford University Press, 1950;

• Paolo Baldi, Calcolo delle probabilità e statistica, Milano, McGraw-Hill, 1998;

• EPA, United States Environmental Protection Agency, Radionuclide Basics: Strontium-90https://www.epa.gov/radiation/radionuclide-basics-strontium-90;

• Nudat 2, sito con tavola dei nuclidi e dati sperimentali

– https://www.nndc.bnl.gov/nudat2/reCenter.jsp?z=39&n=51



• Arkadiy Kheyfits,Yttrium-90 Radioembolization, Radiology Today, Vol. 11 No. 9 P.20

https://www.radiologytoday.net/archive/rt0910p20.shtml

• Joseph Ralph Kallini, Ahmed Gabr, Riad Salem, and Robert J. Lewandowski, Tran-sarterial Radioembolization with Yttrium-90 for the Treatment of Hepatocellular Car-cinoma, 2 aprile 2016

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4882351/

Sommario si ha Q = m nc 2 m pc 2 m e c 2 (1.6) 1.3.2 Decadimento + Nel seguente processo avviene un...

Documents

Transcript of Sommario si ha Q = m nc 2 m pc 2 m e c 2 (1.6) 1.3.2 Decadimento + Nel seguente processo avviene un...