Don't select scalable reflections during post refinement

4 files changed, 61 insertions, 43 deletions

diff --git a/src/geometry.c b/src/geometry.c
index 9b937b24..e2454adf 100644
--- a/src/geometry.c
+++ b/src/geometry.c
@@ -275,45 +275,10 @@ RefList *find_intersections(struct image *image, UnitCell *cell)
 }
 
 
-/* Decide which reflections can be scaled */
-static void select_scalable_reflections(RefList *list, ReflItemList *sc_l)
-{
-	Reflection *refl;
-	RefListIterator *iter;
-
-	for ( refl = first_refl(list, &iter);
-	      refl != NULL;
-	      refl = next_refl(refl, iter) ) {
-
-		int scalable = 1;
-		double v;
-
-		if ( get_partiality(refl) < 0.1 ) scalable = 0;
-		v = fabs(get_intensity(refl));
-		if ( v < 0.1 ) scalable = 0;
-
-		set_scalable(refl, scalable);
-		if ( scalable && (sc_l != NULL) ) {
-
-			signed int h, k, l;
-
-			get_indices(refl, &h, &k, &l);  /* Should already be
-			                                 * asymmetric */
-			if ( (sc_l != NULL) && (!find_item(sc_l, h, k, l)) ) {
-				add_item(sc_l, h, k, l);
-			}
-
-		}
-
-	}
-}
-
-
 /* Calculate partialities and apply them to the image's raw_reflections,
  * while adding to a ReflItemList of the currentl scalable (asymmetric)
  * reflections. */
 void update_partialities(struct image *image, const char *sym,
-                         ReflItemList *scalable,
                          int *n_expected, int *n_found, int *n_notfound)
 {
 	Reflection *refl;
@@ -380,5 +345,4 @@ void update_partialities(struct image *image, const char *sym,
 	}
 
 	reflist_free(predicted);
-	select_scalable_reflections(image->reflections, scalable);
 }
diff --git a/src/geometry.h b/src/geometry.h
index 15c6e552..892ff747 100644
--- a/src/geometry.h
+++ b/src/geometry.h
@@ -22,7 +22,6 @@
 extern RefList *find_intersections(struct image *image, UnitCell *cell);
 
 extern void update_partialities(struct image *image, const char *sym,
-                                ReflItemList *scalable,
                                 int *n_expected, int *n_found, int *n_notfound);
 
 #endif	/* GEOMETRY_H */
diff --git a/src/partialator.c b/src/partialator.c
index 367d3a6e..167f736d 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -151,6 +151,46 @@ static void refine_all(struct image *images, int n_total_patterns,
 }
 
 
+/* Decide which reflections can be scaled */
+static int select_scalable_reflections(RefList *list, ReflItemList *sc_l)
+{
+	Reflection *refl;
+	RefListIterator *iter;
+	int nobs = 0;
+
+	for ( refl = first_refl(list, &iter);
+	      refl != NULL;
+	      refl = next_refl(refl, iter) ) {
+
+		int scalable = 1;
+		double v;
+
+		if ( get_partiality(refl) < 0.1 ) scalable = 0;
+		v = fabs(get_intensity(refl));
+		if ( v < 0.1 ) scalable = 0;
+		set_scalable(refl, scalable);
+
+		if ( scalable ) {
+
+			signed int h, k, l;
+
+			nobs++;
+
+			/* Add (asymmetric) indices to list */
+			get_indices(refl, &h, &k, &l);
+
+			if ( !find_item(sc_l, h, k, l) ) {
+				add_item(sc_l, h, k, l);
+			}
+
+		}
+
+	}
+
+	return nobs;
+}
+
+
 int main(int argc, char *argv[])
 {
 	int c;
@@ -173,6 +213,7 @@ int main(int argc, char *argv[])
 	int n_notfound = 0;
 	char *cref;
 	int n_usable_patterns = 0;
+	int nobs;
 	char *reference_file = NULL;
 	double *reference = NULL;
 	RefList *reference_list = NULL;
@@ -313,6 +354,7 @@ int main(int argc, char *argv[])
 	/* Fill in what we know about the images so far */
 	rewind(fh);
 	scalable = new_items();
+	nobs = 0;
 	for ( i=0; i<n_total_patterns; i++ ) {
 
 		RefList *as;
@@ -355,9 +397,11 @@ int main(int argc, char *argv[])
 		reflist_free(cur->reflections);
 		cur->reflections = as;
 
-		update_partialities(cur, sym, scalable,
+		update_partialities(cur, sym,
 		                    &n_expected, &n_found, &n_notfound);
 
+		nobs += select_scalable_reflections(cur->reflections, scalable);
+
 		progress_bar(i, n_total_patterns-1, "Loading pattern data");
 		n_usable_patterns++;
 
@@ -366,7 +410,7 @@ int main(int argc, char *argv[])
 	STATUS("Found %5.2f%% of the expected peaks (missed %i of %i).\n",
 	       100.0 * (double)n_found / n_expected, n_notfound, n_expected);
 	STATUS("Mean measurements per scalable unique reflection: %5.2f\n",
-	       (double)n_found / num_items(scalable));
+	       (double)nobs / num_items(scalable));
 
 	cref = find_common_reflections(images, n_usable_patterns);
 
@@ -438,6 +482,18 @@ int main(int argc, char *argv[])
 		refine_all(images, n_usable_patterns, det, sym, scalable,
 		           reference_list, nthreads, fhg, fhp);
 
+		nobs = 0;
+		clear_items(scalable);
+		for ( i=0; i<n_usable_patterns; i++ ) {
+
+			struct image *cur = &images[i];
+			nobs += select_scalable_reflections(cur->reflections,
+			                                    scalable);
+
+		}
+		STATUS("Mean measurements per scalable unique "
+		       "reflection: %5.2f\n", (double)nobs/num_items(scalable));
+
 		/* Re-estimate all the full intensities */
 		reflist_free(full);
 		full = scale_intensities(images, n_usable_patterns,
diff --git a/src/post-refinement.c b/src/post-refinement.c
index 46e320e6..7eef5dfa 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -486,7 +486,7 @@ void pr_refine(struct image *image, const RefList *full, const char *sym)
 	const int verbose = 1;
 	int nexp, nfound, nnotfound;
 
-	update_partialities(image, sym, NULL, &nexp, &nfound, &nnotfound);
+	update_partialities(image, sym, &nexp, &nfound, &nnotfound);
 
 	if ( verbose ) {
 		dev = mean_partial_dev(image, full, sym);
@@ -516,8 +516,7 @@ void pr_refine(struct image *image, const RefList *full, const char *sym)
 
 		max_shift = pr_iterate(image, full, sym);
 
-		update_partialities(image, sym, NULL,
-		                    &nexp, &nfound, &nnotfound);
+		update_partialities(image, sym, &nexp, &nfound, &nnotfound);
 
 		if ( verbose ) {
 			dev = mean_partial_dev(image, full, sym);
@@ -536,7 +535,7 @@ void pr_refine(struct image *image, const RefList *full, const char *sym)
 			cell_set_reciprocal(image->indexed_cell, asx, asy, asz,
 				                  bsx, bsy, bsz, csx, csy, csz);
 
-			update_partialities(image, sym, NULL,
+			update_partialities(image, sym,
 		                            &nexp, &nfound, &nnotfound);
 
 			image->pr_dud = 1;