Processing EIGER images in CrystFEL

This has nothing to do with SACLA but may be useful for serial crystallography at synchrotrons.

Make a pixel mask

Save the following script as make_mask.py.

# 2024-Mar-01: updated for Python 3 and latest h5py
import h5py
import sys

if len(sys.argv) != 3 and len(sys.argv) != 4:
    print("Usage: %s master.h5 mask.h5 [dials_mask.pickle]" % sys.argv)
    sys.exit(-1)

f = h5py.File(sys.argv[1], "r")
mask = f["/entry/instrument/detector/detectorSpecific/pixel_mask"][()]
f.close()

if (len(sys.argv) == 4):
    import numpy as np
    import pickle

    m = pickle.load(open(sys.argv[3], "rb"))[0].as_numpy_array()
    mask[m == False] = 1

# Add more masks (horizontal lines in the middle of ASICs)
# these coordinates are for EIGER 16M

mask[255:261,:] = 1
mask[806:812,:] = 1
mask[1357:1363,:] = 1
mask[1908:1914,:] = 1
mask[2459:2465,:] = 1
mask[3010:3016,:] = 1
mask[3561:3566,:] = 1
mask[4112:4118,:] = 1

out = h5py.File(sys.argv[2], "w")
out["pixel_mask"] = mask
out.close()

Run this script on master.h5 to create a mask h5 file.

python make_mask.py col1_master.h5 mask.h5

If you need more complicated mask (e.g. rings), dials.generate_mask is useful. See Processing EIGER images by DIALS to learn how to use it. You can convert its output to the HDF5 format by

python make_mask.py col1_master.h5 mask.h5 mask.pickle

CrystFEL geometry file

Save the following as eiger.geom. Don't forget to replace photon_energy, adu_per_eV, clen, 0/corner_x and 0/corner_y with appropriate values.

photon_energy   = 12398 ; 12398 / (/entry/instrument/beam/incident_wavelength)
adu_per_eV      =  0.00008065 ; 1 / photon_energy
; adu_per_photon = 1 ; in future versions, you can specify like this
clen            = 0.500 ; /entry/instrument/detector/detector_distance
res             = 13333.3 ; 1 m / 75 um
max_adu         = 245900 ; /entry/instrument/detector/detectorSpecific/countrate_correction_count_cutoff

; used by geoptimiser
rigid_group_0 = 0
rigid_group_collection_0 = 0

mask_good = 0x00
mask_bad  = 0xFF

0/min_fs        = 0
0/max_fs        = 4149 ; /entry/instrument/detector/detectorSpecific/x_pixels_in_detector - 1
0/min_ss        = 0
0/max_ss        = 4370 ; /entry/instrument/detector/detectorSpecific/y_pixels_in_detector - 1
0/corner_x      = -1899.33 ; - /entry/instrument/detector/beam_center_x
0/corner_y      = -2187.37 ; - /entry/instrument/detector/beam_center_y
0/fs            = x
0/ss            = y

0/data = /entry/data/data
0/dim0 = %
0/dim1 = ss
0/dim2 = fs
0/mask = /pixel_mask
0/mask_file = mask.h5

If you have many datasets, the following script generate_eiger_geom.py might be useful:

# 2024-Oct-24: updated for Python 3, reordered items for newer CrystFEL
import h5py
import sys

if len(sys.argv) != 2:
    print("Usage: dials.python generate_eiger_geom.py master.h5 > eiger.geom")
    sys.exit(-1)

f = h5py.File(sys.argv[1], "r")

metadata = {}
metadata["photon_energy"] = 12398.0 / f["/entry/instrument/beam/incident_wavelength"][()]
metadata["adu_per_eV"] = 1 / metadata["photon_energy"]
metadata["clen"] = f["/entry/instrument/detector/detector_distance"][()]
metadata["max_adu"] = f["/entry/instrument/detector/detectorSpecific/countrate_correction_count_cutoff"][()]
metadata["det_x"] = int(f["/entry/instrument/detector/detectorSpecific/x_pixels_in_detector"][()]) - 1
metadata["det_y"] = int(f["/entry/instrument/detector/detectorSpecific/y_pixels_in_detector"][()]) - 1
metadata["center_x"] = - f["/entry/instrument/detector/beam_center_x"][()]
metadata["center_y"] = - f["/entry/instrument/detector/beam_center_y"][()]

f.close()


template = """\
photon_energy   = {photon_energy:.1f} ; 12398 / (/entry/instrument/beam/incident_wavelength)
adu_per_eV      = {adu_per_eV:.8f} ; 1 / photon_energy
; adu_per_photon = 1 ; in future versions, you can specify like this
clen            = {clen:.4f}; /entry/instrument/detector/detector_distance
res             = 13333.3 ; 1 m / 75 um
max_adu         = {max_adu:.1f} ; /entry/instrument/detector/detectorSpecific/countrate_correction_count_cutoff

; used by geoptimiser
rigid_group_0 = 0
rigid_group_collection_0 = 0

mask_good = 0x00
mask_bad  = 0xFF

0/min_fs        = 0
0/max_fs        = {det_x:d} ; /entry/instrument/detector/detectorSpecific/x_pixels_in_detector - 1
0/min_ss        = 0
0/max_ss        = {det_y:d} ; /entry/instrument/detector/detectorSpecific/y_pixels_in_detector - 1
0/corner_x      = {center_x:.2f} ; - /entry/instrument/detector/beam_center_x
0/corner_y      = {center_y:.2f} ; - /entry/instrument/detector/beam_center_y
0/fs            = x
0/ss            = y

0/data = /entry/data/data
0/dim0 = %
0/dim1 = ss
0/dim2 = fs
0/mask = /pixel_mask
0/mask_file = mask.h5
"""

print(template.format(**metadata))

CrystFEL input file list

Save the input file list as test.lst. Note that you have to specify data_NNNNNN.h5 files, not the master.h5 file.

col_1_data_000001.h5
col_1_data_000002.h5
col_1_data_000003.h5

Index and Integrate

Make sure that HDF5_PLUGIN_PATH points to the bitshuffle plugin. The plugin must be linked to the same version of HDF5 library used for CrystFEL. The easiest way is to install DIALS, build and run CrystFEL using HDF5 library that came with DIALS.

export LD_LIBRARY_PATH=$HOME/prog/dials/base/lib:$LD_LIBRARY_PATH
export HDF5_PLUGIN_PATH=$HOME/prog/dials/base/lib

indexamajig -g eiger.geom -i test.lst -o test.stream --indexing=mosflm \
 --int-radius=3,4,7 --peaks=zaef --min-gradient=30 --threshold=10 --min-snr=5

Best parameters depend on the dataset, of course. See Running CrystFEL to learn how to optimize parameters. Since the background was cleaner and spot intensites were weaker in this fine-sliced EIGER dataset than in MPCCD SFX images, the threshold and min-gradient are smaller here.

Indexing only hits

In addition to file names, you can specify image numbers in the input list.

col_1_data_000001.h5 //0
col_1_data_000001.h5 //30
col_1_data_000001.h5 //99
col_1_data_000002.h5 //0
col_1_data_000002.h5 //45

Note that image numbers start from 0 in CrystFEL. Since we have 100 frames in a file, this input instructs CrystFEL to process images 1, 31, 100, 101 and 146 (counting from 1).

We can write a script to make the list from Cheetah's log file frames.txt.

Why cannot we use master.h5 instead?

In principle, setting

0/data = /entry/data/%

in the geometry file should work. Unfortunately, CrystFEL complains:

All panels' data and mask entries must have the same number of placeholders.
Error: Data path ends with a / symbol

This is fixable, but we took the simplest approach in this tutorial.