From 89bf2f6de2be2a276d8c349a1fa6ff6c63b5e7fe Mon Sep 17 00:00:00 2001
From: Thomas White <taw@physics.org>
Date: Fri, 16 Nov 2018 15:07:30 +0100
Subject: partialator: Handle partial reflections properly in deltaCChalf

---
 src/rejection.c | 32 ++++++++++++++++++++++++--------
 1 file changed, 24 insertions(+), 8 deletions(-)

(limited to 'src/rejection.c')

diff --git a/src/rejection.c b/src/rejection.c
index c07287b3..05acb0db 100644
--- a/src/rejection.c
+++ b/src/rejection.c
@@ -40,6 +40,7 @@
 #include "rejection.h"
 #include "cell-utils.h"
 #include "post-refinement.h"
+#include "merge.h"
 
 
 static double mean_intensity(RefList *list)
@@ -235,15 +236,19 @@ static double calculate_cchalf(RefList *template, RefList *full,
 			G = crystal_get_osf(exclude);
 			B = crystal_get_Bfac(exclude);
 
-			/* Total (multiplicative) correction factor */
-			Ii *= 1.0/G * exp(B*res*res) * get_lorentz(exrefl)
-			          / get_partiality(exrefl);
+			if ( get_partiality(exrefl) > MIN_PART_MERGE ) {
+
+				/* Total (multiplicative) correction factor */
+				Ii *= 1.0/G * exp(B*res*res) * get_lorentz(exrefl)
+				          / get_partiality(exrefl);
+
+				/* Remove contribution of this reflection */
+				Ex -= Ii - K;
+				Ex2 -= (Ii - K)*(Ii - K);
 
-			/* Remove contribution of this reflection */
-			Ex -= Ii - K;
-			Ex2 -= (Ii - K)*(Ii - K);
+				n_removed++;
 
-			n_removed++;
+			}
 
 			exrefl = next_found_refl(exrefl);
 		}
@@ -282,6 +287,7 @@ static void check_deltacchalf(Crystal **crystals, int n, RefList *full)
 	double *vals;
 	double mean, sd;
 	int nref = 0;
+	int nnan = 0;
 
 	calculate_refl_mean_var(full);
 
@@ -305,15 +311,25 @@ static void check_deltacchalf(Crystal **crystals, int n, RefList *full)
 		//STATUS("  Delta = %f  ", (cchalf - cchalfi)*100.0);
 		//STATUS("(nref = %i)\n", nref);
 		vals[i] = cchalf - cchalfi;
+		if ( isnan(vals[i]) || isinf(vals[i]) ) {
+			vals[i] = 0.0;
+			nnan++;
+		}
 		progress_bar(i, n-1, "Calculating deltaCChalf");
 	}
+	if ( nnan > 0 ) {
+		STATUS("WARNING: %i NaN or inf deltaCChalf values were "
+		       "replaced with zero\n", nnan);
+	}
 
 	mean = gsl_stats_mean(vals, 1, n);
 	sd = gsl_stats_sd_m(vals, 1, n, mean);
 	STATUS("deltaCChalf = %f ± %f %%\n", mean*100.0, sd*100.0);
 
 	for ( i=0; i<n; i++ ) {
-		if ( (vals[i]<0.0) && (vals[i] < mean-2.0*sd) ) {
+		if ( isnan(vals[i]) || isinf(vals[i])
+		  || ((vals[i]<0.0) && (vals[i] < mean-2.0*sd)) )
+		{
 			crystal_set_user_flag(crystals[i], PRFLAG_DELTACCHALF);
 		}
 	}
-- 
cgit v1.2.3