diff --git a/.cspell/custom-dictionary.txt b/.cspell/custom-dictionary.txt index aaab614b..f634cba4 100644 --- a/.cspell/custom-dictionary.txt +++ b/.cspell/custom-dictionary.txt @@ -172,6 +172,7 @@ levelname linalg linekwds linesegkwds +linestyles linewidth linspace literalinclude @@ -325,6 +326,7 @@ scandir scatt scatterkwds scicat +scipy SDIAG sdir segs @@ -341,6 +343,7 @@ stepsize subdir subdirs subfolders +suptitle symscores targcenter termorder @@ -367,6 +370,7 @@ TZCYXS tzoffset ubid UDLD +ufunc unbinned uncategorised undoc @@ -375,6 +379,7 @@ varnames venv verts viewcode +vlines vmax voxels VTOF diff --git a/docs/scripts/build_flash_parquets.py b/docs/scripts/build_flash_parquets.py index d3e45c26..8a679969 100644 --- a/docs/scripts/build_flash_parquets.py +++ b/docs/scripts/build_flash_parquets.py @@ -27,3 +27,25 @@ system_config=config_file, collect_metadata=False, ) + +dataset.get("W110", root_dir="./tutorial") +data_path = dataset.dir + + +config_override = { + "core": { + "paths": { + "data_raw_dir": data_path, + "data_parquet_dir": data_path + "/processed/", + }, + }, +} + +runs = ["44498", "44455"] +for run in runs: + sp = SedProcessor( + runs=run, + config=config_override, + system_config=config_file, + collect_metadata=False, + ) diff --git a/docs/scripts/download_data.py b/docs/scripts/download_data.py index cf1e042e..1bbd2e5b 100644 --- a/docs/scripts/download_data.py +++ b/docs/scripts/download_data.py @@ -6,3 +6,4 @@ dataset.get("Gd_W110", remove_zip=True, root_dir=root_dir) dataset.get("TaS2", remove_zip=True, root_dir=root_dir) dataset.get("Au_Mica", remove_zip=True, root_dir=root_dir) +dataset.get("W110", remove_zip=True, root_dir=root_dir) diff --git a/docs/user_guide/advanced_topics.md b/docs/user_guide/advanced_topics.md index ebd773c8..d35c19ce 100644 --- a/docs/user_guide/advanced_topics.md +++ b/docs/user_guide/advanced_topics.md @@ -3,4 +3,6 @@ ../tutorial/6_binning_with_time-stamped_data ../tutorial/7_correcting_orthorhombic_symmetry ../tutorial/8_jittering_tutorial +../tutorial/10_hextof_workflow_trXPS_bam_correction +../tutorial/11_hextof_workflow_trXPS_energy_calibration_using_SB ``` diff --git a/docs/workflows/index.md b/docs/workflows/index.md index 5acb9c04..a1219085 100644 --- a/docs/workflows/index.md +++ b/docs/workflows/index.md @@ -8,5 +8,6 @@ myst: ```{toctree} ../tutorial/4_hextof_workflow -../tutorial/5_sxp_workflow.ipynb +../tutorial/5_sxp_workflow +../tutorial/9_hextof_workflow_trXPD ``` diff --git a/sed/core/processor.py b/sed/core/processor.py index 882c5370..016ad19e 100644 --- a/sed/core/processor.py +++ b/sed/core/processor.py @@ -1343,14 +1343,15 @@ def calibrate_energy_axis( **kwds, ) if verbose: - print("Quality of Calibration:") - self.ec.view( - traces=self.ec.traces_normed, - xaxis=self.ec.calibration["axis"], - align=True, - energy_scale=energy_scale, - backend="bokeh", - ) + if self.ec.traces_normed is not None: + print("Quality of Calibration:") + self.ec.view( + traces=self.ec.traces_normed, + xaxis=self.ec.calibration["axis"], + align=True, + energy_scale=energy_scale, + backend="bokeh", + ) print("E/TOF relationship:") self.ec.view( traces=self.ec.calibration["axis"][None, :], diff --git a/sed/dataset/datasets.json b/sed/dataset/datasets.json index 914213a5..749a3098 100644 --- a/sed/dataset/datasets.json +++ b/sed/dataset/datasets.json @@ -14,6 +14,14 @@ ], "rearrange_files": true }, + "W110": { + "url": "https://zenodo.org/records/12609441/files/single_event_data.zip", + "subdirs": [ + "analysis_data", + "calibration_data" + ], + "rearrange_files": true + }, "TaS2": { "url": "https://zenodo.org/records/10160182/files/TaS2.zip", "subdirs": [ diff --git a/tutorial/10_hextof_workflow_trXPS_bam_correction.ipynb b/tutorial/10_hextof_workflow_trXPS_bam_correction.ipynb new file mode 100644 index 00000000..1b7e0d54 --- /dev/null +++ b/tutorial/10_hextof_workflow_trXPS_bam_correction.ipynb @@ -0,0 +1,561 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "9192f0aa", + "metadata": {}, + "source": [ + "# Tutorial for trXPS for the HEXTOF instrument at FLASH: t0, cross-correlation and BAM correction" + ] + }, + { + "cell_type": "markdown", + "id": "1e583d01", + "metadata": {}, + "source": [ + "## Preparation" + ] + }, + { + "cell_type": "markdown", + "id": "303d7dfb", + "metadata": {}, + "source": [ + "### Import necessary libraries" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e73f4a19", + "metadata": {}, + "outputs": [], + "source": [ + "%load_ext autoreload\n", + "%autoreload 2\n", + "\n", + "from pathlib import Path\n", + "import os\n", + "\n", + "from sed import SedProcessor\n", + "from sed.dataset import dataset\n", + "import numpy as np\n", + "\n", + "%matplotlib widget\n", + "import matplotlib.pyplot as plt\n", + "\n", + "# For peak fitting\n", + "from lmfit.models import GaussianModel" + ] + }, + { + "cell_type": "markdown", + "id": "07076e2f", + "metadata": {}, + "source": [ + "### Get data paths\n", + "\n", + "If it is your beamtime, you can read the raw data and write to the processed directory. For the public data, you can not write to the processed directory.\n", + "\n", + "The paths are such that if you are on Maxwell, it uses those. Otherwise, data is downloaded in the current directory from Zenodo:\n", + "https://zenodo.org/records/12609441" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f9713ebb", + "metadata": {}, + "outputs": [], + "source": [ + "beamtime_dir = \"/asap3/flash/gpfs/pg2/2023/data/11019101\" # on Maxwell\n", + "if os.path.exists(beamtime_dir) and os.access(beamtime_dir, os.R_OK):\n", + " path = beamtime_dir + \"/raw/hdf/offline/fl1user3\"\n", + " buffer_path = beamtime_dir + \"/processed/tutorial/\"\n", + "else:\n", + " # data_path can be defined and used to store the data in a specific location\n", + " dataset.get(\"W110\") # Put in Path to a storage of at least 10 Byte free space.\n", + " path = dataset.dir\n", + " buffer_path = path + \"/processed/\"" + ] + }, + { + "cell_type": "markdown", + "id": "5cf71d69", + "metadata": {}, + "source": [ + "### Config setup\n", + "Here, we get the path to the config file and set up the relevant directories. This can also be done directly in the config file." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "195d7ff5", + "metadata": {}, + "outputs": [], + "source": [ + "# pick the default configuration file for hextof@FLASH\n", + "config_file = Path('../sed/config/flash_example_config.yaml')\n", + "assert config_file.exists()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "18c3db1b", + "metadata": {}, + "outputs": [], + "source": [ + "# here we setup a dictionary that will be used to override the path configuration\n", + "config_override = {\n", + " \"core\": {\n", + " \"beamtime_id\": 11019101,\n", + " \"paths\": {\n", + " \"data_raw_dir\": path,\n", + " \"data_parquet_dir\": buffer_path\n", + " },\n", + " },\n", + "}" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3145b03c", + "metadata": {}, + "outputs": [], + "source": [ + "energy_cal = {\n", + " \"energy\": {\n", + " \"calibration\": {\n", + " \"E0\": -53.96145014592986,\n", + " \"creation_date\": 1732056868.029444,\n", + " \"d\": 0.8096677233434938,\n", + " \"energy_scale\": \"kinetic\",\n", + " \"t0\": 4.0148196718030886e-07,\n", + " },\n", + " \"offsets\":{\n", + " \"constant\": -77.5,\n", + " \"creation_date\": 1732056874.060922,\n", + " \"monochromatorPhotonEnergy\": {\n", + " \"preserve_mean\": True,\n", + " \"weight\": -1,\n", + " },\n", + " \"sampleBias\": {\n", + " \"preserve_mean\": False,\n", + " \"weight\": 1,\n", + " },\n", + " \"tofVoltage\": {\n", + " \"preserve_mean\": True,\n", + " \"weight\": -1,\n", + " },\n", + " },\n", + " },\n", + "}" + ] + }, + { + "cell_type": "markdown", + "id": "29a56a9e", + "metadata": {}, + "source": [ + "### We use the stored energy calibration parameters and load trXPS data set to define:\n", + "* t0 position with respect to delay stage values;\n", + "* correct accordingly delay stage offset\n", + "* fit cross-correlation \n", + "* apply BAM correction and see its effect on cross-correlation" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "678d23c4", + "metadata": {}, + "outputs": [], + "source": [ + "run_number = 44498\n", + "sp_44498 = SedProcessor(runs=[run_number], config=config_override, folder_config=energy_cal, system_config=config_file, verbose=True)\n", + "\n", + "sp_44498.add_jitter()\n", + "sp_44498.align_dld_sectors()\n", + "sp_44498.append_energy_axis()\n", + "sp_44498.add_energy_offset()" + ] + }, + { + "cell_type": "markdown", + "id": "ba9bb7f7", + "metadata": {}, + "source": [ + "Check which channels are included in the dataframe" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c17be152", + "metadata": {}, + "outputs": [], + "source": [ + "sp_44498.dataframe.head()" + ] + }, + { + "cell_type": "markdown", + "id": "db246a19", + "metadata": {}, + "source": [ + "## Data w/o BAM correction" + ] + }, + { + "cell_type": "markdown", + "id": "1349568e", + "metadata": {}, + "source": [ + "First, we take a look at our sideband measurement before any corrections.\n", + "The sidebands on the W4f core levels can be used as a measure of the pump and probe cross-correlation,\n", + "and hence our temporal resolution.\n", + "We plot the data delay stage position vs Energy data, normalized by acquisition time." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3df8530e", + "metadata": {}, + "outputs": [], + "source": [ + "axes = ['energy', 'delayStage']\n", + "ranges = [[-37.5,-27.5], [1446.75,1449.15]]\n", + "bins = [200,40]\n", + "res = sp_44498.compute(bins=bins, axes=axes, ranges=ranges, normalize_to_acquisition_time=\"delayStage\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "59764b73", + "metadata": {}, + "outputs": [], + "source": [ + "fig,ax = plt.subplots(1,2,figsize=(8,3), layout='constrained')\n", + "res.plot(robust=True, ax=ax[0], cmap='terrain')\n", + "fig.suptitle(f\"Run {run_number}: W 4f, side bands\")\n", + "ax[0].set_title('raw')\n", + "bg = res.sel(delayStage=slice(1448.7,1449.1)).mean('delayStage')\n", + "(res.sel(delayStage=slice(1446.8,1449.3))-bg).plot(robust=True, ax=ax[1])\n", + "ax[1].set_title('difference')" + ] + }, + { + "cell_type": "markdown", + "id": "43191608", + "metadata": {}, + "source": [ + "Now we make fit to determine precise t$_0$ position and cross-correlation using lmfit fit models" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4de472b3", + "metadata": {}, + "outputs": [], + "source": [ + "Gauss_mod = GaussianModel()\n", + "\n", + "#first order sideband:\n", + "x1=res['delayStage']\n", + "y1=res.sel(energy=slice(-30.5,-29.5)).sum('energy')\n", + "y1=y1-np.mean(y1.sel(delayStage=slice(1448.7,1449.1)))\n", + "\n", + "pars1 = Gauss_mod.make_params(amplitude=0.1, center=1447.8, sigma=0.02)\n", + "out1 = Gauss_mod.fit(y1, pars1, x=x1)\n", + "\n", + "#second order sideband\n", + "x2=res['delayStage']\n", + "y2=res.sel(energy=slice(-29.5,-28.5)).sum('energy')\n", + "y2=y2-np.mean(y2.sel(delayStage=slice(1448.7,1449.1)))\n", + "\n", + "pars2 = Gauss_mod.make_params(amplitude=0.1, center=1447.8, sigma=0.02)\n", + "out2 = Gauss_mod.fit(y2, pars2, x=x2)\n", + "\n", + "plt.figure()\n", + "plt.plot(x1,y1,'rx', label='$1^{st}$ order sideband')\n", + "plt.plot(x1,out1.best_fit,'r', label=\"FWHM = {:.3f} ps\".format(out1.values['fwhm']))\n", + "plt.legend(loc=\"best\")\n", + "plt.title('run44498, W4f, sidebands comparison')\n", + "plt.plot(x2,y2,'bx', label='$2^{nd}$ order sideband')\n", + "plt.plot(x2,out2.best_fit,'b', label=\"FWHM = {:.3f} ps\".format(out2.values['fwhm']))\n", + "plt.legend(loc=\"best\")\n", + "plt.xlabel(\"delayStage [ps]\")\n", + "plt.ylabel(\"Intensity [cts/s]\")\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "f6a20bf2", + "metadata": {}, + "source": [ + "As we see the sidebands are quite broad and one of the possible reasons for this could be long or short-term drifts (jitter) of the FEL arrival time with respect to e.g. optical laser or differences in the intra-bunch arrival time. To check and correct for this we can look at beam arrival monitor (BAM). The BAM gives a pulse-resolved measure of the FEL arrival time with respect to a master clock." + ] + }, + { + "cell_type": "markdown", + "id": "2742bb17", + "metadata": {}, + "source": [ + "## Check BAM versus pulse and train IDs" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8feb38c7", + "metadata": {}, + "outputs": [], + "source": [ + "axes = ['trainId', 'pulseId', 'bam']\n", + "ranges = [[1628022640,1628046700], [0,500], [-6400,100]]\n", + "bins = [250, 100, 1000]\n", + "res_bam = sp_44498.compute(bins=bins, axes=axes, ranges=ranges)" + ] + }, + { + "cell_type": "markdown", + "id": "de89978e", + "metadata": {}, + "source": [ + "As we can see, jitter between FEL and pump laser is quite significant withing a pulse train as well as over the whole measurement period." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8203d22d", + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "fig,ax = plt.subplots(1,2,figsize=(8,3), layout='constrained')\n", + "res_bam.sel(bam=slice(-6400,-5100)).sum('trainId').plot(ax=ax[0],robust=True, cmap='terrain')\n", + "res_bam.sel(bam=slice(-6400,-5100)).sum('pulseId').plot(ax=ax[1],robust=True, cmap='terrain')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "babbc574", + "metadata": {}, + "source": [ + "## Apply BAM correction" + ] + }, + { + "cell_type": "markdown", + "id": "b0c70ad0", + "metadata": {}, + "source": [ + "To correct the SASE jitter, using information from the bam column and to calibrate the pump-probe delay axis, we need to shift the delay stage values to centre the pump-probe-time overlap time zero." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "90b702e4", + "metadata": {}, + "outputs": [], + "source": [ + "sp_44498.add_delay_offset(\n", + " constant=-1448, # this is time zero position determined from side band fit\n", + " flip_delay_axis=True, # invert the direction of the delay axis\n", + " columns=['bam'], # use the bam to offset the values\n", + " weights=[-0.001], # bam is in fs, delay in ps\n", + " preserve_mean=True # preserve the mean of the delay axis to keep t0 position\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "9f83b41e", + "metadata": {}, + "source": [ + "### bin in the corrected delay axis" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "89a72e3a", + "metadata": {}, + "outputs": [], + "source": [ + "axes = ['energy', 'delayStage']\n", + "ranges = [[-37.5,-27.5], [-1.5,1.5]]\n", + "bins = [200,60]\n", + "res_corr = sp_44498.compute(bins=bins, axes=axes, ranges=ranges, normalize_to_acquisition_time=\"delayStage\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8443ddcb", + "metadata": {}, + "outputs": [], + "source": [ + "fig,ax = plt.subplots(1,2,figsize=(8,3), layout='constrained')\n", + "fig.suptitle(f\"Run {run_number}: W 4f, side bands\")\n", + "res_corr.plot(robust=True, ax=ax[0], cmap='terrain')\n", + "ax[0].set_title('raw')\n", + "bg = res_corr.sel(delayStage=slice(-1.3,-1.0)).mean('delayStage')\n", + "(res_corr-bg).plot(robust=True, ax=ax[1])\n", + "ax[1].set_title('difference')" + ] + }, + { + "cell_type": "markdown", + "id": "202404c8", + "metadata": {}, + "source": [ + "We clearly see an effect of BAM corrections - side bands are visible much nicer and width became smaller." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3141cc76", + "metadata": {}, + "outputs": [], + "source": [ + "sp_44498.save_delay_offsets()" + ] + }, + { + "cell_type": "markdown", + "id": "67deb9f9", + "metadata": {}, + "source": [ + "Now we can repeat fit procedure to determine true cross-correlation value." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3257d2bb", + "metadata": {}, + "outputs": [], + "source": [ + "Gauss_mod = GaussianModel()\n", + "\n", + "#first order sideband:\n", + "x5=res_corr['delayStage'].sel(delayStage=slice(-1.6,1.5))\n", + "y5=res_corr.sel(energy=slice(-30.4,-29.5),delayStage=slice(-1.6,1.5)).sum('energy')\n", + "y5=y5-np.mean(y5.sel(delayStage=slice(-1.4,-1.0)))\n", + "\n", + "pars5 = Gauss_mod.make_params(amplitude=0.1, center=0.0, sigma=0.02)\n", + "out5 = Gauss_mod.fit(y5, pars5, x=x5)\n", + "\n", + "print(out5.fit_report())\n", + "\n", + "#second order sideband\n", + "x6=res_corr['delayStage'].sel(delayStage=slice(-1.6,1.5))\n", + "y6=res_corr.sel(energy=slice(-29.5,-27.5),delayStage=slice(-1.6,1.5)).sum('energy')\n", + "y6=y6-np.mean(y6.sel(delayStage=slice(-1.4,-1.0)))\n", + "\n", + "pars6 = Gauss_mod.make_params(amplitude=0.1, center=0.0, sigma=0.02)\n", + "out6 = Gauss_mod.fit(y6, pars6, x=x6)\n", + "\n", + "print(out6.fit_report())\n", + "\n", + "#comparison plot\n", + "plt.figure()\n", + "plt.plot(x5,y5,'rx', label='$1^{st}$ order sideband')\n", + "plt.plot(x5,out5.best_fit,'r', label=\"FWHM = {:.3f} ps\".format(out5.values['fwhm']))\n", + "plt.legend(loc=\"best\")\n", + "plt.title('run44498, W4f, sidebands comparison')\n", + "plt.plot(x6,y6,'bx', label='$2^{nd}$ order sideband')\n", + "plt.plot(x6,out6.best_fit,'b', label=\"FWHM = {:.3f} ps\".format(out6.values['fwhm']))\n", + "plt.legend(loc=\"best\")\n", + "plt.xlabel(\"pump probe delay [ps]\")\n", + "plt.ylabel(\"Intensity [cts/s]\")\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "8da4db1f", + "metadata": {}, + "source": [ + "## Comparison of the BAM correction effect" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5c0ce10c", + "metadata": {}, + "outputs": [], + "source": [ + "fig,ax=plt.subplots(2,2,figsize=(9,7),layout=\"constrained\")\n", + "\n", + "plt.axes(ax[0,0])\n", + "res.plot(cmap='terrain', robust=True)\n", + "plt.title(\"W4f, no bam correction\")\n", + "\n", + "plt.axes(ax[0,1])\n", + "plt.plot(x1,y1,'rx',label='integrated intensity 1. order')\n", + "plt.plot(x1,out1.best_fit,'r',label='1. order fit, FWHM = {:.3f} ps'.format(out1.values['fwhm']))\n", + "plt.plot(x2,y2,'bx',label='integrated intensity 2. order')\n", + "plt.plot(x2,out2.best_fit,'b',label='2. order fit, FWHM = {:.3f} ps'.format(out2.values['fwhm']))\n", + "plt.legend(loc=1) \n", + "plt.title(\"Sidebands without bam correction\")\n", + "\n", + "plt.axes(ax[1,0])\n", + "res_corr.sel(delayStage=slice(-1.6,1.5)).plot(robust=True,cmap='terrain')\n", + "plt.title(\"W4f, with bam correction\")\n", + "\n", + "plt.axes(ax[1,1])\n", + "plt.plot(x5,y5,'rx',label='integrated intensity 1. order')\n", + "plt.plot(x5,out5.best_fit,'r',label='1. order fit, FWHM = {:.3f} ps'.format(out5.values['fwhm']))\n", + "plt.plot(x6,y6,'bx',label='integrated intensity 2. order')\n", + "plt.plot(x6,out6.best_fit,'b',label='2. order fit, FWHM = {:.3f} ps'.format(out6.values['fwhm']))\n", + "plt.legend(loc=1)\n", + "plt.title(\"Sidebands with bam correction\")\n", + "\n", + "fig.suptitle(f'Run {run_number}: Effect of BAM correction',fontsize='22')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "16fba27b", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "python3", + "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.9.19" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/tutorial/11_hextof_workflow_trXPS_energy_calibration_using_SB.ipynb b/tutorial/11_hextof_workflow_trXPS_energy_calibration_using_SB.ipynb new file mode 100644 index 00000000..3c89547c --- /dev/null +++ b/tutorial/11_hextof_workflow_trXPS_energy_calibration_using_SB.ipynb @@ -0,0 +1,511 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "916b6dc6", + "metadata": {}, + "source": [ + "# Tutorial for trXPS for energy calibration using core level side-bands" + ] + }, + { + "cell_type": "markdown", + "id": "61c9f989", + "metadata": {}, + "source": [ + "## Preparation" + ] + }, + { + "cell_type": "markdown", + "id": "aacad72f", + "metadata": {}, + "source": [ + "### Import necessary libraries" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5e6dcd73", + "metadata": {}, + "outputs": [], + "source": [ + "%load_ext autoreload\n", + "%autoreload 2\n", + "\n", + "from pathlib import Path\n", + "import os\n", + "\n", + "from sed import SedProcessor\n", + "from sed.dataset import dataset\n", + "import numpy as np\n", + "\n", + "%matplotlib widget\n", + "import matplotlib.pyplot as plt\n", + "\n", + "### for automatic peak finding\n", + "from scipy.signal import find_peaks" + ] + }, + { + "cell_type": "markdown", + "id": "2c147a5f-d84d-448a-9685-056035e464ed", + "metadata": {}, + "source": [ + "### Get data paths\n", + "\n", + "If it is your beamtime, you can read the raw data and write to the processed directory. For the public data, you can not write to the processed directory.\n", + "\n", + "The paths are such that if you are on Maxwell, it uses those. Otherwise, data is downloaded in the current directory from Zenodo:\n", + "https://zenodo.org/records/12609441" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7435190c-02b4-4e4a-82ca-3a8569751eb8", + "metadata": {}, + "outputs": [], + "source": [ + "beamtime_dir = \"/asap3/flash/gpfs/pg2/2023/data/11019101\" # on Maxwell\n", + "if os.path.exists(beamtime_dir) and os.access(beamtime_dir, os.R_OK):\n", + " path = beamtime_dir + \"/raw/hdf/offline/fl1user3\"\n", + " buffer_path = beamtime_dir + \"/processed/tutorial/\"\n", + "else:\n", + " # data_path can be defined and used to store the data in a specific location\n", + " dataset.get(\"W110\") # Put in Path to a storage of at least 10 GByte free space.\n", + " path = dataset.dir\n", + " buffer_path = path + \"/processed/\"" + ] + }, + { + "cell_type": "markdown", + "id": "92b884c6-416c-43dc-9ea4-59aeb04cd045", + "metadata": {}, + "source": [ + "### Config setup\n", + "\n", + "Here, we get the path to the config file and set up the relevant directories. This can also be done directly in the config file." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3cfe124d-4758-47be-afb2-d9d735575c0a", + "metadata": {}, + "outputs": [], + "source": [ + "# pick the default configuration file for hextof@FLASH\n", + "config_file = Path('../sed/config/flash_example_config.yaml')\n", + "assert config_file.exists()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "25135c58-41ae-467d-97fa-b726075e687d", + "metadata": {}, + "outputs": [], + "source": [ + "# here we setup a dictionary that will be used to override the path configuration\n", + "config_override = {\n", + " \"core\": {\n", + " \"beamtime_id\": 11019101,\n", + " \"paths\": {\n", + " \"data_raw_dir\": path,\n", + " \"data_parquet_dir\": buffer_path\n", + " },\n", + " },\n", + "}" + ] + }, + { + "cell_type": "markdown", + "id": "a47990f8", + "metadata": {}, + "source": [ + "## Reference calibration from a bias series" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "df0520c7", + "metadata": {}, + "outputs": [], + "source": [ + "sp_44455 = SedProcessor(runs=[44455], config=config_override, system_config=config_file)\n", + "sp_44455.add_jitter()\n", + "sp_44455.align_dld_sectors()" + ] + }, + { + "cell_type": "markdown", + "id": "3728fdd2", + "metadata": {}, + "source": [ + "### find calibration parameters\n", + "We now will fit the tof-energy relation. This is done by finding the maxima of a peak in the tof spectrum, and then fitting the square root relation to obtain the calibration parameters." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c7156708", + "metadata": {}, + "outputs": [], + "source": [ + "axes = ['sampleBias','dldTimeSteps']\n", + "bins = [4, 250]\n", + "ranges = [[77.5,81.5], [4050,4500]]\n", + "res = sp_44455.compute(bins=bins, axes=axes, ranges=ranges)\n", + "sp_44455.load_bias_series(binned_data=res)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e30cfea9", + "metadata": {}, + "outputs": [], + "source": [ + "ranges=(4120, 4200)\n", + "ref_id=0\n", + "sp_44455.find_bias_peaks(ranges=ranges, ref_id=ref_id, apply=True)" + ] + }, + { + "cell_type": "markdown", + "id": "d3e4f30d", + "metadata": {}, + "source": [ + "We offset the reference energy by the different in bias voltages between this run and run 44498" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bfe5f039", + "metadata": {}, + "outputs": [], + "source": [ + "sp_44455.calibrate_energy_axis(\n", + " ref_id=0,\n", + " ref_energy=-34.9,\n", + " method=\"lmfit\",\n", + " energy_scale='kinetic',\n", + " d={'value':1.0,'min': .7, 'max':1.0, 'vary':True},\n", + " t0={'value':5e-7, 'min': 1e-7, 'max': 1e-6, 'vary':True},\n", + " E0={'value': 0., 'min': -200, 'max': 100, 'vary': True},\n", + " verbose=True,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "11e375bd", + "metadata": {}, + "source": [ + "Now that we have the calibration parameters, we can generate the energy axis for each spectrum" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "626093ec", + "metadata": {}, + "outputs": [], + "source": [ + "sp_44455.save_energy_calibration(\"reference_calib.yaml\")" + ] + }, + { + "cell_type": "markdown", + "id": "98266c62-ab48-4746-96c8-2d47cf92c0e9", + "metadata": {}, + "source": [ + "### Now we can use those parameters and load our trXPS data using the additional config file" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2153366b-f27e-43cc-bd27-af005e7ea5d6", + "metadata": {}, + "outputs": [], + "source": [ + "run_number = 44498\n", + "sp_44498 = SedProcessor(runs=[run_number], config=config_override, folder_config=\"reference_calib.yaml\", system_config=config_file, verbose=True)\n", + "sp_44498.add_jitter()\n", + "sp_44498.append_energy_axis()" + ] + }, + { + "cell_type": "markdown", + "id": "ca21d0dd", + "metadata": {}, + "source": [ + "And bin an energy spectrum for reference" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "217395ee", + "metadata": {}, + "outputs": [], + "source": [ + "axes = ['energy']\n", + "ranges = [[-37.5,-27.5]]\n", + "bins = [200]\n", + "res_ref = sp_44498.compute(bins=bins, axes=axes, ranges=ranges)\n", + "\n", + "plt.figure()\n", + "res_ref.plot()" + ] + }, + { + "cell_type": "markdown", + "id": "0d735518", + "metadata": {}, + "source": [ + "## Energy calibration using side-band peaks" + ] + }, + { + "cell_type": "markdown", + "id": "6c2c8918-322b-4429-8ee9-cc8d032c82d3", + "metadata": {}, + "source": [ + "### Visualize trXPS data bin in the dldTimeSteps and the corrected delay axis to prepare for energy calibration using SB\n", + "We now prepare for an alternative energy calibration based on the side-bands of the time-dependent dataset. This is e.g. helpful if no bias series has been obtained." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1260af8e", + "metadata": {}, + "outputs": [], + "source": [ + "run_number = 44498\n", + "sp_44498 = SedProcessor(runs=[run_number], config=config_override, system_config=config_file, verbose=True)\n", + "sp_44498.add_jitter()" + ] + }, + { + "cell_type": "markdown", + "id": "b0c90a6f", + "metadata": {}, + "source": [ + "### We correct delay stage, t0 position and BAM (see previous tutorial)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ed13e712", + "metadata": {}, + "outputs": [], + "source": [ + "sp_44498.add_delay_offset(\n", + " constant=-1448, # this is time zero position determined from side band fit\n", + " flip_delay_axis=True, # invert the direction of the delay axis\n", + " columns=['bam'], # use the bam to offset the values\n", + " weights=[-0.001], # bam is in fs, delay in ps\n", + " preserve_mean=True # preserve the mean of the delay axis to keep t0 position\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "10c45742-f7f4-4c80-bdcf-fc04a5158b77", + "metadata": {}, + "outputs": [], + "source": [ + "axes = ['dldTimeSteps', 'delayStage']\n", + "ranges = [[3900,4200], [-1.5,1.5]]\n", + "bins = [100,60]\n", + "res_corr = sp_44498.compute(bins=bins, axes=axes, ranges=ranges, normalize_to_acquisition_time=\"delayStage\")\n", + "\n", + "fig,ax = plt.subplots(1,2,figsize=(8,3), layout='constrained')\n", + "fig.suptitle(f\"Run {run_number}: W 4f, side bands\")\n", + "res_corr.plot(robust=True, ax=ax[0], cmap='terrain')\n", + "ax[0].set_title('raw')\n", + "bg = res_corr.sel(delayStage=slice(-1.3,-1.0)).mean('delayStage')\n", + "(res_corr-bg).plot(robust=True, ax=ax[1])\n", + "ax[1].set_title('difference')" + ] + }, + { + "cell_type": "markdown", + "id": "423a4936-38e6-4432-bbee-b52bd655fe03", + "metadata": {}, + "source": [ + "### Automatically extract number and position of peaks in the ROI around t0" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "310b0737-4d28-4b76-a725-7d84a16355f7", + "metadata": {}, + "outputs": [], + "source": [ + "# binned data\n", + "roi = slice(3980, 4160)\n", + "delay = slice(-0.5,0.5)\n", + "data = res_corr.sel(dldTimeSteps = roi, delayStage=delay).sum('delayStage')\n", + "distance = 7\n", + "peaks, _ = find_peaks(data, height=None, distance=distance)\n", + "\n", + "p1SB = data[peaks]['dldTimeSteps'][0]\n", + "W4f5 = data[peaks]['dldTimeSteps'][1]\n", + "m1SB = data[peaks]['dldTimeSteps'][2]\n", + "W4f7 = data[peaks]['dldTimeSteps'][3]\n", + "mm1SB = data[peaks]['dldTimeSteps'][4]\n", + "plt.figure()\n", + "data.plot()\n", + "plt.scatter(data[peaks]['dldTimeSteps'], data[peaks], c='r')#, \"x\")\n", + "plt.vlines([p1SB-7,p1SB+7], 0, 150, color='violet', linestyles='dashed', label='$1^{st}$ order SB')\n", + "plt.vlines([W4f5-7,W4f5+7], 0, 150, color='b', linestyles='dashed', label='W 4f 7/2')\n", + "plt.vlines([m1SB-7,m1SB+7], 0, 150, color='g', linestyles='dashed', label='$-1^{st}$ order SB')\n", + "plt.vlines([W4f7-7,W4f7+7], 0, 150, color='r', linestyles='dashed', label='W 4f 5/2')\n", + "plt.vlines([mm1SB-7,mm1SB+7], 0, 150, color='orange', linestyles='dashed', label='$2nd -1^{st}$ order SB')\n", + "plt.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "f2bee3d0-0ad9-4ed7-8356-16333cf48774", + "metadata": {}, + "source": [ + "### find calibration parameters\n", + "We now will fit the tof-energy relation. This is done using the maxima of a peak in the ToF spectrum and the known kinetic energy of those peaks (kinetic energy of e.g. W4f peaks (-31.4 and -33.6 eV) and their SB of different orders accounting energy of pump beam of 1030 nm = 1.2 eV. The calibration parameters are obtained by fitting the square root relation." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "10752d1f-53c1-4871-b2fd-65a0101aba9a", + "metadata": {}, + "outputs": [], + "source": [ + "### Kinetic energy of w4f peaks and their SB\n", + "ref_energy = -31.4\n", + "ref_id = 1\n", + "sp_44498.ec.biases = -1*np.array([-30.2,-31.4,-32.6,-33.6,-34.8])\n", + "sp_44498.ec.peaks = np.expand_dims(data[peaks]['dldTimeSteps'].data,1)\n", + "sp_44498.ec.tof = res_corr.dldTimeSteps.data\n", + "\n", + "sp_44498.calibrate_energy_axis(\n", + " ref_id=ref_id,\n", + " ref_energy=ref_energy,\n", + " method=\"lmfit\",\n", + " d={'value':1.0,'min': .8, 'max':1.0, 'vary':True},\n", + " t0={'value':5e-7, 'min': 1e-7, 'max': 1e-6, 'vary':True},\n", + " E0={'value': -100., 'min': -200, 'max': 15, 'vary': True},\n", + " labels=\"\",\n", + " verbose=True)" + ] + }, + { + "cell_type": "markdown", + "id": "4052d629-1178-4248-a945-d60a6ff34bf3", + "metadata": {}, + "source": [ + "### Append energy axis into a data frame, bin and visualize data in the calibrated energy and corrected delay axis " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d37a8e16-d91c-4b93-be79-4d4f5ddd2dd0", + "metadata": {}, + "outputs": [], + "source": [ + "sp_44498.append_energy_axis()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4b78aab7-7a04-46f8-9e68-05f595758b7c", + "metadata": {}, + "outputs": [], + "source": [ + "axes = ['energy', 'delayStage']\n", + "ranges = [[-37.5,-27.5], [-1.5,1.5]]\n", + "bins = [200,60]\n", + "res_corr = sp_44498.compute(bins=bins, axes=axes, ranges=ranges, normalize_to_acquisition_time=\"delayStage\")\n", + "\n", + "fig,ax = plt.subplots(1,2,figsize=(8,3), layout='constrained')\n", + "fig.suptitle(f\"Run {run_number}: W 4f, side bands\")\n", + "res_corr.plot(robust=True, ax=ax[0], cmap='terrain')\n", + "ax[0].set_title('raw')\n", + "bg = res_corr.sel(delayStage=slice(-1.3,-1.0)).mean('delayStage')\n", + "(res_corr-bg).plot(robust=True, ax=ax[1])\n", + "ax[1].set_title('difference')" + ] + }, + { + "cell_type": "markdown", + "id": "a66fac5e", + "metadata": {}, + "source": [ + "## Compare to reference\n", + "While this calibration methods gives a reasonable approximation to the energy axis, there are some deviations to the bias method, so it should be used with care" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d83d0eb6-34f2-4eec-a568-3a1767d8e705", + "metadata": {}, + "outputs": [], + "source": [ + "axes = ['energy']\n", + "ranges = [[-37.5,-27.5]]\n", + "bins = [200]\n", + "res_1D = sp_44498.compute(bins=bins, axes=axes, ranges=ranges)\n", + "\n", + "plt.figure()\n", + "(res_ref/res_ref.max()).plot(label=\"bias series calibration\")\n", + "(res_1D/res_1D.max()).plot(label=\"side band calibration\")\n", + "plt.legend()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a26d270d", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "python3", + "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.9.19" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/tutorial/9_hextof_workflow_trXPD.ipynb b/tutorial/9_hextof_workflow_trXPD.ipynb new file mode 100644 index 00000000..90f3ef63 --- /dev/null +++ b/tutorial/9_hextof_workflow_trXPD.ipynb @@ -0,0 +1,555 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "733ed22f", + "metadata": {}, + "source": [ + "# Tutorial for trXPD for the HEXTOF instrument at FLASH with background normalization" + ] + }, + { + "cell_type": "markdown", + "id": "d7eaa93d", + "metadata": {}, + "source": [ + "## Preparation" + ] + }, + { + "cell_type": "markdown", + "id": "fc871acf", + "metadata": {}, + "source": [ + "### Import necessary libraries" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "368cf206", + "metadata": {}, + "outputs": [], + "source": [ + "%load_ext autoreload\n", + "%autoreload 2\n", + "\n", + "from pathlib import Path\n", + "import os\n", + "\n", + "from sed import SedProcessor\n", + "from sed.dataset import dataset\n", + "import xarray as xr\n", + "\n", + "%matplotlib widget\n", + "import matplotlib.pyplot as plt\n", + "from scipy.ndimage import gaussian_filter" + ] + }, + { + "cell_type": "markdown", + "id": "d3ed214a", + "metadata": {}, + "source": [ + "### Get data paths\n", + "\n", + "If it is your beamtime, you can access both read the raw data and write to processed directory. For the public data, you can not write to processed directory.\n", + "\n", + "The paths are such that if you are on Maxwell, it uses those. Otherwise data is downloaded in current directory from Zenodo:\n", + "https://zenodo.org/records/12609441" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "195737a6", + "metadata": {}, + "outputs": [], + "source": [ + "beamtime_dir = \"/asap3/flash/gpfs/pg2/2023/data/11019101\" # on Maxwell\n", + "if os.path.exists(beamtime_dir) and os.access(beamtime_dir, os.R_OK):\n", + " path = beamtime_dir + \"/raw/hdf/offline/fl1user3\"\n", + " buffer_path = beamtime_dir + \"/processed/tutorial/\"\n", + "else:\n", + " # data_path can be defined and used to store the data in a specific location\n", + " dataset.get(\"W110\") # Put in Path to a storage of at least 10 GByte free space.\n", + " path = dataset.dir\n", + " buffer_path = path + \"/processed/\"" + ] + }, + { + "cell_type": "markdown", + "id": "30effaac", + "metadata": {}, + "source": [ + "### Config setup\n", + "\n", + "Here we get the path to the config file and setup the relevant directories. This can also be done directly in the config file." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2ca5745a", + "metadata": {}, + "outputs": [], + "source": [ + "# pick the default configuration file for hextof@FLASH\n", + "config_file = Path('../sed/config/flash_example_config.yaml')\n", + "assert config_file.exists()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "347d338c", + "metadata": {}, + "outputs": [], + "source": [ + "# here we setup a dictionary that will be used to override the path configuration\n", + "config_override = {\n", + " \"core\": {\n", + " \"beamtime_id\": 11019101,\n", + " \"paths\": {\n", + " \"data_raw_dir\": path,\n", + " \"data_parquet_dir\": buffer_path\n", + " },\n", + " },\n", + "}" + ] + }, + { + "cell_type": "markdown", + "id": "ee483947", + "metadata": {}, + "source": [ + "### Prepare Energy Calibration\n", + "Instead of making completely new energy calibration we can take existing values from the calibration made in the previous tutorial. This allows us to calibrate the conversion between the digital values of the dld and the energy." + ] + }, + { + "cell_type": "markdown", + "id": "7417ab20", + "metadata": {}, + "source": [ + " For this we need to add all those parameters as a dictionary and use them during creation of the processor object." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cdcbd374", + "metadata": {}, + "outputs": [], + "source": [ + "en_cal_config = {\n", + " 'energy': {\n", + " 'calibration': {\n", + " 'E0': -54.971004271795664,\n", + " 'creation_date': 1718801358.232129,\n", + " 'd': 0.8096677238144319,\n", + " 'energy_scale': 'kinetic',\n", + " 't0': 4.0148196706891397e-07,\n", + " 'calib_type': 'fit',\n", + " 'fit_function': '(a0/(x0-a1))**2 + a2',\n", + " 'coefficients': ([ 8.09667724e-01, 4.01481967e-07, -5.49710043e+01]),\n", + " 'axis': 0},\n", + " 'tof': None,\n", + " 'offsets': {\n", + " 'constant': -76.5,\n", + " 'creation_date': 1718801360.817963,\n", + " 'monochromatorPhotonEnergy': {'preserve_mean': True,'reduction': None,'weight': -1},\n", + " 'sampleBias': {'preserve_mean': False, 'reduction': None, 'weight': 1},\n", + " 'tofVoltage': {'preserve_mean': True, 'reduction': None, 'weight': -1}}}}" + ] + }, + { + "cell_type": "markdown", + "id": "f070b3bb", + "metadata": {}, + "source": [ + "## Read data\n", + "Now we can use those parameters and load our trXPD data using additional config file" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "39b316d8", + "metadata": {}, + "outputs": [], + "source": [ + "run_number = 44498\n", + "sp_44498 = SedProcessor(runs=[run_number], folder_config=en_cal_config, config=config_override, system_config=config_file, verbose=True)\n", + "sp_44498.add_jitter()" + ] + }, + { + "cell_type": "markdown", + "id": "a2afd119", + "metadata": {}, + "source": [ + "We can inspect dataframe right after data readout" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "167aaf25", + "metadata": {}, + "outputs": [], + "source": [ + "sp_44498.dataframe.head()" + ] + }, + { + "cell_type": "markdown", + "id": "d56f2f83", + "metadata": {}, + "source": [ + "Now we will do energy calibration, add energy offset, jittering and dld sectors alignment" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c1d12db4", + "metadata": {}, + "outputs": [], + "source": [ + "sp_44498.align_dld_sectors()\n", + "sp_44498.append_energy_axis()\n", + "sp_44498.add_energy_offset()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c7addaa4", + "metadata": {}, + "outputs": [], + "source": [ + "sp_44498.attributes.metadata['energy_calibration']" + ] + }, + { + "cell_type": "markdown", + "id": "c97643f5", + "metadata": {}, + "source": [ + "We can do the SASE jitter correction, using information from the bam column and do calibration of the pump-probe delay axis, we need to shift the delay stage values to center the pump-probe-time overlap time zero." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "514e9956", + "metadata": {}, + "outputs": [], + "source": [ + "sp_44498.add_delay_offset(\n", + " constant=-1448, # this is time zero position determined from side band fit\n", + " flip_delay_axis=True, # invert the direction of the delay axis\n", + " columns=['bam'], # use the bam to offset the values\n", + " weights=[-0.001], # bam is in fs, delay in ps\n", + " preserve_mean=True # preserve the mean of the delay axis to keep t0 position\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "6551542d", + "metadata": {}, + "source": [ + "### bin in the calibrated energy and corrected delay axis\n", + "Visualize trXPS data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "88055c6d", + "metadata": {}, + "outputs": [], + "source": [ + "axes = ['energy', 'delayStage']\n", + "ranges = [[-37.5,-27.5], [-1.5,1.5]]\n", + "bins = [200,60]\n", + "res_corr = sp_44498.compute(bins=bins, axes=axes, ranges=ranges, normalize_to_acquisition_time=\"delayStage\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b02865c8", + "metadata": {}, + "outputs": [], + "source": [ + "fig,ax = plt.subplots(1,2,figsize=(8,3), layout='constrained')\n", + "fig.suptitle(f\"Run {run_number}: W 4f, side bands\")\n", + "res_corr.plot(robust=True, ax=ax[0], cmap='terrain')\n", + "ax[0].set_title('raw')\n", + "bg = res_corr.sel(delayStage=slice(-1.3,-1.0)).mean('delayStage')\n", + "(res_corr-bg).plot(robust=True, ax=ax[1])\n", + "ax[1].set_title('difference')" + ] + }, + { + "cell_type": "markdown", + "id": "43e84eff", + "metadata": {}, + "source": [ + "## XPD from W4f core level\n", + "\n", + "Now we can bin not only in energy but also in both momentum directions to get XPD patterns of different core level line of tungsten." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b768cfc7", + "metadata": {}, + "outputs": [], + "source": [ + "axes = ['energy', 'dldPosX', 'dldPosY']\n", + "ranges = [[-38,-28], [420,900], [420,900]]\n", + "bins = [100,240,240]\n", + "res_kx_ky = sp_44498.compute(bins=bins, axes=axes, ranges=ranges)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "73aed405", + "metadata": {}, + "outputs": [], + "source": [ + "## EDC and integration region for XPD\n", + "plt.figure()\n", + "res_kx_ky.mean(('dldPosX', 'dldPosY')).plot()\n", + "plt.vlines([-30.3,-29.9], 0, 2.4, color='r', linestyles='dashed')\n", + "plt.vlines([-31.4,-31.2], 0, 2.4, color='orange', linestyles='dashed')\n", + "plt.vlines([-33.6,-33.4], 0, 2.4, color='g', linestyles='dashed')\n", + "plt.vlines([-37.0,-36.0], 0, 2.4, color='b', linestyles='dashed')\n", + "plt.title('EDC and integration regions for XPD')\n", + "plt.show()\n", + "\n", + "## XPD plots\n", + "fig,ax = plt.subplots(2,2,figsize=(9,7), layout='constrained')\n", + "res_kx_ky.sel(energy=slice(-30.3,-29.9)).mean('energy').plot(robust=True, ax=ax[0,0], cmap='terrain')\n", + "ax[0,0].set_title(\"XPD of $1^{st}$ order sidebands\")\n", + "res_kx_ky.sel(energy=slice(-31.4,-31.2)).mean('energy').plot(robust=True, ax=ax[0,1], cmap='terrain')\n", + "ax[0,1].set_title(\"XPD of W4f 7/2 peak\")\n", + "res_kx_ky.sel(energy=slice(-33.6,-33.4)).mean('energy').plot(robust=True, ax=ax[1,0], cmap='terrain')\n", + "ax[1,0].set_title(\"XPD of W4f 5/2 peak\")\n", + "res_kx_ky.sel(energy=slice(-37.0,-36.0)).mean('energy').plot(robust=True, ax=ax[1,1], cmap='terrain')\n", + "ax[1,1].set_title(\"XPD of W5p 3/2 peak\")" + ] + }, + { + "cell_type": "markdown", + "id": "c30cdbd0", + "metadata": {}, + "source": [ + "As we can see there is some structure visible, but it looks very similar to each other.\n", + "We probably have to do some normalization to remove the detector structure/artefacts.\n", + "The best option is to divide by a flat-field image. The flat-field image can be obtained from a sample that shows no structure under identical measurement conditions.\n", + "Unfortunately, we don't have such a flat-field image.\n", + "\n", + "In this case, we can make a flat-field image from the actual dataset using several different approaches.\n", + "\n", + "As a first option, we can integrate in energy over the whole region and use this image as a background.\n", + "Additionally, we introduce the Gaussian Blur for comparison." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2e4a7faf", + "metadata": {}, + "outputs": [], + "source": [ + "## Background image\n", + "bgd = res_kx_ky.mean(('energy'))\n", + "\n", + "## Apply Gaussian Blur to background image\n", + "bgd_blur = xr.apply_ufunc(gaussian_filter, bgd, 15)\n", + "\n", + "fig,ax = plt.subplots(1,2,figsize=(9,4), layout='constrained')\n", + "bgd.plot(robust=True, cmap='terrain', ax=ax[0])\n", + "ax[0].set_title('Background image')\n", + "bgd_blur.plot(cmap='terrain', ax=ax[1])\n", + "ax[1].set_title('Gaussian Blur of background image')\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6fa01c0a", + "metadata": {}, + "outputs": [], + "source": [ + "## XPD normalized by background image\n", + "fig,ax = plt.subplots(2,2,figsize=(9,7), layout='constrained')\n", + "(res_kx_ky/bgd).sel(energy=slice(-30.3,-29.9)).mean('energy').plot(robust=True, ax=ax[0,0], cmap='terrain')\n", + "(res_kx_ky/bgd).sel(energy=slice(-31.4,-31.2)).mean('energy').plot(robust=True, ax=ax[0,1], cmap='terrain')\n", + "(res_kx_ky/bgd).sel(energy=slice(-33.6,-33.4)).mean('energy').plot(robust=True, ax=ax[1,0], cmap='terrain')\n", + "(res_kx_ky/bgd).sel(energy=slice(-37.0,-36.0)).mean('energy').plot(robust=True, ax=ax[1,1], cmap='terrain')\n", + "fig.suptitle(f'Run {run_number}: XPD patterns after background normalization',fontsize='18')\n", + "\n", + "## XPD normalized by Gaussian-blurred background image\n", + "fig,ax = plt.subplots(2,2,figsize=(9,7), layout='constrained')\n", + "(res_kx_ky/bgd_blur).sel(energy=slice(-30.3,-29.9)).mean('energy').plot(robust=True, ax=ax[0,0], cmap='terrain')\n", + "(res_kx_ky/bgd_blur).sel(energy=slice(-31.4,-31.2)).mean('energy').plot(robust=True, ax=ax[0,1], cmap='terrain')\n", + "(res_kx_ky/bgd_blur).sel(energy=slice(-33.6,-33.4)).mean('energy').plot(robust=True, ax=ax[1,0], cmap='terrain')\n", + "(res_kx_ky/bgd_blur).sel(energy=slice(-37.0,-36.0)).mean('energy').plot(robust=True, ax=ax[1,1], cmap='terrain')\n", + "fig.suptitle(f'Run {run_number}: XPD patterns after Gaussian-blurred background normalization',fontsize='18')\n", + "\n", + "## XPD normalized by Gaussian-blurred background image and blurred to improve contrast\n", + "fig,ax = plt.subplots(2,2,figsize=(9,7), layout='constrained')\n", + "(xr.apply_ufunc(gaussian_filter, res_kx_ky/bgd_blur, 1)).sel(energy=slice(-30.3,-29.9)).mean('energy').plot(robust=True, ax=ax[0,0], cmap='terrain')\n", + "(xr.apply_ufunc(gaussian_filter, res_kx_ky/bgd_blur, 1)).sel(energy=slice(-31.4,-31.2)).mean('energy').plot(robust=True, ax=ax[0,1], cmap='terrain')\n", + "(xr.apply_ufunc(gaussian_filter, res_kx_ky/bgd_blur, 1)).sel(energy=slice(-33.6,-33.4)).mean('energy').plot(robust=True, ax=ax[1,0], cmap='terrain')\n", + "(xr.apply_ufunc(gaussian_filter, res_kx_ky/bgd_blur, 1)).sel(energy=slice(-37.0,-36.0)).mean('energy').plot(robust=True, ax=ax[1,1], cmap='terrain')\n", + "fig.suptitle(f'Run {run_number}: resulting Gaussian-blurred XPD patterns',fontsize='18')" + ] + }, + { + "cell_type": "markdown", + "id": "6f415a21", + "metadata": {}, + "source": [ + "Sometimes, after this division, you may not be happy with intensity distribution. Thus, other option for background correction is to duplicate the XPD pattern, apply large Gaussian blurring that eliminates the fine structures in the XPD pattern. Then divide the XPD pattern by its blurred version. This process sometimes enhances the visibility of the fine structures a lot." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5e5607ea", + "metadata": {}, + "outputs": [], + "source": [ + "## XPD normalized by Gaussian-blurred background image\n", + "\n", + "### Define integration regions for XPD\n", + "SB = res_kx_ky.sel(energy=slice(-30.3,-29.9)).mean('energy')\n", + "W_4f_7 = res_kx_ky.sel(energy=slice(-31.4,-31.2)).mean('energy')\n", + "W_4f_5 = res_kx_ky.sel(energy=slice(-33.6,-33.4)).mean('energy')\n", + "W_5p = res_kx_ky.sel(energy=slice(-37.0,-36.0)).mean('energy')\n", + "\n", + "### Make corresponding Gaussian Blur background\n", + "SB_blur = xr.apply_ufunc(gaussian_filter, SB, 15)\n", + "W_4f_7_blur = xr.apply_ufunc(gaussian_filter, W_4f_7, 15)\n", + "W_4f_5_blur = xr.apply_ufunc(gaussian_filter, W_4f_5, 15)\n", + "W_5p_blur = xr.apply_ufunc(gaussian_filter, W_5p, 15)\n", + "\n", + "### Visualize results\n", + "fig,ax = plt.subplots(2,2,figsize=(9,7), layout='constrained')\n", + "(SB/SB_blur).plot(robust=True, ax=ax[0,0], cmap='terrain')\n", + "(W_4f_7/W_4f_7_blur).plot(robust=True, ax=ax[0,1], cmap='terrain')\n", + "(W_4f_5/W_4f_5_blur).plot(robust=True, ax=ax[1,0], cmap='terrain')\n", + "(W_5p/W_5p_blur).plot(robust=True, ax=ax[1,1], cmap='terrain')\n", + "fig.suptitle(f'Run {run_number}: XPD patterns after Gaussian Blur normalization',fontsize='18')\n", + "\n", + "### Apply Gaussian Blur to resulted images to improve contrast\n", + "SB_norm = xr.apply_ufunc(gaussian_filter, SB/SB_blur, 1)\n", + "W_4f_7_norm = xr.apply_ufunc(gaussian_filter, W_4f_7/W_4f_7_blur, 1)\n", + "W_4f_5_norm = xr.apply_ufunc(gaussian_filter, W_4f_5/W_4f_5_blur, 1)\n", + "W_5p_norm = xr.apply_ufunc(gaussian_filter, W_5p/W_5p_blur, 1)\n", + "\n", + "### Visualize results\n", + "fig,ax = plt.subplots(2,2,figsize=(9,7), layout='constrained')\n", + "SB_norm.plot(robust=True, ax=ax[0,0], cmap='terrain')\n", + "W_4f_7_norm.plot(robust=True, ax=ax[0,1], cmap='terrain')\n", + "W_4f_5_norm.plot(robust=True, ax=ax[1,0], cmap='terrain')\n", + "W_5p_norm.plot(robust=True, ax=ax[1,1], cmap='terrain')\n", + "fig.suptitle(f'Run {run_number}: XPD patterns after Gauss Blur normalization',fontsize='18') " + ] + }, + { + "cell_type": "markdown", + "id": "94d41b4f", + "metadata": {}, + "source": [ + "Third option for background normalization is to use the simultaneously acquired pre-core level region.\n", + "As an example for W4f 7/2 peak, we define a region on the high energy side of it and integrate in energy to use as a background" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6c8513b5", + "metadata": {}, + "outputs": [], + "source": [ + "### Define peak and background region on the high energy side of the peak\n", + "W_4f_7 = res_kx_ky.sel(energy=slice(-31.4,-31.2)).mean('energy')\n", + "W_4f_7_bgd = res_kx_ky.sel(energy=slice(-32.0,-31.8)).mean('energy')\n", + "\n", + "### Make normalization by background, add Gaussian Blur to the resulting image\n", + "W_4f_7_nrm1 = W_4f_7/(W_4f_7_bgd+W_4f_7_bgd.max()*0.00001)\n", + "W_4f_7_nrm1_blur = xr.apply_ufunc(gaussian_filter, W_4f_7_nrm1, 1)\n", + "\n", + "### Add Gaussian Blur to the background image, normalize by it and add Gaussian Blur to the resulting image\n", + "W_4f_7_bgd_blur = xr.apply_ufunc(gaussian_filter, W_4f_7_bgd, 15)\n", + "W_4f_7_nrm2 = W_4f_7/W_4f_7_bgd_blur\n", + "W_4f_7_nrm2_blur = xr.apply_ufunc(gaussian_filter, W_4f_7_nrm2, 1)\n", + "\n", + "### Visualize all steps\n", + "fig,ax = plt.subplots(4,2,figsize=(9,10), layout='constrained')\n", + "W_4f_7.plot(robust=True, ax=ax[0,0], cmap='terrain')\n", + "W_4f_7_bgd.plot(robust=True, ax=ax[0,1], cmap='terrain')\n", + "W_4f_7_nrm1.plot(robust=True, ax=ax[1,0], cmap='terrain')\n", + "W_4f_7_nrm1_blur.plot(robust=True, ax=ax[1,1], cmap='terrain')\n", + "W_4f_7_bgd_blur.plot(robust=True, ax=ax[2,0], cmap='terrain')\n", + "W_4f_7_nrm2.plot(robust=True, ax=ax[2,1], cmap='terrain')\n", + "W_4f_7_nrm2_blur.plot(robust=True, ax=ax[3,0], cmap='terrain')\n", + "fig.suptitle(f'Run {run_number}: XPD patterns of W4f7/2 with pre-core level normalization',fontsize='18') " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "807883fd", + "metadata": {}, + "outputs": [], + "source": [ + "fig,ax = plt.subplots(1,3,figsize=(9,3), layout='constrained')\n", + "(xr.apply_ufunc(gaussian_filter, res_kx_ky/bgd_blur, 1)).sel(energy=slice(-31.4,-31.2)).mean('energy').plot(robust=True, ax=ax[0], cmap='terrain')\n", + "W_4f_7_norm.plot(robust=True, ax=ax[1], cmap='terrain')\n", + "W_4f_7_nrm2_blur.plot(robust=True, ax=ax[2], cmap='terrain')\n", + "fig.suptitle(f'Run {run_number}: comparison of different normalizations\\nof XPD pattern for W4f 7/2 peak with Gaussian Blur',fontsize='18')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7654516c", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "python3", + "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.9.19" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}