Merge branch 'detwin'

16 files changed, 1323 insertions, 21 deletions

diff --git a/.gitignore b/.gitignore
index 40000069..5091a898 100644
--- a/.gitignore
+++ b/.gitignore
@@ -34,6 +34,7 @@ src/partialator
 src/check_hkl
 src/partial_sim
 src/cell_explorer
+src/ambigator
 src/.dirstamp
 *~
 doc/reference/libcrystfel/*
diff --git a/AUTHORS b/AUTHORS
index 91876b85..4b24d960 100644
--- a/AUTHORS
+++ b/AUTHORS
@@ -42,3 +42,6 @@
 
 * Alexandra Tolstikova
   SASE spectrum simulation
+
+* Wolfgang Brehm <wolfgang.brehm@gmail.com>
+  "Detwinning" algorithm (ambigator)
diff --git a/Makefile.am b/Makefile.am
index ad39c1ee..cc5c6485 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -4,7 +4,8 @@ SUBDIRS = lib doc/reference/libcrystfel libcrystfel
 ACLOCAL_AMFLAGS = -I m4
 
 bin_PROGRAMS = src/pattern_sim src/process_hkl src/get_hkl src/indexamajig \
-               src/compare_hkl src/partialator src/check_hkl src/partial_sim
+               src/compare_hkl src/partialator src/check_hkl src/partial_sim \
+               src/ambigator
 
 noinst_PROGRAMS = tests/list_check tests/integration_check \
                   tests/pr_p_gradient_check tests/symmetry_check \
@@ -85,6 +86,8 @@ endif
 src_partialator_SOURCES = src/partialator.c src/post-refinement.c \
                           src/hrs-scaling.c
 
+src_ambigator_SOURCES = src/ambigator.c
+
 if HAVE_CAIRO
 if HAVE_PANGO
 if HAVE_PANGOCAIRO
@@ -137,7 +140,7 @@ man_MANS = doc/man/crystfel.7 doc/man/check_hkl.1 doc/man/compare_hkl.1 \
            doc/man/crystfel_geometry.5 doc/man/get_hkl.1 doc/man/hdfsee.1 \
            doc/man/indexamajig.1 doc/man/partialator.1 doc/man/partial_sim.1 \
            doc/man/pattern_sim.1 doc/man/process_hkl.1 doc/man/render_hkl.1 \
-           doc/man/cell_explorer.1
+           doc/man/cell_explorer.1 doc/man/ambigator.1
 
 
 EXTRA_DIST += $(man_MANS)
@@ -173,7 +176,7 @@ script_DATA = scripts/alternate-stream scripts/cell-please \
               scripts/random-image scripts/README scripts/sequence-image \
               scripts/split-indexed scripts/stream_grep scripts/zone-axes \
               scripts/Rsplit_surface scripts/Rsplit_surface.py \
-              scripts/clean-stream.py
+              scripts/clean-stream.py scripts/fg-graph
 
 EXTRA_DIST += $(script_DATA)
 
diff --git a/README b/README
index e5d3b5db..83f44b38 100644
--- a/README
+++ b/README
@@ -19,6 +19,7 @@ Authors:
   Cornelius Gati <cornelius.gati@desy.de>
   Fedor Chervinskii <fedor.chervinskii@gmail.com>
   Alexandra Tolstikova
+  Wolfgang Brehm <wolfgang.brehm@gmail.com>
 
 Please read the AUTHORS file for a full list of contributions and contributors.
 See the COPYING file for licence conditions.  Summary: GPLv3+.
diff --git a/doc/man/ambigator.1 b/doc/man/ambigator.1
new file mode 100644
index 00000000..b0978fb7
--- /dev/null
+++ b/doc/man/ambigator.1
@@ -0,0 +1,103 @@
+.\"
+.\" ambigator man page
+.\"
+.\" Copyright © 2014 Deutsches Elektronen-Synchrotron DESY,
+.\"                  a research centre of the Helmholtz Association.
+.\"
+.\" Part of CrystFEL - crystallography with a FEL
+.\"
+
+.TH AMBIGATOR 1
+.SH NAME
+ambigator \- Resolve indexing ambiguities
+.SH SYNOPSIS
+.PP
+.B ambigator \fIinput.stream\fR \fB[-o\fR \fIoutput.stream\fR\fB] [options]
+
+.B ambigator --help
+
+.SH DESCRIPTION
+This program resolves indexing ambiguities using a simplified variant of the clustering algorithm described by (FIXME: Reference).
+
+The algorithm works by calculating the individual correlation coefficients between the intensities from one crystal and those from each of the other crystals in turn.  The mean \fIf\fR is then taken over all crystals with the same indexing assignment, and separately the mean \fIg\fR is taken over all crystals with the opposite indexing assignment.  The indexing assignment for a crystal is changed if \fIg\fR > \fIf\fR.  Every crystal is visited once in turn, and the pass over all the crystals repeated several times.
+
+Only one indexing ambiguity can be resolved at once.  In other words, each crystal is considered to be indexable in one of only two ways.  If the true indexing ambiguity has more possibilities than this, the resolution must be performed by running \fBambigator\fR multiple times with a different ambiguity operator each time.
+
+If the ambiguity operator is known (or, equivalently, the actual and apparent symmetries are both known), then the algorithm can be enhanced by including in \fIf\fR the correlation coefficients of all the patterns with the opposite indexing assignment to the current one, but after reindexing the other pattern first.  Likewise, \fg\fR includes the correlation coefficients of the patterns with the same indexing assignment after reindexing.  This enhances the algorithm to an extent roughly equivalent to doubling the number of crystals.
+
+The default behaviour is to compare each crystal to every other crystal.  This leads to a computational complexity proportional to the square of the number of crystals.  If the number of crystals is large, the number of comparisons can be limited without compromising the algorithm much.  In this case, the crystals to correlate against will be selected randomly.
+
+
+.SH OPTIONS
+.PD 0
+.IP "\fB-o\fR \fIfilename\fR"
+.IP \fB--output=\fR\fIfilename\fR
+.PD
+Write a re-indexed version of the input stream to \fIfilename\fR.  This stream can then be merged normally using \fBprocess_hkl\fR or \fBpartialator\fR, but using the actual symmetry instead of the apparent one.
+.IP
+\fBWARNING\fR: There is no default filename.  The default behaviour is not to output any reindexed stream!
+
+.PD 0
+.IP "\fB-y\fR \fIpg\fR"
+.IP \fB--symmetry=\fR\fIpg\fR
+.PD
+Set the actual symmetry of the crystals.  If you're not sure, set this to the highest symmetry which you want to assume, which might be \fB-1\fR to assume Friedel's Law alone or \fB1\fR (the default) for no symmetry at all.  The algorithm will work significantly better if you can use a higher symmetry here.
+
+.PD 0
+.IP "\fB-w\fR \fIpg\fR"
+.PD
+Set the apparent symmetry of the crystals.  The ambiguity operator will be determined by comparing this to the actual symmetry.  If there is more than one operator, only the first will be used.
+.IP
+Using this option improves the algorithm to an extent roughly equivalent to doubling the number of crystals.
+
+.PD 0
+.IP "\fB-n\fR \fIn\fR"
+.IP \fB--iterations=\fR\fIn\fR
+The number of passes through the data to make.  Extra iterations are not expensive once the initial correlation calculation has been performed, so set this value quite high.  Two or three iterations are normally sufficient unless the number of correlations (see \fB--ncorr\fR) is small compared to the number of crystals.  The default is \fB--iterations=6\fR.
+
+.PD 0
+.IP "\fB-j\fR \fIn\fR"
+Number of threads to use for the CC calculation.
+
+.PD 0
+.IP \fB--highres=\fR\fId\fR
+High resolution cutoff in Angstroms.
+
+.PD 0
+.IP \fB--lowres=\fR\fId\fR
+Low resolution cutoff in Angstroms.
+
+.PD 0
+.IP \fB--end-assignments=\fR\fIfilename\fR
+Write the end assignments to \fIfilename\fR.  The file will be a list of 0 or 1, one value per line, in the same order as the crystals appear in the input stream.  1 means that the pattern should be reindexed according to the ambiguity operator.
+
+.PD 0
+.IP \fB--fg-graph=\fR\fIfilename\fR
+Write f and g values to \fIfilename\fR, one line per crystal, repeating all crystals as they are visited by the algorithm.  Plot these using \fBfg-graph\fR from the CrystFEL script folder to evaluate the ambiguity resolution.
+
+.PD 0
+.IP \fB--ncorr=\fR\fIn\fR
+Use \fIn\fR correlations per crystal.  The default is to correlate against every crystal.  If the CC calculation is too slow, try \fB--ncorr=1000\fR.  Note that this option sets the maximum number of correlations, and some crystals might not have enough common reflections to correlate to the number requested.  The mean number of actual correlations per crystal will be output by the program after the CC calculation, and if this number is much smaller than \fIn\fR then this option will not have a significant effect.
+
+
+.SH AUTHOR
+This page was written by Thomas White.
+
+.SH REPORTING BUGS
+Report bugs to <taw@physics.org>, or visit <http://www.desy.de/~twhite/crystfel>.
+
+.SH COPYRIGHT AND DISCLAIMER
+Copyright © 2014 Deutsches Elektronen-Synchrotron DESY, a research centre of the Helmholtz Association.
+
+ambigator, and this manual, are part of CrystFEL.
+
+CrystFEL is free software: you can redistribute it and/or modify it under the
+terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
+
+CrystFEL is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
+
+You should have received a copy of the GNU General Public License along with CrystFEL.  If not, see <http://www.gnu.org/licenses/>.
+
+.SH SEE ALSO
+.BR crystfel (7),
+.BR indexamajig (1).
diff --git a/libcrystfel/src/stream.c b/libcrystfel/src/stream.c
index db8a4f72..80256829 100644
--- a/libcrystfel/src/stream.c
+++ b/libcrystfel/src/stream.c
@@ -133,6 +133,61 @@ static void write_peaks(struct image *image, FILE *ofh)
 }
 
 
+static RefList *read_stream_reflections_2_1(FILE *fh)
+{
+	char *rval = NULL;
+	int first = 1;
+	RefList *out;
+
+	out = reflist_new();
+
+	do {
+
+		char line[1024];
+		signed int h, k, l;
+		float intensity, sigma, fs, ss;
+		char phs[1024];
+		int cts;
+		int r;
+		Reflection *refl;
+
+		rval = fgets(line, 1023, fh);
+		if ( rval == NULL ) continue;
+		chomp(line);
+
+		if ( strcmp(line, REFLECTION_END_MARKER) == 0 ) return out;
+
+		r = sscanf(line, "%i %i %i %f %s %f %i %f %f",
+		           &h, &k, &l, &intensity, phs, &sigma, &cts, &fs, &ss);
+		if ( (r != 9) && (!first) ) {
+			reflist_free(out);
+			return NULL;
+		}
+
+		first = 0;
+		if ( r == 9 ) {
+
+			double ph;
+			char *v;
+
+			refl = add_refl(out, h, k, l);
+			set_intensity(refl, intensity);
+			set_detector_pos(refl, 0.0, fs, ss);
+			set_esd_intensity(refl, sigma);
+			set_redundancy(refl, cts);
+
+			ph = strtod(phs, &v);
+			if ( v != phs ) set_phase(refl, deg2rad(ph));
+
+		}
+
+	} while ( rval != NULL );
+
+	/* Got read error of some kind before finding PEAK_LIST_END_MARKER */
+	return NULL;
+}
+
+
 static RefList *read_stream_reflections(FILE *fh)
 {
 	char *rval = NULL;
@@ -214,7 +269,50 @@ static void write_stream_reflections(FILE *fh, RefList *list)
 		       h, k, l, intensity, esd_i, pk, bg, fs, ss);
 
 	}
