process_hkl: Implement Rick's ESD calculation

1 files changed, 42 insertions, 57 deletions

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);
 
 	}