Dati germanio¶

Dati raccolti col germanio del professor Clemenza, dell'università di Milano Bicocca. Sono stati raccolti campioni di Torio (elettrodi) e di fondo.

In [1]:
import struct
import os
import sys

import numpy as np
from matplotlib import pyplot as plt

from scipy.optimize import curve_fit
In [2]:
myEvent = np.dtype([
    ('timestamp', np.dtype('u8')),
    ('qshort', np.dtype('u2')),
    ('qlong', np.dtype('u2')),
    ('baseline', np.dtype('u2')),
    ('channel', np.dtype('u1')),
    ('PUR', np.dtype('u1')),
])

$^{232}$Th¶

In [3]:
with open("232Th_1_events.ade", "rb") as f:
    myData = np.fromfile(f, dtype = myEvent)
In [4]:
%matplotlib inline

myBin = np.arange(-.5, 60000.5, 1)

h, bins = np.histogram(myData["qlong"], bins = myBin)
binc = bins[:-1] + (bins[1] - bins[0])/2

fig, ax = plt.subplots()
fig.set_size_inches(12,5)

ax.plot(binc, h, ds = "steps-mid", c = "darkgreen", lw = 1, label = "Spettro")
ax.fill_between(binc, h, step = "mid", color = "lime", alpha = 1)



ax.grid()
ax.set_xlabel("Energia [ADC]", fontsize = 14)
ax.set_ylabel("Entries", fontsize = 14)
ax.set_title("$^{232}Th$ - Elettrodi", fontsize = 16)
ax.legend()
ax.set_yscale("log")

plt.show()
In [5]:
%matplotlib inline
#%matplotlib qt


h, bins = np.histogram(myData["qlong"], bins = 1000)
binc = bins[:-1] + (bins[1] - bins[0])/2

fig, ax = plt.subplots()
fig.set_size_inches(12,5)

ax.plot(binc, h, ds = "steps-mid", c = "darkgreen", lw = 1, label = "Spettro")
ax.fill_between(binc, h, step = "mid", color = "lime", alpha = 1)



ax.grid()
ax.set_xlabel("Energia [ADC]", fontsize = 14)
ax.set_ylabel("Entries", fontsize = 14)
ax.set_title("$^{232}Th$ - Elettrodi", fontsize = 16)
ax.legend()
ax.set_yscale("log")

plt.show()

Calibrazione rapida¶

Uso solo i primi 9 picchi di questa figura, degli altri non sono sicuro

In [6]:
%matplotlib inline
# Mi ricreo l'istogramma con i bin desiderati (ora quello sorpa)
h, bins = np.histogram(myData["qlong"], bins = 1000)
binc = bins[:-1] + (bins[1] - bins[0])/2


# Definisco l'intervallo in cui cercare il max per una rapida calibrazione
intervalli = ((7700, 9000), (11300, 12400), (17500, 19260), (19260, 21430), (25000, 26460), 
             (27960, 28840), (30000, 31430), (31430, 33330), (33330, 35910), )

energieVere = np.array((238, 338, 510, 583, 727, 794, 860, 911, 968))

# Lista che conterrà x e y dei picchi
xmax = []
ymax = []
for i in intervalli:
    tmpCond = (binc >= i[0]) & (binc <= i[1])
    tmpXdata = binc[tmpCond]
    tmpYdata = h[tmpCond]
    
    # Indice del max sui sotto vettori appena costruiti
    tmpIdx = np.argmax(tmpYdata)
    
    xmax.append(tmpXdata[tmpIdx])
    ymax.append(tmpYdata[tmpIdx])
    
xmax = np.array(xmax)
ymax = np.array(ymax)



# Plotto spettro e picchi
fig, ax = plt.subplots()
fig.set_size_inches(12,5)

ax.plot(binc, h, ds = "steps-mid", c = "darkgreen", lw = 1, label = "Spettro")
ax.fill_between(binc, h, step = "mid", color = "lime", alpha = 1)
ax.plot(xmax, ymax, marker = "*", c = "tab:orange", ls = "", label = "Found peaks")