+}
+
+
+static void write_stream_reflections_2_1(FILE *fh, RefList *list)
+{
+	Reflection *refl;
+	RefListIterator *iter;
+
+	fprintf(fh, "  h   k   l          I    phase   sigma(I) "
+		     " counts  fs/px  ss/px\n");
+
+	for ( refl = first_refl(list, &iter);
+	      refl != NULL;
+	      refl = next_refl(refl, iter) )
+	{
+
+		signed int h, k, l;
+		double intensity, esd_i, ph;
+		int red;
+		double fs, ss;
+		char phs[16];
+		int have_phase;
+
+		get_indices(refl, &h, &k, &l);
+		get_detector_pos(refl, &fs, &ss);
+		intensity = get_intensity(refl);
+		esd_i = get_esd_intensity(refl);
+		red = get_redundancy(refl);
+		ph = get_phase(refl, &have_phase);
+
+		/* Reflections with redundancy = 0 are not written */
+		if ( red == 0 ) continue;
 
+		if ( have_phase ) {
+			snprintf(phs, 16, "%8.2f", rad2deg(ph));
+		} else {
+			strncpy(phs, "       -", 15);
+		}
+
+		fprintf(fh,
+		       "%3i %3i %3i %10.2f %s %10.2f %7i %6.1f %6.1f\n",
+		       h, k, l, intensity, phs, esd_i, red,  fs, ss);
+
+	}
 }
 
 
