From 6445bd6887ac815d07f38692a98d0b1dd166fed6 Mon Sep 17 00:00:00 2001
From: Thomas White <taw@physics.org>
Date: Wed, 6 Jan 2021 15:03:26 +0100
Subject: partialator: Reject crystals with obviously too large profile radii

The criterion for "too large" is 20% of the 1/d value for the lowest
reflection which is not systematically absent according to the
centering.

A profile radius larger than the 1/d value for a reflection will crash
the xsphere partiality model, and some visualisation shows that this is
a clearly non-physical situation.  The profile radius shouldn't be
anywhere near the inter-Bragg spacing for reasonable data.

However, feedback shows that this is happening quite often in real data,
probably due to bad indexing.
---
 libcrystfel/src/cell-utils.c | 33 ++++++++++++++++++++++++++++++++-
 libcrystfel/src/cell-utils.h |  4 +++-
 2 files changed, 35 insertions(+), 2 deletions(-)

(limited to 'libcrystfel')

diff --git a/libcrystfel/src/cell-utils.c b/libcrystfel/src/cell-utils.c
index 22610ce1..93dea212 100644
--- a/libcrystfel/src/cell-utils.c
+++ b/libcrystfel/src/cell-utils.c
@@ -3,7 +3,7 @@
  *
  * Unit Cell utility functions
  *
- * Copyright © 2012-2020 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2021 Deutsches Elektronen-Synchrotron DESY,
  *                       a research centre of the Helmholtz Association.
  * Copyright © 2012 Lorenzo Galli
  *
@@ -1312,6 +1312,37 @@ double cell_get_volume(UnitCell *cell)
 }
 
 
+/**
+ * \param cell: A %UnitCell
+ *
+ * \returns the value of 1/d for the lowest order reflection
+ * that is not systematically absent according to the centering.
+ *
+ */
+double lowest_reflection(UnitCell *cell)
+{
+	signed int h, k, l;
+	double lowres = INFINITY;
+
+	/* FIXME: Inelegant and nasty.  Anyone want to work out
+	 * all the possible cases? */
+	for ( h=0; h<4; h++ ) {
+		for ( k=0; k<4; k++ ) {
+			for ( l=0; l<4; l++ ) {
+				if ( (h==0) && (k==0) && (l==0) ) continue;
+				if ( !forbidden_reflection(cell, h, k, l) ) {
+					double r = resolution(cell, h, k, l);
+					if ( r < lowres ) {
+						lowres = r;
+					}
+				}
+			}
+		}
+	}
+	return lowres;
+}
+
+
 /* Return true if the two centering symbols are identical,
  * or if they are a pair of R/P, which should be considered the
  * same for the purposes of cell comparison */
diff --git a/libcrystfel/src/cell-utils.h b/libcrystfel/src/cell-utils.h
index 00563ee6..0a110bf2 100644
--- a/libcrystfel/src/cell-utils.h
+++ b/libcrystfel/src/cell-utils.h
@@ -3,7 +3,7 @@
  *
  * Unit Cell utility functions
  *
- * Copyright © 2012-2020 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2021 Deutsches Elektronen-Synchrotron DESY,
  *                       a research centre of the Helmholtz Association.
  * Copyright © 2012 Lorenzo Galli
  *
@@ -79,6 +79,8 @@ extern int forbidden_reflection(UnitCell *cell,
 
 extern double cell_get_volume(UnitCell *cell);
 
+extern double lowest_reflection(UnitCell *cell);
+
 extern int compare_cell_parameters(UnitCell *cell, UnitCell *reference,
                                    const double *tols);
 
-- 
cgit v1.2.3