# Untitled

a guest
Aug 19th, 2019
80
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. from scipy.optimize import curve_fit
2. import numpy as np
3. import matplotlib.pyplot as plt
4.
5. # ---------- DEFINE ---------------
6. datapath = '/Volumes/TOSHIBA EXT/Test-Pulses4(22deg)/'
7. datapath2 = '/Users/marialauraperez/Desktop/SummerCERN/ChunkedSources/'
8. data_path_output = '/Users/marialauraperez/Desktop/'
9. txt_output = '/Users/marialauraperez/Desktop/'
10. first_charges = np.arange(500, 1600, 100)
11. second_charges = np.arange(2000, 12000, 1000)
12. charges = np.concatenate((first_charges, second_charges))
13. ToT_TP = list()
14. ToT_sources = list()
15. energies_TP = charges*3.62/1000
16. del first_charges, second_charges
17. energies_sources = [4.512, 5.9, 6.405, 8.046, 15.775, 17.48, 25.271]#, 59.5]
18.
19. # --------- IMPORT --------------
20. for q in range(len(charges)):
21. matrix = np.genfromtxt(datapath+str(charges[q])+'e.txt_iToT.txt')
22. ToT_TP.append(matrix)
23. for source in ['Ti', 'Mn', 'Fe', 'Cu', 'Zr', 'Mo', 'Sn']:#, 'Am241']:
24. matrix = np.genfromtxt(datapath2+source+'_FinalPeakMatrix.txt').reshape([256, 256])
25. ToT_sources.append(matrix)
26.
27. # -------- FUNCTIONS -----------
28.
29.
30. def surrogate_function_1(x, a, b, c, t, g):
31. y = a*g*x + b - c/(g*x-t)
32. y[np.where(g*x < t)] = 0
33. y *= y > 0
34. return y
35.
36.
37. def surrogate_function_2(A, B, C, T):
38. def surrogate_function_g(x, g):
39. y = g*A*x + B - C/(g*x-T)
40. y[np.where(g * x < T)] = 0
41. y *= y > 0
42. return y
43. return surrogate_function_g
44.
45.
46. def normal_surrogate_function(x, a, b, c, t):
47. y = a*x + b - c/(x-t)
48. y[np.where(x < t)] = 0
49. y *= y > 0
50. return y
51.
52.
53. def get_calibration_matrices(energies_TP, ToT_mats_TP, energies_sources, ToT_mats_sources):
54. energies = [energies_TP, energies_sources]
55. ToTs = [ToT_mats_TP, ToT_mats_sources]
56. a_mat = np.zeros(np.shape(ToT_mats_TP[0]))
57. b_mat = np.zeros(np.shape(ToT_mats_TP[0]))
58. c_mat = np.zeros(np.shape(ToT_mats_TP[0]))
59. t_mat = np.zeros(np.shape(ToT_mats_TP[0]))
60. g_mats = [np.zeros(np.shape(ToT_mats_TP[0])), np.zeros(np.shape(ToT_mats_TP[0]))]
61. px125 = list()
62. for w in range(2):
63. ToT_matrices = ToTs[w]
64. energy_vals = energies[w]
65. g = g_mats[w]
66. for i in range(256):
67. for j in range(256):
68. entry = list()
69. for k in range(len(ToT_matrices)):
70. if w is 0:
71. entry.append(ToT_matrices[k][i, j]/100)
72. else:
73. entry.append(ToT_matrices[k][i, j])
74. if i is 125 and j is 125:
75. px125.append(entry)
76. try:
77. if w is 0:
78. popt, pcov = curve_fit(surrogate_function_1, energy_vals, entry)
79. a_mat[i, j] = popt[0]
80. b_mat[i, j] = popt[1]
81. c_mat[i, j] = popt[2]
82. t_mat[i, j] = popt[3]
83. g[i, j] = popt[4]
84. else:
85. try:
86. popt, pcov = curve_fit(surrogate_function_2(a_mat[i, j], b_mat[i, j], c_mat[i, j],
87. t_mat[i, j]), energy_vals, entry)
88. g[i, j] = popt[0]
89. except ValueError:
90. print('Error: calibration not performed in pixel (', str(i)+', '+str(j)+' )')
91. except RuntimeError:
92. print('Error: calibration not performed in pixel (', str(i)+', '+str(j)+' )')
93. # plt.figure()
94. # plt.scatter(energy_vals, entry)
95. # plt.show()
96. return a_mat, b_mat, c_mat, t_mat, g_mats, px125
97.
98. # ------------- EXCECUTION ------------------
99.
100.
101. a, b, c, t, [g1, g2], [px1, px2] = get_calibration_matrices(energies_TP, ToT_TP, energies_sources,
102. ToT_sources)
103.
104. new_a = g2*a
105. new_b = b
106. new_c = c/g2
107. new_t = t/g2
108.
109. x1 = np.linspace(0, 70, 100)
110. plt.figure()
111. plt.scatter(energies_TP, px1, label='Test pulse points')
112. plt.plot(x1, surrogate_function_1(x1, a[125, 125], b[125, 125], c[125, 125], t[125, 125], g1[125, 125]),
113. label='Test pulse fit: g = '+str(np.round(g1[125, 125], 3)))
114. plt.scatter(energies_sources, px2, label='Source points')
115. plt.plot(x1, surrogate_function_1(x1, a[125, 125], b[125, 125], c[125, 125], t[125, 125], g2[125, 125]),
116. label='Source fit: g = '+str(np.round(g2[125, 125], 3)))
117. plt.legend()
118. plt.title('Energy calibration curve for pixel (125, 125)', fontweight='bold')
119. plt.xlabel('Energy [keV]')
120. plt.ylabel('ToT value')
121. plt.savefig(data_path_output+'Pixel125_CalibGain_NoAm.png')
122. # plt.show()
123.
124. np.savetxt(txt_output+'A_STP_NoAm.txt', new_a, fmt='%.5f')
125. np.savetxt(txt_output+'B_STP_NoAm.txt', new_b, fmt='%.5f')
126. np.savetxt(txt_output+'C_STP_NoAm.txt', new_c, fmt='%.5f')
127. np.savetxt(txt_output+'T_STP_NoAm.txt', new_t, fmt='%.5f')
128.
129. np.savetxt(txt_output+'G_TestPulse_NoAm.txt', g1, fmt='%.5f')
130. np.savetxt(txt_output+'G_Source_NoAm.txt', g2, fmt='%.5f')
RAW Paste Data