-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathQBatchPeakFit.py
69 lines (52 loc) · 2.26 KB
/
QBatchPeakFit.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
# -*- coding: utf-8 -*-
"""
Created on Thu Mar 7 08:14:03 2019
@author: Laura
"""
import sys
import numpy as np
import scipy.optimize as op
from scipy import integrate
from QUtility import *
def fit_peak(specfile, peak):
result = np.zeros(1, dtype=QUtility.peakfit_dtype)
result['name'] = specfile
print('Reading spectrum: {}'.format(specfile))
spectrum = QUtility.read_spec(specfile)
aoi_spec = QUtility.spectrum_slice(spectrum, peak-50, peak+50)
print('Fitting baseline to spectrum.')
bg_spec_l = QUtility.spectrum_slice(spectrum, peak-50, peak-20)
bg_spec_r = QUtility.spectrum_slice(spectrum, peak+20, peak+50)
bg_spec = np.vstack((bg_spec_l, bg_spec_r))
aoi_bg_params = np.polyfit(bg_spec[:,0], bg_spec[:,1], 3) # fit 3rd order polynomial baseline
aoi_bg = np.polyval(aoi_bg_params, aoi_spec[:,0])
(bg_a, bg_b, bg_c, bg_d) = aoi_bg_params
result['bg_a'] = bg_a
result['bg_b'] = bg_b
result['bg_c'] = bg_c
result['bg_d'] = bg_d
absorp_new = aoi_spec[:,1] - aoi_bg
wav_new = aoi_spec[:,0]
spec_new = np.column_stack((wav_new, absorp_new))
print('Fitting peak at {}.'.format(peak))
wav_inter = np.arange(wav_new[0], wav_new[-1], 0.1)
peak_inter = QUtility.inter(spec_new, wav_inter)
first_guess = (peak, 0, 1, 1, 0.5) #first guess for x0, I, HWHM_left and HWHM_right, sigma
bounds = [(peak-2,peak+2),(0,None),(0.001,10),(0.001,10),(0,1)]
res = op.minimize(QUtility.pseudovoigt, method='SLSQP', args=(wav_inter, peak_inter), x0=first_guess, bounds=bounds)
#fit = utility.pseudovoigt_fit(wav_new, *res.x)
(pos, I, HWHM_l, HWHM_r, sigma) = res.x
area_ana = QUtility.peak_area(I, HWHM_l, HWHM_r, sigma)
area_num = integrate.simps(absorp_new)
result['x0'] = pos
result['I'] = I
result['HWHM_l'] = HWHM_l
result['HWHM_r'] = HWHM_r
result['sigma'] = sigma
result['area_ana'] = area_ana
result['area_num'] = area_num
for item in zip((pos, I, HWHM_l, HWHM_r, sigma), ('pos','I', 'HWHM_l','HWHM_r','sigma')):
print(item[1]+': '+str(item[0])+'\n')
return result #(specfile, pos, I, HWHM_l, HWHM_r, sigma, area_ana, area_num, bg_a, bg_b, bg_c, bg_d)
if __name__ == "__main__":
fit_peak(sys.argv[1], sys.argv[2])