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