Tutorial for trXPD for the HEXTOF instrument at FLASH with background normalization#
Preparation#
Import necessary libraries#
[1]:
%load_ext autoreload
%autoreload 2
from pathlib import Path
import os
from sed import SedProcessor
from sed.dataset import dataset
import xarray as xr
%matplotlib widget
import matplotlib.pyplot as plt
from scipy.ndimage import gaussian_filter
Get data paths#
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.
The paths are such that if you are on Maxwell, it uses those. Otherwise data is downloaded in current directory from Zenodo: https://zenodo.org/records/12609441
[2]:
beamtime_dir = "/asap3/flash/gpfs/pg2/2023/data/11019101" # on Maxwell
if os.path.exists(beamtime_dir) and os.access(beamtime_dir, os.R_OK):
path = beamtime_dir + "/raw/hdf/offline/fl1user3"
buffer_path = beamtime_dir + "/processed/tutorial/"
else:
# data_path can be defined and used to store the data in a specific location
dataset.get("W110") # Put in Path to a storage of at least 10 GByte free space.
path = dataset.dir
buffer_path = path + "/processed/"
INFO - Not downloading W110 data as it already exists at "/home/runner/work/sed/sed/docs/tutorial/datasets/W110".
Set 'use_existing' to False if you want to download to a new location.
INFO - Using existing data path for "W110": "/home/runner/work/sed/sed/docs/tutorial/datasets/W110"
INFO - W110 data is already present.
Config setup#
Here we get the path to the config file and setup the relevant directories. This can also be done directly in the config file.
[3]:
# pick the default configuration file for hextof@FLASH
config_file = Path('../src/sed/config/flash_example_config.yaml')
assert config_file.exists()
[4]:
# here we setup a dictionary that will be used to override the path configuration
config_override = {
"core": {
"beamtime_id": 11019101,
"paths": {
"raw": path,
"processed": buffer_path
},
},
}
Prepare Energy Calibration#
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.
For this we need to add all those parameters as a dictionary and use them during creation of the processor object.
[5]:
energy_cal = {
"energy": {
"calibration": {
"E0": -132.47100427179566,
"creation_date": '2024-11-30T20:47:03.305244',
"d": 0.8096677238144319,
"energy_scale": "kinetic",
"t0": 4.0148196706891397e-07,
},
"offsets":{
"constant": 1,
"creation_date": '2024-11-30T21:17:07.762199',
"columns": {
"monochromatorPhotonEnergy": {
"preserve_mean": True,
"weight": -1,
},
"tofVoltage": {
"preserve_mean": True,
"weight": -1,
},
},
},
},
}
Read data#
Now we can use those parameters and load our trXPD data using additional config file
[6]:
run_number = 44498
sp_44498 = SedProcessor(runs=[run_number], folder_config=energy_cal, config=config_override, system_config=config_file, verbose=True)
sp_44498.add_jitter()
INFO - System config loaded from: [/home/runner/work/sed/sed/docs/src/sed/config/flash_example_config.yaml]
INFO - Default config loaded from: [/opt/hostedtoolcache/Python/3.10.16/x64/lib/python3.10/site-packages/sed/config/default.yaml]
INFO - Reading files: 0 new files of 14 total.
loading complete in 0.13 s
INFO - add_jitter: Added jitter to columns ['dldPosX', 'dldPosY', 'dldTimeSteps'].
We can inspect dataframe right after data readout
[7]:
sp_44498.dataframe.head()
[7]:
| trainId | pulseId | electronId | dldPosX | dldPosY | dldTimeSteps | pulserSignAdc | bam | timeStamp | monochromatorPhotonEnergy | gmdBda | delayStage | sampleBias | tofVoltage | extractorVoltage | extractorCurrent | cryoTemperature | sampleTemperature | dldTimeBinSize | dldSectorID | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 1628022830 | 1 | 0 | 651.154674 | 895.154674 | 4595.154674 | 32919.0 | -6187.96875 | 1.677563e+09 | 116.858299 | 2.521457 | 1448.602051 | 72.995903 | 19.995356 | 6029.299805 | -0.073857 | 49.208 | 78.989998 | 0.020576 | 3 |
| 1 | 1628022830 | 1 | 1 | 650.795029 | 887.795029 | 4595.795029 | 32919.0 | -6187.96875 | 1.677563e+09 | 116.858299 | 2.521457 | 1448.602051 | 72.995903 | 19.995356 | 6029.299805 | -0.073857 | 49.208 | 78.989998 | 0.020576 | 0 |
| 2 | 1628022830 | 5 | 0 | 681.824157 | 671.824157 | 4422.824157 | 32914.0 | -6170.15625 | 1.677563e+09 | 116.858299 | 2.521457 | 1448.602051 | 72.995903 | 19.995356 | 6029.299805 | -0.073857 | 49.208 | 78.989998 | 0.020576 | 6 |
| 3 | 1628022830 | 5 | 1 | 684.877689 | 657.877689 | 4424.877689 | 32914.0 | -6170.15625 | 1.677563e+09 | 116.858299 | 2.521457 | 1448.602051 | 72.995903 | 19.995356 | 6029.299805 | -0.073857 | 49.208 | 78.989998 | 0.020576 | 3 |
| 4 | 1628022830 | 5 | 2 | 669.901034 | 686.901034 | 4423.901034 | 32914.0 | -6170.15625 | 1.677563e+09 | 116.858299 | 2.521457 | 1448.602051 | 72.995903 | 19.995356 | 6029.299805 | -0.073857 | 49.208 | 78.989998 | 0.020576 | 5 |
Now we will do energy calibration, add energy offset, jittering and dld sectors alignment
[8]:
sp_44498.align_dld_sectors()
sp_44498.append_energy_axis()
sp_44498.add_energy_offset()
INFO - Aligning 8s sectors of dataframe
INFO - Dask DataFrame Structure:
trainId pulseId electronId dldPosX dldPosY dldTimeSteps pulserSignAdc bam timeStamp monochromatorPhotonEnergy gmdBda delayStage sampleBias tofVoltage extractorVoltage extractorCurrent cryoTemperature sampleTemperature dldTimeBinSize dldSectorID
npartitions=14
uint32 int64 int64 float64 float64 float32 float32 float32 float64 float32 float32 float32 float32 float32 float32 float32 float32 float32 float32 int8
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
Dask Name: assign, 16 graph layers
INFO - Adding energy column to dataframe:
INFO - Using energy calibration parameters generated on 11/30/2024, 20:47:03
INFO - Dask DataFrame Structure:
trainId pulseId electronId dldPosX dldPosY dldTimeSteps pulserSignAdc bam timeStamp monochromatorPhotonEnergy gmdBda delayStage sampleBias tofVoltage extractorVoltage extractorCurrent cryoTemperature sampleTemperature dldTimeBinSize dldSectorID energy
npartitions=14
uint32 int64 int64 float64 float64 float32 float32 float32 float64 float32 float32 float32 float32 float32 float32 float32 float32 float32 float32 int8 float64
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
Dask Name: assign, 31 graph layers
INFO - Adding energy offset to dataframe:
INFO - Using energy offset parameters generated on 11/30/2024, 21:17:07
INFO - Energy offset parameters:
Constant: 1.0
Column[monochromatorPhotonEnergy]: Weight=-1.0, Preserve Mean: True, Reductions: None.
Column[tofVoltage]: Weight=-1.0, Preserve Mean: True, Reductions: None.
INFO - Dask DataFrame Structure:
trainId pulseId electronId dldPosX dldPosY dldTimeSteps pulserSignAdc bam timeStamp monochromatorPhotonEnergy gmdBda delayStage sampleBias tofVoltage extractorVoltage extractorCurrent cryoTemperature sampleTemperature dldTimeBinSize dldSectorID energy
npartitions=14
uint32 int64 int64 float64 float64 float32 float32 float32 float64 float32 float32 float32 float32 float32 float32 float32 float32 float32 float32 int8 float64
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
Dask Name: assign, 64 graph layers
[9]:
sp_44498.attributes.metadata['energy_calibration']
[9]:
{'applied': True,
'calibration': {'creation_date': datetime.datetime(2024, 11, 30, 20, 47, 3, 305244),
'd': 0.8096677238144319,
't0': 4.0148196706891397e-07,
'E0': -132.47100427179566,
'energy_scale': 'kinetic',
'calib_type': 'fit',
'fit_function': '(a0/(x0-a1))**2 + a2',
'coefficients': array([ 8.09667724e-01, 4.01481967e-07, -1.32471004e+02]),
'axis': 0.0},
'tof': 0.0}
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.
[10]:
sp_44498.add_delay_offset(
constant=-1448, # this is time zero position determined from side band fit
flip_delay_axis=True, # invert the direction of the delay axis
columns=['bam'], # use the bam to offset the values
weights=[-0.001], # bam is in fs, delay in ps
preserve_mean=True # preserve the mean of the delay axis to keep t0 position
)
INFO - Adding delay offset to dataframe:
INFO - Delay offset parameters:
Column[bam]: Weight=-0.001, Preserve Mean: True, Reductions: None.
Constant: -1448
Flip delay axis: True
INFO - Dask DataFrame Structure:
trainId pulseId electronId dldPosX dldPosY dldTimeSteps pulserSignAdc bam timeStamp monochromatorPhotonEnergy gmdBda delayStage sampleBias tofVoltage extractorVoltage extractorCurrent cryoTemperature sampleTemperature dldTimeBinSize dldSectorID energy
npartitions=14
uint32 int64 int64 float64 float64 float32 float32 float32 float64 float32 float32 float64 float32 float32 float32 float32 float32 float32 float32 int8 float64
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
Dask Name: assign, 84 graph layers
bin in the calibrated energy and corrected delay axis#
Visualize trXPS data
[11]:
axes = ['energy', 'delayStage']
ranges = [[-37.5,-27.5], [-1.5,1.5]]
bins = [200,60]
res_corr = sp_44498.compute(bins=bins, axes=axes, ranges=ranges, normalize_to_acquisition_time="delayStage")
INFO - Calculate normalization histogram for axis 'delayStage'...
[12]:
fig,ax = plt.subplots(1,2,figsize=(6,2.25), layout='constrained')
fig.suptitle(f"Run {run_number}: W 4f, side bands")
res_corr.plot(robust=True, ax=ax[0], cmap='terrain')
ax[0].set_title('raw')
bg = res_corr.sel(delayStage=slice(-1.3,-1.0)).mean('delayStage')
(res_corr-bg).plot(robust=True, ax=ax[1])
ax[1].set_title('difference')
[12]:
Text(0.5, 1.0, 'difference')
XPD from W4f core level#
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.
[13]:
axes = ['energy', 'dldPosX', 'dldPosY']
ranges = [[-38,-28], [420,900], [420,900]]
bins = [100,240,240]
res_kx_ky = sp_44498.compute(bins=bins, axes=axes, ranges=ranges)
[14]:
## EDC and integration region for XPD
plt.figure()
res_kx_ky.mean(('dldPosX', 'dldPosY')).plot()
plt.vlines([-30.3,-29.9], 0, 2.4, color='r', linestyles='dashed')
plt.vlines([-31.4,-31.2], 0, 2.4, color='orange', linestyles='dashed')
plt.vlines([-33.6,-33.4], 0, 2.4, color='g', linestyles='dashed')
plt.vlines([-37.0,-36.0], 0, 2.4, color='b', linestyles='dashed')
plt.title('EDC and integration regions for XPD')
plt.show()
## XPD plots
fig,ax = plt.subplots(2,2,figsize=(6,4.7), layout='constrained')
res_kx_ky.sel(energy=slice(-30.3,-29.9)).mean('energy').plot(robust=True, ax=ax[0,0], cmap='terrain')
ax[0,0].set_title("XPD of $1^{st}$ order sidebands")
res_kx_ky.sel(energy=slice(-31.4,-31.2)).mean('energy').plot(robust=True, ax=ax[0,1], cmap='terrain')
ax[0,1].set_title("XPD of W4f 7/2 peak")
res_kx_ky.sel(energy=slice(-33.6,-33.4)).mean('energy').plot(robust=True, ax=ax[1,0], cmap='terrain')
ax[1,0].set_title("XPD of W4f 5/2 peak")
res_kx_ky.sel(energy=slice(-37.0,-36.0)).mean('energy').plot(robust=True, ax=ax[1,1], cmap='terrain')
ax[1,1].set_title("XPD of W5p 3/2 peak")
[14]:
Text(0.5, 1.0, 'XPD of W5p 3/2 peak')
As we can see there is some structure visible, but it looks very similar to each other. We probably have to do some normalization to remove the detector structure/artefacts. 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. Unfortunately, we don’t have such a flat-field image.
In this case, we can make a flat-field image from the actual dataset using several different approaches.
As a first option, we can integrate in energy over the whole region and use this image as a background. Additionally, we introduce the Gaussian Blur for comparison.
[15]:
## Background image
bgd = res_kx_ky.mean(('energy'))
## Apply Gaussian Blur to background image
bgd_blur = xr.apply_ufunc(gaussian_filter, bgd, 15)
fig,ax = plt.subplots(1,2,figsize=(6,2.7), layout='constrained')
bgd.plot(robust=True, cmap='terrain', ax=ax[0])
ax[0].set_title('Background image')
bgd_blur.plot(cmap='terrain', ax=ax[1])
ax[1].set_title('Gaussian Blur of background image')
plt.show()
[16]:
## XPD normalized by background image
fig,ax = plt.subplots(2,2,figsize=(6,4.7), layout='constrained')
(res_kx_ky/bgd).sel(energy=slice(-30.3,-29.9)).mean('energy').plot(robust=True, ax=ax[0,0], cmap='terrain')
(res_kx_ky/bgd).sel(energy=slice(-31.4,-31.2)).mean('energy').plot(robust=True, ax=ax[0,1], cmap='terrain')
(res_kx_ky/bgd).sel(energy=slice(-33.6,-33.4)).mean('energy').plot(robust=True, ax=ax[1,0], cmap='terrain')
(res_kx_ky/bgd).sel(energy=slice(-37.0,-36.0)).mean('energy').plot(robust=True, ax=ax[1,1], cmap='terrain')
fig.suptitle(f'Run {run_number}: XPD patterns after background normalization',fontsize='11')
## XPD normalized by Gaussian-blurred background image
fig,ax = plt.subplots(2,2,figsize=(6,4.7), layout='constrained')
(res_kx_ky/bgd_blur).sel(energy=slice(-30.3,-29.9)).mean('energy').plot(robust=True, ax=ax[0,0], cmap='terrain')
(res_kx_ky/bgd_blur).sel(energy=slice(-31.4,-31.2)).mean('energy').plot(robust=True, ax=ax[0,1], cmap='terrain')
(res_kx_ky/bgd_blur).sel(energy=slice(-33.6,-33.4)).mean('energy').plot(robust=True, ax=ax[1,0], cmap='terrain')
(res_kx_ky/bgd_blur).sel(energy=slice(-37.0,-36.0)).mean('energy').plot(robust=True, ax=ax[1,1], cmap='terrain')
fig.suptitle(f'Run {run_number}: XPD patterns after Gaussian-blurred background normalization',fontsize='11')
## XPD normalized by Gaussian-blurred background image and blurred to improve contrast
fig,ax = plt.subplots(2,2,figsize=(6,4.7), layout='constrained')
(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')
(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')
(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')
(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')
fig.suptitle(f'Run {run_number}: resulting Gaussian-blurred XPD patterns',fontsize='11')
[16]:
Text(0.5, 0.98, 'Run 44498: resulting Gaussian-blurred XPD patterns')
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.
[17]:
## XPD normalized by Gaussian-blurred background image
### Define integration regions for XPD
SB = res_kx_ky.sel(energy=slice(-30.3,-29.9)).mean('energy')
W_4f_7 = res_kx_ky.sel(energy=slice(-31.4,-31.2)).mean('energy')
W_4f_5 = res_kx_ky.sel(energy=slice(-33.6,-33.4)).mean('energy')
W_5p = res_kx_ky.sel(energy=slice(-37.0,-36.0)).mean('energy')
### Make corresponding Gaussian Blur background
SB_blur = xr.apply_ufunc(gaussian_filter, SB, 15)
W_4f_7_blur = xr.apply_ufunc(gaussian_filter, W_4f_7, 15)
W_4f_5_blur = xr.apply_ufunc(gaussian_filter, W_4f_5, 15)
W_5p_blur = xr.apply_ufunc(gaussian_filter, W_5p, 15)
### Visualize results
fig,ax = plt.subplots(2,2,figsize=(6,4.7), layout='constrained')
(SB/SB_blur).plot(robust=True, ax=ax[0,0], cmap='terrain')
(W_4f_7/W_4f_7_blur).plot(robust=True, ax=ax[0,1], cmap='terrain')
(W_4f_5/W_4f_5_blur).plot(robust=True, ax=ax[1,0], cmap='terrain')
(W_5p/W_5p_blur).plot(robust=True, ax=ax[1,1], cmap='terrain')
fig.suptitle(f'Run {run_number}: XPD patterns after Gaussian Blur normalization',fontsize='11')
### Apply Gaussian Blur to resulted images to improve contrast
SB_norm = xr.apply_ufunc(gaussian_filter, SB/SB_blur, 1)
W_4f_7_norm = xr.apply_ufunc(gaussian_filter, W_4f_7/W_4f_7_blur, 1)
W_4f_5_norm = xr.apply_ufunc(gaussian_filter, W_4f_5/W_4f_5_blur, 1)
W_5p_norm = xr.apply_ufunc(gaussian_filter, W_5p/W_5p_blur, 1)
### Visualize results
fig,ax = plt.subplots(2,2,figsize=(6,4.7), layout='constrained')
SB_norm.plot(robust=True, ax=ax[0,0], cmap='terrain')
W_4f_7_norm.plot(robust=True, ax=ax[0,1], cmap='terrain')
W_4f_5_norm.plot(robust=True, ax=ax[1,0], cmap='terrain')
W_5p_norm.plot(robust=True, ax=ax[1,1], cmap='terrain')
fig.suptitle(f'Run {run_number}: XPD patterns after Gauss Blur normalization',fontsize='11')
[17]:
Text(0.5, 0.98, 'Run 44498: XPD patterns after Gauss Blur normalization')
Third option for background normalization is to use the simultaneously acquired pre-core level region. 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
[18]:
### Define peak and background region on the high energy side of the peak
W_4f_7 = res_kx_ky.sel(energy=slice(-31.4,-31.2)).mean('energy')
W_4f_7_bgd = res_kx_ky.sel(energy=slice(-32.0,-31.8)).mean('energy')
### Make normalization by background, add Gaussian Blur to the resulting image
W_4f_7_nrm1 = W_4f_7/(W_4f_7_bgd+W_4f_7_bgd.max()*0.00001)
W_4f_7_nrm1_blur = xr.apply_ufunc(gaussian_filter, W_4f_7_nrm1, 1)
### Add Gaussian Blur to the background image, normalize by it and add Gaussian Blur to the resulting image
W_4f_7_bgd_blur = xr.apply_ufunc(gaussian_filter, W_4f_7_bgd, 15)
W_4f_7_nrm2 = W_4f_7/W_4f_7_bgd_blur
W_4f_7_nrm2_blur = xr.apply_ufunc(gaussian_filter, W_4f_7_nrm2, 1)
### Visualize all steps
fig,ax = plt.subplots(4,2,figsize=(6,8), layout='constrained')
W_4f_7.plot(robust=True, ax=ax[0,0], cmap='terrain')
W_4f_7_bgd.plot(robust=True, ax=ax[0,1], cmap='terrain')
W_4f_7_nrm1.plot(robust=True, ax=ax[1,0], cmap='terrain')
W_4f_7_nrm1_blur.plot(robust=True, ax=ax[1,1], cmap='terrain')
W_4f_7_bgd_blur.plot(robust=True, ax=ax[2,0], cmap='terrain')
W_4f_7_nrm2.plot(robust=True, ax=ax[2,1], cmap='terrain')
W_4f_7_nrm2_blur.plot(robust=True, ax=ax[3,0], cmap='terrain')
fig.suptitle(f'Run {run_number}: XPD patterns of W4f7/2 with pre-core level normalization',fontsize='11')
[18]:
Text(0.5, 0.98, 'Run 44498: XPD patterns of W4f7/2 with pre-core level normalization')
[19]:
fig,ax = plt.subplots(1,3,figsize=(6,2), layout='constrained')
(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')
W_4f_7_norm.plot(robust=True, ax=ax[1], cmap='terrain')
W_4f_7_nrm2_blur.plot(robust=True, ax=ax[2], cmap='terrain')
fig.suptitle(f'Run {run_number}: comparison of different normalizations\nof XPD pattern for W4f 7/2 peak with Gaussian Blur',fontsize='11')
[19]:
Text(0.5, 0.98, 'Run 44498: comparison of different normalizations\nof XPD pattern for W4f 7/2 peak with Gaussian Blur')
[ ]: