From 57031407c10f59e5a5cd58ee185c299d13c88bae Mon Sep 17 00:00:00 2001 From: "Aaron S. Brewster" Date: Tue, 26 Sep 2023 15:59:55 -0700 Subject: [PATCH] Break the DIALS/xfel circular dependency (#872) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Brief summary: * A new top level project in cctbx_project has been created: serialtbx * Any code needed by DIALS and dxtbx that was in xfel has been moved to serialtbx or directly into DIALS/dxtbx. That mostly includes code to read data at XFELs used by dxtbx, and stills-specific indexing/integration methods. * serialtbx has no imports of DIALS or dxtbx and so it can be directly configured by the conda packages cctbx and cctbx-base, @bkpoon. That said, there are no dispatchers or libtbx_config file so nothing specifically needs to be done to bootstrap.py, I believe. * All imports of xfel in DIALS and dxtbx have been updated to use the new serialtbx or DIALS/dxtbx imports * All other known repos that have xfel imports have been updated to use the new serialtbx or DIALS/dxtbx imports * xfel will remain configured in the DIALS binary bundles on dials.github.io and distributed by phenix. For conda, the plan is to make a cctbx.xfel conda package that includes all its dependencies (MPI, MySQL, etc.) Se also cctbx/dxtbx#627 and dials/dials#2404 Detailed listing of code movements: * xfel.util.jungfrau → serialtbx.detector From xfel.cftbx.detector: * cspad_detector → serialtbx.detector.legacy_metrology * generic_detector → serialtbx.detector.legacy_metrology * metrology → serialtbx.detector.legacy_metrology From xfel.mono_simulation: * max_like → serialtbx.mono_simulation * util → serialtbx.mono_simulation From xfel.cftbx.detector.cspad_cbf_tbx: * Constants moved to serialtbx.detector.cspad * get_psana_corrected_data → serialtbx.detector.cspad * cbf_wrapper → merged with dxtbx.format.FormatCBFMultiTile.cbf_wrapper * angle_and_axis → dxtbx.format.FormatCBFMultiTile.cbf_wrapper * center → serialtbx.detector * basis → serialtbx.detector * basis_from_geo → serialtbx.detector.xtc * read_slac_metrology → serialtbx.detector.cspad * add_frame_specific_cbf_tables → dxtbx.format.cbf_writer From xfel.command_line.cspad_detector_congruence: * iterate_detector_at_level → serialtbx.detector * iterate_panels → serialtbx.detector * id_from_name → serialtbx.detector * get_center → serialtbx.detector From xfel.command_line.frame_extractor: * ConstructFrame → serialtbx.util.construct_frame From xfel.cxi.cspad_ana.cspad_tbx: * Constants moved to serialtbx.detector.cspad * address_split → serialtbx.detector.xtc * dpack → serialtbx.detector.cspad * get_ebeam → serialtbx.detector.xtc * env_distance → serialtbx.detector.xtc * evt_wavelength → serialtbx.detector.xtc * env_detz → serialtbx.detector.xtc * old_address_to_new_address → serialtbx.detector.xtc From xfel.cxi.cspad_ana.rayonix_tbx.py: * Constants moved to serialtbx.detector.rayonix * get_rayonix_pixel_size → serialtbx.detector.rayonix * get_rayonix_detector_dimensions → serialtbx.detector.rayonix * get_data_from_psana_event → serialtbx.detector.rayonix C++ code from xfel/ext.cpp: * radial_average → dxtbx.ext From xfel.util: * sublattice_helper → dials.algorithms.integration Commits: * Move xfel.mono_simulation.max_like to serialtbx * Move time functionality from cspad_tbx to serialtbx/util/time.py * Move xfel.mono_simulation.util to serialtbx/mono_simulation * xfel.util.sublattice_helper moved to dials.algorithms.integration * Move xfel.command_line.frame_extractor.ConstructFrame to serialtbx.util * Move radial average c++ code to dxtbx * Move dpack to serialtbx.detector.cspad * Move jiffy functions from cspad_detector_congruence to serialtbx.detector * Move utility functions and constants out of cspad_tbx into serialtbx.detector.cspad and serialtbx.detector.xtc * Move movement of functions into serialtbx * Moving out of cspad_cbf_tbx: - Into serialtbx.detector: basis, center - Into serialtbx.detector.cspad: read_slac_metrology - Into serialtbx.detector.xtc: basis_from_geo * Use dxtbx's cbf_wrapper by adding in the only function in the xfel subclass of cbf_wrapper * Move add_frame_specific_cbf_tables to dxtbx * Move xfel.util.jungfrau to serialtbx.detector * cspad constants and data reading moved to serialtbx * Move rayonix functions and constants into serialtbx * Move old_address_to_new_address to serialtbx * Move legacy image pickle metrology code to serialtbx.detector.legacy_metrology * Missing imports and a few more functions moving into serialtbx * Temporarily check out serial_tbx branches on other repos * Duplicated code * Fix import * Add serialtbx to CCIBuilder configuration list. This should ensure serialtbx is included in the conda cctbx-base package * Revert "Temporarily check out serial_tbx branches on other repos" --------- Co-authored-by: Billy K. Poon --- libtbx/auto_build/bootstrap.py | 1 + rstbx/command_line/slip_viewer.py | 2 +- rstbx/slip_viewer/calibration_frame.py | 2 +- rstbx/slip_viewer/tile_generation.py | 6 +- serialtbx/detector/__init__.py | 124 ++++++ serialtbx/detector/cspad.py | 282 +++++++++++++ {xfel/util => serialtbx/detector}/jungfrau.py | 0 .../legacy_metrology}/cspad_detector.py | 14 +- .../legacy_metrology}/generic_detector.py | 0 .../detector/legacy_metrology}/metrology.py | 0 serialtbx/detector/rayonix.py | 56 +++ serialtbx/detector/xtc.py | 212 ++++++++++ .../mono_simulation/max_like.py | 0 {xfel => serialtbx}/mono_simulation/util.py | 0 serialtbx/util/construct_frame.py | 192 +++++++++ serialtbx/util/time.py | 28 ++ simtbx/command_line/hopper_process.py | 2 +- simtbx/diffBragg/hopper_utils.py | 2 +- xfel/cftbx/detector/cspad_cbf_tbx.py | 374 +----------------- xfel/cftbx/detector/metrology2phil.py | 6 +- xfel/clustering/singleframe.py | 2 +- xfel/command_line/cspad_cbf_metrology.py | 2 +- xfel/command_line/cspad_circular_gain_mask.py | 3 +- .../command_line/cspad_detector_congruence.py | 61 +-- xfel/command_line/cspad_detector_shifts.py | 6 +- xfel/command_line/cxi_calibdir2cbfheader.py | 3 +- xfel/command_line/cxi_display_metrology.py | 2 +- xfel/command_line/cxi_slaccalib2cbfheader.py | 3 +- xfel/command_line/detector_residuals.py | 3 +- xfel/command_line/frame_extractor.py | 189 +-------- xfel/command_line/quadrants_cbf.py | 2 +- xfel/command_line/radial_average.py | 2 +- xfel/cxi/cspad_ana/common_mode.py | 3 +- xfel/cxi/cspad_ana/cspad_tbx.py | 315 +-------------- xfel/cxi/cspad_ana/parse_calib.py | 2 +- xfel/cxi/cspad_ana/rayonix_tbx.py | 63 +-- xfel/ext.cpp | 37 -- xfel/mono_simulation/max_like_pos_neg.py | 2 +- xfel/mono_simulation/mono_treatment.py | 4 +- xfel/poly/recompute_mosaic_params.py | 2 +- xfel/swissfel/jf16m_cxigeom2nexus.py | 2 +- xfel/util/dials_pickle_reader.py | 2 +- xfel/util/sublattice_helper.py | 164 -------- 43 files changed, 965 insertions(+), 1212 deletions(-) create mode 100644 serialtbx/detector/__init__.py create mode 100644 serialtbx/detector/cspad.py rename {xfel/util => serialtbx/detector}/jungfrau.py (100%) rename {xfel/cftbx/detector => serialtbx/detector/legacy_metrology}/cspad_detector.py (96%) rename {xfel/cftbx/detector => serialtbx/detector/legacy_metrology}/generic_detector.py (100%) rename {xfel/cftbx/detector => serialtbx/detector/legacy_metrology}/metrology.py (100%) create mode 100644 serialtbx/detector/rayonix.py create mode 100644 serialtbx/detector/xtc.py rename {xfel => serialtbx}/mono_simulation/max_like.py (100%) rename {xfel => serialtbx}/mono_simulation/util.py (100%) create mode 100644 serialtbx/util/construct_frame.py create mode 100644 serialtbx/util/time.py delete mode 100644 xfel/util/sublattice_helper.py diff --git a/libtbx/auto_build/bootstrap.py b/libtbx/auto_build/bootstrap.py index 142148595c..b8f3c2578d 100644 --- a/libtbx/auto_build/bootstrap.py +++ b/libtbx/auto_build/bootstrap.py @@ -2146,6 +2146,7 @@ class CCTBXLiteBuilder(CCIBuilder): 'cctbx', 'cctbx_website', 'scitbx', + 'serialtbx', 'libtbx', 'iotbx', 'mmtbx', diff --git a/rstbx/command_line/slip_viewer.py b/rstbx/command_line/slip_viewer.py index 1d7d269459..2903f75542 100644 --- a/rstbx/command_line/slip_viewer.py +++ b/rstbx/command_line/slip_viewer.py @@ -69,7 +69,7 @@ def run(argv=None): frame.settings_frame.panel.collect_values() if (work_params.effective_metrology is not None): - from xfel.cftbx.detector.metrology import \ + from serialtbx.detector.legacy_metrology.metrology import \ master_phil, metrology_as_transformation_matrices stream = open(work_params.effective_metrology) diff --git a/rstbx/slip_viewer/calibration_frame.py b/rstbx/slip_viewer/calibration_frame.py index 1066a2569b..976d60114e 100644 --- a/rstbx/slip_viewer/calibration_frame.py +++ b/rstbx/slip_viewer/calibration_frame.py @@ -106,7 +106,7 @@ def OnRestoreMetrology(self, event): if dialog.ShowModal() == wx.ID_OK: path = dialog.GetPath() if (path != ""): - from xfel.cftbx.detector.metrology import \ + from serialtbx.detector.legacy_metrology.metrology import \ master_phil, metrology_as_transformation_matrices from libtbx import phil diff --git a/rstbx/slip_viewer/tile_generation.py b/rstbx/slip_viewer/tile_generation.py index a6a8dd6e6b..da71fe6775 100644 --- a/rstbx/slip_viewer/tile_generation.py +++ b/rstbx/slip_viewer/tile_generation.py @@ -50,8 +50,8 @@ def _get_flex_image( def _get_flex_image_multipanel(panels, raw_data, beam, brightness=1.0, binning=1, show_untrusted=False, color_scheme=0): - # From xfel.cftbx.cspad_detector.readHeader() and - # xfel.cftbx.cspad_detector.get_flex_image(). XXX Is it possible to + # From serialtbx.detector.legacy_metrology.cspad_detector.readHeader() and + # serialtbx.detector.legacy_metrology.cspad_detector.get_flex_image(). XXX Is it possible to # merge this with _get_flex_image() above? XXX Move to dxtbx Format # class (or a superclass for multipanel images)? @@ -61,7 +61,7 @@ def _get_flex_image_multipanel(panels, raw_data, beam, brightness=1.0, from libtbx.test_utils import approx_equal from scitbx.array_family import flex from scitbx.matrix import col, rec, sqr - from xfel.cftbx.detector.metrology import get_projection_matrix + from serialtbx.detector.legacy_metrology.metrology import get_projection_matrix assert len(panels) == len(raw_data), (len(panels), len(raw_data)) diff --git a/serialtbx/detector/__init__.py b/serialtbx/detector/__init__.py new file mode 100644 index 0000000000..91c9436b37 --- /dev/null +++ b/serialtbx/detector/__init__.py @@ -0,0 +1,124 @@ +from __future__ import division +from scitbx import matrix + +def center(coords): + """ Returns the average of a list of vectors + @param coords List of vectors to return the center of + """ + for c in coords: + if 'avg' not in locals(): + avg = c + else: + avg += c + return avg / len(coords) + +class basis(object): + """ Bucket for detector element information """ + def __init__(self, orientation = None, translation = None, panelgroup = None, homogenous_transformation = None, name = None): + """ + Provide only orientation + translation or a panelgroup or a homogenous_transformation. + + @param orientation rotation in the form of a quarternion + @param translation vector translation in relation to the parent frame + @param panelgroup dxtbx panelgroup object whose local d matrix will represent the + basis shift + @param homogenous_transformation 4x4 matrix.sqr object representing a translation + and a rotation. Must not also contain a scale as this won't be decomposed properly. + @param name optional name for this basis shift + """ + self.include_translation = True + self.name = name + + if panelgroup is not None: + d_mat = panelgroup.get_local_d_matrix() + fast = matrix.col((d_mat[0],d_mat[3],d_mat[6])).normalize() + slow = matrix.col((d_mat[1],d_mat[4],d_mat[7])).normalize() + orig = matrix.col((d_mat[2],d_mat[5],d_mat[8])) + + v3 = fast.cross(slow).normalize() + + r3 = matrix.sqr((fast[0],slow[0],v3[0], + fast[1],slow[1],v3[1], + fast[2],slow[2],v3[2])) + + self.orientation = r3.r3_rotation_matrix_as_unit_quaternion() + self.translation = orig + + if not self.name: + self.name = panelgroup.get_name() + + elif orientation is not None or translation is not None: + assert orientation is not None and translation is not None + self.orientation = orientation + self.translation = translation + + else: + # Decompose the homegenous transformation assuming no scale factors were used + h = homogenous_transformation + self.orientation = matrix.sqr((h[0],h[1],h[2], + h[4],h[5],h[6], + h[8],h[9],h[10])).r3_rotation_matrix_as_unit_quaternion() + self.translation = matrix.col((h[3], + h[7], + h[11])) + assert h[12] == h[13] == h[14] == 0 and h[15] == 1 + +def iterate_detector_at_level(item, depth = 0, level = 0): + """ + Iterate through all panel groups or panels of a detector object at a given + hierarchy level + @param item panel group or panel. Use detector.hierarchy(). + @param depth current dept for recursion. Should be 0 for initial call. + @param level iterate groups at this level + @return next panel or panel group object + """ + if level == depth: + yield item + else: + for child in item: + for subitem in iterate_detector_at_level(child, depth+1, level): + yield subitem + +def iterate_panels(panelgroup): + """ + Find and iterate all panels in the given panel group, regardless of the hierarchly level + of this panelgroup + @param panelgroup the panel group of interest + @return the next panel + """ + if panelgroup.is_group(): + for child in panelgroup: + for subitem in iterate_panels(child): + yield subitem + else: + yield panelgroup + +def id_from_name(detector, name): + """ Jiffy function to get the id of a panel using its name + @param detector detector object + @param name panel name + @return index of panel in detector + """ + return [p.get_name() for p in detector].index(name) + +def get_center(pg): + """ Find the center of a panel group pg, projected on its fast/slow plane """ + if pg.is_group(): + # find the average center of all this group's children + children_center = matrix.col((0,0,0)) + count = 0 + for p in iterate_panels(pg): + children_center += get_center(p) + count += 1 + children_center /= count + + # project the children center onto the plane of the panel group + pgf = matrix.col(pg.get_fast_axis()) + pgs = matrix.col(pg.get_slow_axis()) + pgn = matrix.col(pg.get_normal()) + pgo = matrix.col(pg.get_origin()) + + return (pgf.dot(children_center) * pgf) + (pgs.dot(children_center) * pgs) + (pgn.dot(pgo) * pgn) + else: + s = pg.get_image_size() + return matrix.col(pg.get_pixel_lab_coord((s[0]/2, s[1]/2))) diff --git a/serialtbx/detector/cspad.py b/serialtbx/detector/cspad.py new file mode 100644 index 0000000000..e411535e72 --- /dev/null +++ b/serialtbx/detector/cspad.py @@ -0,0 +1,282 @@ +from __future__ import division + +from libtbx.utils import Sorry +from scitbx import matrix +from scitbx.array_family import flex +import serialtbx.util.time +import os + +# The CAMP and CSpad counters are both 14 bits wide (Strüder et al +# 2010; Philipp et al., 2007), which means the physical limit is 2**14 - 1. +# However, in practice, when the pixels are in the low gain mode, after +# correcting by a gain value of around 6.87, the pixels tend to saturate +# around 90000. See xpp experiment xppe0314, run 184 as evidence. +cspad_saturated_value = 90000 + +# The dark average for the CSPAD detector is around 1100-1500. A pixel +# histogram of a minimum projection of an uncorrected (raw) light run shows +# a mostly flat tail up to ~800 ADU with a few bumps in the tail which +# represent true underloads. Assume a dark average of 1200 ADU. After dark +# subtraction, 800 - 1200 gives a minimum trusted value of -400. Reject +# pixels less than this. +cspad_min_trusted_value = -400 + +# The pixel size in mm. The pixel size is fixed and square, with side +# length of 110 µm (Philipp et al., 2007). XXX Should really clarify +# this with Sol and Chris. +# +# XXX Andor: 13.5 µm square, CAMP: 75 µm, square (Strüder et al., +# 2010) +# Apr 14 2023: commenting this out and using the more accurate version below +#pixel_size = 110e-3 + + +# need to define these here since it not defined in SLAC's metrology definitions +asic_dimension = (194,185) +asic_gap = 3 +pixel_size = 0.10992 + +PSANA2_VERSION = 0 +try: + PSANA2_VERSION = os.environ.get('PSANA2_VERSION', 0) +except AttributeError: + pass + +from serialtbx.detector import basis, center +from serialtbx.detector.xtc import basis_from_geo + +def read_slac_metrology(path = None, geometry = None, plot=False, include_asic_offset=False): + if path is None and geometry is None: + raise Sorry("Need to provide a geometry object or a path to a geometry file") + + if path is not None and geometry is not None: + raise Sorry("Cannot provide a geometry object and a geometry file. Ambiguous") + + if geometry is None: + try: + from PSCalib.GeometryAccess import GeometryAccess + geometry = GeometryAccess(path) + except Exception: + raise Sorry("Can't parse this metrology file") + + metro = {} + pixel_size = geometry.get_pixel_scale_size()/1000 + null_ori = matrix.col((0,0,1)).axis_and_angle_as_unit_quaternion(0, deg=True) + + # collapse any transformations above those of the quadrants into one X/Y offset, + # but don't keep Z transformations, as those come from the XTC stream + root = geometry.get_top_geo() + root_basis = basis_from_geo(root, use_z=False) + while len(root.get_list_of_children()) != 4 and len(root.get_list_of_children()) != 32: + assert len(root.get_list_of_children()) == 1 + root = root.get_list_of_children()[0] + root_basis *= basis_from_geo(root, use_z=False) + + metro[(0,)] = root_basis + + + + def add_sensor(quad_id, sensor_id, sensor): + metro[(0,quad_id,sensor_id)] = basis_from_geo(sensor) + + x, y, z = sensor.get_pixel_coords() + x/=1000; y/=1000; z/=1000 + assert x.shape == y.shape == z.shape + sensor_px_slow = x.shape[0] + sensor_px_fast = x.shape[1] + assert sensor_px_fast % 2 == 0 + + a0ul = sul = matrix.col((x[0,0],y[0,0],z[0,0])) + a1ur = sur = matrix.col((x[0,sensor_px_fast-1],y[0,sensor_px_fast-1],z[0,sensor_px_fast-1])) + a1lr = slr = matrix.col((x[sensor_px_slow-1,sensor_px_fast-1],y[sensor_px_slow-1,sensor_px_fast-1],z[sensor_px_slow-1,sensor_px_fast-1])) + a0ll = sll = matrix.col((x[sensor_px_slow-1,0],y[sensor_px_slow-1,0],z[sensor_px_slow-1,0])) + + a0ur = matrix.col((x[0,sensor_px_fast//2-1],y[0,sensor_px_fast//2-1],z[0,sensor_px_fast//2-1])) + a0lr = matrix.col((x[sensor_px_slow-1,sensor_px_fast//2-1],y[sensor_px_slow-1,sensor_px_fast//2-1],z[sensor_px_slow-1,sensor_px_fast//2-1])) + + a1ul = matrix.col((x[0,sensor_px_fast//2],y[0,sensor_px_fast//2],z[0,sensor_px_fast//2])) + a1ll = matrix.col((x[sensor_px_slow-1,sensor_px_fast//2],y[sensor_px_slow-1,sensor_px_fast//2],z[sensor_px_slow-1,sensor_px_fast//2])) + + sensor_center = center([sul,sur,slr,sll]) + asic0_center = center([a0ul,a0ur,a0lr,a0ll]) + asic1_center = center([a1ul,a1ur,a1lr,a1ll]) + + asic_trans0 = (asic0_center-sensor_center).length() + asic_trans1 = (asic1_center-sensor_center).length() + + if include_asic_offset: + rotated_ori = matrix.col((1,0,0)).axis_and_angle_as_unit_quaternion(180.0, deg=True) + offset_fast = -pixel_size*((sensor_px_fast) / 4) # 4 because sensor_px_fast is for sensor + offset_slow = +pixel_size*((sensor_px_slow) / 2) # Sensor is divided into 2 only in fast direction + metro[(0,quad_id,sensor_id,0)] = basis(orientation=rotated_ori,translation=matrix.col((-asic_trans0,0,0))) + metro[(0,quad_id,sensor_id,1)] = basis(orientation=rotated_ori,translation=matrix.col((+asic_trans1,0,0))) + metro[(0,quad_id,sensor_id,0)].translation += matrix.col((offset_fast, offset_slow, 0)) + metro[(0,quad_id,sensor_id,1)].translation += matrix.col((offset_fast, offset_slow, 0)) + else: + metro[(0,quad_id,sensor_id,0)] = basis(orientation=null_ori,translation=matrix.col((-asic_trans0,0,0))) + metro[(0,quad_id,sensor_id,1)] = basis(orientation=null_ori,translation=matrix.col((+asic_trans1,0,0))) + + if len(root.get_list_of_children()) == 4: + for quad_id, quad in enumerate(root.get_list_of_children()): + metro[(0,quad_id)] = basis_from_geo(quad) + for sensor_id, sensor in enumerate(quad.get_list_of_children()): + add_sensor(quad_id, sensor_id, sensor) + elif len(root.get_list_of_children()) == 32: + for quad_id in range(4): + metro[(0,quad_id)] = basis(orientation = null_ori, translation = matrix.col((0,0,0))) + sensors = root.get_list_of_children() + for sensor_id in range(8): + add_sensor(quad_id, sensor_id, sensors[quad_id*4+sensor_id]) + else: + assert False + + return metro + +def get_psana_corrected_data(psana_det, evt, use_default=False, dark=True, common_mode=None, apply_gain_mask=True, + gain_mask_value=None, per_pixel_gain=False, gain_mask=None, additional_gain_factor=None): + """ + Given a psana Detector object, apply corrections as appropriate and return the data from the event + @param psana_det psana Detector object + @param evt psana event + @param use_default If true, apply the default calibration only, using the psana algorithms. Otherise, use the corrections + specified by the rest of the flags and values passed in. + @param dark Whether to apply the detector dark, bool or numpy array + @param common_mode Which common mode algorithm to apply. None: apply no algorithm. Default: use the algorithm specified + in the calib folder. Otherwise should be a list as specified by the psana documentation for common mode customization + @param apply_gain_mask Whether to apply the common mode gain mask correction + @param gain_mask_value Multiplier to apply to the pixels, according to the gain mask + @param per_pixel_gain If available, use the per pixel gain deployed to the calibration folder + @param gain_mask gain mask showing which pixels to apply gain mask value + @param additional_gain_factor Additional gain factor. Pixels counts are divided by this number after all other + corrections. + @return Numpy array corrected as specified. + """ + # order is pedestals, then common mode, then gain mask, then per pixel gain + import numpy as np + + if PSANA2_VERSION: + # in psana2, data are stored as raw, fex, etc so the selection + # has to be given here when the detector interface is used. + # for now, assumes cctbx uses "raw". + psana_det = psana_det.raw + + if use_default: + return psana_det.calib(evt) # applies psana's complex run-dependent calibrations + data = psana_det.raw_data(evt) + if data is None: + return + + data = data.astype(np.float64) + if isinstance(dark, bool): + if dark: + if PSANA2_VERSION: + data -= psana_det.pedestals() + else: + data -= psana_det.pedestals(evt) + elif isinstance( dark, np.ndarray ): + data -= dark + + if common_mode is not None and common_mode != "default": + if common_mode == 'cspad_default': + common_mode = (1,25,25,100,1) # default parameters for CSPAD images + psana_det.common_mode_apply(data, common_mode) + elif common_mode == 'unbonded': + common_mode = (5,0,0,0,0) # unbonded pixels used for correction + psana_det.common_mode_apply(data, common_mode) + else: # this is how it was before.. Though I think common_mode would need to be a tuple.. + psana_det.common_mode_apply(data, common_mode) + if apply_gain_mask: + if gain_mask is None: # TODO: consider try/except here + gain_mask = psana_det.gain_mask(evt) == 1 + if gain_mask_value is None: + try: + gain_mask_value = psana_det._gain_mask_factor + except AttributeError: + print("No gain set for psana detector, using gain value of 1, consider disabling gain in your phil file") + gain_mask_value = 1 + data[gain_mask] = data[gain_mask]*gain_mask_value + if per_pixel_gain: # TODO: test this + data *= psana_det.gain() + if additional_gain_factor is not None: + data /= additional_gain_factor + return data + + + + + +def dpack(active_areas=None, + address=None, + beam_center_x=None, + beam_center_y=None, + ccd_image_saturation=None, + data=None, + distance=None, + pixel_size=pixel_size, + saturated_value=None, + timestamp=None, + wavelength=None, + xtal_target=None, + min_trusted_value=None): + """XXX Check completeness. Should fill in sensible defaults.""" + + # Must have data. + if data is None: + return None + + # Create a time stamp of the current time if none was supplied. + if timestamp is None: + timestamp = serialtbx.util.time.timestamp() + + # For unknown historical reasons, the dictionary must contain both + # CCD_IMAGE_SATURATION and SATURATED_VALUE items. + if ccd_image_saturation is None: + if saturated_value is None: + ccd_image_saturation = cspad_saturated_value + else: + ccd_image_saturation = saturated_value + if saturated_value is None: + saturated_value = ccd_image_saturation + + # Use a minimum value if provided for the pixel range + if min_trusted_value is None: + min_trusted_value = cspad_min_trusted_value + + # By default, the beam center is the center of the image. The slow + # (vertical) and fast (horizontal) axes correspond to x and y, + # respectively. + if beam_center_x is None: + beam_center_x = pixel_size * data.focus()[1] / 2 + if beam_center_y is None: + beam_center_y = pixel_size * data.focus()[0] / 2 + + # By default, the entire detector image is an active area. There is + # no sensible default for distance nor wavelength. XXX But setting + # wavelength to zero may be disastrous? + if active_areas is None: + # XXX Verify order with non-square detector + active_areas = flex.int((0, 0, data.focus()[0], data.focus()[1])) + if distance is None: + distance = 0 + if wavelength is None: + wavelength = 0 + + # The size must match the image dimensions. The length along the + # slow (vertical) axis is SIZE1, the length along the fast + # (horizontal) axis is SIZE2. + return {'ACTIVE_AREAS': active_areas, + 'BEAM_CENTER_X': beam_center_x, + 'BEAM_CENTER_Y': beam_center_y, + 'CCD_IMAGE_SATURATION': ccd_image_saturation, + 'DATA': data, + 'DETECTOR_ADDRESS': address, + 'DISTANCE': distance, + 'PIXEL_SIZE': pixel_size, + 'SATURATED_VALUE': saturated_value, + 'MIN_TRUSTED_VALUE': min_trusted_value, + 'SIZE1': data.focus()[0], + 'SIZE2': data.focus()[1], + 'TIMESTAMP': timestamp, + 'SEQUENCE_NUMBER': 0, # XXX Deprecated + 'WAVELENGTH': wavelength, + 'xtal_target': xtal_target} diff --git a/xfel/util/jungfrau.py b/serialtbx/detector/jungfrau.py similarity index 100% rename from xfel/util/jungfrau.py rename to serialtbx/detector/jungfrau.py diff --git a/xfel/cftbx/detector/cspad_detector.py b/serialtbx/detector/legacy_metrology/cspad_detector.py similarity index 96% rename from xfel/cftbx/detector/cspad_detector.py rename to serialtbx/detector/legacy_metrology/cspad_detector.py index 3263d684a4..9da6b58885 100644 --- a/xfel/cftbx/detector/cspad_detector.py +++ b/serialtbx/detector/legacy_metrology/cspad_detector.py @@ -10,7 +10,7 @@ from libtbx import easy_pickle from scitbx.array_family import flex from scitbx.matrix import col, rec, sqr -from xfel.cftbx.detector.generic_detector import GenericDetector +from serialtbx.detector.legacy_metrology.generic_detector import GenericDetector from six.moves import range import six @@ -43,7 +43,7 @@ def readHeader(self): # However, several code paths still depend on the member variables # created here. - from xfel.cftbx.detector.metrology import metrology_as_transformation_matrices + from serialtbx.detector.legacy_metrology.metrology import metrology_as_transformation_matrices d = easy_pickle.load(self.filename) @@ -117,11 +117,11 @@ def get_flex_image(self, brightness, **kwargs): # rstbx.slip_viewer.tile_generation._get_flex_image_multitile(). # XXX Still used by iotbx/command_line/detector_image_as_png.py #raise DeprecationWarning( - # "xfel.cftbx.cspad_detector.get_flex_image() is deprecated") + # "serialtbx.detector.legacy_metrology.cspad_detector.get_flex_image() is deprecated") # no kwargs supported at present - from xfel.cftbx.detector.metrology import get_projection_matrix + from serialtbx.detector.legacy_metrology.metrology import get_projection_matrix # E maps picture coordinates onto metric Cartesian coordinates, # i.e. [row, column, 1 ] -> [x, y, z, 1]. Both frames share the @@ -183,7 +183,7 @@ def get_panel_fast_slow(self, serial): directly. """ - from xfel.cftbx.detector.metrology import get_projection_matrix + from serialtbx.detector.legacy_metrology.metrology import get_projection_matrix center = col([self._asic_focus[0] / 2, self._asic_focus[1] / 2, 1]) fast, nmemb, slow = 0, 0, 0 @@ -269,7 +269,7 @@ def transformation_matrices_as_metrology(self): """ from libtbx import phil - from xfel.cftbx.detector.metrology import \ + from serialtbx.detector.legacy_metrology.metrology import \ master_phil, regularize_transformation_matrices # XXX Experimental! @@ -323,7 +323,7 @@ def readout_coords_as_detector_coords(self, coords): assert self._pixel_size is not None assert self._asic_focus is not None - from xfel.cftbx.detector.metrology import get_projection_matrix + from serialtbx.detector.legacy_metrology.metrology import get_projection_matrix P_f, P_b = get_projection_matrix(self._pixel_size, self._asic_focus) diff --git a/xfel/cftbx/detector/generic_detector.py b/serialtbx/detector/legacy_metrology/generic_detector.py similarity index 100% rename from xfel/cftbx/detector/generic_detector.py rename to serialtbx/detector/legacy_metrology/generic_detector.py diff --git a/xfel/cftbx/detector/metrology.py b/serialtbx/detector/legacy_metrology/metrology.py similarity index 100% rename from xfel/cftbx/detector/metrology.py rename to serialtbx/detector/legacy_metrology/metrology.py diff --git a/serialtbx/detector/rayonix.py b/serialtbx/detector/rayonix.py new file mode 100644 index 0000000000..20a762af2b --- /dev/null +++ b/serialtbx/detector/rayonix.py @@ -0,0 +1,56 @@ +from __future__ import division + +# given value of rayonix detector saturation xppi6115 +rayonix_saturated_value = 2**16 -1 + +# minimum value for rayonix data +rayonix_min_trusted_value = 0 +rayonix_max_trusted_value = rayonix_saturated_value - 1 + +def get_rayonix_pixel_size(bin_size): + ''' Given a bin size determine a pixel size. + +Michael Blum from Rayonix said The pixel size is recorded in the header, +but can be derived trivially from the overall dimension of the corrected imaging +area (170mm) and the number of pixels. (3840 unbinned). The corrected image is +forced to this size. + +unbinned 170/3840 = 0.04427 + +I believe the accuracy of the MEAN pixel size to be at least as good as 0.1% +which is the limit to which I can measure our calibration plate and exceeds the + parallax error in our calibration station. + +Note, the Rayonix MX340 has the same pixel size as the MX170: + +unbinned 340/7680 = 0.04427 + + @param bin_size rayonix bin size as an integer + ''' + pixel_size=bin_size*170/3840 + return pixel_size + +def get_rayonix_detector_dimensions(env): + ''' Given a psana env object, find the detector dimensions + @param env psana environment object + ''' + import psana + cfgs = env.configStore() + rayonix_cfg = cfgs.get(psana.Rayonix.ConfigV2, psana.Source('Rayonix')) + if not rayonix_cfg: return None, None + return rayonix_cfg.width(), rayonix_cfg.height() + +def get_data_from_psana_event(evt, address): + """ Read the pixel data for a Rayonix image from an event + @param psana event object + @param address old style psana detector address + @return numpy array with raw data""" + from psana import Source, Camera + from serialtbx.detector import xtc + import numpy as np + address = xtc.old_address_to_new_address(address) + src=Source(address) + data = evt.get(Camera.FrameV1,src) + if data is not None: + data = data.data16().astype(np.float64) + return data diff --git a/serialtbx/detector/xtc.py b/serialtbx/detector/xtc.py new file mode 100644 index 0000000000..d25f282a97 --- /dev/null +++ b/serialtbx/detector/xtc.py @@ -0,0 +1,212 @@ +from __future__ import division + +from scitbx import matrix +from serialtbx.detector import basis + +def basis_from_geo(geo, use_z = True): + """ Given a psana GeometryObject, construct a basis object """ + rotx = matrix.col((1,0,0)).axis_and_angle_as_r3_rotation_matrix( + geo.rot_x + geo.tilt_x, deg=True) + roty = matrix.col((0,1,0)).axis_and_angle_as_r3_rotation_matrix( + geo.rot_y + geo.tilt_y, deg=True) + rotz = matrix.col((0,0,1)).axis_and_angle_as_r3_rotation_matrix( + geo.rot_z + geo.tilt_z, deg=True) + + rot = (rotx*roty*rotz).r3_rotation_matrix_as_unit_quaternion() + + if use_z: + trans = matrix.col((geo.x0/1000, geo.y0/1000, geo.z0/1000)) + else: + trans = matrix.col((geo.x0/1000, geo.y0/1000, 0)) + + return basis(orientation = rot, translation = trans) + +def old_address_to_new_address(address): + """ Change between old and new style detector addresses. + I.E. CxiDs1-0|Cspad-0 becomes CxiDs1.0:Cspad.0 + @param address detector address to convert + """ + return address.replace('-','.').replace('|',':') + +def address_split(address, env=None): + """The address_split() function splits an address into its four + components. Address strings are on the form + detector-detectorID|device-deviceID, where the detectors must be in + dir(xtc.DetInfo.Detector) and device must be in + (xtc.DetInfo.Device). + @param address Full data source address of the DAQ device + @param env Optional env to dereference an alias into an address + @return Four-tuple of detector name, detector ID, device, and + device ID + """ + + import re + + # pyana + m = re.match( + r"^(?P\S+)\-(?P\d+)\|(?P\S+)\-(?P\d+)$", address) + if m is not None: + return (m.group('det'), m.group('det_id'), m.group('dev'), m.group('dev_id')) + + # psana + m = re.match( + r"^(?P\S+)\.(?P\d+)\:(?P\S+)\.(?P\d+)$", address) + if m is not None: + return (m.group('det'), m.group('det_id'), m.group('dev'), m.group('dev_id')) + + # psana DetInfo string + m = re.match( + r"^DetInfo\((?P\S+)\.(?P\d+)\:(?P\S+)\.(?P\d+)\)$", address) + if m is not None: + return (m.group('det'), m.group('det_id'), m.group('dev'), m.group('dev_id')) + + if env is not None: + # Try to see if this is a detector alias, and if so, dereference it. Code from psana's Detector/PyDetector.py + amap = env.aliasMap() + alias_src = amap.src(address) # string --> DAQ-style psana.Src + + # if it is an alias, look up the full name + if amap.alias(alias_src) != '': # alias found + address = str(alias_src) + return address_split(address) + + return (None, None, None, None) + +def get_ebeam(evt): + try: + # pyana + ebeam = evt.getEBeam() + except AttributeError: + from psana import Source, Bld + src = Source('BldInfo(EBeam)') + ebeam = evt.get(Bld.BldDataEBeamV6, src) + if ebeam is None: + ebeam = evt.get(Bld.BldDataEBeamV5, src) + if ebeam is None: + ebeam = evt.get(Bld.BldDataEBeamV4, src) + if ebeam is None: + ebeam = evt.get(Bld.BldDataEBeamV3, src) + if ebeam is None: + ebeam = evt.get(Bld.BldDataEBeamV2, src) + if ebeam is None: + ebeam = evt.get(Bld.BldDataEBeamV1, src) + if ebeam is None: + ebeam = evt.get(Bld.BldDataEBeamV0, src) + if ebeam is None: + ebeam = evt.get(Bld.BldDataEBeam, src) # recent version of psana will return a V7 event or higher if this type is asked for + + return ebeam + +def evt_wavelength(evt, delta_k=0): + """The evt_wavelength() function returns the wavelength in Ångström + of the event pointed to by @p evt. From Margaritondo & Rebernik + Ribic (2011): the dimensionless relativistic γ-factor is derived + from beam energy in MeV and the electron rest mass, K is a + dimensionless "undulator parameter", and L is the macroscopic + undulator period in Ångström. See also + https://people.eecs.berkeley.edu/~attwood/srms/2007/Lec10.pdf + + @param evt Event data object, a configure object + @param delta_k Optional K-value correction + @return Wavelength, in Ångström + """ + + if evt is not None: + ebeam = get_ebeam(evt) + + if hasattr(ebeam, 'fEbeamPhotonEnergy') and ebeam.fEbeamPhotonEnergy > 0: + # pyana + return 12398.4187 / ebeam.fEbeamPhotonEnergy + if hasattr(ebeam, 'ebeamPhotonEnergy') and ebeam.ebeamPhotonEnergy() > 0: + # psana + return 12398.4187 / ebeam.ebeamPhotonEnergy() + + if hasattr(ebeam, 'fEbeamL3Energy') and ebeam.fEbeamL3Energy > 0: + # pyana + gamma = ebeam.fEbeamL3Energy / 0.510998910 + elif hasattr(ebeam, 'ebeamL3Energy') and ebeam.ebeamL3Energy() > 0: + # psana + gamma = ebeam.ebeamL3Energy() / 0.510998910 + else: + return None + K = 3.5 + delta_k + L = 3.0e8 + return L / (2 * gamma**2) * (1 + K**2 / 2) + return None + +def env_detz(address, env): + """The env_detz() function returns the position of the detector with + the given address string on the z-axis in mm. The zero-point is as + far away as possible from the sample, and values decrease as the + detector is moved towards the sample. + @param address Full data source address of the DAQ device + @param env Environment object + @return Detector z-position, in mm + """ + + if env is not None: + detector = address_split(address, env)[0] + if detector is None: + return None + elif detector == 'CxiDs1': + pv = env.epicsStore().value('CXI:DS1:MMS:06.RBV') + if pv is None: + # Even though potentially unsafe, fall back on the commanded + # value if the corresponding read-back value cannot be read. + # According to Sébastien Boutet, this particular motor has not + # caused any problem in the past. + pv = env.epicsStore().value('CXI:DS1:MMS:06') + if pv is None: + # Try the other detector. These are sometimes inconsistent + pv = env.epicsStore().value('CXI:DS2:MMS:06.RBV') + elif detector == 'CxiDsd' or detector == 'CxiDs2': + # XXX Note inconsistency in naming: Dsd vs Ds2! + pv = env.epicsStore().value('CXI:DS2:MMS:06.RBV') + if pv is None: + # Try the other detector. These are sometimes inconsistent + pv = env.epicsStore().value('CXI:DS1:MMS:06.RBV') + elif detector == 'XppGon': + # There is no distance recorded for the XPP's CSPAD on the robot + # arm. Always return zero to allow the distance to be set using + # the offset. + return 0 + elif detector == 'XppEndstation' or \ + detector == 'MfxEndstation': + # There is no distance recorded for the XPP's or MFX's Rayonix + # on the robot arm. Always return zero to allow the distance to + # be set using the offset. + return 0 + else: + return None + + if pv is None: + return None + + if hasattr(pv, "values"): + if len(pv.values) == 1: + return pv.values[0] + else: + return None + return pv + + return None + +def env_distance(address, env, offset): + """The env_distance() function returns the distance between the + sample and the detector with the given address string in mm. The + distance between the sample and the the detector's zero-point can + vary by an inch or more between different LCLS runs. According to + Sébastien Boutet the offset should be stable to within ±0.5 mm + during a normal experiment. + + @param address Full data source address of the DAQ device + @param env Environment object + @param offset Detector-sample offset in mm, corresponding to + longest detector-sample distance + @return Detector-sample distance, in mm + """ + + detz = env_detz(address, env) + if detz is not None: + return detz + offset + return None diff --git a/xfel/mono_simulation/max_like.py b/serialtbx/mono_simulation/max_like.py similarity index 100% rename from xfel/mono_simulation/max_like.py rename to serialtbx/mono_simulation/max_like.py diff --git a/xfel/mono_simulation/util.py b/serialtbx/mono_simulation/util.py similarity index 100% rename from xfel/mono_simulation/util.py rename to serialtbx/mono_simulation/util.py diff --git a/serialtbx/util/construct_frame.py b/serialtbx/util/construct_frame.py new file mode 100644 index 0000000000..3fed896624 --- /dev/null +++ b/serialtbx/util/construct_frame.py @@ -0,0 +1,192 @@ +from __future__ import absolute_import, division, print_function +from six.moves import range + +from scitbx.array_family import flex +from cctbx import crystal, miller +from cctbx.crystal_orientation import crystal_orientation +import cctbx + +class ConstructFrame(object): + def get_template_pickle(self): + return { + 'beam_s0': 0, + 'beam_polarization_normal': 0, + 'current_cb_op_to_primitive': 0, + 'current_orientation':0, + 'distance':0, + 'effective_tiling':0, + 'identified_isoform':None, + 'mapped_predictions':[[]], + 'max_signal':0, + 'ML_domain_size_ang':[0], + 'ML_half_mosaicity_deg':[0], + 'mosaicity':0, + 'model_partialities':[None], + 'observations':[0], + 'pointgroup':0, + 'residual':0, + 'sa_parameters':['None'], + 's1_vec':[0], + 'wavelength':0, + 'xbeam':0, + 'ybeam':0} + + def __init__(self, reflections, experiment): + # assemble template and unpack files + self.frame = self.get_template_pickle() + self.pixel_size = experiment.detector[0].get_pixel_size()[0] + + if 'intensity.prf.value' in reflections: + self.method = 'prf' # integration by profile fitting + elif 'intensity.sum.value' in reflections: + self.method = 'sum' # integration by simple summation + self.reflections = reflections.select(reflections['intensity.' + self.method + '.variance'] > 0) # keep only spots with sigmas above zero + + self.xtal = experiment.crystal + self.beam_obj = experiment.beam + self.det = experiment.detector + self.gonio = experiment.goniometer + self.scan = experiment.scan + self.img_sweep = experiment.imageset + + # experiment-dependent components --------------------------------------------------------------------------- + + # get wavelength + def populate_wavelength(self): + assert self.beam_obj.get_wavelength() is not None, "no wavelength" + self.frame['wavelength'] = self.beam_obj.get_wavelength() + self.frame['beam_s0'] = self.beam_obj.get_s0() + self.frame['beam_polarization_normal'] = self.beam_obj.get_polarization_normal() + + # get detector distance in mm + def populate_distance(self): + assert self.det[0].get_distance() is not None, "no detector distance" + self.frame['distance'] = self.det[0].get_distance() + + # get xbeam and ybeam in mm + def populate_beam_dir(self): + assert self.beam_obj.get_s0() is not None, "no beam direction" + self.frame['xbeam'], self.frame['ybeam'] = self.det[0].get_beam_centre(self.beam_obj.get_s0()) + + # get isoform + def populate_isoform(self): + if hasattr(self.xtal, "identified_isoform"): + self.frame['identified_isoform'] = self.xtal.identified_isoform + + # get max signal + def populate_max_signal(self): + pass + + # get effective tiling + def populate_effective_tiling(self): + pass + + # indicate simulated annealing parameters, if present + def populate_sa_params(self): + pass + + # crystal-dependent components ------------------------------------------------------------------------------ + + # generate a crystal orientation object from the A* matrix + def populate_orientation(self): + assert self.xtal.get_A() is not None, "no crystal orientation matrix" + self.frame['current_orientation'] = [crystal_orientation(self.xtal.get_A(), True)] + + # generate change-of-basis operation for current to primitive cell + def populate_op_to_primitive(self): + assert self.xtal.get_space_group() is not None, "no space group" + self.frame['current_cb_op_to_primitive'] = [self.xtal.get_space_group().z2p_op()] + + # fetch the point group associated with the crystal + def populate_point_group(self): + assert self.xtal.get_space_group() is not None, "no space group" + self.frame['pointgroup'] = str(self.xtal.get_space_group().build_derived_point_group().info()) + + # get mosaicity + def populate_mosaicity(self): + try: + assert self.xtal.get_mosaicity() is not None, "no mosaicity" + self.frame['mosaicity'] = self.xtal.get_mosaicity() + except Exception: + pass + + # get any available ML values + def populate_ML_values(self): + try: + self.frame['ML_half_mosaicity_deg'] = [self.xtal.get_half_mosaicity_deg()] + except AttributeError: + pass + try: + self.frame['ML_domain_size_ang'] = [self.xtal.get_domain_size_ang()] + except AttributeError: + pass + + # observations-dependent components ------------------------------------------------------------------------- + + # generate a miller array containing the Miller indices, intensities and variances for one frame + def populate_observations(self): + intensities = self.reflections['intensity.' + self.method + '.value'] + variances = self.reflections['intensity.' + self.method + '.variance'] + space_group = crystal.symmetry(self.xtal.get_unit_cell(), str(self.xtal.get_space_group().info())) + miller_set = miller.set(space_group, self.reflections['miller_index']) + self.frame['observations'][0] = cctbx.miller.array(miller_set, intensities, flex.sqrt(variances)).set_observation_type_xray_intensity() + self.frame['s1_vec'][0] = self.reflections['s1'] + + # collect predicted spot positions + def populate_pixel_positions(self): + assert 'xyzcal.px' in self.reflections, "no calculated spot positions" + self.frame['mapped_predictions'][0] = flex.vec2_double() + for i in range(len(self.reflections['xyzcal.px'])): + self.frame['mapped_predictions'][0].append(tuple(self.reflections['xyzcal.px'][i][1::-1])) # 1::-1 reverses the order taking only the first two elements first. + + # generate a list of dictionaries containing a series of corrections for each predicted reflection + def populate_corrections(self): + assert 'xyzobs.px.value' in self.reflections and 'xyzcal.px' in self.reflections, "no calculated or observed spot positions" + assert self.frame['xbeam'] != 0 and self.frame['ybeam'] != 0, "invalid beam center" + self.frame['correction_vectors'] = [[]] + for idx in range(len(self.reflections['xyzobs.px.value'])): + if self.reflections['xyzcal.px'][idx][0:2] != self.reflections['xyzobs.px.value'][idx][0:2]: + theoret_center = 1765/2, 1765/2 + refined_center = self.frame['xbeam']/self.pixel_size, self.frame['ybeam']/self.pixel_size # px to mm conversion + hkl = self.reflections['miller_index'][idx] + obsspot = tuple(self.reflections['xyzobs.px.value'][idx][0:2]) + predspot = tuple(self.reflections['xyzcal.px'][idx][0:2]) + self.frame['correction_vectors'][0].append({'refinedcenter':refined_center, 'hkl':hkl, 'setting_id':0, 'azimuthal':0, 'radial':0, + 'obsspot':obsspot, 'obscenter':theoret_center, 'predspot':predspot}) + + # get partialities + def populate_partialities(self): + pass + + # produce residuals + def populate_residuals(self): + pass + + # collect kapton corrections if present + def populate_kapton_corrections(self): + if 'kapton_absorption_correction' in self.reflections: + self.frame['fuller_kapton_absorption_correction'] = [self.reflections['kapton_absorption_correction']] + if 'kapton_absorption_correction_sigmas' in self.reflections: + self.frame['fuller_kapton_absorption_correction_sigmas'] = [self.reflections['kapton_absorption_correction_sigmas']] + + # combine all of the above + def make_frame(self): + self.populate_wavelength() + self.populate_distance() + self.populate_beam_dir() + self.populate_isoform() + self.populate_max_signal() + self.populate_effective_tiling() + self.populate_sa_params() + self.populate_orientation() + self.populate_op_to_primitive() + self.populate_point_group() + self.populate_mosaicity() + self.populate_ML_values() + self.populate_observations() + self.populate_pixel_positions() + # self.populate_corrections() # works, but unnecessary + self.populate_partialities() + self.populate_residuals() + self.populate_kapton_corrections() + return self.frame diff --git a/serialtbx/util/time.py b/serialtbx/util/time.py new file mode 100644 index 0000000000..87d2393a82 --- /dev/null +++ b/serialtbx/util/time.py @@ -0,0 +1,28 @@ +from __future__ import division + +import math, time + +def now_s_ms(): + """The now_s_ms() function returns the time since + midnight, 1 January 1970 UTC (Unix time) to millisecond precision. + + @return Unix time as a tuple of seconds and milliseconds + """ + + t = time.time() + s = int(math.floor(t)) + return (s, int(round((t - s) * 1000))) + +def timestamp(t=None): + """The timestamp() function returns a string representation of + an extended human-readable ISO 8601 timestamp. If @p t is @c None + the current time is used. + + @param t Tuple of the time in seconds and milliseconds + @return Human-readable ISO 8601 timestamp in string representation + """ + + if t is None: + t = now_s_ms() + return time.strftime("%Y-%m-%dT%H:%MZ%S", time.gmtime(t[0])) + \ + (".%03d" % t[1]) diff --git a/simtbx/command_line/hopper_process.py b/simtbx/command_line/hopper_process.py index f170c89018..f04d136b72 100644 --- a/simtbx/command_line/hopper_process.py +++ b/simtbx/command_line/hopper_process.py @@ -256,7 +256,7 @@ def integrate(self, experiments, indexed): ) if self.params.dispatch.coset: - from xfel.util.sublattice_helper import integrate_coset + from dials.algorithms.integration.sublattice_helper import integrate_coset integrate_coset(self, experiments, indexed) diff --git a/simtbx/diffBragg/hopper_utils.py b/simtbx/diffBragg/hopper_utils.py index 8d44789a8b..1b29754898 100644 --- a/simtbx/diffBragg/hopper_utils.py +++ b/simtbx/diffBragg/hopper_utils.py @@ -14,7 +14,7 @@ from scipy.ndimage import binary_dilation from dxtbx.model.experiment_list import ExperimentListFactory from dxtbx.model import Spectrum -from xfel.util.jungfrau import get_pedestalRMS_from_jungfrau +from serialtbx.detector.jungfrau import get_pedestalRMS_from_jungfrau from simtbx.nanoBragg.utils import downsample_spectrum from dials.array_family import flex from simtbx.diffBragg import utils diff --git a/xfel/cftbx/detector/cspad_cbf_tbx.py b/xfel/cftbx/detector/cspad_cbf_tbx.py index da6d9a94a0..e288e4df7a 100644 --- a/xfel/cftbx/detector/cspad_cbf_tbx.py +++ b/xfel/cftbx/detector/cspad_cbf_tbx.py @@ -18,304 +18,10 @@ except AttributeError: pass -# need to define these here since it not defined in SLAC's metrology definitions -asic_dimension = (194,185) -asic_gap = 3 -pixel_size = 0.10992 -from xfel.cxi.cspad_ana.cspad_tbx import cspad_saturated_value, cspad_min_trusted_value - -def get_psana_corrected_data(psana_det, evt, use_default=False, dark=True, common_mode=None, apply_gain_mask=True, - gain_mask_value=None, per_pixel_gain=False, gain_mask=None, additional_gain_factor=None): - """ - Given a psana Detector object, apply corrections as appropriate and return the data from the event - @param psana_det psana Detector object - @param evt psana event - @param use_default If true, apply the default calibration only, using the psana algorithms. Otherise, use the corrections - specified by the rest of the flags and values passed in. - @param dark Whether to apply the detector dark, bool or numpy array - @param common_mode Which common mode algorithm to apply. None: apply no algorithm. Default: use the algorithm specified - in the calib folder. Otherwise should be a list as specified by the psana documentation for common mode customization - @param apply_gain_mask Whether to apply the common mode gain mask correction - @param gain_mask_value Multiplier to apply to the pixels, according to the gain mask - @param per_pixel_gain If available, use the per pixel gain deployed to the calibration folder - @param gain_mask gain mask showing which pixels to apply gain mask value - @param additional_gain_factor Additional gain factor. Pixels counts are divided by this number after all other - corrections. - @return Numpy array corrected as specified. - """ - # order is pedestals, then common mode, then gain mask, then per pixel gain - import numpy as np - - if PSANA2_VERSION: - # in psana2, data are stored as raw, fex, etc so the selection - # has to be given here when the detector interface is used. - # for now, assumes cctbx uses "raw". - psana_det = psana_det.raw - - if use_default: - return psana_det.calib(evt) # applies psana's complex run-dependent calibrations - data = psana_det.raw_data(evt) - if data is None: - return - - data = data.astype(np.float64) - if isinstance(dark, bool): - if dark: - if PSANA2_VERSION: - data -= psana_det.pedestals() - else: - data -= psana_det.pedestals(evt) - elif isinstance( dark, np.ndarray ): - data -= dark - - if common_mode is not None and common_mode != "default": - if common_mode == 'cspad_default': - common_mode = (1,25,25,100,1) # default parameters for CSPAD images - psana_det.common_mode_apply(data, common_mode) - elif common_mode == 'unbonded': - common_mode = (5,0,0,0,0) # unbonded pixels used for correction - psana_det.common_mode_apply(data, common_mode) - else: # this is how it was before.. Though I think common_mode would need to be a tuple.. - psana_det.common_mode_apply(data, common_mode) - if apply_gain_mask: - if gain_mask is None: # TODO: consider try/except here - gain_mask = psana_det.gain_mask(evt) == 1 - if gain_mask_value is None: - try: - gain_mask_value = psana_det._gain_mask_factor - except AttributeError: - print("No gain set for psana detector, using gain value of 1, consider disabling gain in your phil file") - gain_mask_value = 1 - data[gain_mask] = data[gain_mask]*gain_mask_value - if per_pixel_gain: # TODO: test this - data *= psana_det.gain() - if additional_gain_factor is not None: - data /= additional_gain_factor - return data - -from dxtbx.format.FormatCBFMultiTile import cbf_wrapper as dxtbx_cbf_wrapper -class cbf_wrapper(dxtbx_cbf_wrapper): - """ Wrapper class that provides convenience functions for working with cbflib""" - - def add_frame_shift(self, basis, axis_settings): - """Add an axis representing a frame shift (a rotation axis with an offset)""" - angle, axis = angle_and_axis(basis) - - if angle == 0: - axis = (0,0,1) - - if basis.include_translation: - translation = basis.translation - else: - translation = (0,0,0) - - self.add_row([basis.axis_name,"rotation","detector",basis.depends_on, - str(axis[0]),str(axis[1]),str(axis[2]), - str(translation[0]), - str(translation[1]), - str(translation[2]), - basis.equipment_component]) - - axis_settings.append([basis.axis_name, "FRAME1", str(angle), "0"]) - -def angle_and_axis(basis): - """Normalize a quaternion and return the angle and axis - @param params metrology object""" - q = matrix.col(basis.orientation).normalize() - return q.unit_quaternion_as_axis_and_angle(deg=True) - -def center(coords): - """ Returns the average of a list of vectors - @param coords List of vectors to return the center of - """ - for c in coords: - if 'avg' not in locals(): - avg = c - else: - avg += c - return avg / len(coords) - -class basis(object): - """ Bucket for detector element information """ - def __init__(self, orientation = None, translation = None, panelgroup = None, homogenous_transformation = None, name = None): - """ - Provide only orientation + translation or a panelgroup or a homogenous_transformation. - - @param orientation rotation in the form of a quarternion - @param translation vector translation in relation to the parent frame - @param panelgroup dxtbx panelgroup object whose local d matrix will represent the - basis shift - @param homogenous_transformation 4x4 matrix.sqr object representing a translation - and a rotation. Must not also contain a scale as this won't be decomposed properly. - @param name optional name for this basis shift - """ - self.include_translation = True - self.name = name - - if panelgroup is not None: - d_mat = panelgroup.get_local_d_matrix() - fast = matrix.col((d_mat[0],d_mat[3],d_mat[6])).normalize() - slow = matrix.col((d_mat[1],d_mat[4],d_mat[7])).normalize() - orig = matrix.col((d_mat[2],d_mat[5],d_mat[8])) - - v3 = fast.cross(slow).normalize() - - r3 = matrix.sqr((fast[0],slow[0],v3[0], - fast[1],slow[1],v3[1], - fast[2],slow[2],v3[2])) - - self.orientation = r3.r3_rotation_matrix_as_unit_quaternion() - self.translation = orig - - if not self.name: - self.name = panelgroup.get_name() - - elif orientation is not None or translation is not None: - assert orientation is not None and translation is not None - self.orientation = orientation - self.translation = translation - - else: - # Decompose the homegenous transformation assuming no scale factors were used - h = homogenous_transformation - self.orientation = matrix.sqr((h[0],h[1],h[2], - h[4],h[5],h[6], - h[8],h[9],h[10])).r3_rotation_matrix_as_unit_quaternion() - self.translation = matrix.col((h[3], - h[7], - h[11])) - assert h[12] == h[13] == h[14] == 0 and h[15] == 1 - - def as_homogenous_transformation(self): - """ Returns this basis change as a 4x4 transformation matrix in homogenous coordinates""" - r3 = self.orientation.normalize().unit_quaternion_as_r3_rotation_matrix() - return matrix.sqr((r3[0],r3[1],r3[2],self.translation[0], - r3[3],r3[4],r3[5],self.translation[1], - r3[6],r3[7],r3[8],self.translation[2], - 0,0,0,1)) - - def __mul__(self, other): - """ Use homogenous matrices to multiply bases together """ - if hasattr(other, 'as_homogenous_transformation'): - return basis(homogenous_transformation = self.as_homogenous_transformation() * other.as_homogenous_transformation()) - elif hasattr(other, 'n'): - if other.n == (3,1): - b = matrix.col((other[0], other[1], other[2], 1)) - elif other.n == (4,1): - b = other - else: - raise TypeError(b, "Incompatible matrices") - p = self.as_homogenous_transformation() * b - if other.n == (3,1): - return matrix.col(p[0:3]) - else: - return p - else: - raise TypeError(b) - - -def basis_from_geo(geo, use_z = True): - """ Given a psana GeometryObject, construct a basis object """ - rotx = matrix.col((1,0,0)).axis_and_angle_as_r3_rotation_matrix( - geo.rot_x + geo.tilt_x, deg=True) - roty = matrix.col((0,1,0)).axis_and_angle_as_r3_rotation_matrix( - geo.rot_y + geo.tilt_y, deg=True) - rotz = matrix.col((0,0,1)).axis_and_angle_as_r3_rotation_matrix( - geo.rot_z + geo.tilt_z, deg=True) - - rot = (rotx*roty*rotz).r3_rotation_matrix_as_unit_quaternion() - - if use_z: - trans = matrix.col((geo.x0/1000, geo.y0/1000, geo.z0/1000)) - else: - trans = matrix.col((geo.x0/1000, geo.y0/1000, 0)) - - return basis(orientation = rot, translation = trans) - -def read_slac_metrology(path = None, geometry = None, plot=False, include_asic_offset=False): - if path is None and geometry is None: - raise Sorry("Need to provide a geometry object or a path to a geometry file") - - if path is not None and geometry is not None: - raise Sorry("Cannot provide a geometry object and a geometry file. Ambiguous") - - if geometry is None: - try: - from PSCalib.GeometryAccess import GeometryAccess - geometry = GeometryAccess(path) - except Exception as e: - raise Sorry("Can't parse this metrology file") - - metro = {} - pixel_size = geometry.get_pixel_scale_size()/1000 - null_ori = matrix.col((0,0,1)).axis_and_angle_as_unit_quaternion(0, deg=True) - - # collapse any transformations above those of the quadrants into one X/Y offset, - # but don't keep Z transformations, as those come from the XTC stream - root = geometry.get_top_geo() - root_basis = basis_from_geo(root, use_z=False) - while len(root.get_list_of_children()) != 4 and len(root.get_list_of_children()) != 32: - assert len(root.get_list_of_children()) == 1 - root = root.get_list_of_children()[0] - root_basis *= basis_from_geo(root, use_z=False) - - metro[(0,)] = root_basis - - def add_sensor(quad_id, sensor_id, sensor): - metro[(0,quad_id,sensor_id)] = basis_from_geo(sensor) - - x, y, z = sensor.get_pixel_coords() - x/=1000; y/=1000; z/=1000 - assert x.shape == y.shape == z.shape - sensor_px_slow = x.shape[0] - sensor_px_fast = x.shape[1] - assert sensor_px_fast % 2 == 0 - - a0ul = sul = matrix.col((x[0,0],y[0,0],z[0,0])) - a1ur = sur = matrix.col((x[0,sensor_px_fast-1],y[0,sensor_px_fast-1],z[0,sensor_px_fast-1])) - a1lr = slr = matrix.col((x[sensor_px_slow-1,sensor_px_fast-1],y[sensor_px_slow-1,sensor_px_fast-1],z[sensor_px_slow-1,sensor_px_fast-1])) - a0ll = sll = matrix.col((x[sensor_px_slow-1,0],y[sensor_px_slow-1,0],z[sensor_px_slow-1,0])) - - a0ur = matrix.col((x[0,sensor_px_fast//2-1],y[0,sensor_px_fast//2-1],z[0,sensor_px_fast//2-1])) - a0lr = matrix.col((x[sensor_px_slow-1,sensor_px_fast//2-1],y[sensor_px_slow-1,sensor_px_fast//2-1],z[sensor_px_slow-1,sensor_px_fast//2-1])) - - a1ul = matrix.col((x[0,sensor_px_fast//2],y[0,sensor_px_fast//2],z[0,sensor_px_fast//2])) - a1ll = matrix.col((x[sensor_px_slow-1,sensor_px_fast//2],y[sensor_px_slow-1,sensor_px_fast//2],z[sensor_px_slow-1,sensor_px_fast//2])) - - sensor_center = center([sul,sur,slr,sll]) - asic0_center = center([a0ul,a0ur,a0lr,a0ll]) - asic1_center = center([a1ul,a1ur,a1lr,a1ll]) - - asic_trans0 = (asic0_center-sensor_center).length() - asic_trans1 = (asic1_center-sensor_center).length() - - if include_asic_offset: - rotated_ori = matrix.col((1,0,0)).axis_and_angle_as_unit_quaternion(180.0, deg=True) - offset_fast = -pixel_size*((sensor_px_fast) / 4) # 4 because sensor_px_fast is for sensor - offset_slow = +pixel_size*((sensor_px_slow) / 2) # Sensor is divided into 2 only in fast direction - metro[(0,quad_id,sensor_id,0)] = basis(orientation=rotated_ori,translation=matrix.col((-asic_trans0,0,0))) - metro[(0,quad_id,sensor_id,1)] = basis(orientation=rotated_ori,translation=matrix.col((+asic_trans1,0,0))) - metro[(0,quad_id,sensor_id,0)].translation += matrix.col((offset_fast, offset_slow, 0)) - metro[(0,quad_id,sensor_id,1)].translation += matrix.col((offset_fast, offset_slow, 0)) - else: - metro[(0,quad_id,sensor_id,0)] = basis(orientation=null_ori,translation=matrix.col((-asic_trans0,0,0))) - metro[(0,quad_id,sensor_id,1)] = basis(orientation=null_ori,translation=matrix.col((+asic_trans1,0,0))) - - if len(root.get_list_of_children()) == 4: - for quad_id, quad in enumerate(root.get_list_of_children()): - metro[(0,quad_id)] = basis_from_geo(quad) - for sensor_id, sensor in enumerate(quad.get_list_of_children()): - add_sensor(quad_id, sensor_id, sensor) - elif len(root.get_list_of_children()) == 32: - for quad_id in range(4): - metro[(0,quad_id)] = basis(orientation = null_ori, translation = matrix.col((0,0,0))) - sensors = root.get_list_of_children() - for sensor_id in range(8): - add_sensor(quad_id, sensor_id, sensors[quad_id*4+sensor_id]) - else: - assert False - - return metro +from dxtbx.format.FormatCBFMultiTile import cbf_wrapper +from dxtbx.format.cbf_writer import add_frame_specific_cbf_tables +from serialtbx.detector import basis, center +from serialtbx.detector.cspad import cspad_saturated_value, cspad_min_trusted_value, read_slac_metrology, asic_dimension, asic_gap, pixel_size def get_calib_file_path(env, address, run): """ Findes the path to the SLAC metrology file stored in a psana environment @@ -729,78 +435,6 @@ def metro_phil_to_basis_dict(metro): return bd -def add_frame_specific_cbf_tables(cbf, wavelength, timestamp, trusted_ranges, diffrn_id = "DS1", is_xfel = True, gain = 1.0, flux = None): - """ Adds tables to cbf handle that won't already exsist if the cbf file is just a header - @ param wavelength Wavelength in angstroms - @ param timestamp String formatted timestamp for the image - @ param trusted_ranges Array of trusted range tuples (min, max), one for each element """ - - """Data items in the DIFFRN_RADIATION category describe - the radiation used for measuring diffraction intensities, - its collimation and monochromatization before the sample. - - Post-sample treatment of the beam is described by data - items in the DIFFRN_DETECTOR category.""" - if flux: - cbf.add_category("diffrn_radiation", ["diffrn_id","wavelength_id","probe","beam_flux"]) - cbf.add_row([diffrn_id,"WAVELENGTH1","x-ray","%f"%flux]) - else: - cbf.add_category("diffrn_radiation", ["diffrn_id","wavelength_id","probe"]) - cbf.add_row([diffrn_id,"WAVELENGTH1","x-ray"]) - - """ Data items in the DIFFRN_RADIATION_WAVELENGTH category describe - the wavelength of the radiation used in measuring the diffraction - intensities. Items may be looped to identify and assign weights - to distinct wavelength components from a polychromatic beam.""" - cbf.add_category("diffrn_radiation_wavelength", ["id","wavelength","wt"]) - cbf.add_row(["WAVELENGTH1",str(wavelength),"1.0"]) - - """Data items in the DIFFRN_MEASUREMENT category record details - about the device used to orient and/or position the crystal - during data measurement and the manner in which the - diffraction data were measured.""" - cbf.add_category("diffrn_measurement",["diffrn_id","id","number_of_axes","method","details"]) - cbf.add_row([diffrn_id, - "INJECTION" if is_xfel else "unknown","0", - "electrospray" if is_xfel else "unknown" - "crystals injected by electrospray" if is_xfel else "unknown"]) - - """ Data items in the DIFFRN_SCAN category describe the parameters of one - or more scans, relating axis positions to frames.""" - cbf.add_category("diffrn_scan",["id","frame_id_start","frame_id_end","frames"]) - cbf.add_row(["SCAN1","FRAME1","FRAME1","1"]) - - """Data items in the DIFFRN_SCAN_FRAME category describe - the relationships of particular frames to scans.""" - cbf.add_category("diffrn_scan_frame",["frame_id","frame_number","integration_time","scan_id","date"]) - cbf.add_row(["FRAME1","1","0.0","SCAN1",timestamp]) - - """ Data items in the ARRAY_INTENSITIES category record the - information required to recover the intensity data from - the set of data values stored in the ARRAY_DATA category.""" - # More detail here: http://www.iucr.org/__data/iucr/cifdic_html/2/cif_img.dic/Carray_intensities.html - array_names = [] - cbf.find_category(b"diffrn_data_frame") - while True: - try: - cbf.find_column(b"array_id") - array_names.append(cbf.get_value().decode()) - cbf.next_row() - except Exception as e: - assert "CBF_NOTFOUND" in str(e) - break - - if not isinstance(gain, list): - gain = [gain] * len(array_names) - - - cbf.add_category("array_intensities",["array_id","binary_id","linearity","gain","gain_esd","overload","underload","undefined_value"]) - for i, array_name in enumerate(array_names): - overload = trusted_ranges[i][1] + 1 - underload = trusted_ranges[i][0] - undefined = underload - 1 - cbf.add_row([array_name,str(i+1),"linear","%f"%gain[i],"0.0",str(overload),str(underload),str(undefined)]) - def add_tiles_to_cbf(cbf, tiles, verbose = False): """ Given a cbf handle, add the raw data and the necessary tables to support it diff --git a/xfel/cftbx/detector/metrology2phil.py b/xfel/cftbx/detector/metrology2phil.py index f1f2ffc364..110826d104 100644 --- a/xfel/cftbx/detector/metrology2phil.py +++ b/xfel/cftbx/detector/metrology2phil.py @@ -27,15 +27,15 @@ def metrology2phil(calib_dir, verbose): return sections2phil(sections, verbose) def sections2phil(sections, verbose): - from xfel.cftbx.detector.metrology import master_phil + from serialtbx.detector.legacy_metrology.metrology import master_phil # Properties of CSPad pixels (Philipp et al., 2007). The counters # are 14 bits wide, and the pixels are square with a side length of # 110 um. Because cspad_tbx depends on pyana, it may fail to import # in which case a hardcoded fallback is provided. try: - from xfel.cxi.cspad_ana.cspad_tbx import cspad_saturated_value as sv - from xfel.cxi.cspad_ana.cspad_tbx import pixel_size as ps + from serialtbx.detector.cspad import cspad_saturated_value as sv + from serialtbx.detector.cspad import pixel_size as ps saturated_value = sv pixel_size = ps * 1e-3 except ImportError: diff --git a/xfel/clustering/singleframe.py b/xfel/clustering/singleframe.py index d4eb462723..a1c3443eed 100644 --- a/xfel/clustering/singleframe.py +++ b/xfel/clustering/singleframe.py @@ -319,7 +319,7 @@ def make_g6(uc): class SingleDialsFrame(SingleFrame): def __init__(self, refl=None, expt=None, id=None, **kwargs): - from xfel.command_line.frame_extractor import ConstructFrame + from serialtbx.util.construct_frame import ConstructFrame frame = ConstructFrame(refl, expt).make_frame() SingleFrame.__init__(self, dicti=frame, path=str(id), **kwargs) self.experiment = expt diff --git a/xfel/command_line/cspad_cbf_metrology.py b/xfel/command_line/cspad_cbf_metrology.py index ee29a6d2a7..3aa29617fa 100644 --- a/xfel/command_line/cspad_cbf_metrology.py +++ b/xfel/command_line/cspad_cbf_metrology.py @@ -496,7 +496,7 @@ def refine_expanding(params, merged_scope, combine_phil): # in this step,and apply that shift to the unrefined panels in this step if params.flat_refinement and params.flat_refinement_with_distance and i > 0: from dxtbx.model.experiment_list import ExperimentListFactory - from xfel.command_line.cspad_detector_congruence import iterate_detector_at_level, iterate_panels + from serialtbx.detector import iterate_detector_at_level, iterate_panels from scitbx.array_family import flex from scitbx.matrix import col from libtbx.test_utils import approx_equal diff --git a/xfel/command_line/cspad_circular_gain_mask.py b/xfel/command_line/cspad_circular_gain_mask.py index 63fa67938d..dfd5d64e98 100644 --- a/xfel/command_line/cspad_circular_gain_mask.py +++ b/xfel/command_line/cspad_circular_gain_mask.py @@ -69,7 +69,8 @@ print("Generating circular gain mask %s angstroms, assuming a distance %s mm and wavelength %s angstroms" % \ (params.resolution, params.distance, params.wavelength)) - from xfel.cftbx.detector.cspad_cbf_tbx import read_slac_metrology, get_cspad_cbf_handle + from serialtbx.detector.cspad import read_slac_metrology + from xfel.cftbx.detector.cspad_cbf_tbx import get_cspad_cbf_handle from dxtbx.format.FormatCBFCspad import FormatCBFCspadInMemory from dxtbx.model import BeamFactory metro = read_slac_metrology(path = params.optical_metrology_path) diff --git a/xfel/command_line/cspad_detector_congruence.py b/xfel/command_line/cspad_detector_congruence.py index 5196db4354..2e72fcca88 100644 --- a/xfel/command_line/cspad_detector_congruence.py +++ b/xfel/command_line/cspad_detector_congruence.py @@ -22,6 +22,7 @@ import libtbx.load_env import math from six.moves import zip +from serialtbx.detector import iterate_detector_at_level, iterate_panels, id_from_name, get_center try: from matplotlib import pyplot as plt @@ -62,66 +63,6 @@ for visualizing tilt. ''') -def iterate_detector_at_level(item, depth = 0, level = 0): - """ - Iterate through all panel groups or panels of a detector object at a given - hierarchy level - @param item panel group or panel. Use detector.hierarchy(). - @param depth current dept for recursion. Should be 0 for initial call. - @param level iterate groups at this level - @return next panel or panel group object - """ - if level == depth: - yield item - else: - for child in item: - for subitem in iterate_detector_at_level(child, depth+1, level): - yield subitem - -def iterate_panels(panelgroup): - """ - Find and iterate all panels in the given panel group, regardless of the hierarchly level - of this panelgroup - @param panelgroup the panel group of interest - @return the next panel - """ - if panelgroup.is_group(): - for child in panelgroup: - for subitem in iterate_panels(child): - yield subitem - else: - yield panelgroup - -def id_from_name(detector, name): - """ Jiffy function to get the id of a panel using its name - @param detector detector object - @param name panel name - @return index of panel in detector - """ - return [p.get_name() for p in detector].index(name) - -def get_center(pg): - """ Find the center of a panel group pg, projected on its fast/slow plane """ - if pg.is_group(): - # find the average center of all this group's children - children_center = col((0,0,0)) - count = 0 - for p in iterate_panels(pg): - children_center += get_center(p) - count += 1 - children_center /= count - - # project the children center onto the plane of the panel group - pgf = col(pg.get_fast_axis()) - pgs = col(pg.get_slow_axis()) - pgn = col(pg.get_normal()) - pgo = col(pg.get_origin()) - - return (pgf.dot(children_center) * pgf) + (pgs.dot(children_center) * pgs) + (pgn.dot(pgo) * pgn) - else: - s = pg.get_image_size() - return col(pg.get_pixel_lab_coord((s[0]/2, s[1]/2))) - def get_bounds(root, pg): """ Find the max extent of the panel group pg, projected onto the fast/slow plane of root """ def panel_bounds(root, panel): diff --git a/xfel/command_line/cspad_detector_shifts.py b/xfel/command_line/cspad_detector_shifts.py index c0f462265e..501e3df467 100644 --- a/xfel/command_line/cspad_detector_shifts.py +++ b/xfel/command_line/cspad_detector_shifts.py @@ -19,7 +19,7 @@ from scitbx.matrix import col from libtbx.phil import parse from libtbx.utils import Sorry -from xfel.command_line.cspad_detector_congruence import get_center +from serialtbx.detector import get_center import libtbx.load_env from six.moves import zip @@ -38,7 +38,7 @@ .help = Maximum hierarchy level to compute shifts to ''', process_includes=True) -from xfel.command_line.cspad_detector_congruence import iterate_detector_at_level, iterate_panels, id_from_name +from serialtbx.detector import iterate_detector_at_level, iterate_panels, id_from_name from xfel.command_line.cspad_detector_congruence import Script as ParentScript class Script(ParentScript): @@ -115,7 +115,7 @@ def run(self): table_header = ["PanelG","BC dist","Delta XY","R Offsets","T Offsets","Z Offsets","dR Norm","dT Norm","Local dNorm","Rot Z","N Refls"] table_header2 = ["ID","(mm)","(microns)","(microns)","(microns)","(microns)","(deg)","(deg)","(deg)","(deg)",""] - from xfel.cftbx.detector.cspad_cbf_tbx import basis + from serialtbx.detector import basis def get_full_basis_shift(pg): """Compute basis shift from pg to lab space""" shift = basis(panelgroup=pg) diff --git a/xfel/command_line/cxi_calibdir2cbfheader.py b/xfel/command_line/cxi_calibdir2cbfheader.py index 2781d067a8..362eaa715f 100644 --- a/xfel/command_line/cxi_calibdir2cbfheader.py +++ b/xfel/command_line/cxi_calibdir2cbfheader.py @@ -92,7 +92,8 @@ def get_asics_center(asics): # In order to match the image pickle metrology, we have to discard the tilts and offets in the # section objects, and instead use the active areas returned by corners_asic as the tile # locations. - from xfel.cftbx.detector.cspad_cbf_tbx import basis, pixel_size, asic_dimension, asic_gap + from serialtbx.detector import basis + from xfel.cftbx.detector.cspad_cbf_tbx import pixel_size, asic_dimension, asic_gap null_ori = col((0,0,1)).axis_and_angle_as_unit_quaternion(0, deg=True) # the detector is rotated 90 degrees from how corners_asic reports things diff --git a/xfel/command_line/cxi_display_metrology.py b/xfel/command_line/cxi_display_metrology.py index 1b1388cbdd..a30a7908c6 100644 --- a/xfel/command_line/cxi_display_metrology.py +++ b/xfel/command_line/cxi_display_metrology.py @@ -99,7 +99,7 @@ ax.set_ylim((0, 2000)) ax.set_ylim(ax.get_ylim()[::-1]) else: - from xfel.cftbx.detector.cspad_cbf_tbx import basis_from_geo + from serialtbx.detector.xtc import basis_from_geo root = geometry.get_top_geo() root_basis = basis_from_geo(root) diff --git a/xfel/command_line/cxi_slaccalib2cbfheader.py b/xfel/command_line/cxi_slaccalib2cbfheader.py index 299bb3fbbd..279f37e303 100644 --- a/xfel/command_line/cxi_slaccalib2cbfheader.py +++ b/xfel/command_line/cxi_slaccalib2cbfheader.py @@ -8,7 +8,8 @@ import sys, os import libtbx.phil from libtbx.utils import Usage, Sorry -from xfel.cftbx.detector.cspad_cbf_tbx import read_slac_metrology, write_cspad_cbf +from serialtbx.detector.cspad import read_slac_metrology +from xfel.cftbx.detector.cspad_cbf_tbx import write_cspad_cbf master_phil = libtbx.phil.parse(""" metrology_file = None diff --git a/xfel/command_line/detector_residuals.py b/xfel/command_line/detector_residuals.py index 8ca4533c46..e192ccdc5f 100644 --- a/xfel/command_line/detector_residuals.py +++ b/xfel/command_line/detector_residuals.py @@ -403,7 +403,8 @@ def trumpet_plot(experiment, reflections, axis = None): axis.plot(two_thetas, tan_outer_deg, "g.") axis.plot(two_thetas, -tan_outer_deg, "g.") -from xfel.command_line.cspad_detector_congruence import iterate_detector_at_level, iterate_panels, id_from_name, get_center, detector_plot_dict +from xfel.command_line.cspad_detector_congruence import detector_plot_dict +from serialtbx.detector import iterate_detector_at_level, iterate_panels, id_from_name, get_center from xfel.command_line.cspad_detector_congruence import Script as DCScript class Script(DCScript): ''' Class to parse the command line options. ''' diff --git a/xfel/command_line/frame_extractor.py b/xfel/command_line/frame_extractor.py index 17dd0974d2..77ebea0164 100644 --- a/xfel/command_line/frame_extractor.py +++ b/xfel/command_line/frame_extractor.py @@ -7,14 +7,12 @@ # # $Id: frame_extractor.py idyoung $ -from dials.array_family import flex from dials.util.options import Importer, flatten_reflections, flatten_experiments, ArgumentParser -from cctbx import crystal, miller -from cctbx.crystal_orientation import crystal_orientation import iotbx.phil import cctbx, os, glob from libtbx import easy_pickle from six.moves import zip +from serialtbx.util.construct_frame import ConstructFrame phil_scope = iotbx.phil.parse(""" input { @@ -35,191 +33,6 @@ } """) -class ConstructFrame(object): - def get_template_pickle(self): - return { - 'beam_s0': 0, - 'beam_polarization_normal': 0, - 'current_cb_op_to_primitive': 0, - 'current_orientation':0, - 'distance':0, - 'effective_tiling':0, - 'identified_isoform':None, - 'mapped_predictions':[[]], - 'max_signal':0, - 'ML_domain_size_ang':[0], - 'ML_half_mosaicity_deg':[0], - 'mosaicity':0, - 'model_partialities':[None], - 'observations':[0], - 'pointgroup':0, - 'residual':0, - 'sa_parameters':['None'], - 's1_vec':[0], - 'wavelength':0, - 'xbeam':0, - 'ybeam':0} - - def __init__(self, reflections, experiment): - # assemble template and unpack files - self.frame = self.get_template_pickle() - self.pixel_size = experiment.detector[0].get_pixel_size()[0] - - if 'intensity.prf.value' in reflections: - self.method = 'prf' # integration by profile fitting - elif 'intensity.sum.value' in reflections: - self.method = 'sum' # integration by simple summation - self.reflections = reflections.select(reflections['intensity.' + self.method + '.variance'] > 0) # keep only spots with sigmas above zero - - self.xtal = experiment.crystal - self.beam_obj = experiment.beam - self.det = experiment.detector - self.gonio = experiment.goniometer - self.scan = experiment.scan - self.img_sweep = experiment.imageset - - # experiment-dependent components --------------------------------------------------------------------------- - - # get wavelength - def populate_wavelength(self): - assert self.beam_obj.get_wavelength() is not None, "no wavelength" - self.frame['wavelength'] = self.beam_obj.get_wavelength() - self.frame['beam_s0'] = self.beam_obj.get_s0() - self.frame['beam_polarization_normal'] = self.beam_obj.get_polarization_normal() - - # get detector distance in mm - def populate_distance(self): - assert self.det[0].get_distance() is not None, "no detector distance" - self.frame['distance'] = self.det[0].get_distance() - - # get xbeam and ybeam in mm - def populate_beam_dir(self): - assert self.beam_obj.get_s0() is not None, "no beam direction" - self.frame['xbeam'], self.frame['ybeam'] = self.det[0].get_beam_centre(self.beam_obj.get_s0()) - - # get isoform - def populate_isoform(self): - if hasattr(self.xtal, "identified_isoform"): - self.frame['identified_isoform'] = self.xtal.identified_isoform - - # get max signal - def populate_max_signal(self): - pass - - # get effective tiling - def populate_effective_tiling(self): - pass - - # indicate simulated annealing parameters, if present - def populate_sa_params(self): - pass - - # crystal-dependent components ------------------------------------------------------------------------------ - - # generate a crystal orientation object from the A* matrix - def populate_orientation(self): - assert self.xtal.get_A() is not None, "no crystal orientation matrix" - self.frame['current_orientation'] = [crystal_orientation(self.xtal.get_A(), True)] - - # generate change-of-basis operation for current to primitive cell - def populate_op_to_primitive(self): - assert self.xtal.get_space_group() is not None, "no space group" - self.frame['current_cb_op_to_primitive'] = [self.xtal.get_space_group().z2p_op()] - - # fetch the point group associated with the crystal - def populate_point_group(self): - assert self.xtal.get_space_group() is not None, "no space group" - self.frame['pointgroup'] = str(self.xtal.get_space_group().build_derived_point_group().info()) - - # get mosaicity - def populate_mosaicity(self): - try: - assert self.xtal.get_mosaicity() is not None, "no mosaicity" - self.frame['mosaicity'] = self.xtal.get_mosaicity() - except Exception: - pass - - # get any available ML values - def populate_ML_values(self): - try: - self.frame['ML_half_mosaicity_deg'] = [self.xtal.get_half_mosaicity_deg()] - except AttributeError: - pass - try: - self.frame['ML_domain_size_ang'] = [self.xtal.get_domain_size_ang()] - except AttributeError: - pass - - # observations-dependent components ------------------------------------------------------------------------- - - # generate a miller array containing the Miller indices, intensities and variances for one frame - def populate_observations(self): - intensities = self.reflections['intensity.' + self.method + '.value'] - variances = self.reflections['intensity.' + self.method + '.variance'] - space_group = crystal.symmetry(self.xtal.get_unit_cell(), str(self.xtal.get_space_group().info())) - miller_set = miller.set(space_group, self.reflections['miller_index']) - self.frame['observations'][0] = cctbx.miller.array(miller_set, intensities, flex.sqrt(variances)).set_observation_type_xray_intensity() - self.frame['s1_vec'][0] = self.reflections['s1'] - - # collect predicted spot positions - def populate_pixel_positions(self): - assert 'xyzcal.px' in self.reflections, "no calculated spot positions" - self.frame['mapped_predictions'][0] = flex.vec2_double() - for i in range(len(self.reflections['xyzcal.px'])): - self.frame['mapped_predictions'][0].append(tuple(self.reflections['xyzcal.px'][i][1::-1])) # 1::-1 reverses the order taking only the first two elements first. - - # generate a list of dictionaries containing a series of corrections for each predicted reflection - def populate_corrections(self): - assert 'xyzobs.px.value' in self.reflections and 'xyzcal.px' in self.reflections, "no calculated or observed spot positions" - assert self.frame['xbeam'] != 0 and self.frame['ybeam'] != 0, "invalid beam center" - self.frame['correction_vectors'] = [[]] - for idx in range(len(self.reflections['xyzobs.px.value'])): - if self.reflections['xyzcal.px'][idx][0:2] != self.reflections['xyzobs.px.value'][idx][0:2]: - theoret_center = 1765/2, 1765/2 - refined_center = self.frame['xbeam']/self.pixel_size, self.frame['ybeam']/self.pixel_size # px to mm conversion - hkl = self.reflections['miller_index'][idx] - obsspot = tuple(self.reflections['xyzobs.px.value'][idx][0:2]) - predspot = tuple(self.reflections['xyzcal.px'][idx][0:2]) - self.frame['correction_vectors'][0].append({'refinedcenter':refined_center, 'hkl':hkl, 'setting_id':0, 'azimuthal':0, 'radial':0, - 'obsspot':obsspot, 'obscenter':theoret_center, 'predspot':predspot}) - - # get partialities - def populate_partialities(self): - pass - - # produce residuals - def populate_residuals(self): - pass - - # collect kapton corrections if present - def populate_kapton_corrections(self): - if 'kapton_absorption_correction' in self.reflections: - self.frame['fuller_kapton_absorption_correction'] = [self.reflections['kapton_absorption_correction']] - if 'kapton_absorption_correction_sigmas' in self.reflections: - self.frame['fuller_kapton_absorption_correction_sigmas'] = [self.reflections['kapton_absorption_correction_sigmas']] - - # combine all of the above - def make_frame(self): - self.populate_wavelength() - self.populate_distance() - self.populate_beam_dir() - self.populate_isoform() - self.populate_max_signal() - self.populate_effective_tiling() - self.populate_sa_params() - self.populate_orientation() - self.populate_op_to_primitive() - self.populate_point_group() - self.populate_mosaicity() - self.populate_ML_values() - self.populate_observations() - self.populate_pixel_positions() - # self.populate_corrections() # works, but unnecessary - self.populate_partialities() - self.populate_residuals() - self.populate_kapton_corrections() - return self.frame - class ConstructFrameFromFiles(ConstructFrame): def __init__(self, refl_name, json_name, outname=None): # load the integration.refl file (reflection table) into memory and diff --git a/xfel/command_line/quadrants_cbf.py b/xfel/command_line/quadrants_cbf.py index 127dca2fff..f585055777 100644 --- a/xfel/command_line/quadrants_cbf.py +++ b/xfel/command_line/quadrants_cbf.py @@ -5,7 +5,7 @@ import sys,os from scitbx.matrix import col -from xfel.cftbx.detector.cspad_cbf_tbx import center +from serialtbx.detector import center import dxtbx import libtbx.load_env from libtbx.utils import Usage diff --git a/xfel/command_line/radial_average.py b/xfel/command_line/radial_average.py index 03c53d39db..6ed8632505 100644 --- a/xfel/command_line/radial_average.py +++ b/xfel/command_line/radial_average.py @@ -59,7 +59,7 @@ h_y = 3 def run (args, source_data = None) : - from xfel import radial_average + from dxtbx.ext import radial_average from scitbx.array_family import flex from iotbx.detectors.cspad_detector_formats import reverse_timestamp from iotbx.detectors.cspad_detector_formats import detector_format_version as detector_format_function diff --git a/xfel/cxi/cspad_ana/common_mode.py b/xfel/cxi/cspad_ana/common_mode.py index d825fe322b..cc15802c5d 100644 --- a/xfel/cxi/cspad_ana/common_mode.py +++ b/xfel/cxi/cspad_ana/common_mode.py @@ -21,6 +21,7 @@ from xfel.cxi.cspad_ana import cspad_tbx from xfel.cxi.cspad_ana import skip_event_flag from xfel.cxi.cspad_ana.mod_event_info import mod_event_info +from serialtbx.detector.xtc import old_address_to_new_address class common_mode_correction(mod_event_info): @@ -380,7 +381,7 @@ def event(self, evt, env): self.address == 'MfxEndstation-0|Rayonix-0': from psana import Source, Camera import numpy as np - address = cspad_tbx.old_address_to_new_address(self.address) + address = old_address_to_new_address(self.address) src=Source('DetInfo(%s)'%address) self.cspad_img = evt.get(Camera.FrameV1,src) if self.cspad_img is not None: diff --git a/xfel/cxi/cspad_ana/cspad_tbx.py b/xfel/cxi/cspad_ana/cspad_tbx.py index a9ebab1ae8..f42232c767 100644 --- a/xfel/cxi/cspad_ana/cspad_tbx.py +++ b/xfel/cxi/cspad_ana/cspad_tbx.py @@ -23,23 +23,13 @@ import six from six.moves import zip -__version__ = "$Revision$" - +import serialtbx.util.time +import serialtbx.detector.cspad +from serialtbx.detector.cspad import pixel_size +from serialtbx.detector.xtc import old_address_to_new_address, get_ebeam, env_detz, address_split # import dependency -# The CAMP and CSpad counters are both 14 bits wide (Strüder et al -# 2010; Philipp et al., 2007), which means the physical limit is 2**14 - 1. -# However, in practice, when the pixels are in the low gain mode, after -# correcting by a gain value of around 6.87, the pixels tend to saturate -# around 90000. See xpp experiment xppe0314, run 184 as evidence. -cspad_saturated_value = 90000 +__version__ = "$Revision$" -# The dark average for the CSPAD detector is around 1100-1500. A pixel -# histogram of a minimum projection of an uncorrected (raw) light run shows -# a mostly flat tail up to ~800 ADU with a few bumps in the tail which -# represent true underloads. Assume a dark average of 1200 ADU. After dark -# subtraction, 800 - 1200 gives a minimum trusted value of -400. Reject -# pixels less than this. -cspad_min_trusted_value = -400 # As long as the mask value is outside of the trusted range, the pixel should # be ignored by any downstream software. @@ -49,14 +39,6 @@ # XXX This should be obsoleted! npix_quad = 850 -# The pixel size in mm. The pixel size is fixed and square, with side -# length of 110 µm (Philipp et al., 2007). XXX Should really clarify -# this with Sol and Chris. -# -# XXX Andor: 13.5 µm square, CAMP: 75 µm, square (Strüder et al., -# 2010) -pixel_size = 110e-3 - # origin of section in quad coordinate system. x-position # correspond to column number. XXX Note/reference the source! # XXX This should be obsoleted! @@ -71,51 +53,6 @@ [ 0, 0, 214, 1, 425, 425, 615, 403]] # 2:5 were not measured -def address_split(address, env=None): - """The address_split() function splits an address into its four - components. Address strings are on the form - detector-detectorID|device-deviceID, where the detectors must be in - dir(xtc.DetInfo.Detector) and device must be in - (xtc.DetInfo.Device). - @param address Full data source address of the DAQ device - @param env Optional env to dereference an alias into an address - @return Four-tuple of detector name, detector ID, device, and - device ID - """ - - import re - - # pyana - m = re.match( - r"^(?P\S+)\-(?P\d+)\|(?P\S+)\-(?P\d+)$", address) - if m is not None: - return (m.group('det'), m.group('det_id'), m.group('dev'), m.group('dev_id')) - - # psana - m = re.match( - r"^(?P\S+)\.(?P\d+)\:(?P\S+)\.(?P\d+)$", address) - if m is not None: - return (m.group('det'), m.group('det_id'), m.group('dev'), m.group('dev_id')) - - # psana DetInfo string - m = re.match( - r"^DetInfo\((?P\S+)\.(?P\d+)\:(?P\S+)\.(?P\d+)\)$", address) - if m is not None: - return (m.group('det'), m.group('det_id'), m.group('dev'), m.group('dev_id')) - - if env is not None: - # Try to see if this is a detector alias, and if so, dereference it. Code from psana's Detector/PyDetector.py - amap = env.aliasMap() - alias_src = amap.src(address) # string --> DAQ-style psana.Src - - # if it is an alias, look up the full name - if amap.alias(alias_src) != '': # alias found - address = str(alias_src) - return address_split(address) - - return (None, None, None, None) - - def cbcaa(config, sections): """The cbcaa() function uses on-disk calibration data to estimate the beam centre and the active detector areas. The beam centre is @@ -448,82 +385,9 @@ def CsPadElement(data3d, qn, config): quadrant = numpy.rot90(quadrant, 4 - qn) return quadrant -def dpack(active_areas=None, - address=None, - beam_center_x=None, - beam_center_y=None, - ccd_image_saturation=None, - data=None, - distance=None, - pixel_size=pixel_size, - saturated_value=None, - timestamp=None, - wavelength=None, - xtal_target=None, - min_trusted_value=None): - """XXX Check completeness. Should fill in sensible defaults.""" - - # Must have data. - if data is None: - return None - - # Create a time stamp of the current time if none was supplied. - if timestamp is None: - timestamp = evt_timestamp() - - # For unknown historical reasons, the dictionary must contain both - # CCD_IMAGE_SATURATION and SATURATED_VALUE items. - if ccd_image_saturation is None: - if saturated_value is None: - ccd_image_saturation = cspad_saturated_value - else: - ccd_image_saturation = saturated_value - if saturated_value is None: - saturated_value = ccd_image_saturation - - # Use a minimum value if provided for the pixel range - if min_trusted_value is None: - min_trusted_value = cspad_min_trusted_value - - # By default, the beam center is the center of the image. The slow - # (vertical) and fast (horizontal) axes correspond to x and y, - # respectively. - if beam_center_x is None: - beam_center_x = pixel_size * data.focus()[1] / 2 - if beam_center_y is None: - beam_center_y = pixel_size * data.focus()[0] / 2 - - # By default, the entire detector image is an active area. There is - # no sensible default for distance nor wavelength. XXX But setting - # wavelength to zero may be disastrous? - if active_areas is None: - # XXX Verify order with non-square detector - active_areas = flex.int((0, 0, data.focus()[0], data.focus()[1])) - if distance is None: - distance = 0 - if wavelength is None: - wavelength = 0 - - # The size must match the image dimensions. The length along the - # slow (vertical) axis is SIZE1, the length along the fast - # (horizontal) axis is SIZE2. - return {'ACTIVE_AREAS': active_areas, - 'BEAM_CENTER_X': beam_center_x, - 'BEAM_CENTER_Y': beam_center_y, - 'CCD_IMAGE_SATURATION': ccd_image_saturation, - 'DATA': data, - 'DETECTOR_ADDRESS': address, - 'DISTANCE': distance, - 'PIXEL_SIZE': pixel_size, - 'SATURATED_VALUE': saturated_value, - 'MIN_TRUSTED_VALUE': min_trusted_value, - 'SIZE1': data.focus()[0], - 'SIZE2': data.focus()[1], - 'TIMESTAMP': timestamp, - 'SEQUENCE_NUMBER': 0, # XXX Deprecated - 'WAVELENGTH': wavelength, - 'xtal_target': xtal_target} - +def dpack(*kwargs): + """ thin wrapper """ + return serialtbx.detector.cspad.dpack(*kwargs) def hdf5pack(hdf5_file, active_areas=None, @@ -827,31 +691,6 @@ def get_value(self, key, args, kwargs_local): subprocess=env.subprocess(), user=getuser()) -def get_ebeam(evt): - try: - # pyana - ebeam = evt.getEBeam() - except AttributeError as e: - from psana import Source, Bld - src = Source('BldInfo(EBeam)') - ebeam = evt.get(Bld.BldDataEBeamV6, src) - if ebeam is None: - ebeam = evt.get(Bld.BldDataEBeamV5, src) - if ebeam is None: - ebeam = evt.get(Bld.BldDataEBeamV4, src) - if ebeam is None: - ebeam = evt.get(Bld.BldDataEBeamV3, src) - if ebeam is None: - ebeam = evt.get(Bld.BldDataEBeamV2, src) - if ebeam is None: - ebeam = evt.get(Bld.BldDataEBeamV1, src) - if ebeam is None: - ebeam = evt.get(Bld.BldDataEBeamV0, src) - if ebeam is None: - ebeam = evt.get(Bld.BldDataEBeam, src) # recent version of psana will return a V7 event or higher if this type is asked for - - return ebeam - def env_laser_status(env, laser_id): """The return value is a bool that indicates whether the laser in question was on for that particular shot. Bear in mind that sample @@ -895,83 +734,9 @@ def env_injector_xyz(env): for i in range(3)]) -def env_detz(address, env): - """The env_detz() function returns the position of the detector with - the given address string on the z-axis in mm. The zero-point is as - far away as possible from the sample, and values decrease as the - detector is moved towards the sample. - @param address Full data source address of the DAQ device - @param env Environment object - @return Detector z-position, in mm - """ - - if env is not None: - detector = address_split(address, env)[0] - if detector is None: - return None - elif detector == 'CxiDs1': - pv = env.epicsStore().value('CXI:DS1:MMS:06.RBV') - if pv is None: - # Even though potentially unsafe, fall back on the commanded - # value if the corresponding read-back value cannot be read. - # According to Sébastien Boutet, this particular motor has not - # caused any problem in the past. - pv = env.epicsStore().value('CXI:DS1:MMS:06') - if pv is None: - # Try the other detector. These are sometimes inconsistent - pv = env.epicsStore().value('CXI:DS2:MMS:06.RBV') - elif detector == 'CxiDsd' or detector == 'CxiDs2': - # XXX Note inconsistency in naming: Dsd vs Ds2! - pv = env.epicsStore().value('CXI:DS2:MMS:06.RBV') - if pv is None: - # Try the other detector. These are sometimes inconsistent - pv = env.epicsStore().value('CXI:DS1:MMS:06.RBV') - elif detector == 'XppGon': - # There is no distance recorded for the XPP's CSPAD on the robot - # arm. Always return zero to allow the distance to be set using - # the offset. - return 0 - elif detector == 'XppEndstation' or \ - detector == 'MfxEndstation': - # There is no distance recorded for the XPP's or MFX's Rayonix - # on the robot arm. Always return zero to allow the distance to - # be set using the offset. - return 0 - else: - return None - - if pv is None: - return None - - if hasattr(pv, "values"): - if len(pv.values) == 1: - return pv.values[0] - else: - return None - return pv - - return None - - -def env_distance(address, env, offset): - """The env_distance() function returns the distance between the - sample and the detector with the given address string in mm. The - distance between the sample and the the detector's zero-point can - vary by an inch or more between different LCLS runs. According to - Sébastien Boutet the offset should be stable to within ±0.5 mm - during a normal experiment. - - @param address Full data source address of the DAQ device - @param env Environment object - @param offset Detector-sample offset in mm, corresponding to - longest detector-sample distance - @return Detector-sample distance, in mm - """ - - detz = env_detz(address, env) - if detz is not None: - return detz + offset - return None +def env_distance(*kwargs): + """ thin wrapper """ + return serialtbx.detector.xtc.env_distance(*kwargs) def env_sifoil(env): @@ -1231,9 +996,7 @@ def evt_time(evt=None): """ if evt is None: - t = time.time() - s = int(math.floor(t)) - return (s, int(round((t - s) * 1000))) + return serialtbx.util.time.now_s_ms() if hasattr(evt, "getTime"): t = evt.getTime() @@ -1255,57 +1018,11 @@ def evt_timestamp(t=None): @return Human-readable ISO 8601 timestamp in string representation """ - if t is None: - t = evt_time(evt=None) - if t is None: - return None - return time.strftime("%Y-%m-%dT%H:%MZ%S", time.gmtime(t[0])) + \ - (".%03d" % t[1]) - - -def evt_wavelength(evt, delta_k=0): - """The evt_wavelength() function returns the wavelength in Ångström - of the event pointed to by @p evt. From Margaritondo & Rebernik - Ribic (2011): the dimensionless relativistic γ-factor is derived - from beam energy in MeV and the electron rest mass, K is a - dimensionless "undulator parameter", and L is the macroscopic - undulator period in Ångström. See also - https://people.eecs.berkeley.edu/~attwood/srms/2007/Lec10.pdf - - @param evt Event data object, a configure object - @param delta_k Optional K-value correction - @return Wavelength, in Ångström - """ - - if evt is not None: - ebeam = get_ebeam(evt) + return serialtbx.util.time.timestamp(t) - if hasattr(ebeam, 'fEbeamPhotonEnergy') and ebeam.fEbeamPhotonEnergy > 0: - # pyana - return 12398.4187 / ebeam.fEbeamPhotonEnergy - if hasattr(ebeam, 'ebeamPhotonEnergy') and ebeam.ebeamPhotonEnergy() > 0: - # psana - return 12398.4187 / ebeam.ebeamPhotonEnergy() - - if hasattr(ebeam, 'fEbeamL3Energy') and ebeam.fEbeamL3Energy > 0: - # pyana - gamma = ebeam.fEbeamL3Energy / 0.510998910 - elif hasattr(ebeam, 'ebeamL3Energy') and ebeam.ebeamL3Energy() > 0: - # psana - gamma = ebeam.ebeamL3Energy() / 0.510998910 - else: - return None - K = 3.5 + delta_k - L = 3.0e8 - return L / (2 * gamma**2) * (1 + K**2 / 2) - return None - -def old_address_to_new_address(address): - """ Change between old and new style detector addresses. - I.E. CxiDs1-0|Cspad-0 becomes CxiDs1.0:Cspad.0 - @param address detector address to convert - """ - return address.replace('-','.').replace('|',':') +def evt_wavelength(*kwargs): + """ thin wrapper """ + return serialtbx.detector.xtc.evt_wavelength(*kwargs) def getConfig(address, env): """ Given a detector address, find the config object in an env object diff --git a/xfel/cxi/cspad_ana/parse_calib.py b/xfel/cxi/cspad_ana/parse_calib.py index be2091dc68..761886f24f 100644 --- a/xfel/cxi/cspad_ana/parse_calib.py +++ b/xfel/cxi/cspad_ana/parse_calib.py @@ -256,7 +256,7 @@ def v2calib2sections(filename): @return Section objects """ - from xfel.cftbx.detector.cspad_cbf_tbx import read_slac_metrology + from serialtbx.detector.cspad import read_slac_metrology from scitbx.matrix import sqr from xfel.cxi.cspad_ana.cspad_tbx import pixel_size diff --git a/xfel/cxi/cspad_ana/rayonix_tbx.py b/xfel/cxi/cspad_ana/rayonix_tbx.py index fc7c46999e..9fe82670af 100644 --- a/xfel/cxi/cspad_ana/rayonix_tbx.py +++ b/xfel/cxi/cspad_ana/rayonix_tbx.py @@ -1,46 +1,10 @@ from __future__ import absolute_import, division, print_function # Utility functions for the Rayonix Detector. +from dxtbx.format.cbf_writer import add_frame_specific_cbf_tables from scitbx.array_family import flex -# given value of rayonix detector saturation xppi6115 -rayonix_saturated_value = 2**16 -1 - -# minimum value for rayonix data -rayonix_min_trusted_value = 0 -rayonix_max_trusted_value = rayonix_saturated_value - 1 - -def get_rayonix_pixel_size(bin_size): - ''' Given a bin size determine a pixel size. - -Michael Blum from Rayonix said The pixel size is recorded in the header, -but can be derived trivially from the overall dimension of the corrected imaging -area (170mm) and the number of pixels. (3840 unbinned). The corrected image is -forced to this size. - -unbinned 170/3840 = 0.04427 - -I believe the accuracy of the MEAN pixel size to be at least as good as 0.1% -which is the limit to which I can measure our calibration plate and exceeds the - parallax error in our calibration station. - -Note, the Rayonix MX340 has the same pixel size as the MX170: - -unbinned 340/7680 = 0.04427 - - @param bin_size rayonix bin size as an integer - ''' - pixel_size=bin_size*170/3840 - return pixel_size - -def get_rayonix_detector_dimensions(env): - ''' Given a psana env object, find the detector dimensions - @param env psana environment object - ''' - import psana - cfgs = env.configStore() - rayonix_cfg = cfgs.get(psana.Rayonix.ConfigV2, psana.Source('Rayonix')) - if not rayonix_cfg: return None, None - return rayonix_cfg.width(), rayonix_cfg.height() +from serialtbx.detector.rayonix import rayonix_min_trusted_value, rayonix_max_trusted_value, get_rayonix_pixel_size +from serialtbx.detector.rayonix import get_data_from_psana_event, rayonix_saturated_value, get_rayonix_detector_dimensions # import dependency def get_rayonix_cbf_handle(tiles, metro, timestamp, cbf_root, wavelength, distance, bin_size, detector_size, verbose = True, header_only = False): # set up the metrology dictionary to include axis names, pixel sizes, and so forth @@ -68,7 +32,7 @@ def get_rayonix_cbf_handle(tiles, metro, timestamp, cbf_root, wavelength, distan basis.axis_name = detector_axes_names[-1] # the data block is the root cbf node - from xfel.cftbx.detector.cspad_cbf_tbx import cbf_wrapper + from dxtbx.format.FormatCBFMultiTile import cbf_wrapper import os cbf=cbf_wrapper() cbf.new_datablock(os.path.splitext(os.path.basename(cbf_root))[0].encode()) @@ -231,7 +195,7 @@ class FormatCBFRayonixInMemory(FormatCBFFullStillInMemory): def get_dxtbx_from_params(params, detector_size): """ Build a dxtbx format object for the Rayonix based on input paramters (beam center and binning) """ - from xfel.cftbx.detector.cspad_cbf_tbx import basis + from serialtbx.detector import basis from scitbx.matrix import col fake_distance = 100 null_ori = col((0,0,1)).axis_and_angle_as_unit_quaternion(0, deg=True) @@ -250,21 +214,6 @@ def get_dxtbx_from_params(params, detector_size): base_dxtbx = FormatCBFRayonixInMemory(cbf) return base_dxtbx -def get_data_from_psana_event(evt, address): - """ Read the pixel data for a Rayonix image from an event - @param psana event object - @param address old style psana detector address - @return numpy array with raw data""" - from psana import Source, Camera - from xfel.cxi.cspad_ana import cspad_tbx - import numpy as np - address = cspad_tbx.old_address_to_new_address(address) - src=Source(address) - data = evt.get(Camera.FrameV1,src) - if data is not None: - data = data.data16().astype(np.float64) - return data - def format_object_from_data(base_dxtbx, data, distance, wavelength, timestamp, address, round_to_int=True): """ Given a preloaded dxtbx format object and raw data, assemble the tiles @@ -288,7 +237,7 @@ def format_object_from_data(base_dxtbx, data, distance, wavelength, timestamp, a data = flex.double(data.astype(np.float64)) n_asics = data.focus()[0] * data.focus()[1] - cspad_cbf_tbx.add_frame_specific_cbf_tables(cbf, wavelength,timestamp, + add_frame_specific_cbf_tables(cbf, wavelength,timestamp, [(rayonix_min_trusted_value, rayonix_max_trusted_value)]*n_asics) # Set the distance, I.E., the length translated along the Z axis diff --git a/xfel/ext.cpp b/xfel/ext.cpp index db658f9794..a19b3834ec 100644 --- a/xfel/ext.cpp +++ b/xfel/ext.cpp @@ -678,37 +678,6 @@ double distance_between_points(scitbx::vec2 const& a, scitbx::vec2 con return std::sqrt((std::pow(double(b[0]-a[0]),2)+std::pow(double(b[1]-a[1]),2))); } -void radial_average(scitbx::af::versa > & data, -scitbx::af::versa > & mask, -scitbx::vec2 const& beam_center, -scitbx::af::shared sums, -scitbx::af::shared sums_sq, -scitbx::af::shared counts, -double pixel_size, double distance, -scitbx::vec2 const& upper_left, -scitbx::vec2 const& lower_right) { - std::size_t extent = sums.size(); - double extent_in_mm = extent * pixel_size; - double extent_two_theta = std::atan(extent_in_mm/distance)*180/scitbx::constants::pi; - - for(std::size_t y = upper_left[1]; y < lower_right[1]; y++) { - for(std::size_t x = upper_left[0]; x < lower_right[0]; x++) { - double val = data(x,y); - if(val > 0 && mask(x,y)) { - scitbx::vec2 point((int)x,(int)y); - double d_in_mm = distance_between_points(point,beam_center) * pixel_size; - double twotheta = std::atan(d_in_mm/distance)*180/scitbx::constants::pi; - std::size_t bin = (std::size_t)std::floor(twotheta*extent/extent_two_theta); - if (bin >= extent) - continue; - sums[bin] += val; - sums_sq[bin] += val*val; - counts[bin]++; - } - } - } -} - void radial_histogram(scitbx::af::versa > & data, scitbx::af::versa > & histogram, scitbx::vec2 const& intensity_extent, @@ -849,12 +818,6 @@ namespace boost_python { namespace { .def("parse_from_line",&column_parser::parse_from_line) ; - def("radial_average", &radial_average, - (arg("data"), arg("beam_center"), arg("sums"), arg("sums_sq"), arg("counts"), - arg("pixel_size"), arg("distance"), - arg("upper_left"), arg("lower_right"))) - ; - def("radial_histogram", &radial_histogram, (arg("data"), arg("histogram"), arg("intensity_extent"), arg("beam_center"), arg("pixel_size"), arg("distance"), diff --git a/xfel/mono_simulation/max_like_pos_neg.py b/xfel/mono_simulation/max_like_pos_neg.py index 80163c4dab..c40578d402 100644 --- a/xfel/mono_simulation/max_like_pos_neg.py +++ b/xfel/mono_simulation/max_like_pos_neg.py @@ -169,7 +169,7 @@ def compute_functional_and_gradients(self): print("pos_neg output Deff/2, eta_deg ",1./M.x[0], M.x[1]*180./pi) - from xfel.mono_simulation.max_like import minimizer as legacy_minimizer + from serialtbx.mono_simulation.max_like import minimizer as legacy_minimizer Q = legacy_minimizer(d_i, psi_i, eta_rad, Deff) print("legacy output Deff/2, eta_deg ",1./Q.x[0], Q.x[1]*180./pi) diff --git a/xfel/mono_simulation/mono_treatment.py b/xfel/mono_simulation/mono_treatment.py index 67820c445f..4861a157de 100644 --- a/xfel/mono_simulation/mono_treatment.py +++ b/xfel/mono_simulation/mono_treatment.py @@ -344,7 +344,7 @@ def refine_rotx_roty2(OO,enable_rotational_target=True): tan_outer_deg = tan_phi_deg + k_degrees if OO.mosaic_refinement_target=="ML": - from xfel.mono_simulation.max_like import minimizer + from serialtbx.mono_simulation.max_like import minimizer print("input", s_ang,2. * solution[1]*180/math.pi) # coerce the estimates to be positive for max-likelihood lower_limit_domain_size = math.pow( @@ -474,7 +474,7 @@ def refine_rotx_roty2(OO,enable_rotational_target=True): OO.parent.show_figure(plt,fig,"psi") plt.close() - from xfel.mono_simulation.util import green_curve_area,ewald_proximal_volume + from serialtbx.mono_simulation.util import green_curve_area,ewald_proximal_volume if OO.mosaic_refinement_target=="ML": OO.parent.green_curve_area = green_curve_area(two_thetas, tan_outer_deg_ML) OO.parent.inputai.setMosaicity(M.x[1]*180./math.pi) # full width, degrees diff --git a/xfel/poly/recompute_mosaic_params.py b/xfel/poly/recompute_mosaic_params.py index 8e4ac48faa..de5762ef78 100644 --- a/xfel/poly/recompute_mosaic_params.py +++ b/xfel/poly/recompute_mosaic_params.py @@ -1,7 +1,7 @@ from __future__ import absolute_import, division, print_function from scitbx.matrix import col, sqr -from xfel.mono_simulation import max_like +from serialtbx.mono_simulation import max_like import numpy as np from scipy import constants ENERGY_CONV = 1e10*constants.c*constants.h / constants.electron_volt diff --git a/xfel/swissfel/jf16m_cxigeom2nexus.py b/xfel/swissfel/jf16m_cxigeom2nexus.py index 4ccb3f044c..7c30b5d269 100644 --- a/xfel/swissfel/jf16m_cxigeom2nexus.py +++ b/xfel/swissfel/jf16m_cxigeom2nexus.py @@ -8,7 +8,7 @@ import six from libtbx.utils import Sorry import datetime -from xfel.util.jungfrau import pad_stacked_format +from serialtbx.detector.jungfrau import pad_stacked_format phil_scope = parse(""" unassembled_file = None diff --git a/xfel/util/dials_pickle_reader.py b/xfel/util/dials_pickle_reader.py index 77abba1fa0..4087560b6b 100644 --- a/xfel/util/dials_pickle_reader.py +++ b/xfel/util/dials_pickle_reader.py @@ -6,7 +6,7 @@ import os from dials.util.options import Importer, flatten_reflections, flatten_experiments -from xfel.command_line.frame_extractor import ConstructFrame +from serialtbx.util.construct_frame import ConstructFrame def find_json(pickle, pickle_ext=None, json_ext=None): """find the matching json file for a given dials-format non-image pickle file""" diff --git a/xfel/util/sublattice_helper.py b/xfel/util/sublattice_helper.py deleted file mode 100644 index 859e8d877e..0000000000 --- a/xfel/util/sublattice_helper.py +++ /dev/null @@ -1,164 +0,0 @@ -from __future__ import division, print_function, absolute_import -from scitbx.matrix import sqr, row -from rstbx.sublattice_support.change_basis import sublattice_change_of_basis -# for debug: from cctbx.crystal_orientation import crystal_orientation -from cctbx.miller import set as miller_set -from cctbx.crystal import symmetry -import copy, logging -from dials.array_family import flex -from dials.algorithms.indexing.stills_indexer import calc_2D_rmsd_and_displacements - -SL = sublattice_change_of_basis(max_modulus=2) -SLT = list(SL.yield_transformations_ascending_modulus()) -# reorder ersatz yield from the toolbox to conform with organized published list -# Set up the phil file so common sense lookup gives the desired transformation -SLT = [SLT[3], SLT[1], SLT[0], SLT[5], SLT[4], SLT[2], SLT[6]] -assert SLT[0].matS().elems == (2,0,0,0,1,0,0,0,1) # Double the a-axis -assert SLT[1].matS().elems == (1,0,0,0,2,0,0,0,1) # Double the b-axis -assert SLT[2].matS().elems == (1,0,0,0,1,0,0,0,2) # Double the c-axis -assert SLT[3].matS().elems == (2,1,0,0,1,0,0,0,1) # C-face centering -assert SLT[4].matS().elems == (2,0,1,0,1,0,0,0,1) # B-face centering -assert SLT[5].matS().elems == (1,0,0,0,2,1,0,0,1) # A-face centering -assert SLT[6].matS().elems == (2,1,1,0,1,0,0,0,1) # Body centering - -logger = logging.getLogger(__name__) - -def integrate_coset(self, experiments, indexed): - TRANS = self.params.integration.coset.transformation - - # here get a deepcopy that we are not afraid to modify: - experiments_local = copy.deepcopy(experiments) - - logger.info("*" * 80) - logger.info("Coset Reflections for modeling or validating the background") - logger.info("*" * 80) - from dials.algorithms.profile_model.factory import ProfileModelFactory - from dials.algorithms.integration.integrator import create_integrator - - # XXX Fixme later implement support for non-primitive lattices NKS - base_set = miller_set( crystal_symmetry = symmetry( - unit_cell = experiments_local[0].crystal.get_unit_cell(), - space_group = experiments_local[0].crystal.get_space_group()), - indices = indexed["miller_index"] - ) - triclinic = base_set.customized_copy( - crystal_symmetry=symmetry(unit_cell = experiments_local[0].crystal.get_unit_cell(),space_group="P1")) - - # ================ - # Compute the profile model - # Predict the reflections - # Create the integrator - # This creates a reference to the experiment, not a copy: - experiments_local = ProfileModelFactory.create(self.params, experiments_local, indexed) - # for debug SLT[TRANS].show_summary() - - for e in experiments_local: - e.crystal.set_space_group(triclinic.space_group()) - Astar = e.crystal.get_A() - # debug OriAstar = crystal_orientation(Astar,True) - # debug OriAstar.show(legend="old ") - Astarprime = sqr(Astar)* ( sqr(SLT[TRANS]._reindex_N).transpose().inverse() ) - e.crystal.set_A(Astarprime) - # debug OriAstarprime = crystal_orientation(Astarprime,True) - # debug OriAstarprime.show(legend="new ") - - logger.info("Predicting coset reflections") - logger.info("") - predicted = flex.reflection_table.from_predictions_multi( - experiments_local, - dmin=self.params.prediction.d_min, - dmax=self.params.prediction.d_max, - margin=self.params.prediction.margin, - force_static=self.params.prediction.force_static, - ) - logger.info("sublattice total predictions %d"%len(predicted)) - - # filter the sublattice, keep only the coset indices - miller = predicted["miller_index"] - # coset of modulus 2, wherein there is a binary choice - # see Sauter & Zwart, Acta D (2009) 65:553, Table 1; select the valid coset using eqn(5). - coset_select_algorithm_2 = flex.bool() - M_mat = SLT[TRANS].matS() # the transformation - M_p = M_mat.inverse() - for idx in miller: - H_row = row(idx) - h_orig_setting = H_row * M_p - on_coset=False - for icom in h_orig_setting.elems: - if icom.denominator() > 1: on_coset=True; break - coset_select_algorithm_2.append(on_coset) - predicted = predicted.select(coset_select_algorithm_2) - logger.info("of which %d are in coset %d"%(len(predicted), TRANS)) - - logger.info("") - integrator = create_integrator(self.params, experiments_local, predicted) - - # Integrate the reflections - integrated = integrator.integrate() - - # Delete the shoeboxes used for intermediate calculations, if requested - if self.params.integration.debug.delete_shoeboxes and "shoebox" in integrated: - del integrated["shoebox"] - - if self.params.output.composite_output: - if ( - self.params.output.coset_experiments_filename - or self.params.output.coset_filename - ): - assert ( - self.params.output.coset_experiments_filename is not None - and self.params.output.coset_filename is not None - ) - n = len(self.all_coset_experiments) - self.all_coset_experiments.extend(experiments_local) - for i, experiment in enumerate(experiments_local): - refls = integrated.select(integrated["id"] == i) - refls["id"] = flex.int(len(refls), n) - del refls.experiment_identifiers()[i] - refls.experiment_identifiers()[n] = experiment.identifier - self.all_coset_reflections.extend(refls) - n += 1 - else: - # Dump experiments to disk - if self.params.output.coset_experiments_filename: - - experiments_local.as_json(self.params.output.coset_experiments_filename) - - if self.params.output.coset_filename: - # Save the reflections - self.save_reflections( - integrated, self.params.output.coset_filename - ) - - rmsd_indexed, _ = calc_2D_rmsd_and_displacements(indexed) - log_str = "coset RMSD indexed (px): %f\n" % (rmsd_indexed) - log_str += "integrated %d\n"%len(integrated) - for i in range(6): - bright_integrated = integrated.select( - ( - integrated["intensity.sum.value"] - / flex.sqrt(integrated["intensity.sum.variance"]) - ) - >= i - ) - if len(bright_integrated) > 0: - rmsd_integrated, _ = calc_2D_rmsd_and_displacements(bright_integrated) - else: - rmsd_integrated = 0 - log_str += ( - "N reflections integrated at I/sigI >= %d: % 4d, RMSD (px): %f\n" - % (i, len(bright_integrated), rmsd_integrated) - ) - - for crystal_model in experiments_local.crystals(): - if hasattr(crystal_model, "get_domain_size_ang"): - log_str += ( - ". Final ML model: domain size angstroms: %f, half mosaicity degrees: %f" - % ( - crystal_model.get_domain_size_ang(), - crystal_model.get_half_mosaicity_deg(), - ) - ) - - logger.info(log_str) - logger.info("")