compare_hkl: Add --shell-file, update documentation

2 files changed, 61 insertions, 29 deletions

diff --git a/doc/man/compare_hkl.1 b/doc/man/compare_hkl.1
index 12f19489..71c3d544 100644
--- a/doc/man/compare_hkl.1
+++ b/doc/man/compare_hkl.1
@@ -16,7 +16,7 @@ compare_hkl \- compare reflection data
 \fBcompare_hkl --help\fR
 
 .SH DESCRIPTION
-compare_hkl compares two sets of reflection data and calculates figures of merit such as R-factors.  Several figures of merit will be displayed on the command line.  Figures of merit can also be calculated in resolution shells.  Reflections will be considered equivalent according to your choice of point group.  The ratios of the reflections can be written to a new file, perhaps to look for systematic effects.
+compare_hkl compares two sets of reflection data and calculates figures of merit such as R-factors.  Reflections will be considered equivalent according to your choice of point group.
 
 .SH OPTIONS
 .PD 0
@@ -32,49 +32,71 @@ Specify the symmetry of the reflections.  The symmetry must be the same for both
 Specify the name of the PDB file containing at least a CRYST1 line describing the unit cell.
 
 .PD 0
-.IP \fB--rmin=\fR\fI1/d\fR
+.IP \fB--fom=\fR\fIFoM\fR
+.PD
+Calculate figure of merit \fIFoM\fR.  Possible figures of merit are:
+.RS
+.IP \fBRsplit\fR
+.PD
+2^(-0.5) * sum(I1-kI2) / [ 0.5*sum(I1+kI2) ]
+.IP \fBR1f\fR
 .PD
-Fix the lower resolution limit for the resolutions shells, as 1/d in m^-1.  Reflections outside the specified resolution range will still be included in the calculation of overall figures of merit.
+sum(sqrt(I1)-sqrt(kI2)) / sum(sqrt(I1))
+.IP \fBR1i\fR
+.PD
+sum(I1-kI2) / sum(I1)
+.IP \fBR2\fR
+.PD
+sqrt(sum[(I1-kI2)^2] / sum(I1^2))
+.IP \fBCC\fR
+.PD
+The Pearson correlation coefficient.
+.IP \fBCC*\fR
+.PD
+The value "CC*", see Karplus and Diederichs, Science 336 (2012) p1030.
+.PP
+I1 and I2 are the intensities of the same reflection in both reflection lists.  The scale factor, k, is given by sum(i1*i2) / sum(i2^2), unless you use \fB-u\fR.
+.RE
 
 .PD 0
-.IP \fB--rmax=\fR\fI1/d\fR
+.IP \fB-u\fR
 .PD
-Fix the upper resolution limit for the resolutions shells, as 1/d in m^-1  Reflections outside the specified resolution range will still be included in the calculation of overall figures of merit.
+Fix scale factors to unity.
 
 .PD 0
-.IP \fB--sigma-cutoff=\fR\fIn\fR
+.IP \fB--nshells=\fR\fIn\fR
 .PD
-Discard reflections with I/sigma(I) < \fIn\fR.  Default: -infinity (no cutoff).
+Use \fIn\fR resolution shells.  Default: 10.
 
 .PD 0
-.IP \fB--shells=\fR\fIFoM\fR
+.IP \fB--shell-file=\fIfilename\fR
 .PD
-Calculate \fIFoM\fR in resolution shells.  Possible figures of merit for \fIFoM\fR are:
-.RS
-.IP \fBRsplit\fR
+Write the figure of merit in resolution shells to \fIfilename\fR.  Default: "shells.dat".
+
+.PD 0
+.IP \fB--ignore-negs\fR
 .PD
-2^(-0.5) * sum(I1-kI2) / [ 0.5*sum(I1+kI2) ]
-.IP \fBR1f\fR
+Ignore reflections with negative intensities.
+
+.PD 0
+.IP \fB--zero-negs\fR
 .PD
-sum(F1-kF2) / sum(f1)
-.IP \fBR1i\fR
+Set to zero the intensities of reflections with negative intensities.
+
+.PD 0
+.IP \fB--sigma-cutoff=\fR\fIn\fR
 .PD
