python应用 曲线拟合03

  1 # -*- coding: utf-8 -*-
  2 """
  3 @author: kurrrr
  4 """
  5 
  6 import numpy as np
  7 import matplotlib.pyplot as plt
  8 import struct
  9 from scipy.optimize import curve_fit
 10 
 11 t_pre = 26  # global
 12 p = 0.6  # global
 13 q = 1.2  # global
 14 
 15 
 16 class Signal(object):
 17     _t = [x/250.0 for x in range(0, 20000)]
 18 
 19     def __init__(self):
 20         self._u = []
 21 
 22     @property
 23     def u(self):
 24         return self._u
 25 
 26     @u.setter
 27     def u(self, u):
 28         self._u = u
 29 
 30     @property
 31     def t(self):
 32         return self._t
 33 
 34     def draw(self):
 35         plt.scatter(self._t, self._u, marker=‘.‘, label="original data")
 36         plt.show()
 37 
 38     def get_baseline(self):
 39         return np.mean(self._u[0:250])  # 250 points to get baseline
 40 
 41     def get_max(self):
 42         return np.max(self._u)
 43 
 44     def get_max_t(self):
 45         return np.argmax(self._u)/250.0
 46 
 47 
 48 ‘‘‘
 49 t_pre :     preamplifier tau
 50 b:          number of fast photons
 51 c:          number of slow photons
 52 p:          tau of fast photons
 53 q:          tau of slow photons
 54 v:          x offset
 55 w:          y offset
 56 ‘‘‘
 57 
 58 
 59 def pre_func_1(x, w):
 60     return w
 61 
 62 
 63 def pre_func_2(x, b, c, p, q, v, w):
 64     global t_pre
 65     return b*t_pre/(t_pre-p) *  66             (np.exp(-(x-v)/t_pre)-np.exp(-(x-v)/p))  67             + c*t_pre/(t_pre-q) *  68             (np.exp(-(x-v)/t_pre)-np.exp(-(x-v)/q)) + w
 69 
 70 
 71 def pre_func(x, b, c, p, q, v, w):
 72     return np.piecewise(x, [x < v, x >= v], [lambda x: pre_func_1(x, w),
 73                         lambda x: pre_func_2(x, b, c, p, q, v, w)])
 74 
 75 
 76 def main():
 77     global t_pre
 78     global p
 79     global q
 80 
 81     data_file = open(‘run0035.bin‘, ‘rb‘)
 82     signal = Signal()
 83 
 84     fast = []
 85     slow = []
 86     u_temp = []
 87     for n in range(10000):
 88         u_temp.clear()
 89         for _ in range(0, 20000):
 90             data_temp = data_file.read(2)
 91             data_dec, = struct.unpack(‘h‘, data_temp)
 92             u_temp.append(-data_dec)
 93         signal.u = u_temp
 94 
 95         u_baseline = signal.get_baseline()
 96         u_max = signal.get_max()
 97         t_max = signal.get_max_t()
 98         b = (u_baseline - u_max) / 2.0
 99         c = (u_baseline - u_max) / 2.0
100         v = t_max
101         w = u_max
102 
103         fit_par = [b, c, p, q, v, w]
104         popt, pcov = curve_fit(pre_func, signal.t,
105                                signal.u, fit_par, maxfev=50000)
106 
107         if popt[2] < popt[3]:
108             fast.append(popt[2])
109             slow.append(popt[3])
110         else:
111             fast.append(popt[3])
112             slow.append(popt[2])
113 
114     fast_slow = list(zip(fast, slow))
115     for xx in fast_slow:
116         if xx[0] > 100 or xx[1] >100:
117             fast_slow.remove(xx)
118 
119     fast_new = []
120     slow_new = []
121     for xx in fast_slow:
122         fast_new.append(xx[0])
123         slow_new.append(xx[1])
124     plt.scatter(fast_new, slow_new, marker=‘.‘)
125     plt.show()
126 
127     data_file.close()
128 
129 
130 if __name__ == ‘__main__‘:
131     main()