• API
• FAQ
• Tools
• Archive
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.
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!

Top