Reading experiment and data

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

        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
        # a rotation image in here i.e. if scan is not None; derive goniometer
        # matrix else set this to identity
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
# Example code for how to load experiments and reflections in the DIALS


from libtbx.phil import parse

import dials.util.options
from dials.util import show_mail_handle_errors

help_message = """

pass experiment.expt indexed.refl
"""

phil_scope = parse(
    """
    png = 'example.png'
      .type = str
      .help = 'Output name for .png'
    """,
    process_includes=True,
)


class Script:
    """A class for running the script."""

    def __init__(self):
        usage = "dials.experiment_data [options] indexed.expt indexed.refl"

        self.parser = dials.util.options.OptionParser(
            usage=usage,
            phil=phil_scope,
            epilog=help_message,
            check_format=False,
            read_experiments=True,
            read_reflections=True,
        )

    def run(self):
        from scitbx import matrix

        params, options = self.parser.parse_args(show_diff_phil=True)

        experiments = dials.util.options.flatten_experiments(params.input.experiments)
        reflections = dials.util.options.flatten_reflections(params.input.reflections)

        if len(experiments) == 0:
            self.parser.print_help()
            return

        if len(reflections) != 1:
            self.parser.print_help()
            return

        reflections = reflections[0]

        print(f"Read {len(reflections)} reflections")

        indexed = reflections.select(reflections.get_flags(reflections.flags.indexed))

        print(f"Kept {len(indexed)} indexed reflections")

        for name in sorted(indexed.keys()):
            print(f"Found column {name}")

        for reflection in indexed[:3]:
            print(reflection)

        # verify that these experiments correspond to exactly one imageset, one
        # detector, one beam (obviously)
        for experiment in experiments[1:]:
            assert experiment.imageset == experiments[0].imageset
            assert experiment.beam == experiments[0].beam
            assert experiment.detector == experiments[0].detector

        # now perform some calculations - the only things different from one
        # experiment to the next will be crystal models
        print("Crystals:")
        for experiment in experiments:
            print(experiment.crystal)
        detector = experiments[0].detector
        beam = experiments[0].beam
        imageset = experiments[0].imageset

        # derived quantities
        wavelength = beam.get_wavelength()
        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
        # a rotation image in here i.e. if scan is not None; derive goniometer
        # matrix else set this to identity

        scan = experiments[0].scan
        goniometer = experiments[0].goniometer

        if scan and goniometer:
            angle = scan.get_angle_from_array_index(
                0.5 * sum(imageset.get_array_range())
            )
            axis = matrix.col(goniometer.get_rotation_axis_datum())
            F = matrix.sqr(goniometer.get_fixed_rotation())
            S = matrix.sqr(goniometer.get_setting_rotation())
            R = S * axis.axis_and_angle_as_r3_rotation_matrix(angle, deg=True) * F
        else:
            R = matrix.sqr((1, 0, 0, 0, 1, 0, 0, 0, 1))

        print(f"Rotation matrix:\n{R}")

        assert len(detector) == 1


if __name__ == "__main__":
    with show_mail_handle_errors():
        script = Script()
        script.run()