ax.grid()
ax.set_xlabel("Energia [ADC]", fontsize = 14)
ax.set_ylabel("Entries", fontsize = 14)
ax.set_title("$^{232}Th$ - Elettrodi", fontsize = 16)
ax.legend()
ax.set_yscale("log")

plt.show()

Retta di calibrazione¶

In [7]:
# Retta di calibrazione
def myLine(x,m,q):
    return m*x+q

popt, pcov = curve_fit(myLine, energieVere, xmax)


fig, ax = plt.subplots()
ax.plot(energieVere, xmax, "*g", label = "Dati")
ax.plot(energieVere, myLine(energieVere, *popt), ":k", label = "Retta fittata")

ax.set_title("Retta di calibrazione", fontsize = 16)
ax.set_xlabel("Energia vera [keV]", fontsize = 14)
ax.set_ylabel("Energia [ADC]", fontsize = 14)

plt.legend(fontsize = 14)
plt.grid()
plt.show()
print(popt)

[35.63868843 31.91914223]

Spettro calibrato¶

In [8]:
# Plotto lo spettro calibrato
m = 1/popt[0]
q = -popt[1]/popt[0]
binckev = myLine(binc, m, q)


fig, ax = plt.subplots()
fig.set_size_inches(12,5)

ax.plot(binckev, h, ds = "steps-mid", c = "darkgreen", lw = 1, label = "Spettro")
ax.fill_between(binckev, h, step = "mid", color = "lime", alpha = 1)
#ax.plot(xmax, ymax, marker = "*", c = "tab:orange", ls = "", label = "Found peaks")


ax.grid()
ax.set_xlabel("Energia [keV]", fontsize = 14)
ax.set_ylabel("Entries", fontsize = 14)
ax.set_title("$^{232}Th$ - Elettrodi", fontsize = 16)
ax.legend()
ax.set_yscale("log")

plt.show()

Calibrazione con fit¶

In accordo con i sacri testi, i picchi sono

Io utilizzo (238, 338, 510, 583, 727, 794, 860, 911, 968, 1460)

In [9]:
%matplotlib inline
#%matplotlib qt

# Mi ricreo l'istogramma con i bin desiderati (ora quello sorpa)
h, bins = np.histogram(myData["qlong"], bins = 1000)
binc = bins[:-1] + (bins[1] - bins[0])/2


# Definisco l'intervallo in cui cercare il max per una rapida calibrazione
intervalli = ((7700, 9000), (11300, 12400), (17500, 19260), (19260, 21430), (25000, 26460), 
             (27960, 28840), (30000, 31430), (31430, 33330), (33330, 35910), (51679, 52277))

energieVere = np.array((238, 338, 510, 583, 727, 794, 860, 911, 968, 1460))

def myGauss(x,a,mu,sigma,m,q):
    return a*np.exp(-(x-mu)**2 / (2*sigma**2)) + m*x+q


# Plotto spettro e picchi
fig, ax = plt.subplots()
fig.set_size_inches(12,5)

ax.plot(binc, h, ds = "steps-mid", c = "darkgreen", lw = 1, label = "Spettro")
ax.fill_between(binc, h, step = "mid", color = "lime", alpha = 1)


xmax1 = []
for i in intervalli:
    tmpCond = (binc >= i[0]) & (binc <= i[1])
    tmpXdata = binc[tmpCond]
    tmpYdata = h[tmpCond]
    
    # Indice del max sui sotto vettori appena costruiti
    tmpIdx = np.argmax(tmpYdata)
    
    xmax = tmpXdata[tmpIdx]
    ymax = tmpYdata[tmpIdx]
    
    
    
    # fitto
    popt, pcov = curve_fit(myGauss, tmpXdata, tmpYdata, 
                           sigma = np.sqrt(tmpXdata), absolute_sigma = True, 
                           p0 = (ymax, xmax, 50, -.03, 2000))
    
    #print(popt)
    ax.plot(tmpXdata, myGauss(tmpXdata, *popt), ":k")
    xmax1.append(popt[1])
    







