Advertisement
Stex6299

CreaDati - Risoluzione - AutoGUI

May 1st, 2021
1,259
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Python 33.43 KB | None | 0 0
  1. # -*- coding: utf-8 -*-
  2. """
  3. Created on Sun Apr 25 10:44:22 2021
  4.  
  5. Autore: Stefano Carsi
  6. Data ultima modifica: 202104271819
  7.  
  8.    
  9. - Cosa ho fatto:
  10.    *
  11.    
  12. - Cosa devo fare:
  13.    *
  14.    
  15. - Per relazione:
  16.    *
  17.    
  18.    
  19.  
  20. - Info per l'uso:
  21.    > Serve un pizzico di manualità
  22.    
  23.    > Runnare la cella, aggiustare le cose e provare a far convergere
  24.      il plot
  25.      
  26.    > Solo dopo che si è soddisfatti, runnare la cella sottostante per memo-
  27.      rizzarsi i parametri ottimizzati dal fit
  28.      
  29.    > Infine esportare il file con tutti i dati
  30.    
  31.    > Le figure dei fit, se serve, meglio esportarle a mano
  32.    
  33.    > Nota bene: il vettore degli starting point non può essere vuoto!
  34.  
  35.  
  36.  
  37.    
  38. """
  39.  
  40.  
  41. #INIZIO STANDARD
  42. #%matplotlib auto
  43. import matplotlib.pyplot as plt
  44. import numpy as np
  45. from scipy.optimize import curve_fit
  46.  
  47. from matplotlib import use as mplUse
  48. mplUse("Qt5Agg")
  49. #%matplotlib inline
  50.  
  51.  
  52. import os
  53. path = r"F:\Corsi\Laboratorio 4 - Fisica nucleare e subnucleare\3 - Spettroscopia gamma\PITONE"
  54. os.chdir(path)
  55.  
  56. Export = False
  57.  
  58.  
  59.  
  60.  
  61. #%% Carico i dati
  62.  
  63. listaMatrici = []
  64. listaNum = []
  65.  
  66. for e in os.scandir(r"..\datadir_voltagescan_2021"):
  67.     if e.name.split(".")[1] != "dat": continue
  68.  
  69.     tmp = np.loadtxt(e.path)
  70.     listaMatrici.append(tmp.copy())
  71.     listaNum.append(int(e.name.split(".")[0][-2:]))
  72.    
  73. #listaNum = np.array(listaNum, dtype = np.int8)    
  74. zipIt = zip(listaNum, range(len(listaNum)))
  75. myDict = dict(zipIt)
  76.    
  77.  
  78.  
  79.  
  80. #%% Inizio - Modelli di fitting
  81.  
  82. from dataphile.statistics.regression.modeling import Parameter, Model, CompositeModel, AutoGUI
  83.  
  84.  
  85. def gauss(x,a,mu,sigma):
  86.     return a/(sigma*np.sqrt(2*np.pi))* np.exp(-(x-mu)**2 / (2*sigma**2))
  87.  
  88. def line(x, m, q):
  89.     return m*x+q
  90.  
  91. def exp(x,a,k):
  92.     return a*np.exp(-k*x)
  93.  
  94.  
  95.  
  96. #%% 52 - funziona, ma non usato
  97. numRun = 52
  98. dati = listaMatrici[myDict[numRun]][:,1]
  99. dati = dati[dati<4095]
  100.  
  101. h, bins = np.histogram(dati, bins=100)
  102. binc = bins[:-1] + (bins[1] - bins[0])
  103.    
  104.    
  105. fig = plt.figure()
  106. fig.set_size_inches(12,5)
  107. fig.suptitle("Titolo", fontsize=16)
  108.  
  109. ax = fig.add_axes([0.15, 0.14, 0.84, 0.70])
  110. ax.set_yscale("log")
  111. ax.grid(True)
  112. ax.set_xlabel("ADC")
  113. ax.set_ylabel("Entries")
  114.  
  115.  
  116.  
  117. # Plotto tutti i dati
  118. ax.plot(binc, h, color='black', lw=1, label='Dati completi', ds="steps-mid")
  119. # Plotto i dati da fittare
  120.  
  121. cond = (binc>1273) & (binc<1560)
  122. xVect = binc[cond]
  123. yVect = h[cond]
  124.  
  125. data_graph, = ax.plot(xVect, yVect, lw=1, label='Dati per fit',\
  126.                       ds="steps-mid", c="tab:green")
  127.  
  128.  
  129.  
  130.  
  131. model = CompositeModel(
  132.     Model(gauss,
  133.           Parameter(value=443026, bounds=(3e5, 3e6),      label='A'),
  134.           Parameter(value=1434.95,   bounds=(0, 4095),   label='mu'),
  135.           Parameter(value=42.4055,   bounds=(0, 300), label='sigma'),
  136.           label='gaussiana'),
  137.     Model(line,
  138.           Parameter(value=-2.54899, bounds=(-10, 10),  label='m'),
  139.           Parameter(value=4205.7, bounds=(-100, 10000), label='q'),
  140.           label='retta'),
  141.     label='gaussian_peaks')
  142.  
  143.  
  144. xsample = np.linspace(0, 4095, 1500)
  145. model_graph, = ax.plot(xsample, model(xsample), color='steelblue', label='model')
  146. ax.legend();
  147.  
  148.  
  149. fig
  150.  
  151.  
  152. ax.set_position([0.15, 0.30, 0.84, 0.56])
  153. gui = AutoGUI(model, [model_graph], bbox=[0.20, 0.07, 0.75, 0.12], figure=fig,
  154.               slider_options={'color': 'steelblue'}, data=(xVect, yVect));
  155.  
  156.  
  157. A_Riassunto = model.summary()
  158.  
  159.  
  160.  
  161. #%% 52 Ufficiale (old)
  162.  
  163. numRun = 52
  164. dati = listaMatrici[myDict[numRun]][:,1]
  165. dati = dati[dati<4095]
  166.  
  167. h, bins = np.histogram(dati, bins=150)
  168. binc = bins[:-1] + (bins[1] - bins[0])
  169.    
  170.    
  171. fig = plt.figure()
  172. fig.set_size_inches(12,5)
  173. fig.suptitle("226Ra - 600V - Gain: 1000 ", fontsize=16)
  174.  
  175. ax = fig.add_axes([0.15, 0.14, 0.84, 0.70])
  176. #ax.set_yscale("log")
  177. ax.grid(True)
  178. ax.set_xlabel("ADC", fontsize=14)
  179. ax.set_ylabel("Entries", fontsize=14)
  180.  
  181.  
  182.  
  183. # Plotto tutti i dati
  184. ax.plot(binc, h, color='black', lw=1, label='Dati completi', ds="steps-mid")
  185. # Plotto i dati da fittare
  186.  
  187. cond = (binc<1560) & (binc>243)
  188. xVect = binc[cond]
  189. yVect = h[cond]
  190.  
  191. data_graph, = ax.plot(xVect, yVect, lw=1, label='Dati per fit',\
  192.                       ds="steps-mid", c="tab:green")
  193.  
  194.  
  195.  
  196.  
  197. model = CompositeModel(
  198.     Model(gauss,
  199.           Parameter(value=463743, bounds=(1e5, 3e6),      label='A'),
  200.           Parameter(value=304.195,   bounds=(0, 4095),   label='mu'),
  201.           Parameter(value=16.2228,   bounds=(0, 100), label='sigma'),
  202.           label='gaussiana 1'),
  203.     Model(gauss,
  204.           Parameter(value=128113, bounds=(1e5, 3e6),      label='A'),
  205.           Parameter(value=536.718,   bounds=(0, 4095),   label='mu'),
  206.           Parameter(value=21.1004,   bounds=(0, 100), label='sigma'),
  207.           label='gaussiana 2'),
  208.     Model(gauss,
  209.           Parameter(value=134248, bounds=(1e5, 3e6),      label='A'),
  210.           Parameter(value=655.77,   bounds=(0, 4095),   label='mu'),
  211.           Parameter(value=21.771,   bounds=(0, 100), label='sigma'),
  212.           label='gaussiana 3'),
  213.     Model(gauss,
  214.           Parameter(value=296409, bounds=(1e5, 3e6),      label='A'),
  215.           Parameter(value=768.839,   bounds=(0, 4095),   label='mu'),
  216.           Parameter(value=25.0682,   bounds=(0, 100), label='sigma'),
  217.           label='gaussiana 4'),
  218.     Model(gauss,
  219.           Parameter(value=474508, bounds=(1e5, 3e6),      label='A'),
  220.           Parameter(value=888.935,   bounds=(0, 4095),   label='mu'),
  221.           Parameter(value=28.5505,   bounds=(0, 100), label='sigma'),
  222.           label='gaussiana 5'),
  223.     Model(gauss,
  224.           Parameter(value=276598, bounds=(1e5, 3e6),      label='A'),
  225.           Parameter(value=1427.27,   bounds=(0, 4095),   label='mu'),
  226.           Parameter(value=40.4266,   bounds=(0, 100), label='sigma'),
  227.           label='gaussiana 6'),
  228.  
  229. #    Model(line,
  230. #          Parameter(value=-4.55851, bounds=(-10, 10),  label='m'),
  231. #          Parameter(value=5881.62, bounds=(-100, 10000), label='q'),
  232. #          label='retta'),
  233.     Model(exp,
  234.           Parameter(value=9142.08, bounds=(0, 1e5),  label='a'),
  235.           Parameter(value=0.00208337, bounds=(0, 1e-2), label='k'),
  236.           label='exp'),
  237.     label='gaussian_peaks')
  238.  
  239.  
  240. xsample = np.linspace(0, 4095, 1500)
  241. model_graph, = ax.plot(xsample, model(xsample), color='steelblue', \
  242.                        label='Curva fittata')
  243. ax.legend();
  244.  
  245.  
  246. fig
  247.  
  248.  
  249. ax.set_position([0.15, 0.30, 0.84, 0.56])
  250. gui = AutoGUI(model, [model_graph], bbox=[0.20, 0.07, 0.75, 0.12], figure=fig,
  251.               slider_options={'color': 'steelblue'}, data=(xVect, yVect));
  252.  
  253. # Eseguire a mano dopo aver premuto il pulsantone ed ispezionare a mano
  254. # è facilmente castabile in un array di numpy
  255. Riassunto_52 = model.summary()
  256. Riassunto52 = np.array(Riassunto_52)
  257.  
  258. if Export: pass
  259.     #plt.savefig(r'../RELAZIONE/FIGURE/risoluzione-52.eps', format='eps')
  260. plt.show()
  261.  
  262.  
  263.  
  264. #%% 52 nuovo ufficiale OK
  265. # Nota, ho alzato il bin150->200, ora converge
  266.  
  267. # Carico i dati
  268. numRun = 52
  269. dati = listaMatrici[myDict[numRun]][:,1]
  270. dati = dati[dati<4095]
  271.  
  272.  
  273. # Costruisco l'istogramma e trovo il centro
  274. h, bins = np.histogram(dati, bins=200)
  275. binc = bins[:-1] + (bins[1] - bins[0])
  276.    
  277.    
  278. # Creo figura e assi su cui disegnare
  279. fig = plt.figure()
  280. fig.set_size_inches(12,5)
  281. fig.suptitle("226Ra - 600V - Gain: 1000 ", fontsize=16)
  282.  
  283. ax = fig.add_axes([0.15, 0.14, 0.84, 0.70])
  284. ax.set_yscale("log")
  285. ax.grid(True)
  286. ax.set_xlabel("ADC", fontsize=14)
  287. ax.set_ylabel("Entries", fontsize=14)
  288.  
  289.  
  290.  
  291. # Plotto tutti i dati
  292. ax.plot(binc, h, color='black', lw=1, label='Dati completi', ds="steps-mid")
  293.  
  294.  
  295. # Seleziono e Plotto i dati da fittare
  296. cond = (binc<2633) & (binc>243)
  297. xVect = binc[cond]
  298. yVect = h[cond]
  299.  
  300. data_graph, = ax.plot(xVect, yVect, lw=1, label='Dati per fit',\
  301.                       ds="steps-mid", c="tab:green")
  302.  
  303.  
  304.    
  305. # Starting points: inseriti copiando quelli di convergenza
  306. startPars = np.array([
  307.        3.14020187e+05, 2.98589322e+02, 1.43257226e+01,
  308.        9.41411082e+04, 5.34214983e+02, 2.02581202e+01,
  309.        1.04459955e+05, 6.51821651e+02, 2.23514168e+01,
  310.        
  311.        2.30961972e+05, 7.65474019e+02, 2.49634128e+01,
  312.        3.65195534e+05, 8.86502718e+02, 2.82249031e+01,
  313.        2.14699579e+05, 1.42224403e+03, 4.12398651e+01,
  314.        
  315.        1.13048417e+05, 1.76150564e+03, 1.49847505e+02,
  316.        6.24648382e+04, 2.11018949e+03, 1.12551885e+02,
  317.        9.82593192e+04, 2.47588684e+03, 1.06711679e+02,
  318.        
  319.        7.81582354e+03, 2.32745796e-03])
  320.  
  321.  
  322. # Movello di fitting
  323. model = CompositeModel(
  324.     Model(gauss,
  325.           Parameter(value=startPars[0], bounds=(1e5, 3e6),      label='A'),
  326.           Parameter(value=startPars[1],   bounds=(0, 4095),   label='mu'),
  327.           Parameter(value=startPars[2],   bounds=(0, 100), label='sigma'),
  328.           label='gaussiana 1'),
  329.     Model(gauss,
  330.           Parameter(value=startPars[3], bounds=(1e5, 3e6),      label='A'),
  331.           Parameter(value=startPars[4],   bounds=(0, 4095),   label='mu'),
  332.           Parameter(value=startPars[5],   bounds=(0, 100), label='sigma'),
  333.           label='gaussiana 2'),
  334.     Model(gauss,
  335.           Parameter(value=startPars[6], bounds=(1e5, 3e6),      label='A'),
  336.           Parameter(value=startPars[7],   bounds=(0, 4095),   label='mu'),
  337.           Parameter(value=startPars[8],   bounds=(0, 100), label='sigma'),
  338.           label='gaussiana 3'),
  339.    
  340.     Model(gauss,
  341.           Parameter(value=startPars[9], bounds=(1e5, 3e6),      label='A'),
  342.           Parameter(value=startPars[10],   bounds=(0, 4095),   label='mu'),
  343.           Parameter(value=startPars[11],   bounds=(0, 100), label='sigma'),
  344.           label='gaussiana 4'),
  345.     Model(gauss,
  346.           Parameter(value=startPars[12], bounds=(1e5, 3e6),      label='A'),
  347.           Parameter(value=startPars[13],   bounds=(0, 4095),   label='mu'),
  348.           Parameter(value=startPars[14],   bounds=(0, 100), label='sigma'),
  349.           label='gaussiana 5'),
  350.     Model(gauss,
  351.           Parameter(value=startPars[15], bounds=(1e5, 3e6),      label='A'),
  352.           Parameter(value=startPars[16],   bounds=(0, 4095),   label='mu'),
  353.           Parameter(value=startPars[17],   bounds=(0, 100), label='sigma'),
  354.           label='gaussiana 6'),
  355.    
  356.     Model(gauss,
  357.           Parameter(value=startPars[18], bounds=(1e4, 3e6),      label='A'),
  358.           Parameter(value=startPars[19],   bounds=(0, 4095),   label='mu'),
  359.           Parameter(value=startPars[20],   bounds=(0, 100), label='sigma'),
  360.           label='gaussiana 7'),
  361.     Model(gauss,
  362.           Parameter(value=startPars[21], bounds=(1e4, 3e6),      label='A'),
  363.           Parameter(value=startPars[22],   bounds=(0, 4095),   label='mu'),
  364.           Parameter(value=startPars[23],   bounds=(0, 100), label='sigma'),
  365.           label='gaussiana 8'),
  366.     Model(gauss,
  367.           Parameter(value=startPars[24], bounds=(1e4, 3e6),      label='A'),
  368.           Parameter(value=startPars[25],   bounds=(0, 4095),   label='mu'),
  369.           Parameter(value=startPars[26],   bounds=(0, 100), label='sigma'),
  370.           label='gaussiana 9'),
  371.  
  372. #    Model(line,
  373. #          Parameter(value=-4.55851, bounds=(-10, 10),  label='m'),
  374. #          Parameter(value=5881.62, bounds=(-100, 10000), label='q'),
  375. #          label='retta'),
  376.     Model(exp,
  377.           Parameter(value=startPars[27], bounds=(0, 1e5),  label='a'),
  378.           Parameter(value=startPars[28], bounds=(0, 1e-2), label='k'),
  379.           label='exp'),
  380.     label='gaussian_peaks')
  381.  
  382.  
  383. # Disegno la curva da muovere
  384. xsample = np.linspace(0, 4095, 1500)
  385. model_graph, = ax.plot(xsample, model(xsample), color='steelblue', \
  386.                        label='Curva fittata')
  387. ax.legend();
  388.  
  389.  
  390. # Mostro la figura
  391. fig
  392.  
  393.  
  394. # Ridimensiono gli assi per fare spazio ai cursori
  395. ax.set_position([0.10, 0.30, 0.84, 0.62])
  396.  
  397. # Chiamo ciò che fa la magia
  398. gui = AutoGUI(model, [model_graph], bbox=[0.20, 0.07, 0.75, 0.16], figure=fig,
  399.               slider_options={'color': 'steelblue'}, data=(xVect, yVect),
  400.               sigma=np.sqrt(yVect));
  401.  
  402. #%%
  403. # Eseguire a mano dopo aver premuto il pulsantone ed ispezionare a mano
  404. # è facilmente castabile in un array di numpy
  405. Riassunto_52 = model.summary()
  406. Riassunto52 = np.array(Riassunto_52)
  407.  
  408. if Export: pass
  409.     #plt.savefig(r'../RELAZIONE/FIGURE/risoluzione-52.eps', format='eps')
  410. #plt.show()
  411.  
  412.  
  413.  
  414.  
  415.  
  416.  
  417. #%% 53 - Fatto OK
  418. # Nota, ho alzato il bin150->200, ora converge
  419.  
  420. # Carico i dati
  421. numRun = 53
  422. dati = listaMatrici[myDict[numRun]][:,1]
  423. dati = dati[dati<4095]
  424.  
  425.  
  426. # Costruisco l'istogramma e trovo il centro
  427. h, bins = np.histogram(dati, bins=200)
  428. binc = bins[:-1] + (bins[1] - bins[0])
  429.    
  430.    
  431. # Creo figura e assi su cui disegnare
  432. fig = plt.figure()
  433. fig.set_size_inches(12,5)
  434. fig.suptitle("226Ra - 650V - Gain: 500 ", fontsize=16)
  435.  
  436. ax = fig.add_axes([0.15, 0.14, 0.84, 0.70])
  437. ax.set_yscale("log")
  438. ax.grid(True)
  439. ax.set_xlabel("ADC", fontsize=14)
  440. ax.set_ylabel("Entries", fontsize=14)
  441.  
  442.  
  443.  
  444. # Plotto tutti i dati
  445. ax.plot(binc, h, color='black', lw=1, label='Dati completi', ds="steps-mid")
  446.  
  447.  
  448. # Seleziono e Plotto i dati da fittare
  449. cond = (binc<2360) & (binc>220)
  450. xVect = binc[cond]
  451. yVect = h[cond]
  452.  
  453. data_graph, = ax.plot(xVect, yVect, lw=1, label='Dati per fit',\
  454.                       ds="steps-mid", c="tab:green")
  455.  
  456.  
  457.  
  458. # Starting points: inseriti copiando quelli di convergenza
  459. startPars = np.array([
  460.        3.42507160e+05, 2.80817395e+02, 1.24616954e+01,
  461.        1.03919529e+05, 4.90186785e+02, 2.03356879e+01,
  462.        1.18326472e+05, 5.99110882e+02, 2.07618554e+01,
  463.        
  464.        2.37704890e+05, 7.02525659e+02, 2.16424925e+01,
  465.        3.85261776e+05, 8.10998553e+02, 2.55724833e+01,
  466.        2.21297335e+05, 1.29723065e+03, 3.59272711e+01,
  467.        
  468.        1.12844278e+05, 1.60399837e+03, 1.34508650e+02,
  469.        7.36554529e+04, 1.92969659e+03, 1.30558647e+02,
  470.        8.67953207e+04, 2.24743304e+03, 7.91527781e+01,
  471.        
  472.        8.00390614e+03, 2.48022403e-03])
  473.  
  474.  
  475. # Movello di fitting
  476. model = CompositeModel(
  477.     Model(gauss,
  478.           Parameter(value=startPars[0], bounds=(1e5, 3e6),      label='A'),
  479.           Parameter(value=startPars[1],   bounds=(0, 4095),   label='mu'),
  480.           Parameter(value=startPars[2],   bounds=(0, 100), label='sigma'),
  481.           label='gaussiana 1'),
  482.     Model(gauss,
  483.           Parameter(value=startPars[3], bounds=(1e5, 3e6),      label='A'),
  484.           Parameter(value=startPars[4],   bounds=(0, 4095),   label='mu'),
  485.           Parameter(value=startPars[5],   bounds=(0, 100), label='sigma'),
  486.           label='gaussiana 2'),
  487.     Model(gauss,
  488.           Parameter(value=startPars[6], bounds=(1e5, 3e6),      label='A'),
  489.           Parameter(value=startPars[7],   bounds=(0, 4095),   label='mu'),
  490.           Parameter(value=startPars[8],   bounds=(0, 100), label='sigma'),
  491.           label='gaussiana 3'),
  492.    
  493.     Model(gauss,
  494.           Parameter(value=startPars[9], bounds=(1e5, 3e6),      label='A'),
  495.           Parameter(value=startPars[10],   bounds=(0, 4095),   label='mu'),
  496.           Parameter(value=startPars[11],   bounds=(0, 100), label='sigma'),
  497.           label='gaussiana 4'),
  498.     Model(gauss,
  499.           Parameter(value=startPars[12], bounds=(1e5, 3e6),      label='A'),
  500.           Parameter(value=startPars[13],   bounds=(0, 4095),   label='mu'),
  501.           Parameter(value=startPars[14],   bounds=(0, 100), label='sigma'),
  502.           label='gaussiana 5'),
  503.     Model(gauss,
  504.           Parameter(value=startPars[15], bounds=(1e5, 3e6),      label='A'),
  505.           Parameter(value=startPars[16],   bounds=(0, 4095),   label='mu'),
  506.           Parameter(value=startPars[17],   bounds=(0, 100), label='sigma'),
  507.           label='gaussiana 6'),
  508.    
  509.     Model(gauss,
  510.           Parameter(value=startPars[18], bounds=(1e4, 3e6),      label='A'),
  511.           Parameter(value=startPars[19],   bounds=(0, 4095),   label='mu'),
  512.           Parameter(value=startPars[20],   bounds=(0, 100), label='sigma'),
  513.           label='gaussiana 7'),
  514.     Model(gauss,
  515.           Parameter(value=startPars[21], bounds=(1e4, 3e6),      label='A'),
  516.           Parameter(value=startPars[22],   bounds=(0, 4095),   label='mu'),
  517.           Parameter(value=startPars[23],   bounds=(0, 100), label='sigma'),
  518.           label='gaussiana 8'),
  519.     Model(gauss,
  520.           Parameter(value=startPars[24], bounds=(1e4, 3e6),      label='A'),
  521.           Parameter(value=startPars[25],   bounds=(0, 4095),   label='mu'),
  522.           Parameter(value=startPars[26],   bounds=(0, 100), label='sigma'),
  523.           label='gaussiana 9'),
  524.  
  525. #    Model(line,
  526. #          Parameter(value=-4.55851, bounds=(-10, 10),  label='m'),
  527. #          Parameter(value=5881.62, bounds=(-100, 10000), label='q'),
  528. #          label='retta'),
  529.     Model(exp,
  530.           Parameter(value=startPars[27], bounds=(0, 1e5),  label='a'),
  531.           Parameter(value=startPars[28], bounds=(0, 1e-2), label='k'),
  532.           label='exp'),
  533.     label='gaussian_peaks')
  534.  
  535.  
  536. # Disegno la curva da muovere
  537. xsample = np.linspace(0, 4095, 1500)
  538. model_graph, = ax.plot(xsample, model(xsample), color='steelblue', \
  539.                        label='Curva fittata')
  540. ax.legend();
  541.  
  542.  
  543. # Mostro la figura
  544. fig
  545.  
  546.  
  547. # Ridimensiono gli assi per fare spazio ai cursori
  548. ax.set_position([0.10, 0.30, 0.84, 0.62])
  549.  
  550. # Chiamo ciò che fa la magia
  551. gui = AutoGUI(model, [model_graph], bbox=[0.20, 0.07, 0.75, 0.16], figure=fig,
  552.               slider_options={'color': 'steelblue'}, data=(xVect, yVect),
  553.               sigma=np.sqrt(yVect));
  554.  
  555. #%%
  556. # Eseguire a mano dopo aver premuto il pulsantone ed ispezionare a mano
  557. # è facilmente castabile in un array di numpy
  558. Riassunto_53 = model.summary()
  559. Riassunto53 = np.array(Riassunto_53)
  560.  
  561.  
  562.  
  563.  
  564.  
  565.  
  566. #%% 54 - Fatto
  567. # I parametri soo fatti sui 200 bin
  568.  
  569. # Carico i dati
  570. numRun = 54
  571. dati = listaMatrici[myDict[numRun]][:,1]
  572. dati = dati[dati<4095]
  573.  
  574.  
  575. # Costruisco l'istogramma e trovo il centro
  576. h, bins = np.histogram(dati, bins=200)
  577. binc = bins[:-1] + (bins[1] - bins[0])
  578.    
  579.    
  580. # Creo figura e assi su cui disegnare
  581. fig = plt.figure()
  582. fig.set_size_inches(12,5)
  583. fig.suptitle("226Ra - 700V - Gain: 200 ", fontsize=16)
  584.  
  585. ax = fig.add_axes([0.15, 0.14, 0.84, 0.70])
  586. ax.set_yscale("log")
  587. ax.grid(True)
  588. ax.set_xlabel("ADC", fontsize=14)
  589. ax.set_ylabel("Entries", fontsize=14)
  590.  
  591.  
  592.  
  593. # Plotto tutti i dati
  594. ax.plot(binc, h, color='black', lw=1, label='Dati completi', ds="steps-mid")
  595.  
  596.  
  597. # Seleziono e Plotto i dati da fittare
  598. cond = (binc<1705) & (binc>166)
  599. xVect = binc[cond]
  600. yVect = h[cond]
  601.  
  602. data_graph, = ax.plot(xVect, yVect, lw=1, label='Dati per fit',\
  603.                       ds="steps-mid", c="tab:green")
  604.  
  605.    
  606.  
  607. # Starting points: inseriti copiando quelli di convergenza
  608. startPars = np.array([
  609.        4.62943107e+05, 2.26198954e+02, 1.42027020e+01,
  610.        4.30064153e+05, 3.51438064e+02, 6.29233950e+01,
  611.        1.10434659e+05, 4.57722754e+02, 1.32710721e+01,
  612.        
  613.        2.80277831e+05, 5.27133855e+02, 1.74496915e+01,
  614.        4.17788000e+05, 6.04375499e+02, 1.88606343e+01,
  615.        2.22374631e+05, 9.48896242e+02, 2.56032992e+01,
  616.        
  617.        7.15420610e+04, 1.16368668e+03, 7.36926220e+01,
  618.        6.46185002e+04, 1.38351571e+03, 8.84068589e+01,
  619.        8.02610692e+04, 1.61800918e+03, 5.50063439e+01,
  620.        
  621.        6.00904050e+03, 2.55138106e-03])
  622.  
  623.  
  624. # Movello di fitting
  625. model = CompositeModel(
  626.     Model(gauss,
  627.           Parameter(value=startPars[0], bounds=(1e5, 3e6),      label='A'),
  628.           Parameter(value=startPars[1],   bounds=(0, 4095),   label='mu'),
  629.           Parameter(value=startPars[2],   bounds=(0, 100), label='sigma'),
  630.           label='gaussiana 1'),
  631.     Model(gauss,
  632.           Parameter(value=startPars[3], bounds=(1e5, 3e6),      label='A'),
  633.           Parameter(value=startPars[4],   bounds=(0, 4095),   label='mu'),
  634.           Parameter(value=startPars[5],   bounds=(0, 100), label='sigma'),
  635.           label='gaussiana 2'),
  636.     Model(gauss,
  637.           Parameter(value=startPars[6], bounds=(1e5, 3e6),      label='A'),
  638.           Parameter(value=startPars[7],   bounds=(0, 4095),   label='mu'),
  639.           Parameter(value=startPars[8],   bounds=(0, 100), label='sigma'),
  640.           label='gaussiana 3'),
  641.    
  642.     Model(gauss,
  643.           Parameter(value=startPars[9], bounds=(1e5, 3e6),      label='A'),
  644.           Parameter(value=startPars[10],   bounds=(0, 4095),   label='mu'),
  645.           Parameter(value=startPars[11],   bounds=(0, 100), label='sigma'),
  646.           label='gaussiana 4'),
  647.     Model(gauss,
  648.           Parameter(value=startPars[12], bounds=(1e5, 3e6),      label='A'),
  649.           Parameter(value=startPars[13],   bounds=(0, 4095),   label='mu'),
  650.           Parameter(value=startPars[14],   bounds=(0, 100), label='sigma'),
  651.           label='gaussiana 5'),
  652.     Model(gauss,
  653.           Parameter(value=startPars[15], bounds=(1e5, 3e6),      label='A'),
  654.           Parameter(value=startPars[16],   bounds=(0, 4095),   label='mu'),
  655.           Parameter(value=startPars[17],   bounds=(0, 100), label='sigma'),
  656.           label='gaussiana 6'),
  657.    
  658.     Model(gauss,
  659.           Parameter(value=startPars[18], bounds=(1e4, 3e6),      label='A'),
  660.           Parameter(value=startPars[19],   bounds=(0, 4095),   label='mu'),
  661.           Parameter(value=startPars[20],   bounds=(0, 100), label='sigma'),
  662.           label='gaussiana 7'),
  663.     Model(gauss,
  664.           Parameter(value=startPars[21], bounds=(1e4, 3e6),      label='A'),
  665.           Parameter(value=startPars[22],   bounds=(0, 4095),   label='mu'),
  666.           Parameter(value=startPars[23],   bounds=(0, 100), label='sigma'),
  667.           label='gaussiana 8'),
  668.     Model(gauss,
  669.           Parameter(value=startPars[24], bounds=(1e4, 3e6),      label='A'),
  670.           Parameter(value=startPars[25],   bounds=(0, 4095),   label='mu'),
  671.           Parameter(value=startPars[26],   bounds=(0, 100), label='sigma'),
  672.           label='gaussiana 9'),
  673.  
  674. #    Model(line,
  675. #          Parameter(value=-4.55851, bounds=(-10, 10),  label='m'),
  676. #          Parameter(value=5881.62, bounds=(-100, 10000), label='q'),
  677. #          label='retta'),
  678.     Model(exp,
  679.           Parameter(value=startPars[27], bounds=(0, 1e5),  label='a'),
  680.           Parameter(value=startPars[28], bounds=(0, 1e-2), label='k'),
  681.           label='exp'),
  682.     label='gaussian_peaks')
  683.  
  684.  
  685. # Disegno la curva da muovere
  686. xsample = np.linspace(0, 4095, 1500)
  687. model_graph, = ax.plot(xsample, model(xsample), color='steelblue', \
  688.                        label='Curva fittata')
  689. ax.legend();
  690.  
  691.  
  692. # Mostro la figura
  693. fig
  694.  
  695.  
  696. # Ridimensiono gli assi per fare spazio ai cursori
  697. ax.set_position([0.10, 0.30, 0.84, 0.62])
  698.  
  699. # Chiamo ciò che fa la magia
  700. gui = AutoGUI(model, [model_graph], bbox=[0.20, 0.07, 0.75, 0.16], figure=fig,
  701.               slider_options={'color': 'steelblue'}, data=(xVect, yVect),
  702.               sigma=np.sqrt(yVect));
  703.  
  704. #%%
  705. # Eseguire a mano dopo aver premuto il pulsantone ed ispezionare a mano
  706. # è facilmente castabile in un array di numpy
  707. Riassunto_54 = model.summary()
  708. Riassunto54 = np.array(Riassunto_54)
  709.  
  710.  
  711.  
  712.  
  713.  
  714.  
  715. #%% 55 - Fatto OK
  716. # Nota, ho alzato il bin150->200, ora converge
  717.  
  718. # Carico i dati
  719. numRun = 55
  720. dati = listaMatrici[myDict[numRun]][:,1]
  721. dati = dati[dati<4095]
  722.  
  723.  
  724. # Costruisco l'istogramma e trovo il centro
  725. h, bins = np.histogram(dati, bins=200)
  726. binc = bins[:-1] + (bins[1] - bins[0])
  727.    
  728.    
  729. # Creo figura e assi su cui disegnare
  730. fig = plt.figure()
  731. fig.set_size_inches(12,5)
  732. fig.suptitle("226Ra - 750V - Gain: 200 ", fontsize=16)
  733.  
  734. ax = fig.add_axes([0.15, 0.14, 0.84, 0.70])
  735. ax.set_yscale("log")
  736. ax.grid(True)
  737. ax.set_xlabel("ADC", fontsize=14)
  738. ax.set_ylabel("Entries", fontsize=14)
  739.  
  740.  
  741.  
  742. # Plotto tutti i dati
  743. ax.plot(binc, h, color='black', lw=1, label='Dati completi', ds="steps-mid")
  744.  
  745.  
  746. # Seleziono e Plotto i dati da fittare
  747. cond = (binc<2803) & (binc>243)
  748. xVect = binc[cond]
  749. yVect = h[cond]
  750.  
  751. data_graph, = ax.plot(xVect, yVect, lw=1, label='Dati per fit',\
  752.                       ds="steps-mid", c="tab:green")
  753.  
  754.  
  755.  
  756. # Starting points: inseriti copiando quelli di convergenza
  757. startPars = np.array([
  758.        3.27007341e+05, 3.13976549e+02, 1.51769090e+01,
  759.        1.10288407e+05, 5.66753954e+02, 2.90578040e+01,
  760.        1.16014769e+05, 6.98831507e+02, 2.36358556e+01,
  761.        
  762.        2.39377275e+05, 8.23293555e+02, 2.62860390e+01,
  763.        3.79262612e+05, 9.53982266e+02, 2.88812872e+01,
  764.        2.18900196e+05, 1.53509437e+03, 4.27045957e+01,
  765.        
  766.        1.04806678e+05, 1.88950564e+03, 1.57225121e+02,
  767.        7.48170187e+04, 2.27747309e+03, 1.51051278e+02,
  768.        8.90622408e+04, 2.65905490e+03, 9.61389581e+01,
  769.        
  770.        6.20122878e+03, 2.03132197e-03])
  771.  
  772.  
  773. # Movello di fitting
  774. model = CompositeModel(
  775.     Model(gauss,
  776.           Parameter(value=startPars[0], bounds=(1e5, 3e6),      label='A'),
  777.           Parameter(value=startPars[1],   bounds=(0, 4095),   label='mu'),
  778.           Parameter(value=startPars[2],   bounds=(0, 100), label='sigma'),
  779.           label='gaussiana 1'),
  780.     Model(gauss,
  781.           Parameter(value=startPars[3], bounds=(1e5, 3e6),      label='A'),
  782.           Parameter(value=startPars[4],   bounds=(0, 4095),   label='mu'),
  783.           Parameter(value=startPars[5],   bounds=(0, 100), label='sigma'),
  784.           label='gaussiana 2'),
  785.     Model(gauss,
  786.           Parameter(value=startPars[6], bounds=(1e5, 3e6),      label='A'),
  787.           Parameter(value=startPars[7],   bounds=(0, 4095),   label='mu'),
  788.           Parameter(value=startPars[8],   bounds=(0, 100), label='sigma'),
  789.           label='gaussiana 3'),
  790.    
  791.     Model(gauss,
  792.           Parameter(value=startPars[9], bounds=(1e5, 3e6),      label='A'),
  793.           Parameter(value=startPars[10],   bounds=(0, 4095),   label='mu'),
  794.           Parameter(value=startPars[11],   bounds=(0, 100), label='sigma'),
  795.           label='gaussiana 4'),
  796.     Model(gauss,
  797.           Parameter(value=startPars[12], bounds=(1e5, 3e6),      label='A'),
  798.           Parameter(value=startPars[13],   bounds=(0, 4095),   label='mu'),
  799.           Parameter(value=startPars[14],   bounds=(0, 100), label='sigma'),
  800.           label='gaussiana 5'),
  801.     Model(gauss,
  802.           Parameter(value=startPars[15], bounds=(1e5, 3e6),      label='A'),
  803.           Parameter(value=startPars[16],   bounds=(0, 4095),   label='mu'),
  804.           Parameter(value=startPars[17],   bounds=(0, 100), label='sigma'),
  805.           label='gaussiana 6'),
  806.    
  807.     Model(gauss,
  808.           Parameter(value=startPars[18], bounds=(1e4, 3e6),      label='A'),
  809.           Parameter(value=startPars[19],   bounds=(0, 4095),   label='mu'),
  810.           Parameter(value=startPars[20],   bounds=(0, 100), label='sigma'),
  811.           label='gaussiana 7'),
  812.     Model(gauss,
  813.           Parameter(value=startPars[21], bounds=(1e4, 3e6),      label='A'),
  814.           Parameter(value=startPars[22],   bounds=(0, 4095),   label='mu'),
  815.           Parameter(value=startPars[23],   bounds=(0, 100), label='sigma'),
  816.           label='gaussiana 8'),
  817.     Model(gauss,
  818.           Parameter(value=startPars[24], bounds=(1e4, 3e6),      label='A'),
  819.           Parameter(value=startPars[25],   bounds=(0, 4095),   label='mu'),
  820.           Parameter(value=startPars[26],   bounds=(0, 100), label='sigma'),
  821.           label='gaussiana 9'),
  822.  
  823. #    Model(line,
  824. #          Parameter(value=-4.55851, bounds=(-10, 10),  label='m'),
  825. #          Parameter(value=5881.62, bounds=(-100, 10000), label='q'),
  826. #          label='retta'),
  827.     Model(exp,
  828.           Parameter(value=startPars[27], bounds=(0, 1e5),  label='a'),
  829.           Parameter(value=startPars[28], bounds=(0, 1e-2), label='k'),
  830.           label='exp'),
  831.     label='gaussian_peaks')
  832.  
  833.  
  834. # Disegno la curva da muovere
  835. xsample = np.linspace(0, 4095, 1500)
  836. model_graph, = ax.plot(xsample, model(xsample), color='steelblue', \
  837.                        label='Curva fittata')
  838. ax.legend();
  839.  
  840.  
  841. # Mostro la figura
  842. fig
  843.  
  844.  
  845. # Ridimensiono gli assi per fare spazio ai cursori
  846. ax.set_position([0.10, 0.30, 0.84, 0.62])
  847.  
  848. # Chiamo ciò che fa la magia
  849. gui = AutoGUI(model, [model_graph], bbox=[0.20, 0.07, 0.75, 0.16], figure=fig,
  850.               slider_options={'color': 'steelblue'}, data=(xVect, yVect),
  851.               sigma=np.sqrt(yVect));
  852.  
  853. #%%
  854. # Eseguire a mano dopo aver premuto il pulsantone ed ispezionare a mano
  855. # è facilmente castabile in un array di numpy
  856. Riassunto_55 = model.summary()
  857. Riassunto55 = np.array(Riassunto_55)
  858.  
  859.  
  860.  
  861.  
  862.  
  863.  
  864. #%% 56 - Fatto OK
  865. # Nota, ho alzato il bin150->200, ora converge
  866.  
  867. # Carico i dati
  868. numRun = 56
  869. dati = listaMatrici[myDict[numRun]][:,1]
  870. dati = dati[dati<4095]
  871.  
  872.  
  873. # Costruisco l'istogramma e trovo il centro
  874. h, bins = np.histogram(dati, bins=200)
  875. binc = bins[:-1] + (bins[1] - bins[0])
  876.    
  877.    
  878. # Creo figura e assi su cui disegnare
  879. fig = plt.figure()
  880. fig.set_size_inches(12,5)
  881. fig.suptitle("226Ra - 800V - Gain: 100 ", fontsize=16)
  882.  
  883. ax = fig.add_axes([0.15, 0.14, 0.84, 0.70])
  884. ax.set_yscale("log")
  885. ax.grid(True)
  886. ax.set_xlabel("ADC", fontsize=14)
  887. ax.set_ylabel("Entries", fontsize=14)
  888.  
  889.  
  890.  
  891. # Plotto tutti i dati
  892. ax.plot(binc, h, color='black', lw=1, label='Dati completi', ds="steps-mid")
  893.  
  894.  
  895. # Seleziono e Plotto i dati da fittare
  896. cond = (binc<2243) & (binc>213)
  897. xVect = binc[cond]
  898. yVect = h[cond]
  899.  
  900. data_graph, = ax.plot(xVect, yVect, lw=1, label='Dati per fit',\
  901.                       ds="steps-mid", c="tab:green")
  902.  
  903.  
  904.  
  905. # Starting points: inseriti copiando quelli di convergenza
  906. startPars = np.array([
  907.        3.41057648e+05, 2.74009366e+02, 1.21576501e+01,
  908.        9.65121696e+04, 4.77133343e+02, 1.87943137e+01,
  909.        1.23293444e+05, 5.84143244e+02, 1.98466608e+01,
  910.        
  911.        2.36748059e+05, 6.82714027e+02, 2.20560395e+01,
  912.        3.80932808e+05, 7.88170951e+02, 2.36810707e+01,
  913.        2.29361831e+05, 1.25312868e+03, 3.60282529e+01,
  914.        
  915.        9.47249003e+04, 1.53170346e+03, 1.09290881e+02,
  916.        8.10669238e+04, 1.83305278e+03, 1.16567200e+02,
  917.        8.98130756e+04, 2.13912681e+03, 7.72823389e+01,
  918.        
  919.        8.22077398e+03, 2.52499607e-03])
  920.  
  921.  
  922. # Movello di fitting
  923. model = CompositeModel(
  924.     Model(gauss,
  925.           Parameter(value=startPars[0], bounds=(1e5, 3e6),      label='A'),
  926.           Parameter(value=startPars[1],   bounds=(0, 4095),   label='mu'),
  927.           Parameter(value=startPars[2],   bounds=(0, 100), label='sigma'),
  928.           label='gaussiana 1'),
  929.     Model(gauss,
  930.           Parameter(value=startPars[3], bounds=(1e5, 3e6),      label='A'),
  931.           Parameter(value=startPars[4],   bounds=(0, 4095),   label='mu'),
  932.           Parameter(value=startPars[5],   bounds=(0, 100), label='sigma'),
  933.           label='gaussiana 2'),
  934.     Model(gauss,
  935.           Parameter(value=startPars[6], bounds=(1e5, 3e6),      label='A'),
  936.           Parameter(value=startPars[7],   bounds=(0, 4095),   label='mu'),
  937.           Parameter(value=startPars[8],   bounds=(0, 100), label='sigma'),
  938.           label='gaussiana 3'),
  939.    
  940.     Model(gauss,
  941.           Parameter(value=startPars[9], bounds=(1e5, 3e6),      label='A'),
  942.           Parameter(value=startPars[10],   bounds=(0, 4095),   label='mu'),
  943.           Parameter(value=startPars[11],   bounds=(0, 100), label='sigma'),
  944.           label='gaussiana 4'),
  945.     Model(gauss,
  946.           Parameter(value=startPars[12], bounds=(1e5, 3e6),      label='A'),
  947.           Parameter(value=startPars[13],   bounds=(0, 4095),   label='mu'),
  948.           Parameter(value=startPars[14],   bounds=(0, 100), label='sigma'),
  949.           label='gaussiana 5'),
  950.     Model(gauss,
  951.           Parameter(value=startPars[15], bounds=(1e5, 3e6),      label='A'),
  952.           Parameter(value=startPars[16],   bounds=(0, 4095),   label='mu'),
  953.           Parameter(value=startPars[17],   bounds=(0, 100), label='sigma'),
  954.           label='gaussiana 6'),
  955.    
  956.     Model(gauss,
  957.           Parameter(value=startPars[18], bounds=(1e4, 3e6),      label='A'),
  958.           Parameter(value=startPars[19],   bounds=(0, 4095),   label='mu'),
  959.           Parameter(value=startPars[20],   bounds=(0, 100), label='sigma'),
  960.           label='gaussiana 7'),
  961.     Model(gauss,
  962.           Parameter(value=startPars[21], bounds=(1e4, 3e6),      label='A'),
  963.           Parameter(value=startPars[22],   bounds=(0, 4095),   label='mu'),
  964.           Parameter(value=startPars[23],   bounds=(0, 100), label='sigma'),
  965.           label='gaussiana 8'),
  966.     Model(gauss,
  967.           Parameter(value=startPars[24], bounds=(1e4, 3e6),      label='A'),
  968.           Parameter(value=startPars[25],   bounds=(0, 4095),   label='mu'),
  969.           Parameter(value=startPars[26],   bounds=(0, 100), label='sigma'),
  970.           label='gaussiana 9'),
  971.  
  972. #    Model(line,
  973. #          Parameter(value=-4.55851, bounds=(-10, 10),  label='m'),
  974. #          Parameter(value=5881.62, bounds=(-100, 10000), label='q'),
  975. #          label='retta'),
  976.     Model(exp,
  977.           Parameter(value=startPars[27], bounds=(0, 1e5),  label='a'),
  978.           Parameter(value=startPars[28], bounds=(0, 1e-2), label='k'),
  979.           label='exp'),
  980.     label='gaussian_peaks')
  981.  
  982.  
  983. # Disegno la curva da muovere
  984. xsample = np.linspace(0, 4095, 1500)
  985. model_graph, = ax.plot(xsample, model(xsample), color='steelblue', \
  986.                        label='Curva fittata')
  987. ax.legend();
  988.  
  989.  
  990. # Mostro la figura
  991. fig
  992.  
  993.  
  994. # Ridimensiono gli assi per fare spazio ai cursori
  995. ax.set_position([0.10, 0.30, 0.84, 0.62])
  996.  
  997. # Chiamo ciò che fa la magia
  998. gui = AutoGUI(model, [model_graph], bbox=[0.20, 0.07, 0.75, 0.16], figure=fig,
  999.               slider_options={'color': 'steelblue'}, data=(xVect, yVect),
  1000.               sigma=np.sqrt(yVect));
  1001.  
  1002. #%%
  1003. # Eseguire a mano dopo aver premuto il pulsantone ed ispezionare a mano
  1004. # è facilmente castabile in un array di numpy
  1005. Riassunto_56 = model.summary()
  1006. Riassunto56 = np.array(Riassunto_56)
  1007.  
  1008.  
  1009.  
  1010. #%% Esporto i dati
  1011. np.savez("Dati", rias52 = Riassunto52, rias53 = Riassunto53, \
  1012.                  rias54 = Riassunto54, rias55 = Riassunto55, \
  1013.                  rias56 = Riassunto56, )
  1014.  
  1015.  
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement