Files
crystfel_tools/geom_files/jf8m_cxigeom2nexus.py
2024-02-13 16:46:10 +01:00

495 lines
25 KiB
Python

from __future__ import absolute_import, division, print_function
from six.moves import range
import os
import h5py
import numpy as np
from read_geom import read_geom
from libtbx.phil import parse
import six
from libtbx.utils import Sorry
import datetime
from jungfrau import pad_stacked_format
phil_scope = parse("""
unassembled_file = None
.type = path
.help = hdf5 file used to read in image data.
geom_file = None
.type = path
.help = geometry file to be read in for detector (.geom).
output_file = None
.type = path
.help = output file path
detector_distance = None
.type = float
.help = Detector distance
wavelength = None
.type = float
.help = If not provided, try to find wavelength in unassembled file.
beam_file = None
.type = path
.help = Overrides wavelength. Reads the pulse IDs in the provided file \
to get a list of wavelengths for the master.
include_spectra = False
.type = bool
.help = If true, 2D spectral data will be included in the master file, \
as read from the beam_file.
energy_offset = None
.type = float
.help = If set, add this offset (in eV) to the energy axis in the \
spectra in the beam file and to the per-shot wavelength.
mask_file = None
.type = str
.help = Path to file with external bad pixel mask.
split_modules_into_asics = True
.type = bool
.help = Whether to split the 4x2 modules into indivdual asics \
accounting for borders and gaps.
trusted_range = None
.type = floats(size=2)
.help = Set the trusted range
raw = False
.type = bool
.help = Whether the data being analyzed is raw data from the JF16M or has \
been corrected and padded.
unassembled_data_key = None
.type = str
.expert_level = 2
.help = Override hdf5 key name in unassembled file
spectral_data_key = None
.type = str
.expert_level = 2
.help = Override beam hdf5 file key for spectral data
pedestal_file = None
.type = str
.help = path to Jungfrau pedestal file
gain_file = None
.type = str
.help = path to Jungfrau gain file
raw_file = None
.type = str
.help = path to Jungfrau raw file
recalculate_wavelength = False
.type = bool
.help = If True, the incident_wavelength is recalculated from the \
spectrum. The original wavelength is saved as \
incident_wavelength_fitted_spectrum_average.
nexus_details {
instrument_name = SwissFEL ARAMIS BEAMLINE ESC
.type = str
.help = Name of instrument
instrument_short_name = Cristallina
.type = str
.help = short name for instrument, perhaps the acronym
source_name = SwissFEL ARAMIS
.type = str
.help = Name of the neutron or x-ray storage ring/facility
source_short_name = SwissFEL ARAMIS
.type = str
.help = short name for source, perhaps the acronym
start_time = None
.type = str
.help = ISO 8601 time/date of the first data point collected in UTC, \
using the Z suffix to avoid confusion with local time
end_time = None
.type = str
.help = ISO 8601 time/date of the last data point collected in UTC, \
using the Z suffix to avoid confusion with local time. \
This field should only be filled when the value is accurately \
observed. If the data collection aborts or otherwise prevents \
accurate recording of the end_time, this field should be omitted
end_time_estimated = None
.type = str
.help = ISO 8601 time/date of the last data point collected in UTC, \
using the Z suffix to avoid confusion with local time. \
This field may be filled with a value estimated before an \
observed value is avilable.
sample_name = None
.type = str
.help = Descriptive name of sample
total_flux = None
.type = float
.help = flux incident on beam plane in photons per second
}
""")
'''
This script creates a master nexus file by taking in as input a) an hdf5 file and b) a .geom file
The hd5f file is generated by the JF16M after processing the raw images and doing appropriate gain corrections
The assumed parameters for the detector can be seen in the __init__ function and should be changed
if they are modified at in the future
'''
class jf16m_cxigeom2nexus(object):
def __init__(self, args):
self.params_from_phil(args)
if self.params.detector_distance == None:
self.params.detector_distance = 100.0 # Set detector distance arbitrarily if nothing is provided
self.hierarchy = read_geom(self.params.geom_file)
self.n_quads = 4
self.n_modules = 4
def params_from_phil(self, args):
user_phil = []
for arg in args:
if os.path.isfile(arg):
user_phil.append(parse(file_name=arg))
else:
try:
user_phil.append(parse(arg))
except Exception as e:
raise Sorry("Unrecognized argument: %s"%arg)
self.params = phil_scope.fetch(sources=user_phil).extract()
def _create_scalar(self, handle,path,dtype,value):
dataset = handle.create_dataset(path, (),dtype=dtype)
dataset[()] = value
def create_vector(self,handle, name, value, **attributes):
handle.create_dataset(name, (1,), data = [value], dtype='f')
for key,attribute in six.iteritems(attributes):
handle[name].attrs[key] = attribute
def create_nexus_master_file(self):
'''
Hierarchical structure of master nexus file. Format information available here
http://download.nexusformat.org/sphinx/classes/base_classes/NXdetector_module.html#nxdetector-module
--> entry
--> data
--> definition (leaf)
--> instrument
--> sample
'''
output_file_name = self.params.output_file if self.params.output_file is not None else os.path.splitext(self.params.unassembled_file)[0]+'_master.h5'
f = h5py.File(output_file_name, 'w')
f.attrs['NX_class'] = 'NXroot'
f.attrs['file_name'] = os.path.basename(output_file_name)
f.attrs['file_time'] = datetime.datetime.utcnow().strftime("%Y-%m-%dT%H:%M:%SZ")
f.attrs['HDF5_Version'] = h5py.version.hdf5_version
entry = f.create_group('entry')
entry.attrs['NX_class'] = 'NXentry'
if self.params.nexus_details.start_time: entry['start_time'] = self.params.nexus_details.start_time
if self.params.nexus_details.end_time: entry['end_time'] = self.params.nexus_details.end_time
if self.params.nexus_details.end_time_estimated: entry['end_time_estimated'] = self.params.nexus_details.end_time_estimated
# --> definition
self._create_scalar(entry, 'definition', 'S4', np.string_('NXmx'))
# --> data
data = entry.create_group('data')
data.attrs['NX_class'] = 'NXdata'
data_key = 'data'
if self.params.unassembled_data_key:
unassembled_data_key = self.params.unassembled_data_key
else:
if self.params.raw:
unassembled_data_key = "data/JF17T16V01"
else:
unassembled_file = h5py.File(self.params.unassembled_file, "r")
unassembled_data_key = "data"
if not unassembled_file.get(unassembled_data_key):
unassembled_data_key = "data/JF17T16V01"
if not unassembled_file.get(unassembled_data_key):
raise Sorry("couldn't find unassembled data in the indicated file")
unassembled_file.close()
data[data_key] = h5py.ExternalLink(self.params.unassembled_file, unassembled_data_key+"/data")
if self.params.raw_file is not None:
assert not self.params.raw
with h5py.File(self.params.pedestal_file, "r") as pedh5:
print("Padding raw pedestal data")
mean_pedestal = [pad_stacked_format(raw) for raw in pedh5["gains"]]
print("Padding raw pedestal RMS data")
sigma_pedestal = [pad_stacked_format(raw) for raw in pedh5["gainsRMS"]]
data.create_dataset("pedestal", data=mean_pedestal, dtype=np.float32)
data.create_dataset('pedestalRMS', data=sigma_pedestal, dtype=np.float32)
with h5py.File(self.params.gain_file, "r") as gainh5:
print("Padding gains")
gains = [pad_stacked_format(raw) for raw in gainh5["gains"]]
data.create_dataset("gains", data=gains, dtype=np.float32)
data.attrs['signal'] = 'data'
raw_file_handle = h5py.File(self.params.raw_file, "r")
res_file_handle = h5py.File(self.params.unassembled_file, "r")
raw_dset = raw_file_handle["data/JF17T16V01/data"]
raw_shape = raw_dset.shape
_, raw_slowDim, raw_fastDim = raw_shape
raw_type = raw_dset.dtype
num_imgs = res_file_handle['data/data'].shape[0]
raw_layout = h5py.VirtualLayout(shape=(num_imgs, raw_slowDim, raw_fastDim), dtype=raw_type)
raw_pulses = raw_file_handle['data/JF17T16V01/pulse_id'][()][:, 0]
assert np.all(raw_pulses == np.sort(raw_pulses)) # NOTE; this is quick, however I think this is always the case
res_pulses = h5py.File(self.params.unassembled_file, 'r')['data/pulse_id'][()]
raw_source = h5py.VirtualSource(self.params.raw_file, 'data/JF17T16V01/data', shape=raw_shape)
for res_imgnum, raw_imgnum in enumerate(np.searchsorted(raw_pulses, res_pulses)):
raw_layout[res_imgnum] = raw_source[raw_imgnum]
data.create_virtual_dataset('raw', raw_layout)
if self.params.raw:
if self.params.pedestal_file:
# named gains instead of pedestal in JF data files
data['pedestal'] = h5py.ExternalLink(self.params.pedestal_file, 'gains')
data['pedestalRMS'] = h5py.ExternalLink(self.params.pedestal_file, 'gainsRMS')
if self.params.gain_file:
data['gains'] = h5py.ExternalLink(self.params.gain_file, 'gains')
if self.params.pedestal_file or self.params.gain_file:
data.attrs['signal'] = 'data'
#--> sample
sample = entry.create_group('sample')
sample.attrs['NX_class'] = 'NXsample'
if self.params.nexus_details.sample_name: sample['name'] = self.params.nexus_details.sample_name
sample['depends_on'] = '.' # This script does not support scans/gonios
# --> source
source = entry.create_group('source')
source.attrs['NX_class'] = 'NXsource'
source['name'] = self.params.nexus_details.source_name
source['name'].attrs['short_name'] = self.params.nexus_details.source_short_name
# --> instrument
instrument = entry.create_group('instrument')
instrument.attrs['NX_class'] = 'NXinstrument'
instrument["name"] = self.params.nexus_details.instrument_name
instrument["name"].attrs["short_name"] = self.params.nexus_details.instrument_short_name
beam = instrument.create_group('beam')
beam.attrs['NX_class'] = 'NXbeam'
if self.params.nexus_details.total_flux:
self._create_scalar(beam, 'total_flux', 'f', self.params.nexus_details.total_flux)
beam['total_flux'].attrs['units'] = 'Hz'
if self.params.wavelength is None and self.params.beam_file is None:
wavelengths = h5py.File(self.params.unassembled_file, 'r')['instrument/photon_wavelength_A']
beam.create_dataset('incident_wavelength', (1,), data=np.mean(wavelengths), dtype='f8')
elif self.params.beam_file is not None:
# data file has pulse ids, need to match those to the beam file, which may have more pulses
if self.params.raw:
data_pulse_ids = h5py.File(self.params.unassembled_file, 'r')[unassembled_data_key+'/pulse_id'][()]
else:
data_pulse_ids = h5py.File(self.params.unassembled_file, 'r')[unassembled_data_key+'/pulse_id'][()]
beam_h5 = h5py.File(self.params.beam_file, 'r')
if self.params.spectral_data_key:
spectral_data_key = self.params.spectral_data_key
else:
spectral_data_key = "SARFE10-PSSS059"
if not beam_h5.get(spectral_data_key+":SPECTRUM_CENTER/pulse_id"):
spectral_data_key = "data/SARFE10-PSSS059"
if not beam_h5.get(spectral_data_key+":SPECTRUM_CENTER/pulse_id"):
raise Sorry("couldn't find spectral data at the indicated address")
beam_pulse_ids = beam_h5[spectral_data_key+':SPECTRUM_CENTER/pulse_id'][()]
beam_energies = beam_h5[spectral_data_key+':SPECTRUM_CENTER/data'][()]
energies = np.ndarray((len(data_pulse_ids),), dtype='f8')
if self.params.recalculate_wavelength:
orig_energies = np.ndarray((len(data_pulse_ids),), dtype='f8')
if self.params.include_spectra:
beam_spectra_x = beam_h5[spectral_data_key+':SPECTRUM_X/data'][()]
beam_spectra_y = beam_h5[spectral_data_key+':SPECTRUM_Y/data'][()]
spectra_x = np.ndarray((len(data_pulse_ids),beam_spectra_x.shape[1]), dtype='f8')
spectra_y = np.ndarray((len(data_pulse_ids),beam_spectra_y.shape[1]), dtype='f8')
for i, pulse_id in enumerate(data_pulse_ids):
where = np.where(beam_pulse_ids==pulse_id)[0][0]
if self.params.recalculate_wavelength:
assert self.params.beam_file is not None, "Must provide beam file to recalculate wavelength"
x = beam_spectra_x[where]
y = beam_spectra_y[where].astype(float)
ycorr = y - np.percentile(y, 10)
e = sum(x*ycorr) / sum(ycorr)
energies[i] = e
orig_energies[i] = beam_energies[where]
else:
energies[i] = beam_energies[where]
if self.params.include_spectra:
spectra_x[i] = beam_spectra_x[where]
spectra_y[i] = beam_spectra_y[where]
if self.params.energy_offset:
energies += self.params.energy_offset
wavelengths = 12398.4187/energies
if self.params.recalculate_wavelength:
orig_wavelengths = 12398.4187/orig_energies
beam.create_dataset('incident_wavelength', wavelengths.shape, data=wavelengths, dtype=wavelengths.dtype)
if self.params.recalculate_wavelength:
beam.create_dataset(
'incident_wavelength_fitted_spectrum_average',
orig_wavelengths.shape,
data=orig_wavelengths,
dtype=orig_wavelengths.dtype
)
beam['incident_wavelength_fitted_spectrum_average'].attrs['variant'] = 'incident_wavelength'
if self.params.include_spectra:
if self.params.energy_offset:
spectra_x += self.params.energy_offset
if (spectra_x == spectra_x[0]).all(): # check if all rows are the same
spectra_x = spectra_x[0]
beam.create_dataset('incident_wavelength_1Dspectrum', spectra_x.shape, data=12398.4187/spectra_x, dtype=spectra_x.dtype)
beam.create_dataset('incident_wavelength_1Dspectrum_weight', spectra_y.shape, data=spectra_y, dtype=spectra_y.dtype)
beam['incident_wavelength_1Dspectrum'].attrs['units'] = 'angstrom'
beam['incident_wavelength'].attrs['variant'] = 'incident_wavelength_1Dspectrum'
else:
assert self.params.wavelength is not None, "Provide a wavelength"
beam.create_dataset('incident_wavelength', (1,), data=self.params.wavelength, dtype='f8')
beam['incident_wavelength'].attrs['units'] = 'angstrom'
jf16m = instrument.create_group('JF8M')
jf16m.attrs['NX_class'] = 'NXdetector_group'
jf16m.create_dataset('group_index', data = list(range(1,3)), dtype='i')
data = [np.string_('JF8M'),np.string_('ELE_D0')]
jf16m.create_dataset('group_names',(2,), data=data, dtype='S12')
jf16m.create_dataset('group_parent',(2,), data=[-1,1], dtype='i')
jf16m.create_dataset('group_type', (2,), data=[1,2], dtype='i')
detector = instrument.create_group('ELE_D0')
detector.attrs['NX_class'] = 'NXdetector'
detector['description'] = 'JUNGFRAU 8M'
detector['depends_on'] = '/entry/instrument/ELE_D0/transformations/AXIS_RAIL'
detector['gain_setting'] = 'auto'
detector['sensor_material'] = 'Si'
self._create_scalar(detector, 'sensor_thickness', 'f', 320.)
self._create_scalar(detector, 'bit_depth_readout', 'i', 16)
self._create_scalar(detector, 'count_time', 'f', 10)
self._create_scalar(detector, 'frame_time', 'f', 40)
detector['sensor_thickness'].attrs['units'] = 'microns'
detector['count_time'].attrs['units'] = 'us'
detector['frame_time'].attrs['units'] = 'us'
transformations = detector.create_group('transformations')
transformations.attrs['NX_class'] = 'NXtransformations'
# Create AXIS leaves for RAIL, D0 and different hierarchical levels of detector
self.create_vector(transformations, 'AXIS_RAIL', self.params.detector_distance, depends_on='.', equipment='detector', equipment_component='detector_arm',transformation_type='translation', units='mm', vector=(0., 0., 1.), offset=(0.,0.,0.))
self.create_vector(transformations, 'AXIS_D0', 0.0, depends_on='AXIS_RAIL', equipment='detector', equipment_component='detector_arm',transformation_type='rotation', units='degrees', vector=(0., 0., -1.), offset=self.hierarchy.local_origin, offset_units = 'mm')
# Add 4 quadrants
# Nexus coordiate system, into the board JF16M detector
# o --------> (+x) Q3=(12,13,14,15) Q0=(0,1,2,3)
# | o
# | Q2=(8,9,10,11) Q1=(4,5,6,7)
# v
# (+y)
panels = []
for q, quad in six.iteritems(self.hierarchy):
panels.extend([quad[key] for key in quad])
pixel_size = panels[0]['pixel_size']
assert [pixel_size == panels[i+1]['pixel_size'] for i in range(len(panels)-1)].count(False) == 0
if self.params.raw:
fast = 1030; slow = 8224
else:
fast = max([int(panel['max_fs']) for panel in panels])+1
slow = max([int(panel['max_ss']) for panel in panels])+1
quad_fast = fast
quad_slow = slow // self.n_quads
module_fast = quad_fast
module_slow = quad_slow // self.n_modules
if self.params.split_modules_into_asics:
n_fast_asics = 4; n_slow_asics = 2
if self.params.raw:
border = 1; data_gap = 0; real_gap = 2
asic_fast = (module_fast - (n_fast_asics-1)*data_gap) / n_fast_asics # includes border but not gap
asic_slow = (module_slow - (n_slow_asics-1)*data_gap) / n_slow_asics # includes border but not gap
else:
border = 1; gap = 2
asic_fast = (module_fast - (n_fast_asics-1)*gap) / n_fast_asics # includes border but not gap
asic_slow = (module_slow - (n_slow_asics-1)*gap) / n_slow_asics # includes border but not gap
array_name = 'ARRAY_D0'
if self.params.mask_file is not None:
self._create_scalar(detector, 'pixel_mask_applied', 'b', True)
#detector['pixel_mask'] = h5py.ExternalLink(self.params.mask_file, "mask") # If mask was formatted right, could use this
mask = h5py.File(self.params.mask_file, 'r')['pixel_mask'][()].astype(np.int32)
detector.create_dataset('pixel_mask', mask.shape, data=mask==0, dtype=mask.dtype)
if self.params.trusted_range is not None:
underload, overload = self.params.trusted_range
detector.create_dataset('underload_value', (1,), data=[underload], dtype='int32')
detector.create_dataset('saturation_value', (1,), data=[overload], dtype='int32')
alias = 'data'
data_name = 'data'
detector[alias] = h5py.SoftLink('/entry/data/%s'%data_name)
for q_key in sorted(self.hierarchy.keys()):
quad = int(q_key.lstrip('quad'))
q_name = 'AXIS_D0Q%d'%quad
quad_vector = self.hierarchy[q_key].local_origin.elems
self.create_vector(transformations, q_name, 0.0, depends_on='AXIS_D0', equipment='detector', equipment_component='detector_quad',transformation_type='rotation', units='degrees', vector=(0., 0., -1.), offset = quad_vector, offset_units = 'mm')
for m_key in sorted(self.hierarchy[q_key].keys()):
module_num = int(m_key.lstrip('m'))
m_name = 'AXIS_D0Q%dM%d'%(quad, module_num)
module_vector = self.hierarchy[q_key][m_key]['local_origin']
fast = self.hierarchy[q_key][m_key]['local_fast']
slow = self.hierarchy[q_key][m_key]['local_slow']
if self.params.split_modules_into_asics:
# module_vector points to the corner of the module. Instead, point to the center of it.
offset = (module_fast/2 * pixel_size * fast) + (module_slow/2 * pixel_size * slow)
module_vector = module_vector + offset
self.create_vector(transformations, m_name, 0.0, depends_on=q_name, equipment='detector', equipment_component='detector_module',transformation_type='rotation', units='degrees', vector=(0., 0., -1.), offset = module_vector.elems, offset_units = 'mm')
for asic_fast_number in range(n_fast_asics):
for asic_slow_number in range(n_slow_asics):
asic_num = asic_fast_number + (asic_slow_number * n_fast_asics)
a_name = 'AXIS_D0Q%dM%dA%d'%(quad, module_num, asic_num)
# Modules look like this:
# bbbbbggbbbbb Assuming a 3x3 asic (X), a 1 pixel border (b) and a 2 pixel gap (g).
# bXXXbggbXXXb This is as if there were only two asics in a module
# bXXXbggbXXXb The math below skips past the borders and gaps to the first real pixel.
# bXXXbggbXXXb
# bbbbbggbbbbb
if self.params.raw:
asic_vector = -offset + (fast * pixel_size * (asic_fast * asic_fast_number + border + (asic_fast_number * real_gap))) + \
(slow * pixel_size * (asic_slow * asic_slow_number + border + (asic_slow_number * real_gap)))
else:
asic_vector = -offset + (fast * pixel_size * (asic_fast * asic_fast_number + border + (asic_fast_number * gap))) + \
(slow * pixel_size * (asic_slow * asic_slow_number + border + (asic_slow_number * gap)))
self.create_vector(transformations, a_name, 0.0, depends_on=m_name, equipment='detector', equipment_component='detector_asic',
transformation_type='rotation', units='degrees', vector=(0., 0., -1.), offset = asic_vector.elems, offset_units = 'mm')
asicmodule = detector.create_group(array_name+'Q%dM%dA%d'%(quad,module_num,asic_num))
asicmodule.attrs['NX_class'] = 'NXdetector_module'
if self.params.raw:
asicmodule.create_dataset('data_origin', (2,), data=[module_slow * module_num + asic_slow * asic_slow_number + border + data_gap*asic_slow_number,
asic_fast * asic_fast_number + border + data_gap*asic_fast_number], dtype='i')
else:
asicmodule.create_dataset('data_origin', (2,), data=[module_slow * module_num + asic_slow * asic_slow_number + border + gap*asic_slow_number,
asic_fast * asic_fast_number + border + gap*asic_fast_number], dtype='i')
asicmodule.create_dataset('data_size', (2,), data=[asic_slow - border*2, asic_fast - border*2], dtype='i')
self.create_vector(asicmodule, 'fast_pixel_direction',pixel_size,
depends_on=transformations.name+'/AXIS_D0Q%dM%dA%d'%(quad,module_num,asic_num),
transformation_type='translation', units='mm', vector=fast.elems, offset=(0. ,0., 0.))
self.create_vector(asicmodule, 'slow_pixel_direction',pixel_size,
depends_on=transformations.name+'/AXIS_D0Q%dM%dA%d'%(quad,module_num,asic_num),
transformation_type='translation', units='mm', vector=slow.elems, offset=(0., 0., 0.))
else:
module_vector = self.hierarchy[q_key][m_key]['local_origin'].elems
self.create_vector(transformations, m_name, 0.0, depends_on=q_name,
equipment='detector', equipment_component='detector_module',
transformation_type='rotation', units='degrees', vector=(0. ,0., -1.),
offset = module_vector, offset_units = 'mm')
modulemodule = detector.create_group(array_name+'Q%dM%d'%(quad,module_num))
modulemodule.attrs['NX_class'] = 'NXdetector_module'
modulemodule.create_dataset('data_origin', (2,), data=[module_slow * module_num, 0],
dtype='i')
modulemodule.create_dataset('data_size', (2,), data=[module_slow, module_fast], dtype='i')
fast = self.hierarchy[q_key][m_key]['local_fast'].elems
slow = self.hierarchy[q_key][m_key]['local_slow'].elems
self.create_vector(modulemodule, 'fast_pixel_direction',pixel_size,
depends_on=transformations.name+'/AXIS_D0Q%dM%d'%(quad,module_num),
transformation_type='translation', units='mm', vector=fast, offset=(0. ,0., 0.))
self.create_vector(modulemodule, 'slow_pixel_direction',pixel_size,
depends_on=transformations.name+'/AXIS_D0Q%dM%d'%(quad,module_num),
transformation_type='translation', units='mm', vector=slow, offset=(0., 0., 0.))
f.close()
if __name__ == '__main__':
import sys
nexus_helper = jf16m_cxigeom2nexus(sys.argv[1:])
nexus_helper.create_nexus_master_file()