ax.grid()
ax.set_xlabel("Energia [ADC]", fontsize = 14)
ax.set_ylabel("Entries", fontsize = 14)
ax.set_title("$^{232}Th$ - Elettrodi", fontsize = 16)
ax.legend()
ax.set_yscale("log")
ax.set_ylim((20, 1e4))

plt.show()

xmax1 = np.array(xmax1)

Retta di calibrazione¶

In [10]:
# Retta di calibrazione
def myLine(x,m,q):
    return m*x+q

popt, pcov = curve_fit(myLine, energieVere, xmax1, p0 = (35, 31))


fig, ax = plt.subplots()
ax.plot(energieVere, xmax1, "*g", label = "Dati")
ax.plot(energieVere, myLine(energieVere, *popt), ":k", label = "Retta fittata")

ax.set_title("Retta di calibrazione", fontsize = 16)
ax.set_xlabel("Energia vera [keV]", fontsize = 14)
ax.set_ylabel("Energia [ADC]", fontsize = 14)

plt.legend(fontsize = 14)
plt.grid()
plt.show()

Verifica non linearità - Residui¶

Plotto $\dfrac{E_\text{fit}-E_\text{vera}}{E_\text{vera}}$

In [11]:
m = 1/popt[0]
q = -popt[1]/popt[0]


residui = (myLine(xmax1, m, q) - energieVere) / energieVere

fig, ax = plt.subplots()

ax.plot(energieVere, 100*residui, "*g")

ax.set_title("Residui percentuali", fontsize = 16)
ax.set_xlabel("Energia vera [keV]", fontsize = 14)
ax.set_ylabel("Residui [%]", fontsize = 14)

ax.grid()

plt.show()

Spettro calibrato¶

In [12]:
# Plotto lo spettro calibrato

binckev = myLine(binc, m, q)


fig, ax = plt.subplots()
fig.set_size_inches(12,5)

ax.plot(binckev, h, ds = "steps-mid", c = "darkgreen", lw = 1, label = "Spettro")
ax.fill_between(binckev, h, step = "mid", color = "lime", alpha = 1)
#ax.plot(xmax, ymax, marker = "*", c = "tab:orange", ls = "", label = "Found peaks")


ax.grid()
ax.set_xlabel("Energia [keV]", fontsize = 14)
ax.set_ylabel("Entries", fontsize = 14)
ax.set_title("$^{232}Th$ - Elettrodi", fontsize = 16)
ax.legend()
ax.set_yscale("log")

plt.show()

Per poster - [FOTO SPETTRO GERMANIO]¶

In [17]:
# Plotto lo spettro calibrato

binckev = myLine(binc, m, q)


fig, ax = plt.subplots()
fig.set_size_inches(12,5)

ax.plot(binckev[20:-1], h[20:-1], ds = "steps-mid", c = "darkgreen", lw = 1, label = "Spectrum")
ax.fill_between(binckev[20:-1], h[20:-1], step = "mid", color = "lime", alpha = 1)
#ax.plot(xmax, ymax, marker = "*", c = "tab:orange", ls = "", label = "Found peaks")


ax.grid()
ax.set_xlabel("Energy [keV]", fontsize = 14)
ax.set_ylabel("Entries", fontsize = 14)
ax.set_title("$^{232}$Th spectrum ", fontsize = 16)
ax.legend()
ax.set_yscale("log")

plt.show()

Misura della risoluzione¶

In [13]:
%matplotlib inline

# Creo l'istogramma a step di 1ADC
myBin = np.arange(-.5, 60000.5, 1)

h, bins = np.histogram(myData["qlong"], bins = myBin)
binc = bins[:-1] + (bins[1] - bins[0])/2
binckev = myLine(binc, m, q) # binc in kev, estratta dal fit


lstRes = []

for i in energieVere:
    mu = i
    sigma = 8
    
  
    cond = (binckev >= (mu-sigma)) &  (binckev <= (mu+sigma)) & (h>0)
    xData = binckev[cond]
    yData = h[cond]
    popt, pcov = curve_fit(myGauss, xData, yData, sigma = np.sqrt(yData), 
                           absolute_sigma = True, maxfev = 1000000,
                           p0 = (np.max(yData), mu, 2, -.03, 2000),)




    fig, ax = plt.subplots()
    fig.set_size_inches(12,5)

    ax.plot(binckev, h, ds = "steps-mid", c = "darkgreen", lw = 1, label = "Spettro")
    ax.fill_between(binckev, h, step = "mid", color = "lime", alpha = 1)
    ax.plot(xData, myGauss(xData, *popt), ":k")


    ax.grid()
    ax.set_xlabel("Energia [keV]", fontsize = 14)
    ax.set_ylabel("Entries", fontsize = 14)
    ax.set_title(f"Risoluzione a {mu} keV", fontsize = 16)
    ax.legend()
    ax.set_yscale("log")


    ax.set_xlim((mu-sigma, mu+sigma))

    plt.show()
    
    tmpRes = 235*np.abs(popt[2])/popt[1]
    lstRes.append(tmpRes)

    print(f"La risoluzione per il picco a {mu} keV vale {tmpRes:.4f} %")
    #print(popt)
    
lstRes = np.array(lstRes)


for mu, tmpRes in zip(energieVere, lstRes):
    print(f"La risoluzione per il picco a {mu} keV vale\t {tmpRes:.4f} %")

La risoluzione per il picco a 238 keV vale 0.9870 %

La risoluzione per il picco a 338 keV vale 0.6837 %

La risoluzione per il picco a 510 keV vale 0.5532 %

La risoluzione per il picco a 583 keV vale 0.4353 %

La risoluzione per il picco a 727 keV vale 0.3467 %

La risoluzione per il picco a 794 keV vale 0.3392 %

La risoluzione per il picco a 860 keV vale 0.3007 %

La risoluzione per il picco a 911 keV vale 0.3017 %

La risoluzione per il picco a 968 keV vale 0.2932 %

La risoluzione per il picco a 1460 keV vale 0.2141 %
La risoluzione per il picco a 238 keV vale	 0.9870 %
La risoluzione per il picco a 338 keV vale	 0.6837 %
La risoluzione per il picco a 510 keV vale	 0.5532 %
La risoluzione per il picco a 583 keV vale	 0.4353 %
La risoluzione per il picco a 727 keV vale	 0.3467 %
La risoluzione per il picco a 794 keV vale	 0.3392 %
La risoluzione per il picco a 860 keV vale	 0.3007 %
La risoluzione per il picco a 911 keV vale	 0.3017 %
La risoluzione per il picco a 968 keV vale	 0.2932 %
La risoluzione per il picco a 1460 keV vale	 0.2141 %
In [14]:
#voglio un errore

Fondo¶

In [15]:
with open("Fondo_events.ade", "rb") as f:
    myData = np.fromfile(f, dtype = myEvent)
In [16]:
%matplotlib inline

myBin = np.arange(-.5, 60000.5, 1)

h, bins = np.histogram(myData["qlong"], bins = myBin)
binc = bins[:-1] + (bins[1] - bins[0])/2

fig, ax = plt.subplots()
fig.set_size_inches(12,5)

ax.plot(binc, h, ds = "steps-mid", c = "darkgreen", lw = 1, label = "Spettro")
ax.fill_between(binc, h, step = "mid", color = "lime", alpha = 1)



ax.grid()
ax.set_xlabel("Energia [ADC]", fontsize = 14)
ax.set_ylabel("Entries", fontsize = 14)
ax.set_title(f"Fondo", fontsize = 16)
ax.legend()
ax.set_yscale("log")

plt.show()
In [17]:
%matplotlib inline

h, bins = np.histogram(myData["qlong"], bins = 1000)
binc = bins[:-1] + (bins[1] - bins[0])/2

fig, ax = plt.subplots()
fig.set_size_inches(12,5)

ax.plot(binc, h, ds = "steps-mid", c = "darkgreen", lw = 1, label = "Spettro")
ax.fill_between(binc, h, step = "mid", color = "lime", alpha = 1)



ax.grid()
ax.set_xlabel("Energia [ADC]", fontsize = 14)
ax.set_ylabel("Entries", fontsize = 14)
ax.set_title(f"Fondo", fontsize = 16)
ax.legend()
ax.set_yscale("log")

plt.show()