Guest User

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