Scan in pressione

import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit
import os

# Questa lista contiene le varie matrici di dati, ogni matrice è un file
listaVect = []
# Questa lista invece contiene 
listaNum = []

for e in os.scandir(r"./datiScan"):
    print(e.name)
    
    if e.name.split(".")[1] != "txt": continue
    
    _, _, _, ql, _ = np.loadtxt(e, unpack = True)
    listaVect.append(ql.copy())
    listaNum.append(int(e.name.split("_")[1][:3]))
    
listaNum

.ipynb_checkpoints
241Am_000mbar.txt
241Am_000mbar_events.ade
241Am_050mbar.txt
241Am_100mbar.txt
241Am_150mbar.txt
241Am_200mbar.txt
241Am_250mbar.txt
241Am_300mbar.txt
241Am_350mbar.txt
241Am_400mbar.txt
241Am_450mbar.txt
241Am_500mbar.txt
241Am_550mbar.txt
241Am_600mbar.txt
241Am_650mbar.txt
241Am_700mbar.txt
241Am_750mbar.txt

[0, 50, 100, 150, 200, 250, 300, 350, 400, 450, 500, 550, 600, 650, 700, 750]

def gausC(x, a, mu, sigma, c):
    return a/(sigma*np.sqrt(2*np.pi)) * np.exp( -(x-mu)**2 / (2*sigma**2)) + c

myRange = ((1200, 1600), (1100, 1550), (1000, 1400), (900, 1400),
           (800, 1200), (700, 1200), (600, 1100), (500, 900),
           (300, 800), (250, 700), (240, 500), (170, 350))

lstPopt = []
lstPcov = []
lstPres = [] # Pressioni dei dati fittati

fig, ax = plt.subplots(8,2)
fig.set_size_inches(35, 50)
fig.subplots_adjust(wspace = .1)
ax = ax.flatten()

for i in range(len(listaVect)):
    h, bins = np.histogram(listaVect[i], bins = 200, range = (0, 1700))
    binc = bins[:-1] + (bins[1] - bins[0])/2
    
    ax[i].plot(binc, h, ds = "steps-mid", c = "darkgreen", lw = 2, label = "Spettro")
    ax[i].fill_between(binc, h, step = "mid", color = "lime", alpha = 1)
    
    # fittto fino a 400 mbar
    if listaNum[i] <= 550:
        cond = (binc > myRange[i][0]) & (binc < myRange[i][1])
        xData = binc[cond]
        yData = h[cond]
        popt, pcov = curve_fit(gausC, xData, yData, p0 = (100*np.max(yData), xData[np.argmax(yData)], 80, 10), 
                               absolute_sigma = True, sigma = np.sqrt(yData), maxfev = 1000000)
        ax[i].plot(xData, gausC(xData, *popt), "--k", lw = 3, label = "Gaussiana fittata")
        
        #print(popt)
        lstPopt.append(popt)
        lstPcov.append(pcov)
        lstPres.append(listaNum[i])
        
    
    ax[i].axvline(x = 100, ls = "--", c = "tab:red")
    ax[i].grid()
    ax[i].set_xlabel("Energia [ADC]", fontsize = 14)
    ax[i].set_ylabel("Entries", fontsize = 14)
    ax[i].set_title(f"Pressione: {listaNum[i]} mbar", fontsize = 16)
    ax[i].legend()

plt.show()

I vettori che uso

• listaVect: vettori dei dati
• listaNum: tutte le pressioni
• lstPopt: parametri fittati
• lstPcov: covarianze fit
• lstPres: pressioni solo dei dataset fittati
• lstMu: valori medi fittati

lstMu = [i[1] for i in lstPopt]
lstErrMu = [np.diag(i)[1]**.5 for i in lstPcov]

fig, ax = plt.subplots(figsize = (12,5))

ax.plot(lstPres, lstMu, "*g", label = "Dati fit")
#ax.errorbar(lstPres, lstMu, yerr = lstErrMu, label = "err", ls = "", marker="s")

ax.grid()
ax.set_xlabel("Pressione [mbar]", fontsize = 14)
ax.set_ylabel("$\mu$ fittato [ADC]", fontsize = 14)
ax.set_title("Mu vs Pressione", fontsize = 16)
ax.legend()


plt.show()


