<commit_before># LIBTBX_SET_DISPATCHER_NAME dials.background

from __future__ import division

import iotbx.phil

help_message = '''

Examples::

  dials.background image_*.cbf

'''

phil_scope = iotbx.phil.parse("""\
bins = 100
  .type = int
frames = None
  .type = int
  .multiple = True
""", process_includes=True)

def main():
  import sys
  run(sys.argv[1:])

def run(args):

  from dials.util.options import OptionParser
  from dials.util.options import flatten_datablocks
  import libtbx.load_env

  usage = "%s [options] image_*.cbf" % (
    libtbx.env.dispatcher_name)

  parser = OptionParser(
    usage=usage,
    phil=phil_scope,
    read_datablocks=True,
    read_datablocks_from_images=True,
    epilog=help_message)

  params, options = parser.parse_args(show_diff_phil=True)
  datablocks = flatten_datablocks(params.input.datablock)

  if len(datablocks) == 0 and len(experiments) == 0 and len(reflections) == 0:
    parser.print_help()
    exit()

  assert(len(datablocks) == 1)

  datablock = datablocks[0]
  imagesets = datablock.extract_imagesets()

  assert(len(imagesets) == 1)

  imageset = imagesets[0]

  images = imageset.indices()
  if params.frames:
    images = params.frames

  for indx in images:
    print 'For frame %d:' % indx
    background(imageset, indx, params)

def background(imageset, indx, params):
  from dials.array_family import flex
  from libtbx.phil import parse
  from scitbx import matrix
  import math

  detector = imageset.get_detector()
  beam = imageset.get_beam()
  assert(len(detector) == 1)
  detector = detector[0]
  trusted = detector.get_trusted_range()

  n = matrix.col(detector.get_normal()).normalize()
  b = matrix.col(beam.get_s0()).normalize()
  wavelength = beam.get_wavelength()

  if math.fabs(b.dot(n)) < 0.95:
    from libtbx.utils import Sorry
    raise Sorry('Detector not perpendicular to beam')

  data = imageset.get_raw_data(indx)
  assert(len(data) == 1)
  data = data[0]
  negative = (data < 0)
  hot = (data > int(round(trusted[1])))
  bad = negative | hot

  from dials.algorithms.spot_finding.factory import SpotFinderFactory
  from dials.algorithms.spot_finding.factory import phil_scope

  data = data.as_double()

  from dxtbx import datablock

  spot_params = phil_scope.fetch(source=parse("")).extract()
  threshold_function = SpotFinderFactory.configure_threshold(
    spot_params, datablock.DataBlock([imageset]))
  peak_pixels = threshold_function.compute_threshold(data, ~bad)
  signal = data.select(peak_pixels.iselection())
  background = data.select((~bad & ~peak_pixels).iselection())

  # print some summary information
  print 'Mean background: %.3f' % (flex.sum(background) / background.size())
  print 'Max/total signal pixels: %.0f / %.0f' % (flex.max(signal),
                                                 flex.sum(signal))
  print 'Peak/background/hot pixels: %d / %d / %d' % (peak_pixels.count(True),
                                                      background.size(),
                                                      hot.count(True))

  # compute histogram of two-theta values, then same weighted
  # by pixel values, finally divide latter by former to get
  # the radial profile out, need to set the number of bins
  # sensibly; flex.histogram does not allow weights so use
  # numpy.histogram to get the same effect... inspired by
  # method in PyFAI

  data = data.as_1d()
  two_theta_array = detector.get_two_theta_array(beam.get_s0())
  two_theta_array.set_selected((bad | peak_pixels).iselection(), 0.0)
  data.set_selected((bad | peak_pixels).iselection(), 0.0)

  # new fangled flex.weighted_histogram :-)
  n_bins = params.bins
  h0 = flex.weighted_histogram(two_theta_array, n_slots=n_bins)
  h1 = flex.weighted_histogram(two_theta_array, data, n_slots=n_bins)
  h2 = flex.weighted_histogram(two_theta_array, data * data, n_slots=n_bins)

  d0 = h0.slots()
  d1 = h1.slots()
  d2 = h2.slots()

  I = d1 / d0
  I2 = d2 / d0

  tt_range = 0, flex.max(two_theta_array)

  print '%8s %8s %8s' % ('d', 'I', 'sig')
  for j in range(0, len(I)):
    tt = tt_range[0] + (j + 0.5) * (1.0 / n_bins) * (tt_range[1] - tt_range[0])
    d = wavelength / (2.0 * math.sin(0.5 * tt))
    print '%8.3f %8.3f %8.3f' % (d, I[j], math.sqrt(I2[j] - I[j]**2))

if __name__ == '__main__':
  main()<commit_msg>no need to recalculate bin tt values as already calculated by histogram code<commit_after># LIBTBX_SET_DISPATCHER_NAME dials.background

