Fitto gli spettri illuminati, quindi plotto la delta pp al variare del bias.

I files sono in punto3_bis

Usiamo 3 dataset per ogni bias. Il fatto che si vedano le linee tratteggiate che fanno cose strane dipende dal fatto che non ho corretto a mano i files, ma ho usato pandas che non si lamenta. Semplicemente non mettiamo quei plot nella presentazione...

import numpy as np
from matplotlib import pyplot as plt

from scipy.optimize import curve_fit
from scipy.signal import find_peaks
import os

import pandas as pd

def gaus5Gaus(x,*args):
    return  args[0] * np.exp( -(x-args[1])**2 / (2*args[2]**2)) + \
            args[3] * np.exp( -(x-args[4])**2 / (2*args[5]**2)) + \
            args[6] * np.exp( -(x-args[7])**2 / (2*args[8]**2)) + \
            args[9] * np.exp( -(x-args[10])**2 / (2*args[11]**2)) + \
            args[12] * np.exp( -(x-args[13])**2 / (2*args[14]**2)) + \
            args[15] * np.exp( -(x-args[16])**2 / (2*args[17]**2))


def myLine(x, m, q):
    return m*x + q


def fittaSpettri(fileName):
    #x,y = np.loadtxt(fileName, unpack=True)
    df = pd.read_csv(fileName, sep = "\t", header = None, names = ("x", "y"))
    x = df.x.to_numpy()
    y = df.y.to_numpy()
    #print(fileName)
    
    
    cond = (x > -500) & (x < 6921) & (y > 0)
    cond = (x > -500) & (x < 5500) & (y > 0)
    xData = x[cond]
    yData = y[cond]
    
    # Trovo i massimi
    idx, kkk = find_peaks(yData, width = 3, prominence = 50, distance = 30)
    """
    if len(idx) > 5:
        # li prendo in ordine di altezza
        tmpy = yData[idx]
        perm = np.flip(np.argsort(tmpy))
        idx = idx[perm[:5]]
        idx = np.sort(idx)
    """
    if len(idx) > 5:
        idx = idx[:5]

    
    """
    print(idx)
    print(len(idx))
    print(kkk)
    """
    
    # Definisco gli startPars
    tmp = []
    for i in idx:
        tmp.append(yData[i])
        tmp.append(xData[i])
        tmp.append(30)
        
    tmp.append(yData[idx[2]]/2)
    tmp.append(xData[idx[2]])
    tmp.append(750)
    tmp = np.array(tmp)
    #print(len(tmp))
    
    
    # Nuove condizioni
    cond = (x > -500) & (x < (xData[idx[4]] + 100)) & (y > 0)
    xxData = x[cond]
    yyData = y[cond]
    
    
    # Effettuo il fit
    popt, pcov = curve_fit(gaus5Gaus, xxData, yyData, sigma = np.sqrt(yyData), p0 = tmp, maxfev = 1000000)
    #print(popt)
    
    lstPopt.append(popt)
    lstPcov.append(pcov)
    #print(popt)
    

    
    # Plotto spettro, fit e massimi trovati
    fig, ax = plt.subplots(figsize=(8,6))
    
    ax.set_title(f"{fileName}")
    
    ax.plot(x, y, ds = "steps-mid", c = "tab:green", lw = 3, ls = "", marker = ".")
    ax.plot(xxData, gaus5Gaus(xxData, *popt), "--k")
    ax.set_xlabel('Energia [ADC]', fontsize = 14)
    ax.set_ylabel('Entries', fontsize = 14)
    ax.plot(xData[idx], yData[idx], "*r", ms = 8)

    ax.set_xlim(-500, 3000)
    ax.grid()

    plt.show()
    

%matplotlib inline

plt.close("all")
# Parametri fit
lstPopt = []
lstPcov = []
# Bias
lstBias = []


for e in os.scandir(".\punto3_bis"):
    if e.name.split(".")[-1] != "txt":  continue
    if e.name.split("_")[-1][:3] == "old":  continue
    if e.name.split("_")[-2] in ("4", "5", "0"):  continue
    
    
    """
    if e.path == r".\punto3_bis\MPH_l_55_3_histo.txt":
        %matplotlib qt
        fittaSpettri(e.path)
    """

    
    #print(e.path)
    fittaSpettri(e.path)
    #print(lstPopt,lstPcov)
    lstBias.append(e.name.split("_")[2])
    print(e.name.split("_")[0][4:])

lstBias = np.array(lstBias, dtype = float)
lstPopt = np.array(lstPopt)
lstPcov = np.array(lstPcov)

C:\Users\steca\anaconda3\lib\site-packages\scipy\optimize\minpack.py:833: OptimizeWarning: Covariance of the parameters could not be estimated
  warnings.warn('Covariance of the parameters could not be estimated',

C:\Users\steca\anaconda3\lib\site-packages\scipy\optimize\minpack.py:833: OptimizeWarning: Covariance of the parameters could not be estimated
  warnings.warn('Covariance of the parameters could not be estimated',

var = []
for i in np.arange(0,36):
    var1 = np.sqrt(np.diag(lstPcov[i]))
    a = np.divide(1,var1)
    if np.sum(a) == 0:
        var1 = np.sqrt(np.diag(lstPcov[i-1]))
    var.append(var1[1]); var.append(var1[4])  ; var.append(var1[7]) ; var.append(var1[10])
var = np.array(var) ; var = np.reshape(var,  (36, 4))
print(var, lstPcov.shape)

[[0.25648602 0.23725136 0.31307344 0.42509318]
 [0.27077933 0.25194246 0.31922942 0.39393445]
 [0.29671669 0.2803174  0.36725708 0.47254018]
 [0.27208825 0.22240223 0.25400812 0.3549257 ]
 [0.31758412 0.25392219 0.27335925 0.34521088]
 [0.30035907 0.24800875 0.27810633 0.44156789]
 [0.31503492 0.23555512 0.24067386 0.30157728]
 [0.31488688 0.23454178 0.23692316 0.29349178]
 [0.32929118 0.24165647 0.24800234 0.31765363]
 [0.32850086 0.22305249 0.21299896 0.24077332]
 [0.37459751 0.25275861 0.23583236 0.2748778 ]
 [0.38136788 0.25847864 0.24720306 0.4119007 ]
 [0.3785258  0.25578311 0.22702216 0.24451843]
 [0.41909482 0.27274963 0.24565624 0.26060392]
 [0.41909482 0.27274963 0.24565624 0.26060392]
 [0.63776735 0.3952531  0.34691153 0.31848464]
 [1.11164616 0.40638932 0.35669384 0.35792235]
 [1.11164616 0.40638932 0.35669384 0.35792235]
 [1.12144159 0.69310178 0.60366995 0.56285552]
 [0.53487239 0.452854   0.42859712 0.47701343]
 [0.53229668 0.43441214 0.43091953 0.4691444 ]
 [0.17914069 0.14925685 0.16237337 0.20617724]
 [0.19199238 0.16102151 0.1754827  0.22300351]
 [0.63844237 0.35506475 0.38996764 0.49841769]
 [0.25267414 0.19723874 0.2058513  0.24790683]
 [0.26875178 0.21119031 0.22010083 0.26581496]
 [0.24283316 0.19173891 0.19842901 0.24334749]
 [0.33616138 0.26633646 0.27424609 0.33227872]
 [0.38344349 0.27456542 0.29258719 0.35331462]
 [0.31219389 0.24480627 0.25678459 0.30875507]
 [0.39060653 0.30956742 0.3326798  0.41158033]
 [0.41776535 0.32802065 0.3526231  0.42961568]
 [0.43312014 0.33452512 0.36098847 0.43739186]
 [0.47640439 0.40295199 0.44409527 0.55646455]
 [0.46791048 0.38942112 0.43889624 0.55412888]
 [0.48250448 0.37801438 0.4387166  0.54255701]] (36, 18, 18)

lstBias = np.array([52.7 , 52.75, 52.8, 52.85, 52.9, 52.95, 53, 53.5 , 54, 54.5, 55, 55.5]) 

dpp21 = []
dpp43 = []

for i in lstPopt:
    dpp21.append((i[7]-i[4]))   
    dpp43.append((i[13]-i[10]))

errp1 = var[:,0] ; errp1 =  np.reshape(errp1,  (12,3))
errp2 = var[:,1] ; errp2 =  np.reshape(errp2,  (12,3))
errp3 = var[:,2] ; errp3 =  np.reshape(errp3,  (12,3))
errp4 = var[:,3] ; errp4 =  np.reshape(errp4,  (12,3))


errdpp21 = np.sqrt(errp1**2+errp2**2) ; errdpp43 = np.sqrt(errp3**2+errp4**2)

dpp21 = np.array(dpp21)
dpp21 =  np.reshape(dpp21,  (int(len(lstBias)), 3))

t = []

for j in np.arange(0,12):
    v0 = dpp21[j,:]
    v1 = errdpp21[j,:]
    t.append(np.average(v0,weights=1/v1))
# dpp21_mean = np.mean(dpp21, axis = 1)
# dpp21_std = np.std(dpp21, axis = 1) 
#print(dpp21)

dpp43 = np.array(dpp43)
dpp43 =  np.reshape(dpp43,  (int(len(lstBias)), 3))

s = []

for j in np.arange(0,12):
    v0 = dpp43[j,:]
    v2 = errdpp43[j,:]
    s.append(np.average(v0,weights=1/v2))
s = np.array(s)
t = np.array(t)

err_t = []
for j in np.arange(0,12):
    w = 1/np.sqrt(np.sum(1/errdpp21[j,:]**2))
    err_t.append(w)
print(err_t)
err_s = []
for j in np.arange(0,12):
    w = 1/np.sqrt(np.sum(1/errdpp43[j,:]**2))
    err_s.append(w)

[0.21564718598932855, 0.21959884512947808, 0.22975785034518828, 0.25038500194106006, 0.279613799064912, 0.5586739856756243, 0.459812886338669, 0.16622175325391225, 0.18653991304969425, 0.2477943609374534, 0.3027649618913991, 0.3551423875444803]

def myLine(x,m,q):
    return m*x+q

fig, ax = plt.subplots(1, 2,figsize=(15,8))

popt1, pcov1 = curve_fit(myLine,lstBias, t, sigma = np.sqrt(2)*np.ones(t.size), absolute_sigma = True)
popt2, pcov2 = curve_fit(myLine,lstBias, s, sigma = np.sqrt(2)*np.ones(s.size), absolute_sigma = True)

myChi2_21 = np.sqrt(np.sum(((t-myLine(lstBias, *popt1))/np.sqrt(2) )**2)) / (len(t) - len(popt1) )
print(myChi2_21)

myChi2_43 = np.sqrt(np.sum(((s-myLine(lstBias, *popt2))/np.sqrt(2) )**2)) / (len(s) - len(popt2) )
print(myChi2_43)

ax[0].errorbar(lstBias, t, yerr =  np.sqrt(2)*np.ones(t.size), marker = "o", ls = "")
ax[1].errorbar(lstBias, s, yerr =  np.sqrt(2)*np.ones(s.size), marker = "o", ls = "")

ax[0].plot(lstBias, myLine(lstBias, *popt1), "-.r", ms = 5)
ax[1].plot(lstBias, myLine(lstBias, *popt2), "-.r", ms= 5)

ax[0].set_title("$\Delta_{pp}$ 21 vs Voltage", fontsize = 16)
ax[1].set_title("$\Delta_{pp}$ 43 vs Voltage", fontsize = 16)

for i in ax:
    i.grid()
    i.set_xlabel("Bias [mV]", fontsize = 14)
    i.set_ylabel("$\Delta_{pp}$", fontsize = 14)
    
plt.show()

0.763640982398013
1.0198864082868924

r1 = -popt1[1]/popt1[0] ; r2 = -popt2[1]/popt2[0]

err_r1 = np.sqrt (pcov1[1,1] / popt1[0]**2 + popt1[1]**2/popt1[0]**4 * pcov1[0,0]  + 2*popt1[1]/(popt1[0]**3) * pcov1[1,0])
err_r2 = np.sqrt (pcov2[1,1] / popt2[0]**2 + popt2[1]**2/popt2[0]**4 * pcov2[0,0]  + 2*popt2[1]/(popt2[0]**3) * pcov2[1,0])

r_tot = ((r1*(1/err_r1**2))+(r2*(1/err_r2)**2))/((1/err_r1)**2+(1/err_r2)**2)
err_mean = 1/(np.sqrt(((1/err_r1)**2+(1/err_r2)**2)))

print(f"Il breakdown {round(r_tot,3)} V +- {round(err_mean,3)} V")

Il breakdown 52.167 V +- 0.176 V

1

1