Per queste misure sono state effettuate delle single acquisition sull'oscilloscopio, in assenza di luce ed in una finestra lunga 50us, triggerando sul led (che non era collegato, usato solo per il trigger).
Sono quindi stati registrati il numero di conteggi per ciascuna finestra. L'istogramma è stato quindi fittato con una poissoniana
import numpy as np
from matplotlib import pyplot as plt
from scipy.optimize import curve_fit
import scipy
Conteggi in una finestra lunga 50us
counts = np.array([
7, 7, 12, 9, 7,
10, 6, 3, 8, 11,
6, 5, 4, 7, 11,
3, 7, 11, 13, 8,
4, 12, 3, 3, 5,
4, 10, 5, 10, 9,
8, 10, 11, 10, 9,
5, 7, 3, 7, 8,
3, 7, 17, 12, 6,
2, 10, 8, 10, 9,
5, 10, 4, 3, 19,
11, 6, 5, 6, 5,
4, 7, 5, 6, 7,
12, 8, 6, 5, 11
])
# plt.hist(counts, bins=17)
%matplotlib inline
h, bins = np.histogram(counts, bins = np.arange(2,20))
binc = bins[:-1] + (bins[1] - bins[0])/2
fig, ax = plt.subplots()
ax.plot(bins[:-1], h, ds = "steps-mid", c = "tab:green")
ax.errorbar(bins[:-1], h, yerr = np.sqrt(h), c = "k", ls = "")
plt.show()
# Probabilità di avere k conteggi con mu valor medio (da stimare)
binsVect = bins[:-1].astype(int)
def myPoiss(k, mu):
return np.exp(-mu) * (mu **k) / scipy.special.factorial(k)
def myChi2(_bins, _h, _mu):
return np.sqrt(np.sum( ((_h - np.sum(_h) * myPoiss(_bins, _mu))/np.sqrt(_h))**2 )) / (len(_bins)-1)
valoriScan = np.arange(6, 9, .05)
lstchi = []
cond = h>0
for i in valoriScan:
lstchi.append(myChi2(binsVect[cond], h[cond], i))
lstchi = np.array(lstchi)
fig, ax = plt.subplots()
ax.plot(valoriScan, lstchi, "*g")
ax.set_title("Reduced chi2 vs $\mu$", fontsize = 16)
ax.set_xlabel("$\mu$", fontsize = 14)
ax.set_ylabel("Red. chi2", fontsize = 14)
ax.grid()
plt.show()
optMu = valoriScan[np.argmin(lstchi)]
#print(optMu)
print(f"Il valore ottimo di mu vale {optMu:.4f}. Poiché la finestra era lunga 50us, il rate di DCR vale {(optMu/50e-6):.4f} Hz")
print(lstchi.min())
Il valore ottimo di mu vale 7.2500. Poiché la finestra era lunga 50us, il rate di DCR vale 145000.0000 Hz 0.28015898215285395
Per l’errore sul rate di DCR usiamo la formula Rdcr = \dfrac{ \mu{ott} }{ gate } poiché entrambi hanno un errore la propagazione diventa: \err{Rdcr} = \sqrt{ (1/gate \err_{mu_ott} )^2 + (mu_ott / (gate)^2 err_{gate})^2 }
#gio
tmpGate = 50e-6
tmpErrMu = .05
tmpErrGate = .01e-6
tmpDCR = optMu/tmpGate
tmpErrDCR = np.sqrt( (tmpErrMu/tmpGate)**2 + (optMu*tmpErrGate/(tmpGate**2))**2 )
print(f"Il valore ottimo di mu vale {optMu:.4f}. Poiché la finestra era lunga 50us, il rate di DCR vale {tmpDCR:.4f} ± {tmpErrDCR:.4f} Hz")
Il valore ottimo di mu vale 7.2500. Poiché la finestra era lunga 50us, il rate di DCR vale 145000.0000 ± 1000.4204 Hz
%matplotlib inline
h, bins = np.histogram(counts, bins = np.arange(2,20))
binc = bins[:-1] + (bins[1] - bins[0])/2
fig, ax = plt.subplots()
ax.plot(bins[:-1], h, ds = "steps-mid", c = "tab:green", label = "Data")
ax.errorbar(bins[:-1], h, yerr = np.sqrt(h), c = "k", ls = "")
xDenso = np.arange(binsVect.min(), binsVect.max(), 100)
ax.plot(binsVect[cond], np.sum(h) * myPoiss(binsVect[cond], optMu), c = "tab:orange", label = f"Fit poiss - Opt. mu = {optMu:.2f}" )
ax.set_title("Istogramma di conteggi", fontsize = 16)
ax.set_xlabel("# Peaks per window (50 us)", fontsize = 14)
ax.set_ylabel("Entries", fontsize = 14)
ax.grid(True)
ax.legend()
plt.show()
Tempi di interarrivo tra impulsi per ciascuna finestra
times = ((-13.4, -11.35, -9.1, -5.35, -0.2, 0.9, 7.45),
(-11.7, -7.9, -7.1, 1.75, 7.95, 11.7, 21.85),
(-24.5, -23.7, -19.3, -14.75, -6.9, 0.75, 3.2, 4.55, 6.15, 8.05, 13.9, 17.9),
(-13.9, -4.95, 8.15, 11, 13.75, 15.45, 15.6, 16.55, 18.8), #4
(-23.7, -19.6, -18.45, -14, -12.65, -8.45, 24.75),
(-24, -14.05, -12.9, -12.05, -7.85, -4.25, 5.9, 6.85, 16.7, 17.25),
(-5.5, 2.5, 6.2, 11.95, 15.45, 18.2),
(-23.45, -22.05, -14.9), #8
(-16.25, -8.8, -8.4, -2.85, -0.25, 14.55, 16.15, 19.85),
(-19.25, -18.35, -14.65, -8.1, 1.8, 3.05, 3.95, 11.75, 18, 20.60, 24.85),
(-22.65, -11.8, -7.45, 6, 11.8, 23.85),
(-21.1, -19.6, -3.7, -0.35, 0.05), #12
(-7.45, 1.1, 4.05, 17.8),
(-22.9, -17.25, 13.15, 14.95, 17.65, 20.55, 22.5),
(-22.9, -8.35, -6.85, 8.05, 8.4, 8.9, 11.5, 12.6, 14.4, 19.95, 24.25))
lstDeltat = []
for i in times:
tmp = np.array(i)
#lstDeltat = np.vstack((lstDeltat, np.diff(tmp)))
#print(np.diff(tmp))
lstDeltat.append(np.diff(tmp))
lstDeltat = np.hstack(lstDeltat)
#print(lstDeltat.shape, type(lstDeltat[0]), lstDeltat)
def myExp(x, a, tau):
return a * np.exp(-x/tau)
h, bins = np.histogram(lstDeltat, bins = 32, range = (0,32))
binc = bins[:-1] + (bins[1] - bins[0])/2
cond = h>0
xData = bins[:-1][cond]
yData = h[cond]
popt, pcov = curve_fit(myExp, xData, yData, sigma = np.sqrt(yData),
p0 = (20, 32), absolute_sigma = True, maxfev = 1000000)
fig, ax = plt.subplots()
ax.plot(bins[:-1], h, ds = "steps-mid",c = "tab:green", label = "Isto")
ax.errorbar(bins[:-1], h, yerr = np.sqrt(h), c = "k", ls = "")
ax.plot(xData, myExp(xData, *popt), ":k", label = "Fit")
ax.set_title("Tempi di interarrivo", fontsize = 16)
ax.set_xlabel("$\Delta_{time}$", fontsize = 14)
ax.set_ylabel("Entries", fontsize = 14)
ax.grid(True)
ax.legend()
plt.show()
print(f"Il tau fittato vale {popt[1]:.4f}")
myChi2 = np.sqrt(np.sum(((h[cond]-myExp(xData, *popt))/np.sqrt(h[cond]) )**2)) / (len(h[cond]) - len(popt) )
myChi2
Il tau fittato vale 5.6720
0.17221647034088483
h[cond].shape
(16,)
questoTau = popt[1]#*1e-6
newDCR = 1/(questoTau) # tempo in us
print(f"Il DCR calcolato dai tempi di interarrivo vale {newDCR*1e6:.4f} Hz")
questoErrore = np.sqrt(np.diag(pcov)[1])/(questoTau**2)
# 1e6 poichè ho usato il tempo in us
print(f"Il DCR calcolato dai tempi di interarrivo vale {newDCR*1e6:.4f} ± {questoErrore*1e6} Hz")
print(f"Il tau fittato vale {questoTau:.4f} +- {np.sqrt(np.diag(pcov)[1]):.4f} us")
Il DCR calcolato dai tempi di interarrivo vale 176303.8226 Hz Il DCR calcolato dai tempi di interarrivo vale 176303.8226 ± 27578.010656196373 Hz Il tau fittato vale 5.6720 +- 0.8872 us
Z = np.abs(tmpDCR - newDCR*1e6) / np.sqrt(tmpErrDCR**2 + (questoErrore*1e6)**2)
Z
1.1343546928862969
Copio da file "alcuni calcoli"
Testo la compatibilità di quello dell'osci con questi due
lstDCR = np.array((149.07, 150.56, 150.25, 151.39, 151.93, 147.12, 152.85, 151.76, 149.49, 148.97))*1e3
print(lstDCR)
myMean = np.mean(lstDCR)
myStd = np.std(lstDCR)
print (myMean, myStd)
print(f"Il DCR misurato sull'osci vale {myMean:.4f} ± {myStd:.4f}")
[149070. 150560. 150250. 151390. 151930. 147120. 152850. 151760. 149490. 148970.] 150339.0 1630.4689509463221 Il DCR misurato sull'osci vale 150339.0000 ± 1630.4690
Z = np.abs(myMean - newDCR*1e6) / np.sqrt(myStd**2 + (questoErrore*1e6)**2)
Z
0.9398633790471964
Z = np.abs(tmpDCR - myMean) / np.sqrt(tmpErrDCR**2 + myStd**2)
Z
2.7910183592024285