process_hkl: Implement Rick's ESD calculation

6 files changed, 87 insertions, 73 deletions

diff --git a/src/compare_hkl.c b/src/compare_hkl.c
index 489e2b43..29dfe2a5 100644
--- a/src/compare_hkl.c
+++ b/src/compare_hkl.c
@@ -456,7 +456,8 @@ int main(int argc, char *argv[])
 
 	if ( outfile != NULL ) {
 
-		write_reflections(outfile, icommon, out, NULL, NULL, cell, 1.0);
+		write_reflections(outfile, icommon, out, NULL,
+		                  NULL, NULL, cell);
 		STATUS("Sigma(I) values in output file are not meaningful.\n");
 
 	}
diff --git a/src/facetron.c b/src/facetron.c
index 53fb00f2..f4148ae1 100644
--- a/src/facetron.c
+++ b/src/facetron.c
@@ -449,8 +449,7 @@ int main(int argc, char *argv[])
 	}
 
 	/* Output results */
-	write_reflections(outfile, obs, i_full, NULL, cts, NULL, 1.0);
-	STATUS("Sigma(I) values in output file are not (yet) meaningful.\n");
+	write_reflections(outfile, obs, i_full, NULL, NULL, cts, NULL);
 
 	/* Clean up */
 	free(i_full);
diff --git a/src/get_hkl.c b/src/get_hkl.c
index 8a8a71c2..7fb14f7e 100644
--- a/src/get_hkl.c
+++ b/src/get_hkl.c
@@ -429,8 +429,8 @@ int main(int argc, char *argv[])
 		adu_per_photon = beam->adu_per_photon;
 	}
 
-	write_reflections(output, write_items, ideal_ref, phases, NULL, cell,
-	                  adu_per_photon);
+	write_reflections(output, write_items, ideal_ref, NULL, phases,
+	                  NULL, cell);
 
 	delete_items(input_items);
 	delete_items(write_items);
diff --git a/src/process_hkl.c b/src/process_hkl.c
index c05adbe9..2533102a 100644
--- a/src/process_hkl.c
+++ b/src/process_hkl.c
@@ -220,7 +220,7 @@ static void merge_pattern(double *model, ReflItemList *observed,
                           double *hist_vals,
                           signed int hist_h, signed int hist_k,
                           signed int hist_l, int *hist_n, double *devs,
-                          double *tots, double *means, FILE *outfh)
+                          double *means, FILE *outfh)
 {
 	int i;
 	int twin;
@@ -275,11 +275,11 @@ static void merge_pattern(double *model, ReflItemList *observed,
 			}
 		}
 
-		if ( tots != NULL ) {
+		if ( devs != NULL ) {
 			double m;
 			m = lookup_intensity(means, h, k, l);
-			integrate_intensity(tots, h, k, l, intensity);
-			integrate_intensity(devs, h, k, l, fabs(intensity-m));
+			integrate_intensity(devs, h, k, l,
+			                    pow(intensity-m, 2.0));
 		}
 
 		if ( !find_item(sym_items, h, k, l) ) {
@@ -399,8 +399,7 @@ static void merge_all(FILE *fh, double **pmodel, ReflItemList **pobserved,
                       ReflItemList *reference, double *reference_i,
                       double *hist_vals,
                       signed int hist_h, signed int hist_k, signed int hist_l,
-                      int *hist_i, double *devs, double *tots, double *means,
-                      FILE *outfh)
+                      int *hist_i, double *devs, double *means, FILE *outfh)
 {
 	char *rval;
 	float f0;
@@ -470,7 +469,7 @@ static void merge_all(FILE *fh, double **pmodel, ReflItemList **pobserved,
 			              twins, holo, mero,
 			              reference, reference_i,
 			              hist_vals, hist_h, hist_k, hist_l,
-			              hist_i, devs, tots, means, outfh);
+			              hist_i, devs, means, outfh);
 
 			n_patterns++;
 			if ( n_patterns == config_stopafter ) break;
@@ -570,6 +569,8 @@ int main(int argc, char *argv[])
 	FILE *outfh = NULL;
 	struct beam_params *beam = NULL;
 	int space_for_hist = 0;
+	double *devs;
+	double *esds;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -771,7 +772,7 @@ int main(int argc, char *argv[])
 	          config_startafter, config_stopafter,
                   twins, in_sym, sym, n_total_patterns,
                   reference_items, reference_i,
-                  hist_vals, hist_h, hist_k, hist_l, &hist_i, NULL, NULL, NULL,
+                  hist_vals, hist_h, hist_k, hist_l, &hist_i, NULL, NULL,
                   outfh);
 	rewind(fh);
 	if ( space_for_hist && (hist_i >= space_for_hist) ) {
@@ -784,67 +785,51 @@ int main(int argc, char *argv[])
 		plot_histogram(hist_vals, hist_i);
 	}
 
-	if ( output != NULL ) {
-
-		double adu_per_photon;
-
-		if ( beam == NULL ) {
-			adu_per_photon = 167.0;
-			STATUS("No beam parameters file provided (use -b), "
-			       "so I'm assuming 167.0 ADU per photon.\n");
-		} else {
-			adu_per_photon = beam->adu_per_photon;
-		}
-
-		write_reflections(output, observed, model, NULL, counts, cell,
-		                  adu_per_photon);
-	}
-
-	if ( config_rmerge ) {
+	STATUS("Extra pass to calculate ESDs...\n");
+	devs = new_list_intensity();
+	esds = new_list_intensity();
+	rewind(fh);
+	merge_all(fh, &model, &observed, &counts,
+	          config_maxonly, config_scale, 0,
+	          config_startafter, config_stopafter,
+	          twins, in_sym, sym,
+	          n_total_patterns, reference_items, reference_i,
+	          NULL, 0, 0, 0, NULL, devs, model, NULL);
 
-		double *devs = new_list_intensity();
-		double *tots = new_list_intensity();
-		double total_dev = 0.0;
-		double total_tot = 0.0;
-		double total_pdev = 0.0;
-		double total_rdev = 0.0;
-		int i;
+	for ( i=0; i<num_items(observed); i++ ) {
 
-		STATUS("Extra pass to calculate figures of merit...\n");
+		struct refl_item *it;
+		signed int h, k, l;
+		double dev, esd;
+		unsigned int count;
 
-		rewind(fh);
-		merge_all(fh, &model, &observed, &counts,
-		          config_maxonly, config_scale, 0,
-		          config_startafter, config_stopafter,
-		          twins, in_sym, sym,
-		          n_total_patterns, reference_items, reference_i,
-		          NULL, 0, 0, 0, NULL, devs, tots, model, NULL);
+		it = get_item(observed, i);
+		h = it->h;  k = it->k,  l = it->l;
 
-		for ( i=0; i<num_items(observed); i++ ) {
+		dev = lookup_intensity(devs, h, k, l);
+		count = lookup_count(counts, h, k, l);
 
-			struct refl_item *it;
-			signed int h, k, l;
-			double dev;
-			unsigned int count;
+		if ( count < 2 ) continue;
 
-			it = get_item(observed, i);
-			h = it->h;  k = it->k,  l = it->l;
+		esd = sqrt(dev) / (double)count;
+		set_intensity(esds, h, k, l, esd);
 
-			dev = lookup_intensity(devs, h, k, l);
-			count = lookup_count(counts, h, k, l);
+	}
 
-			if ( count < 2 ) continue;
+	if ( output != NULL ) {
 
-			total_dev += dev;
-			total_pdev += sqrt(1.0/(count-1.0))*dev;
-			total_rdev += sqrt(count/(count-1.0))*dev;
-			total_tot += lookup_intensity(tots, h, k, l);
+		double adu_per_photon;
 
+		if ( beam == NULL ) {
+			adu_per_photon = 167.0;
+			STATUS("No beam parameters file provided (use -b), "
+			       "so I'm assuming 167.0 ADU per photon.\n");
+		} else {
+			adu_per_photon = beam->adu_per_photon;
 		}
 
-		STATUS("Rmerge = %f%%\n", 100.0*total_dev/total_tot);
-		STATUS("  Rpim = %f%%\n", 100.0*total_pdev/total_tot);
-		STATUS(" Rmeas = %f%%\n", 100.0*total_rdev/total_tot);
+		write_reflections(output, observed, model, esds,
+		                  NULL, counts, cell);
 
 	}
 
diff --git a/src/reflections.c b/src/reflections.c
index 0474d517..98fc9bab 100644
--- a/src/reflections.c
+++ b/src/reflections.c
@@ -22,10 +22,40 @@
 #include "beam-parameters.h"
 
 
+double *poisson_esds(double *intensities, ReflItemList *items,
+                     double adu_per_photon)
+{
+	int i;
+	double *esds = new_list_intensity();
+
+	for ( i=0; i<num_items(items); i++ ) {
+
+		struct refl_item *it;
+		signed int h, k, l;
+		double intensity, sigma;
+
+		it = get_item(items, i);
+		h = it->h;  k = it->k;  l = it->l;
+
+		intensity = lookup_intensity(intensities, h, k, l);
+
+		if ( intensity > 0.0 ) {
+			sigma = adu_per_photon * sqrt(intensity/adu_per_photon);
+		} else {
+			sigma = 0.0;
+		}
+
+		set_intensity(esds, h, k, l, sigma);
+
+	}
+
+	return esds;
+}
+
+
 void write_reflections(const char *filename, ReflItemList *items,
-                       double *intensities, double *phases,
-                       unsigned int *counts, UnitCell *cell,
-                       double adu_per_photon)
+                       double *intensities, double *esds, double *phases,
+                       unsigned int *counts, UnitCell *cell)
 {
 	FILE *fh;
 	int i;
@@ -79,8 +109,8 @@ void write_reflections(const char *filename, ReflItemList *items,
 			s = 0.0;
 		}
 
-		if ( intensity > 0.0 ) {
-			sigma = adu_per_photon * sqrt(intensity/adu_per_photon);
+		if ( esds != NULL ) {
+			sigma = lookup_intensity(esds, h, k, l);
 		} else {
 			sigma = 0.0;
 		}
@@ -90,9 +120,6 @@ void write_reflections(const char *filename, ReflItemList *items,
 		        h, k, l, intensity, ph, sigma, s/1.0e9, N);
 
 	}
-	STATUS("Warning: Errors have been estimated from Poisson distribution"
-	       " assuming %5.2f ADU per photon.\n", adu_per_photon);
-	STATUS("This is not necessarily a useful estimate.\n");
 	fclose(fh);
 }
 
diff --git a/src/reflections.h b/src/reflections.h
index 7b5e5e42..20060f21 100644
--- a/src/reflections.h
+++ b/src/reflections.h
@@ -21,10 +21,12 @@
 #include "utils.h"
 
 
+extern double *poisson_esds(double *intensities, ReflItemList *items,
+                            double adu_per_photon);
+
 extern void write_reflections(const char *filename, ReflItemList *items,
-                              double *intensities, double *phases,
-                              unsigned int *counts, UnitCell *cell,
-                              double adu_per_photon);
+                              double *intensities, double *esds, double *phases,
+                              unsigned int *counts, UnitCell *cell);
 
 extern ReflItemList *read_reflections(const char *filename,
                                       double *intensities, double *phases,