SHARE
TWEET

Untitled

a guest Aug 19th, 2019 70 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
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