Reading experiment and data

The following example code demonstrates how to load experiments and reflections in the DIALS framework:

        s0 = matrix.col(beam.get_s0())

        print(f"Wavelength: {wavelength:g}Å")
        print(f"Beam vector s0:\n{s0}")

        # in here do some jiggery-pokery to cope with this being interpreted as

  1# Example code for how to load experiments and reflections in the DIALS
  2
  3
  4from __future__ import annotations
  5
  6from libtbx.phil import parse
  7
  8import dials.util.options
  9from dials.util import show_mail_handle_errors
 10
 11help_message = """
 12
 13pass experiment.expt indexed.refl
 14"""
 15
 16phil_scope = parse(
 17    """
 18    png = 'example.png'
 19      .type = str
 20      .help = 'Output name for .png'
 21    """,
 22    process_includes=True,
 23)
 24
 25
 26class Script:
 27    """A class for running the script."""
 28
 29    def __init__(self):
 30        usage = "dials.experiment_data [options] indexed.expt indexed.refl"
 31
 32        self.parser = dials.util.options.ArgumentParser(
 33            usage=usage,
 34            phil=phil_scope,
 35            epilog=help_message,
 36            check_format=False,
 37            read_experiments=True,
 38            read_reflections=True,
 39        )
 40
 41    def run(self):
 42        from scitbx import matrix
 43
 44        params, options = self.parser.parse_args(show_diff_phil=True)
 45
 46        experiments = dials.util.options.flatten_experiments(params.input.experiments)
 47        reflections = dials.util.options.flatten_reflections(params.input.reflections)
 48
 49        if len(experiments) == 0:
 50            self.parser.print_help()
 51            return
 52
 53        if len(reflections) != 1:
 54            self.parser.print_help()
 55            return
 56
 57        reflections = reflections[0]
 58
 59        print(f"Read {len(reflections)} reflections")
 60
 61        indexed = reflections.select(reflections.get_flags(reflections.flags.indexed))
 62
 63        print(f"Kept {len(indexed)} indexed reflections")
 64
 65        for name in sorted(indexed.keys()):
 66            print(f"Found column {name}")
 67
 68        for reflection in indexed[:3]:
 69            print(reflection)
 70
 71        # verify that these experiments correspond to exactly one imageset, one
 72        # detector, one beam (obviously)
 73        for experiment in experiments[1:]:
 74            assert experiment.imageset == experiments[0].imageset
 75            assert experiment.beam == experiments[0].beam
 76            assert experiment.detector == experiments[0].detector
 77
 78        # now perform some calculations - the only things different from one
 79        # experiment to the next will be crystal models
 80        print("Crystals:")
 81        for experiment in experiments:
 82            print(experiment.crystal)
 83        detector = experiments[0].detector
 84        beam = experiments[0].beam
 85        imageset = experiments[0].imageset
 86
 87        # derived quantities
 88        wavelength = beam.get_wavelength()
 89        s0 = matrix.col(beam.get_s0())
 90
 91        print(f"Wavelength: {wavelength:g}Å")
 92        print(f"Beam vector s0:\n{s0}")
 93
 94        # in here do some jiggery-pokery to cope with this being interpreted as
 95        # a rotation image in here i.e. if scan is not None; derive goniometer
 96        # matrix else set this to identity
 97
 98        scan = experiments[0].scan
 99        goniometer = experiments[0].goniometer
100
101        if scan and goniometer:
102            angle = scan.get_angle_from_array_index(
103                0.5 * sum(imageset.get_array_range())
104            )
105            axis = matrix.col(goniometer.get_rotation_axis_datum())
106            F = matrix.sqr(goniometer.get_fixed_rotation())
107            S = matrix.sqr(goniometer.get_setting_rotation())
108            R = S * axis.axis_and_angle_as_r3_rotation_matrix(angle, deg=True) * F
109        else:
110            R = matrix.sqr((1, 0, 0, 0, 1, 0, 0, 0, 1))
111
112        print(f"Rotation matrix:\n{R}")
113
114        assert len(detector) == 1
115
116
117if __name__ == "__main__":
118    with show_mail_handle_errors():
119        script = Script()
120        script.run()