SHARE
TWEET

Untitled

a guest Aug 19th, 2019 60 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.  
  6. def surrogate_function(x, a, b, c, t, g):
  7.     y = a*g*x + b - c/(g*x-t)
  8.     y[np.where(g*x < t)] = 0
  9.     y *= y > 0
  10.     return y
  11.  
  12.  
  13. datapath = '/Volumes/TOSHIBA EXT/23-07-19 - Test Pulses (22deg)/'
  14. data_path_output = '/Users/marialauraperez/Desktop/SummerCERN/Results/TestPulses-22deg/'
  15. first_charges = np.arange(500, 1600, 100)
  16. second_charges = np.arange(2000, 12000, 1000)
  17. charges = np.concatenate((first_charges, second_charges))
  18. ToT_data = list()
  19. mean_ToT = list()
  20. energies = charges*3.62/1000
  21.  
  22. for i in range(len(charges)):
  23.     matrix = np.genfromtxt(datapath+str(charges[i])+'e.txt_iToT.txt')
  24.     ToT_data.append(matrix)
  25.     mean_ToT.append(np.mean(matrix)/100)
  26.  
  27. theoretical_XRF = [4.512, 6.405, 8.046, 15.775, 17.48, 25.271, 59.5]
  28. ToT_XRF = [12.05, 16.91, 20.98, 38.83, 42.64, 58.97, 118.54]
  29.  
  30. poptmean, pcovn = curve_fit(surrogate_function, energies, mean_ToT)
  31.  
  32.  
  33. def surrogate_function_ABCT(A, B, C, T):
  34.     def surrogate_function_g(x, g):
  35.         y = g*A*x + B - C/(g*x-T)
  36.         y[np.where(g * x < T)] = 0
  37.         y *= y > 0
  38.         return y
  39.     return surrogate_function_g
  40.  
  41.  
  42. popt1, pcov1 = curve_fit(surrogate_function_ABCT(poptmean[0], poptmean[1], poptmean[2], poptmean[3]),
  43.                          theoretical_XRF, ToT_XRF)
  44.  
  45. x = np.linspace(2, 60, 100)
  46. plt.figure()
  47. plt.title('Test pulse and source calibrations (whole detector)', fontweight='bold')
  48. plt.xlabel('Energy [keV]')
  49. plt.ylabel('ToT value')
  50. plt.plot(x, surrogate_function(x, *poptmean), label='TP Fit: a = '+str(np.round(poptmean[0], 3))
  51.                                                     + '; b = ' + str(np.round(poptmean[1], 3))
  52.                                                     + '; c = '+str(np.round(poptmean[2], 3))
  53.                                                     + '; t = '+str(np.round(poptmean[3], 3))
  54.                                                     + '; g = '+str(np.round(poptmean[4], 3)))
  55. plt.scatter(energies, mean_ToT, label='TP points', s=10)
  56. plt.plot(x, surrogate_function_ABCT(poptmean[0], poptmean[1], poptmean[2], poptmean[3])(x, popt1[0]),
  57.          label='Source fit: g = '+str(np.round(popt1[0], 3)))
  58. plt.scatter(theoretical_XRF, ToT_XRF, label='Source points', s=10)
  59. plt.text(5.25, 6.94, 'Ti', fontsize=9)
  60. plt.text(7.4, 13.5, 'Fe', fontsize=9)
  61. plt.text(9.4, 19.9, 'Cu', fontsize=9)
  62. plt.text(16.7, 34.7, 'Zr', fontsize=9)
  63. plt.text(18.5, 41.4, 'Mo', fontsize=9)
  64. plt.text(26.16, 57, 'Sn', fontsize=9)
  65. plt.text(56.2, 110.4, 'Am 241', fontsize=9)
  66. plt.legend(fontsize='x-small', loc=2)
  67. plt.savefig('/Users/marialauraperez/Desktop/MeanCalib_Source+TP_22deg.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. OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!
 
Top