From 0347566aeb16fc7c84d10d01dc9b289c69c110b1 Mon Sep 17 00:00:00 2001 From: JakHolzer <53743814+JakHolzer@users.noreply.github.com> Date: Thu, 22 Oct 2020 15:16:40 +0200 Subject: [PATCH 1/4] Update comm_export.py fixed the order of hkls --- pyzebra/comm_export.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyzebra/comm_export.py b/pyzebra/comm_export.py index e178bb0..326792e 100644 --- a/pyzebra/comm_export.py +++ b/pyzebra/comm_export.py @@ -67,8 +67,8 @@ def export_comm(data, path, lorentz=False): line = ( scan_number_str + h_str - + l_str + k_str + + l_str + int_str + sigma_str + angle_str1 From 6099df650b4368aeff922ebc345e544a887d3288 Mon Sep 17 00:00:00 2001 From: JakHolzer <53743814+JakHolzer@users.noreply.github.com> Date: Fri, 23 Oct 2020 10:23:46 +0200 Subject: [PATCH 2/4] Update param_study_moduls.py Updated parametric study module with merging, adding etc... --- pyzebra/param_study_moduls.py | 211 +++++++++++++++++++++++++++++++--- 1 file changed, 195 insertions(+), 16 deletions(-) diff --git a/pyzebra/param_study_moduls.py b/pyzebra/param_study_moduls.py index 2593039..928fb12 100644 --- a/pyzebra/param_study_moduls.py +++ b/pyzebra/param_study_moduls.py @@ -1,5 +1,4 @@ from load_1D import load_1D -from ccl_dict_operation import add_dict import pandas as pd from mpl_toolkits.mplot3d import Axes3D # dont delete, otherwise waterfall wont work import matplotlib.pyplot as plt @@ -7,6 +6,17 @@ import matplotlib as mpl import numpy as np import pickle import scipy.io as sio +import uncertainties as u + + +def create_tuples(x, y, y_err): + """creates tuples for sorting and merginng of the data + Counts need to be normalized to monitor before""" + t = list() + for i in range(len(x)): + tup = (x[i], y[i], y_err[i]) + t.append(tup) + return t def load_dats(filepath): @@ -38,10 +48,10 @@ def load_dats(filepath): dict1 = add_dict(dict1, load_1D(file_list[i][0])) else: dict1 = add_dict(dict1, load_1D(file_list[i])) - dict1["scan"][i + 1]["params"] = {} + dict1["meas"][i + 1]["params"] = {} if data_type == "txt": for x in range(len(col_names) - 1): - dict1["scan"][i + 1]["params"][col_names[x + 1]] = file_list[i][x + 1] + dict1["meas"][i + 1]["params"][col_names[x + 1]] = file_list[i][x + 1] return dict1 @@ -53,7 +63,7 @@ def create_dataframe(dict1): # create dictionary to which we pull only wanted items before transforming it to pd.dataframe pull_dict = {} pull_dict["filenames"] = list() - for key in dict1["scan"][1]["params"]: + for key in dict1["meas"][1]["params"]: pull_dict[key] = list() pull_dict["temperature"] = list() pull_dict["mag_field"] = list() @@ -63,19 +73,19 @@ def create_dataframe(dict1): pull_dict["Counts"] = list() # populate the dict - for keys in dict1["scan"]: - if "file_of_origin" in dict1["scan"][keys]: - pull_dict["filenames"].append(dict1["scan"][keys]["file_of_origin"].split("/")[-1]) + for keys in dict1["meas"]: + if "file_of_origin" in dict1["meas"][keys]: + pull_dict["filenames"].append(dict1["meas"][keys]["file_of_origin"].split("/")[-1]) else: pull_dict["filenames"].append(dict1["meta"]["original_filename"].split("/")[-1]) - for key in dict1["scan"][keys]["params"]: - pull_dict[str(key)].append(float(dict1["scan"][keys]["params"][key])) - pull_dict["temperature"].append(dict1["scan"][keys]["temperature"]) - pull_dict["mag_field"].append(dict1["scan"][keys]["mag_field"]) - pull_dict["fit_area"].append(dict1["scan"][keys]["fit"]["fit_area"]) - pull_dict["int_area"].append(dict1["scan"][keys]["fit"]["int_area"]) - pull_dict["om"].append(dict1["scan"][keys]["om"]) - pull_dict["Counts"].append(dict1["scan"][keys]["Counts"]) + for key in dict1["meas"][keys]["params"]: + pull_dict[str(key)].append(float(dict1["meas"][keys]["params"][key])) + pull_dict["temperature"].append(dict1["meas"][keys]["temperature"]) + pull_dict["mag_field"].append(dict1["meas"][keys]["mag_field"]) + pull_dict["fit_area"].append(dict1["meas"][keys]["fit"]["fit_area"]) + pull_dict["int_area"].append(dict1["meas"][keys]["fit"]["int_area"]) + pull_dict["om"].append(dict1["meas"][keys]["om"]) + pull_dict["Counts"].append(dict1["meas"][keys]["Counts"]) return pd.DataFrame(data=pull_dict) @@ -144,7 +154,7 @@ def make_graph(data, sorting_parameter, style): def save_dict(obj, name): - """ saves dictionary as pickle file in binary format + """saves dictionary as pickle file in binary format :arg obj - object to save :arg name - name of the file NOTE: path should be added later""" @@ -200,3 +210,172 @@ def save_table(data, filetype, name, path=None): hdf.close() if filetype == "json": data.to_json((path + name + ".json")) + + def normalize(dict, key, monitor): + """Normalizes the measurement to monitor, checks if sigma exists, otherwise creates it + :arg dict : dictionary to from which to tkae the scan + :arg key : which scan to normalize from dict1 + :arg monitor : final monitor + :return counts - normalized counts + :return sigma - normalized sigma""" + + counts = np.array(dict["meas"][key]["Counts"]) + sigma = np.sqrt(counts) if "sigma" not in dict["meas"][key] else dict["meas"][key]["sigma"] + monitor_ratio = monitor / dict["meas"][key]["monitor"] + scaled_counts = counts * monitor_ratio + scaled_sigma = np.array(sigma) * monitor_ratio + + return scaled_counts, scaled_sigma + + def merge(dict1, dict2, scand_dict_result, auto=True, monitor=100000): + """merges the two tuples and sorts them, if om value is same, Counts value is average + averaging is propagated into sigma if dict1 == dict2, key[1] is deleted after merging + :arg dict1 : dictionary to which measurement will be merged + :arg dict2 : dictionary from which measurement will be merged + :arg keys : tuple with key to dict1 and dict2 + :arg auto : if true, when monitors are same, does not change it, if flase, takes monitor always + :arg monitor : final monitor after merging + note: dict1 and dict2 can be same dict + :return dict1 with merged scan""" + for keys in scand_dict_result: + for j in range(len(scand_dict_result[keys])): + first, second = scand_dict_result[keys][j][0], scand_dict_result[keys][j][1] + print(first, second) + if auto: + if dict1["meas"][first]["monitor"] == dict2["meas"][second]["monitor"]: + monitor = dict1["meas"][first]["monitor"] + + # load om and Counts + x1, x2 = dict1["meas"][first]["om"], dict2["meas"][second]["om"] + cor_y1, y_err1 = normalize(dict1, first, monitor=monitor) + cor_y2, y_err2 = normalize(dict2, second, monitor=monitor) + # creates touples (om, Counts, sigma) for sorting and further processing + tuple_list = create_tuples(x1, cor_y1, y_err1) + create_tuples(x2, cor_y2, y_err2) + # Sort the list on om and add 0 0 0 tuple to the last position + sorted_t = sorted(tuple_list, key=lambda tup: tup[0]) + sorted_t.append((0, 0, 0)) + om, Counts, sigma = [], [], [] + seen = list() + for i in range(len(sorted_t) - 1): + if sorted_t[i][0] not in seen: + if sorted_t[i][0] != sorted_t[i + 1][0]: + om = np.append(om, sorted_t[i][0]) + Counts = np.append(Counts, sorted_t[i][1]) + sigma = np.append(sigma, sorted_t[i][2]) + else: + om = np.append(om, sorted_t[i][0]) + counts1, counts2 = sorted_t[i][1], sorted_t[i + 1][1] + sigma1, sigma2 = sorted_t[i][2], sorted_t[i + 1][2] + count_err1 = u.ufloat(counts1, sigma1) + count_err2 = u.ufloat(counts2, sigma2) + avg = (count_err1 + count_err2) / 2 + Counts = np.append(Counts, avg.n) + sigma = np.append(sigma, avg.s) + seen.append(sorted_t[i][0]) + else: + continue + + if dict1 == dict2: + del dict1["meas"][second] + + note = ( + f"This measurement was merged with measurement {second} from " + f'file {dict2["meta"]["original_filename"]} \n' + ) + if "notes" not in dict1["meas"][first]: + dict1["meas"][first]["notes"] = note + else: + dict1["meas"][first]["notes"] += note + + dict1["meas"][first]["om"] = om + dict1["meas"][first]["Counts"] = Counts + dict1["meas"][first]["sigma"] = sigma + dict1["meas"][first]["monitor"] = monitor + print("merging done") + return dict1 + + +def add_dict(dict1, dict2): + """adds two dictionaries, meta of the new is saved as meata+original_filename and + measurements are shifted to continue with numbering of first dict + :arg dict1 : dictionarry to add to + :arg dict2 : dictionarry from which to take the measurements + :return dict1 : combined dictionary + Note: dict1 must be made from ccl, otherwise we would have to change the structure of loaded + dat file""" + if dict1["meta"]["zebra_mode"] != dict2["meta"]["zebra_mode"]: + print("You are trying to add scans measured with different zebra modes") + return + max_measurement_dict1 = max([keys for keys in dict1["meas"]]) + new_filenames = np.arange( + max_measurement_dict1 + 1, max_measurement_dict1 + 1 + len(dict2["meas"]) + ) + new_meta_name = "meta" + str(dict2["meta"]["original_filename"]) + if new_meta_name not in dict1: + for keys, name in zip(dict2["meas"], new_filenames): + dict2["meas"][keys]["file_of_origin"] = str(dict2["meta"]["original_filename"]) + dict1["meas"][name] = dict2["meas"][keys] + + dict1[new_meta_name] = dict2["meta"] + else: + raise KeyError( + str( + "The file %s has alredy been added to %s" + % (dict2["meta"]["original_filename"], dict1["meta"]["original_filename"]) + ) + ) + return dict1 + + +def auto(dict): + """takes just unique tuples from all tuples in dictionary returend by scan_dict + intendet for automatic merge if you doesent want to specify what scans to merge together + args: dict - dictionary from scan_dict function + :return dict - dict without repetitions""" + for keys in dict: + tuple_list = dict[keys] + new = list() + for i in range(len(tuple_list)): + if tuple_list[0][0] == tuple_list[i][0]: + new.append(tuple_list[i]) + dict[keys] = new + return dict + + +def scan_dict(dict, precision=0.5): + """scans dictionary for duplicate angles indexes + :arg dict : dictionary to scan + :arg precision : in deg, sometimes angles are zero so its easier this way, instead of + checking zero division + :return dictionary with matching scans, if there are none, the dict is empty + note: can be checked by "not d", true if empty + """ + if dict["meta"]["zebra_mode"] == "bi": + angles = ["twotheta_angle", "omega_angle", "chi_angle", "phi_angle"] + elif dict["meta"]["zebra_mode"] == "nb": + angles = ["gamma_angle", "omega_angle", "nu_angle"] + else: + print("Unknown zebra mode") + return + + d = {} + for i in dict["meas"]: + for j in dict["meas"]: + if dict["meas"][i] != dict["meas"][j]: + itup = list() + for k in angles: + itup.append(abs(abs(dict["meas"][i][k]) - abs(dict["meas"][j][k]))) + + if all(i <= precision for i in itup): + if str([np.around(dict["meas"][i][k], 1) for k in angles]) not in d: + d[str([np.around(dict["meas"][i][k], 1) for k in angles])] = list() + d[str([np.around(dict["meas"][i][k], 1) for k in angles])].append((i, j)) + else: + d[str([np.around(dict["meas"][i][k], 1) for k in angles])].append((i, j)) + + else: + pass + + else: + continue + return d From 6ff1b2b54f4b37586e1b13fcd91b93d6279036b3 Mon Sep 17 00:00:00 2001 From: JakHolzer <53743814+JakHolzer@users.noreply.github.com> Date: Fri, 23 Oct 2020 15:19:08 +0200 Subject: [PATCH 3/4] Update param_study_moduls.py meas to scan --- pyzebra/param_study_moduls.py | 92 ++++++++++++++++++----------------- 1 file changed, 47 insertions(+), 45 deletions(-) diff --git a/pyzebra/param_study_moduls.py b/pyzebra/param_study_moduls.py index 928fb12..a8cfb56 100644 --- a/pyzebra/param_study_moduls.py +++ b/pyzebra/param_study_moduls.py @@ -48,10 +48,10 @@ def load_dats(filepath): dict1 = add_dict(dict1, load_1D(file_list[i][0])) else: dict1 = add_dict(dict1, load_1D(file_list[i])) - dict1["meas"][i + 1]["params"] = {} + dict1["scan"][i + 1]["params"] = {} if data_type == "txt": for x in range(len(col_names) - 1): - dict1["meas"][i + 1]["params"][col_names[x + 1]] = file_list[i][x + 1] + dict1["scan"][i + 1]["params"][col_names[x + 1]] = file_list[i][x + 1] return dict1 @@ -63,7 +63,7 @@ def create_dataframe(dict1): # create dictionary to which we pull only wanted items before transforming it to pd.dataframe pull_dict = {} pull_dict["filenames"] = list() - for key in dict1["meas"][1]["params"]: + for key in dict1["scan"][1]["params"]: pull_dict[key] = list() pull_dict["temperature"] = list() pull_dict["mag_field"] = list() @@ -73,19 +73,19 @@ def create_dataframe(dict1): pull_dict["Counts"] = list() # populate the dict - for keys in dict1["meas"]: - if "file_of_origin" in dict1["meas"][keys]: - pull_dict["filenames"].append(dict1["meas"][keys]["file_of_origin"].split("/")[-1]) + for keys in dict1["scan"]: + if "file_of_origin" in dict1["scan"][keys]: + pull_dict["filenames"].append(dict1["scan"][keys]["file_of_origin"].split("/")[-1]) else: pull_dict["filenames"].append(dict1["meta"]["original_filename"].split("/")[-1]) - for key in dict1["meas"][keys]["params"]: - pull_dict[str(key)].append(float(dict1["meas"][keys]["params"][key])) - pull_dict["temperature"].append(dict1["meas"][keys]["temperature"]) - pull_dict["mag_field"].append(dict1["meas"][keys]["mag_field"]) - pull_dict["fit_area"].append(dict1["meas"][keys]["fit"]["fit_area"]) - pull_dict["int_area"].append(dict1["meas"][keys]["fit"]["int_area"]) - pull_dict["om"].append(dict1["meas"][keys]["om"]) - pull_dict["Counts"].append(dict1["meas"][keys]["Counts"]) + for key in dict1["scan"][keys]["params"]: + pull_dict[str(key)].append(float(dict1["scan"][keys]["params"][key])) + pull_dict["temperature"].append(dict1["scan"][keys]["temperature"]) + pull_dict["mag_field"].append(dict1["scan"][keys]["mag_field"]) + pull_dict["fit_area"].append(dict1["scan"][keys]["fit"]["fit_area"]) + pull_dict["int_area"].append(dict1["scan"][keys]["fit"]["int_area"]) + pull_dict["om"].append(dict1["scan"][keys]["om"]) + pull_dict["Counts"].append(dict1["scan"][keys]["Counts"]) return pd.DataFrame(data=pull_dict) @@ -219,21 +219,22 @@ def save_table(data, filetype, name, path=None): :return counts - normalized counts :return sigma - normalized sigma""" - counts = np.array(dict["meas"][key]["Counts"]) - sigma = np.sqrt(counts) if "sigma" not in dict["meas"][key] else dict["meas"][key]["sigma"] - monitor_ratio = monitor / dict["meas"][key]["monitor"] + counts = np.array(dict["scan"][key]["Counts"]) + sigma = np.sqrt(counts) if "sigma" not in dict["scan"][key] else dict["scan"][key]["sigma"] + monitor_ratio = monitor / dict["scan"][key]["monitor"] scaled_counts = counts * monitor_ratio scaled_sigma = np.array(sigma) * monitor_ratio return scaled_counts, scaled_sigma - def merge(dict1, dict2, scand_dict_result, auto=True, monitor=100000): + def merge(dict1, dict2, scand_dict_result, keep=True, monitor=100000): """merges the two tuples and sorts them, if om value is same, Counts value is average averaging is propagated into sigma if dict1 == dict2, key[1] is deleted after merging :arg dict1 : dictionary to which measurement will be merged :arg dict2 : dictionary from which measurement will be merged - :arg keys : tuple with key to dict1 and dict2 - :arg auto : if true, when monitors are same, does not change it, if flase, takes monitor always + :arg scand_dict_result : result of scan_dict after auto function + :arg keep : if true, when monitors are same, does not change it, if flase, takes monitor + always :arg monitor : final monitor after merging note: dict1 and dict2 can be same dict :return dict1 with merged scan""" @@ -241,12 +242,12 @@ def save_table(data, filetype, name, path=None): for j in range(len(scand_dict_result[keys])): first, second = scand_dict_result[keys][j][0], scand_dict_result[keys][j][1] print(first, second) - if auto: - if dict1["meas"][first]["monitor"] == dict2["meas"][second]["monitor"]: - monitor = dict1["meas"][first]["monitor"] + if keep: + if dict1["scan"][first]["monitor"] == dict2["scan"][second]["monitor"]: + monitor = dict1["scan"][first]["monitor"] # load om and Counts - x1, x2 = dict1["meas"][first]["om"], dict2["meas"][second]["om"] + x1, x2 = dict1["scan"][first]["om"], dict2["scan"][second]["om"] cor_y1, y_err1 = normalize(dict1, first, monitor=monitor) cor_y2, y_err2 = normalize(dict2, second, monitor=monitor) # creates touples (om, Counts, sigma) for sorting and further processing @@ -276,21 +277,21 @@ def save_table(data, filetype, name, path=None): continue if dict1 == dict2: - del dict1["meas"][second] + del dict1["scan"][second] note = ( f"This measurement was merged with measurement {second} from " f'file {dict2["meta"]["original_filename"]} \n' ) - if "notes" not in dict1["meas"][first]: - dict1["meas"][first]["notes"] = note + if "notes" not in dict1["scan"][first]: + dict1["scan"][first]["notes"] = note else: - dict1["meas"][first]["notes"] += note + dict1["scan"][first]["notes"] += note - dict1["meas"][first]["om"] = om - dict1["meas"][first]["Counts"] = Counts - dict1["meas"][first]["sigma"] = sigma - dict1["meas"][first]["monitor"] = monitor + dict1["scan"][first]["om"] = om + dict1["scan"][first]["Counts"] = Counts + dict1["scan"][first]["sigma"] = sigma + dict1["scan"][first]["monitor"] = monitor print("merging done") return dict1 @@ -306,15 +307,15 @@ def add_dict(dict1, dict2): if dict1["meta"]["zebra_mode"] != dict2["meta"]["zebra_mode"]: print("You are trying to add scans measured with different zebra modes") return - max_measurement_dict1 = max([keys for keys in dict1["meas"]]) + max_measurement_dict1 = max([keys for keys in dict1["scan"]]) new_filenames = np.arange( - max_measurement_dict1 + 1, max_measurement_dict1 + 1 + len(dict2["meas"]) + max_measurement_dict1 + 1, max_measurement_dict1 + 1 + len(dict2["scan"]) ) new_meta_name = "meta" + str(dict2["meta"]["original_filename"]) if new_meta_name not in dict1: - for keys, name in zip(dict2["meas"], new_filenames): - dict2["meas"][keys]["file_of_origin"] = str(dict2["meta"]["original_filename"]) - dict1["meas"][name] = dict2["meas"][keys] + for keys, name in zip(dict2["scan"], new_filenames): + dict2["scan"][keys]["file_of_origin"] = str(dict2["meta"]["original_filename"]) + dict1["scan"][name] = dict2["scan"][keys] dict1[new_meta_name] = dict2["meta"] else: @@ -350,6 +351,7 @@ def scan_dict(dict, precision=0.5): :return dictionary with matching scans, if there are none, the dict is empty note: can be checked by "not d", true if empty """ + if dict["meta"]["zebra_mode"] == "bi": angles = ["twotheta_angle", "omega_angle", "chi_angle", "phi_angle"] elif dict["meta"]["zebra_mode"] == "nb": @@ -359,19 +361,19 @@ def scan_dict(dict, precision=0.5): return d = {} - for i in dict["meas"]: - for j in dict["meas"]: - if dict["meas"][i] != dict["meas"][j]: + for i in dict["scan"]: + for j in dict["scan"]: + if dict["scan"][i] != dict["scan"][j]: itup = list() for k in angles: - itup.append(abs(abs(dict["meas"][i][k]) - abs(dict["meas"][j][k]))) + itup.append(abs(abs(dict["scan"][i][k]) - abs(dict["scan"][j][k]))) if all(i <= precision for i in itup): - if str([np.around(dict["meas"][i][k], 1) for k in angles]) not in d: - d[str([np.around(dict["meas"][i][k], 1) for k in angles])] = list() - d[str([np.around(dict["meas"][i][k], 1) for k in angles])].append((i, j)) + if str([np.around(dict["scan"][i][k], 1) for k in angles]) not in d: + d[str([np.around(dict["scan"][i][k], 1) for k in angles])] = list() + d[str([np.around(dict["scan"][i][k], 1) for k in angles])].append((i, j)) else: - d[str([np.around(dict["meas"][i][k], 1) for k in angles])].append((i, j)) + d[str([np.around(dict["scan"][i][k], 1) for k in angles])].append((i, j)) else: pass From 430ffc2caa45d007f1d10a78797586dd2a3a91d0 Mon Sep 17 00:00:00 2001 From: JakHolzer <53743814+JakHolzer@users.noreply.github.com> Date: Fri, 23 Oct 2020 16:45:04 +0200 Subject: [PATCH 4/4] Generalized fitting function This is first idea how the function could work. Use should be the same as previous one but we need to find a way how to pass parameters to the function. There is a new parameter called variable, which should choose the x coordinate, since "om" might not be the only axis here. Function does not change the initial dictionary yet, but process will be the same as in the first one. It is still not clear how the peaks should be reported, more so what to report in case of two overlapping peaks (same goes for numerical integration), but the process will be similar to the fitvol2. The function can be used, but is posted here for a reason of discussion and finding the best way of passing the parameters. --- pyzebra/fitvol3.py | 167 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 167 insertions(+) create mode 100644 pyzebra/fitvol3.py diff --git a/pyzebra/fitvol3.py b/pyzebra/fitvol3.py new file mode 100644 index 0000000..e3cf6e4 --- /dev/null +++ b/pyzebra/fitvol3.py @@ -0,0 +1,167 @@ +import numpy as np +from lmfit import Model, Parameters +from scipy.integrate import simps +import matplotlib.pyplot as plt +import uncertainties as u +from lmfit.models import GaussianModel +from lmfit.models import VoigtModel +from lmfit.models import PseudoVoigtModel + + +def bin_data(array, binsize): + if isinstance(binsize, int) and 0 < binsize < len(array): + return [ + np.mean(array[binsize * i : binsize * i + binsize]) + for i in range(int(np.ceil(len(array) / binsize))) + ] + else: + print("Binsize need to be positive integer smaller than lenght of array") + return array + + +def create_uncertanities(y, y_err): + # create array with uncertanities for error propagation + combined = np.array([]) + for i in range(len(y)): + part = u.ufloat(y[i], y_err[i]) + combined = np.append(combined, part) + return combined + + +def find_nearest(array, value): + # find nearest value and return index + array = np.asarray(array) + idx = (np.abs(array - value)).argmin() + return idx + + +# predefined peak positions +# peaks = [6.2, 8.1, 9.9, 11.5] +peaks = [23.5, 24.5] +# peaks = [24] +def fitccl(scan, variable="om", peak_type="gauss", binning=None): + + x = list(scan[variable]) + y = list(scan["Counts"]) + peak_centre = np.mean(x) + if binning is None or binning == 0 or binning == 1: + x = list(scan["om"]) + y = list(scan["Counts"]) + y_err = list(np.sqrt(y)) if scan.get("sigma", None) is None else list(scan["sigma"]) + print(scan["peak_indexes"]) + if not scan["peak_indexes"]: + peak_centre = np.mean(x) + else: + centre = x[int(scan["peak_indexes"])] + else: + x = list(scan["om"]) + if not scan["peak_indexes"]: + peak_centre = np.mean(x) + else: + peak_centre = x[int(scan["peak_indexes"])] + x = bin_data(x, binning) + y = list(scan["Counts"]) + y_err = list(np.sqrt(y)) if scan.get("sigma", None) is None else list(scan["sigma"]) + combined = bin_data(create_uncertanities(y, y_err), binning) + y = [combined[i].n for i in range(len(combined))] + y_err = [combined[i].s for i in range(len(combined))] + + def background(x, slope, intercept): + """background""" + return slope * (x - peak_centre) + intercept + + def gaussian(x, center, g_sigma, amplitude): + """1-d gaussian: gaussian(x, amp, cen, wid)""" + return (amplitude / (np.sqrt(2.0 * np.pi) * g_sigma)) * np.exp( + -((x - center) ** 2) / (2 * g_sigma ** 2) + ) + + def lorentzian(x, center, l_sigma, amplitude): + """1d lorentzian""" + return (amplitude / (1 + ((1 * x - center) / l_sigma) ** 2)) / (np.pi * l_sigma) + + def pseudoVoigt1(x, center, g_sigma, amplitude, l_sigma, fraction): + """PseudoVoight peak with different widths of lorenzian and gaussian""" + return (1 - fraction) * gaussian(x, center, g_sigma, amplitude) + fraction * ( + lorentzian(x, center, l_sigma, amplitude) + ) + + mod = Model(background) + params = Parameters() + params.add_many( + ("slope", 0, True, None, None, None, None), ("intercept", 0, False, None, None, None, None) + ) + for i in range(len(peaks)): + if peak_type == "gauss": + mod = mod + GaussianModel(prefix="p%d_" % (i + 1)) + params.add(str("p%d_" % (i + 1) + "amplitude"), 20, True, 0, None, None) + params.add(str("p%d_" % (i + 1) + "center"), peaks[i], True, None, None, None) + params.add(str("p%d_" % (i + 1) + "sigma"), 0.2, True, 0, 5, None) + elif peak_type == "voigt": + mod = mod + VoigtModel(prefix="p%d_" % (i + 1)) + params.add(str("p%d_" % (i + 1) + "amplitude"), 20, True, 0, None, None) + params.add(str("p%d_" % (i + 1) + "center"), peaks[i], True, None, None, None) + params.add(str("p%d_" % (i + 1) + "sigma"), 0.2, True, 0, 3, None) + params.add(str("p%d_" % (i + 1) + "gamma"), 0.2, True, 0, 5, None) + elif peak_type == "pseudovoigt": + mod = mod + PseudoVoigtModel(prefix="p%d_" % (i + 1)) + params.add(str("p%d_" % (i + 1) + "amplitude"), 20, True, 0, None, None) + params.add(str("p%d_" % (i + 1) + "center"), peaks[i], True, None, None, None) + params.add(str("p%d_" % (i + 1) + "sigma"), 0.2, True, 0, 5, None) + params.add(str("p%d_" % (i + 1) + "fraction"), 0.5, True, -5, 5, None) + elif peak_type == "pseudovoigt1": + mod = mod + Model(pseudoVoigt1, prefix="p%d_" % (i + 1)) + params.add(str("p%d_" % (i + 1) + "amplitude"), 20, True, 0, None, None) + params.add(str("p%d_" % (i + 1) + "center"), peaks[i], True, None, None, None) + params.add(str("p%d_" % (i + 1) + "g_sigma"), 0.2, True, 0, 5, None) + params.add(str("p%d_" % (i + 1) + "l_sigma"), 0.2, True, 0, 5, None) + params.add(str("p%d_" % (i + 1) + "fraction"), 0.5, True, 0, 1, None) + # add parameters + + result = mod.fit( + y, params, weights=[np.abs(1 / y_err[i]) for i in range(len(y_err))], x=x, calc_covar=True + ) + + comps = result.eval_components() + + reportstring = list() + for keys in result.params: + if result.params[keys].value is not None: + str2 = np.around(result.params[keys].value, 3) + else: + str2 = 0 + if result.params[keys].stderr is not None: + str3 = np.around(result.params[keys].stderr, 3) + else: + str3 = 0 + reportstring.append("%s = %2.3f +/- %2.3f" % (keys, str2, str3)) + + reportstring = "\n".join(reportstring) + + plt.figure(figsize=(20, 10)) + plt.plot(x, result.best_fit, "k-", label="Best fit") + + plt.plot(x, y, "b-", label="Original data") + plt.plot(x, comps["background"], "g--", label="Line component") + for i in range(len(peaks)): + plt.plot( + x, + comps[str("p%d_" % (i + 1))], + "r--", + ) + plt.fill_between(x, comps[str("p%d_" % (i + 1))], alpha=0.4, label=str("p%d_" % (i + 1))) + plt.legend() + plt.text( + np.min(x), + np.max(y), + reportstring, + fontsize=9, + verticalalignment="top", + ) + plt.title(str(peak_type)) + + plt.xlabel("Omega [deg]") + plt.ylabel("Counts [a.u.]") + plt.show() + + print(result.fit_report())