Geiger 2

Programma base

Analizza tutti i file (restanti) di dati a diverse distanze, trova i conteggi e fitta i picchi

%matplotlib inline

import numpy as np
import matplotlib.pyplot as plt
import scipy.optimize as opt

def gaus(x, A, mu, sigma):
    return A*np.exp(-(x-mu)**2/(2*sigma**2))

# prepara per memorizzare i risultati
mean  = []
mean_err = []

# figura
fig, axs = plt.subplots(4,4, sharex=True, figsize=(15,10))

# init fit
par0 = [20, 100,10]
bin0 = 0
for i in range(31,47):
    datafile = open('dati/run%6.6d.dat' % (i))
    a,b,c,counts,time = np.loadtxt(datafile,unpack=True)
    rate = 1000*counts/time
    h, bins, _ = axs.reshape(16)[i-31].hist(rate,190, range=(0,800), histtype='step', color='k', linewidth=2, label='Run %d'%i)
    p, cov = opt.curve_fit(gaus,bins[bin0:-1],h[bin0:],p0=[100,rate.mean(),rate.mean()*0.1] )
    #par0 = p
    errp = np.sqrt(np.diag(cov))
    axs.reshape(16)[i-31].plot(bins[bin0:],gaus(bins[bin0:],*p),'--r', linewidth=2, label='fit')
    axs.reshape(16)[i-31].legend()
    mean.append(p[1])
    mean_err.append(errp[1])
    
plt.show()

Continuano i fit (file con range un po' diversi)

# figura
fig, axs = plt.subplots(1,3, sharex=True, figsize=(15,4))

# init fit
par0 = [20, 100,10]
bin0 = 0
for i in range(47,50):
    datafile = open('dati/run%6.6d.dat' % (i))
    a,b,c,counts,time = np.loadtxt(datafile,unpack=True)
    rate = 1000*counts/time
    h, bins, _ = axs[i-47].hist(rate,170, range=(700,3200), histtype='step', color='k', linewidth=2, label='Run %d'%i)
    p, cov = opt.curve_fit(gaus,bins[bin0:-1],h[bin0:],p0=[100,rate.mean(),rate.mean()*0.1] )
    #par0 = p
    errp = np.sqrt(np.diag(cov))
    axs[i-47].plot(bins[bin0:],gaus(bins[bin0:],*p),'--r', linewidth=2, label='fit')
    axs[i-47].legend()
    mean.append(p[1])
    mean_err.append(errp[2])
    
plt.show()

Grafico e fit del rate in funzione della distanza

def fun(x,p1,p2):
    return p1/(p2-x)**2


dist = np.array([120,100,80,75,70,65,60,55,50,45,40,35,30,25,20,15,10,5,0])
stp  = np.array([15,20,25,35,40,50,55,60,65,90,110,137,176,250,350,520,840,1400,2900])

plt.errorbar(dist[:12],mean[:12],mean_err[:12], label='dati')
#p = np.polyfit(dist,mean,1)
#plt.plot(dist,p[1]+p[0]*dist,'r--', label='y=%.3f*x+%.3f' % (p[0],p[1]))
fpar, cov = opt.curve_fit(fun,dist[:12],mean[:12], p0=[200000,0.4])
plt.plot(dist[:12],fun(dist[:12],*fpar),'--r', linewidth=2, label='fit con y = %.0f / (%.1f-x)**2' % (fpar[0],fpar[1]))
plt.title('Rate in funzione della distanza')
plt.xlabel('distanza [mm]')
plt.ylabel('Rate [Hz]')
plt.legend()
plt.show()

Grafico e fit del rate in funzione dell'inverso del quadrato della distanza, fit lineare della parte iniziale del plot

invdist_iniz = 1/dist[:12]**2
invdist_fine = 1/dist[13:18]**2
invdist_all = 1/dist[:16]**2

plt.plot(invdist_iniz,mean[:12],lw=3, label='lontano')
plt.plot(invdist_fine,mean[13:18], lw=2, label='vicino')
p = np.polyfit(invdist_iniz,mean[:12],1)
plt.plot(invdist_all,p[0]*invdist_all+p[1],'--r', label='y = %.0f x +%.1f'%(p[0],p[1]))
plt.title('Conteggi in funzione dell\'inverso della distanza', fontsize=14)
plt.xlabel('1/d$^2$ [mm$^{-2}$]', fontsize=14)
plt.ylabel('Conteggi', fontsize=14)
plt.legend(fontsize=16)
plt.show()

Calcolo finale (rate vero vs misurato + fit con funzione che contiene il parametro del tempo morto)

#sigma rtrue=car(1)/((car(2)-dist)**2)
rtrue = fpar[0] / (fpar[1]-dist)**2
plt.errorbar(mean[:18],rtrue[:18], mean_err[:18], label='dati', lw=2)
plt.title('Rate vero vs misurato', fontsize=14)
plt.xlabel('misurato [Hz]', fontsize=14)
plt.ylabel('vero [Hz]', fontsize=14)

def fun2(x,tau):
    return x/(1-x*tau)

p, cov = opt.curve_fit(fun2,mean[:18],rtrue[:18],p0=2e-4)
x=np.linspace(0,1400,100)
plt.plot(x,fun2(x,p),'--r', linewidth=2, label='fit con x/(x-x*%.1e)' % p)
plt.legend(fontsize=16)

plt.show()