from __future__ import division

import iotbx.phil

help_message = '''

Examples::

  dials.background image_*.cbf

'''

phil_scope = iotbx.phil.parse("""\
bins = 100
  .type = int
frames = None
  .type = int
  .multiple = True
""", process_includes=True)

def main():
  import sys
  run(sys.argv[1:])

def run(args):

  from dials.util.options import OptionParser
  from dials.util.options import flatten_datablocks
  import libtbx.load_env

  usage = "%s [options] image_*.cbf" % (
    libtbx.env.dispatcher_name)

  parser = OptionParser(
    usage=usage,
    phil=phil_scope,
    read_datablocks=True,
    read_datablocks_from_images=True,
    epilog=help_message)

  params, options = parser.parse_args(show_diff_phil=True)
  datablocks = flatten_datablocks(params.input.datablock)

  if len(datablocks) == 0 and len(experiments) == 0 and len(reflections) == 0:
    parser.print_help()
    exit()

  assert(len(datablocks) == 1)

  datablock = datablocks[0]
  imagesets = datablock.extract_imagesets()

  assert(len(imagesets) == 1)

  imageset = imagesets[0]

  images = imageset.indices()
  if params.frames:
    images = params.frames

  for indx in images:
    print 'For frame %d:' % indx
    background(imageset, indx, params)

def background(imageset, indx, params):
  from dials.array_family import flex
  from libtbx.phil import parse
  from scitbx import matrix
  import math

  detector = imageset.get_detector()
  beam = imageset.get_beam()
  assert(len(detector) == 1)
  detector = detector[0]
  trusted = detector.get_trusted_range()

  n = matrix.col(detector.get_normal()).normalize()
  b = matrix.col(beam.get_s0()).normalize()
  wavelength = beam.get_wavelength()

  if math.fabs(b.dot(n)) < 0.95:
    from libtbx.utils import Sorry
    raise Sorry('Detector not perpendicular to beam')

  data = imageset.get_raw_data(indx)
  assert(len(data) == 1)
  data = data[0]
  negative = (data < 0)
  hot = (data > int(round(trusted[1])))
  bad = negative | hot

  from dials.algorithms.spot_finding.factory import SpotFinderFactory
  from dials.algorithms.spot_finding.factory import phil_scope

  data = data.as_double()

  from dxtbx import datablock

  spot_params = phil_scope.fetch(source=parse("")).extract()
  threshold_function = SpotFinderFactory.configure_threshold(
    spot_params, datablock.DataBlock([imageset]))
  peak_pixels = threshold_function.compute_threshold(data, ~bad)
  signal = data.select(peak_pixels.iselection())
  background = data.select((~bad & ~peak_pixels).iselection())

  # print some summary information
  print 'Mean background: %.3f' % (flex.sum(background) / background.size())
  print 'Max/total signal pixels: %.0f / %.0f' % (flex.max(signal),
                                                 flex.sum(signal))
  print 'Peak/background/hot pixels: %d / %d / %d' % (peak_pixels.count(True),
                                                      background.size(),
                                                      hot.count(True))

  # compute histogram of two-theta values, then same weighted
  # by pixel values, finally divide latter by former to get
  # the radial profile out, need to set the number of bins
  # sensibly; flex.histogram does not allow weights so use
  # numpy.histogram to get the same effect... inspired by
  # method in PyFAI

  data = data.as_1d()
  two_theta_array = detector.get_two_theta_array(beam.get_s0())
  two_theta_array.set_selected((bad | peak_pixels).iselection(), 0.0)
  data.set_selected((bad | peak_pixels).iselection(), 0.0)

  # new fangled flex.weighted_histogram :-)
  n_bins = params.bins
  h0 = flex.weighted_histogram(two_theta_array, n_slots=n_bins)
  h1 = flex.weighted_histogram(two_theta_array, data, n_slots=n_bins)
  h2 = flex.weighted_histogram(two_theta_array, data * data, n_slots=n_bins)

  d0 = h0.slots()
  d1 = h1.slots()
  d2 = h2.slots()

  I = d1 / d0
  I2 = d2 / d0
  sig = flex.sqrt(I2 - flex.pow2(I))

  tt = h0.slot_centers()

  print '%8s %8s %8s' % ('d', 'I', 'sig')
  for j in range(len(I)):
    d = wavelength / (2.0 * math.sin(0.5 * tt[j]))
    print '%8.3f %8.3f %8.3f' % (d, I[j], sig[j])

if __name__ == '__main__':
  main()
{'fraction_non_alphanumeric': 0.12440469276338716, 'fraction_numerical': 0.014868161226623301, 'mean_word_length': 3.4153846153846152, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '2613324', 'n_tokens_mistral': 3151, 'n_tokens_neox': 2989, 'n_words': 689}