SHARE
TWEET

Untitled

a guest Aug 19th, 2019 63 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. import numpy as np
  2. import matplotlib.pyplot as plt
  3. import Fit
  4.  
  5. data_path_disk = '/Users/marialauraperez/Desktop/SummerCERN/npz histos 22 deg/'
  6. output_data_path = '/Users/marialauraperez/Desktop/SummerCERN/ChunkedSources/'
  7. sources = ['Ti', 'Mn', 'Fe', 'Cu', 'Zr', 'Mo', 'Sn']#, 'Am241']
  8. colors = ['red', 'orange', 'yellowgreen', 'green', 'dodgerblue', 'navy', 'purple', 'pink']
  9.  
  10. # for i in range(4, 7):
  11. #     print('processing data from ', sources[i])
  12. #     matrix_path = output_data_path + ('%s_FinalPeakMatrix.txt' % sources[i])
  13. #     data = np.load(data_path_disk+sources[i]+'XRF_cluster_sp_hist2.npz')
  14. #     data_hist = data['hist']  # 65536 x 200 array. make a histogram for every index
  15. #     data_tot = data['ToT']
  16. #     p0 = None
  17. #     f = Fit.FitGauss()
  18. #     peaks = np.zeros(len(data_hist))
  19. #     for j in range(len(data_hist)):
  20. #         try:
  21. #             gauss = f.fit(x=data_tot, y=data_hist[j], p0=p0)
  22. #             p0 = gauss[0]
  23. #             mu = gauss[0][1]
  24. #             if gauss[2].R2 < 0.8:
  25. #                 p0 = None
  26. #         except:
  27. #             mu = 0.0
  28. #         peaks[j] = mu
  29. #         if j%1000 is 0:
  30. #             print(j)
  31. #     np.savetxt(matrix_path, peaks, fmt='%.4f', delimiter='\n')
  32.  
  33. # ----------- FOR AMERICIUM, GADOLINIUM AND TIN ---------------------------------
  34. # print('processing data from ', sources[7])
  35. # matrix_path = output_data_path + ('%s_FinalPeakMatrix.txt' % sources[7])
  36. # data = np.load(data_path_disk+sources[7]+'XRF_cluster_sp_hist2.npz')
  37. # data_hist = data['hist']  # 65536 x 200 array. make a histogram for every index
  38. # data_tot = data['ToT']
  39. # p0 = None
  40. # p2 = None
  41. # f = Fit.FitGauss()
  42. # peaks_am = np.zeros(len(data_hist))
  43. # peaks_sn = np.zeros(len(data_hist))
  44. # for j in range(len(data_hist)):
  45. #     try:
  46. #         argmax_am = np.argmax(data_hist[j][100:200])+100
  47. #         argmax_sn = np.argmax(data_hist[j][40:90])+40
  48. #
  49. #         gauss_am = f.fit(x=data_tot[argmax_am-20:argmax_am+20], y=data_hist[j][argmax_am-20:argmax_am+20], p0=p0)
  50. #         gauss_sn = f.fit(x=data_tot[argmax_sn-20:argmax_sn+20], y=data_hist[j][argmax_sn-20:argmax_sn+20], p0=p2)
  51. #         mu_am = gauss_am[0][1]
  52. #         mu_sn = gauss_sn[0][1]
  53. #         p0 = gauss_am[0]
  54. #         p2 = gauss_sn[0]
  55. #         if gauss_am[2].R2 < 0.8:
  56. #             p0 = None
  57. #         if gauss_sn[2].R2 < 0.8:
  58. #             p2 = None
  59. #     except:
  60. #         mu_am = 0.0
  61. #         mu_sn = 0.0
  62. #     peaks_am[j] = mu_am
  63. #     peaks_sn[j] = mu_sn
  64. #     if j % 1000 is 0:
  65. #         print(j)
  66. # np.savetxt(output_data_path+('%s_FinalPeakMatrix.txt' % sources[7]), peaks_am, fmt='%.4f', delimiter='\n')
  67. # np.savetxt(output_data_path+'(Check)Sn_FinalPeakMatrix.txt', peaks_sn, fmt='%.4f', delimiter='\n')
  68.  
  69.  
  70. def normalize(v):
  71.     norm = np.linalg.norm(v)
  72.     if norm == 0:
  73.         return v
  74.     return v / norm
  75.  
  76.  
  77. mus = np.zeros(len(sources))
  78. plt.figure(figsize=(9,5))
  79. plt.xlabel('ToT value')
  80. plt.ylabel('Normalized counts')
  81. plt.title('Total spectra with Timepix3 (T=22ºC)', fontweight='bold')
  82. #plt.ticklabel_format(style='sci', axis='y', scilimits=(0, 5))
  83. plt.xlim(0, 90)
  84. for i in range(len(sources)):
  85.     data = np.load(data_path_disk + sources[i] + 'XRF_cluster_sp_hist2.npz')
  86.     data_hist = data['hist']  # 65536 x 200 array. make a histogram for every index
  87.     data_tot = data['ToT']
  88.     mean = np.zeros(np.shape(data_tot))
  89.     for j in range(len(data_tot)):
  90.         mean[j] = np.sum(np.transpose(data_hist)[:][j])
  91.     mean = normalize(mean)
  92.     f = Fit.FitGauss()
  93.     if sources[i] != 'Am241':
  94.         gauss = f.fit(x=data_tot, y=mean)
  95.         y1 = Fit.gauss(data_tot, gauss[0])
  96.         p0 = gauss[0]
  97.         mu = gauss[0][1]
  98.         mus[i] = mu
  99.         plt.scatter(data_tot, mean, label='_nolegend_', s=10, color=colors[i])
  100.         plt.plot(data_tot, y1, label=sources[i]+r': $\mu=$ ' + str(np.round(mu, 3)), linestyle='--', color=colors[i])
  101.     else:
  102.         gauss_am = f.fit(x=data_tot[100:150], y=mean[100:150])
  103.         gauss_gd = f.fit(x=data_tot[75:100], y=mean[75:100])
  104.         gauss_sn = f.fit(x=data_tot[40:75], y=mean[40:75])
  105.         y1 = Fit.gauss(data_tot, gauss_am[0])
  106.         y2 = Fit.gauss(data_tot, gauss_gd[0])
  107.         y3 = Fit.gauss(data_tot, gauss_sn[0])
  108.         mu_am = gauss_am[0][1]
  109.         mu_gd = gauss_gd[0][1]
  110.         mu_sn = gauss_sn[0][1]
  111.         mus[i] = mu_am
  112.         plt.scatter(data_tot, mean, label='_nolegend_', s=10, color=colors[i])
  113.         plt.plot(data_tot[80:150], y1[80:150], label=r'Am: $\mu=$ ' + str(np.round(mu_am, 3)), linestyle='--', color='hotpink')
  114.         #plt.plot(data_tot, y2, label='_nolegend_', linestyle='--')
  115.         #plt.plot(data_tot, y3, label='_nolegend_', linestyle='--')
  116. plt.legend()
  117. plt.savefig('/Users/marialauraperez/Desktop/totalfitssources_noam.png')
  118.  
  119. # -------- FOR CALCIUM AND COPPER --------------
  120.  
  121. # data = np.load(data_path_disk + 'CaXRF_cluster_sp_hist2.npz')
  122. # data_hist = data['hist']  # 65536 x 200 array. make a histogram for every index
  123. # data_tot = data['ToT']
  124. # mean = np.zeros(np.shape(data_tot))
  125. # for j in range(len(data_tot)):
  126. #     mean[j] = np.sum(np.transpose(data_hist)[:][j])
  127. # f = Fit.FitGauss()
  128. # gauss_ca = f.fit(x=data_tot[0:13], y=mean[0:13])
  129. # gauss_cu = f.fit(x=data_tot[13:25], y=mean[13:25])
  130. # y1 = Fit.gauss(data_tot, gauss_ca[0])
  131. # y2 = Fit.gauss(data_tot, gauss_cu[0])
  132. # mu_ca = gauss_ca[0][1]
  133. # mu_cu = gauss_cu[0][1]
  134. # plt.figure()
  135. # plt.plot(data_tot, mean, label='Data')
  136. # plt.plot(data_tot, y1, label=r'Fit Ca: $\mu=$ ' + str(np.round(mu_ca, 3)), linestyle='--')
  137. # plt.plot(data_tot, y2, label=r'Fit Cu: $\mu=$ ' + str(np.round(mu_cu, 3)), linestyle='--')
  138. # plt.xlabel('ToT value')
  139. # plt.ylabel('Total counts')
  140. # plt.title('Ca total spectrum with Timepix3 (T=22ºC)', fontweight='bold')
  141. # plt.ticklabel_format(style='sci', axis='y', scilimits=(0, 5))
  142. # plt.legend()
  143. # plt.savefig('/Users/marialauraperez/Desktop/' + 'Ca_totalfit.png')
RAW Paste Data
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!
 
Top