@@ -276,7 +374,9 @@ static void write_crystal(Stream *st, Crystal *cr, int include_reflections)
 			if ( AT_LEAST_VERSION(st, 2, 2) ) {
 				write_stream_reflections(st->fh, reflist);
 			} else {
-				write_reflections_to_file(st->fh, reflist);
+				/* This function writes like a normal reflection
+				 * list was written in stream 2.1 */
+				write_stream_reflections_2_1(st->fh, reflist);
 			}
 			fprintf(st->fh, REFLECTION_END_MARKER"\n");
 
@@ -465,7 +565,7 @@ void read_crystal(Stream *st, struct image *image)
 			if ( AT_LEAST_VERSION(st, 2, 2) ) {
 				reflist = read_stream_reflections(st->fh);
 			} else {
-				reflist = read_reflections_from_file(st->fh);
+				reflist = read_stream_reflections_2_1(st->fh);
 			}
 			if ( reflist == NULL ) {
 				ERROR("Failed while reading reflections\n");
diff --git a/libcrystfel/src/symmetry.c b/libcrystfel/src/symmetry.c
index 81d87b25..dec8a520 100644
--- a/libcrystfel/src/symmetry.c
+++ b/libcrystfel/src/symmetry.c
@@ -1110,8 +1110,6 @@ void get_equiv(const SymOpList *ops, const SymOpMask *m, int idx,
 
 	}
 
-
-
 	if ( idx >= n ) {
 
 		ERROR("Index %i out of range for point group '%s'\n", idx,
diff --git a/scripts/fg-graph b/scripts/fg-graph
new file mode 100755
index 00000000..30ca522a
--- /dev/null
+++ b/scripts/fg-graph
@@ -0,0 +1,61 @@
+#!/bin/bash
+
+# Copyright © 2014 Wolfgang Brehm
+# Copyright © 2014 Deutsches Elektronen-Synchrotron DESY,
+#                  a research centre of the Helmholtz Association.
+#
+# Authors:
+#   2013-2014 Wolfgang Brehm <wolfgang.brehm@gmail.com>
+#   2014 Thomas White <taw@physics.org>
+#
+# This file is part of CrystFEL.
+#
+# CrystFEL is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# CrystFEL is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with CrystFEL.  If not, see <http://www.gnu.org/licenses/>.
+#
+
+FILE=fg.dat
+NPATT=4000
+NITER=4
+CMIN=-0.1
+CMAX=0.3
+
+gnuplot -p << EOF
+set terminal pngcairo size 1600,1000 enhanced monochrome font 'Helvetica,20' linewidth 2
+set output "correlation.png"
+set xlabel "Number of crystals"
+set ylabel "Correlation"
+rnd(x) = x - floor(x) < 0.5 ? floor(x) : ceil(x)
+round(x1,x2)=rnd(10**x2*x1)/10.0**x2
+set xzeroaxis lc rgb "black" lt 1
+set key top left
+set xrange [0:${NITER}*${NPATT}]
+set yrange [${CMIN}:${CMAX}]
+set ytics nomirror
+set xtics nomirror
+set x2tics nomirror
+
+set x2range [0:${NITER}]
+set x2label "Number of passes over all crystals"
+set arrow from ${NPATT},${CMIN} to ${NPATT},${CMAX} nohead lc rgb "black"
+
+plot\
+ 1 lc rgb "black" lt 1 notitle,\
+ "${FILE}" u 0:1 ps 0.2 pt 7 lc rgb "orange" notitle,\
+ "${FILE}" u 0:2 ps 0.2 pt 7 lc rgb "#0088FF" notitle,\
+ "${FILE}" u 0:(rand(0)>0.5?\$1:1/0) ps 0.2 pt 7 lc rgb "orange" notitle ,\
+ "${FILE}" u (round(\$0, -3)):(\$1>\$2?\$1:\$2) lw 3 lc rgb "black" smooth unique notitle,\
+ "${FILE}" u (round(\$0, -3)):(\$1<\$2?\$1:\$2) lw 3 lc rgb "black" smooth unique notitle,\
+ "${FILE}" u (round(\$0, -3)):1 lw 3 lc rgb "red" smooth unique notitle,\
+ "${FILE}" u (round(\$0, -3)):2 lw 3 lc rgb "blue" smooth unique notitle
+EOF
diff --git a/src/ambigator.c b/src/ambigator.c
new file mode 100644
index 00000000..da7e0e11
--- /dev/null
+++ b/src/ambigator.c
@@ -0,0 +1,1018 @@
+/*
+ * ambigator.c
+ *
+ * Resolve indexing ambiguities
+ *
+ * Copyright © 2014 Deutsches Elektronen-Synchrotron DESY,
+ *                  a research centre of the Helmholtz Association.
+ * Copyright © 2014 Wolfgang Brehm
+ *
+ * Authors:
+ *   2014 Thomas White <taw@physics.org>
+ *   2014 Wolfgang Brehm <wolfgang.brehm@gmail.com>
+ *
+ * This file is part of CrystFEL.
+ *
+ * CrystFEL is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * CrystFEL is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with CrystFEL.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ */
+
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdarg.h>
+#include <stdlib.h>
+#include <stdio.h>
+#include <string.h>
+#include <unistd.h>
+#include <getopt.h>
+#include <assert.h>
+#include <gsl/gsl_rng.h>
+#include <gsl/gsl_permutation.h>
+#include <gsl/gsl_randist.h>
+
+#include <image.h>
+#include <utils.h>
+#include <symmetry.h>
+#include <stream.h>
+#include <reflist.h>
+#include <reflist-utils.h>
+#include <cell.h>
+#include <cell-utils.h>
+#include <thread-pool.h>
+
+
+static void show_help(const char *s)
+{
+	printf("Syntax: %s [options] input.stream\n\n", s);
+	printf(
+"Resolve indexing ambiguities.\n"
+"\n"
+"  -h, --help                 Display this help message.\n"
+"\n"
+"  -o, --output=<filename>    Output stream.\n"
+"  -y, --symmetry=<sym>       Actual (\"target\") symmetry.\n"
+"  -w <sym>                   Apparent (\"source\" or \"twinned\") symmetry.\n"
+"  -n, --iterations=<n>       Iterate <n> times.\n"
+"      --highres=<n>          High resolution cutoff in A.\n"
+"      --lowres=<n>           Low resolution cutoff in A.\n"
+"      --end-assignments=<fn> Save end assignments to file <fn>.\n"
+"      --fg-graph=<fn>        Save f and g correlation values to file <fn>.\n"
+"      --ncorr=<n>            Use <n> correlations per crystal.  Default 1000\n"
+"      --stop-after=<n>       Use at most the first <n> crystals.\n"
+"  -j <n>                     Use <n> threads for CC calculation.\n"
+);
+}
+
+#define SERIAL(h, k, l) ((((h)+512)<<20) + (((k)+512)<<10) + ((l)+512))
+#define GET_H(serial) ((((serial) & 0x3ff00000)>>20)-512)
+#define GET_K(serial) ((((serial) & 0x000ffc00)>>10)-512)
+#define GET_L(serial) (((serial) & 0x000003ff)-512)
+
+struct flist
+{
+	int n;
+
+	unsigned int *s;
+	float *i;
+
+	unsigned int *s_reidx;
+	float *i_reidx;
+};
+
+
+static struct flist *asymm_and_merge(RefList *in, const SymOpList *sym,
+                                     UnitCell *cell, double rmin, double rmax,
+                                     SymOpList *amb)
+{
+	Reflection *refl;
+	RefListIterator *iter;
+	RefList *asym;
+	struct flist *f;
+	int n;
+
+	asym = reflist_new();
+	if ( asym == NULL ) return NULL;
+
+	for ( refl = first_refl(in, &iter);
+	      refl != NULL;
+	      refl = next_refl(refl, iter) )
+	{
+		signed int h, k, l;
+		signed int ha, ka, la;
+		Reflection *cr;
+		double res;
+
+		get_indices(refl, &h, &k, &l);
+
+		res = 2.0*resolution(cell, h, k, l);
+		if ( res < rmin ) continue;
+		if ( res > rmax ) continue;
+
+		get_asymm(sym, h, k, l, &ha, &ka, &la);
+
+		if ( amb != NULL ) {
+
+			signed int hr, kr, lr;
+			signed int hra, kra, lra;
+
+			get_equiv(amb, NULL, 1, ha, ka, la, &hr, &kr, &lr);
+			get_asymm(sym, hr, kr, lr, &hra, &kra, &lra);
+
+			/* Skip twin-proof reflections */
+			if ( (ha==hra) && (ka==kra) && (la==lra) ) {
+				//STATUS("%i %i %i is twin proof\n", h, k, l);
+				continue;
+			}
+
+		}
+
+		cr = find_refl(asym, ha, ka, la);
+		if ( cr == NULL ) {
+			cr = add_refl(asym, ha, ka, la);
+			assert(cr != NULL);
+			copy_data(cr, refl);
+		} else {
+			const double i = get_intensity(cr);
+			const int r = get_redundancy(cr);
+			set_intensity(cr, (r*i + get_intensity(refl))/(r+1));
+			set_redundancy(cr, r+1);
+		}
+	}
+
+	f = malloc(sizeof(struct flist));
+	if ( f == NULL ) {
+		ERROR("Failed to allocate flist\n");
+		return NULL;
+	}
+
+	n = num_reflections(asym);
+	f->s = malloc(n*sizeof(unsigned int));
+	f->s_reidx = malloc(n*sizeof(unsigned int));
+	f->i = malloc(n*sizeof(float));
+	f->i_reidx = malloc(n*sizeof(float));
+	if ( (f->s == NULL) || (f->i == NULL)
+	   || (f->s_reidx == NULL) || (f->i_reidx == NULL) ) {
+		ERROR("Failed to allocate flist\n");
+		return NULL;
+	}
+
+	f->n = 0;
+	for ( refl = first_refl(asym, &iter);
+	      refl != NULL;
+	      refl = next_refl(refl, iter) )
+	{
+		signed int h, k, l;
+
+		get_indices(refl, &h, &k, &l);
+		f->s[f->n] = SERIAL(h, k, l);
+
+		f->i[f->n] = get_intensity(refl);
+		f->n++;
+	}
+	assert(f->n == n);
+
+	if ( amb != NULL ) {
+
+		RefList *reidx = reflist_new();
+		if ( reidx == NULL ) return NULL;
+
+		for ( refl = first_refl(asym, &iter);
+		      refl != NULL;
+		      refl = next_refl(refl, iter) )
+		{
+			signed int h, k, l;
+			signed int hr, kr, lr;
+			signed int hra, kra, lra;
+			Reflection *cr;
+
+			get_indices(refl, &h, &k, &l);
+			get_equiv(amb, NULL, 1, h, k, l, &hr, &kr, &lr);
+			get_asymm(sym, hr, kr, lr, &hra, &kra, &lra);
+
+			cr = add_refl(reidx, hra, kra, lra);
+			copy_data(cr, refl);
+		}
+
+		n = 0;
+		for ( refl = first_refl(reidx, &iter);
+		      refl != NULL;
+		      refl = next_refl(refl, iter) )
+		{
+			signed int h, k, l;
+			get_indices(refl, &h, &k, &l);
+			f->s_reidx[n] = SERIAL(h, k, l);
+			f->i_reidx[n++] = get_intensity(refl);
+		}
+
+		reflist_free(reidx);
+	}
+
+	reflist_free(asym);
+
+	return f;
+}
+
+
+static float corr(struct flist *a, struct flist *b, int *pn, int a_reidx)
+{
+	float s_xy = 0.0;
+	float s_x = 0.0;
+	float s_y = 0.0;
+	float s_x2 = 0.0;
+	float s_y2 = 0.0;
+	int n = 0;
+	float t1, t2;
+	int ap = 0;
+	int bp = 0;
+	int done = 0;
+	unsigned int *sa;
+	float *ia;
+
+	if ( a_reidx ) {
+		sa = a->s_reidx;
+		ia = a->i_reidx;
+	} else {
+		sa = a->s;
+		ia = a->i;
+	}
+
+	if ( (a->n == 0) || (b->n == 0) ) {
+		*pn = 0;
+		return 0.0;
+	}
+
+	while ( 1 ) {
+
+		while ( sa[ap] > b->s[bp] ) {
+			if ( ++bp == b->n ) {
+				done = 1;
+				break;
+			}
+		}
+		if ( done ) break;
+
+		while ( sa[ap] < b->s[bp] ) {
+			if ( ++ap == a->n ) {
+				done = 1;
+				break;
+			}
+		}
+		if ( done ) break;
+
+		if ( sa[ap] == b->s[bp] ) {
+
+			float aint, bint;
+
+			aint = ia[ap];
+			bint = b->i[bp];
+
+			s_xy += aint*bint;
+			s_x += aint;
+			s_y += bint;
+			s_x2 += aint*aint;
+			s_y2 += bint*bint;
+			n++;
+
+			if ( ++ap == a->n ) break;
+			if ( ++bp == b->n ) break;
+
+		}
+
+	}
+	*pn = n;
+
+	t1 = s_x2 - s_x*s_x / n;
+	t2 = s_y2 - s_y*s_y / n;
+
+	if ( (t1 <= 0.0) || (t2 <= 0.0) ) return 0.0;
+
+	return (s_xy - s_x*s_y/n) / sqrt(t1*t2);
+}
+
+
+struct cc_list
+{
+	signed int *ind;
+	float *cc;
+
+	signed int *ind_reidx;
+	float *cc_reidx;
+};
+
+
+struct queue_args
+{
+	int n_started;
+	int n_finished;
+	int n_to_do;
+	int mean_nac;
+	int nmean_nac;
+
+	struct cc_list *ccs;
+	struct flist **crystals;
+	int n_crystals;
+	int ncorr;
+	SymOpList *amb;
+	gsl_rng **rngs;
+};
+
+
+struct cc_job
+{
+	struct cc_list *ccs;
+	int i;
+	int mean_nac;
+	int nmean_nac;
+
+	struct flist **crystals;
+	int n_crystals;
+	int ncorr;
+	SymOpList *amb;
+	gsl_rng **rngs;
+};
+
+
+static void *get_task(void *vp)
+{
+	struct queue_args *qargs = vp;
+	struct cc_job *job;
+
+	if ( qargs->n_started == qargs->n_to_do ) return NULL;
+
+	job = malloc(sizeof(struct cc_job));
+	if ( job == NULL ) return NULL;
+
+	job->ccs = qargs->ccs;
+	job->i = qargs->n_started++;
+
+	job->crystals = qargs->crystals;
+	job->n_crystals = qargs->n_crystals;
+	job->ncorr = qargs->ncorr;
+	job->amb = qargs->amb;
+	job->rngs = qargs->rngs;
+
+	return job;
+}
+
+
+static void final(void *qp, void *wp)
+{
+	struct queue_args *qargs = qp;
+	struct cc_job *job = wp;
+
+	qargs->mean_nac += job->mean_nac;
+	qargs->nmean_nac += job->nmean_nac;
+
+	free(job);
+
+	qargs->n_finished++;
+	progress_bar(qargs->n_finished, qargs->n_to_do, "Calculating CCs");
+}
+
+
+static void work(void *wp, int cookie)
+{
+	struct cc_job *job = wp;
+	int i = job->i;
+	int k, l;
+	struct cc_list *ccs = job->ccs;
+	struct flist **crystals = job->crystals;
+	int n_crystals = job->n_crystals;
+	int ncorr = job->ncorr;
+	SymOpList *amb = job->amb;
+	int mean_nac = 0;
+	int nmean_nac = 0;
+	gsl_permutation *p;
+
+	p = gsl_permutation_alloc(n_crystals);
+	gsl_permutation_init(p);
+	gsl_ran_shuffle(job->rngs[cookie], p->data, n_crystals, sizeof(size_t));
+
+	ccs[i].ind = malloc(ncorr*sizeof(int));
+	ccs[i].cc = malloc(ncorr*sizeof(float));
+	ccs[i].ind_reidx = calloc(ncorr, sizeof(int));
+	ccs[i].cc_reidx = calloc(ncorr, sizeof(float));
+	if ( (ccs[i].ind==NULL) || (ccs[i].cc==NULL) ||
+	     (ccs[i].ind_reidx==NULL) ||  (ccs[i].cc_reidx==NULL) ) {
+		return;
+	}
+
+	k = 0;
+	for ( l=0; l<n_crystals; l++ ) {
+
+		int n;
+		int j;
+		float cc;
+
+		j = gsl_permutation_get(p, l);
+		if ( i == j ) continue;
+
+		cc = corr(crystals[i], crystals[j], &n, 0);
+
+		if ( n < 4 ) continue;
+
+		ccs[i].ind[k] = j+1;
+		ccs[i].cc[k] = cc;
+		k++;
+
+		if ( k == ncorr-1 ) break;
+
+	}
+	ccs[i].ind[k] = 0;
+	mean_nac += k;
+	nmean_nac++;
+
+	if ( amb != NULL ) {
+
+		k = 0;
+		for ( l=0; l<n_crystals; l++ ) {
+
+			int n;
+			int j;
+			float cc;
+
+			j = gsl_permutation_get(p, l);
+			if ( i == j ) continue;
+
+			cc = corr(crystals[i], crystals[j], &n, 1);
+
+			if ( n < 4 ) continue;
+
+			ccs[i].ind_reidx[k] = j+1;
+			ccs[i].cc_reidx[k] = cc;
+			k++;
+
+			if ( k == ncorr-1 ) break;
+
+		}
+		ccs[i].ind_reidx[k] = 0;
+		mean_nac += k;
+		nmean_nac++;
+
+	}
+
+	gsl_permutation_free(p);
+
+	job->mean_nac = mean_nac;
+	job->nmean_nac = nmean_nac;
+}
+
+
+static gsl_rng **setup_random(gsl_rng *rng, int n)
+{
+	gsl_rng **rngs;
+	int i;
+
+	rngs = malloc(n * sizeof(gsl_rng *));
+	if ( rngs == NULL ) return NULL;
+
+	for ( i=0; i<n; i++ ) {
+		rngs[i] = gsl_rng_alloc(gsl_rng_mt19937);
+		gsl_rng_set(rngs[i], gsl_rng_get(rng));
+	}
+
+	return rngs;
+}
+
+
+static struct cc_list *calc_ccs(struct flist **crystals, int n_crystals,
+                                int ncorr, SymOpList *amb, gsl_rng *rng,
+                                float *pmean_nac, int nthreads)
+{
+	struct cc_list *ccs;
+	struct queue_args qargs;
+	int i;
+
+	assert(n_crystals >= ncorr);
+	ncorr++;  /* Extra value at end for sentinel */
+
+	qargs.rngs = setup_random(rng, nthreads);
+
+	ccs = malloc(n_crystals*sizeof(struct cc_list));
+	if ( ccs == NULL ) return NULL;
+
+	qargs.n_started = 0;
+	qargs.n_finished = 0;
+	qargs.n_to_do = n_crystals;
+	qargs.ccs = ccs;
+	qargs.mean_nac = 0;
+	qargs.nmean_nac = 0;
+
+	qargs.crystals = crystals;
+	qargs.n_crystals = n_crystals;
+	qargs.ncorr = ncorr;
+	qargs.amb = amb;
+
+	run_threads(nthreads, work, get_task, final, &qargs, n_crystals,
+	            0, 0, 0);
+
+	for ( i=0; i<nthreads; i++ ) {
+		gsl_rng_free(qargs.rngs[i]);
+	}
+
+	*pmean_nac = (float)qargs.mean_nac/qargs.nmean_nac;
+
+	return ccs;
+}
+
+
+static void detwin(struct cc_list *ccs, int n_crystals, int *assignments,
+                   FILE *fh, struct flist **crystals)
+{
+	int i;
+	int nch = 0;
+	float mf = 0.0;
+	int nmf = 0;
+	int ndud = 0;
+
+	for ( i=0; i<n_crystals; i++ ) {
+
+		int k;
+		float f = 0.0;
+		float g = 0.0;;
+		int p = 0;
+		int q = 0;
+
+		for ( k=0; (ccs[i].ind[k] != 0); k++ ) {
+
+			int j = ccs[i].ind[k]-1;
+			float cc = ccs[i].cc[k];
+
+			if ( assignments[i] == assignments[j] ) {
+				f += cc;
+				p++;
+			} else {
+				g += cc;
+				q++;
+			}
+
+		}
+
+		for ( k=0; (ccs[i].ind_reidx[k] != 0); k++ ) {
+
+			int j = ccs[i].ind_reidx[k]-1;
+			float cc = ccs[i].cc_reidx[k];
+
+			if ( assignments[i] == assignments[j] ) {
+				g += cc;
+				q++;
+			} else {
+				f += cc;
+				p++;
+			}
+
+		}
+
+		if ( (p==0) || (q==0) ) {
+			ndud++;
+			continue;
+		}
+
+		f /= p;
+		g /= q;
+
+		if ( fh != NULL ) fprintf(fh, "%5.3f %5.3f\n", f, g);
+
+		mf += f;
+		nmf++;
+
+		if ( f < g ) {
+			assignments[i] = 1 - assignments[i];
+			nch++;
+		}
+
+	}
+
+	if ( ndud > 0 ) {
+		STATUS("Warning: %i crystals had no correlation\n", ndud);
+	}
+
+	STATUS("Mean f = %f, changed %i assignments this time.\n", mf/nmf, nch);
+}
+
+
+static void reindex_reflections(FILE *fh, FILE *ofh, int assignment,
+                                SymOpList *amb)
+{
+	int first = 1;
+
+	do {
+
+		char *rval;
+		char line[1024];
+		int n;
+		signed int h, k, l;
+		int r;
+
+		rval = fgets(line, 1023, fh);
+		if ( rval == NULL ) break;
+
+		if ( strcmp(line, REFLECTION_END_MARKER"\n") == 0 ) {
+			fputs(line, ofh);
+			return;
+		}
+
+		if ( first ) {
+			fputs(line, ofh);
+			first = 0;
+			continue;
+		}
+
+		r = sscanf(line, "%i %i %i%n", &h, &k, &l, &n);
+
+		/* See scanf() manual page about %n to see why <3 is used */
+		if ( (r < 3) && !first ) return;
+
+		get_equiv(amb, NULL, assignment, h, k, l, &h, &k, &l);
+
+		fprintf(ofh, "%4i %4i %4i%s", h, k, l, line+n);
+
+	} while ( 1 );
+}
+
+
+/* This is nasty, but means the output includes absolutely everything in the
+ * input, even stuff ignored by read_chunk() */
+static void write_reindexed_stream(const char *infile, const char *outfile,
+                                   int *assignments, SymOpList *amb)
+{
+	FILE *fh;
+	FILE *ofh;
+	int i;
+
+	fh = fopen(infile, "r");
+	if ( fh == NULL ) {
+		ERROR("Failed to open '%s'\n", infile);
+		return;
+	}
+
+	ofh = fopen(outfile, "w");
+	if ( ofh == NULL ) {
+		ERROR("Failed to open '%s'\n", outfile);
+		return;
+	}
+
+	i = 0;
+	do {
+
+		char *rval;
+		char line[1024];
+
+		rval = fgets(line, 1023, fh);
+		if ( rval == NULL ) break;
+
+		fputs(line, ofh);
+
+		if ( strcmp(line, REFLECTION_START_MARKER"\n") == 0 ) {
+			reindex_reflections(fh, ofh, assignments[i++], amb);
+		}
+
+	} while ( 1 );
+
+	if ( !feof(fh) ) {
+		ERROR("Error reading stream.\n");
+	}
+
+	fclose(fh);
+	fclose(ofh);
+}
+
+
+int main(int argc, char *argv[])
+{
+	int c;
+	const char *infile;
+	char *outfile = NULL;
+	char *end_ass_fn = NULL;
+	char *fg_graph_fn = NULL;
+	char *s_sym_str = NULL;
+	SymOpList *s_sym;
+	char *w_sym_str = NULL;
+	SymOpList *w_sym;
+	SymOpList *amb;
+	int n_iter = 6;
+	int n_crystals, n_chunks, max_crystals;
+	int n_dif;
+	struct flist **crystals;
+	Stream *st;
+	int i;
+	int *assignments;
+	int *orig_assignments;
+	gsl_rng *rng;
+	float highres, lowres;
+	double rmin = 0.0;  /* m^-1 */
+	double rmax = INFINITY;  /* m^-1 */
+	FILE *fgfh = NULL;
+	struct cc_list *ccs;
+	int ncorr;
+	int ncorr_set = 0;
+	float mean_nac;
+	int n_threads = 1;
+
+	/* Long options */
+	const struct option longopts[] = {
+		{"help",               0, NULL,               'h'},
+		{"output",             1, NULL,               'o'},
+		{"symmetry",           1, NULL,               'y'},
+		{"iterations",         1, NULL,               'n'},
+
+		{"highres",            1, NULL,                2},
+		{"lowres",             1, NULL,                3},
+		{"end-assignments",    1, NULL,                4},
+		{"fg-graph",           1, NULL,                5},
+		{"ncorr",              1, NULL,                6},
+
+		{0, 0, NULL, 0}
+	};
+
+	/* Short options */
+	while ((c = getopt_long(argc, argv, "ho:y:n:w:j:",
+	                        longopts, NULL)) != -1)
+	{
+
+		switch (c) {
+
+			case 'h' :
+			show_help(argv[0]);
+			return 0;
+
+			case 'o' :
+			outfile = strdup(optarg);
+			break;
+
+			case 'y' :
+			s_sym_str = strdup(optarg);
+			break;
+
+			case 'w' :
+			w_sym_str = strdup(optarg);
+			break;
+
+			case 'n' :
+			n_iter = atoi(optarg);
+			break;
+
+			case 'j' :
+			if ( sscanf(optarg, "%i", &n_threads) != 1 ) {
+				ERROR("Invalid value for -j\n");
+				return 1;
+			}
+			break;
+
+			case 2 :
+			if ( sscanf(optarg, "%e", &highres) != 1 ) {
+				ERROR("Invalid value for --highres\n");
+				return 1;
+			}
+			rmax = 1.0 / (highres/1e10);
+			break;
+
+			case 3 :
+			if ( sscanf(optarg, "%e", &lowres) != 1 ) {
+				ERROR("Invalid value for --lowres\n");
+				return 1;
+			}
+			rmin = 1.0 / (lowres/1e10);
+			break;
+
+			case 4 :
+			end_ass_fn = strdup(optarg);
+			break;
+
+			case 5 :
+			fg_graph_fn = strdup(optarg);
+			break;
+
+			case 6 :
+			if ( sscanf(optarg, "%i", &ncorr) != 1 ) {
+				ERROR("Invalid value for --ncorr\n");
+				return 1;
+			} else {
+				ncorr_set = 1;
+			}
+			break;
+
+			case 0 :
+			break;
+
+			case '?' :
+			break;
+
+			default :
+			ERROR("Unhandled option '%c'\n", c);
+			break;
+
+		}
+
+	}
+
+	if ( argc != (optind+1) ) {
+		ERROR("Please provide exactly one stream filename.\n");
+		return 1;
+	}
+
+	infile = argv[optind++];
+	st = open_stream_for_read(infile);
+	if ( st == NULL ) {
+		ERROR("Failed to open input stream '%s'\n", infile);
+		return 1;
+	}
+
+	/* Sanitise output filename */
+	if ( outfile == NULL ) {
+		outfile = strdup("partialator.hkl");
+	}
+
+	if ( s_sym_str == NULL ) {
+		ERROR("You must specify the input symmetry (with -y)\n");
+		return 1;
+	}
+	s_sym = get_pointgroup(s_sym_str);
+	if ( s_sym == NULL ) return 1;
+	free(s_sym_str);
+
+	if ( w_sym_str == NULL ) {
+		w_sym = NULL;
+		amb = NULL;
+	} else {
+		w_sym = get_pointgroup(w_sym_str);
+		free(w_sym_str);
+		if ( w_sym == NULL ) return 1;
+		amb = get_ambiguities(w_sym, s_sym);
+		if ( amb == NULL ) {
+			ERROR("Couldn't find ambiguity operator.\n");
+			ERROR("Check that your values for -y and -w are "
+			      "correct.\n");
+			return 1;
+		}
+		STATUS("Ambiguity operations:\n");
+		describe_symmetry(amb);
+		if ( num_equivs(amb, NULL) != 2 ) {
+			ERROR("There must be only one ambiguity operator.\n");
+			ERROR("Try again with a different value for -w.\n");
+			return 1;
+		}
+	}
+
+	crystals = NULL;
+	n_crystals = 0;
+	max_crystals = 0;
+	n_chunks = 0;
+	do {
+
+		struct image cur;
+		int i;
+
+		cur.det = NULL;
+
+		if ( read_chunk(st, &cur) != 0 ) {
+			break;
+		}
+
+		image_feature_list_free(cur.features);
+
+		for ( i=0; i<cur.n_crystals; i++ ) {
+
+			Crystal *cr;
+			RefList *list;
+			UnitCell *cell;
+
+			cr = cur.crystals[i];
+			cell = crystal_get_cell(cr);
+
+			if ( n_crystals == max_crystals ) {
+
+				struct flist **crystals_new;
+				size_t nsz;
+
+				nsz = (max_crystals+1024)*sizeof(struct flist *);
+				crystals_new = realloc(crystals, nsz);
+				if ( crystals_new == NULL ) {
+					fprintf(stderr, "Failed to allocate "
+					        "memory for crystals.\n");
+					break;
+				}
+
+				max_crystals += 1024;
+				crystals = crystals_new;
+
+			}
+
+			list = crystal_get_reflections(cr);
+			crystals[n_crystals] = asymm_and_merge(list, s_sym,
+			                                       cell,
+			                                       rmin, rmax,
+			                                       amb);
+			cell_free(cell);
+			n_crystals++;
+			reflist_free(list);
+
+		}
+
+		fprintf(stderr, "Loaded %i crystals from %i chunks\r",
+		        n_crystals, ++n_chunks);
+
+	} while ( 1 );
+	fprintf(stderr, "\n");
+
+	close_stream(st);
+
+	assignments = malloc(n_crystals*sizeof(int));
+	if ( assignments == NULL ) {
+		ERROR("Couldn't allocate memory for assignments.\n");
+		return 1;
+	}
+
+	orig_assignments = malloc(n_crystals*sizeof(int));
+	if ( orig_assignments == NULL ) {
+		ERROR("Couldn't allocate memory for original assignments.\n");
+		return 1;
+	}
+
+	rng = gsl_rng_alloc(gsl_rng_mt19937);
+	for ( i=0; i<n_crystals; i++ ) {
+		assignments[i] = (random_flat(rng, 1.0) > 0.5);
+		orig_assignments[i] = assignments[i];
+	}
+
+	if ( fg_graph_fn != NULL ) {
+		fgfh = fopen(fg_graph_fn, "w");
+		if ( fgfh == NULL ) {
+			ERROR("Failed to open '%s'\n", fg_graph_fn);
+		}
+	}
+
+	if ( !ncorr_set || (ncorr > n_crystals) ) {
+		ncorr = n_crystals;
+	}
+
+	ccs = calc_ccs(crystals, n_crystals, ncorr, amb, rng, &mean_nac,
+	               n_threads);
+	if ( ccs == NULL ) {
+		ERROR("Failed to allocate CCs\n");
+		return 1;
+	}
+	STATUS("Mean number of correlations per crystal: %.1f\n", mean_nac);
+
+	for ( i=0; i<n_crystals; i++ ) {
+		free(crystals[i]->s);
+		free(crystals[i]->i);
+		free(crystals[i]->s_reidx);
+		free(crystals[i]->i_reidx);
+		free(crystals[i]);
+	}
+	free(crystals);
+
+	for ( i=0; i<n_iter; i++ ) {
+		detwin(ccs, n_crystals, assignments, fgfh, crystals);
+	}
+
+	if ( fgfh != NULL ) {
+		fclose(fgfh);
+	}
+
+	if ( end_ass_fn != NULL ) {
+		FILE *fh = fopen(end_ass_fn, "w");
+		if ( fh == NULL ) {
+			ERROR("Failed to open '%s'\n", end_ass_fn);
+		} else {
+			for ( i=0; i<n_crystals; i++ ) {
+				fprintf(fh, "%i\n", assignments[i]);
+			}
+		}
+		fclose(fh);
+	}
+
+	n_dif = 0;
+	for ( i=0; i<n_crystals; i++ ) {
+		if ( orig_assignments[i] != assignments[i] ) n_dif++;
+	}
+	STATUS("%i assignments are different from their starting values.\n",
+	       n_dif);
+
+	if ( amb != NULL ) {
+		write_reindexed_stream(infile, outfile, assignments, amb);
+	} else {
+		ERROR("Can only write stream with known ambiguity operator.\n");
+		ERROR("Try again with -w\n");
+	}
+
+	free(assignments);
+	gsl_rng_free(rng);
+
+	return 0;
+}
diff --git a/tests/first_merge_check b/tests/first_merge_check
index d649e284..5cf7cc70 100755
--- a/tests/first_merge_check
+++ b/tests/first_merge_check
@@ -35,8 +35,10 @@ End of reflections
 EOF
 
 cat > first_merge_check_ans.hkl << EOF
-  h   k   l          I    phase   sigma(I)  counts  fs/px  ss/px
-  1   0   0     150.03        -      35.36       2    0.0    0.0
+CrystFEL reflection list version 2.0
+Symmetry: 1
+   h    k    l          I    phase   sigma(I)   nmeas
+   1    0    0     150.03        -      35.36       2
 End of reflections
 EOF
 
diff --git a/tests/fourth_merge_check b/tests/fourth_merge_check
index e2d2487f..bbe32b9b 100755
--- a/tests/fourth_merge_check
+++ b/tests/fourth_merge_check
@@ -20,7 +20,9 @@ End of reflections
 EOF
 
 cat > fourth_merge_check_ans.hkl << EOF
-  h   k   l          I    phase   sigma(I)  counts  fs/px  ss/px
+CrystFEL reflection list version 2.0
+Symmetry: 1
+   h    k    l          I    phase   sigma(I)   nmeas
 End of reflections
 EOF
 
diff --git a/tests/partialator_merge_check_1 b/tests/partialator_merge_check_1
index 05b74e97..f1a6adf5 100755
--- a/tests/partialator_merge_check_1
+++ b/tests/partialator_merge_check_1
@@ -32,8 +32,10 @@ EOF
 # We merge two patterns, without scaling or partiality, the result should just
 # be an average.
 cat > partialator_merge_check_1_ans.hkl << EOF
-  h   k   l          I    phase   sigma(I)  counts  fs/px  ss/px
-  1   0   0     150.00        -      35.36       2    0.0    0.0
+CrystFEL reflection list version 2.0
+Symmetry: unknown
+   h    k    l          I    phase   sigma(I)   nmeas
+   1    0    0     150.00        -      35.36       2
 End of reflections
 EOF
 
diff --git a/tests/partialator_merge_check_2 b/tests/partialator_merge_check_2
index be48249c..5513f5ec 100755
--- a/tests/partialator_merge_check_2
+++ b/tests/partialator_merge_check_2
@@ -33,8 +33,10 @@ EOF
 # be the mean but with the standard deviation should be zero because the scaling
 # factor can absorb the difference.
 cat > partialator_merge_check_2_ans.hkl << EOF
-  h   k   l          I    phase   sigma(I)  counts  fs/px  ss/px
-  1   0   0     150.00        -       0.00       2    0.0    0.0
+CrystFEL reflection list version 2.0
+Symmetry: unknown
+   h    k    l          I    phase   sigma(I)   nmeas
+   1    0    0     150.00        -       0.00       2
 End of reflections
 EOF
 
diff --git a/tests/partialator_merge_check_3 b/tests/partialator_merge_check_3
index 9a2293fd..274d7c2d 100755
--- a/tests/partialator_merge_check_3
+++ b/tests/partialator_merge_check_3
@@ -34,9 +34,11 @@ EOF
 
 # W
 cat > partialator_merge_check_3_ans.hkl << EOF
-  h   k   l          I    phase   sigma(I)  counts  fs/px  ss/px
-  1   0   0      78.40        -       3.15       2    0.0    0.0
-  2   0   0     145.51        -       1.70       2    0.0    0.0
+CrystFEL reflection list version 2.0
+Symmetry: unknown
+   h    k    l          I    phase   sigma(I)   nmeas
+   1    0    0      78.40        -       3.15       2
+   2    0    0     145.51        -       1.70       2
 End of reflections
 EOF
 
diff --git a/tests/second_merge_check b/tests/second_merge_check
index aec5b120..d34cc75e 100755
--- a/tests/second_merge_check
+++ b/tests/second_merge_check
@@ -35,8 +35,10 @@ End of reflections
 EOF
 
 cat > second_merge_check_ans.hkl << EOF
-  h   k   l          I    phase   sigma(I)  counts  fs/px  ss/px
-  1   0   0     150.03        -      35.36       2    0.0    0.0
+CrystFEL reflection list version 2.0
+Symmetry: -1
+   h    k    l          I    phase   sigma(I)   nmeas
+   1    0    0     150.03        -      35.36       2
 End of reflections
 EOF
 
diff --git a/tests/third_merge_check b/tests/third_merge_check
index 569d5ddf..81f9ad78 100755
--- a/tests/third_merge_check
+++ b/tests/third_merge_check
@@ -50,8 +50,10 @@ End of reflections
 EOF
 
 cat > third_merge_check_ans.hkl << EOF
-  h   k   l          I    phase   sigma(I)  counts  fs/px  ss/px
-  1   0   0     133.36        -      27.22       3    0.0    0.0
+CrystFEL reflection list version 2.0
+Symmetry: 1
+   h    k    l          I    phase   sigma(I)   nmeas
+   1    0    0     133.36        -      27.22       3
 End of reflections
 EOF