# Derivata
fig, ax = plt.subplots(figsize = (12,5))

ax.plot(lstPres[:-1], np.abs(np.diff(lstMu)), "*r", label = "Derivata")

ax.grid()
ax.set_xlabel("Pressione [mbar]", fontsize = 14)
#ax.set_ylabel("$\mu$ fittato", fontsize = 14)
ax.set_title("Derivata", fontsize = 16)
ax.legend()


plt.show()

lstValMedio = [i[i>100].mean() for i in listaVect] # Salvo il valor medio dei dati >50


fig, ax = plt.subplots(figsize = (12,5))

ax.plot(listaNum, lstValMedio, "*g", label = "Dati medi")

ax.grid()
ax.set_xlabel("Pressione [mbar]", fontsize = 14)
ax.set_ylabel("Val medio (>50) [ADC]", fontsize = 14)
ax.set_title("Val medio (>50) vs Pressione", fontsize = 16)
ax.legend()


plt.show()


# Derivata
fig, ax = plt.subplots(figsize = (12,5))

ax.plot(listaNum[:-1], np.abs(np.diff(lstValMedio)), "*r", label = "Derivata")

ax.grid()
ax.set_xlabel("Pressione [mbar]", fontsize = 14)
#ax.set_ylabel("Val medio (>50)", fontsize = 14)
ax.set_title("Derivata", fontsize = 16)
ax.legend()


plt.show()

lstRefP0 = [(lstMu[0] - i) for i in lstMu]
    
    
fig, ax = plt.subplots(figsize = (12,5))

ax.plot(lstPres, lstRefP0, "*g", label = "Dati")

ax.grid()
ax.set_xlabel("Pressione [mbar]", fontsize = 14)
ax.set_ylabel("fit0 - fit [ADC]", fontsize = 14)
ax.set_title("(fit0 - fit) vs Pressione", fontsize = 16)
ax.legend()


plt.show()

RIPROVIAMOCI - scan fine

import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit
import os

# Questa lista contiene le varie matrici di dati, ogni matrice è un file
listaVect = []
# Questa lista invece contiene 
listaNum = []

for e in os.scandir(r"./scanPres"):
    print(e.name)
    
    if e.name.split(".")[1] != "txt": continue
    
    _, _, _, ql, _ = np.loadtxt(e, unpack = True)
    listaVect.append(ql.copy())
    listaNum.append(int(e.name.split("_")[2][:3]))
    
listaNum = np.array(listaNum)
listaNum

241Am_scan_000mbar.txt
241Am_scan_050mbar.txt
241Am_scan_100mbar.txt
241Am_scan_150mbar.txt
241Am_scan_200mbar.txt
241Am_scan_250mbar.txt
241Am_scan_300mbar.txt
241Am_scan_350mbar.txt
241Am_scan_360mbar.txt
241Am_scan_370mbar.txt
241Am_scan_380mbar.txt
241Am_scan_390mbar.txt
241Am_scan_400mbar.txt
241Am_scan_410mbar.txt
241Am_scan_420mbar.txt
241Am_scan_430mbar.txt
241Am_scan_440mbar.txt
241Am_scan_450mbar.txt
241Am_scan_460mbar.txt
241Am_scan_470mbar.txt
241Am_scan_480mbar.txt
241Am_scan_490mbar.txt
241Am_scan_500mbar.txt
241Am_scan_510mbar.txt
241Am_scan_520mbar.txt
241Am_scan_530mbar.txt
241Am_scan_540mbar.txt
241Am_scan_550mbar.txt

array([  0,  50, 100, 150, 200, 250, 300, 350, 360, 370, 380, 390, 400,
       410, 420, 430, 440, 450, 460, 470, 480, 490, 500, 510, 520, 530,
       540, 550])

%matplotlib qt
%matplotlib inline

fig, ax = plt.subplots(14,2)
fig.set_size_inches(35, 80)
fig.subplots_adjust(wspace = .1)
ax = ax.flatten()

