In [1]:
import numpy as np
import matplotlib.pyplot as plt
import scipy.optimize as opt

Definisco un'unica funzione per operare su un certo file. La funzione può ricevere diversi input:

  • il nome del file
  • il valore per cui moltiplicare i valori misurati (perché i dati sono acquisiti con diverse impostazioni di amplificazione)
  • i parametri dei range e di inizializzazione del fit sui 6 picchi
In [2]:
def gaus(x, A, mu, sigma):
    return A*np.exp(-(x-mu)**2/(2*sigma**2))

def fitpicco(x,m,q,A,mu,sigma):
    return m*x+q+gaus(x,A,mu,sigma)

def apri_file_e_fit_picchi(nomefile, ampli, ranges, initpars):
    time, a, b, c, d = np.loadtxt(nomefile, unpack=True)
    # i dati sono nella colonna "a" a meno di un fattore moltplicativo 
    dati = a * ampli
    
    #plt.figure(figsize=(15,5))
    h, bins, _ = plt.hist(dati[(dati<2000*ampli)],400, histtype='step', ec='k')
    hy = np.array(h)
    hx = np.array([(bins[i+1]+bins[i])/2 for i in range(len(bins)-1)])
    # eseguo il fit più volte memorizzando gli output
    optPars_picchi=[]
    pCov_picchi=[]
    plt.close()
    for i in range(6):
        # definisco la zona da fittare
        hxfit = hx[ (hx> ranges[i][0]) & (hx<ranges[i][1]) ]
        hyfit = hy[ (hx>ranges[i][0]) & (hx<ranges[i][1]) ]
        optimizedPars, pCov = opt.curve_fit(gaus,hxfit,hyfit, p0=initpars[i])
        optPars_picchi.append(optimizedPars)
        pCov_picchi.append(pCov)
        #plt.plot(hxfit,gaus(hxfit,*optimizedPars),'-',color='C%d'%i, linewidth=3, label='%d: $\mu$ %.1f, $\sigma$ %.1f'%(i,optimizedPars[1],optimizedPars[2]))
    
    #plt.grid()
    #plt.legend(fontsize=16)
    #plt.show()
    plt.close()
    
    optPars_picchi2=[]
    pCov_picchi2=[]
    plt.figure(figsize=(20,10))
    for i in range(6):
        plt.subplot(2,3,1+i)
        hxfit = hx[ (hx> ranges[i][0]-30) & (hx<ranges[i][1]+30) ]
        hyfit = hy[ (hx>ranges[i][0]-30) & (hx<ranges[i][1]+30) ]
        optimizedPars, pCov = opt.curve_fit(fitpicco,hxfit,hyfit, p0=[0,0]+list(optPars_picchi[i]))
        optPars_picchi2.append(optimizedPars[2:]) # memorizzo solo dal 3 in poi (i primi due sono la retta)
        pCov_picchi2.append(pCov[2:])
        plt.plot(hx,hy,'-k', label='dati')
        plt.plot(hxfit,fitpicco(hxfit,*optimizedPars),'-r', linewidth=3, label='fit ($\mu=%.1f$, $\sigma=%.1f$)'%(optimizedPars[3],np.abs(optimizedPars[4])))
        plt.legend(fontsize=16)
    
    plt.savefig(nomefile.split('/')[-1][:-4]+'_all.pdf')
    plt.close()
    plt.figure(figsize=(10,10))
    picchi_KeV = [77, 185, 241, 295, 351, 609]
    picchi_ADC = np.array(optPars_picchi2)[:,1]
    errori_ADC = np.absolute(np.array(optPars_picchi2)[:,2])
    plt.errorbar(picchi_KeV, picchi_ADC, errori_ADC, marker='o', color='k', label='Dati', linewidth=0.5)
    m, q = np.polyfit(picchi_KeV, picchi_ADC, 1, w=errori_ADC)
    x = np.linspace(picchi_KeV[0], picchi_KeV[-1],100)
    plt.plot(x,m*x+q,'--r',label='Fit y={0:.1f} x + {1:.1f}'.format(m,q), linewidth=2)
    plt.legend()
    plt.savefig(nomefile.split('/')[-1][:-4]+'_fit.pdf')
    plt.close()
    return m, picchi_ADC[-1], errori_ADC[-1]
    #plt.show()
    
    
    
gain = list()
means = list()
sigmas =list()


range_picchi = np.array([(260, 310), (480, 560), (600, 680), (720, 800), (820, 920), (1350, 1500)])
initpars     = np.array([[5000, 300, 50], [1000, 550, 50], [1000, 650, 50],  [1000, 750, 50], [2000, 900, 50], [1000, 1400, 50]])
g, m, s = apri_file_e_fit_picchi('./datadir/run020752.dat', 1, range_picchi, initpars)
gain.append(g)
means.append(m)
sigmas.append(s)

range_picchi = np.array([(470, 600), (870, 1050), (1100, 1200), (1300, 1480), (1500, 1700), (2400, 2800)])
initpars     = np.array([[3000, 600, 50], [1000, 950, 50], [1000, 1150, 50],  [1000, 1400, 100], [2000, 1600, 100], [1000, 2600, 100]])
g, m, s = apri_file_e_fit_picchi('./datadir/run020753.dat', 2, range_picchi, initpars)
gain.append(g)
means.append(m)
sigmas.append(s)

range_picchi = np.array([[ 940, 1200],[1740, 2000],[2050, 2300],[2500, 2700],[2850, 3200],[4400, 5000]])
initpars = np.array([[3000, 1100, 50], [1000, 1800, 50], [1000, 2150, 50],  [1000, 2600, 100], [2000, 3200, 100], [1000, 4800, 100]])
g, m, s = apri_file_e_fit_picchi('./datadir/run020754.dat', 5, range_picchi, initpars)
gain.append(g)
means.append(m)
sigmas.append(s)

range_picchi = np.array([[ 1300, 1700],[2500, 3000],[3150, 3700],[3900, 4250],[4350, 5000],[7300, 8100]])
initpars = np.array([[3000, 1700, 150], [1000, 2800, 150], [1000, 3500, 150],  [1000, 4000, 150], [2000, 4700, 150], [1000, 7800, 150]])
g, m, s = apri_file_e_fit_picchi('./datadir/run020755.dat', 5, range_picchi, initpars)
gain.append(g)
means.append(m)
sigmas.append(s)

range_picchi = np.array([[ 2450, 2800],[4500, 5000],[5200, 6100],[6300, 7200],[7350, 8400],[11800, 13000]])
initpars = np.array([[4000, 2600, 150], [1000, 4700, 150], [1000, 5500, 150],  [1000, 7000, 150], [2000, 7500, 150], [1000, 12000, 150]])
g, m, s = apri_file_e_fit_picchi('./datadir/run020756.dat', 10, range_picchi, initpars)
gain.append(g)
means.append(m)
sigmas.append(s)
In [3]:
def gainfunction(v, C, gamma):
    return C*v**gamma

v = [600, 650, 700, 750, 800]
In [4]:
plt.plot(v,gain,'o')
optpars, cov = opt.curve_fit(gainfunction,v,gain,p0=[0.01,4])
x=np.linspace(500,850,100)
plt.plot(x,gainfunction(x,*optpars),'--r')
#plt.show()
optpars
Out[4]:
array([3.17193497e-21, 7.49618149e+00])
In [ ]: