Spettri

In queste misure sono stati acquisiti multi photon spectra al variare della lunghezza della finestra di integrazione, con il led acceso.

Gli spettri sono quindi stati fittati con delle gaussiane, e si è plottata la delta peak-to-peak in funzione della larghezza della finestra di integrazione, per trovare il plateau

import numpy as np
from matplotlib import pyplot as plt

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

 
def gaus10Gaus(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)) + \
            args[18] * np.exp( -(x-args[19])**2 / (2*args[20]**2)) + \
            args[21] * np.exp( -(x-args[22])**2 / (2*args[23]**2)) + \
            args[24] * np.exp( -(x-args[25])**2 / (2*args[26]**2)) + \
            args[27] * np.exp( -(x-args[28])**2 / (2*args[29]**2)) 
            #args[30] * np.exp( -(x-args[31])**2 / (2*args[32]**2)) Senno col 72 non funziona
            #args[33] * np.exp( -(x-args[34])**2 / (2*args[35]**2)) 
        
def myLine(x, m, q):
    return m*x + q
    
"""
lstPopt = []
lstPcov = []
lstInt = [] # Integrale manuale a pm 3sigma
lstA = [] # Aree fittate

lstMu0 = [] # mu picco a zero
lstMu = [] # valor medio istogramma
"""
def fittaSpettri(fileName):
    x,y = np.loadtxt(fileName, unpack=True)
    #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) > 9:
        # li prendo in ordine di altezza
        tmpy = yData[idx]
        perm = np.flip(np.argsort(tmpy))
        idx = idx[perm[:9]]
        idx = np.sort(idx)
    """
    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(80)
        
    tmp.append(100)
    tmp.append(3500)
    tmp.append(1500)
    tmp = np.array(tmp)
    #print(len(tmp))
    
    
    

    
    
    # Effettuo il fit
    popt, pcov = curve_fit(gaus10Gaus, xData, yData, sigma = np.sqrt(yData), p0 = tmp, maxfev = 100000)
    #print(popt)
    
    lstPopt.append(popt)
    lstPcov.append(pcov)
    #print(popt)

    
    # Plotto spettro, fit e massimi trovati
    fig, ax = plt.subplots()
    
    ax.set_title(f"{fileName}")
    ax.set_xlabel("Energia [ADC]", fontsize = 14)
    ax.set_ylabel("Entries", fontsize = 14)
    
    
    ax.plot(x, y, ds = "steps-mid", c = "tab:green", lw = 3)
    ax.plot(xData, gaus10Gaus(xData, *popt), "--k")

    ax.plot(xData[idx], yData[idx], "*r")

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

    plt.show()
    
    
    # Integrali a mano a $\pm 3 \sigma$
    tmpLstInt = []
    tmpLstA = []
    for i in range(9): # Considero solo le 9 gaussiane fisiche
        myMu = popt[3*i+1]
        mySigma = np.abs(popt[3*i+2])
        
        inizio = myMu - 3 * mySigma
        fine = myMu + 3 * mySigma
        
        tmpCond = (x >= inizio) & (x <= fine)
        
        # Integrale manuale
        tmpInt = np.sum(y[tmpCond])        
        tmpLstInt.append(tmpInt)
        
        # Aree fittate/largh. bin
        tmpLstA.append(popt[3*i] * np.sqrt(2*np.pi) * mySigma / (xData[1] - xData[0]))
        
    lstInt.append(tmpLstInt.copy())
    lstA.append(tmpLstA.copy())
    
    
    
    # Devo fare una retta di calibrazione
    xxData = np.array([popt[3*i+1] for i in range(9)])
    yyData =  np.arange(0,9, dtype = float)
    
    ppopt, ppcov = curve_fit(myLine, xxData, yyData)
    
    fig, ax = plt.subplots()
    ax.plot(xxData, yyData, "*g", label = "Data", zorder = 2)
    ax.plot(xxData, myLine(xxData, *ppopt), c = "tab:orange", label = "Fit line", zorder = 1)
    
    ax.set_xlabel("ADC", fontsize = 14)
    ax.set_ylabel("# pe", fontsize = 14)
    
    ax.grid()
    ax.legend()
    
    plt.show()
    
    
    
    
    # Calcolo mu0 (mu riferito al picco a zero)
    # mu0 = area0/area tot
    A0 = tmpLstInt[0]
    A = np.sum(y)
    #mu0 = -np.log(tmpLstInt[0] / np.sum(y))
    mu0 = -np.log(A0/A)
    lstMu0.append(mu0)
    # FORMULA SBAGLIATA PER MOTIVI LEGATI ALLA DISLESSIA
    #errMu0 = np.sqrt( ( np.sqrt(A0) / (A0*A )  )**2 + (np.log(A0) * np.sqrt(A) / (A**2))**2)
    errMu0 = np.sqrt( 1/A0 + 1/A)
    
    
    
    # Calcolo mu (mu riferito al valor medio, ma affetto dal crosstalk)
    xCalib = myLine(x, *ppopt)
    #mu = np.sum(myLine(x, *ppopt) * y)/np.sum(y)
    mu = np.sum(xCalib * y)/np.sum(y)
    lstMu.append(mu)
    errMu = np.sqrt( (1/ np.sum(y)) * np.sum(y * (xCalib-mu)**2)) / np.sqrt(len(y))
    
    # Calcolo la probabilità di crosstalk
    xtalk = (mu-mu0)/mu
    lstXtalk.append(xtalk)
    errXtalk = np.sqrt( (mu0/mu**2 + errMu**2)**2 + (errMu0/mu)**2 )
    lstErrXtalk.append(errXtalk)
    print(mu0, errMu0, mu, errMu)

    
    

%matplotlib inline

plt.close("all")
lstPopt = []
lstPcov = []
lstGate = []
lstInt = [] # Integrale manuale a pm 3sigma
lstA = [] # Aree fittate
lstMu0 = [] # mu picco a zero
lstMu = [] # valor medio istogramma
lstXtalk = []
lstErrXtalk= []

for e in os.scandir(".\scanGate"):
    if e.name.split(".")[1] != "txt":  continue
    #print(e.path)
    fittaSpettri(e.path)
    lstGate.append(e.name.split("_")[0][4:])
    #print(e.name.split("_")[0][4:])

lstGate = np.array(lstGate, dtype = int)
lstGate

4.613178965598897 0.019212571192631887 4.721214541164828 0.04380026065659773

4.697351345472215 0.01852983211518994 4.7286131399410705 0.04343796931924798

4.66775367228986 0.018721212589633895 4.724735815654573 0.04341691236337501

4.700146330982884 0.01925403424694673 4.737429012112496 0.043409423222724994

4.7305798347239065 0.02056231053848219 4.80304747146964 0.043952376094796235

4.7352698513464695 0.020915516667433288 4.8095417177985444 0.044010029620817674

4.745082042163096 0.020851427932925038 4.806947729296337 0.044087445714784755

4.751040186986844 0.021378021444031368 4.804955513441819 0.04402848702122323

4.719169487169115 0.021241883494597984 4.813433569844982 0.044067790301908626

4.739929051846016 0.021225719293695086 4.816974387398447 0.04408508643647014

4.713670276688899 0.0197183617775058 4.80690924749114 0.044243533690762595

4.731202495888604 0.020075923246923484 4.796239420959825 0.044122916106965025

4.729323290055958 0.02054950882894474 4.791547955630241 0.0440554815505488

4.724590260688497 0.020271822360273056 4.789857136197606 0.04406699684283584

4.703231399164306 0.020608033620194085 4.78164736801757 0.04396577242560318

4.73081063793336 0.020758978941468876 4.776206380200284 0.04401127931087194

4.719133993577479 0.020671928622944032 4.781339043523081 0.04398484113037629

4.721671798335672 0.020676076065078163 4.767353660040763 0.043910072226321595

4.7296604065644665 0.020670971415976517 4.769137103701632 0.04394822303189723

4.71346302877214 0.019745028992229522 4.778523359559374 0.043878282123794896

4.667362206244539 0.01733250648447646 4.740772073450372 0.04383460583094334

4.657530311384126 0.016664748935599043 4.750956949096718 0.043921636681136654

4.620480947225183 0.02621014056400242 4.657170750579318 0.042747830662397644

array([256,  80,  88,  96, 104, 112, 120, 128, 136, 144, 152, 160, 168,
       176, 184, 192, 200, 208, 216, 224, 240, 248,  72])

def expDesc(x, a, tau, c):
    return a * (1 - c * np.exp(-x/tau))
                

#dpp21 = np.array([(i[7]-i[4]) for i in lstPopt])
dpp21 = []
for i in lstPopt:
    dpp21.append((i[7]-i[4]))
dpp21 = np.array(dpp21)
    
errDpp21 = []
for i in lstPcov:
    err2 = np.sqrt(np.diag(i)[7])
    err1 = np.sqrt(np.diag(i)[4])
    
    errDpp21.append(np.sqrt(err1**2 + err2**2))
errDpp21 = np.array(errDpp21)
    
    
cond = (lstGate >= 72) & (lstGate <= 232)
    


xData = lstGate[cond]
yData = dpp21[cond]
yErr = errDpp21[cond]

popt21, pcov21 = curve_fit(expDesc, xData, yData, sigma = yErr, 
                       absolute_sigma = True, p0 = (660, 5, 1)) 



xFine = np.linspace(xData.min(), xData.max(), 100)

fig, ax = plt.subplots()

ax.errorbar(xData, yData, yerr = yErr, ls = "", marker = "*" , c = "tab:green", label = "Dati")
ax.plot(xFine, expDesc(xFine, *popt21), ":k", label = "Fit")


ax.set_title("Stima gate ottimale", fontsize = 16)
ax.set_xlabel("Gate [ns]", fontsize = 14)
ax.set_ylabel("$\Delta pp_{21}$ [ADC]", fontsize = 14)

ax.grid(True)
ax.legend()

plt.show()


print(popt21)

chi2 = np.sqrt(np.sum(((yData - expDesc(xData, *popt21))/yErr ) **2)) / (len(xData) - 3)
print(f"il chi2 ridotto vale {chi2:.4f}")

print (f"Il tau fittato vale {popt21[1]:.4f} pm {np.sqrt(np.diag(pcov21))[1]:.4f}")

[649.79021619  32.33439057   1.86598064]
il chi2 ridotto vale 0.6731
Il tau fittato vale 32.3344 pm 0.4417

def expDesc(x, a, tau, c):
    return a * (1 - c * np.exp(-x/tau))
                

#dpp21 = np.array([(i[7]-i[4]) for i in lstPopt])
dpp21 = []
for i in lstPopt:
    dpp21.append((i[13]-i[10]))
dpp21 = np.array(dpp21)
    
errDpp21 = []
for i in lstPcov:
    err2 = np.sqrt(np.diag(i)[13])
    err1 = np.sqrt(np.diag(i)[10])
    
    errDpp21.append(np.sqrt(err1**2 + err2**2))
errDpp21 = np.array(errDpp21)
    
    
cond = (lstGate >= 72) & (lstGate <= 232)
    


xData = lstGate[cond]
yData = dpp21[cond]
yErr = errDpp21[cond]

popt43, pcov43 = curve_fit(expDesc, xData, yData, sigma = yErr, 
                       absolute_sigma = True, p0 = (660, 5, 1)) 



xFine = np.linspace(xData.min(), xData.max(), 100)

fig, ax = plt.subplots()

ax.errorbar(xData, yData, yerr = yErr, ls = "", marker = "*" , c = "tab:green", label = "Dati")
ax.plot(xFine, expDesc(xFine, *popt43), ":k", label = "Fit")


ax.set_title("Stima gate ottimale", fontsize = 16)
ax.set_xlabel("Gate [ns]", fontsize = 14)
ax.set_ylabel("$\Delta pp_{43}$ [ADC]", fontsize = 14)

ax.grid(True)
ax.legend()

plt.show()


print(popt43)

chi2 = np.sqrt(np.sum(((yData - expDesc(xData, *popt43))/yErr ) **2)) / (len(xData) - 3)
print(f"il chi2 ridotto vale {chi2:.4f}")

print (f"Il tau fittato vale {popt43[1]:.4f} pm {np.sqrt(np.diag(pcov43))[1]:.4f}")

[647.15198087  31.3778212    2.04922845]
il chi2 ridotto vale 0.6262
Il tau fittato vale 31.3778 pm 0.4092

Copio da file alcuni calcoli

# Inizio - picco - 1/e picco
times = np.array(((-9, 6.8, 43.6), 
            (-8.8, 7,  40.60), 
            (-9.6, 6.2, 39), 
            (-9.2, 6.8, 39.4), 
            (-8.2, 7, 43), 
            (-8.6, 7, 39.8)) )

timeToPeak = np.mean(times[:,1] - times[:,0])
errToPeak = np.std(times[:,1] - times[:,0])

timeToE = np.mean(times[:,2] - times[:,1])
errToE = np.std(times[:,2] - times[:,1])
      
print(f"Il tempo dall'inizo al picco vale {timeToPeak:.4f} ± {errToPeak:.4f}\nIl tempo dal picco a 1/e, ovvero tau, vale {timeToE:.4f} ± {errToE:.4f}")

Il tempo dall'inizo al picco vale 15.7000 ± 0.2517
Il tempo dal picco a 1/e, ovvero tau, vale 34.1000 ± 1.6723

Compatibilità tra le tau

def myZ(x1, x2, err1, err2):
    return np.abs(x1-x2)/np.sqrt(err1**2 + err2**2)

# Tau 21 vs tau 43

Z = myZ(popt21[1], popt43[1], np.sqrt(np.diag(pcov21)[1]), np.sqrt(np.diag(pcov43)[1]) )
print(f"Dpp21 vs Dpp43:\t {Z}")

Z = myZ(popt21[1], timeToE, np.sqrt(np.diag(pcov21)[1]), errToE )
print(f"Dpp21 vs osci:\t {Z}")

Z = myZ(popt43[1], timeToE, np.sqrt(np.diag(pcov43)[1]), errToE )
print(f"Dpp43 vs osci:\t {Z}")

Dpp21 vs Dpp43:	 1.5885840248464287
Dpp21 vs osci:	 1.0207739407451217
Dpp43 vs osci:	 1.581129661924065

np.sqrt(np.diag(pcov21)[1])

0.4417200444679869

# print(lstInt) # print(lstA) for idx, (i, j) in enumerate(zip(lstInt, lstA)): print("\n") for k,l in zip(i, j): print(f"Gaus{idx} - Integrale a +- 3sigma: {k:.2f} - Area fittata/largh. bin: {l:.2f}")

whos

Variable       Type           Data/Info
---------------------------------------
Z              float64        1.581129661924065
ax             AxesSubplot    AxesSubplot(0.125,0.125;0.775x0.755)
chi2           float64        0.6262087066682316
cond           ndarray        23: 23 elems, type `bool`, 23 bytes
curve_fit      function       <function curve_fit at 0x000001C6D68E1940>
dpp21          ndarray        23: 23 elems, type `float64`, 184 bytes
e              DirEntry       <DirEntry 'Gate72_histo.txt'>
err1           float64        1.0858375592553784
err2           float64        1.3739966588059904
errDpp21       ndarray        23: 23 elems, type `float64`, 184 bytes
errToE         float64        1.6723237326147924
errToPeak      float64        0.2516611478423587
expDesc        function       <function expDesc at 0x000001C6D9532D30>
fig            Figure         Figure(432x288)
find_peaks     function       <function find_peaks at 0x000001C6D720E310>
fittaSpettri   function       <function fittaSpettri at 0x000001C6D7210310>
gaus10Gaus     function       <function gaus10Gaus at 0x000001C6D3FD20D0>
i              ndarray        30x30: 900 elems, type `float64`, 7200 bytes
lstA           list           n=23
lstErrXtalk    list           n=23
lstGate        ndarray        23: 23 elems, type `int32`, 92 bytes
lstInt         list           n=23
lstMu          list           n=23
lstMu0         list           n=23
lstPcov        list           n=23
lstPopt        list           n=23
lstXtalk       list           n=23
myLine         function       <function myLine at 0x000001C6D7210280>
myZ            function       <function myZ at 0x000001C6D9356EE0>
np             module         <module 'numpy' from 'C:\<...>ges\\numpy\\__init__.py'>
os             module         <module 'os' from 'C:\\Us<...>\\anaconda3\\lib\\os.py'>
pcov21         ndarray        3x3: 9 elems, type `float64`, 72 bytes
pcov43         ndarray        3x3: 9 elems, type `float64`, 72 bytes
plt            module         <module 'matplotlib.pyplo<...>\\matplotlib\\pyplot.py'>
popt21         ndarray        3: 3 elems, type `float64`, 24 bytes
popt43         ndarray        3: 3 elems, type `float64`, 24 bytes
timeToE        float64        34.1
timeToPeak     float64        15.700000000000001
times          ndarray        6x3: 18 elems, type `float64`, 144 bytes
xData          ndarray        20: 20 elems, type `int32`, 80 bytes
xFine          ndarray        100: 100 elems, type `float64`, 800 bytes
yData          ndarray        20: 20 elems, type `float64`, 160 bytes
yErr           ndarray        20: 20 elems, type `float64`, 160 bytes

for i,j in zip(lstXtalk, lstErrXtalk):
    print(f"{i:.4f} ± {j:.4f}")

0.0229 ± 0.2089
0.0066 ± 0.2120
0.0121 ± 0.2110
0.0079 ± 0.2113
0.0151 ± 0.2070
0.0154 ± 0.2067
0.0129 ± 0.2073
0.0112 ± 0.2078
0.0196 ± 0.2057
0.0160 ± 0.2063
0.0194 ± 0.2060
0.0136 ± 0.2077
0.0130 ± 0.2080
0.0136 ± 0.2079
0.0164 ± 0.2077
0.0095 ± 0.2094
0.0130 ± 0.2084
0.0096 ± 0.2097
0.0083 ± 0.2099
0.0136 ± 0.2084
0.0155 ± 0.2096
0.0197 ± 0.2083
0.0079 ± 0.2149

perm = np.argsort(lstGate)

_lstGate = lstGate[perm]
_lstXtalk = np.array(lstXtalk)[perm]
_lstErrXtalk = np.array(lstErrXtalk)[perm]

for i,j in zip(_lstXtalk, _lstErrXtalk):
    print(f"{i:.4f} ± {j:.4f}")
    
plt.errorbar(_lstGate, _lstXtalk)#, yerr = _lstErrXtalk)
plt.show()

0.0079 ± 0.2149
0.0066 ± 0.2120
0.0121 ± 0.2110
0.0079 ± 0.2113
0.0151 ± 0.2070
0.0154 ± 0.2067
0.0129 ± 0.2073
0.0112 ± 0.2078
0.0196 ± 0.2057
0.0160 ± 0.2063
0.0194 ± 0.2060
0.0136 ± 0.2077
0.0130 ± 0.2080
0.0136 ± 0.2079
0.0164 ± 0.2077
0.0095 ± 0.2094
0.0130 ± 0.2084
0.0096 ± 0.2097
0.0083 ± 0.2099
0.0136 ± 0.2084
0.0155 ± 0.2096
0.0197 ± 0.2083
0.0229 ± 0.2089

Faccio la media dei cross talk. L'errore è la media / N

meanXtalk = np.mean(_lstXtalk)
errxTalk = np.sqrt(np.sum(_lstErrXtalk**2))/len(_lstErrXtalk)

print(f"Il cross talk medio vale {100*meanXtalk:.2f} pm {100*errxTalk:.2f} %")

Il cross talk medio vale 1.36 pm 4.35 %