diff --git a/CMakeLists.txt b/CMakeLists.txt index c286d39e2..96c6b76d5 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -3,7 +3,7 @@ cmake_minimum_required(VERSION 3.20 FATAL_ERROR) project(dxtbx LANGUAGES C CXX) # Add the included modules -set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_SOURCE_DIR}/cmake/Modules/") +set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_CURRENT_SOURCE_DIR}/cmake/Modules/") # General cmake environment configuration include(SetDefaultBuildRelWithDebInfo) # Default builds to release with debug info diff --git a/newsfragments/621.feature b/newsfragments/621.feature new file mode 100644 index 000000000..394ff8e74 --- /dev/null +++ b/newsfragments/621.feature @@ -0,0 +1 @@ +Add new Beam class "PolychromaticBeam" for polychromatic/multi-wavelength/wide bandpass experiments. diff --git a/newsfragments/626.feature b/newsfragments/626.feature new file mode 100644 index 000000000..c28cc4293 --- /dev/null +++ b/newsfragments/626.feature @@ -0,0 +1 @@ +Update Format handling to reflect move of Eiger detector from PETRA P14 to P13. diff --git a/newsfragments/633.bugfix b/newsfragments/633.bugfix new file mode 100644 index 000000000..7ce5110e5 --- /dev/null +++ b/newsfragments/633.bugfix @@ -0,0 +1 @@ +Slicing of imageset objects is now consistently 0-based, including for the sliced data accessor. Previously, the data accessor had to be accessed with the original index offsets. diff --git a/newsfragments/647.feature b/newsfragments/647.feature new file mode 100644 index 000000000..0ff6b0846 --- /dev/null +++ b/newsfragments/647.feature @@ -0,0 +1 @@ +The ``Beam`` model now has a ``probe`` value to keep track of the type of radiation. diff --git a/newsfragments/652.bugfix b/newsfragments/652.bugfix new file mode 100644 index 000000000..ec7635b59 --- /dev/null +++ b/newsfragments/652.bugfix @@ -0,0 +1 @@ +``dxtbx``: add fix for Eiger / NXmx data from i19-2 to correctly assign the image bit depth diff --git a/newsfragments/653.feature b/newsfragments/653.feature new file mode 100644 index 000000000..22c7d8c0e --- /dev/null +++ b/newsfragments/653.feature @@ -0,0 +1 @@ +``FormatROD``: include support for multi-axis goniometers and faster decompression. diff --git a/src/dxtbx/boost_python/compression.cc b/src/dxtbx/boost_python/compression.cc index ff5d91c5f..6957d1bd1 100644 --- a/src/dxtbx/boost_python/compression.cc +++ b/src/dxtbx/boost_python/compression.cc @@ -164,3 +164,129 @@ unsigned int dxtbx::boost_python::cbf_decompress(const char *packed, return values - original; } + +inline uint32_t read_uint32_from_bytearray(const char *buf) { + // `char` can be signed or unsigned depending on the platform. + // For bit shift operations, we need unsigned values. + // If `char` on the platform is signed, converting directly to "unsigned int" can + // produce huge numbers because modulo 2^n is taken by the integral conversion + // rules. Thus, we have to explicitly cast to `unsigned char` first. + // Then the automatic integral promotion converts them to `int`. + // Note that the unsigned to signed conversion is implementation-dependent + // and might not produce the intended result if two's complement is not used. + // Fortunately, DIALS targets only two's complement. + // https://github.com/cctbx/dxtbx/issues/11#issuecomment-1657809645 + // Moreover, C++20 standarized this: + // https://stackoverflow.com/questions/54947427/going-from-signed-integers-to-unsigned-integers-and-vice-versa-in-c20 + + return ((unsigned char)buf[0]) | (((unsigned char)buf[1]) << 8) + | (((unsigned char)buf[2]) << 16) | (((unsigned char)buf[3]) << 24); +} + +inline uint16_t read_uint16_from_bytearray(const char *buf) { + return ((unsigned char)buf[0]) | ((unsigned char)buf[1] << 8); +} + +void dxtbx::boost_python::rod_TY6_decompress(int *const ret, + const char *const buf_data, + const char *const buf_offsets, + const int slow, + const int fast) { + const size_t BLOCKSIZE = 8; // Codes below assume this is at most 8 + const signed int SHORT_OVERFLOW = 127; // after 127 is subtracted + const signed int LONG_OVERFLOW = 128; + + const size_t nblock = (fast - 1) / (BLOCKSIZE * 2); + const size_t nrest = (fast - 1) % (BLOCKSIZE * 2); + + for (size_t iy = 0; iy < slow; iy++) { + size_t ipos = read_uint32_from_bytearray(buf_offsets + iy * sizeof(uint32_t)); + size_t opos = fast * iy; + + // Values from -127 to +126 (inclusive) are stored with an offset of 127 + // as 0 to 253. 254 and 255 mark short and long overflows. + // Other values ("overflows") are represented in two's complement. + + int firstpx = (unsigned char)buf_data[ipos++] - 127; + if (firstpx == LONG_OVERFLOW) { + // See comments in read_uint32_from_bytearray() about + // the safety of the unsigned to signed conversion. + firstpx = (signed int)read_uint32_from_bytearray(buf_data + ipos); + ipos += 4; + } else if (firstpx == SHORT_OVERFLOW) { + firstpx = (signed short)read_uint16_from_bytearray(buf_data + ipos); + ipos += 2; + } + ret[opos++] = firstpx; + + // For every two blocks + for (int k = 0; k < nblock; k++) { + const size_t bittypes = buf_data[ipos++]; + const size_t nbits[2] = {bittypes & 15, (bittypes >> 4) & 15}; + + // One pixel is stored using `nbit` bits. + // Although `nbit` itself is stored using 4 bits, + // only values 1 (0001b) to 8 (1000b) are allowed. + // Negative values are encoded as follows. (Not 2's complement!) + // - When nbit = 1, the pixel value is 0 or 1 + // - When nbit = 2, the pixel value is -1, 0, 1, 2 + // - When nbit = 3, the pixel value is -3, -2, 1, 0, 1, 2, 3, 4 + // - When nbit - 8, the pixel value is -127, -126, ..., + // 127 (== // SHORT_OVERFLOW), 128 (== LONG_OVERFLOW) + + // Load values + for (int i = 0; i < 2; i++) { + const size_t nbit = nbits[i]; + assert(nbit >= 0 && nbit <= 8); + + int zero_at = 0; + if (nbit > 1) { + zero_at = (1 << (nbit - 1)) - 1; + } + + // Since nbit is at most 8, 8 * 8 (= BLOCKSIZE) = 64 bits are sufficient. + unsigned long long v = 0; + for (int j = 0; j < nbit; j++) { + // Implicit promotion is only up to 32 bits, not 64 bits so we have to be + // explicit. + v |= (long long)((unsigned char)buf_data[ipos++]) << (BLOCKSIZE * j); + } + + const unsigned long long mask = (1 << nbit) - 1; + for (int j = 0; j < BLOCKSIZE; j++) { + ret[opos++] = ((v >> (nbit * j)) & mask) - zero_at; + } + } + + // Apply differences. Load more values when overflown. + for (size_t i = opos - 2 * BLOCKSIZE; i < opos; i++) { + int offset = ret[i]; + + if (offset == LONG_OVERFLOW) { + offset = (signed int)read_uint32_from_bytearray(buf_data + ipos); + ipos += 4; + } else if (offset == SHORT_OVERFLOW) { + offset = (signed short)read_uint16_from_bytearray(buf_data + ipos); + ipos += 2; + } + + ret[i] = offset + ret[i - 1]; + } + } + + for (int i = 0; i < nrest; i++) { + int offset = (unsigned char)buf_data[ipos++] - 127; + + if (offset == LONG_OVERFLOW) { + offset = (signed int)read_uint32_from_bytearray(buf_data + ipos); + ipos += 4; + } else if (offset == SHORT_OVERFLOW) { + offset = (signed short)read_uint16_from_bytearray(buf_data + ipos); + ipos += 2; + } + + ret[opos] = ret[opos - 1] + offset; + opos++; + } + } +} diff --git a/src/dxtbx/boost_python/compression.h b/src/dxtbx/boost_python/compression.h index d046b3965..b50b5e73f 100644 --- a/src/dxtbx/boost_python/compression.h +++ b/src/dxtbx/boost_python/compression.h @@ -6,6 +6,12 @@ namespace dxtbx { namespace boost_python { unsigned int cbf_decompress(const char*, std::size_t, int*, const std::size_t); std::vector cbf_compress(const int*, const std::size_t&); + // Decompress Rigaku Oxford diffractometer TY6 compression + void rod_TY6_decompress(int* const, + const char* const, + const char* const, + const int, + const int); }} // namespace dxtbx::boost_python #endif diff --git a/src/dxtbx/boost_python/ext.cpp b/src/dxtbx/boost_python/ext.cpp index 5273eb69e..70c69c06b 100644 --- a/src/dxtbx/boost_python/ext.cpp +++ b/src/dxtbx/boost_python/ext.cpp @@ -193,6 +193,24 @@ namespace dxtbx { namespace boost_python { return PyBytes_FromStringAndSize(&*packed.begin(), packed.size()); } + // Python entry point to decompress Rigaku Oxford Diffractometer TY6 compression + scitbx::af::flex_int uncompress_rod_TY6(const boost::python::object &data, + const boost::python::object &offsets, + const int &slow, + const int &fast) { + // Cannot I extract const char* directly? + std::string str_data = boost::python::extract(data); + std::string str_offsets = boost::python::extract(offsets); + + scitbx::af::flex_int z((scitbx::af::flex_grid<>(slow, fast)), + scitbx::af::init_functor_null()); + + dxtbx::boost_python::rod_TY6_decompress( + z.begin(), str_data.c_str(), str_offsets.c_str(), slow, fast); + + return z; + } + void init_module() { using namespace boost::python; def("read_uint8", read_uint8, (arg("file"), arg("count"))); @@ -206,6 +224,9 @@ namespace dxtbx { namespace boost_python { def("is_big_endian", is_big_endian); def("uncompress", &uncompress, (arg_("packed"), arg_("slow"), arg_("fast"))); def("compress", &compress); + def("uncompress_rod_TY6", + &uncompress_rod_TY6, + (arg_("data"), arg_("offsets"), arg_("slow"), arg_("fast"))); } BOOST_PYTHON_MODULE(dxtbx_ext) { diff --git a/src/dxtbx/dxtbx_model_ext.pyi b/src/dxtbx/dxtbx_model_ext.pyi index 4cbee00fb..8a4a85934 100644 --- a/src/dxtbx/dxtbx_model_ext.pyi +++ b/src/dxtbx/dxtbx_model_ext.pyi @@ -22,6 +22,8 @@ from scitbx.array_family import shared as flex_shared # Attempt to use the stub typing for flex-inheritance from scitbx.array_family.flex import FlexPlain +from dxtbx_model_ext import Probe # type: ignore + # TypeVar for the set of Experiment models that can be joint-accepted # - profile, imageset and scalingmodel are handled as 'object' TExperimentModel = TypeVar( @@ -113,6 +115,37 @@ class Beam(BeamBase): @staticmethod def from_dict(data: Dict) -> Beam: ... def to_dict(self) -> Dict: ... + @staticmethod + def get_probe_from_name(name: str) -> Probe: ... + +class PolychromaticBeam(Beam): + @overload + def __init__(self, beam: PolychromaticBeam) -> None: ... + @overload + def __init__(self, direction: Vec3Float) -> None: ... + @overload + def __init__( + self, + direction: Vec3Float, + divergence: float, + sigma_divergence: float, + deg: bool = ..., + ) -> None: ... + @overload + def __init__( + self, + direction: Vec3Float, + divergence: float, + sigma_divergence: float, + polarization_normal: Vec3Float, + polarization_fraction: float, + flux: float, + transmission: float, + deg: bool = ..., + ) -> None: ... + @staticmethod + def from_dict(data: Dict) -> PolychromaticBeam: ... + def to_dict(self) -> Dict: ... class CrystalBase: @property diff --git a/src/dxtbx/format/FormatCBFMini.py b/src/dxtbx/format/FormatCBFMini.py index 4561665fd..39060f84a 100644 --- a/src/dxtbx/format/FormatCBFMini.py +++ b/src/dxtbx/format/FormatCBFMini.py @@ -8,6 +8,7 @@ from __future__ import annotations import binascii +import datetime import os import pathlib import sys @@ -73,6 +74,20 @@ def __init__(self, image_file, **kwargs): self._raw_data = None super().__init__(image_file, **kwargs) + @staticmethod + def _get_timestamp_from_raw_header( + header: str | list[str], + ) -> datetime.datetime | None: + """Given a raw header, or lines from, attempt to extract the timestamp field""" + if isinstance(header, str): + header = header.splitlines() + timestamp = None + for record in header: + if len(record[1:].split()) <= 2 and record.count(":") == 2: + timestamp = datetime.datetime.fromisoformat(record[1:].strip()) + break + return timestamp + def _start(self): """Open the image file, read the image header, copy it into a dictionary for future reference.""" diff --git a/src/dxtbx/format/FormatCBFMiniEigerPetraP14.py b/src/dxtbx/format/FormatCBFMiniEigerPetraP14.py index 530948590..36b5cce91 100644 --- a/src/dxtbx/format/FormatCBFMiniEigerPetraP14.py +++ b/src/dxtbx/format/FormatCBFMiniEigerPetraP14.py @@ -3,6 +3,7 @@ from __future__ import annotations +import datetime import sys from dxtbx.format.FormatCBFMiniEiger import FormatCBFMiniEiger @@ -19,11 +20,21 @@ def understand(image_file): header = FormatCBFMiniEiger.get_cbf_header(image_file) + # Valid from 22nd May 2021 + expected_serial = "E-32-0129" + if timestamp := FormatCBFMiniEiger._get_timestamp_from_raw_header(header): + # We have a timestamp. Let's see what detector we should expect + + # Before 22nd May 2021 + if timestamp < datetime.datetime(2021, 5, 22): + expected_serial = "E-32-0107" + + # Find the line recording detector serial, and check for record in header.split("\n"): if ( "# detector" in record.lower() and "eiger" in record.lower() - and "E-32-0107" in record + and expected_serial in record ): return True diff --git a/src/dxtbx/format/FormatGatanDM4.py b/src/dxtbx/format/FormatGatanDM4.py index d5cb31181..acb93f1ef 100644 --- a/src/dxtbx/format/FormatGatanDM4.py +++ b/src/dxtbx/format/FormatGatanDM4.py @@ -21,6 +21,7 @@ ) from dxtbx.format.Format import Format from dxtbx.format.FormatMultiImage import FormatMultiImage +from dxtbx.model.beam import Probe def read_tag(f, byteorder): @@ -358,6 +359,7 @@ def _beam(self): wavelength=wavelength, polarization=(0, 1, 0), polarization_fraction=0.5, + probe=Probe.electron, ) def _scan(self): diff --git a/src/dxtbx/format/FormatMRC.py b/src/dxtbx/format/FormatMRC.py index 0b9b06226..b7d9b30bb 100644 --- a/src/dxtbx/format/FormatMRC.py +++ b/src/dxtbx/format/FormatMRC.py @@ -19,6 +19,7 @@ from dxtbx.format.Format import Format from dxtbx.format.FormatMultiImage import FormatMultiImage from dxtbx.model import ScanFactory +from dxtbx.model.beam import Probe logger = logging.getLogger("dials") @@ -226,6 +227,7 @@ def _beam(self): wavelength=wavelength, polarization=(0, 1, 0), polarization_fraction=0.5, + probe=Probe.electron, ) diff --git a/src/dxtbx/format/FormatNXmxDLSI19_2.py b/src/dxtbx/format/FormatNXmxDLSI19_2.py index 6d7f0c7ad..255330eb4 100644 --- a/src/dxtbx/format/FormatNXmxDLSI19_2.py +++ b/src/dxtbx/format/FormatNXmxDLSI19_2.py @@ -8,7 +8,7 @@ from libtbx import Auto from dxtbx.format.FormatNXmx import FormatNXmx -from dxtbx.format.FormatNXmxDLS import FormatNXmxDLS +from dxtbx.format.FormatNXmxDLS import FormatNXmxDLS, get_bit_depth_from_meta from dxtbx.masking import GoniometerMaskerFactory @@ -41,6 +41,11 @@ def __init__(self, image_file, **kwargs): """Initialise the image structure from the given file.""" self._dynamic_shadowing = self.has_dynamic_shadowing(**kwargs) super().__init__(image_file, **kwargs) + try: + self._bit_depth_readout = get_bit_depth_from_meta(self._meta) + except Exception: + self._bit_depth_readout = 16 + def get_goniometer_shadow_masker(self, goniometer=None): """Apply the dynamic mask for a diamond anvil cell with a 76° aperture.""" diff --git a/src/dxtbx/format/FormatROD.py b/src/dxtbx/format/FormatROD.py index 417f99f1a..5ce6bd814 100644 --- a/src/dxtbx/format/FormatROD.py +++ b/src/dxtbx/format/FormatROD.py @@ -12,9 +12,7 @@ from __future__ import annotations __author__ = "David Waterman, Takanori Nakane" -__copyright__ = ( - "Copyright 2018 United Kingdom Research and Innovation & 2022 Takanori Nakane" -) +__copyright__ = "Copyright 2018-2023 United Kingdom Research and Innovation & 2022-2023 Takanori Nakane" __license__ = "BSD 3-clause" import re @@ -23,7 +21,9 @@ import numpy as np from scitbx.array_family import flex +from scitbx.math import r3_rotation_axis_and_angle_as_matrix +from dxtbx.ext import uncompress_rod_TY6 from dxtbx.format.Format import Format @@ -153,7 +153,7 @@ def _read_binary_header( f.seek(offset + general_nbytes + special_nbytes + 640) # detector rotation in degrees along e1, e2, e3 detector_rotns = struct.unpack("> 4) & 15) ipos += 1 - # ipos_bit = ipos * 8 for i in range(2): nbit = nbits[i] @@ -455,5 +488,12 @@ def decode_TY6_oneline(self, linedata, w): reader = FormatROD(arg) print(reader._txt_header) print(reader._bin_header) + print() + print( + "Starting angles in degrees (OMEGA, THETA, KAPPA=CHI, PHI, OMEGA PRIME, THETA PRIME)" + ) + print(reader._gonio_start_angles) + print("Ending angles in degrees") + print(reader._gonio_end_angles) else: print("Unsupported format") diff --git a/src/dxtbx/format/FormatSMVTimePix_SU.py b/src/dxtbx/format/FormatSMVTimePix_SU.py index 223061836..b3492ff56 100644 --- a/src/dxtbx/format/FormatSMVTimePix_SU.py +++ b/src/dxtbx/format/FormatSMVTimePix_SU.py @@ -12,6 +12,7 @@ from scitbx import matrix from dxtbx.format.FormatSMV import FormatSMV +from dxtbx.model.beam import Probe from dxtbx.model.detector import Detector @@ -76,6 +77,7 @@ def _beam(self): wavelength=wavelength, polarization=(0, 1, 0), polarization_fraction=0.5, + probe=Probe.electron, ) def _scan(self): diff --git a/src/dxtbx/imageset.py b/src/dxtbx/imageset.py index 64c51e354..564ced681 100644 --- a/src/dxtbx/imageset.py +++ b/src/dxtbx/imageset.py @@ -303,29 +303,12 @@ def __getitem__(self, item): if not isinstance(item, slice): return self.get_corrected_data(item) else: - offset = self.get_scan().get_batch_offset() if item.step is not None: raise IndexError("Sequences must be sequential") - # nasty workaround for https://github.com/dials/dials/issues/1153 - # slices with -1 in them are meaningful :-/ so grab the original - # constructor arguments of the slice object. - # item.start and item.stop may have been compromised at this point. - if offset < 0: - start, stop, step = item.__reduce__()[1] - if start is None: - start = 0 - else: - start -= offset - if stop is None: - stop = len(self) - else: - stop -= offset - else: - start = item.start or 0 - stop = item.stop or (len(self) + offset) - start -= offset - stop -= offset + start = item.start or 0 + stop = item.stop or len(self) + if self.data().has_single_file_reader(): reader = self.reader().copy(self.reader().paths(), stop - start) else: diff --git a/src/dxtbx/model/__init__.py b/src/dxtbx/model/__init__.py index 26a2fc4d8..a7e3ada72 100644 --- a/src/dxtbx/model/__init__.py +++ b/src/dxtbx/model/__init__.py @@ -45,6 +45,7 @@ OffsetPxMmStrategy, Panel, ParallaxCorrectedPxMmStrategy, + PolychromaticBeam, PxMmStrategy, Scan, ScanBase, @@ -80,6 +81,7 @@ OffsetPxMmStrategy, Panel, ParallaxCorrectedPxMmStrategy, + PolychromaticBeam, PxMmStrategy, Scan, ScanBase, @@ -97,6 +99,7 @@ __all__ = ( "Beam", "BeamBase", + "PolychromaticBeam", "BeamFactory", "Crystal", "CrystalBase", diff --git a/src/dxtbx/model/beam.h b/src/dxtbx/model/beam.h index 9903341ec..6df5325d1 100644 --- a/src/dxtbx/model/beam.h +++ b/src/dxtbx/model/beam.h @@ -13,6 +13,7 @@ #include #include +#include #include #include #include @@ -24,6 +25,9 @@ namespace dxtbx { namespace model { using scitbx::vec3; + // probe type enumeration + enum Probe { xray = 1, electron = 2, neutron = 3 }; + /** Base class for beam objects */ class BeamBase { public: @@ -44,6 +48,8 @@ namespace dxtbx { namespace model { virtual std::size_t get_num_scan_points() const = 0; virtual scitbx::af::shared > get_s0_at_scan_points() const = 0; virtual vec3 get_s0_at_scan_point(std::size_t index) const = 0; + virtual Probe get_probe() const = 0; + virtual std::string get_probe_name() const = 0; virtual void set_direction(vec3 direction) = 0; virtual void set_wavelength(double wavelength) = 0; @@ -59,6 +65,7 @@ namespace dxtbx { namespace model { virtual void set_transmission(double transmission) = 0; virtual void set_s0_at_scan_points( const scitbx::af::const_ref > &s0) = 0; + virtual void set_probe(Probe probe) = 0; virtual void reset_scan_points() = 0; virtual bool is_similar_to(const BeamBase &rhs, @@ -82,7 +89,8 @@ namespace dxtbx { namespace model { polarization_normal_(0.0, 1.0, 0.0), polarization_fraction_(0.999), flux_(0), - transmission_(1.0) {} + transmission_(1.0), + probe_(Probe::xray) {} /** * @param s0 The incident beam vector. @@ -93,7 +101,8 @@ namespace dxtbx { namespace model { polarization_normal_(0.0, 1.0, 0.0), polarization_fraction_(0.999), flux_(0), - transmission_(1.0) { + transmission_(1.0), + probe_(Probe::xray) { DXTBX_ASSERT(s0.length() > 0); wavelength_ = 1.0 / s0.length(); direction_ = -s0.normalize(); @@ -110,7 +119,8 @@ namespace dxtbx { namespace model { polarization_normal_(0.0, 1.0, 0.0), polarization_fraction_(0.999), flux_(0), - transmission_(1.0) { + transmission_(1.0), + probe_(Probe::xray) { DXTBX_ASSERT(direction.length() > 0); direction_ = direction.normalize(); } @@ -126,7 +136,8 @@ namespace dxtbx { namespace model { polarization_normal_(0.0, 1.0, 0.0), polarization_fraction_(0.999), flux_(0), - transmission_(1.0) { + transmission_(1.0), + probe_(Probe::xray) { DXTBX_ASSERT(s0.length() > 0); wavelength_ = 1.0 / s0.length(); direction_ = -s0.normalize(); @@ -148,7 +159,8 @@ namespace dxtbx { namespace model { polarization_normal_(0.0, 1.0, 0.0), polarization_fraction_(0.999), flux_(0), - transmission_(1.0) { + transmission_(1.0), + probe_(Probe::xray) { DXTBX_ASSERT(direction.length() > 0); direction_ = direction.normalize(); } @@ -162,6 +174,7 @@ namespace dxtbx { namespace model { * @param polarization_fraction The polarization fraction * @param flux The beam flux * @param transmission The beam transmission + * @param probe The probe value */ Beam(vec3 direction, double wavelength, @@ -170,14 +183,16 @@ namespace dxtbx { namespace model { vec3 polarization_normal, double polarization_fraction, double flux, - double transmission) + double transmission, + Probe probe) : wavelength_(wavelength), divergence_(divergence), sigma_divergence_(sigma_divergence), polarization_normal_(polarization_normal), polarization_fraction_(polarization_fraction), flux_(flux), - transmission_(transmission) { + transmission_(transmission), + probe_(probe) { DXTBX_ASSERT(direction.length() > 0); direction_ = direction.normalize(); } @@ -188,7 +203,7 @@ namespace dxtbx { namespace model { return direction_; } - double get_wavelength() const { + virtual double get_wavelength() const { return wavelength_; } @@ -207,16 +222,16 @@ namespace dxtbx { namespace model { direction_ = direction.normalize(); } - void set_wavelength(double wavelength) { + virtual void set_wavelength(double wavelength) { wavelength_ = wavelength; } - vec3 get_s0() const { + virtual vec3 get_s0() const { DXTBX_ASSERT(wavelength_ != 0.0); return -direction_ * 1.0 / wavelength_; } - void set_s0(vec3 s0) { + virtual void set_s0(vec3 s0) { DXTBX_ASSERT(s0.length() > 0); direction_ = -s0.normalize(); wavelength_ = 1.0 / s0.length(); @@ -289,11 +304,49 @@ namespace dxtbx { namespace model { return s0_at_scan_points_[index]; } + Probe get_probe() const { + return probe_; + } + + std::string get_probe_name() const { + // Return a name that matches NeXus definitions from + // https://manual.nexusformat.org/classes/base_classes/NXsource.html + switch (probe_) { + case xray: + return std::string("x-ray"); + case electron: + return std::string("electron"); + case neutron: + return std::string("neutron"); + default: + DXTBX_ERROR("Unknown probe type"); + } + } + + static Probe get_probe_from_name(const std::string probe) { + // Return a Probe matched to NeXus definitions from + // https://manual.nexusformat.org/classes/base_classes/NXsource.html + + if (probe == "x-ray") { + return Probe::xray; + } else if (probe == "electron") { + return Probe::electron; + } else if (probe == "neutron") { + return Probe::neutron; + } + + DXTBX_ERROR("Unknown probe " + probe); + } + + void set_probe(Probe probe) { + probe_ = probe; + } + void reset_scan_points() { s0_at_scan_points_.clear(); } - bool operator==(const BeamBase &rhs) const { + virtual bool operator==(const BeamBase &rhs) const { double eps = 1.0e-6; // scan-varying model checks @@ -324,14 +377,15 @@ namespace dxtbx { namespace model { angle_safe(polarization_normal_, rhs.get_polarization_normal())) <= eps && std::abs(polarization_fraction_ - rhs.get_polarization_fraction()) - <= eps; + <= eps + && (probe_ == rhs.get_probe()); } - bool is_similar_to(const BeamBase &rhs, - double wavelength_tolerance, - double direction_tolerance, - double polarization_normal_tolerance, - double polarization_fraction_tolerance) const { + virtual bool is_similar_to(const BeamBase &rhs, + double wavelength_tolerance, + double direction_tolerance, + double polarization_normal_tolerance, + double polarization_fraction_tolerance) const { // scan varying model checks if (get_num_scan_points() != rhs.get_num_scan_points()) { return false; @@ -361,7 +415,8 @@ namespace dxtbx { namespace model { angle_safe(polarization_normal_, rhs.get_polarization_normal())) <= polarization_normal_tolerance && std::abs(polarization_fraction_ - rhs.get_polarization_fraction()) - <= polarization_fraction_tolerance; + <= polarization_fraction_tolerance + && (probe_ == rhs.get_probe()); } bool operator!=(const BeamBase &rhs) const { @@ -375,8 +430,7 @@ namespace dxtbx { namespace model { friend std::ostream &operator<<(std::ostream &os, const Beam &b); - private: - double wavelength_; + protected: vec3 direction_; double divergence_; double sigma_divergence_; @@ -384,12 +438,17 @@ namespace dxtbx { namespace model { double polarization_fraction_; double flux_; double transmission_; + Probe probe_; + + private: + double wavelength_; scitbx::af::shared > s0_at_scan_points_; }; /** Print beam information */ inline std::ostream &operator<<(std::ostream &os, const Beam &b) { os << "Beam:\n"; + os << " probe: " << b.get_probe_name() << "\n"; os << " wavelength: " << b.get_wavelength() << "\n"; os << " sample to source direction : " << b.get_sample_to_source_direction().const_ref() << "\n"; @@ -403,6 +462,180 @@ namespace dxtbx { namespace model { return os; } + class PolychromaticBeam : public Beam { + public: + PolychromaticBeam() { + set_direction(vec3(0.0, 0.0, 1.0)); + set_divergence(0.0); + set_sigma_divergence(0.0); + set_polarization_normal(vec3(0.0, 1.0, 0.0)); + set_polarization_fraction(0.5); + set_flux(0); + set_transmission(1.0); + set_probe(Probe::xray); + } + + /** + * @param direction The beam direction pointing sample to source + */ + PolychromaticBeam(vec3 direction) { + DXTBX_ASSERT(direction.length() > 0); + direction_ = direction.normalize(); + set_divergence(0.0); + set_sigma_divergence(0.0); + set_polarization_normal(vec3(0.0, 1.0, 0.0)); + set_polarization_fraction(0.5); + set_flux(0); + set_transmission(1.0); + set_probe(Probe::xray); + } + + /** + * @param direction The beam direction pointing sample to source + * @param divergence The beam divergence + * @param sigma_divergence The standard deviation of the beam divergence + */ + PolychromaticBeam(vec3 direction, + double divergence, + double sigma_divergence) { + DXTBX_ASSERT(direction.length() > 0); + direction_ = direction.normalize(); + set_divergence(divergence); + set_sigma_divergence(sigma_divergence); + set_polarization_normal(vec3(0.0, 1.0, 0.0)); + set_polarization_fraction(0.5); + set_flux(0); + set_transmission(1.0); + set_probe(Probe::xray); + } + + /** + * @param direction The beam direction pointing sample to source + * @param divergence The beam divergence + * @param sigma_divergence The standard deviation of the beam divergence + * @param polarization_normal The polarization plane + * @param polarization_fraction The polarization fraction + * @param flux The beam flux + * @param transmission The beam transmission + * @param probe The probe value + */ + PolychromaticBeam(vec3 direction, + double divergence, + double sigma_divergence, + vec3 polarization_normal, + double polarization_fraction, + double flux, + double transmission, + Probe probe) { + DXTBX_ASSERT(direction.length() > 0); + direction_ = direction.normalize(); + set_divergence(divergence); + set_sigma_divergence(sigma_divergence); + set_polarization_normal(polarization_normal); + set_polarization_fraction(polarization_fraction); + set_flux(flux); + set_transmission(transmission); + set_probe(probe); + } + + double get_wavelength() const { + throw DXTBX_ERROR("PolychromaticBeam has no fixed wavelength"); + return -1.; + } + + void set_wavelength(double wavelength) { + throw DXTBX_ERROR("PolychromaticBeam has no fixed wavelength"); + } + + vec3 get_s0() const { + throw DXTBX_ERROR("PolychromaticBeam has no fixed s0"); + return vec3(0., 0., 0.); + } + + void set_s0(vec3 s0) { + throw DXTBX_ERROR("PolychromaticBeam has no fixed s0"); + } + + std::size_t get_num_scan_points() const { + throw DXTBX_ERROR("PolychromaticBeam has no fixed s0"); + return 1; + } + + void set_s0_at_scan_points(const scitbx::af::const_ref > &s0) { + throw DXTBX_ERROR("PolychromaticBeam has no fixed s0"); + } + + scitbx::af::shared > get_s0_at_scan_points() const { + throw DXTBX_ERROR("PolychromaticBeam has no fixed s0"); + return scitbx::af::shared >(1, (0., 0., 0.)); + } + + vec3 get_s0_at_scan_point(std::size_t index) const { + throw DXTBX_ERROR("PolychromaticBeam has no fixed s0"); + return vec3(0., 0., 0.); + } + + void reset_scan_points() { + throw DXTBX_ERROR("PolychromaticBeam has no fixed s0"); + } + + bool operator==(const BeamBase &rhs) const { + double eps = 1.0e-6; + + return std::abs(angle_safe(direction_, rhs.get_sample_to_source_direction())) + <= eps + && std::abs(divergence_ - rhs.get_divergence()) <= eps + && std::abs(sigma_divergence_ - rhs.get_sigma_divergence()) <= eps + && std::abs( + angle_safe(polarization_normal_, rhs.get_polarization_normal())) + <= eps + && std::abs(polarization_fraction_ - rhs.get_polarization_fraction()) + <= eps + && (probe_ == rhs.get_probe()); + } + + bool is_similar_to(const BeamBase &rhs, + double wavelength_tolerance, + double direction_tolerance, + double polarization_normal_tolerance, + double polarization_fraction_tolerance) const { + return is_similar_to(rhs, + direction_tolerance, + polarization_normal_tolerance, + polarization_fraction_tolerance); + } + + bool is_similar_to(const BeamBase &rhs, + double direction_tolerance, + double polarization_normal_tolerance, + double polarization_fraction_tolerance) const { + return std::abs(angle_safe(direction_, rhs.get_sample_to_source_direction())) + <= direction_tolerance + && std::abs( + angle_safe(polarization_normal_, rhs.get_polarization_normal())) + <= polarization_normal_tolerance + && std::abs(polarization_fraction_ - rhs.get_polarization_fraction()) + <= polarization_fraction_tolerance + && (probe_ == rhs.get_probe()); + } + }; + + /** Print beam information */ + inline std::ostream &operator<<(std::ostream &os, const PolychromaticBeam &b) { + os << "Beam:\n"; + os << " probe: " << b.get_probe_name() << "\n"; + os << " sample to source direction : " + << b.get_sample_to_source_direction().const_ref() << "\n"; + os << " divergence: " << b.get_divergence() << "\n"; + os << " sigma divergence: " << b.get_sigma_divergence() << "\n"; + os << " polarization normal: " << b.get_polarization_normal().const_ref() + << "\n"; + os << " polarization fraction: " << b.get_polarization_fraction() << "\n"; + os << " flux: " << b.get_flux() << "\n"; + os << " transmission: " << b.get_transmission() << "\n"; + return os; + } + }} // namespace dxtbx::model #endif // DXTBX_MODEL_BEAM_H diff --git a/src/dxtbx/model/beam.py b/src/dxtbx/model/beam.py index 870dd9c4b..a7360e9e5 100644 --- a/src/dxtbx/model/beam.py +++ b/src/dxtbx/model/beam.py @@ -1,15 +1,18 @@ from __future__ import annotations import math +from typing import Tuple import pycbf import libtbx.phil try: - from ..dxtbx_model_ext import Beam + from ..dxtbx_model_ext import Beam, PolychromaticBeam, Probe except ModuleNotFoundError: - from dxtbx_model_ext import Beam # type: ignore + from dxtbx_model_ext import Beam, PolychromaticBeam, Probe # type: ignore + +Vec3Float = Tuple[float, float, float] beam_phil_scope = libtbx.phil.parse( """ @@ -17,6 +20,15 @@ .expert_level = 1 .short_caption = "Beam overrides" { + type = *monochromatic polychromatic + .type = choice + .help = "Override the beam type" + .short_caption = "beam_type" + + probe = *x-ray electron neutron + .type = choice + .help = "Override the beam probe" + .short_caption = "beam_probe" wavelength = None .type = float @@ -27,6 +39,14 @@ .help = "Override the sample to source direction" .short_caption = "Sample to source direction" + divergence = None + .type = float + .help = "Override the beam divergence" + + sigma_divergence = None + .type = float + .help = "Override the beam sigma divergence" + polarization_normal = None .type = floats(size=3) .help = "Override the polarization normal" @@ -57,62 +77,81 @@ class BeamFactory: will be used, otherwise simplified descriptions can be applied.""" @staticmethod - def from_phil(params, reference=None): + def from_phil( + params: libtbx.phil.scope_extract, + reference: Beam | PolychromaticBeam | None = None, + ) -> Beam | PolychromaticBeam: """ Convert the phil parameters into a beam model - """ + # Check the input if reference is None: - beam = Beam() + beam = ( + PolychromaticBeam() if params.beam.type == "polychromatic" else Beam() + ) else: beam = reference # Set the parameters - if params.beam.wavelength is not None: - beam.set_wavelength(params.beam.wavelength) - elif reference is None: - raise RuntimeError("No wavelength set") + if params.beam.type == "monochromatic": + if params.beam.wavelength is not None: + beam.set_wavelength(params.beam.wavelength) + elif reference is None: + raise RuntimeError("No wavelength set") if params.beam.direction is not None: beam.set_direction(params.beam.direction) elif reference is None: raise RuntimeError("No beam direction set") + + if params.beam.divergence is not None: + beam.set_divergence(params.beam.divergence) + if params.beam.sigma_divergence is not None: + beam.set_sigma_divergence(params.beam.sigma_divergence) if params.beam.polarization_normal is not None: beam.set_polarization_normal(params.beam.polarization_normal) if params.beam.polarization_fraction is not None: beam.set_polarization_fraction(params.beam.polarization_fraction) + if params.beam.transmission is not None: + beam.set_transmission(params.beam.transmission) + if params.beam.flux is not None: + beam.set_flux(params.beam.flux) + beam.set_probe(Beam.get_probe_from_name(params.beam.probe)) return beam @staticmethod - def from_dict(d, t=None): - """Convert the dictionary to a beam model + def from_dict(dict: dict, template: dict = None) -> Beam | PolychromaticBeam: + """Convert the dictionary to a beam model""" - Params: - d The dictionary of parameters - t The template dictionary to use + if template is not None: + if "__id__" in dict and "__id__" in template: + assert ( + dict["__id__"] == template["__id__"] + ), "Beam and template dictionaries are not the same type." - Returns: - The beam model - """ - if d is None and t is None: + if dict is None and template is None: return None - joint = t.copy() if t else {} - joint.update(d) + joint = template.copy() if template else {} + joint.update(dict) # Create the model from the joint dictionary + if "probe" not in joint: + joint["probe"] = "x-ray" + if joint.get("__id__") == "polychromatic": + return PolychromaticBeam.from_dict(joint) + return Beam.from_dict(joint) @staticmethod def make_beam( - sample_to_source=None, - wavelength=None, - s0=None, - unit_s0=None, - divergence=None, - sigma_divergence=None, - ): - + sample_to_source: Vec3Float = None, + wavelength: float = None, + s0: Vec3Float = None, + unit_s0: Vec3Float = None, + divergence: float = None, + sigma_divergence: float = None, + ) -> Beam: if divergence is None or sigma_divergence is None: divergence = 0.0 sigma_divergence = 0.0 @@ -137,19 +176,44 @@ def make_beam( assert s0 return Beam(tuple(map(float, s0))) + @staticmethod + def make_polychromatic_beam( + direction: Vec3Float, + divergence: float = 0.0, + sigma_divergence: float = 0.0, + polarization_normal: Vec3Float = (0.0, 1.0, 0.0), + polarization_fraction: float = 0.5, + flux: float = 0.0, + transmission: float = 1.0, + probe: Probe = Probe.xray, + deg: bool = True, + ) -> PolychromaticBeam: + return PolychromaticBeam( + tuple(map(float, direction)), + float(divergence), + float(sigma_divergence), + tuple(map(float, polarization_normal)), + float(polarization_fraction), + float(flux), + float(transmission), + probe, + bool(deg), + ) + @staticmethod def make_polarized_beam( - sample_to_source=None, - wavelength=None, - s0=None, - unit_s0=None, - polarization=None, - polarization_fraction=None, - divergence=None, - sigma_divergence=None, - flux=None, - transmission=None, - ): + sample_to_source: Vec3Float = None, + wavelength: float = None, + s0: Vec3Float = None, + unit_s0: Vec3Float = None, + polarization: Vec3Float = None, + polarization_fraction: float = None, + divergence: float = None, + sigma_divergence: float = None, + flux: float = None, + transmission: float = None, + probe: Probe = Probe.xray, + ) -> Beam: assert polarization assert 0.0 <= polarization_fraction <= 1.0 @@ -173,6 +237,7 @@ def make_polarized_beam( float(polarization_fraction), float(flux), float(transmission), + probe, ) elif unit_s0: assert wavelength @@ -185,21 +250,27 @@ def make_polarized_beam( float(polarization_fraction), float(flux), float(transmission), + probe, ) else: assert s0 + sum_sq_s0 = s0[0] ** 2 + s0[1] ** 2 + s0[2] ** 2 + assert sum_sq_s0 > 0 + wavelength = 1.0 / math.sqrt(sum_sq_s0) return Beam( tuple(map(float, s0)), + wavelength, float(divergence), float(sigma_divergence), tuple(map(float, polarization)), float(polarization_fraction), float(flux), float(transmission), + probe, ) @staticmethod - def simple(wavelength): + def simple(wavelength: float) -> Beam: """Construct a beam object on the principle that the beam is aligned with the +z axis, as is quite normal. Also assume the beam has polarization fraction 0.999 and is polarized in the x-z plane, unless @@ -219,7 +290,7 @@ def simple(wavelength): ) @staticmethod - def simple_directional(sample_to_source, wavelength): + def simple_directional(sample_to_source: Vec3Float, wavelength: float) -> Beam: """Construct a beam with direction and wavelength.""" if wavelength > 0.05: @@ -236,8 +307,11 @@ def simple_directional(sample_to_source, wavelength): @staticmethod def complex( - sample_to_source, polarization_fraction, polarization_plane_normal, wavelength - ): + sample_to_source: Vec3Float, + polarization_fraction: float, + polarization_plane_normal: Vec3Float, + wavelength: float, + ) -> Beam: """Full access to the constructor for cases where we do know everything that we need...""" @@ -249,7 +323,7 @@ def complex( ) @staticmethod - def imgCIF(cif_file): + def imgCIF(cif_file: str) -> Beam: """Initialize a detector model from an imgCIF file. N.B. the definition of the polarization plane is not completely helpful in this - it is the angle between the polarization plane and the @@ -263,7 +337,7 @@ def imgCIF(cif_file): return result @staticmethod - def imgCIF_H(cbf_handle): + def imgCIF_H(cbf_handle: pycbf.cbf_handle_struct) -> Beam: """Initialize a detector model from an imgCIF file. N.B. the definition of the polarization plane is not completely helpful in this - it is the angle between the polarization plane and the diff --git a/src/dxtbx/model/boost_python/beam.cc b/src/dxtbx/model/boost_python/beam.cc index bf709d526..86cdab927 100644 --- a/src/dxtbx/model/boost_python/beam.cc +++ b/src/dxtbx/model/boost_python/beam.cc @@ -17,6 +17,7 @@ #include #include #include +#include namespace dxtbx { namespace model { namespace boost_python { namespace beam_detail { @@ -25,6 +26,8 @@ namespace dxtbx { namespace model { namespace boost_python { using scitbx::deg_as_rad; using scitbx::rad_as_deg; + // Beam + std::string beam_to_string(const Beam &beam) { std::stringstream ss; ss << beam; @@ -40,7 +43,8 @@ namespace dxtbx { namespace model { namespace boost_python { obj.get_polarization_normal(), obj.get_polarization_fraction(), obj.get_flux(), - obj.get_transmission()); + obj.get_transmission(), + obj.get_probe()); } static boost::python::tuple getstate(boost::python::object obj) { @@ -107,6 +111,7 @@ namespace dxtbx { namespace model { namespace boost_python { double polarization_fraction, double flux, double transmission, + Probe probe, bool deg) { Beam *beam = NULL; if (deg) { @@ -117,7 +122,8 @@ namespace dxtbx { namespace model { namespace boost_python { polarization_normal, polarization_fraction, flux, - transmission); + transmission, + probe); } else { beam = new Beam(sample_to_source, wavelength, @@ -126,7 +132,8 @@ namespace dxtbx { namespace model { namespace boost_python { polarization_normal, polarization_fraction, flux, - transmission); + transmission, + probe); } return beam; } @@ -180,6 +187,7 @@ namespace dxtbx { namespace model { namespace boost_python { template <> boost::python::dict to_dict(const Beam &obj) { boost::python::dict result; + result["__id__"] = "monochromatic"; result["direction"] = obj.get_sample_to_source_direction(); result["wavelength"] = obj.get_wavelength(); result["divergence"] = rad_as_deg(obj.get_divergence()); @@ -188,6 +196,7 @@ namespace dxtbx { namespace model { namespace boost_python { result["polarization_fraction"] = obj.get_polarization_fraction(); result["flux"] = obj.get_flux(); result["transmission"] = obj.get_transmission(); + result["probe"] = obj.get_probe_name(); if (obj.get_num_scan_points() > 0) { boost::python::list l; scitbx::af::shared > s0_at_scan_points = obj.get_s0_at_scan_points(); @@ -212,7 +221,9 @@ namespace dxtbx { namespace model { namespace boost_python { obj.get("polarization_normal", vec3(0.0, 1.0, 0.0))), boost::python::extract(obj.get("polarization_fraction", 0.999)), boost::python::extract(obj.get("flux", 0)), - boost::python::extract(obj.get("transmission", 1))); + boost::python::extract(obj.get("transmission", 1)), + Beam::get_probe_from_name( + boost::python::extract(obj.get("probe", "x-ray")))); if (obj.has_key("s0_at_scan_points")) { boost::python::list s0_at_scan_points = boost::python::extract(obj["s0_at_scan_points"]); @@ -221,9 +232,117 @@ namespace dxtbx { namespace model { namespace boost_python { return b; } + // PolychromaticBeam + + std::string PolychromaticBeam_to_string(const PolychromaticBeam &beam) { + std::stringstream ss; + ss << beam; + return ss.str(); + } + + struct PolychromaticBeamPickleSuite : boost::python::pickle_suite { + static boost::python::tuple getinitargs(const PolychromaticBeam &obj) { + return boost::python::make_tuple(obj.get_sample_to_source_direction(), + obj.get_divergence(), + obj.get_sigma_divergence(), + obj.get_polarization_normal(), + obj.get_polarization_fraction(), + obj.get_flux(), + obj.get_transmission(), + obj.get_probe()); + } + }; + + static PolychromaticBeam *make_PolychromaticBeam(vec3 direction) { + return new PolychromaticBeam(direction); + } + + static PolychromaticBeam *make_PolychromaticBeam_w_divergence(vec3 direction, + double divergence, + double sigma_divergence, + bool deg) { + PolychromaticBeam *beam = NULL; + if (deg) { + beam = new PolychromaticBeam( + direction, deg_as_rad(divergence), deg_as_rad(sigma_divergence)); + } else { + beam = new PolychromaticBeam(direction, divergence, sigma_divergence); + } + return beam; + } + + static PolychromaticBeam *make_PolychromaticBeam_w_all( + vec3 direction, + double divergence, + double sigma_divergence, + vec3 polarization_normal, + double polarization_fraction, + double flux, + double transmission, + Probe probe, + bool deg) { + PolychromaticBeam *beam = NULL; + if (deg) { + beam = new PolychromaticBeam(direction, + deg_as_rad(divergence), + deg_as_rad(sigma_divergence), + polarization_normal, + polarization_fraction, + flux, + transmission, + probe); + } else { + beam = new PolychromaticBeam(direction, + divergence, + sigma_divergence, + polarization_normal, + polarization_fraction, + flux, + transmission, + probe); + } + return beam; + } + + template <> + boost::python::dict to_dict(const PolychromaticBeam &obj) { + boost::python::dict result; + result["__id__"] = "polychromatic"; + result["direction"] = obj.get_sample_to_source_direction(); + result["divergence"] = rad_as_deg(obj.get_divergence()); + result["sigma_divergence"] = rad_as_deg(obj.get_sigma_divergence()); + result["polarization_normal"] = obj.get_polarization_normal(); + result["polarization_fraction"] = obj.get_polarization_fraction(); + result["flux"] = obj.get_flux(); + result["transmission"] = obj.get_transmission(); + result["probe"] = obj.get_probe_name(); + return result; + } + + template <> + PolychromaticBeam *from_dict(boost::python::dict obj) { + PolychromaticBeam *b = new PolychromaticBeam( + boost::python::extract >(obj["direction"]), + deg_as_rad(boost::python::extract(obj.get("divergence", 0.0))), + deg_as_rad(boost::python::extract(obj.get("sigma_divergence", 0.0))), + boost::python::extract >( + obj.get("polarization_normal", vec3(0.0, 1.0, 0.0))), + boost::python::extract(obj.get("polarization_fraction", 0.999)), + boost::python::extract(obj.get("flux", 0)), + boost::python::extract(obj.get("transmission", 1)), + Beam::get_probe_from_name( + boost::python::extract(obj.get("probe", "x-ray")))); + return b; + } + void export_beam() { using namespace beam_detail; + enum_("Probe") + .value("xray", xray) + .value("electron", electron) + .value("neutron", neutron); + class_("BeamBase", no_init) .def("get_sample_to_source_direction", &BeamBase::get_sample_to_source_direction) .def("set_direction", &BeamBase::set_direction) @@ -254,6 +373,9 @@ namespace dxtbx { namespace model { namespace boost_python { .def("set_s0_at_scan_points", &Beam_set_s0_at_scan_points_from_list) .def("get_s0_at_scan_points", &BeamBase::get_s0_at_scan_points) .def("get_s0_at_scan_point", &BeamBase::get_s0_at_scan_point) + .def("get_probe", &BeamBase::get_probe) + .def("get_probe_name", &BeamBase::get_probe_name) + .def("set_probe", &BeamBase::set_probe) .def("reset_scan_points", &BeamBase::reset_scan_points) .def("rotate_around_origin", &rotate_around_origin, @@ -298,13 +420,49 @@ namespace dxtbx { namespace model { namespace boost_python { arg("polarization_fraction"), arg("flux"), arg("transmission"), + arg("probe") = Probe::xray, arg("deg") = true))) .def("__str__", &beam_to_string) .def("to_dict", &to_dict) .def("from_dict", &from_dict, return_value_policy()) .staticmethod("from_dict") + .def("get_probe_from_name", &Beam::get_probe_from_name) + .staticmethod("get_probe_from_name") .def_pickle(BeamPickleSuite()); + class_, bases >( + "PolychromaticBeam") + .def(init()) + .def("__init__", + make_constructor( + &make_PolychromaticBeam, default_call_policies(), (arg("direction")))) + .def("__init__", + make_constructor(&make_PolychromaticBeam_w_divergence, + default_call_policies(), + (arg("direction"), + arg("divergence"), + arg("sigma_divergence"), + arg("deg") = true))) + .def("__init__", + make_constructor(&make_PolychromaticBeam_w_all, + default_call_policies(), + (arg("direction"), + arg("divergence"), + arg("sigma_divergence"), + arg("polarization_normal"), + arg("polarization_fraction"), + arg("flux"), + arg("transmission"), + arg("probe") = Probe::xray, + arg("deg") = true))) + .def("__str__", &PolychromaticBeam_to_string) + .def("to_dict", &to_dict) + .def("from_dict", + &from_dict, + return_value_policy()) + .staticmethod("from_dict") + .def_pickle(PolychromaticBeamPickleSuite()); + scitbx::af::boost_python::flex_wrapper::plain("flex_Beam"); } diff --git a/src/dxtbx/model/experiment_list.py b/src/dxtbx/model/experiment_list.py index d415e1ccd..ec836614f 100644 --- a/src/dxtbx/model/experiment_list.py +++ b/src/dxtbx/model/experiment_list.py @@ -625,8 +625,7 @@ def from_sequence_and_crystal(imageset, crystal, load_models=True): # if imagesequence is still images, make one experiment for each # all referencing into the same image set if imageset.get_scan().is_still(): - start, end = imageset.get_scan().get_array_range() - for j in range(start, end): + for j in range(len(imageset)): subset = imageset[j : j + 1] experiments.append( Experiment( diff --git a/tests/model/test_beam.py b/tests/model/test_beam.py index 7f1ef5b9d..40659f4d4 100644 --- a/tests/model/test_beam.py +++ b/tests/model/test_beam.py @@ -5,8 +5,8 @@ from libtbx.phil import parse from scitbx import matrix -from dxtbx.model import Beam -from dxtbx.model.beam import BeamFactory, beam_phil_scope +from dxtbx.model import Beam, PolychromaticBeam +from dxtbx.model.beam import BeamFactory, Probe, beam_phil_scope def test_setting_direction_and_wavelength(): @@ -102,6 +102,18 @@ def test_from_phil(): assert b3.get_polarization_fraction() == 0.5 assert b3.get_polarization_normal() == (1.0, 0.0, 0.0) + params3 = beam_phil_scope.fetch( + parse( + """ + beam { + probe = electron + } + """ + ) + ).extract() + b4 = BeamFactory.from_phil(params3, reference) + assert b4.get_probe() == Probe.electron + def test_scan_varying(): direction = matrix.col((0.013142, 0.002200, 1.450476)) @@ -176,3 +188,98 @@ def test_beam_object_comparison(): def test_beam_self_serialization(): beam = Beam() assert beam == BeamFactory.from_dict(beam.to_dict()) + + +def test_polychromatic_beam_from_phil(): + params = beam_phil_scope.fetch( + parse( + """ + beam { + type = polychromatic + direction = (0., 0., 1.) + divergence = 0.2 + sigma_divergence = 0.3 + polarization_normal = (0., -1., 0.) + polarization_fraction = .65 + transmission = .5 + flux = .75 + } + """ + ) + ).extract() + + beam = BeamFactory.from_phil(params) + assert isinstance(beam, PolychromaticBeam) + + assert beam.get_sample_to_source_direction() == pytest.approx((0.0, 0.0, 1.0)) + assert beam.get_divergence() == pytest.approx(0.2) + assert beam.get_sigma_divergence() == pytest.approx(0.3) + assert beam.get_polarization_normal() == pytest.approx((0.0, -1.0, 0.0)) + assert beam.get_polarization_fraction() == pytest.approx(0.65) + assert beam.get_transmission() == pytest.approx(0.5) + assert beam.get_flux() == pytest.approx(0.75) + + +def test_polychromatic_beam_from_dict(): + beam = PolychromaticBeam() + assert beam == BeamFactory.from_dict(beam.to_dict()) + + +def test_make_polychromatic_beam(): + + direction = (0.0, 0.0, 1.0) + divergence = 0.2 + sigma_divergence = 0.3 + polarization_normal = (0.0, -1.0, 0.0) + polarization_fraction = 0.65 + transmission = 0.5 + flux = 0.75 + probe = Probe.neutron + + beam = BeamFactory.make_polychromatic_beam( + direction=direction, + divergence=divergence, + sigma_divergence=sigma_divergence, + polarization_normal=polarization_normal, + polarization_fraction=polarization_fraction, + transmission=transmission, + flux=flux, + probe=probe, + ) + + assert beam.get_sample_to_source_direction() == pytest.approx((0.0, 0.0, 1.0)) + assert beam.get_divergence() == pytest.approx(0.2) + assert beam.get_sigma_divergence() == pytest.approx(0.3) + assert beam.get_polarization_normal() == pytest.approx((0.0, -1.0, 0.0)) + assert beam.get_polarization_fraction() == pytest.approx(0.65) + assert beam.get_transmission() == pytest.approx(0.5) + assert beam.get_flux() == pytest.approx(0.75) + assert beam.get_probe() == Probe.neutron + + +def test_polychromatic_beam_wavelength_guards(): + beam = PolychromaticBeam() + with pytest.raises(RuntimeError): + _ = beam.get_wavelength() + with pytest.raises(RuntimeError): + _ = beam.get_s0() + with pytest.raises(RuntimeError): + _ = beam.get_num_scan_points() + with pytest.raises(RuntimeError): + _ = beam.get_s0_at_scan_points() + with pytest.raises(RuntimeError): + _ = beam.get_s0_at_scan_point(0) + with pytest.raises(RuntimeError): + beam.reset_scan_points() + with pytest.raises(RuntimeError): + beam.set_wavelength(1.0) + with pytest.raises(RuntimeError): + beam.set_s0((0.0, 0.0, 0.1)) + + +def test_polychromatic_beam_str(): + beam = PolychromaticBeam() + assert ( + beam.__str__() + == "Beam:\n probe: x-ray\n sample to source direction : {0,0,1}\n divergence: 0\n sigma divergence: 0\n polarization normal: {0,1,0}\n polarization fraction: 0.5\n flux: 0\n transmission: 1\n" + ) diff --git a/tests/test_imageset.py b/tests/test_imageset.py index 9285db95d..7c54d65b7 100644 --- a/tests/test_imageset.py +++ b/tests/test_imageset.py @@ -387,6 +387,11 @@ def tst_get_item(self, sequence): with pytest.raises(RuntimeError): _ = sequence2[5] + # Check data access matches expected images from the slice + panel_data1 = sequence[3][0] + panel_data2 = sequence2[0][0] + assert panel_data1.all_eq(panel_data2) + assert len(sequence2) == 4 assert_can_get_detectorbase(sequence2, range(0, 4), 5) self.tst_get_models(sequence2, range(0, 4), 5) @@ -396,11 +401,12 @@ def tst_get_item(self, sequence): with pytest.raises(IndexError): _ = sequence[3:7:2] + # Batch offset should not affect slicing of imagesequence # Simulate a scan starting from image 0 sequence_ = copy.deepcopy(sequence) sequence_.get_scan().set_batch_offset(-1) sequence3 = sequence_[3:7] - assert sequence3.get_array_range() == (4, 8) + assert sequence3.get_array_range() == (3, 7) @staticmethod def tst_paths(sequence, filenames1):