diff --git a/pyzebra/fit2.py b/pyzebra/fit2.py index 967bb41..720f95b 100644 --- a/pyzebra/fit2.py +++ b/pyzebra/fit2.py @@ -5,7 +5,9 @@ import scipy as sc from scipy import integrate import numpy as np -def fitccl(data, keys, guess, vary, constraints_min, constraints_max, numfit_min=None, numfit_max=None): +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 @@ -29,100 +31,110 @@ def fitccl(data, keys, guess, vary, constraints_min, constraints_max, numfit_min 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) + if len(data["Measurements"][str(keys)]["peak_indexes"]) != 1: + print("NO PEAK or more than 1 peak") + return - print('after', constraints_min) + 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 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 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 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 + 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)) + 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) + 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)) - 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 + 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: - return print('NO PEAK or more than 1 peak') + 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), + ) + 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