Reading experiment and dataΒΆ

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

    # now perform some calculations - the only things different from one
    # experiment to the next will be crystal models
    crystals = [experiment.crystal for experiment in experiments]
    detector = experiments[0].detector
    beam = experiments[0].beam
    imageset = experiments[0].imageset
  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
119
120
121
122
123
124
125
126
#!/usr/bin/env dials.python
#
# example_experiment_data.py
#
#  Copyright (C) 2017 Diamond Light Source
#
#  Author: Graeme Winter
#
#  This code is distributed under the BSD license, a copy of which is
#  included in the root directory of this package.
#
# Example code for how to load experiments and reflections in the DIALS
# framework

from __future__ import division
from libtbx.phil import parse

help_message = '''

pass experiment.json indexed.pickle

'''

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

class Script(object):
  '''A class for running the script.'''

  def __init__(self):
    from dials.util.options import OptionParser
    import libtbx.load_env

    usage = 'usage: %s [options] experiments.json indexed.pickle' \
            % libtbx.env.dispatcher_name

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

  def run(self):
    from dials.array_family import flex # import dependency
    from scitbx import matrix
    from dials.util.options import flatten_experiments
    from dials.util.options import flatten_reflections

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

    experiments = flatten_experiments(params.input.experiments)
    reflections = 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 'Read %d reflections' % len(reflections)

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

    print 'Kept %d indexed reflections' % len(indexed)

    for name in sorted(indexed.keys()):
      print 'Found column %s' % 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
    crystals = [experiment.crystal for experiment in experiments]
    detector = experiments[0].detector
    beam = experiments[0].beam
    imageset = experiments[0].imageset

    # derived quantities
    wavelength = beam.get_wavelength()
    s0 = matrix.col(beam.get_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))

    assert(len(detector) == 1)


if __name__ == '__main__':
  from dials.util import halraiser
  try:
    script = Script()
    script.run()
  except Exception as e:
    halraiser(e)