#!/usr/bin/env python
# -*- coding: utf-8 -*-
#
# Generate "peak powder" from CrystFEL stream
#
# Copyright © 2017 Deutsches Elektronen-Synchrotron DESY,
#                  a research centre of the Helmholtz Association.
#
# Author:
#    2017 Thomas White <taw@physics.org>
#

import numpy as np
import h5py
import sys
import re

f = open(sys.argv[1], 'r')
powder = np.zeros((512,1024), dtype=float)
peaks = []

prog1 = re.compile("^\s*([\d\-\.]+)\s+([\d\-\.]+)\s+[\d\-\.]+\s+([\d\-\.]+)\s+\S+$")

while True:

    fline = f.readline()
    if not fline:
        break

    if fline == '----- End chunk -----\n':
        if len(peaks) > 10:
            for p in peaks:
                powder[p[1],p[0]] += p[2]
        peaks = []

    match = prog1.match(fline)
    if match:
        fs = int(float(match.group(1)))
        ss = int(float(match.group(2)))
        intensity = float(match.group(3))
        peaks.append((fs,ss,intensity))

f.close()

fh = h5py.File("summed-peaks.h5", "w")
fh.create_dataset("/data", (512,1024), data=powder)
fh.close()