diff --git a/docs/instruments/loki/nurf/ill_auxilliary_funcs.py b/docs/instruments/loki/nurf/ill_auxilliary_funcs.py new file mode 100644 index 000000000..31ee85a85 --- /dev/null +++ b/docs/instruments/loki/nurf/ill_auxilliary_funcs.py @@ -0,0 +1,438 @@ +import errno +import os +import h5py +import numpy as np +import copy + +def complete_fname(scan_numbers): + """Converts a list of input numbers to a filename uses at ILL. + + Parameters + ---------- + scan_numbers: list of int or a single int + List of filenumbers or one filenumnber. + + Returns: + ---------- + flist_num: list of str or one str + List of filenames following ILL style or string following ILL style. + + """ + if isinstance(scan_numbers, int): + flist_num = f"{str(scan_numbers).zfill(6)}.nxs" + + if isinstance(scan_numbers, list): + # convert a list of input numbers to real filename + flist_num = [str(i).zfill(6) + ".nxs" for i in scan_numbers] + + return flist_num + +# This cell contains function to convert ILL raw data to the forseen Loki file format + +def loki_file_creator(loki_path, loki_filename): + """ + This function creates a dummy loki_basic.nxs that has the expected file structure for ESS Loki. + Afterwards, I will append the nurf data and rename the file. + """ + with h5py.File(os.path.join(loki_path,loki_filename), 'w') as hf: + nxentry = hf.create_group("entry") + nxentry.attrs['NX_class'] = 'NXentry' + nxinstrument=nxentry.create_group("instrument") + nxinstrument.attrs['NX_class'] = 'NXinstrument' + + +def load_one_spectro_file(file_handle, path_rawdata): + """ + This function loads one .nxs file containing spectroscopy data (fluo, uv). Data is stored in multiple np.ndarrays. + + In: + file_handle: file_handle is a file number, one element is expected otherwise an error is raised + type: list + + path_rawdata: Path to the raw data + type: str + + Out: + data: contains all relevant HDF5 entries and their content for the Nurf project (keys and values) + type: dict + + """ + + # create path to file, convert file_number to string + file_path_spectro = os.path.join(path_rawdata, file_handle + '.nxs') + + # check if file exists + if not os.path.isfile(file_path_spectro): + raise FileNotFoundError(errno.ENOENT, os.strerror(errno.ENOENT), + file_path_spectro) + + # we are ready to load the data set + # open the .nxs file and read the values + with h5py.File(file_path_spectro, "r") as f: + # access nurf sub-group + nurf_group = '/entry0/D22/nurf/' + + # access keys in sub-group + nurf_keys = list(f[nurf_group].keys()) + # how many keys exist + len_nurf_keys = len(nurf_keys) + + # print(nurf_keys) + # print(len_nurf_keys) + + # this returns a list with HDF5datasets (I think) + # data_spectro_file=list(f[nurf_group].values()) + + # data_spectro_file=f[nurf_group].values() + # data_spectro_file=f[nurf_group] + + # extract all data of the nurf subgroup and store it in a new dict + + # initialise an empty dict + data = {} + + for key in f[nurf_group].keys(): + # print(key) + # this is how I get string giving the full path to this dataset + path_dataset = f[nurf_group][key].name + # print(path_dataset) + # print(f[nurf_group][key].name) #this prints the full path to the dataset + # print(type(f[path_dataset][:])) #this gives me access to the content of the data set, I could use : or () inside [], out: np.ndarray + + # This gives a dict with full path name as dict entry followed by value. No, I don't want this, but good to know. + # data[f[nurf_group][key].name]=f[path_dataset][:] + + # This gives a dict where the keys corresponds to the key names of the h5 file. + data[key] = f[path_dataset][:] + + # print(f[nurf_group].get('Fluo_spectra')) + + # print a hierachical view of the file (simple) + # like this is would only go through subgroups + # f[nurf_group].visititems(lambda x, y: print(x)) + + # walk through the whole file and show the attributes (found on github as an explanation for this function) + # def print_attrs(name, obj): + # print(name) + # for key, val in obj.attrs.items(): + # print("{0}: {1}".format(key, val)) + # f[nurf_group].visititems(print_attrs) + + # print(data_spectro_file) + + # file_handle is returned as np.ndarray and its an np.array, elements correspond to row indices + return data + + +def nurf_file_creator(loki_file, path_to_loki_file, data): + """ + Appends NUrF group to LOKI NeXus file for ESS + + Args: + loki_file (str): filename of NeXus file for Loki + path_to_loki_file (str): File path where the NeXus file for LOKI is stored + data (dict): Dictionary with dummy data for Nurf + """ + + # change directory where the loki.nxs is located + os.chdir(path_to_loki_file) + + # open the file and append + with h5py.File(loki_file, 'a') as hf: + + #comment on names + #UV/FL_Background is the dark + #UV/FL_Intensity0 is the reference + #UV/FL_Spectra is the sample + + # image_key: number of frames (nFrames) given indirectly as part of the shape of the arrays + # TODO: keep in mind what happens if multiple dark or reference frames are taken + + # remove axis=2 of length one from this array + data['UV_spectra']=np.squeeze(data['UV_spectra'], axis=2) #this removes the third axis in this array, TODO: needs later to be verified with real data from hardware + + # assemble all spectra in one variable + uv_all_data=np.row_stack((data['UV_spectra'],data['UV_background'], data['UV_intensity0'])) + + dummy_vec=np.full(np.shape(uv_all_data)[0], False) + + #create boolean masks for data, dark, reference + uv_nb_spectra=np.shape(data['UV_spectra'])[0] + + # data mask, copy and replace first entries + uv_data_mask=copy.copy(dummy_vec) + uv_data_mask[0:uv_nb_spectra]=True + + # dark mask + # find out how many darks exist + # TODO: Is there always a background or do we need to catch this case if there isn't? + if data['UV_background'].ndim==1: + uv_nb_darks=1 + else: + uv_nb_darks=np.shape(data['UV_background'])[1] #TODO: needs to be verified with real data from Judith's setup + + uv_dark_mask=copy.copy(dummy_vec) + uv_dark_mask[uv_nb_spectra:uv_nb_spectra+uv_nb_darks]=True + + # reference + # how many references where taken? + if data['UV_intensity0'].ndim==1: + uv_nb_ref=1 + else: + uv_nb_ref=np.shape(data['UV_intensity0'])[1] #TODO: needs to be verified with real data from Judith's setup + + # reference mask, copy and replace first entries + uv_ref_mask=copy.copy(dummy_vec) + uv_ref_mask[uv_nb_spectra+uv_nb_darks:uv_nb_spectra+uv_nb_darks+uv_nb_ref]=True + + + + # UV subgroup + grp_uv = hf.create_group("/entry/instrument/uv") + grp_uv.attrs["NX_class"] = 'NXdata' + + # uv spectra + uv_signal_data=grp_uv.create_dataset('data', data=uv_all_data, dtype=np.float32) + uv_signal_data.attrs['long name']= 'all_data' + uv_signal_data.attrs['units']= 'counts' + grp_uv.attrs['signal']= 'data' #indicate that the main signal is data + grp_uv.attrs['axes']= [ "spectrum", "wavelength" ] #time is here the first axis, i.e axis=0, wavelength is axis=1 + + # define the AXISNAME_indices + grp_uv.attrs['time_indices'] = 0 + grp_uv.attrs['integration_time_indices'] = 0 + grp_uv.attrs['wavelength_indices'] = 1 + + # introducing a key that is interpretable for sample, dark, and reference + grp_uv.attrs['is_data_indices'] = 0 + grp_uv.attrs['is_dark_indices'] = 0 + grp_uv.attrs['is_reference_indices'] = 0 + + grp_uv.create_dataset('is_data', data=uv_data_mask, dtype=bool) + grp_uv.create_dataset('is_dark', data=uv_dark_mask, dtype=bool) + grp_uv.create_dataset('is_reference', data=uv_ref_mask, dtype=bool) + + # uv_time + # dummy timestamps for uv_time + # TODO: Codes will have to change later for the real hardware. + uv_time = np.empty(np.shape(uv_all_data)[0], dtype='datetime64[us]') + for i in range(0, np.shape(uv_time)[0]): + uv_time[i]=np.datetime64('now') + + # see https://stackoverflow.com/questions/23570632/store-datetimes-in-hdf5-with-h5py + # suggested work around because h5py does not support time types + uv_time_data=grp_uv.create_dataset('time', data=uv_time.view(' 066150.nxs 066150\n", + "066017 \n", + "066020 \n", + "066023 \n", + "066026 \n", + "066029 \n", + "066032 \n", + "066034 \n", + "066037 \n", + "066040 \n", + "066043 \n", + "066046 \n", + "066050 \n", + "066053 \n", + "066056 \n", + "066059 \n", + "066062 \n", + "066065 \n", + "066068 \n", + "066071 \n", + "066074 \n", + "066077 \n", + "066080 \n", + "066083 \n", + "066086 \n", + "066089 \n", + "066092 \n", + "066095 \n", + "066098 \n", + "066101 \n", + "066104 \n", + "066107 \n", + "066110 \n", + "066113 \n", + "066116 \n", + "066119 \n", + "066122 \n", + "066125 \n", + "066128 \n", + "066131 \n", + "066134 \n", + "066137 \n", + "066140 \n", + "066143 \n", + "066146 \n", + "065925 \n", + "065927 \n", + "065930 \n", + "065933 \n", + "065936 \n", + "065939 \n", + "065942 \n", + "065945 \n", + "065948 \n", + "065951 \n", + "065954 \n", + "065957 \n", + "065962 \n", + "065965 \n", + "065968 \n", + "065971 \n", + "065974 \n", + "065977 \n", + "065980 \n", + "065983 \n", + "065986 \n", + "065989 \n", + "065992 \n" + ] + } + ], + "source": [ + "# Location processed files\n", + "process_folder='/Users/gudlo523/Library/CloudStorage/OneDrive-LundUniversity/UU/ILL_947_June_2021/212/d22/exp_9-13-947/processed/ess_version'\n", + "\n", + "# Location for LOKI_basic.nxs\n", + "path_to_loki_file=process_folder\n", + "\n", + "# create LOKI_basic.nxs\n", + "loki_file='LOKI_basic.nxs'\n", + "loki_file_creator(path_to_loki_file, loki_file)\n", + "\n", + "# Location raw data files\n", + "path_rawdata='/Users/gudlo523/Library/CloudStorage/OneDrive-LundUniversity/UU/ILL_947_June_2021/212/d22/exp_9-13-947/rawdata'\n", + "os.chdir(path_rawdata)\n", + "\n", + "\n", + "flist=os.listdir(path_rawdata)\n", + "print(type(flist), flist[0], os.path.splitext(flist[0])[0])\n", + "# get only filenumbers\n", + "flist_num=[os.path.splitext(fnumber)[0] for fnumber in flist]\n", + "\n", + "# file numbers for ILL experiment 947\n", + "#scan numbers in exp 5 \n", + "exp5=[66017, 66020, 66023, 66026, 66029, 66032, 66034, 66037, 66040, 66043, 66046]\n", + "#scan numbers in exp 6\n", + "exp6= [66050, 66053, 66056, 66059, 66062, 66065, 66068, 66071, 66074, 66077, 66080]\n", + "scan_numbers=[exp5, exp6]\n", + "for i in scan_numbers:\n", + " convert_ill2loki(i, path_to_loki_file, loki_file, path_rawdata, process_folder)\n", + "\n", + "\n", + "# scan number in exp 7 and exp8\n", + "exp7= [66083, 66086, 66089, 66092, 66095, 66098, 66101, 66104, 66107, 66110, 66113]\n", + "exp8= [66116, 66119, 66122, 66125, 66128, 66131, 66134, 66137, 66140, 66143, 66146]\n", + "\n", + "scan_numbers= [exp7, exp8]\n", + "\n", + "for i in scan_numbers:\n", + " convert_ill2loki(i, path_to_loki_file, loki_file, path_rawdata, process_folder)\n", + "\n", + "#exp2 \n", + "exp2=[65925, 65927, 65930, 65933, 65936, 65939, 65942, 65945, 65948, 65951, 65954, 65957]\n", + "exp3= [65962, 65965, 65968, 65971, 65974, 65977, 65980, 65983, 65986, 65989, 65992]\n", + "\n", + "scan_numbers= [exp2, exp3]\n", + "\n", + "for i in scan_numbers:\n", + " convert_ill2loki(i, path_to_loki_file, loki_file, path_rawdata, process_folder)\n", + "\n", + "\n", + "\n", + "\n", + "\n", + " " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3.8.12 ('scippneutron')", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.8.12" + }, + "orig_nbformat": 4, + "vscode": { + "interpreter": { + "hash": "410bfb39d85b3112e66f48ab3e53bb74ad7e7b3fb364f756ca860b5a3cf79ca2" + } + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/docs/instruments/loki/nurf/nurf.py b/docs/instruments/loki/nurf/nurf.py new file mode 100644 index 000000000..bc3206b56 --- /dev/null +++ b/docs/instruments/loki/nurf/nurf.py @@ -0,0 +1,1953 @@ +# standard library imports +import itertools +import os +from typing import Optional + +# related third party imports +import numpy as np +import matplotlib.pyplot as plt +from matplotlib import cm +import matplotlib.gridspec as gridspec +from IPython.display import display, HTML +# possibilites for median filters, I did not benchmark, apparently median_filter +# could be the faster and medfilt2d is faster than medfilt +from scipy.ndimage import median_filter +from scipy.signal import medfilt +from scipy.optimize import leastsq # needed for fitting of turbidity + +# local application imports +import scippneutron as scn +import scippnexus as snx +import scipp as sc + + +def split_sample_dark_reference(da): + """Separate incoming dataarray into the three contributions: sample, dark, reference. + + Parameters + ---------- + da: scipp.DataArray + sc.DataArray that contains spectroscopy contributions sample, dark, reference + + Returns: + ---------- + da_dict: dict + Dictionary that contains spectroscopy data signal (data) from the sample, + the reference, and the dark measurement. + Keys: sample, reference, dark + + """ + assert isinstance(da, sc.DataArray) + + dark = da[da.coords["is_dark"]].squeeze() + ref = da[da.coords["is_reference"]].squeeze() + sample = da[da.coords["is_data"]] + + #TODO Instead of a dict a sc.Dataset? + return {"sample": sample, "reference": ref, "dark": dark} + + +def load_uv(name): + """Loads the UV data from the corresponding entry in the LoKI.nxs filename. + + Parameters + ---------- + name: str + Filename, e.g. 066017.nxs + + Returns: + ---------- + uv_dict: dict + Dictionary that contains UV data signal (data) from the sample, the reference, + and the dark measurement. + Keys: sample, reference, dark + + """ + + # load the nexus and extract the uv entry + with snx.File(name) as f: + uv = f["entry/instrument/uv"][()] + + # separation + uv_dict = split_sample_dark_reference(uv) + + return uv_dict + + +def normalize_uv( + *, sample: sc.DataArray, reference: sc.DataArray, dark: sc.DataArray +) -> sc.DataArray: + """Calculates the absorbance of the UV signal. + + Parameters + ---------- + sample: sc.DataArray + DataArray containing sample UV signal, one spectrum or multiple spectra. + reference: sc.DataArray + DataArray containing reference UV signal, one spectrum expected. + dark: sc.DataArray + DataArray containing dark UV signal, one spectrum expected. + + Returns + ---------- + normalized: sc.DataArray + DataArray that contains the normalized UV signal, one spectrum or mulitple spectra. + + """ + + normalized = sc.log10( + (reference - dark) / (sample - dark) + ) # results in DataArrays with multiple spectra + + return normalized + + +def load_and_normalize_uv(name): + """Loads the UV data from the corresponding entry in the LoKI.nxs filename and calculates the absorbance of each UV spectrum. + For an averaged spectrum based on all UV spectra in the file, use process_uv. + + Parameters + ---------- + name: str + Filename, e.g. 066017.nxs + + Returns: + ---------- + normalized: sc.DataArray + DataArray that contains the normalized UV signal, one spectrum or mulitple spectra. + + """ + uv_dict = load_uv(name) + normalized = normalize_uv(**uv_dict) # results in DataArrays with multiple spectra + + return normalized + + +def process_uv(name): + """Processses all UV spectra in a single LoKI.nxs and averages them to one corrected + UV spectrum. + + Parameters + ---------- + name: str + Filename for a LoKI.nxs file containting UV entry. + + Returns + ---------- + normalized: + One averaged UV spectrum. Averaged over all UV spectra contained in the file + under UV entry data. + + """ + + uv_dict = load_uv(name) + normalized = normalize_uv(**uv_dict) + + # returns averaged uv spectrum + return normalized.mean("spectrum") + + +def plot_uv(name): + """Plots multiple graphs of UV data contained in a LoKI.nxs file. + First graph: Raw UV spectra + Second graph: Plot of dark and reference + Third graph: Individual normalized UV spectra + Fourth graph: Averaged UV spectrum over all found data in LoKI.nxs file. + + Parameters + ---------- + name: str + Filename for a LoKI.nxs file containting UV entry. + + """ + # extracting and preparing uv data + uv_dict = load_uv(name) + normalized = normalize_uv(**uv_dict) # results in DataArrays with multiple spectra + + # legend property + legend_props = {"show": True, "loc": 1} + figs, axs = plt.subplots(1, 4, figsize=(16, 3)) + + # How to plot raw data spectra in each file? + out1 = sc.plot( + sc.collapse(uv_dict["sample"], keep="wavelength"), + linestyle="dashed", + marker=".", + grid=True, + legend=legend_props, + title=f"{name} - raw data", + ax=axs[0], + ) + + # How to plot dark and reference? + out4 = sc.plot( + {"dark": uv_dict["dark"], "reference": uv_dict["reference"]}, + linestyle="dashed", + marker=".", + grid=True, + legend=legend_props, + title=f"{name} - dark and reference", + ax=axs[1], + ) + + # How to plot individual calculated spectra from one .nxs file + out2 = sc.plot( + sc.collapse(normalized, keep="wavelength"), + linestyle="dashed", + marker=".", + grid=True, + legend=legend_props, + title=f"{name} - calculated", + ax=axs[2], + ) + + # How to plot averaged spectra (.mean and NOT .sum) + out3 = sc.plot( + sc.collapse(normalized.mean("spectrum"), keep="wavelength"), + linestyle="dashed", + marker=".", + grid=True, + legend=legend_props, + title=f"{name} - averaged", + ax=axs[3], + ) + # display(figs) + + +def gather_uv_set(flist_num): + """Creates a sc.DataSet for set of given filenames for an experiment composed of + multiple, separated UV measurements over time. + Parameters + ---------- + flist_num: list of int + List of filenames as numbers (ILL style) containing UV data + + Returns + ---------- + uv_spectra_set: sc.DataSet + DataSet of multiple UV DataArrays, where the UV signal for each experiment was averaged + + """ + + uv_spectra_set = sc.Dataset({name: process_uv(name) for name in flist_num}) + return uv_spectra_set + + +def plot_uv_set(flist_num, lambda_min=None, lambda_max=None, vmin=None, vmax=None): + """Plots a set of averaged UV spectra + + Parameters + ---------- + flist_num: list of int + List of filenames as numbers (ILL style) containing UV data + lambda_min: float + Minimum wavelength + lambda_max: float + Maximum wavelength + vmin: float + Minimum y-value + vmax: float + Maximum y-value + + Returns + ---------- + fig_all_spectra: scipp.plotting.objects.Plot + + """ + + # creates a dataset based on flist_num, plots all of them + uv_spectra_set = gather_uv_set(flist_num) + # How to plot multiple UV spectra and zoom in. + # set figure size and legend props + figure_size = (7, 4) + legend_props = {"show": True, "loc": (1.04, 0)} + fig_all_spectra = uv_spectra_set.plot( + linestyle="dashed", + marker=".", + grid=True, + legend=legend_props, + figsize=figure_size, + ) + if vmin is not None and vmax is not None: + fig_all_spectra.ax.set_ylim(vmin, vmax) + if lambda_min is not None and lambda_max is not None: + fig_all_spectra.ax.set_xlim(lambda_min, lambda_max) + display(fig_all_spectra) + return fig_all_spectra + + +def load_fluo(name): + """Loads the data contained in the fluo entry of a LoKI.nxs file + + Parameters + ---------- + name: str + Filename, e.g. 066017.nxs + + Returns + ---------- + fluo_dict: dict + Dictionary of sc.DataArrays. Keys: data, reference, dark. Data contains the fluo signals of the sample. + + """ + + with snx.File(name) as f: + fluo = f["entry/instrument/fluorescence"][()] + + # separation + fluo_dict = split_sample_dark_reference(fluo) + + return fluo_dict + + +def normalize_fluo( + *, sample: sc.DataArray, reference: sc.DataArray, dark: sc.DataArray +) -> sc.DataArray: + """Calculates the corrected fluo signal for each given fluo spectrum in a given sc.DataArray + + Parameters + ---------- + sample: sc.DataArray + DataArray containing sample fluo signal, one spectrum or multiple spectra. + reference: sc.DataArray + DataArray containing reference fluo signal, one spectrum expected. + dark: sc.DataArray + DataArray containing dark fluo signal, one spectrum expected. + + Returns + ---------- + final_fluo: sc.DataArray + DataArray that contains the calculated fluo signal, one spectrum or mulitple spectra. + + """ + + # all spectra in this file are converted to final_fluo + # More explanation on the dark here. + # We keep dark here for consistency reasons. + # The dark measurement is necessary for this type of detector. Sometimes, + # one can use the fluorescence emission without the reference. In that case + # having the dark is important. + # In most cases the reference should have no fluorescence emission, + # basically flat. In more complex solvent the reference may have some + # intrinsic fluorescence that would need to be substracted. + + final_fluo = (sample - dark) - (reference - dark) + + return final_fluo + + +def plot_fluo(name): + """Plots all fluo spectra contained in a LoKI.nxs + Currently, we separate between good and bad fluo spectra collected during a ILL experiment. This differentiation can later go for LoKI. + First graph: Plot of all raw fluo spectra + Second graph: Plot dark and reference for the fluo path + Third graph: All bad fluo spectra recorded during the ILL experiment + Fourth graph: All good fluo spectra recorded during the ILL experiment + Fifth graph: All good fluo spectra for wavelength range 250nm - 600nm + + + Parameters + ---------- + name: str + Filename of a LoKI.nxs file + + """ + # Plots fluo spectra in LOKI.nxs of file name + # Input: name, str + + fluo_dict = load_fluo(name) # this is dictionary with scipp.DataArrays in it + final_fluo = normalize_fluo(**fluo_dict) # this is one scipp.DataArray + + # set figure size and legend props + figure_size = (8, 4) + legend_props = {"show": True, "loc": (1.04, 0)} + + # plot all fluo raw spectra + out1 = sc.plot( + sc.collapse(fluo_dict["sample"]["spectrum", :], keep="wavelength"), + linestyle="dashed", + grid=True, + legend=legend_props, + figsize=figure_size, + title=f"{name}, raw fluo spectrum - all", + ) # legend={"show": True, "loc": (1.0, 1.0)} #figsize=(width, height) + display(out1) + + # plot raw and dark spectra for fluo part + out2 = sc.plot( + {"dark": fluo_dict["dark"], "reference": fluo_dict["reference"]}, + linestyle="dashed", + grid=True, + legend=legend_props, + figsize=figure_size, + title=f"{name}, dark and reference", + ) # legend={"show": True, "loc": (1.0, 1.0)} #figsize=(width, height)) + display(out2) + + # specific for ILL data, every second sppectrum good, pay attention to range() where + # the selection takes place + only_bad_spectra = {} # make empty dict + for i in range(0, fluo_dict["sample"].sizes["spectrum"], 2): + only_bad_spectra[f"spectrum-{i}"] = fluo_dict["sample"]["spectrum", i] + out3 = sc.plot( + only_bad_spectra, + linestyle="dashed", + grid=True, + legend=legend_props, + figsize=figure_size, + title=f"{name} - all bad raw spectra", + ) + display(out3) + + only_good_spectra = {} # make empty dict + for i in range(1, fluo_dict["sample"].sizes["spectrum"], 2): + only_good_spectra[f"spectrum-{i}"] = fluo_dict["sample"]["spectrum", i] + out4 = sc.plot( + only_good_spectra, + linestyle="dashed", + grid=True, + legend=legend_props, + figsize=figure_size, + title=f"{name} - all good raw spectra", + ) + display(out4) + + # plot only good spectra with wavelength in legend name + only_good_fspectra = {} # make empty dict + for i in range(1, fluo_dict["sample"].sizes["spectrum"], 2): + mwl = str(final_fluo.coords["monowavelengths"][i].value) + "nm" + only_good_fspectra[f"spectrum-{i}, {mwl}"] = final_fluo["spectrum", i] + out5 = sc.plot( + only_good_fspectra, + linestyle="dashed", + grid=True, + legend=legend_props, + figsize=figure_size, + title=f"{name} - all good final spectra", + ) + display(out5) + + # zoom in final spectra selection + out6 = sc.plot( + only_good_fspectra, + linestyle="dashed", + grid=True, + legend=legend_props, + figsize=figure_size, + title=f"{name} - all good final spectra", + ) + lambda_min = 250 + lambda_max = 600 + out6.ax.set_xlim(lambda_min, lambda_max) + display(out6) + + fig_handles = [out1, out2, out3, out4, out5, out5, out6] + modify_plt_app(fig_handles) + + +def markers(): + """Creates a list of markers for plots""" + return ["+", ".", "o", "*", "1", "2", "x", "P", "v", "X", "^", "d"] + + +def line_colors(number_of_lines): + """ + Creates a number of colors based on the number of spectra/lines to plot. + + Parameters + ---------- + number_of_lines: int + number of available spectra in a plot + + + Returns + ---------- + colors: list of tuples + """ + + # create better line colors + start = 0.0 + stop = 1.0 + number_of_lines = number_of_lines + cm_subsection = np.linspace(start, stop, number_of_lines) + colors = [cm.jet(x) for x in cm_subsection] + return colors + + +def modify_plt_app(fig_handles): + """Modifies scipp plots. Nicer markers and rainbow colors for multiple curves in a scipp plot object. + + Parameters + ---------- + fig_handles: list of scipp.plotting.objects.Plot + + """ + # Modify scipp plots afterwards + for pl1 in fig_handles: + pl1.fig.tight_layout() + number_lines = len(pl1.ax.get_lines()) + colors = line_colors(number_lines) + for i in range(0, number_lines): + pl1.ax.get_lines()[i].set_color(colors[i]) + pl1.ax.get_legend().legendHandles[i].set_color(colors[i]) + pl1.ax.get_lines()[i].set_marker(markers()[i]) + pl1.ax.get_lines()[i].set_markersize(5) + pl1.ax.get_lines()[i].set_markevery(5) + pl1.ax.get_legend().legendHandles[i].set_marker(markers()[i]) + pl1.ax.get_lines()[i].set_linewidth(1) + + +def export_uv(name, path_output): + """Export normalized all uv data and an averaged uv spectrum in a LoKI.nxs file to .dat file + + Attention: Current output format follows custom format for an individual user, not + any other software. + + Parameters + ---------- + name: str + Filename of LoKI.nxs file that contains the UV data for export + + path_output: str + Absolute path to output folder + + Returns + ---------- + Tab-limited .dat file with columns wavelength, dark, reference, (multiple) raw uv spectra, (multiple) corrected uv spectra, one averaged uv spectrum + + """ + + uv_dict = load_uv(name) + normalized = normalize_uv( + **uv_dict + ) + + normalized_avg = normalized.mean("spectrum") + + # prepare for export as .dat files + output_filename = f"{name}_uv.dat" + + # puzzle the header together + l = "".join( + ["dark_{0}\t".format(i) for i in range(uv_dict['dark'].ndim) + ] + ) + m="".join( + ["reference_{0}\t".format(i) for i in range(uv_dict['reference'].ndim) + ] + ) + + n="".join( + [ + "uv_raw_spectra_{0}\t".format(i) + for i, x in enumerate(range(uv_dict['sample'].sizes["spectrum"])) + ] + ) + o="".join( + [ + "uv_norm_spectra_{0}\t".format(i) + for i, x in enumerate(range(normalized.sizes["spectrum"])) + ] + ) + + p = "uv_spectra_avg\t" + + + hdrtxt = "wavelength [nm]\t" + final_header = hdrtxt + l + m + n + o + p + + + data_to_save = np.column_stack( + ( + normalized.coords["wavelength"].values.transpose(), + + # raw data + uv_dict['dark'].values.transpose(), + uv_dict['reference'].values.transpose(), + uv_dict['sample'].values.transpose(), + + # reduced data + normalized.values.transpose(), + normalized_avg.values.transpose(), + + ) + ) + path_to_save = os.path.join(path_output, output_filename) + + # dump the content + with open(path_to_save, "w") as f: + np.savetxt(f, data_to_save, fmt="%.10f", delimiter="\t", header=final_header) + + +def export_fluo(name, path_output): + """Export corrected fluo data contained in a LoKI.nxs file to .dat file + + Attention: Current output format follows custom format for an individual user, not + any other software. + + Parameters + ---------- + name: str + Filename of LoKI.nxs file that contains the fluo data for export + + path_output: str + Absolute path to output folder + + Returns + ---------- + Tab-limited .dat file with columns wavelength, dark, reference, multiple raw fluo + spectra, and multiple normalized fluo spectra. + Header of each fluo spectrum column contains the incident excitation energy. + + """ + # export of all calculated fluo data in a LOKI.nxs name to .dat + # input: filename of LOKI.nxs: name, str, path_output: absolut path to output folder, str + + fluo_dict = load_fluo(name) + final_fluo = normalize_fluo(**fluo_dict) + + # prepare for export as .dat files + output_filename = f"{name}_fluo.dat" + path_to_save = os.path.join(path_output, output_filename) + + + + l = "".join( + ["dark_{0}\t".format(i) for i in range(fluo_dict['dark'].ndim) + ] + ) + m = "".join( + ["reference_{0}\t".format(i) for i in range(fluo_dict['reference'].ndim) + ] + ) + + n= "".join([f"raw_{i}nm\t" for i in final_fluo.coords["monowavelengths"].values]) + + o= "".join([f"norm_{i}nm\t" for i in final_fluo.coords["monowavelengths"].values]) + + + + hdrtxt = "wavelength [nm]\t" + final_header = hdrtxt + l + m + n + o + + data_to_save = np.column_stack( + ( + final_fluo.coords["wavelength"].values.transpose(), + # dark + fluo_dict['dark'].values.transpose(), + # reference + fluo_dict['reference'].values.transpose(), + # sample + fluo_dict['sample'].values.transpose(), + # final fluo spectra + final_fluo.data["spectrum", :].values.transpose(), + ) + ) + + # dump the content + with open(path_to_save, "w") as f: + np.savetxt( + f, + data_to_save, + fmt="".join(["%.5f\t"] + ["%.5e\t"] * (fluo_dict['dark'].ndim + + fluo_dict['reference'].ndim + 2*final_fluo.sizes["spectrum"])), + delimiter="\t", + header=final_header, + ) + + +def fluo_maxint_max_wavelen( + flist_num, + wllim=300, + wulim=400, + wl_unit=None, + medfilter=True, + kernel_size=15, +): + """For a given list of files this function extracts for each file the maximum intensity and the corresponding + wavelength of each fluo spectrum. + + Parameters + ---------- + flist_num: list + List of filename for a LoKI.nxs file containting Fluo entry. + wllim: float + Lower wavelength limit where the search for the maximum should begin + wulim: float + Upper wavelength limit where the search for the maximum should end + wl_unit: sc.Unit + Unit of the wavelength + medfilter: bool + A medfilter is applied to the fluo spectra as fluo is often more noisy + kernel_size: int, uneven, min>=3 + The width of the median filter, uneven, and at least a value of 3. + + Returns + ---------- + fluo_int_dict: dict + A dictionary of nested dictionaries. For each found monowavelength, there are nested dictionaries + for each file containing the maximum intensity "intensity_max" and the corresponding wavelength "wavelength_max" + + """ + + from collections import defaultdict + + fluo_int_dict = defaultdict(dict) + + for name in flist_num: + fluo_dict = load_fluo(name) + fluo_da = normalize_fluo(**fluo_dict) + # check for the unit + if wl_unit is None: + wl_unit = fluo_da.coords["wavelength"].unit + + + print(f'Number of fluo spectra in {name}: {fluo_da.sizes["spectrum"]}') + print(f"This is the fluo dataarray for {name}.") + display(fluo_da) + + fluo_filt_max = fluo_peak_int( + fluo_da, + wllim=wllim, + wulim=wulim, + wl_unit=wl_unit, + medfilter=medfilter, + kernel_size=kernel_size, + ) + print(f"This is the fluo filt max intensity and wavelength dataset for {name}.") + display(fluo_filt_max) + + # print(f'{name} Unique mwl:', np.unique(fluo_filt_max.coords['monowavelengths'].values)) + unique_monowavelen = np.unique(fluo_filt_max.coords["monowavelengths"].values) + + # create entries into dict + for mwl in unique_monowavelen: + fluo_int_max = fluo_filt_max["intensity_max"][ + fluo_filt_max.coords["monowavelengths"] == mwl * wl_unit + ].values + fluo_wavelen_max = fluo_filt_max["wavelength_max"][ + fluo_filt_max.coords["monowavelengths"] == mwl * wl_unit + ].values + + # I collect the values in a dict with nested dicts, separated by wavelength + fluo_int_dict[f"{mwl}{wl_unit}"][f"{name}"] = { + "intensity_max": fluo_int_max, + "wavelength_max": fluo_wavelen_max, + } + + return fluo_int_dict + + +def fluo_plot_maxint_max_wavelen(fluo_int_dict: dict): + """Plots for fluorescence maximum intensity and corressponding wavelength as function of pertubation. + Pertubation parameter: currently filename as unique identifier + + Parameters + ---------- + fluo_int_dict: dict + Contains per found monowavelength and per measurement the value of the max fluo intensity for a selected wavelength range and the corresponding maximum wavelength. + + Returns + ---------- + + """ + + # separate information + max_int_dict = {} + max_int_wavelen_dict = {} + + monowavelen = fluo_int_dict.keys() + + # prepare the plots + monowavelen = len(fluo_int_dict.keys()) + print("Number of keys", monowavelen) + + fig, ax = plt.subplots(nrows=monowavelen, ncols=2, figsize=(10, 7)) + plt.subplots_adjust( + left=None, bottom=0.05, right=None, top=0.95, wspace=0.2, hspace=0.5 + ) + + for j, mwl in enumerate(fluo_int_dict): + + for name in fluo_int_dict[mwl]: + max_int_dict[name] = fluo_int_dict[mwl][name]["intensity_max"] + max_int_wavelen_dict[name] = fluo_int_dict[mwl][name]["wavelength_max"] + + for i, fname in enumerate(max_int_dict): + + x_val = [i] * len(max_int_dict[fname]) + ax[j, 0].set_title(f"monowavelength: {mwl}: max. intensity") + ax[j, 0].plot(x_val, max_int_dict[fname], marker="o", linestyle="None") + ax[j, 0].grid(True) + ax[j, 0].set_ylabel("Max. peak intensity [counts]") + print(mwl, fname, "max int", max_int_dict[fname]) + + ax[j, 1].set_title(f"monowavelength: {mwl}: wavelength") + ax[j, 1].plot( + x_val, max_int_wavelen_dict[fname], marker="o", linestyle="None" + ) + ax[j, 1].set_ylabel("Peak wavelength [nm]") + ax[j, 1].grid(True) + # ax[1].set_ylabel(f'Wavelength [{unit_str}]') + # ax[j,1].set_ylim(bottom=0.95*int(mwl[:-4])) + + ax[j, 1].plot( + np.linspace( + 0, + len(max_int_wavelen_dict.keys()) - 1, + num=len(max_int_wavelen_dict.keys()), + ), + [int(mwl[:-4])] * len(max_int_wavelen_dict.keys()), + "--", + ) + # overwrite xticks in ax0 + ax[j, 0].set_xticks( + np.arange(0, len(max_int_dict.keys()), 1), + labels=[f"{name}" for name in max_int_dict], + rotation=90, + ) + # overwrite xticks in ax1 + ax[j, 1].set_xticks( + np.arange(0, len(max_int_wavelen_dict.keys()), 1), + labels=[f"{name}" for name in max_int_wavelen_dict], + rotation=90, + ) + + +def apply_medfilter( + da: sc.DataArray, kernel_size: Optional[int] = None +) -> sc.DataArray: + #TODO: Rewrite this function with median filters offered by scipp. + """Applies a median filter to a sc.DataArray that contains fluo or uv spectra to remove spikes + Filter used: from scipy.ndimage import median_filter. This filter function is maybe subject to change for another scipy function. + Caution: The scipy mean (median?) filter will just ignore errorbars on counts according to Neil Vaytet. + + Parameters + ---------- + da: sc.DataArray + DataArray that contains already fluo or uv spectra in data. + + kernel_size: int + Scalar giving the size of the median filter window. Elements of kernel_size should be odd. Default size is 3 + + Returns + ---------- + da_mfilt: sc.DataArray + New sc.DataArray with median filtered data + + """ + if kernel_size is None: + kernel_size = 3 + else: + assert type(kernel_size) is int, "kernel_size must be an integer" + + # check if kernel_size is odd + if (kernel_size % 2) != 1: + raise ValueError("kernel_size should be odd.") + + # extract all spectreal values + data_spectrum = da.values + + # apply a median filter + # yes, we can apply a median filter to an array (vectorize it), but the kernel needs to be adapted to the ndim off the array + # no filterin in the first dimension (spectrum), but in the second (wavelength) + # magic suggested by Simon + # this will make [kernel_size] for 1d data, and [1, kernel_size] for 2d data + ksize = np.ones_like(da.shape) + ksize[-1] = kernel_size + data_spectrum_filt = median_filter(data_spectrum, size=ksize) + + # create a new sc.DataArray where the data is replaced by the median filtered data + # make a deep copy + da_mfilt = da.copy() + + # replace original data with filtered data + da_mfilt.values = data_spectrum_filt + + # graphical comparison + legend_props = {"show": True, "loc": (1.04, 0)} + figs, axs = plt.subplots(1, 2, figsize=(16, 5)) + # out1=sc.plot(da['spectrum',7], ax=axs[0]) + # out2=sc.plot(da_mfilt['spectrum',7], ax=axs[1]) + out1 = sc.plot( + sc.collapse(da, keep="wavelength"), + linestyle="dashed", + marker=".", + grid=True, + legend=legend_props, + title=f"Before median filter", + ax=axs[0], + ) + out2 = sc.plot( + sc.collapse(da_mfilt, keep="wavelength"), + linestyle="dashed", + marker=".", + grid=True, + legend=legend_props, + title=f"After median filter - kernel_size: {kernel_size}", + ax=axs[1], + ) + # display(figs) + + return da_mfilt + + +def uv_peak_int(uv_da: sc.DataArray, wavelength=None, wl_unit=None, tol=None): + """Extract uv peak intensity for a given wavelength [wl_unit] and a given interval [wl_unit]. + First version: interval around wavelength of width 2*tol, values are averaged and then we need interpolation to get value at requested wavelength??? not sure yet how to realise this + If no wavelength is given, 280 is chosen as value + If no wl_unit is given, the unit of the wavelength coordinate of the sc.DataArray is chosen, e.g. [nm]. Other option to generate a unit: wl_unit=sc.Unit('nm') + + Parameters + ---------- + uv_da: sc.DataArray + DataArray containing uv spectra + wavelength: float + Wavelength + wl_unit: sc.Unit + Unit of the wavelength + tol: float + Tolerance, 2*tol defines the interval around the given wavelength + + Returns + ---------- + uv_peak_int: dict + Dictionary that contains the peak intensity for the requested wavelength, the peak intensity averaged over the requested interval, the requested wavelength with its unit, and the tolerance + + """ + assert ( + "wavelength" in uv_da.dims + ), "sc.DataArray is missing the wavelength dimension" # assert that 'wavelength' is a dimension in the uv_da sc.DataArray + + # obtain unit of wavelength: + if wl_unit is None: + wl_unit = uv_da.coords["wavelength"].unit + else: + if not isinstance(wl_unit, sc.Unit): + raise TypeError + assert ( + wl_unit == uv_da.coords["wavelength"].unit + ) # we check that the given unit corresponds to the unit for the wavelength + + # set default value for wavelength: + if wavelength is None: + wavelength = 280 + + # set default value for tolerance: + if tol is None: + tol = 0.5 + + # filter spectrum values for the specified interval, filtered along the wavelength dimension + uv_da_filt = uv_da[ + "wavelength", + (wavelength - tol) * wl_unit : (wavelength + tol) * wl_unit, + ] + + # try out interpolation + from scipp.interpolate import interp1d + + uv_interp = interp1d(uv_da, "wavelength") + # get new x values, in particular I want the value at one wavelength + x2 = sc.linspace( + dim="wavelength", start=wavelength, stop=wavelength, num=1, unit=wl_unit + ) + # alternative + # x2 = sc.array(dims=['wavelength'], values=[wavelength], unit=wl_unit) + # this gives us the intensity value for one specific wavelength + uv_int_one_wl = uv_interp(x2) + + # now we want to have as well the wavelength for the interval + uv_int_mean_interval = uv_da_filt.mean(dim="wavelength") + + # prepare a dict for output + uv_peak_int = { + "one_wavelength": uv_int_one_wl, + "wl_interval": uv_int_mean_interval, + "wavelength": wavelength, + "unit": wl_unit, + "tol": tol, + } + + return uv_peak_int + + +def plot_uv_peak_int( + uv_da: sc.DataArray, + name, + wavelength=None, + wl_unit=None, + tol=None, + medfilter=False, + kernel_size=None, +): + """Plotting of extracted uv peak intensity for a given wavelength [wl_unit] and a given interval [wl_unit] in file name + First version: interval around wavelength of width 2*tol, values are averaged and then we need interpolation to get value at requested wavelength??? not sure yet how to realise this + If no wavelength is given, 280 is chosen as value + If no wl_unit is given, the unit of the wavelength coordinate of the sc.DataArray is chosen, e.g. [nm]. Other option to generate a unit: wl_unit=sc.Unit('nm') + Kernel size determined the window for the medfilter + + Parameters + ---------- + uv_da: sc.DataArray + DataArray containing uv spectra + name: str + Filename of a file containting uv data + wavelength: float + Wavelength + wl_unit: sc.Unit + Unit of the wavelength + tol: float + Tolerance, 2*tol defines the interval around the given wavelength + medfilter: boolean + If medfilter==False, not medfilter is applied + kernel_size: int + kernel for medianfilter + + Returns + ---------- + uv_peak_int: dict + Dictionary that contains the peak intensity for the requested wavelength, the peak intensity averaged over the requested interval, the requested wavelength with its unit, and the tolerance + + """ + if not isinstance(uv_da, sc.DataArray): + raise TypeError + + if medfilter is not False: + # apply medianfilter + uv_da = apply_medfilter(uv_da, kernel_size=kernel_size) + + # process + uv_peak_int_dict = uv_peak_int( + uv_da, wavelength=wavelength, wl_unit=wl_unit, tol=tol + ) + # append name to dict + uv_peak_int_dict["filename"] = name + + print(uv_peak_int_dict["one_wavelength"]) + print(uv_peak_int_dict["wl_interval"]["spectrum", :]) + + # set figure size and legend props + # figs, axs = plt.subplots(1, 2, figsize=(10 ,3)) + fig = plt.figure(figsize=(12, 4)) + + gs = gridspec.GridSpec(nrows=1, ncols=2, width_ratios=[1, 1], wspace=0.4) + ax0 = fig.add_subplot(gs[0, 0]) + out1 = sc.plot( + uv_peak_int_dict["one_wavelength"]["spectrum", :].squeeze(), + linestyle="none", + grid=True, + title=f"UV peak intensity for {uv_peak_int_dict['wavelength']} {uv_peak_int_dict['unit']} in {name} ", + ax=ax0, + ) + ax0.set_ylabel("UV int") + ax1 = fig.add_subplot(gs[0, 1]) + out2 = sc.plot( + uv_peak_int_dict["wl_interval"]["spectrum", :], + linestyle="none", + grid=True, + title=f"UV peak intensity for intervall [{uv_peak_int_dict['wavelength']-uv_peak_int_dict['tol']} , {uv_peak_int_dict['wavelength']+uv_peak_int_dict['tol']}]{uv_peak_int_dict['unit']} in {name} ", + ax=ax1 + ) + ax1.set_ylabel("UV int") + + +def uv_quick_data_check( + filelist, + wavelength=None, + wl_unit=None, + tol=None, + medfilter=False, + kernel_size=None, +): + """Plots uv peak intensity for a given wavelength with unit wl_unit for a given + filelist in order of items in filelist. Per recorded spectrum in the file a data + point is shown. + + """ + # Currently I cannot rely on having the same number of spectra in each file currently, sc.Dataset will not work + # The following does not work if the number of spectra in each Loki file is not the same :() + # out= sc.Dataset({name:load_and_normalize_uv(name) for name in filelist}) + + # Let's try dictionary comprehension ... yeap that works + uv_dict = {name: load_and_normalize_uv(name) for name in filelist} + + # apply medianfilter + if medfilter is not False: + # apply_medfilter expects a sc.DataArray + uv_dict = { + name: apply_medfilter(uv_dict[name], kernel_size=kernel_size) + for name in filelist + } + + # extract peak intensities + # function: uv_peak_int(uv_da: sc.DataArray, wavelength=wavelength, wl_unit=wl_unit, tol=tol) + # this process will return a dictionary of dictionaries + uv_peak_int_dict = { + name: uv_peak_int( + uv_dict[name], wavelength=wavelength, wl_unit=wl_unit, tol=tol + ) + for name in filelist + } + + # prepare plots + figure_size = (20, 10) + fig = plt.figure(figsize=figure_size) + gs = gridspec.GridSpec( + nrows=2, + ncols=3, + width_ratios=[1, 1, 1], + wspace=0.3, + height_ratios=[1, 1], + ) + ax0 = fig.add_subplot(gs[0, 0]) + ax1 = fig.add_subplot(gs[0, 1]) + + # Testing but so far unsuccessful + # print(uv_peak_int_dict['066017.nxs']['one_wavelength'].data) + # display(sc.plot(sc.collapse(uv_peak_int_dict['066017.nxs']['one_wavelength'].data,keep='spectrum'))) + + uv_peak_one_wl = {} # make empty dict + for name in filelist: + # uv_peak_one_wl[f'{name}']=uv_peak_int_dict[f'{name}']['one_wavelength']['spectrum',:] + # print(uv_peak_one_wl[f'{name}']) + + ax0.plot( + uv_peak_int_dict[f"{name}"]["one_wavelength"]["spectrum", :].values, + "o", + label=f"{name}", + ) + ax0.legend() + + ax1.plot( + uv_peak_int_dict[f"{name}"]["wl_interval"]["spectrum", :].values, + "o", + label=f"{name}", + ) + ax1.legend() + + ax0.grid(visible=True) + ax0.set_title( + f"UV peak intensity for {uv_peak_int_dict[filelist[0]]['wavelength']} {uv_peak_int_dict[filelist[0]]['unit']} " + ) + ax0.set_xlabel("spectrum") + + ax1.plot( + uv_peak_int_dict[filelist[0]]["wl_interval"]["spectrum", :].values, + "o", + label=f"{name}", + ) + ax1.legend() + ax1.grid(visible=True) + ax1.set_title( + f"UV peak intensity for intervall [{uv_peak_int_dict[filelist[0]]['wavelength']-uv_peak_int_dict[filelist[0]]['tol']} , {uv_peak_int_dict[filelist[0]]['wavelength']+uv_peak_int_dict[filelist[0]]['tol']}]{uv_peak_int_dict[filelist[0]]['unit']}" + ) + print( + uv_peak_int_dict[filelist[0]]["wl_interval"]["spectrum", :].values, + uv_peak_int_dict[filelist[0]]["wl_interval"]["spectrum", :].values, + ) + ax1.set_xlabel("spectrum") + + display(fig) + + +def plot_multiple_uv_peak_int( + filelist, + wavelength=None, + wl_unit=None, + tol=None, + medfilter=False, + kernel_size=None, +): + """Plots uv peak intensity for a given wavelength with unit wl_unit for a given filelist in order of items in filelist + Explorative state. I fall back to matplotlib because it is complicated. Try later with sc.plot. + Currently, I cannot rely on all files having the same number of spectra. I cannot create a sc.Dataset. I am unsure how to catch this with Python. + No averaging in uv spectra when calculating the uv peak intensity within one file. + Medfilter can be activated + + #TODO: Get rid of the for loop and and make filename an attribute when loading nexus files, filename as new dimension? + """ + + # currently I cannot rely on having the same number of spectra in each file currently, sc.Dataset will not work + fig = plt.figure(figsize=(15, 5)) + gs = gridspec.GridSpec( + nrows=1, ncols=2, width_ratios=[1, 1], wspace=0.3, height_ratios=[1] + ) + ax0 = fig.add_subplot(gs[0, 0]) + ax1 = fig.add_subplot(gs[0, 1]) + + for i, name in enumerate(filelist): + uv_dict = load_uv(name) + uv_da = normalize_uv(**uv_dict) + + # check for medfilter + if medfilter is not False: + # apply medianfilter + uv_da = apply_medfilter(uv_da, kernel_size=kernel_size) + + # process + uv_peak_int_dict = uv_peak_int( + uv_da, wavelength=wavelength, wl_unit=wl_unit, tol=tol + ) + # append name to dict + uv_peak_int_dict["filename"] = name + + num_spec = uv_peak_int_dict["wl_interval"]["spectrum", :].values.size + + # prints number of spectra in each file and artifical x-value + # print(num_spec,i*np.ones(num_spec)) + + # to offset each spectrum on the x-axis, I create an artifical x-axis at the moment, TODO: could be later replaced by the pertubation parmeter (time, concentration) + ax0.plot( + i * np.ones(num_spec), + uv_peak_int_dict["one_wavelength"]["spectrum", :].values, + "o", + label=f"{name}", + ) + ax0.legend() + ax0.grid(visible=True) + ax0.set_title( + f"UV peak intensity for {uv_peak_int_dict['wavelength']} {uv_peak_int_dict['unit']} " + ) + + ax1.plot( + i * np.ones(num_spec), + uv_peak_int_dict["wl_interval"]["spectrum", :].values, + "o", + label=f"{name}", + ) + ax1.legend() + ax1.grid(visible=True) + ax1.set_title( + f"UV peak intensity for intervall [{uv_peak_int_dict['wavelength']-uv_peak_int_dict['tol']} , {uv_peak_int_dict['wavelength']+uv_peak_int_dict['tol']}]{uv_peak_int_dict['unit']}" + ) + + # overwrite xticks in ax0 + ax0.set_xticks( + np.arange(0, len(filelist), 1), + labels=[f"{name}" for name in filelist], + rotation=90, + ) + # overwrite xticks in ax1 + ax1.set_xticks( + np.arange(0, len(filelist), 1), + labels=[f"{name}" for name in filelist], + rotation=90, + ) + + display(fig) + + +def turbidity(wl, b, m): + """Function describing turbidity tau. tau = b* lambda **(-m) + Fitting parameters: b, m. b corresponds to the baseline found for higher wavelengths (flat line in UV spectrum), m corresponds to the slope. + lambda: wavelength + + Parameters + ---------- + b: np.ndarray + Offset, baseline + + m: np.ndarray + Slope + + wl: np.ndarray + UV wavelengths + + Returns + ---------- + y: np.ndarray + Turbidity + + """ + + y = b * wl ** (-m) + return y + + +def residual(p, x, y): + """Calculates the residuals between fitted turbidity and measured UV data + + Parameters + ---------- + p: list + Fit parameters for turbidity + + x: np.ndarray + x values, here: UV wavelength + + y: np.ndarray + y values, here: UV intensity + + Returns + ---------- + y - turbidity(x, *p): np.ndarray + Difference between measured UV intensity values and fitted turbidity + + """ + + return y - turbidity(x, *p) + + +def uv_turbidity_fit( + uv_da: sc.DataArray, + wl_unit=None, + fit_llim=None, + fit_ulim=None, + b_llim=None, + b_ulim=None, + m=None, + plot_corrections=True, +): + """Fit turbidity to the experimental data. Turbidity: tau=b * wavelength^(-m) Parameters of interest: b, m. + b is the baseline and m is the slope. b can be obtained by averaging over the flat range of the UV spectrum + in the higher wavelength range. + m: make an educated guess. Advice: Limit fitting range to wavelengths after spectroscopy peaks. + + Parameters + ---------- + uv_da: sc.DataArray + UV sc.DataArray containing one or more normalized UV spectra + wl_unit: sc.Unit + Wavelength unit + fit_llim: int + Lower wavelength limit of fit range for turbidity + fit_ulim: int + Upper wavelength limit of fit range for turbidity + b_llim: int + Lower wavelength limit of fit range for b + b_ulim: int + Upper wavelength limit of fit range for b + m: int + Educated guess start value of slope parameter in turbidity + plot_corrections: bool + If true, plot single contribitions for turbidity corrections + + Returns: + ---------- + uv_da_turbcorr: sc.DataArray + uv_da dataarray where each spectrum was corrected for a fitted turbidity, export for all wavelength values + + """ + # obtain unit of wavelength: + if wl_unit is None: + wl_unit = uv_da.coords["wavelength"].unit + else: + if not isinstance(wl_unit, sc.Unit): + raise TypeError + assert ( + wl_unit == uv_da.coords["wavelength"].unit + ) # we check that the given unit corresponds to the unit for the wavelength + + if fit_llim is None: + fit_llim = 400 + if fit_ulim is None: + fit_ulim = 800 + + if fit_llim is not None and fit_ulim is not None: + assert fit_llim < fit_ulim, "fit_llim < fit_ulim" + if b_llim is not None and b_ulim is not None: + assert b_llim < b_ulim, "b_llim dict: + """This function extracts a specific fluo spectrum from all given files. Ideally, the user provides the number index. + A median filter can be applied to the input data. A lower and upper wavelength range will be used for a zoomed-in image. + Selected spectra are all plotted in one graph. + + Parameters + ---------- + flist_num: list + List of LoKI.nxs file containting fluo entry. + spectral_idx: int + Index of spectrum for selection + kernel_size: int + Scalar giving the size of the median filter window. Elements of kernel_size should be odd. Default size is 3 + wllim: float + Lower wavelength limit + wulim: float + Upper wavelength limit + wl_unit: sc.Unit + Wavelength unit of a fluo spectum + + Returns + ---------- + fluo_spec_idx_ds: sc.Dataset + sc.Dataset containing selected spectra for each input LoKI.nxs + + """ + # This is maybe not necessary, because scipp catches it later, but I can catch it earlier. + if not isinstance(spectral_idx, int): + raise TypeError("Spectral index should be of type int.") + + if wllim is None: + wllim = 300 + if wulim is None: + wulim = 400 + + # obtain unit of wavelength, I extract it from the first element in the flist_num, + # but I have to load it + fluo_dict = load_fluo(flist_num[0]) + final_fluo = normalize_fluo(**fluo_dict) + if wl_unit is None: + wl_unit = final_fluo.coords["wavelength"].unit + else: + if not isinstance(wl_unit, sc.Unit): + raise TypeError("wl_unit should be of type sc.Unit.") + assert ( + wl_unit == final_fluo.coords["wavelength"].unit + ) # we check that the given unit corresponds to the unit for the wavelength + + # prepare data + fluo_spec_idx_ds = sc.Dataset() + for name in flist_num: + fluo_dict = load_fluo(name) # this is dictionary with scipp.DataArrays in it + fluo_normalized = normalize_fluo(**fluo_dict) # this is one scipp.DataArray + final_fluo = apply_medfilter( + fluo_normalized, kernel_size=kernel_size + ) # apply medfilter + fluo_spec_idx_ds[name] = final_fluo[ + "spectrum", spectral_idx + ] # append data array to dataset + + # Plotting, TODO: should it be decoupled form the loading? + legend_props = {"show": True, "loc": (1.04, 0)} + fig = plt.figure(figsize=(17, 5)) + gs = gridspec.GridSpec(nrows=1, ncols=2, width_ratios=[1, 1], wspace=0.5) + ax0 = fig.add_subplot(gs[0, 0]) + ax1 = fig.add_subplot(gs[0, 1]) + out0 = sc.plot( + fluo_spec_idx_ds, + grid=True, + linestyle="dashed", + title=f"Normalized fluo spectra, selected spectrum #{spectral_idx}", + legend=legend_props, + ax=ax0, + ) # would also work with fluo_spec_idx_dict + out1 = sc.plot( + fluo_spec_idx_ds["wavelength", wllim * wl_unit : wulim * wl_unit], + grid=True, + linestyle="dashed", + title=f"Normalized fluo spectra, selected spectrum #{spectral_idx}", + legend=legend_props, + ax=ax1, + ) + display(fig) + + return fluo_spec_idx_ds diff --git a/docs/instruments/loki/nurf/scipp_947.ipynb b/docs/instruments/loki/nurf/scipp_947.ipynb new file mode 100644 index 000000000..a268c2bd9 --- /dev/null +++ b/docs/instruments/loki/nurf/scipp_947.ipynb @@ -0,0 +1,879 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Demonstration of NUrF function module\n", + "\n", + "The Nurf Python module contains functions to extract UV and fluo data from a LoKI.nxs file.
\n", + "It is capable of correcting the measured spectra for reference and dark current. It calculates the final fluo and UV spectra and plots them. \n", + "A median filter can be applied to the spectra.
\n", + "For UV the following parameters can be extracted:
\n", + "- Absorbance at 280nm\n", + "- Turbidity and m, b factor\n", + "\n", + "For fluo the following parameters can be extracted:
\n", + "- Maximum peak intensity\n", + "- Peak position corresponding to the maximum peak intensity\n", + "\n", + "Besides standard Python libraries like matplotlib, os, scipy, and numpy, it is relies on scipp, scippnexus, and scippneutron. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Prepare the stage and load the required tools." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "from nurf import *\n", + "from ill_auxilliary_funcs import *\n", + "from scipp.signal import butter, sosfiltfilt" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The code below uses UV and fluo data that was colleced in a multi-modal SANS, UV, Fluo experiment at ILL. One experiment is composed of a series of measurement. In each measurement, UV and fluo spectra were recorded. Alongside the sample, dark and reference measurements were recorded. Here, exp5 and exp6, exp2 and exp3, exp7 and exp8 form complementary datasets. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Provide information on where the relevant files are and which files belong to which experimental datasets." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "# Prepare for export to .dat for uv and fluo\n", + "\n", + "# path to LOKI-like files\n", + "process_folder='/Users/gudlo523/Library/CloudStorage/OneDrive-LundUniversity/UU/ILL_947_June_2021/212/d22/exp_9-13-947/processed/ess_version'\n", + "# change to folder\n", + "os.chdir(process_folder)\n", + "\n", + "# export path for .dat files\n", + "path_output='/Users/gudlo523/Library/CloudStorage/OneDrive-LundUniversity/UU/ILL_947_June_2021/212/d22/exp_9-13-947/processed/ess_version/dat-files'\n", + "\n", + "# experimental data sets\n", + "exp5= [66017, 66020, 66023, 66026, 66029, 66032, 66034, 66037, 66040, 66043, 66046]\n", + "exp6= [66050, 66053, 66056, 66059, 66062, 66065, 66068, 66071, 66074, 66077, 66080]\n", + "\n", + "exp2= [65925, 65927, 65930, 65933, 65936, 65939, 65942, 65945, 65948, 65951, 65954, 65957]\n", + "exp3= [65962, 65965, 65968, 65971, 65974, 65977, 65980, 65983, 65986, 65989, 65992]\n", + "\n", + "exp7= [66083, 66086, 66089, 66092, 66095, 66098, 66101, 66104, 66107, 66110, 66113]\n", + "exp8= [66116, 66119, 66122, 66125, 66128, 66131, 66134, 66137, 66140, 66143, 66146]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The cell below shows how to export from LoKI.nxs to a generic Ascii file .dat if the users would like to use .dat files." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# List comprehensions for export to .dat\n", + "all_data=[exp2, exp3, exp5, exp6, exp7, exp8]\n", + "\n", + "# Create one list only with all entries of the exp sets\n", + "all_data_flat = sum(all_data, [])\n", + "\n", + "# convert numbers to file names\n", + "all_file_names = [complete_fname(i) for i in all_data_flat]\n", + "\n", + "# export uv data for all_data to .dat\n", + "[export_uv(m, path_output) for m in all_file_names]\n", + "\n", + "#export fluo data for all_data to .dat\n", + "[export_fluo(m, path_output) for m in all_file_names]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Plotting of UV spectra for exp5 and exp6 " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Plot two uv sets\n", + "\n", + "# path to LOKI-like files\n", + "process_folder='/Users/gudlo523/Library/CloudStorage/OneDrive-LundUniversity/UU/ILL_947_June_2021/212/d22/exp_9-13-947/processed/ess_version'\n", + "os.chdir(process_folder)\n", + "\n", + "# exp 5 and exp 6\n", + "exp5=[66017, 66020, 66023, 66026, 66029, 66032, 66034, 66037, 66040, 66043, 66046]\n", + "exp6= [66050, 66053, 66056, 66059, 66062, 66065, 66068, 66071, 66074, 66077, 66080]\n", + "\n", + "# Create filenames \n", + "flist_num=complete_fname(exp5)\n", + "\n", + "# Plotting exp5\n", + "main_title='Exp 5'\n", + "pl1=plot_uv_set(flist_num, lambda_min=None, lambda_max=None, vmin=None, vmax=None)\n", + "pl1.ax.set_title(main_title)\n", + "\n", + "pl2=plot_uv_set(flist_num, lambda_min=200,lambda_max=400, vmin=-1, vmax=6)\n", + "pl2.ax.set_title(main_title)\n", + "#modify plots\n", + "fig_handles=[pl1, pl2]\n", + "modify_plt_app(fig_handles)\n", + "\n", + "#Example to plot one UV file. \n", + "plot_uv('066017.nxs')\n", + "\n", + "# We can gather all normalized uv data in a sc.Dataset\n", + "uv_exp5_set=gather_uv_set(flist_num)\n", + "\n", + "\n", + "# Repeating all steps above for exp6.\n", + "flist_num=complete_fname(exp6)\n", + "main_title='Exp 6'\n", + "pl1=plot_uv_set(flist_num, lambda_min=None, lambda_max=None, vmin=None, vmax=None)\n", + "pl1.ax.set_title(main_title)\n", + "plt.tight_layout()\n", + "\n", + "pl2=plot_uv_set(flist_num, lambda_min=200,lambda_max=400, vmin=-1, vmax=6)\n", + "pl2.ax.set_title(main_title)\n", + "plt.tight_layout()\n", + "#modify plots\n", + "fig_handles=[pl1, pl2]\n", + "modify_plt_app(fig_handles)\n", + "uv_exp6_set=gather_uv_set(flist_num)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Back to some basic code functionality \n", + "This cell shows how to load a Loki.nxs files, what is returned and how to separate the components in the .nxs file (sample, darf, reference)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "name='066017.nxs'\n", + "with snx.File(name) as f:\n", + " uv = f[\"entry/instrument/uv\"][()]\n", + "\n", + "print('Output of .nxs is of type ', type(uv))\n", + "display(uv)\n", + "da=uv\n", + "dark = da[da.coords[\"is_dark\"]].squeeze()\n", + "ref = da[da.coords[\"is_reference\"]].squeeze()\n", + "sample = da[da.coords[\"is_data\"]] \n", + "\n", + "\n", + "print('Sample is', type(sample))\n", + "display(sample)\n", + "print('Dark is ', type(dark))\n", + "display(dark)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Plot fluo data for two single measurements" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# path to LOKI-like files\n", + "process_folder='/Users/gudlo523/Library/CloudStorage/OneDrive-LundUniversity/UU/ILL_947_June_2021/212/d22/exp_9-13-947/processed/ess_version'\n", + "os.chdir(process_folder)\n", + "\n", + "# Plot two fluo examples\n", + "name='066017.nxs' #exp5\n", + "#name='066050.nxs' #exp6\n", + "#name='066053.nxs' #exp6\n", + "#name='066083.nxs' #exp7\n", + "#name='066116.nxs' #exp8 \n", + "#name='065925.nxs' #exp2\n", + "#name='065962.nxs' #exp3\n", + "plot_fluo(name)\n", + "\n", + "exp2=[65925, 65927, 65930, 65933, 65936, 65939, 65942, 65945, 65948, 65951, 65954, 65957]\n", + "exp3= [65962, 65965, 65968, 65971, 65974, 65977, 65980, 65983, 65986, 65989, 65992]\n", + "\n", + "# converts numbers to full filenames\n", + "#flist_num=complete_fname(exp2)\n", + "\n", + "name='065925.nxs'\n", + "plot_fluo(name)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Plot UV data for exp7 and exp8" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# path to LOKI-like files\n", + "process_folder='/Users/gudlo523/Library/CloudStorage/OneDrive-LundUniversity/UU/ILL_947_June_2021/212/d22/exp_9-13-947/processed/ess_version'\n", + "os.chdir(process_folder)\n", + "\n", + "exp7= [66083, 66086, 66089, 66092, 66095, 66098, 66101, 66104, 66107, 66110, 66113]\n", + "exp8= [66116, 66119, 66122, 66125, 66128, 66131, 66134, 66137, 66140, 66143, 66146]\n", + "\n", + "\n", + "flist_num=complete_fname(exp7)\n", + "main_title='Exp 7'\n", + "pl1=plot_uv_set(flist_num, lambda_min=None, lambda_max=None, vmin=None, vmax=None)\n", + "pl1.ax.set_title(main_title)\n", + "\n", + "pl2=plot_uv_set(flist_num, lambda_min=200,lambda_max=400, vmin=-1, vmax=7)\n", + "pl2.ax.set_title(main_title)\n", + "#modify plots,\n", + "fig_handles=[pl1, pl2]\n", + "modify_plt_app(fig_handles)\n", + "\n", + "\n", + "uv_exp7_set=gather_uv_set(flist_num)\n", + "\n", + "flist_num=complete_fname(exp8)\n", + "main_title='Exp 8'\n", + "pl1=plot_uv_set(flist_num, lambda_min=None, lambda_max=None, vmin=None, vmax=None)\n", + "pl1.ax.set_title(main_title)\n", + "\n", + "pl2=plot_uv_set(flist_num, lambda_min=200,lambda_max=400, vmin=-1, vmax=10)\n", + "pl2.ax.set_title(main_title)\n", + "#modify plots\n", + "fig_handles=[pl1, pl2]\n", + "modify_plt_app(fig_handles)\n", + "\n", + "\n", + "uv_exp8_set=gather_uv_set(flist_num)\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Plot UV data for exp2 and exp3. In addition UV for two single measurements are shown." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "exp2=[65925, 65927, 65930, 65933, 65936, 65939, 65942, 65945, 65948, 65951, 65954, 65957]\n", + "exp3= [65962, 65965, 65968, 65971, 65974, 65977, 65980, 65983, 65986, 65989, 65992]\n", + "\n", + "\n", + "flist_num=complete_fname(exp2)\n", + "main_title='Exp 2'\n", + "pl1=plot_uv_set(flist_num, lambda_min=None, lambda_max=None, vmin=None, vmax=None)\n", + "pl1.ax.set_title(main_title)\n", + "\n", + "\n", + "pl2=plot_uv_set(flist_num, lambda_min=200,lambda_max=400, vmin=-2, vmax=4)\n", + "pl2.ax.set_title(main_title)\n", + "#modify plots\n", + "fig_handles=[pl1, pl2]\n", + "modify_plt_app(fig_handles)\n", + "\n", + "\n", + "uv_exp2_set=gather_uv_set(flist_num)\n", + "\n", + "name_s='065957.nxs'\n", + "plot_uv(name_s)\n", + "name_s='065925.nxs'\n", + "plot_uv(name_s)\n", + "\n", + "\n", + "flist_num=complete_fname(exp3)\n", + "main_title='Exp 3'\n", + "pl1=plot_uv_set(flist_num, lambda_min=None, lambda_max=None, vmin=None, vmax=None)\n", + "pl1.ax.set_title(main_title)\n", + "\n", + "\n", + "pl2=plot_uv_set(flist_num, lambda_min=200,lambda_max=400, vmin=-2, vmax=5)\n", + "pl2.ax.set_title(main_title)\n", + "#modify plots\n", + "fig_handles=[pl1, pl2]\n", + "modify_plt_app(fig_handles)\n", + "\n", + "\n", + "uv_exp3_set=gather_uv_set(flist_num)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Exploring filter possibilities\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Explore median filter\n", + "\n", + "# path to LOKI-like files\n", + "process_folder='/Users/gudlo523/Library/CloudStorage/OneDrive-LundUniversity/UU/ILL_947_June_2021/212/d22/exp_9-13-947/processed/ess_version'\n", + "os.chdir(process_folder)\n", + "\n", + "name='066017.nxs' #exp5\n", + "#name='066050.nxs' #exp6\n", + "#name='066053.nxs' #exp6\n", + "#name='066083.nxs' #exp7\n", + "#name='066116.nxs' #exp8 \n", + "#name='065925.nxs' #exp2\n", + "#name='065962.nxs' #exp3\n", + "\n", + "#plot_fluo(name)\n", + "\n", + "\n", + "fluo_dict=load_fluo(name)\n", + "fluo_da=normalize_fluo(**fluo_dict)\n", + "\n", + "display(fluo_da)\n", + "\n", + "spectra_number=7 \n", + "\n", + "def explore_medfilt(name, spectra_number=0, lmin=250, lmax=500):\n", + "\n", + " #load fluo data\n", + " fluo_dict=load_fluo(name)\n", + " fluo_da=normalize_fluo(**fluo_dict)\n", + "\n", + " # prepare a curve\n", + " yfluo=fluo_da['spectrum',spectra_number].values\n", + " xfluo=fluo_da.coords['wavelength'].values\n", + " \n", + " #slicing values\n", + " # lidx\n", + " lidx=np.where(np.logical_and((xfluo>=lmin), (xfluo<=lmax)))\n", + " xfluo_filt=xfluo[lidx]\n", + " yfluo_filt=yfluo[lidx]\n", + "\n", + " nsp=3\n", + " fig, ax = plt.subplots(1,nsp, figsize=(30,7))\n", + "\n", + "\n", + " ax[0].plot(xfluo,yfluo)\n", + " #yes, we want the zero line\n", + " ax[1].plot(xfluo,yfluo-yfluo)\n", + " ax[2].plot(xfluo_filt, yfluo_filt)\n", + " #example medfilt\n", + " kernel=[]\n", + " kernel_range=range(3,20,2)\n", + "\n", + "\n", + " for i in kernel_range:\n", + "\n", + " yfluomed=medfilt(yfluo, i)\n", + " ax[0].plot(xfluo, yfluomed)\n", + " ax[1].plot(xfluo,yfluo-yfluomed)\n", + " kernel.append(str(i))\n", + "\n", + " # zommed in\n", + " yfluomed_filt=medfilt(yfluo_filt, i)\n", + " ax[2].plot(xfluo_filt, yfluomed_filt)\n", + "\n", + " \n", + " colors=line_colors(len(kernel_range)+1)\n", + " #print('Colors contains,' ,len(colors))\n", + " ax[0].set_ylabel('Fluo intensity')\n", + " ax[1].set_ylabel('Difference between medfilter and signal')\n", + " ax[2].set_ylabel('Fluo intensity')\n", + "\n", + " title_string=f\"Influence of a median filter, {name}, spectrum {str(spectra_number)}\"\n", + "\n", + " for m in range(0,nsp):\n", + " ax[m].legend(['None', *kernel])\n", + " ax[m].set_xlabel('Wavelength [nm]')\n", + " ax[m].set_title(title_string)\n", + " ax[m].grid()\n", + " [ax[m].get_lines()[i].set_color(colors[i]) for i in range(0,len(kernel_range)+1)]\n", + " [ax[m].get_legend().legendHandles[i].set_color(colors[i]) for i in range(0,len(kernel_range)+1)]\n", + " [ax[m].get_lines()[i].set_marker(markers()[i]) for i in range(0,len(kernel_range)+1)]\n", + " [ax[m].get_lines()[i].set_markersize(2) for i in range(0,len(kernel_range)+1)]\n", + " #[ax[m].get_lines()[i].set_markevery(5) for i in range(0,len(kernel_range)+1)]\n", + " [ax[m].get_legend().legendHandles[i].set_marker(markers()[i]) for i in range(0,len(kernel_range)+1)]\n", + " [ax[m].get_lines()[i].set_linewidth(1) for i in range(0,len(kernel_range)+1)]\n", + "\n", + "\n", + "explore_medfilt(name, spectra_number=7, lmin=300, lmax=400)\n", + "\n", + "\n", + "# scipy.signal.medfilt not yet supported by scipp\n", + "# Simon suggests to check out\n", + "#https://scipp.github.io/generated/modules/scipp.signal.html\n", + "\n", + "\n", + "\n", + "# Check if fluo wavelength is regular. But it is not regular. \n", + "# Calculate the difference between two adjacent wavelength entries \n", + "xw=fluo_da.coords['wavelength'].values\n", + "xdiff = [xw[n]-xw[n-1] for n in range(1,len(xw))]\n", + "#print(xdiff)\n", + "\n", + "\n", + "def explore_butter(name, spectra_number=0, lmin=250, lmax=500):\n", + " from scipp.signal import butter, sosfiltfilt\n", + "\n", + " #load fluo data\n", + " fluo_dict=load_fluo(name)\n", + " fluo_da=normalize_fluo(**fluo_dict)\n", + "\n", + " # prepare a curve\n", + " yfluo=fluo_da['spectrum',spectra_number].values\n", + " x = fluo_da.coords['wavelength']\n", + " # conversion to float64 due to a bug\n", + " fluo_da.coords['wavelength'] = sc.linspace(x.dim, x.values[0], x.values[-1], num=len(x), unit=x.unit, dtype='float64')\n", + " #out=butter(fluo_da.coords['wavelength'], N=12, Wn=0.04/ x.unit).filtfilt(fluo_da.data, 'wavelength') #this applies the butter filter to the whole fluo dataarray\n", + " #out=butter(fluo_da['spectrum',spectra_number].coords['wavelength'], N=12, Wn=0.04/ x.unit).filtfilt(fluo_da['spectrum',spectra_number].data, 'wavelength') # I want to apply the butter filter only to one spectrum\n", + "\n", + " #slicing values\n", + " # lidx\n", + " lidx=np.where(np.logical_and((fluo_da.coords['wavelength'].values>=lmin), (fluo_da.coords['wavelength'].values<=lmax)))\n", + " xfluo_filt=x.values[lidx]\n", + " yfluo_filt=yfluo[lidx]\n", + "\n", + "\n", + " \n", + " # butter parameters\n", + " Wn_range=[0.5,0.3,0.04]\n", + " N_range=[12,4]\n", + " \n", + " nsp=3\n", + " fig, ax = plt.subplots(1,nsp, figsize=(30,7))\n", + " ax[0].plot(x.values,fluo_da['spectrum',spectra_number].values)\n", + " ax[1].plot(x.values, fluo_da['spectrum',spectra_number].values-fluo_da['spectrum',spectra_number].values)\n", + " ax[2].plot(xfluo_filt,yfluo_filt)\n", + "\n", + " ax[0].set_ylabel('Fluo intensity')\n", + " ax[2].set_ylabel('Fluo intensity')\n", + " ax[1].set_ylabel('Difference between butter and signal')\n", + "\n", + " para_list=[]\n", + " for i in Wn_range:\n", + " for m in N_range:\n", + " out=butter(fluo_da['spectrum',spectra_number].coords['wavelength'], N=m, Wn=i/ x.unit).filtfilt(fluo_da['spectrum',spectra_number].data, 'wavelength') \n", + " ax[0].plot(x.values, out.values)\n", + " para_out=[m,i]\n", + " para_list.append(para_out)\n", + " ax[1].plot(x.values,fluo_da['spectrum',spectra_number].values-out.values)\n", + " ax[2].plot(xfluo_filt,out.values[lidx])\n", + "\n", + " kernel_length=len(Wn_range)*len(N_range)\n", + " kernel_range=para_list\n", + "\n", + " #ax[0].plot(x.values, out['spectrum',spectra_number].values)\n", + " #ax[0].plot(x.values, out.values)\n", + "\n", + " title_string=f\"Influence of a butter filter, {name}, spectrum {str(spectra_number)}\"\n", + "\n", + "\n", + " colors=line_colors(len(kernel_range)+1)\n", + " for m in range(0,nsp):\n", + " ax[m].legend(['None', *kernel_range])\n", + " ax[m].set_xlabel('Wavelength [nm]')\n", + " ax[m].set_title(title_string)\n", + " ax[m].grid()\n", + " [ax[m].get_lines()[i].set_color(colors[i]) for i in range(0,len(kernel_range)+1)]\n", + " [ax[m].get_legend().legendHandles[i].set_color(colors[i]) for i in range(0,len(kernel_range)+1)]\n", + " [ax[m].get_lines()[i].set_marker(markers()[i]) for i in range(0,len(kernel_range)+1)]\n", + " [ax[m].get_lines()[i].set_markersize(2) for i in range(0,len(kernel_range)+1)]\n", + " #[ax[m].get_lines()[i].set_markevery(5) for i in range(0,len(kernel_range)+1)]\n", + " [ax[m].get_legend().legendHandles[i].set_marker(markers()[i]) for i in range(0,len(kernel_range)+1)]\n", + " [ax[m].get_lines()[i].set_linewidth(1) for i in range(0,len(kernel_range)+1)]\n", + "\n", + "\n", + "\n", + "explore_butter(name, spectra_number=7, lmin=300, lmax=400)\n", + "\n", + "\n", + "# there is a bug, hence conversion to float64\n", + "#fluo_da.coords['wavelength'] = fluo_da.coords.pop('wavelength').to(dtype='float64')\n", + "#out = butter(fluo_da.coords['wavelength'], N=4, Wn=20 / fluo_da.coords['wavelength'].unit).filtfilt(fluo_da,'wavelength')\n", + "# Above line fails currently, because the coords for wavelength are not regular\n", + "\n", + "#Simon suggestion\n", + "\n", + "#x = fluo_da.coords['wavelength']\n", + "#fluo_da.coords['wavelength'] = sc.linspace(x.dim, x.values[0], x.values[-1], num=len(x), unit=x.unit, dtype='float64')\n", + "#out=butter(fluo_da.coords['wavelength'], N=12, Wn=0.04/ x.unit).filtfilt(fluo_da.data, 'wavelength')\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Developing functions for a data analysis pipeline for UV, fluo spectroscopy\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# path to LOKI-like files\n", + "process_folder='/Users/gudlo523/Library/CloudStorage/OneDrive-LundUniversity/UU/ILL_947_June_2021/212/d22/exp_9-13-947/processed/ess_version'\n", + "os.chdir(process_folder)\n", + "\n", + "name='065992.nxs'\n", + "#load uv data\n", + "uv_dict=load_uv(name)\n", + "#display(uv_dict['data'],uv_dict['reference'],uv_dict['dark'])\n", + "uv_da=normalize_uv(**uv_dict) #returns sc.DataArray with all uv spectra present in file\n", + "#display(uv_da)\n", + "\n", + "plot_uv(name)\n", + "uv_turbidity_fit(uv_da, wl_unit=sc.Unit('nm'), fit_llim=500, fit_ulim=850, b_llim=300, b_ulim=400,m=0.1, plot_corrections=True)\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### UV pipeline" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# path to LOKI-like files\n", + "process_folder='/Users/gudlo523/Library/CloudStorage/OneDrive-LundUniversity/UU/ILL_947_June_2021/212/d22/exp_9-13-947/processed/ess_version'\n", + "os.chdir(process_folder)\n", + "\n", + "\"\"\" What functions are available for uv data? Overview on only the function names\n", + "#How to apply a medilter to fluo or uv data? \n", + "apply_medfilter\n", + "\n", + "#This is for uv\n", + "#Extract uv peak intensity\n", + "uv_peak_int\n", + "\n", + "#Plot uv peak intensity\n", + "plot_uv_peak_int\n", + "\n", + "#Quickly check uv data\n", + "uv_quick_data_check\n", + "\n", + "#plot mulitple uv peak int\n", + "plot_multiple_uv_peak_int\n", + "\n", + "#turbidity\n", + "turbidity\n", + "\n", + "#residuals\n", + "residual\n", + "\n", + "#UV turbidity fit\n", + "uv_turbidity_fit\n", + "\n", + "#multi turbidity fit\n", + "multi_uv_turbidity_fit\n", + "\"\"\"\n", + "\n", + "#%%%%%%%%%%%%%%%%%%%%%%%%%%% Code examples %%%%%%%%%%%%%%%%%%%%%%%%%%\n", + "\n", + "name='066017.nxs'\n", + "#load uv data\n", + "uv_dict=load_uv(name)\n", + "#display(uv_dict['data'],uv_dict['reference'],uv_dict['dark'])\n", + "uv_da=normalize_uv(**uv_dict) #returns sc.DataArray with all uv spectra present in a LoKI.nxs file\n", + "#display(uv_da)\n", + "\n", + "#How to apply a medfilter to a uv dataarray\n", + "#apply_medfilter(uv_da,kernel_size=9) #applies median filter to multiple spectra in uv_da\n", + "#apply_medfilter(process_uv(name),kernel_size=9) #applies median filter to normalised spectra resulting from process_uv\n", + "\n", + "# Plots uv peak intensity for a given nexus file\n", + "uv_280_result=uv_peak_int(uv_da, wavelength=None, wl_unit=sc.Unit('nm'), tol=None)\n", + "plot_uv_peak_int(uv_da, name, wavelength=None, wl_unit=sc.Unit('nm'), tol=1)\n", + "\n", + "# exp 5 and exp 6\n", + "exp5=[66017, 66020, 66023, 66026, 66029, 66032, 66034, 66037, 66040, 66043, 66046]\n", + "#exp6= [66050, 66053, 66056, 66059, 66062, 66065, 66068, 66071, 66074, 66077, 66080]\n", + "\n", + "# Create complete filenames \n", + "#flist_num=complete_fname(exp5)\n", + "#flist_num=complete_fname(exp6)\n", + "\n", + "# a quic visualisation check if all files have equivalent number of spectra\n", + "#uv_quick_data_check(flist_num, wavelength=None, wl_unit=None, tol=None, medfilter=False, kernel_size=None)\n", + "\n", + "# select peak intensities at 280nm and plot them\n", + "plot_multiple_uv_peak_int(flist_num, wavelength=280, wl_unit=None, tol=None, medfilter=True, kernel_size=None)\n", + "\n", + "# Plot the spectra in a single LoKI.nxs file\n", + "#plot_uv('066020.nxs')\n", + "\n", + "# Plot multiple UV spectra for a given set of files\n", + "#plot_uv_set(flist_num)\n", + "\n", + "# Example for tubidity\n", + "#uv_turbidity_fit(uv_da, wl_unit=sc.Unit('nm'), fit_llim=300, fit_ulim=850, b_llim=450, b_ulim=700,m=0.1, plot_corrections=True)\n", + "\n", + "# Performs for a given set of files turbidity correction\n", + "uv_da_multi=uv_multi_turbidity_fit(flist_num, wl_unit=sc.Unit('nm'), fit_llim=300, fit_ulim=850, b_llim=450, b_ulim=700,m=0.1, plot_corrections=True)\n", + "\n", + "\n", + "\n", + "\n", + " " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Fluo pipeline" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# path to LOKI-like files\n", + "process_folder='/Users/gudlo523/Library/CloudStorage/OneDrive-LundUniversity/UU/ILL_947_June_2021/212/d22/exp_9-13-947/processed/ess_version'\n", + "os.chdir(process_folder)\n", + "\n", + "\"\"\" What functions are available for fluo data? Overview on only the function names\n", + "#How to apply a medilter to fluo or uv data? \n", + "apply_medfilter\n", + "\n", + "#Fluo peak intensity\n", + "fluo_peak_int\n", + "\n", + "#Plot fluo peak intensity\n", + "plot_fluo_peak_int\n", + "\n", + "#Plot multiple peak intensities\n", + "plot_fluo_multiple_peak_int\n", + "\"\"\"\n", + "\n", + "#Example of how to load one nexusfile and extract the fluo dataarray\n", + "#name='066017.nxs' #exp5\n", + "#fluo_dict=load_fluo(name)\n", + "#fluo_da=normalize_fluo(**fluo_dict)\n", + "#display(fluo_da) \n", + "\n", + "\n", + "# Example of how to apply the med filter to a fluo dataarray fluo_da\n", + "#apply_medfilter(fluo_da,kernel_size=9)\n", + "#apply_medfilter(fluo_da['spectrum', 7:9], kernel_size=9)\n", + "#apply_medfilter(fluo_da['spectrum', 7], kernel_size=9) \n", + "#plot_fluo(name)\n", + "#fluo_filt_max=fluo_peak_int(fluo_da, wllim=None, wulim=None, wl_unit=None, medfilter=True, kernel_size=15)\n", + "#display(fluo_filt_max)\n", + "#plot_fluo_peak_int(fluo_da,name, wllim=300, wulim=400, wl_unit=None, medfilter=True, kernel_size=15)\n", + "\n", + "exp5=[66017, 66020, 66023, 66026, 66029, 66032, 66034, 66037, 66040, 66043, 66046]\n", + "flist_num5=complete_fname(exp5)\n", + "exp6= [66050, 66053, 66056, 66059, 66062, 66065, 66068, 66071, 66074, 66077, 66080]\n", + "flist_num6=complete_fname(exp6)\n", + "\n", + "\n", + "def plot_fluo_multiple_peak_int(filelist, wllim=None, wulim=None, wl_unit=None, medfilter=True, kernel_size=None):\n", + " \"\"\" Plots multiple peak intensities for fluo spectr\n", + " \n", + " \"\"\"\n", + " import itertools\n", + " marker = itertools.cycle(markers()) \n", + "\n", + "\n", + " print(filelist)\n", + "\n", + " figure_size=(15,5)\n", + " fig, ax=plt.subplots(nrows=1,ncols=2,figsize=figure_size, constrained_layout=True)\n", + "\n", + " \n", + " unique_mwl=[]\n", + " ds_list=[]\n", + " for name in filelist:\n", + " fluo_dict=load_fluo(name)\n", + " fluo_da=normalize_fluo(**fluo_dict)\n", + " #extract max int value and corresponding wavelength position, median filter is applied\n", + " fluo_filt_max=fluo_peak_int(fluo_da, wllim=wllim, wulim=wulim, wl_unit=wl_unit, medfilter=medfilter, kernel_size=kernel_size) \n", + " # attach filename as attribute to dataarray\n", + " #fluo_filt_max.attrs['filename'] = sc.scalar(name)\n", + " #display(fluo_filt_max)\n", + " ds_list.append(fluo_filt_max)\n", + " unique_mwl.append(np.unique(fluo_filt_max.coords['monowavelengths'].values))\n", + " #print(fluo_filt_max)\n", + "\n", + " #same marker for both plots for the same file\n", + " markerchoice=next(marker)\n", + "\n", + " ax[0].plot(fluo_filt_max.coords['monowavelengths'].values, fluo_filt_max['intensity_max'].values, label=f'{name}', linestyle=\"None\", marker=markerchoice, markersize=10)\n", + " ax[0].set_ylabel('Max. Intensity')\n", + " ax[0].set_title('Fluo - max. intensity')\n", + "\n", + " ax[1].plot(fluo_filt_max.coords['monowavelengths'].values, fluo_filt_max['wavelength_max'].values, label=f'{name}', linestyle=\"None\", marker=markerchoice, markersize=10)\n", + " unit_str=str(fluo_filt_max['wavelength_max'].unit)\n", + "\n", + " ax[1].set_ylabel(f'Wavelength [{unit_str}]')\n", + " ax[1].set_title(f'Fluo - corresponding wavelength')\n", + "\n", + " # show the lowest monowavelength as lower boundary on the y-axis\n", + " ax[1].set_ylim(bottom=0.9*np.min(fluo_filt_max.coords['monowavelengths'].values))\n", + "\n", + " # plot the found monowavelengths as additional visual information on the y-axis \n", + " for mwl in np.unique(unique_mwl):\n", + " ax[1].plot(np.unique(unique_mwl), np.full(np.shape(np.unique(unique_mwl)), mwl) , '--', label=f\"{mwl}{sc.Unit('nm')}\")\n", + " #ax[1].legend(loc='upper right', bbox_to_anchor=(1.05, 1.05)) \n", + " \n", + " for axes in ax:\n", + " #axes.legend(loc='upper right', bbox_to_anchor=(1.1, 1.00)) \n", + " axes.legend( bbox_to_anchor=(1.04, 1)) \n", + " axes.grid(True)\n", + " axes.set_xlabel('Monowavelengths')\n", + " \n", + " \n", + " display(fig)\n", + "\n", + "plot_fluo_multiple_peak_int(flist_num5,wllim=300, wulim=400, wl_unit=None, medfilter=True, kernel_size=15)\n", + "plot_fluo_multiple_peak_int(flist_num6,wllim=300, wulim=400, wl_unit=None, medfilter=True, kernel_size=15)\n", + "\n", + "\n", + "#If I understood scipp correctly, I can only have a Dataset, where each DataArray is the same. So if one DataArray has 2 spectra, one DataArray only one, it should not be possible to put them together in a Dataset.\n", + "#Correct. But you could have an 1-D data array and a 2-D data array in the same dataset, if they have the same length along the shared dimension.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#How to access for a given filename the fluo dataarray.\n", + "name='066017.nxs'\n", + "fluo_dict=load_fluo(name)\n", + "fluo_da=normalize_fluo(**fluo_dict)\n", + "display(fluo_da)\n", + "print(f'Number of fluo spectra in {name}: {fluo_da.sizes[\"spectrum\"]}' )\n", + "\n", + "#apply medfilter to fluo sc.DataArray\n", + "fluo_da_filt=apply_medfilter(fluo_da, kernel_size=15)\n", + "\n", + "#showcasing scipp dataarray graphical presentation\n", + "sc.show(fluo_da_filt)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "#### Extract from fluorescence measurements the maximum intensity and the corresponding wavelength.\n", + "\n", + "#Fluo Max int, max wavelength\n", + "\n", + "exp5=[66017, 66020, 66023, 66026, 66029, 66032, 66034, 66037, 66040, 66043, 66046]\n", + "flist_num=complete_fname(exp5)\n", + "\n", + "fluo_int_dict=fluo_maxint_max_wavelen(flist_num,wllim=300, wulim=400, wl_unit=None, medfilter=True, kernel_size=15)\n", + "\n", + "fluo_plot_maxint_max_wavelen(fluo_int_dict) " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Plot a selected fluo spectrum across a series of measurements" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Here we plot a certain fluo spectrum for all measurements in a whole experimental series. Here spectrum #1 corresponds to the white light spectrum.\n", + "\n", + "exp5=[66017, 66020, 66023, 66026, 66029, 66032, 66034, 66037, 66040, 66043, 66046]\n", + "flist_num5=complete_fname(exp5)\n", + "exp6= [66050, 66053, 66056, 66059, 66062, 66065, 66068, 66071, 66074, 66077, 66080]\n", + "flist_num6=complete_fname(exp6)\n", + "\n", + "fwls=plot_fluo_spectrum_selection(flist_num5, spectral_idx=1, kernel_size=15, wllim=400, wulim=500)\n", + "fwls=plot_fluo_spectrum_selection(flist_num6, spectral_idx=1, kernel_size=15, wllim=400, wulim=500)\n", + "\n" + ] + } + ], + "metadata": { + "interpreter": { + "hash": "410bfb39d85b3112e66f48ab3e53bb74ad7e7b3fb364f756ca860b5a3cf79ca2" + }, + "kernelspec": { + "display_name": "Python 3.8.12 ('scippneutron')", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.8.12" + }, + "orig_nbformat": 4 + }, + "nbformat": 4, + "nbformat_minor": 2 +}