API tools faq
paste
Advertisement
Stex6299

analisi cobalto

Stex6299
May 14th, 2021
684
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Python 21.34 KB | None | 0 0
raw download clone embed print report
  1. # -*- coding: utf-8 -*-
  2. """
  3. Created on Fri May  7 09:37:02 2021
  4.  
  5. Autore: Stefano Carsi
  6. Data ultima modifica: 202105111551
  7.  
  8.    
  9. - Cosa ho fatto:
  10.    *
  11.    
  12. - Cosa devo fare:
  13.    *
  14.    
  15. - Per relazione:
  16.    *
  17.    
  18.    > 1 = Spesso
  19.    > 2 = Sottile
  20.    
  21.    > L'ADC è a 14 bit:  va da 0 a 16384-1
  22.    
  23.    
  24.    + Perchè non vediamo il 136 in quello spesso? Saturazione positiva?
  25. """
  26.  
  27.  
  28. #INIZIO STANDARD
  29. %matplotlib inline
  30. import matplotlib.pyplot as plt
  31. import numpy as np
  32. import matplotlib as mpl
  33. from scipy.optimize import curve_fit
  34.  
  35. import copy
  36.  
  37. from matplotlib.patches import Rectangle
  38.  
  39. import os
  40. path = r"F:\Corsi\Laboratorio 4 - Fisica nucleare e subnucleare\4 - Co-57\PITONE"
  41. os.chdir(path)
  42.  
  43. Export = False
  44.  
  45.  
  46. #%% Carico i dati
  47. nev, time, ADC1, time1, ADC2, time2 = \
  48. np.loadtxt(r"..\DATI\run020006.dat", unpack=True)#, dtype=np.int32)
  49.  
  50. #dati = np.loadtxt(r"..\DATI\run020006.dat")
  51.  
  52. # Il tempo è in tick: 1tick=2ns ==> Converto il tempo in ns
  53. time1 *= 2
  54. time2 *= 2
  55.  
  56.  
  57. #%% Plotto gli istogrammi brutti
  58.  
  59. range1 = (2500,14000)
  60. range2 = (200,14000)
  61.  
  62. binnn1 = 200 # In origine era 250
  63. binnn2 = 250
  64.  
  65. # Sopra i 14000 c'è la saturazione?
  66. h1, bins1 = np.histogram(ADC1, bins = binnn1, range = range1)
  67. h2, bins2 = np.histogram(ADC2, bins = binnn2, range = range2)
  68.  
  69. binc1 = bins1[:-1] + bins1[1] - bins1[0]
  70. binc2 = bins2[:-1] + bins2[1] - bins2[0]
  71.  
  72.  
  73.  
  74. fig, ax = plt.subplots(2, 1)#, sharex=True)
  75. #fig.subplots_adjust(hspace=0)
  76. fig.set_size_inches(12,5)
  77. fig.suptitle("Istogrammi individuali", fontsize=16)
  78.  
  79.  
  80. ax[0].plot(binc1, h1, ds="steps-mid", label = "Istogramma spesso")
  81. ax[1].plot(binc2, h2, ds="steps-mid", label = "Istogramma sottile")
  82.  
  83.  
  84.  
  85.  
  86. for i in range(2):
  87.     ax[i].grid()
  88.     ax[i].legend(fontsize=14)
  89.     ax[i].set_ylabel("Entries", fontsize=14)
  90.     ax[i].set_xlabel("ADC", fontsize=14)
  91.  
  92. fig.subplots_adjust(hspace=.4)
  93.  
  94.  
  95. if Export:
  96.     plt.savefig('../RELAZIONE/FIGURE/istogrammi-brutti.eps', format='eps')
  97. plt.show()
  98.  
  99.  
  100.  
  101. #%% Plotto l'istogramma bidimensionale
  102.  
  103. # Esistono problemi con i pixels zero: li risolvo
  104. my_cmap = copy.copy(mpl.cm.viridis) # copy the default cmap
  105. my_cmap.set_bad(my_cmap.colors[0])
  106.  
  107. fig, ax = plt.subplots()
  108. fig.set_size_inches(12,5)
  109. ax.set_title("Istogramma bidimensionale", fontsize=16)
  110.  
  111. myPlot = ax.hist2d(ADC1, ADC2, bins = (binnn1, binnn2), range = (range1, range2),
  112.                     norm = mpl.colors.LogNorm(), cmap = my_cmap)
  113.  
  114. ax.set_xlabel("Rivelatore spesso (ADC)", fontsize = 14)
  115. ax.set_ylabel("Rivelatore sottile (ADC)", fontsize = 14)
  116.  
  117. fig.colorbar(myPlot[3])
  118.  
  119.  
  120.  
  121. if Export:
  122.     plt.savefig('../RELAZIONE/FIGURE/istogramma-2d.eps', format='eps')
  123. plt.show()
  124.  
  125.  
  126.  
  127. #%% Proviamo il contour
  128.  
  129. xcentr = myPlot[1][:-1] + myPlot[1][1] - myPlot[1][0]
  130. ycentr = myPlot[2][:-1] + myPlot[2][1] - myPlot[2][0]
  131.  
  132.  
  133. fig, ax = plt.subplots()
  134. fig.set_size_inches(12,5)
  135. ax.set_title("Curve di livello (Contour)", fontsize=16)
  136.  
  137. ax.set_xlabel("Rivelatore spesso (ADC)", fontsize = 14)
  138. ax.set_ylabel("Rivelatore sottile (ADC)", fontsize = 14)
  139.  
  140.  
  141. myContour = ax.contour(*np.meshgrid(xcentr, ycentr, indexing = "ij"),
  142.                        myPlot[0], levels = range(20,200,20))#np.logspace(1,2.3,5))
  143.  
  144.  
  145. rect1 = Rectangle((2840,10850), width = 900, height = 2400,
  146.                   ec='tab:green', fc='none', lw = 3)
  147. rect2 = Rectangle((8663,3470), width = 1200, height = 1400,
  148.                   ec='tab:green', fc='none', lw = 3)
  149.  
  150. rect3 = Rectangle((10500,490), width = 3000, height = 800,
  151.                   ec='tab:green', fc='none', lw = 3)
  152. rect4 = Rectangle((10500,1600), width = 3000, height = 1050,
  153.                   ec='tab:red', fc='none', lw = 3)
  154.  
  155. rectList = (rect1, rect2, rect3, rect4)
  156.  
  157.  
  158. ax.add_patch(rect1)
  159. ax.add_patch(rect2)
  160. ax.add_patch(rect3)
  161. ax.add_patch(rect4)
  162.  
  163. fig.colorbar(myContour)
  164.  
  165.  
  166.  
  167. # Stampo i rettangoloni
  168. plt.annotate("1", xytext=(5500, 9500), xy=(4000, 12000), xycoords='data',
  169.              arrowprops=dict(arrowstyle="fancy",
  170.                             connectionstyle="arc3, rad=.3", color="tab:green"),
  171.              c="tab:green", size=50)
  172.  
  173. plt.annotate("2", xytext=(10500, 7000), xy=(9000, 5000), xycoords='data',
  174.              arrowprops=dict(arrowstyle="fancy",
  175.                             connectionstyle="arc3, rad=.6", color="tab:green"),
  176.              c="tab:green", size=50)
  177.  
  178. plt.annotate("3", xytext=(6000, 4000), xy=(10300, 1000), xycoords='data',
  179.              arrowprops=dict(arrowstyle="fancy",
  180.                             connectionstyle="arc3, rad=.3", color="tab:green"),
  181.              c="tab:green", size=50)
  182.  
  183. plt.annotate("4", xytext=(13000, 4000), xy=(12000, 3000), xycoords='data',
  184.              arrowprops=dict(arrowstyle="fancy",
  185.                             connectionstyle="arc3, rad=.6", color="tab:red"),
  186.              c="tab:red", size=50)
  187.  
  188.  
  189. if Export:
  190.     plt.savefig('../RELAZIONE/FIGURE/contour.eps', format='eps')
  191. plt.show()
  192.  
  193.  
  194.  
  195.  
  196.  
  197. #%%
  198. """
  199. in ogni rettangolo devo proiettare sui due assi e fittare con una gaussiana
  200.  
  201. """
  202.  
  203.  
  204. #%% Fit gaussiani
  205.  
  206. def gauss(x,a,mu,sigma):
  207.     return a/(sigma*np.sqrt(2*np.pi)) * np.exp(-(x-mu)**2/(2*(sigma**2)))
  208.  
  209.  
  210. lista_popt = []
  211. lista_pcov = []
  212. lista_chi2 = []
  213.  
  214. def calibra(idxRect, startingPts):
  215.    
  216.     myRect = rectList[idxRect]
  217.    
  218.     startx = myRect.get_x()
  219.     endx = startx + myRect.get_width()
  220.    
  221.     idx_startx = np.where(xcentr >= startx)[0][0]
  222.     idx_endx = np.where(xcentr >= endx)[0][0]
  223.    
  224.    
  225.     starty = myRect.get_y()
  226.     endy = starty + myRect.get_height()
  227.    
  228.     idx_starty = np.where(ycentr >= starty)[0][0]
  229.     idx_endy = np.where(ycentr >= endy)[0][0]
  230.    
  231.    
  232.    
  233.     subdata = myPlot[0][ idx_startx:idx_endx , idx_starty:idx_endy]
  234.     subdata_sqrt = np.sqrt(subdata)
  235.    
  236.    
  237.    
  238.    
  239.     # *** EFFETTUO LE PROIEZIONI ***
  240.     # Alla terza riga proietto per tutta la colonna/riga, di cui il
  241.     # rettangolo è intersezione
  242.    
  243.     # Quello che vede il rivelatore sottile, proietto a sx
  244.     proiez_sottile = np.sum(subdata, axis = 0)
  245.     proiez_sottile_sqrt = np.sqrt(np.sum(subdata_sqrt**2, axis = 0))
  246.    
  247.     proiez_sottile_completo = np.sum( myPlot[0][idx_startx:idx_endx , :], axis = 0)
  248.    
  249.     # Quello che vede il rivelatore spesso, proietto sotto
  250.     proiez_spesso = np.sum(subdata, axis = 1)
  251.     proiez_spesso_sqrt = np.sqrt(np.sum(subdata_sqrt**2, axis = 1))
  252.    
  253.     proiez_spesso_completo = np.sum( myPlot[0][:, idx_starty:idx_endy], axis = 1)
  254.  
  255.    
  256.    
  257.    
  258.     # Ottengo le x delle proiezioni
  259.     vect_spesso = np.array([i for i in xcentr if (i >= startx) & (i<=endx)])
  260.     vect_sottile = np.array([i for i in ycentr if (i >= starty) & (i<=endy)])
  261.  
  262.    
  263.  
  264.    
  265.     # fitto con gaussiane
  266.     popt_spesso, pcov_spesso = curve_fit(gauss, vect_spesso, proiez_spesso,
  267.                                          sigma = proiez_spesso_sqrt,
  268.                                          p0 = startingPts[0])
  269.    
  270.     popt_sottile, pcov_sottile = curve_fit(gauss, vect_sottile, proiez_sottile,
  271.                                          sigma = proiez_sottile_sqrt,
  272.                                          p0 = startingPts[1])
  273.    
  274.    
  275.     lista_popt.append([popt_spesso, popt_sottile])
  276.     lista_pcov.append([pcov_spesso, pcov_sottile])
  277.    
  278.    
  279.    
  280.    
  281.     # *** LO PLOTTO ***
  282.    
  283.    
  284.     fig, ax = plt.subplots(2, 1)#, sharex=True)
  285.     #fig.subplots_adjust(hspace=0)
  286.     fig.set_size_inches(12,5)
  287.     fig.subplots_adjust(hspace=.4)
  288.  
  289.    
  290.     fig.suptitle(f"Regione {idxRect+1}", fontsize=16)
  291.    
  292.     # Plotto le proiezioni (tanto dopo si vede già tutto)
  293.     #ax[0].plot(vect_spesso, proiez_spesso, label = "Proiezione rivelatore spesso")
  294.     #ax[1].plot(vect_sottile, proiez_sottile, label = "Proiezione rivelatore sottile")
  295.    
  296.    
  297.     # Plotto la somma su tutta la striscia
  298.     ax[0].plot(xcentr, proiez_spesso_completo,  
  299.                label = "Proiezione rivelatore spesso")
  300.     ax[1].plot(ycentr, proiez_sottile_completo,  
  301.                label = "Proiezione rivelatore sottile")
  302.    
  303.    
  304.     # Plotto le curve fittate
  305.     ax[0].plot(vect_spesso, gauss(vect_spesso, *popt_spesso), ":k",
  306.                label = "Curva fittata")
  307.     ax[1].plot(vect_sottile, gauss(vect_sottile, *popt_sottile), ":k",
  308.                label = "Curva fittata")
  309.    
  310.     # delimito le regioni di fit
  311.     ax[0].axvline(x=startx, ls="--", c="grey")
  312.     ax[0].axvline(x=endx, ls="--", c="grey")
  313.     ax[1].axvline(x=starty, ls="--", c="grey")
  314.     ax[1].axvline(x=endy, ls="--", c="grey")
  315.    
  316.    
  317.     for i in range(2):
  318.         ax[i].grid()
  319.         ax[i].legend(fontsize=14)
  320.         ax[i].set_ylabel("Entries", fontsize=14)
  321.         ax[i].set_xlabel("ADC", fontsize=14)
  322.    
  323.  
  324.     if Export:
  325.         plt.savefig(f'../RELAZIONE/FIGURE/fit-gauss-{idxRect}.eps', format='eps')
  326.     plt.show()
  327.    
  328.     mychi2_spesso = np.sum(((proiez_spesso-gauss(vect_spesso, *popt_spesso))/proiez_spesso_sqrt)**2) / (len(proiez_spesso)-3)
  329.     mychi2_sottile = np.sum(((proiez_sottile-gauss(vect_sottile, *popt_sottile))/proiez_sottile_sqrt)**2) / (len(proiez_sottile)-3)
  330.    
  331.     lista_chi2.append([mychi2_spesso, mychi2_sottile])
  332.  
  333.  
  334.  
  335.  
  336. calibra(0, ((100, 3300, 400), (60, 12000, 2000))) # Regione 1
  337. calibra(1, ((100, 12000, 2000), (60, 2100, 700))) # Regione 2
  338. calibra(2, ((100, 12000, 2000), (60, 2100, 700))) # Regione 3
  339. calibra(3, ((100, 12000, 2000), (60, 2100, 700))) # Regione 4
  340.  
  341.  
  342.  
  343. #%% Rette di calibrazione
  344.  
  345. def line(x,a,b):
  346.     return a*x+b
  347.  
  348.  
  349. picchiVeriSpesso = np.array((29, 93, 122, 122))
  350. picchiVeriSottile = np.array((93, 29, 6.4, 14.4))
  351.  
  352. picchiSpesso = np.array([i[0][1] for i in lista_popt])
  353. sigmaSpesso = np.array([i[0][2] for i in lista_popt])
  354. errPicchiSpesso = np.array([np.sqrt(i[0][1,1]) for i in lista_pcov])
  355.  
  356. picchiSottile = np.array([i[1][1] for i in lista_popt])
  357. sigmaSottile = np.array([i[1][2] for i in lista_popt])
  358. errPicchiSottile = np.array([np.sqrt(i[1][1,1]) for i in lista_pcov])
  359.  
  360.  
  361. # fit lineari
  362. popt1, pcov1 = curve_fit(line, picchiVeriSpesso, picchiSpesso,
  363.                          sigma = sigmaSpesso)
  364. popt2, pcov2 = curve_fit(line, picchiVeriSottile, picchiSottile,
  365.                          sigma = sigmaSottile)
  366.  
  367.  
  368. fig, ax = plt.subplots(1,2)#, sharex=True)
  369. fig.subplots_adjust(wspace=.3)
  370. fig.set_size_inches(12,5)
  371.  
  372.  
  373. # Plotto i punti
  374. ax[0].errorbar(picchiVeriSpesso, picchiSpesso, fmt=".",    
  375.                yerr = sigmaSpesso, label = "Dati" )
  376. ax[1].errorbar(picchiVeriSottile, picchiSottile, fmt=".",
  377.                yerr = sigmaSottile, label = "Dati" )
  378.  
  379. # Plotto le curve fittate
  380. ax[0].plot(picchiVeriSpesso, line(picchiVeriSpesso, *popt1), ":k")
  381. ax[1].plot(picchiVeriSottile, line(picchiVeriSottile, *popt2), ":k")
  382.  
  383. ax[0].set_title("Rivelatore spesso", fontsize = 16)
  384. ax[1].set_title("Rivelatore sottile", fontsize = 16)
  385.  
  386.  
  387.    
  388. for i in range(2):
  389.     ax[i].grid()
  390.     ax[i].legend(fontsize=14)
  391.     ax[i].set_ylabel("Picchi fittati (ADC)", fontsize=14)
  392.     ax[i].set_xlabel("Picchi veri (keV)", fontsize=14)
  393.  
  394.  
  395. if Export:
  396.     plt.savefig(f'../RELAZIONE/FIGURE/retta-cal.eps', format='eps')
  397. plt.show()
  398.  
  399.  
  400. myChi2_Spesso = np.sum((( picchiSpesso - line(picchiVeriSpesso, *popt1)) /sigmaSpesso )**2) /  (len(picchiVeriSpesso)-2)
  401. myChi2_Sottile = np.sum((( picchiSottile - line(picchiVeriSottile, *popt1)) /sigmaSottile )**2) /  (len(picchiVeriSottile)-2)
  402.  
  403.  
  404. print(f"Chi2 spesso: {myChi2_Spesso}")
  405. print(f"Chi2 sottile: {myChi2_Sottile}")
  406.  
  407.  
  408. #%% Stampo la tabella con i parametri delle gaussiane
  409. # Spesso
  410. for i in range(len(lista_popt)):
  411.     print(f"\t& {i+1}\t& {lista_popt[i][0][1]:.2f}\t& {np.sqrt(lista_pcov[i][0][1,1]):.2f}\t& {lista_popt[i][0][2]:.2f}\t& {np.sqrt(lista_pcov[i][0][2,2]):.2f}\t& {lista_chi2[i][0]:.2f} \\\\")
  412.  
  413. print("\n")
  414.  
  415. # Sottile
  416. for i in range(len(lista_popt)):
  417.     print(f"\t& {i+1}\t& {lista_popt[i][1][1]:.2f}\t& {np.sqrt(lista_pcov[i][1][1,1]):.2f}\t& {lista_popt[i][1][2]:.2f}\t& {np.sqrt(lista_pcov[i][1][2,2]):.2f}\t& {lista_chi2[i][1]:.2f} \\\\")
  418.  
  419.  
  420. #%%  NON LA METTIAMO Stampo una tabella con i picchi veri, quelli fittati in kev e errore
  421. # picco fittato kev, errore, picco vero
  422.  
  423. err_retta = lambda x,dx,m,dm,q,dq : np.sqrt((x*dm)**2 + (m*dx)**2 + dq**2)
  424.  
  425. # Spesso
  426. m = 1/popt1[0]
  427. q = -popt1[1]/popt1[0]
  428. errm = (1/(m**2)) * np.sqrt(pcov1[0,0])
  429. errq = np.sqrt((np.sqrt(pcov1[1,1]/m))**2 + ((q/(m**2))*np.sqrt(pcov1[0,0]))**2)
  430. for i in range(len(lista_popt)):
  431.     print(f"\t& {picchiVeriSpesso[i]}\t& {line(lista_popt[i][0][1], m, q):.2f}\t& {err_retta(picchiSpesso[i],errPicchiSpesso[i],m,errm,q,errq):.2f} \\\\")
  432.    
  433.  
  434. # Sottile
  435. m = 1/popt2[0]
  436. q = -popt2[1] / popt2[0]
  437. errm = (1/(m**2)) * np.sqrt(pcov2[0,0])
  438. errq = np.sqrt((np.sqrt(pcov2[1,1]/m))**2 + ((q/(m**2))*np.sqrt(pcov2[0,0]))**2)
  439.  
  440. for i in range(len(lista_popt)):
  441.     print(f"\t& {picchiVeriSottile[i]}\t& {line(lista_popt[i][1][1], m, q):.2f}\t& {err_retta(picchiSottile[i],errPicchiSottile[i],m,errm,q,errq):.2f} \\\\")
  442.  
  443.  
  444.  
  445. #%% Disegno le regioni considerate
  446.  
  447.  
  448.  
  449.  
  450. myRect = rectList[3]
  451.  
  452. startx = myRect.get_x()
  453. endx = startx + myRect.get_width()
  454.  
  455. idx_startx = np.where(xcentr >= startx)[0][0]
  456. idx_endx = np.where(xcentr >= endx)[0][0]
  457.  
  458.  
  459. starty = myRect.get_y()
  460. endy = starty + myRect.get_height()
  461.  
  462. idx_starty = np.where(ycentr >= starty)[0][0]
  463. idx_endy = np.where(ycentr >= endy)[0][0]
  464.  
  465.  
  466.  
  467. subdata = myPlot[0][ idx_startx:idx_endx , idx_starty:idx_endy]
  468. subdata_sqrt = np.sqrt(subdata)
  469.  
  470.  
  471.  
  472. # Quello che vede il rivelatore sottile, proietto a sx
  473. proiez_sottile = np.sum(subdata, axis = 0)
  474. proiez_sottile_sqrt = np.sqrt(np.sum(subdata_sqrt**2, axis = 0))
  475.  
  476. # Quello che vede il rivelatore spesso, proietto sotto
  477. proiez_spesso = np.sum(subdata, axis = 1)
  478. proiez_spesso_sqrt = np.sqrt(np.sum(subdata_sqrt**2, axis = 1))
  479.  
  480.  
  481.  
  482. vect_spesso = np.array([i for i in xcentr if (i >= startx) & (i<=endx)])
  483. vect_sottile = np.array([i for i in ycentr if (i >= starty) & (i<=endy)])
  484.  
  485.  
  486.  
  487. # fitto con gaussiane
  488. popt_spesso, pcov_spesso = curve_fit(gauss, vect_spesso, proiez_spesso,
  489.                                      sigma = proiez_spesso_sqrt,
  490.                                      p0 = lista_popt[3][0])
  491.  
  492. popt_sottile, pcov_sottile = curve_fit(gauss, vect_sottile, proiez_sottile,
  493.                                      sigma = proiez_sottile_sqrt,
  494.                                      p0 = lista_popt[3][1])
  495.  
  496.  
  497.  
  498.  
  499.  
  500.  
  501.  
  502. fig, ax = plt.subplots(2, 1)#, sharex=True)
  503. fig.subplots_adjust(hspace=.4)
  504. fig.set_size_inches(12,5)
  505.  
  506.  
  507. fig.suptitle(f"Regione 4 - aree per calcolo lifetime", fontsize=16)
  508.  
  509. # Plotto le proiezioni
  510. ax[0].plot(vect_spesso, proiez_spesso, "k", label = "Proiezione rivelatore spesso")
  511. ax[1].plot(vect_sottile, proiez_sottile, "k", label = "Proiezione rivelatore sottile")
  512.  
  513.  
  514. # Plotto le curve fittate
  515. ax[0].plot(vect_spesso, gauss(vect_spesso, *popt_spesso), ":k",
  516.            label = "Curva fittata")
  517. ax[1].plot(vect_sottile, gauss(vect_sottile, *popt_sottile), ":k",
  518.            label = "Curva fittata")
  519.  
  520.  
  521.  
  522.  
  523.  
  524. # Parte nuova
  525. # Disegno delle linee ad 1 sigma
  526. ax[0].axvline(lista_popt[3][0][1], ls = "--", c="grey")
  527.  
  528. ax[0].axvline(lista_popt[3][0][1] + .75*lista_popt[3][0][2], ls = "--", c="tab:blue", label = "0.75 Sigma")
  529. ax[0].axvline(lista_popt[3][0][1] - .75*lista_popt[3][0][2], ls = "--", c="tab:blue")
  530. ax[0].axvline(lista_popt[3][0][1] + 1*lista_popt[3][0][2], ls = "--", c="tab:orange", label = "1 Sigma")
  531. ax[0].axvline(lista_popt[3][0][1] - 1*lista_popt[3][0][2], ls = "--", c="tab:orange")
  532. ax[0].axvline(lista_popt[3][0][1] + 1.25*lista_popt[3][0][2], ls = "--", c="tab:green", label = "1.25 Sigma")
  533. ax[0].axvline(lista_popt[3][0][1] - 1.25*lista_popt[3][0][2], ls = "--", c="tab:green")
  534. ax[0].axvline(lista_popt[3][0][1] + 1.5*lista_popt[3][0][2], ls = "--", c="tab:red", label = "1.5 Sigma")
  535. ax[0].axvline(lista_popt[3][0][1] - 1.5*lista_popt[3][0][2], ls = "--", c="tab:red")
  536.  
  537.  
  538.  
  539. ax[1].axvline(lista_popt[3][1][1], ls = "--", c="grey")
  540.  
  541. ax[1].axvline(lista_popt[3][1][1] + .75*lista_popt[3][1][2], ls = "--", c="tab:blue", label = "0.75 Sigma")
  542. ax[1].axvline(lista_popt[3][1][1] - .75*lista_popt[3][1][2], ls = "--", c="tab:blue")
  543. ax[1].axvline(lista_popt[3][1][1] + 1*lista_popt[3][1][2], ls = "--", c="tab:orange", label = "1 Sigma")
  544. ax[1].axvline(lista_popt[3][1][1] - 1*lista_popt[3][1][2], ls = "--", c="tab:orange")
  545. ax[1].axvline(lista_popt[3][1][1] + 1.25*lista_popt[3][1][2], ls = "--", c="tab:green", label = "1.25 Sigma")
  546. ax[1].axvline(lista_popt[3][1][1] - 1.25*lista_popt[3][1][2], ls = "--", c="tab:green")
  547. ax[1].axvline(lista_popt[3][1][1] + 1.5*lista_popt[3][1][2], ls = "--", c="tab:red", label = "1.5 Sigma")
  548. ax[1].axvline(lista_popt[3][1][1] - 1.5*lista_popt[3][1][2], ls = "--", c="tab:red")
  549.  
  550.  
  551. for i in range(2):
  552.     ax[i].grid()
  553.     ax[i].legend(fontsize=10)
  554.     ax[i].set_ylabel("Conteggi", fontsize=14)
  555.     ax[i].set_xlabel("ADC", fontsize=14)
  556.  
  557.  
  558.  
  559. if Export:
  560.     plt.savefig(f'../RELAZIONE/FIGURE/aree.eps', format='eps')
  561. plt.show()
  562.  
  563.  
  564. #%% =================================================================
  565.  
  566.  
  567. #%% Da qui si calcola la lifetime
  568.  
  569. def myExp(x,a,tau):
  570.     return a*np.exp(-x/tau)
  571.  
  572. listaSigme = []
  573. listaLifetime = []
  574. listaErrLifetime = []
  575. nuovalistachi2 = []
  576.  
  577.  
  578.  
  579.  
  580.  
  581. def calcolaLifetime(nsigma, altezza, numBins):
  582.    
  583.     ADC1_min = lista_popt[3][0][1] - nsigma*lista_popt[3][0][2]
  584.     ADC1_max = lista_popt[3][0][1] + nsigma*lista_popt[3][0][2]
  585.     ADC2_min = lista_popt[3][1][1] - nsigma*lista_popt[3][1][2]
  586.     ADC2_max = lista_popt[3][1][1] + nsigma*lista_popt[3][1][2]
  587.  
  588.     listaTempi = []
  589.     for i in range(len(time)):
  590.         if (ADC1[i] >= ADC1_min) & (ADC1[i] <= ADC1_max) & \
  591.            (ADC2[i] >= ADC2_min) & (ADC2[i] <= ADC2_max) :
  592.                
  593.                listaTempi.append(np.abs(time1[i] - time2[i]))
  594.                
  595.              
  596.     listaTempi = np.array(listaTempi)
  597.  
  598.      
  599.  
  600.  
  601.    
  602.     h, bins = np.histogram(listaTempi, bins=numBins)
  603.     binc = bins[:-1] + bins[1] - bins[0]
  604.    
  605.     plt.plot(binc, h, ds = "steps-mid", label = f"Sigma = {nsigma}")
  606.    
  607.    
  608.     # fitto con un exp
  609.     condizione = binc > 480
  610.     xVect = binc[condizione]
  611.     yVect = h[condizione]
  612.    
  613.     popt, pcov = curve_fit(myExp, xVect, yVect, sigma = np.sqrt(yVect),
  614.                            p0 = [altezza, 1.41920784e+02])
  615.    
  616.     plt.plot(xVect, myExp(xVect,*popt), ls = ":", c="k")
  617.    
  618.  
  619.     print(f"Sigma: {nsigma} - Vita media: {popt[1]} più o meno {np.sqrt(pcov[1,1])}")
  620.  
  621.     listaSigme.append(nsigma)
  622.     listaLifetime.append(popt[1])  
  623.     listaErrLifetime.append(np.sqrt(pcov[1,1]))  
  624.    
  625.     # calcolo i chi2
  626.     tmpchi2 = np.sum(( (yVect - myExp(xVect,*popt)) / np.sqrt(yVect)  )**2)/(len(yVect)-2)
  627.     nuovalistachi2.append(tmpchi2)
  628.    
  629.    
  630.    
  631.    
  632. # Calcolo le lifetime  
  633.    
  634.    
  635. plt.figure(figsize=(12,5))
  636. plt.title(f"Lifetime", fontsize=16)
  637.  
  638.  
  639. calcolaLifetime(.75, 1.59507942e+05, 45)
  640. calcolaLifetime(1, 1.59507942e+05, 45)
  641. calcolaLifetime(1.25, 1.59507942e+05, 40)
  642. calcolaLifetime(1.5, 1.59507942e+05, 45)
  643.  
  644. plt.xlabel("Delta time (ns)", fontsize=14)
  645. plt.ylabel("Entries", fontsize=14)
  646. plt.legend(fontsize=14)
  647. plt.grid()
  648. if Export:
  649.     plt.savefig(r'../RELAZIONE/FIGURE/lifetime.eps', format='eps')
  650. plt.show()
  651.  
  652.  
  653.  
  654. listaSigme = np.array(listaSigme)
  655. listaLifetime = np.array(listaLifetime)
  656. listaErrLifetime = np.array(listaErrLifetime)
  657. # Testo la compatibilità con 141
  658. Zvect = np.abs(141-listaLifetime)/listaErrLifetime
  659.  
  660.  
  661. #%% Stampatabella lifetime
  662. for i in range(len(listaSigme)):
  663.     print(f"{listaSigme[i]}\t& {listaLifetime[i]:.2f}\t& {listaErrLifetime[i]:.2f}\t& {Zvect[i]:.2f}\t& {nuovalistachi2[i]:.2f} \\\\")
  664.  
  665.  
  666.  
  667.  
  668. #%% Calcolo la lifetime in modo bruttocon anche tutti i dati
  669. """
  670. Al plot di prima sovrapponiamo anche l'istogramma complessivo dei tempi
  671. """
  672.  
  673.  
  674. def calcolaLifetimebrutta(nsigma, altezza, numBins, primo=False):
  675.    
  676.     ADC1_min = lista_popt[3][0][1] - nsigma*lista_popt[3][0][2]
  677.     ADC1_max = lista_popt[3][0][1] + nsigma*lista_popt[3][0][2]
  678.     ADC2_min = lista_popt[3][1][1] - nsigma*lista_popt[3][1][2]
  679.     ADC2_max = lista_popt[3][1][1] + nsigma*lista_popt[3][1][2]
  680.  
  681.     listaTempi = []
  682.     listaTempiBrutto = []
  683.    
  684.     for i in range(len(time)):
  685.         if (ADC1[i] >= ADC1_min) & (ADC1[i] <= ADC1_max) & \
  686.            (ADC2[i] >= ADC2_min) & (ADC2[i] <= ADC2_max) :
  687.                
  688.                listaTempi.append(np.abs(time1[i] - time2[i]))
  689.                
  690.     if primo:
  691.          for i in range(len(time)):
  692.              listaTempiBrutto.append(np.abs(time1[i] - time2[i]))
  693.                
  694.              
  695.     listaTempi = np.array(listaTempi)
  696.  
  697.      
  698.  
  699.  
  700.    
  701.    
  702.    
  703.     if primo:
  704.         hb, binsb = np.histogram(listaTempiBrutto, bins=100)
  705.         bincb = binsb[:-1] + binsb[1] - binsb[0]
  706.         plt.plot(bincb, hb, ds = "steps-mid", label = f"Complessivo", c="k")
  707.    
  708.    
  709.  
  710.  
  711.     h, bins = np.histogram(listaTempi, bins=numBins)
  712.     binc = bins[:-1] + bins[1] - bins[0]
  713.    
  714.     plt.plot(binc, h, ds = "steps-mid", label = f"Sigma = {nsigma}")
  715.    
  716.    
  717.    
  718. # Calcolo le lifetime  
  719.    
  720.    
  721. plt.figure(figsize=(12,5))
  722. plt.title(f"Lifetime", fontsize=16)
  723.  
  724.  
  725. calcolaLifetimebrutta(.75, 1.59507942e+05, 45, True)
  726. calcolaLifetimebrutta(1, 1.59507942e+05, 45)
  727. calcolaLifetimebrutta(1.25, 1.59507942e+05, 40)
  728. calcolaLifetimebrutta(1.5, 1.59507942e+05, 45)
  729.  
  730. plt.xlabel("Delta time (ns)", fontsize=14)
  731. plt.ylabel("Entries", fontsize=14)
  732. plt.legend(fontsize=14)
  733. plt.yscale("log")
  734. plt.grid()
  735.  
  736.  
  737.  
  738. if Export:
  739.     plt.savefig(r'../RELAZIONE/FIGURE/lifetime-brutto.eps', format='eps')
  740. plt.show()
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!