Tutorial for trXPD for the HEXTOF instrument at FLASH with background normalization

Preparation

Import necessary libraries

%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

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.

# pick the default configuration file for hextof@FLASH
config_file = Path('../sed/config/flash_example_config.yaml')
assert config_file.exists()

# here we setup a dictionary that will be used to override the path configuration
config_override = {
    "core": {
        "beamtime_id": 11019101,
        "paths": {
            "data_raw_dir": path,
            "data_parquet_dir": 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.

en_cal_config = {
    'energy': {
        'calibration': {
            'E0': -54.971004271795664,
            'creation_date': 1718801358.232129,
            'd': 0.8096677238144319,
            'energy_scale': 'kinetic',
            't0': 4.0148196706891397e-07,
            'calib_type': 'fit',
            'fit_function': '(a0/(x0-a1))**2 + a2',
            'coefficients': ([ 8.09667724e-01,  4.01481967e-07, -5.49710043e+01]),
            'axis': 0},
        'tof': None,
        'offsets': {
            'constant': -76.5,
            'creation_date': 1718801360.817963,
            'monochromatorPhotonEnergy': {'preserve_mean': True,'reduction': None,'weight': -1},
            'sampleBias': {'preserve_mean': False, 'reduction': None, 'weight': 1},
            'tofVoltage': {'preserve_mean': True, 'reduction': None, 'weight': -1}}}}

Read data

Now we can use those parameters and load our trXPD data using additional config file

run_number = 44498
sp_44498 = SedProcessor(runs=[run_number], folder_config=en_cal_config, config=config_override, system_config=config_file, verbose=True)
sp_44498.add_jitter()

System config loaded from: [/home/runner/work/sed/sed/docs/sed/config/flash_example_config.yaml]
Default config loaded from: [/home/runner/work/sed/sed/sed/config/default.yaml]
Reading files: 0 new files of 14 total.
All files converted successfully!
Filling nan values...
loading complete in  0.10 s

We can inspect dataframe right after data readout

sp_44498.dataframe.head()

	trainId	pulseId	electronId	timeStamp	dldPosX	dldPosY	dldTimeSteps	cryoTemperature	dldTimeBinSize	extractorCurrent	extractorVoltage	sampleBias	sampleTemperature	tofVoltage	pulserSignAdc	monochromatorPhotonEnergy	gmdBda	bam	delayStage	dldSectorID
0	1628022830	1	0	1.677563e+09	651.227338	895.227338	4595.227338	49.208	0.020576	-0.073857	6029.299805	72.995903	78.989998	19.995356	32919.0	NaN	NaN	-6187.96875	NaN	3
1	1628022830	1	1	1.677563e+09	651.370549	888.370549	4596.370549	49.208	0.020576	-0.073857	6029.299805	72.995903	78.989998	19.995356	32919.0	NaN	NaN	-6187.96875	NaN	0
2	1628022830	5	0	1.677563e+09	681.871577	671.871577	4422.871577	49.208	0.020576	-0.073857	6029.299805	72.995903	78.989998	19.995356	32914.0	NaN	NaN	-6170.15625	NaN	6
3	1628022830	5	1	1.677563e+09	685.191383	658.191383	4425.191383	49.208	0.020576	-0.073857	6029.299805	72.995903	78.989998	19.995356	32914.0	NaN	NaN	-6170.15625	NaN	3
4	1628022830	5	2	1.677563e+09	670.380648	687.380648	4424.380648	49.208	0.020576	-0.073857	6029.299805	72.995903	78.989998	19.995356	32914.0	NaN	NaN	-6170.15625	NaN	5

Now we will do energy calibration, add energy offset, jittering and dld sectors alignment

sp_44498.align_dld_sectors()
sp_44498.append_energy_axis()
sp_44498.add_energy_offset()

Aligning 8s sectors of dataframe
Dask DataFrame Structure:
               trainId pulseId electronId timeStamp  dldPosX  dldPosY dldTimeSteps cryoTemperature dldTimeBinSize extractorCurrent extractorVoltage sampleBias sampleTemperature tofVoltage pulserSignAdc monochromatorPhotonEnergy   gmdBda      bam delayStage dldSectorID
npartitions=14
                uint32   int64      int64   float64  float64  float64      float32         float32        float32          float32          float32    float32           float32    float32       float32                   float32  float32  float32    float32        int8
                   ...     ...        ...       ...      ...      ...          ...             ...            ...              ...              ...        ...               ...        ...           ...                       ...      ...      ...        ...         ...
...                ...     ...        ...       ...      ...      ...          ...             ...            ...              ...              ...        ...               ...        ...           ...                       ...      ...      ...        ...         ...
                   ...     ...        ...       ...      ...      ...          ...             ...            ...              ...              ...        ...               ...        ...           ...                       ...      ...      ...        ...         ...
                   ...     ...        ...       ...      ...      ...          ...             ...            ...              ...              ...        ...               ...        ...           ...                       ...      ...      ...        ...         ...
Dask Name: assign, 9 graph layers
Adding energy column to dataframe:
Using energy calibration parameters generated on 06/19/2024, 12:49:18
Dask DataFrame Structure:
               trainId pulseId electronId timeStamp  dldPosX  dldPosY dldTimeSteps cryoTemperature dldTimeBinSize extractorCurrent extractorVoltage sampleBias sampleTemperature tofVoltage pulserSignAdc monochromatorPhotonEnergy   gmdBda      bam delayStage dldSectorID   energy
npartitions=14
                uint32   int64      int64   float64  float64  float64      float32         float32        float32          float32          float32    float32           float32    float32       float32                   float32  float32  float32    float32        int8  float64
                   ...     ...        ...       ...      ...      ...          ...             ...            ...              ...              ...        ...               ...        ...           ...                       ...      ...      ...        ...         ...      ...
...                ...     ...        ...       ...      ...      ...          ...             ...            ...              ...              ...        ...               ...        ...           ...                       ...      ...      ...        ...         ...      ...
                   ...     ...        ...       ...      ...      ...          ...             ...            ...              ...              ...        ...               ...        ...           ...                       ...      ...      ...        ...         ...      ...
                   ...     ...        ...       ...      ...      ...          ...             ...            ...              ...              ...        ...               ...        ...           ...                       ...      ...      ...        ...         ...      ...
Dask Name: assign, 19 graph layers
Adding energy offset to dataframe:
Using energy offset parameters generated on 06/19/2024, 12:49:20
Energy offset parameters:
   Constant: -76.5
   Column[monochromatorPhotonEnergy]: Weight=-1, Preserve Mean: True,  Reductions: None.
   Column[sampleBias]: Weight=1, Preserve Mean: False,  Reductions: None.
   Column[tofVoltage]: Weight=-1, Preserve Mean: True,  Reductions: None.
Dask DataFrame Structure:
               trainId pulseId electronId timeStamp  dldPosX  dldPosY dldTimeSteps cryoTemperature dldTimeBinSize extractorCurrent extractorVoltage sampleBias sampleTemperature tofVoltage pulserSignAdc monochromatorPhotonEnergy   gmdBda      bam delayStage dldSectorID   energy
npartitions=14
                uint32   int64      int64   float64  float64  float64      float32         float32        float32          float32          float32    float32           float32    float32       float32                   float32  float32  float32    float32        int8  float64
                   ...     ...        ...       ...      ...      ...          ...             ...            ...              ...              ...        ...               ...        ...           ...                       ...      ...      ...        ...         ...      ...
...                ...     ...        ...       ...      ...      ...          ...             ...            ...              ...              ...        ...               ...        ...           ...                       ...      ...      ...        ...         ...      ...
                   ...     ...        ...       ...      ...      ...          ...             ...            ...              ...              ...        ...               ...        ...           ...                       ...      ...      ...        ...         ...      ...
                   ...     ...        ...       ...      ...      ...          ...             ...            ...              ...              ...        ...               ...        ...           ...                       ...      ...      ...        ...         ...      ...
Dask Name: assign, 43 graph layers

sp_44498.attributes.metadata['energy_calibration']

{'applied': True,
 'calibration': {'E0': -54.971004271795664,
  'creation_date': 1718801358.232129,
  'd': 0.8096677238144319,
  'energy_scale': 'kinetic',
  't0': 4.0148196706891397e-07,
  'calib_type': 'fit',
  'fit_function': '(a0/(x0-a1))**2 + a2',
  'coefficients': array([ 8.09667724e-01,  4.01481967e-07, -5.49710043e+01]),
  'axis': 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.

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
)

Adding delay offset to dataframe:
Delay offset parameters:
   Column[bam]: Weight=-0.001, Preserve Mean: True,  Reductions: None.
   Constant: -1448
   Flip delay axis: True
Dask DataFrame Structure:
               trainId pulseId electronId timeStamp  dldPosX  dldPosY dldTimeSteps cryoTemperature dldTimeBinSize extractorCurrent extractorVoltage sampleBias sampleTemperature tofVoltage pulserSignAdc monochromatorPhotonEnergy   gmdBda      bam delayStage dldSectorID   energy
npartitions=14
                uint32   int64      int64   float64  float64  float64      float32         float32        float32          float32          float32    float32           float32    float32       float32                   float32  float32  float32    float64        int8  float64
                   ...     ...        ...       ...      ...      ...          ...             ...            ...              ...              ...        ...               ...        ...           ...                       ...      ...      ...        ...         ...      ...
...                ...     ...        ...       ...      ...      ...          ...             ...            ...              ...              ...        ...               ...        ...           ...                       ...      ...      ...        ...         ...      ...
                   ...     ...        ...       ...      ...      ...          ...             ...            ...              ...              ...        ...               ...        ...           ...                       ...      ...      ...        ...         ...      ...
                   ...     ...        ...       ...      ...      ...          ...             ...            ...              ...              ...        ...               ...        ...           ...                       ...      ...      ...        ...         ...      ...
Dask Name: assign, 62 graph layers

bin in the calibrated energy and corrected delay axis

Visualize trXPS data

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")

Calculate normalization histogram for axis 'delayStage'...

fig,ax = plt.subplots(1,2,figsize=(8,3), 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')

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.

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)

## 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=(9,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")

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.

## 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=(9,4), 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()

## XPD normalized by background image
fig,ax = plt.subplots(2,2,figsize=(9,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='18')

## XPD normalized by Gaussian-blurred background image
fig,ax = plt.subplots(2,2,figsize=(9,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='18')

## XPD normalized by Gaussian-blurred background image and blurred to improve contrast
fig,ax = plt.subplots(2,2,figsize=(9,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='18')

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.

## 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=(9,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='18')

### 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=(9,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='18')

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

### 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=(9,10), 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='18')

Text(0.5, 0.98, 'Run 44498: XPD patterns of W4f7/2 with pre-core level normalization')

fig,ax = plt.subplots(1,3,figsize=(9,3), 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='18')

Text(0.5, 0.98, 'Run 44498: comparison of different normalizations\nof XPD pattern for W4f 7/2 peak with Gaussian Blur')