SHARE
TWEET

Untitled

a guest Aug 19th, 2019 67 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. from scipy.optimize import curve_fit
  4.  
  5.  
  6. datapath = '/Volumes/TOSHIBA EXT/23-07-19 - Test Pulses (22deg)/'
  7. data_path_output = '/Users/marialauraperez/Desktop/SummerCERN/Results/TestPulses-22deg/'
  8. first_charges = np.arange(500, 1600, 100)
  9. second_charges = np.arange(2000, 12000, 1000)
  10. charges = np.concatenate((first_charges, second_charges))
  11. ToT_data = list()
  12. mean_ToT = list()
  13. energies = charges*3.62/1000
  14.  
  15. for i in range(len(charges)):
  16.     matrix = np.genfromtxt(datapath+str(charges[i])+'e.txt_iToT.txt')
  17.     ToT_data.append(matrix)
  18.     mean_ToT.append(np.mean(matrix)/100)
  19.  
  20. # isotope_energies = [5.89, 22.16, 59.54]
  21. # isotope_ToTs = [14.21, 48.76, 118.28]
  22. theoretical_XRF = [4.512, 5.9, 6.405, 8.046, 15.775, 17.48, 25.271, 59.5]
  23. ToT_XRF = [12.05, 15.82, 16.91, 20.98, 38.83, 42.64, 58.97, 118.54]
  24.  
  25.  
  26. def surrogate_function(x, a, b, c, t):
  27.     y = a*x + b - c/(x-t)
  28.     y[np.where(x < t)] = 0
  29.     y *= y > 0
  30.     return y
  31.  
  32.  
  33. poptmean, pcovn = curve_fit(surrogate_function, energies, mean_ToT)
  34.  
  35. x = np.linspace(1.5, 65, 100)
  36.  
  37. plt.figure()
  38. plt.scatter(energies, mean_ToT, s=10, label='TP Data', color='navy')
  39. plt.plot(x, surrogate_function(x, *poptmean), label='TP Fit: a = '+str(np.round(poptmean[0], 3))
  40.                                                     + '; b = ' + str(np.round(poptmean[1], 3))
  41.                                                     + '; c = '+str(np.round(poptmean[2], 3))
  42.                                                     + '; t = '+str(np.round(poptmean[3], 3)), color='blue')
  43. plt.scatter(theoretical_XRF, ToT_XRF, s=10, label='Source Data', color='red')
  44. plt.xlabel('Energy [keV]')
  45. plt.ylabel('ToT value')
  46. plt.title('Source and test pulse calibrations (whole detector)', fontweight='bold')
  47. plt.legend()
  48. plt.savefig('/Users/marialauraperez/Desktop/total_per_source/Source+TestPulse_Calibs22deg.png')
  49. plt.show()
  50.  
  51. # plt.figure()
  52. # plt.scatter(energies, mean_ToT, s=10, label='Data', color='navy')
  53. # plt.plot(x, surrogate_function(x, *poptmean), label='Fit: a = '+str(np.round(poptmean[0], 3))
  54. #                                                     + '; b = ' + str(np.round(poptmean[1], 3))
  55. #                                                     + '; c = '+str(np.round(poptmean[2], 3))
  56. #                                                     + '; t = '+str(np.round(poptmean[3], 3)), color='blue')
  57. # plt.xlabel('Test pulse height [keV]')
  58. # plt.ylabel('Mean ToT value')
  59. # plt.title('Calibration with test pulses', fontweight='bold')
  60. # plt.legend()
  61. # plt.savefig(data_path_output+'22degTestPulse_AvgCalib.png')
  62.  
  63.  
  64. def get_calibration_matrices(energy_vals, ToT_matrices):
  65.     nrows = len(ToT_matrices[0][:, 0])
  66.     ncols = len(ToT_matrices[0][0, :])
  67.     a = np.zeros(np.shape(ToT_matrices[0]))
  68.     b = a.copy()
  69.     c = a.copy()
  70.     t = a.copy()
  71.     test_vals_y = list()
  72.     pixel_i = list()
  73.     pixel_j = list()
  74.     test_popt = list()
  75.     for i in range(nrows):
  76.         for j in range(ncols):
  77.             entry = list()
  78.             for k in range(len(ToT_matrices)):
  79.                 entry.append(ToT_matrices[k][i, j]/100)
  80.             try:
  81.                 popt, pcov = curve_fit(surrogate_function, energy_vals, entry, p0=poptmean)
  82.                 a[i, j] = popt[0]
  83.                 b[i, j] = popt[1]
  84.                 c[i, j] = popt[2]
  85.                 t[i, j] = popt[3]
  86.                 if (i is 20 or i is 148 or i is 50) and (j is 130 or j is 210 or j is 90):
  87.                     test_vals_y.append(entry)
  88.                     pixel_i.append(i)
  89.                     pixel_j.append(j)
  90.                     test_popt.append(popt)
  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.  
  97.     return a, b, c, t, pixel_i, pixel_j, test_vals_y, test_popt
  98.  
  99.  
  100. # a_matrix, b_matrix, c_matrix, t_matrix, pixels_x, pixels_y, test_vals, test_params = \
  101. #     get_calibration_matrices(energies, ToT_data)
  102. #
  103. # plt.figure()
  104. # for k in range(len(pixels_x)):
  105. #     plt.scatter(energies, test_vals[k], s=10, label='('+str(pixels_x[k])+', '+str(pixels_y[k])+')')
  106. #     plt.plot(x, surrogate_function(x, *test_params[k]), label='Fit')
  107. # plt.xlabel('Test pulse height [keV]')
  108. # plt.ylabel('ToT value')
  109. # plt.title('Calibration per pixel with test pulses', fontweight='bold')
  110. # plt.legend(fontsize='small', ncol=2)
  111. # plt.savefig(data_path_output+'22degTestPulse_CalibsPixels.png')
  112. #
  113. # np.savetxt(data_path_output+'A_TestPulse22deg.txt', a_matrix, fmt='%.5f')
  114. # np.savetxt(data_path_output+'B_TestPulse22deg.txt', b_matrix, fmt='%.5f')
  115. # np.savetxt(data_path_output+'C_TestPulse22deg.txt', c_matrix, fmt='%.5f')
  116. # np.savetxt(data_path_output+'T_TestPulse22deg.txt', t_matrix, fmt='%.5f')
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
 
Top