for i in range(len(listaVect)):
    h, bins = np.histogram(listaVect[i], bins = 200)#, range = (0, 1700))
    binc = bins[:-1] + (bins[1] - bins[0])/2
    
    ax[i].plot(binc, h, ds = "steps-mid", c = "darkgreen", lw = 2, label = "Spettro")
    ax[i].fill_between(binc, h, step = "mid", color = "lime", alpha = 1)
    
    
    
    ax[i].axvline(x = 100, ls = "--", c = "tab:red")
    ax[i].grid()
    ax[i].set_xlabel("Energia [ADC]", fontsize = 14)
    ax[i].set_ylabel("Entries", fontsize = 14)
    ax[i].set_title(f"Pressione: {listaNum[i]} mbar", fontsize = 16)
    ax[i].legend()

plt.show()

%matplotlib inline

def gausC(x, a, mu, sigma, c):
    return a/(sigma*np.sqrt(2*np.pi)) * np.exp( -(x-mu)**2 / (2*sigma**2)) + c

myRange = ((1200, 1600), (1100, 1550), (1000, 1400), (900, 1400), #000-150
           (800, 1200), (700, 1200), (600, 1100), (500, 900),     #200-350
           (450, 900), (450, 900), (450, 800), (400, 800),        #360-390
           
           (300, 800), (300, 800), (250, 800), (250, 800),        #400-430
           (200, 700), (200, 700), (200, 600), (200, 600),        #440-470
           (200, 600), (200, 600), (150, 600), (150, 600),        #480-510
           (200, 700), (200, 700), (200, 600), (200, 600),        #420-5500
          )        

lstPopt = []
lstPcov = []
lstPres = [] # Pressioni dei dati fittati

fig, ax = plt.subplots(14,2)
fig.set_size_inches(35, 80)
fig.subplots_adjust(wspace = .1)
ax = ax.flatten()

for i in range(len(listaVect)):
    h, bins = np.histogram(listaVect[i], bins = 200)#, range = (0, 1700))
    binc = bins[:-1] + (bins[1] - bins[0])/2
    
    ax[i].plot(binc, h, ds = "steps-mid", c = "darkgreen", lw = 2, label = "Spettro")
    ax[i].fill_between(binc, h, step = "mid", color = "lime", alpha = 1)
    
    # fittto fino a 400 mbar
    if listaNum[i] <= 510 :
        cond = (binc > myRange[i][0]) & (binc < myRange[i][1])
        xData = binc[cond]
        yData = h[cond]
        popt, pcov = curve_fit(gausC, xData, yData, p0 = (100*np.max(yData), xData[np.argmax(yData)], 80, 10), 
                               absolute_sigma = True, sigma = np.sqrt(yData), maxfev = 1000000)
        ax[i].plot(xData, gausC(xData, *popt), "--k", lw = 3, label = "Gaussiana fittata")
        
        #print(popt)
        lstPopt.append(popt)
        lstPcov.append(pcov)
        lstPres.append(listaNum[i])
        
    
    ax[i].axvline(x = 100, ls = "--", c = "tab:red")
    ax[i].grid()
    ax[i].set_xlabel("Energia [ADC]", fontsize = 14)
    ax[i].set_ylabel("Entries", fontsize = 14)
    ax[i].set_title(f"Pressione: {listaNum[i]} mbar", fontsize = 16)
    ax[i].legend()

plt.show()

lstPres = np.array(lstPres)

I vettori che uso

• listaVect: vettori dei dati
• listaNum: tutte le pressioni
• lstPopt: parametri fittati
• lstPcov: covarianze fit
• lstPres: pressioni solo dei dataset fittati
• lstMu: valori medi fittati

lstMu = np.array([i[1] for i in lstPopt])
lstErrMu = [np.diag(i)[1]**.5 for i in lstPcov]

fig, ax = plt.subplots(figsize = (12,5))

ax.plot(lstPres, lstMu, "*g", label = "Dati fit")
#ax.errorbar(lstPres, lstMu, yerr = lstErrMu, label = "err", ls = "", marker="s")

ax.grid()
ax.set_xlabel("Pressione [mbar]", fontsize = 14)
ax.set_ylabel("$\mu$ fittato [ADC]", fontsize = 14)
ax.set_title("Mu vs Pressione", fontsize = 16)
ax.legend()


plt.show()



tmpDer = (lstMu[1:] - lstMu[:-1]) / (lstPres[1:] - lstPres[:-1])


# Derivata
fig, ax = plt.subplots(figsize = (12,5))

