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/299.bugfix b/newsfragments/299.bugfix new file mode 100644 index 000000000..5a8abfc77 --- /dev/null +++ b/newsfragments/299.bugfix @@ -0,0 +1 @@ +Panel geometry definitions in PHIL are merged by panel id before constructing panels diff --git a/newsfragments/439.bugfix b/newsfragments/439.bugfix new file mode 100644 index 000000000..a97cc5116 --- /dev/null +++ b/newsfragments/439.bugfix @@ -0,0 +1 @@ +``flumpy``: Fix case where incorrect ``flex.vec2``, ``flex.vec3`` could be generated. 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/645.feature b/newsfragments/645.feature new file mode 100644 index 000000000..fb38a5121 --- /dev/null +++ b/newsfragments/645.feature @@ -0,0 +1 @@ +Add partial support for the Rigaku Oxford Diffraction file format. 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/650.misc b/newsfragments/650.misc new file mode 100644 index 000000000..94f0ed982 --- /dev/null +++ b/newsfragments/650.misc @@ -0,0 +1 @@ +Format classes are now tested against invalid binary data with dials-data, for when dials-regression is not present. diff --git a/src/dxtbx/boost_python/flumpy.cc b/src/dxtbx/boost_python/flumpy.cc index 12e5ebf94..d7ce34b20 100644 --- a/src/dxtbx/boost_python/flumpy.cc +++ b/src/dxtbx/boost_python/flumpy.cc @@ -520,6 +520,13 @@ py::object vec_from_numpy(py::array np_array) { static_assert(VecType::fixed_size == 2 || VecType::fixed_size == 3, "Only vec2/vec3 supported"); + + // Only accept arrays that have a dimension higher than 1 - we want + // numpy.array([1,2,3]) to fail but numpy.array([[1,2,3]]) to work + if (np_array.ndim() == 1) { + throw std::invalid_argument("Array for conversion to vec must be multidimensional"); + } + // Only accept arrays whose last dimension is the size of this object if (np_array.shape(np_array.ndim() - 1) != VecType::fixed_size) { throw std::invalid_argument("Input array last dimension is not size " 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/FormatROD.py b/src/dxtbx/format/FormatROD.py new file mode 100644 index 000000000..417f99f1a --- /dev/null +++ b/src/dxtbx/format/FormatROD.py @@ -0,0 +1,459 @@ +"""Support for the Rigaku Oxford Diffraction image format + +NB: Rigaku datasets may use a non-zero-padded image incremental serial number. +At present, this is not compatible with assumptions in dxtbx. In order to +import these datasets, the images should be renumbered first. + +See https://github.com/cctbx/dxtbx/issues/646 for details.""" + +# Takanori Nakane took David Waterman's code to parse headers from +# https://github.com/cctbx/cctbx_project/commit/b95467f3b2a70a37eeb820ea294128a32551700c +# and heavily modified it. The original commit in the cctbx_project repository is orphan now. +from __future__ import annotations + +__author__ = "David Waterman, Takanori Nakane" +__copyright__ = ( + "Copyright 2018 United Kingdom Research and Innovation & 2022 Takanori Nakane" +) +__license__ = "BSD 3-clause" + +import re +import struct + +import numpy as np + +from scitbx.array_family import flex + +from dxtbx.format.Format import Format + + +class FormatROD(Format): + """An image reading class for Rigaku Oxford Diffraction images. + + The Rigaku Oxford Diffraction format is used by CrysAlis(Pro). + """ + + def __init__(self, image_file, **kwargs): + """Initialise the image structure from the given file.""" + + from dxtbx import IncorrectFormatError + + if not self.understand(image_file): + raise IncorrectFormatError(self, image_file) + + Format.__init__(self, image_file, **kwargs) + + @staticmethod + def understand(image_file): + with FormatROD.open_file(image_file, "rb") as f: + try: + hdr = f.read(256).decode("ascii") + except UnicodeDecodeError: + return False + + lines = hdr.splitlines() + if len(lines) < 2: + return False + + vers = lines[0].split() + if len(vers) < 2 or vers[0] != "OD" or vers[1] != "SAPPHIRE": + return False + + compression = lines[1].split("=") + if compression[0] != "COMPRESSION": + return False + + return True + + @staticmethod + def _read_ascii_header(image_file): + """Read the ASCII header comprising the first 256 bytes of the file""" + + hd = {} + with FormatROD.open_file(image_file, "rb") as f: + hdr = f.read(256).decode("ascii") + lines = hdr.splitlines() + + vers = lines[0].split() + if len(vers) < 2 or vers[0] != "OD" or vers[1] != "SAPPHIRE": + raise ValueError("Wrong header format") + hd["version"] = float(vers[-1]) + + compression = lines[1].split("=") + if compression[0] != "COMPRESSION": + raise ValueError("Wrong header format") + hd["compression"] = compression[1] + + # Extract definitions from the 3rd - 5th line + defn = re.compile(r"([A-Z]+=[ 0-9]+)") + for line in lines[2:5]: + sizes = defn.findall(line) + for s in sizes: + n, v = s.split("=") + hd[n] = int(v) + + hd["time"] = lines[5].split("TIME=")[-1].strip("\x1a").rstrip() + + return hd + + @staticmethod + def _read_binary_header( + image_file, + offset=256, + nbytes=5120, + general_nbytes=512, + special_nbytes=768, + km4gonio_nbytes=1024, + statistics_nbytes=512, + history_nbytes=2048, + ): + """Read the most relevant parameters from the binary header""" + + with FormatROD.open_file(image_file, "rb") as f: + # General section + f.seek(offset) + bin_x, bin_y = struct.unpack("> 4) & 15) + ipos += 1 + # ipos_bit = ipos * 8 + + for i in range(2): + nbit = nbits[i] + # 0: 0 + # 1: 1 or 0 + # 2: -1, 0, 1, 2 + # 3: -3, -2, -1, 0, 1, 2, 3, 4 (= b111 -3 = 7 - 3) + zero_at = 0 + if nbit > 1: + zero_at = (1 << (nbit - 1)) - 1 + + v = 0 + for j in range(nbit): + v |= linedata[ipos] << (8 * j) + ipos += 1 + + mask = (1 << nbit) - 1 + for j in range(BLOCKSIZE): + ret[opos] = ((v >> (nbit * j)) & mask) - zero_at + opos += 1 + + for i in range(opos - BLOCKSIZE * 2, opos): + offset = ret[i] + + # surprisingly, this IF is very slow! + if offset >= SHORT_OVERFLOW_SIGNED: + if offset >= LONG_OVERFLOW_SIGNED: + offset = linedata[ipos : (ipos + 4)].view(np.int32)[0] + ipos += 4 + else: + offset = linedata[ipos : (ipos + 2)].view(np.int16)[0] + ipos += 2 + + ret[i] = offset + ret[i - 1] + + for i in range(nrest): + px = linedata[ipos] + ipos += 1 + if px < SHORT_OVERFLOW: + ret[opos] = ret[opos - 1] + px - 127 + elif px == LONG_OVERFLOW: + ret[opos] = ( + ret[opos - 1] + linedata[ipos : (ipos + 4)].view(np.int32)[0] + ) + ipos += 4 + else: + ret[opos] = ( + ret[opos - 1] + linedata[ipos : (ipos + 2)].view(np.int16)[0] + ) + ipos += 2 + opos += 1 + + return ret + + +if __name__ == "__main__": + import sys + + for arg in sys.argv[1:]: + if FormatROD.understand(arg): + reader = FormatROD(arg) + print(reader._txt_header) + print(reader._bin_header) + 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 574779956..85fd2200b 100644 --- a/src/dxtbx/imageset.py +++ b/src/dxtbx/imageset.py @@ -300,29 +300,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..fe2fda864 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,50 @@ 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: + return std::string("unknown"); + } + } + + 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); + return Probe::xray; + } + + 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 +378,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 +416,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 +431,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 +439,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 +463,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/detector.py b/src/dxtbx/model/detector.py index c98828c72..c046ca8d9 100644 --- a/src/dxtbx/model/detector.py +++ b/src/dxtbx/model/detector.py @@ -1,6 +1,7 @@ from __future__ import annotations import os +from copy import deepcopy import pycbf @@ -212,6 +213,42 @@ ) +def merge_panel_scope_extracts_by_id(panel_params): + + id_to_params = {} + for (i, params) in enumerate(panel_params): + if params.id not in id_to_params: + id_to_params[params.id] = [ + i, + ] + else: + id_to_params[params.id].append(i) + + merged_params = [] + for params_set in id_to_params.values(): + params0 = deepcopy(panel_params[params_set[0]]) + for i in params_set[1:]: + params1 = deepcopy(panel_params[i]) + for key in params1.__dict__: + if ( + key.startswith("_") + or key == "id" + or params0.__dict__[key] == params1.__dict__[key] + ): + continue + if ( + params0.__dict__[key] is not None + and params1.__dict__[key] is not None + ): + raise RuntimeError( + f"Multiple definitions for {key} for panel id={params0.id}" + ) + if params0.__dict__[key] is None: + params0.__dict__[key] = params1.__dict__[key] + merged_params.append(params0) + return merged_params + + class DetectorFactory: """A factory class for detector objects, which will encapsulate standard detector designs to make it a little easier to get started with these. In @@ -229,7 +266,11 @@ def generate_from_phil(params, beam=None): # Create a list of panels panel_list = {} - for panel_params in params.detector.panel: + + # merge panel params by id first + merged = merge_panel_scope_extracts_by_id(params.detector.panel) + + for panel_params in merged: panel = Panel() if panel_params.name is not None: panel.set_name(panel_params.name) @@ -335,8 +376,11 @@ def overwrite_from_phil(params, detector, beam=None): Overwrite from phil parameters """ + # merge panel params by id first + merged = merge_panel_scope_extracts_by_id(params.detector.panel) + # Override any panel parameters - for panel_params in params.detector.panel: + for panel_params in merged: panel = detector[panel_params.id] if panel_params.name is not None: panel.set_name(panel_params.name) 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/format/test_format.py b/tests/format/test_format.py new file mode 100644 index 000000000..3b951ba8d --- /dev/null +++ b/tests/format/test_format.py @@ -0,0 +1,12 @@ +from __future__ import annotations + +from dxtbx.format import Registry + + +def test_reading_refl_failure(dials_data): + test_data = dials_data("centroid_test_data", pathlib=True) + + # Without dials_regression, none of the dials-data tests check for this "invalid binary data" case + assert Registry.get_format_class_for_file(test_data / "indexed.refl") is None + # Check .expt while here + assert Registry.get_format_class_for_file(test_data / "indexed.expt") is None 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_flumpy.py b/tests/test_flumpy.py index cd09e3c09..88ef9f8cc 100644 --- a/tests/test_flumpy.py +++ b/tests/test_flumpy.py @@ -136,7 +136,7 @@ def test_reverse_numeric_2d(flex_numeric): def test_numeric_4d(flex_numeric): - #  Check that we can think fourth-dimnesionally + # Check that we can think fourth-dimensionally grid = flex.grid(1, 9, 8, 5) fo = flex_numeric(grid) assert fo.nd() == 4 @@ -402,3 +402,11 @@ def test_int_long_degeneracy(): npo[0] = 42 assert all(fo[x] == npo[x] for x in range(4)) assert fo[0] == 42 + + +def test_single_entry_vec_from_numpy(): + with pytest.raises(ValueError): + flumpy.vec_from_numpy(np.array([2.0, 3.0, 4.0])) + ao = flumpy.vec_from_numpy(np.array([[2.0, 3.0, 4.0]])) + assert len(ao) == 1 + assert ao[0] == pytest.approx([2.0, 3.0, 4.0]) 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):