compare_hkl: Reject reflections with low SNR

src/compare_hkl.c

@@ -142,6 +142,10 @@ int main(int argc, char *argv[])
 	ReflItemList *i1, *i2, *icommon;
 	int config_luzzati = 0;
 	char *pdb = NULL;
+	double *esd1;
+	double *esd2;
+	int rej1 = 0;
+	int rej2 = 0;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -202,13 +206,15 @@ int main(int argc, char *argv[])
 	free(pdb);
 
 	ref1 = new_list_intensity();
-	i1 = read_reflections(afile, ref1, NULL, NULL);
+	esd1 = new_list_sigma();
+	i1 = read_reflections(afile, ref1, NULL, NULL, esd1);
 	if ( i1 == NULL ) {
 		ERROR("Couldn't open file '%s'\n", afile);
 		return 1;
 	}
 	ref2 = new_list_intensity();
-	i2 = read_reflections(bfile, ref2, NULL, NULL);
+	esd2 = new_list_sigma();
+	i2 = read_reflections(bfile, ref2, NULL, NULL, esd2);
 	if ( i2 == NULL ) {
 		ERROR("Couldn't open file '%s'\n", bfile);
 		return 1;
@@ -226,6 +232,8 @@ int main(int argc, char *argv[])
 		signed int h, k, l;
 		signed int he, ke, le;
 		double val1, val2;
+		double sig1, sig2;
+		int ig = 0;
 
 		it = get_item(i1, i);
 		h = it->h;  k = it->k;  l = it->l;
@@ -238,12 +246,29 @@ int main(int argc, char *argv[])
 
 		val1 = lookup_intensity(ref1, h, k, l);
 		val2 = lookup_intensity(ref2, he, ke, le);
+		sig1 = lookup_sigma(esd1, h, k, l);
+		sig2 = lookup_sigma(esd2, h, k, l);
+
+		if ( val1 < 3.0 * sig1 ) {
+			rej1++;
+			ig = 1;
+		}
+		if ( val2 < 3.0 * sig2 ) {
+			rej2++;
+			ig = 1;
+		}
+		if ( ig ) continue;
+
 		set_intensity(ref2_transformed, h, k, l, val2);
 		set_intensity(out, h, k, l, val1/val2);
 		add_item(icommon, h, k, l);
 
 	}
 	ncom = num_items(icommon);
+	STATUS("%i reflections with I < 3.0*sigma(I) rejected from '%s'\n",
+	       rej1, afile);
+	STATUS("%i reflections with I < 3.0*sigma(I) rejected from '%s'\n",
+	       rej2, bfile);
 
 	STATUS("%i,%i reflections: %i in common\n",
 	       num_items(i1), num_items(i2), ncom);

src/get_hkl.c

@@ -238,7 +238,8 @@ int main(int argc, char *argv[])
 		                            phases, input_items);
 	} else {
 		ideal_ref = new_list_intensity();
-		input_items = read_reflections(input, ideal_ref, phases, NULL);
+		input_items = read_reflections(input, ideal_ref, phases,
+		                               NULL, NULL);
 		free(input);
 	}
 
@@ -275,7 +276,8 @@ int main(int argc, char *argv[])
 		/* Write out only reflections which are in the template
 		 * (and which we have in the input) */
 		ReflItemList *template_items;
-		template_items = read_reflections(template, NULL, NULL, NULL);
+		template_items = read_reflections(template,
+		                                  NULL, NULL, NULL, NULL);
 		write_items = intersection_items(input_items, template_items);
 		delete_items(template_items);
 	} else {

src/indexamajig.c

@@ -597,7 +597,8 @@ int main(int argc, char *argv[])
 
 	if ( intfile != NULL ) {
 		ReflItemList *items;
-		items = read_reflections(intfile, intensities, NULL, NULL);
+		items = read_reflections(intfile, intensities,
+		                         NULL, NULL, NULL);
 		delete_items(items);
 	} else {
 		intensities = NULL;

src/pattern_sim.c

@@ -332,7 +332,8 @@ int main(int argc, char *argv[])
 		}
 		intensities = new_list_intensity();
 		phases = new_list_phase();
-		items = read_reflections(intfile, intensities, phases, NULL);
+		items = read_reflections(intfile, intensities, phases,
+		                         NULL, NULL);
 		free(intfile);
 		delete_items(items);
 	}

src/process_hkl.c

@@ -696,7 +696,7 @@ int main(int argc, char *argv[])
 	if ( reference != NULL ) {
 		reference_i = new_list_intensity();
 		reference_items = read_reflections(reference, reference_i,
-		                                   NULL, NULL);
+		                                   NULL, NULL, NULL);
 		if ( reference_items == NULL ) {
 			ERROR("Couldn't read '%s'\n", reference);
 			return 1;

src/reflections.c

@@ -102,7 +102,7 @@ void write_reflections(const char *filename, ReflItemList *items,
  */
 ReflItemList *read_reflections(const char *filename,
                                double *intensities, double *phases,
-                               unsigned int *counts)
+                               unsigned int *counts, double *esds)
 {
 	FILE *fh;
 	char *rval;
@@ -163,6 +163,9 @@ ReflItemList *read_reflections(const char *filename,
 		if ( counts != NULL ) {
 			set_count(counts, h, k, l, cts);
 		}
+		if ( esds != NULL ) {
+			set_sigma(esds, h, k, l, sigma);
+		}
 
 	} while ( rval != NULL );
 

src/reflections.h

@@ -27,7 +27,7 @@ extern void write_reflections(const char *filename, ReflItemList *items,
 
 extern ReflItemList *read_reflections(const char *filename,
                                       double *intensities, double *phases,
-                                      unsigned int *counts);
+                                      unsigned int *counts, double *esds);
 
 extern double *ideal_intensities(double complex *sfac);
 

src/render_hkl.c

@@ -609,7 +609,7 @@ int main(int argc, char *argv[])
 	}
 	ref = new_list_intensity();
 	cts = new_list_count();
-	ReflItemList *items = read_reflections(infile, ref, NULL, cts);
+	ReflItemList *items = read_reflections(infile, ref, NULL, cts, NULL);
 	if ( ref == NULL ) {
 		ERROR("Couldn't open file '%s'\n", infile);
 		return 1;

src/utils.h

@@ -175,6 +175,11 @@ static inline double angle_between(double x1, double y1, double z1,
 #define TYPE unsigned int
 #include "list_tmp.h"
 
+/* As above, but for sigmas */
+#define LABEL(x) x##_sigma
+#define TYPE double
+#include "list_tmp.h"
+
 
 /* ----------- Reflection lists indexed by sequence (not indices) ----------- */