#ax.plot(lstPres[:-1], np.abs(np.diff(lstMu)), "*r", label = "Derivata")
ax.plot(lstPres[:-1], np.abs(tmpDer), "*r", label = "Derivata")

ax.grid()
ax.set_xlabel("Pressione [mbar]", fontsize = 14)
#ax.set_ylabel("$\mu$ fittato", fontsize = 14)
ax.set_title("Derivata", fontsize = 16)
ax.legend()


plt.show()

print(lstMu)
print(lstPres)
print(tmpDer)

[1413.05401731 1326.5398212  1242.41820439 1143.37890442 1049.3546068
  931.39729657  824.61257177  699.95199106  666.98565919  645.51095294
  625.74161034  595.7115545   557.91133495  527.34090609  490.20213608
  459.71274306  423.29937422  395.17604674  366.69049454  331.19073455
  306.07834966  279.71235841  229.53889781  177.37628543]
[  0  50 100 150 200 250 300 350 360 370 380 390 400 410 420 430 440 450
 460 470 480 490 500 510]
[-1.73028392 -1.68243234 -1.980786   -1.88048595 -2.3591462  -2.1356945
 -2.49321161 -3.29663319 -2.14747062 -1.97693426 -3.00300558 -3.78002195
 -3.05704289 -3.713877   -3.0489393  -3.64133688 -2.81233275 -2.84855522
 -3.549976   -2.51123849 -2.63659913 -5.01734606 -5.21626124]

lstValMedio = np.array([i[i>100].mean() for i in listaVect]) # Salvo il valor medio dei dati >50


fig, ax = plt.subplots(figsize = (12,5))

ax.plot(listaNum, lstValMedio, "*g", label = "Dati medi")

ax.grid()
ax.set_xlabel("Pressione [mbar]", fontsize = 14)
ax.set_ylabel("Val medio (>100) [ADC]", fontsize = 14)
ax.set_title("Val medio (>100) vs Pressione", fontsize = 16)
ax.legend()


plt.show()



tmpDer = (lstValMedio[1:] - lstValMedio[:-1]) / (listaNum[1:] - listaNum[:-1])


# Derivata
fig, ax = plt.subplots(figsize = (12,5))

#ax.plot(listaNum[:-1], np.abs(np.diff(lstValMedio)), "*r", label = "Derivata")
ax.plot(listaNum[:-1]/1000*5, np.abs(tmpDer), "*r", label = "Derivata")

ax.grid()
ax.set_xlabel("Pressione [cm eq]", fontsize = 14)
#ax.set_ylabel("Val medio (>50)", fontsize = 14)
ax.set_title("Derivata", fontsize = 16)
ax.legend()


plt.show()

print(lstValMedio)
print(listaNum)
print(tmpDer)

[1346.71781148 1267.40401764 1185.2006111  1098.28909627  998.63617001
  893.19293611  788.78292754  664.41489156  634.25931662  609.44317157
  584.59399008  562.49749632  531.77417127  501.02528009  470.90473396
  441.68550296  414.875       387.30025424  361.78851684  340.07206693
  315.11211837  294.41438356  273.36766024  256.34938756  241.14041812
  225.67483266  214.9046337   199.81911003]
[  0  50 100 150 200 250 300 350 360 370 380 390 400 410 420 430 440 450
 460 470 480 490 500 510 520 530 540 550]
[-1.58627588 -1.64406813 -1.7382303  -1.99305853 -2.10886468 -2.08820017
 -2.48736072 -3.01555749 -2.48161451 -2.48491815 -2.20964938 -3.0723325
 -3.07488912 -3.01205461 -2.9219231  -2.6810503  -2.75747458 -2.55117374
 -2.17164499 -2.49599486 -2.06977348 -2.10467233 -1.70182727 -1.52089694
 -1.54655855 -1.0770199  -1.50855237]

lstRefP0 = [(lstMu[0] - i) for i in lstMu]
    
    
fig, ax = plt.subplots(figsize = (12,5))

ax.plot(lstPres, lstRefP0, "*g", label = "Dati")

ax.grid()
ax.set_xlabel("Pressione [mbar]", fontsize = 14)
ax.set_ylabel("fit0 - fit [ADC]", fontsize = 14)
ax.set_title("(fit0 - fit) vs Pressione", fontsize = 16)
ax.legend()


plt.show()