zebra/pyzebra
3
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

peak fitting

Browse Source 
Fit one peak from ccl as combination of gaussian, lorentian and backgroud
This commit is contained in:
JakHolzer 2020-09-10 15:27:53 +02:00 committed by Ivan Usov
parent 208c165922
commit 087fe0fe2f
1 changed files with 151 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

151
pyzebra/fit2.py Normal file
View File

@ -0,0 +1,151 @@
from lmfit import minimize, Parameters, Model
from lmfit.models import LinearModel, LorentzianModel, GaussianModel
import matplotlib.pyplot as plt
from scipy.integrate import simps
import scipy as sc
from scipy import integrate
import numpy as np
from time import sleep
def fitccl(data, keys, guess, vary, constraints_min, constraints_max, numfit_min=None, numfit_max=None):
"""Made for fitting of ccl date where 1 peak is expected. Allows for combination of gaussian, lorentzian and linear model combination
:param data: dictionary after peak fining
:param keys: name of the measurement in the data dict (i.e. M123)
:param guess: initial guess for the fitting, if none, some values are added automatically in order (see below)
:param vary: True if parameter can vary during fitting, False if it to be fixed
:param numfit_min: minimal value on x axis for numerical integration - if none is centre of gaussian minus 3 sigma
:param numfit_max: maximal value on x axis for numerical integration - if none is centre of gaussian plus 3 sigma
:param constraints_min: min constranits value for fit
:param constraints_max: max constranits value for fit
:return data dict with additional values
order for guess, vary, constraints_min, constraints_max
[Gaussian centre, Gaussian sigma, Gaussian amplitude, Lorentzian centre, Lorentzian sigma, Lorentzian amplitude, background slope, background intercept]
examples:
guess = [None, None, 100, None, None, None, 0, None]
vary = [True, True, True, True, False, True, True, True]
constraints_min = [23, None, 50, None, None, None, 0, 0]
constraints_min = [80, None, 1000, None, None, None, 0, 100]
"""
if len(data["Measurements"][str(keys)]["peak_indexes"]) == 1:
x = list(data["Measurements"][str(keys)]["omega"])
y = list(data["Measurements"][str(keys)]["counts"])
peak_index = data["Measurements"][str(keys)]["peak_indexes"]
peak_height = data["Measurements"][str(keys)]["peak_heights"]
print('before', constraints_min)
guess[0] = x[int(peak_index)] if guess[0] is None else guess[0]
guess[1] = 0.1 if guess[1] is None else guess[1]
guess[2] = float(peak_height/10) if guess[2]is None else float(guess[2])
guess[3] = x[int(peak_index)] if guess[3] is None else guess[3]
guess[4] = 2*guess[1] if guess[4] is None else guess[4]
guess[5] = float(peak_height/10) if guess[5] is None else float(guess[5])
guess[6] = 0 if guess[6] is None else guess[6]
guess[7] = np.median(x) if guess[7] is None else guess[7]
constraints_min[0] = np.min(x) if constraints_min[0] is None else constraints_min[0]
constraints_min[3] = np.min(x) if constraints_min[3] is None else constraints_min[3]
constraints_max[0] = np.max(x) if constraints_max[0] is None else constraints_max[0]
constraints_max[3] = np.max(x) if constraints_max[3] is None else constraints_max[3]
print('key', keys)
print('after', constraints_min)
def find_nearest(array, value):
array = np.asarray(array)
idx = (np.abs(array - value)).argmin()
return idx
def gaussian(x, g_cen, g_width, g_amp):
"""1-d gaussian: gaussian(x, amp, cen, wid)"""
return (g_amp / (np.sqrt(2.0 * np.pi) * g_width)) * np.exp(-(x - g_cen) ** 2 / (2 * g_width ** 2))
def lorentzian(x, l_cen, l_width, l_amp):
"""1d lorentzian"""
return (l_amp / (1 + ((1 * x - l_cen) / l_width) ** 2)) / (np.pi * l_width)
def background(x, slope, intercept):
"""background"""
return slope*x + intercept
mod = Model(gaussian) + Model(lorentzian) + Model(background)
params = Parameters()
params.add_many(('g_cen', x[int(peak_index)], bool(vary[0]), np.min(x), np.max(x), None, None),
('g_width', guess[1], bool(vary[1]), constraints_min[1], constraints_max[1], None, None),
('g_amp', guess[2], bool(vary[2]), constraints_min[2], constraints_max[2], None, None),
('l_cen', guess[3], bool(vary[3]), np.min(x), np.max(x), None, None),
('l_width', guess[4], bool(vary[4]), constraints_min[4], constraints_max[4], None, None),
('l_amp', guess[5], bool(vary[5]), constraints_min[5], constraints_max[5], None, None),
('slope', guess[6], bool(vary[6]), constraints_min[6], constraints_max[6], None, None),
('intercept', guess[7], bool(vary[7]), constraints_min[7], constraints_max[7], None, None))
result = mod.fit(y, params, x=x)
print('Chi-sqr: ', result.chisqr)
comps = result.eval_components()
gauss_3sigmamin = find_nearest(x, result.params['g_cen'].value-3*result.params['g_width'].value)
gauss_3sigmamax = find_nearest(x, result.params['g_cen'].value+3*result.params['g_width'].value)
numfit_min = gauss_3sigmamin if numfit_min is None else find_nearest(x, numfit_min)
numfit_max = gauss_3sigmamax if numfit_max is None else find_nearest(x, numfit_max)
print(numfit_max, numfit_min)
if x[numfit_min] < np.min(x):
numfit_min = gauss_3sigmamin
print('Minimal integration value outside of x range')
elif x[numfit_min] >= x[numfit_max]:
numfit_min = gauss_3sigmamin
print('Minimal integration value higher than maximal')
else:
pass
if x[numfit_max] > np.max(x):
numfit_max = gauss_3sigmamax
print('Maximal integration value outside of x range')
elif x[numfit_max] <= x[numfit_min]:
numfit_max = gauss_3sigmamax
print('Maximal integration value lower than minimal')
else:
pass
print(result.params['g_width'].value)
print(result.params['g_cen'].value)
num_int_area = simps(y[numfit_min:numfit_max], x[numfit_min:numfit_max])
num_int_bacground = integrate.quad(background, numfit_min, numfit_max, args=(result.params['slope'].value,result.params['intercept'].value))
plt.plot(x, y, 'b', label='Original data')
plt.plot(x, comps['gaussian'], 'r--', label='Gaussian component')
plt.fill_between(x, comps['gaussian'], facecolor="red", alpha=0.4)
plt.plot(x, comps['lorentzian'], 'b--', label='Lorentzian component')
plt.fill_between(x, comps['lorentzian'], facecolor="blue", alpha=0.4)
plt.plot(x, comps['background'], 'g--', label='Line component')
plt.fill_between(x, comps['background'], facecolor="green", alpha=0.4)
#plt.plot(x[numfit_min:numfit_max],y[numfit_min:numfit_max], 'vy', markersize=7)
plt.fill_between(x[numfit_min:numfit_max], y[numfit_min:numfit_max], facecolor="yellow", alpha=0.4, label='Integrated area')
#plt.plot(x, result.init_fit, 'k--', label='initial fit')
plt.plot(x, result.best_fit, 'k-', label='Best fit')
plt.title('%s \n Gaussian: centre = %9.4f, width = %9.4f, amp = %9.4f \n'
'Lorentzian: centre, %9.4f, width = %9.4f, amp = %9.4f \n'
'background: slope = %9.4f, intercept = %9.4f, int_area %9.4f' % (keys, result.params['g_cen'].value, result.params['g_width'].value,
result.params['g_amp'].value, result.params['l_cen'].value, result.params['l_width'].value,
result.params['l_amp'].value, result.params['slope'].value, result.params['intercept'].value, num_int_area))
plt.legend(loc='best')
plt.show()
d = {}
for pars in result.params:
d[str(pars)] = (result.params[str(pars)].value, result.params[str(pars)].vary)
d["int_area"] = num_int_area
d["int_background"] = num_int_bacground
d["full_report"] = result.fit_report()
data["Measurements"][str(keys)]["fit"] = d
return data
else:
return print('NO PEAK or more than 1 peak')