-sum(I1-kI2) / sum(I1)
-.PP
-I1 and I2 are the intensities of the same reflection in both reflection lists.  F1 and F2 are the square roots of the intensities.  Negative intensity values are discarded if R1f is used.  The scale factor, k, is chosen to minimise the R-factor in all cases.
-.RE
+Discard reflections with I/sigma(I) < \fIn\fR.  Default: -infinity (no cutoff).
 
 .PD 0
-.IP \fB-p\fR \fIunitcell.pdb\fR
-.IP \fB--pdb=\fR\fIunitcell.pdb\fR
+.IP \fB--rmin=\fR\fI1/d\fR
 .PD
-Specify the name of the PDB file containing at least a CRYST1 line describing the unit cell.  This is required if \fB--shells\fR is used.
+Discard reflections with lower than this value of 1/d in m^-1.
 
 .PD 0
-.IP \fB-o\fR \fIfilename\fR
-.IP \fB--ratio=\fR\fIfilename\fR
+.IP \fB--rmax=\fR\fI1/d\fR
 .PD
-Write the ratios of the reflection intensities to \fIfilename\fR.  The ratio file will also be a CrystFEL reflection list (suitable for use with, for example, \fBrender_hkl\fR), but the intensities column in the file will actually contain the ratios of intensities.
+Discard reflections with higher than this value of 1/d in m^-1.
 
 .SH AUTHOR
 This page was written by Thomas White.
diff --git a/src/compare_hkl.c b/src/compare_hkl.c
index e25b4746..c6f46be4 100644
--- a/src/compare_hkl.c
+++ b/src/compare_hkl.c
@@ -86,6 +86,7 @@ static void show_help(const char *s)
 "                              CCF, CCI, CCFstar, CCIstar.\n"
 "      --nshells=<n>          Use <n> resolution shells.\n"
 "  -u                         Force scale factor to 1.\n"
+"      --shell-file=<file>    Write resolution shells to <file>.\n"
 "\n"
 "You can control which reflections are included in the calculation:\n"
 "\n"
@@ -316,7 +317,7 @@ static double fom_shell(struct fom_context *fctx, int i)
 
 static void do_fom(RefList *list1, RefList *list2, UnitCell *cell,
                    double rmin, double rmax, enum fom fom,
-                   int config_unity, int nshells)
+                   int config_unity, int nshells, const char *filename)
 {
 	int *cts;
 	double *rmins;
@@ -436,9 +437,9 @@ static void do_fom(RefList *list1, RefList *list2, UnitCell *cell,
 
 	}
 
-	fh = fopen("shells.dat", "w");
+	fh = fopen(filename, "w");
 	if ( fh == NULL ) {
-		ERROR("Couldn't open 'shells.dat'\n");
+		ERROR("Couldn't open '%s'\n", filename);
 		return;
 	}
 
@@ -530,6 +531,7 @@ int main(int argc, char *argv[])
 	int config_zeronegs = 0;
 	int config_unity = 0;
 	int nshells = 10;
+	char *shell_file = NULL;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -541,6 +543,7 @@ int main(int argc, char *argv[])
 		{"fom",                1, NULL,                4},
 		{"sigma-cutoff",       1, NULL,                5},
 		{"nshells",            1, NULL,                6},
+		{"shell-file",         1, NULL,                7},
 		{"ignore-negs",        0, &config_ignorenegs,  1},
 		{"zero-negs",          0, &config_zeronegs,    1},
 		{0, 0, NULL, 0}
@@ -604,6 +607,10 @@ int main(int argc, char *argv[])
 			}
 			break;
 
+			case 7 :
+			shell_file = strdup(optarg);
+			break;
+
 			default :
 			ERROR("Unhandled option '%c'\n", c);
 			break;
@@ -652,6 +659,8 @@ int main(int argc, char *argv[])
 		exit(1);
 	}
 
+	if ( shell_file == NULL ) shell_file = strdup("shells.dat");
+
 	cell = load_cell_from_pdb(pdb);
 	free(pdb);
 
@@ -812,8 +821,9 @@ int main(int argc, char *argv[])
 	reflist_free(list2);
 
 	do_fom(list1_acc, list2_acc, cell, rmin, rmax, fom, config_unity,
-	       nshells);
+	       nshells, shell_file);
 
+	free(shell_file);
 	reflist_free(list1_acc);
 	reflist_free(list2_acc);