From e1cea7059ba00610f8eac196df9dab6870f27645 Mon Sep 17 00:00:00 2001 From: Aaron Brewster Date: Tue, 25 Aug 2020 12:54:13 -0700 Subject: [PATCH] Support Epix10k detector at MFX endstation at LCLS https://lcls.slac.stanford.edu/detectors/ePix10k Nearly a copy of FormatXTCJungfrau, but taking account the different shape of the data and the deeper metrology with quadrants --- format/FormatXTCEpix.py | 195 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 195 insertions(+) create mode 100644 format/FormatXTCEpix.py diff --git a/format/FormatXTCEpix.py b/format/FormatXTCEpix.py new file mode 100644 index 000000000..c1cd2549c --- /dev/null +++ b/format/FormatXTCEpix.py @@ -0,0 +1,195 @@ +from __future__ import absolute_import, division, print_function + +import sys +from builtins import range + +import numpy as np +import psana + +from libtbx.phil import parse +from scitbx.array_family import flex +from scitbx.matrix import col + +from dxtbx.format.FormatXTC import FormatXTC, locator_str +from dxtbx.model import Detector + +try: + from xfel.cxi.cspad_ana import cspad_tbx +except ImportError: + # xfel not configured + pass + +epix_locator_str = """ +""" + +epix_locator_scope = parse(epix_locator_str + locator_str, process_includes=True) + + +class FormatXTCEpix(FormatXTC): + def __init__(self, image_file, **kwargs): + super(FormatXTCEpix, self).__init__( + image_file, locator_scope=epix_locator_scope, **kwargs + ) + self._ds = FormatXTC._get_datasource(image_file, self.params) + self._env = self._ds.env() + self.populate_events() + self.n_images = len(self.times) + self._cached_detector = {} + self._cached_psana_detectors = {} + + @staticmethod + def understand(image_file): + try: + params = FormatXTC.params_from_phil(epix_locator_scope, image_file) + except Exception: + return False + return any(["epix" in src.lower() for src in params.detector_address]) + + def get_raw_data(self, index): + d = FormatXTCEpix.get_detector(self, index) + evt = self._get_event(index) + run = self.get_run_from_index(index) + if run.run() not in self._cached_psana_detectors: + assert len(self.params.detector_address) == 1 + self._cached_psana_detectors[run.run()] = psana.Detector( + self.params.detector_address[0], self._env + ) + det = self._cached_psana_detectors[run.run()] + data = det.calib(evt) + data = data.astype(np.float64) + # the shape of the epix 10k data is (16, 352, 384) + + self._raw_data = [] + for quad_count, quad in enumerate(d.hierarchy()): + for module_count, module in enumerate(quad): + for asic_count, asic in enumerate(module): + fdim, sdim = asic.get_image_size() + twobytwo_id = (4 * quad_count) + module_count + sensor_id = asic_count // 2 # There are 2X2 asics per module + asic_in_sensor_id = asic_count % 2 # this number will be 0 or 1 + asic_data = data[twobytwo_id][ + sensor_id * sdim : (sensor_id + 1) * sdim, + asic_in_sensor_id * fdim : (asic_in_sensor_id + 1) * fdim, + ] # 8 sensors per module + self._raw_data.append(flex.double(np.array(asic_data))) + assert len(d) == len(self._raw_data), (len(d), len(self._raw_data)) + return tuple(self._raw_data) + + def get_num_images(self): + return self.n_images + + def get_detector(self, index=None): + return FormatXTCEpix._detector(self, index) + + def get_beam(self, index=None): + return self._beam(index) + + def _beam(self, index=None): + """Returns a simple model for the beam""" + if index is None: + index = 0 + evt = self._get_event(index) + wavelength = cspad_tbx.evt_wavelength(evt) + if wavelength is None: + return None + return self._beam_factory.simple(wavelength) + + def get_goniometer(self, index=None): + return None + + def get_scan(self, index=None): + return None + + def _detector(self, index=None): + from xfel.cftbx.detector.cspad_cbf_tbx import basis_from_geo + from PSCalib.SegGeometryStore import sgs + + run = self.get_run_from_index(index) + if run.run() in self._cached_detector: + return self._cached_detector[run.run()] + + if index is None: + index = 0 + self._env = self._ds.env() + assert len(self.params.detector_address) == 1 + self._det = psana.Detector(self.params.detector_address[0], self._env) + geom = self._det.pyda.geoaccess(self._get_event(index).run()) + pixel_size = ( + self._det.pixel_size(self._get_event(index)) / 1000.0 + ) # convert to mm + d = Detector() + pg0 = d.hierarchy() + # first deal with D0 + det_num = 0 + root = geom.get_top_geo() + root_basis = basis_from_geo(root) + while len(root.get_list_of_children()) == 1: + sub = root.get_list_of_children()[0] + sub_basis = basis_from_geo(sub) + root = sub + root_basis = root_basis * sub_basis + t = root_basis.translation + root_basis.translation = col((t[0], t[1], -t[2])) + + origin = col((root_basis * col((0, 0, 0, 1)))[0:3]) + fast = col((root_basis * col((1, 0, 0, 1)))[0:3]) - origin + slow = col((root_basis * col((0, 1, 0, 1)))[0:3]) - origin + pg0.set_local_frame(fast.elems, slow.elems, origin.elems) + pg0.set_name("D%d" % (det_num)) + + # Now deal with quads + for quad_num in range(len(root.get_list_of_children())): + pg1 = pg0.add_group() + quad = root.get_list_of_children()[quad_num] + quad_basis = basis_from_geo(quad) + origin = col((quad_basis * col((0, 0, 0, 1)))[0:3]) + fast = col((quad_basis * col((1, 0, 0, 1)))[0:3]) - origin + slow = col((quad_basis * col((0, 1, 0, 1)))[0:3]) - origin + pg1.set_local_frame(fast.elems, slow.elems, origin.elems) + pg1.set_name("D%dQ%d" % (det_num, quad_num)) + + # Modules next + for module_num in range(len(quad.get_list_of_children())): + pg2 = pg1.add_group() + module = quad.get_list_of_children()[module_num] + module_basis = basis_from_geo(module) + origin = col((module_basis * col((0, 0, 0, 1)))[0:3]) + fast = col((module_basis * col((1, 0, 0, 1)))[0:3]) - origin + slow = col((module_basis * col((0, 1, 0, 1)))[0:3]) - origin + pg2.set_local_frame(fast.elems, slow.elems, origin.elems) + pg2.set_name("D%dQ%dM%d" % (det_num, quad_num, module_num)) + + # Read the known layout of the Epix 2x2 module + sg = sgs.Create(segname=module.oname) + xx, yy = sg.get_seg_xy_maps_um() + xx = xx / 1000 + yy = yy / 1000 + + # Now deal with ASICs + for asic_num in range(4): + val = "ARRAY_D0Q%dM%dA%d" % (module_num, quad_num, asic_num) + p = pg2.add_panel() + dim_slow = xx.shape[0] + dim_fast = xx.shape[1] + sensor_id = asic_num // 2 # There are 2X2 asics per module + asic_in_sensor_id = asic_num % 2 # this number will be 0 or 1 + id_slow = sensor_id * (dim_slow // 2) + id_fast = asic_in_sensor_id * (dim_fast // 2) + origin = col((xx[id_slow][id_fast], yy[id_slow][id_fast], 0,)) + fp = col((xx[id_slow][id_fast + 1], yy[id_slow][id_fast + 1], 0)) + sp = col((xx[id_slow + 1][id_fast], yy[id_slow + 1][id_fast], 0)) + fast = (fp - origin).normalize() + slow = (sp - origin).normalize() + p.set_local_frame(fast.elems, slow.elems, origin.elems) + p.set_pixel_size((pixel_size, pixel_size)) + p.set_image_size((dim_fast // 2, dim_slow // 2)) + p.set_trusted_range((-1, 2e6)) + p.set_name(val) + self._cached_detector[run.run()] = d + return d + + +if __name__ == "__main__": + for arg in sys.argv[1:]: + # Bug, should call this part differently for understand method to work + print(FormatXTCEpix.understand(arg))