aboutsummaryrefslogtreecommitdiff
path: root/libcrystfel/src
diff options
context:
space:
mode:
Diffstat (limited to 'libcrystfel/src')
-rw-r--r--libcrystfel/src/fom.c220
-rw-r--r--libcrystfel/src/fom.h7
2 files changed, 227 insertions, 0 deletions
diff --git a/libcrystfel/src/fom.c b/libcrystfel/src/fom.c
index e10df5b8..35d5c5e3 100644
--- a/libcrystfel/src/fom.c
+++ b/libcrystfel/src/fom.c
@@ -789,3 +789,223 @@ struct fom_context *fom_calculate(RefList *list1, RefList *list2, UnitCell *cell
return fctx;
}
+
+int fom_select_reflections(RefList *list1, RefList *list2,
+ RefList *list1_acc, RefList *list2_acc,
+ UnitCell *cell, SymOpList *sym,
+ int anom, double rmin_fix, double rmax_fix,
+ double sigma_cutoff, int ignore_negs,
+ int zero_negs, int mul_cutoff,
+ double *pmin_I, double *pmax_I)
+{
+ Reflection *refl1;
+ RefListIterator *iter;
+ int ncom, nrej, nmul, nneg, nres, nbij, ncen;
+ double min_I = +INFINITY;
+ double max_I = -INFINITY;
+
+ /* Select reflections to be used */
+ ncom = 0;
+ nrej = 0;
+ nmul = 0;
+ nneg = 0;
+ nres = 0;
+ nbij = 0;
+ ncen = 0;
+ for ( refl1 = first_refl(list1, &iter);
+ refl1 != NULL;
+ refl1 = next_refl(refl1, iter) )
+ {
+ signed int h, k, l;
+ double val1, val2;
+ double esd1, esd2;
+ int mul1, mul2;
+ Reflection *refl2;
+ Reflection *refl1_acc;
+ Reflection *refl2_acc;
+
+ get_indices(refl1, &h, &k, &l);
+
+ refl2 = find_refl(list2, h, k, l);
+ if ( refl2 == NULL ) continue;
+
+ val1 = get_intensity(refl1);
+ val2 = get_intensity(refl2);
+
+ esd1 = get_esd_intensity(refl1);
+ esd2 = get_esd_intensity(refl2);
+
+ mul1 = get_redundancy(refl1);
+ mul2 = get_redundancy(refl2);
+
+ if ( (val1 < sigma_cutoff * esd1)
+ || (val2 < sigma_cutoff * esd2) )
+ {
+ nrej++;
+ continue;
+ }
+
+ if ( ignore_negs && ((val1 < 0.0) || (val2 < 0.0)) ) {
+ nneg++;
+ continue;
+ }
+
+ if ( (mul1 < mul_cutoff) || (mul2 < mul_cutoff) ) {
+ nmul++;
+ continue;
+ }
+
+ if ( zero_negs ) {
+ int d = 0;
+ if ( val1 < 0.0 ) {
+ val1 = 0.0;
+ d = 1;
+ }
+ if ( val2 < 0.0 ) {
+ val2 = 0.0;
+ d = 1;
+ }
+ if ( d ) nneg++;
+ }
+
+ if ( rmin_fix > 0.0 ) {
+ double res = 2.0*resolution(cell, h, k, l);
+ if ( res < rmin_fix ) {
+ nres++;
+ continue;
+ }
+ }
+
+ if ( rmax_fix > 0.0 ) {
+ double res = 2.0*resolution(cell, h, k, l);
+ if ( res > rmax_fix ) {
+ nres++;
+ continue;
+ }
+ }
+
+ refl1_acc = add_refl(list1_acc, h, k, l);
+ copy_data(refl1_acc, refl1);
+ set_intensity(refl1_acc, val1);
+
+ refl2_acc = add_refl(list2_acc, h, k, l);
+ copy_data(refl2_acc, refl2);
+ set_intensity(refl2_acc, val2);
+
+ if ( val1 > max_I ) max_I = val1;
+ if ( val1 < min_I ) min_I = val1;
+
+ ncom++;
+
+ }
+
+ /* For anomalous figures of merit, we additionally require that we have
+ * all the Bijvoet pairs after the above rejection tests */
+ if ( anom ) {
+
+ list1 = list1_acc;
+ list2 = list2_acc;
+ list1_acc = reflist_new();
+ list2_acc = reflist_new();
+
+ min_I = +INFINITY;
+ max_I = -INFINITY;
+ ncom = 0;
+
+ for ( refl1 = first_refl(list1, &iter);
+ refl1 != NULL;
+ refl1 = next_refl(refl1, iter) )
+ {
+ Reflection *refl1_bij = NULL;
+ Reflection *refl2_bij = NULL;
+ signed int h, k, l;
+ signed int hb, kb, lb;
+ Reflection *refl1_acc;
+ Reflection *refl2_acc;
+ Reflection *refl2;
+ double val1, val2;
+
+ get_indices(refl1, &h, &k, &l);
+
+ refl2 = find_refl(list2, h, k, l);
+ assert(refl2 != NULL);
+
+ val1 = get_intensity(refl1);
+ val2 = get_intensity(refl2);
+
+ if ( is_centric(h, k, l, sym) ) {
+ ncen++;
+ continue;
+ }
+
+ if ( find_equiv_in_list(list1, -h, -k, -l, sym,
+ &hb, &kb, &lb) )
+ {
+ refl1_bij = find_refl(list1, hb, kb, lb);
+ }
+
+ if ( find_equiv_in_list(list2, -h, -k, -l, sym,
+ &hb, &kb, &lb) )
+ {
+ refl2_bij = find_refl(list2, hb, kb, lb);
+ }
+
+ if ( (refl1_bij == NULL) || (refl2_bij == NULL) ) {
+ nbij++;
+ continue;
+ }
+
+ refl1_acc = add_refl(list1_acc, h, k, l);
+ copy_data(refl1_acc, refl1);
+ set_intensity(refl1_acc, val1);
+
+ refl2_acc = add_refl(list2_acc, h, k, l);
+ copy_data(refl2_acc, refl2);
+ set_intensity(refl2_acc, val2);
+
+ if ( val1 > max_I ) max_I = val1;
+ if ( val1 < min_I ) min_I = val1;
+
+ ncom++;
+ }
+ }
+
+ if ( nrej > 0 ) {
+ STATUS("Discarded %i reflection pairs because either or both"
+ " versions had I/sigma(I) < %f.\n", nrej, sigma_cutoff);
+ }
+
+ if ( ignore_negs && (nneg > 0) ) {
+ STATUS("Discarded %i reflection pairs because either or both"
+ " versions had negative intensities.\n", nneg);
+ }
+
+ if ( zero_negs && (nneg > 0) ) {
+ STATUS("For %i reflection pairs, either or both versions had"
+ " negative intensities which were set to zero.\n", nneg);
+ }
+
+ if ( nmul > 0 ) {
+ STATUS("%i reflection pairs rejected because either or both"
+ " versions had too few measurements.\n", nmul);
+ }
+
+ if ( nres > 0 ) {
+ STATUS("%i reflection pairs rejected because either or both"
+ " versions were outside the resolution range.\n", nres);
+ }
+
+ if ( nbij > 0 ) {
+ STATUS("%i reflection pairs rejected because either or both"
+ " versions did not have Bijvoet partners.\n", nres);
+ }
+
+ if ( ncen > 0 ) {
+ STATUS("%i reflection pairs rejected because they were"
+ " centric.\n", ncen);
+ }
+
+ *pmin_I = min_I;
+ *pmax_I = max_I;
+ return ncom;
+}
diff --git a/libcrystfel/src/fom.h b/libcrystfel/src/fom.h
index 8af434e8..1adb940f 100644
--- a/libcrystfel/src/fom.h
+++ b/libcrystfel/src/fom.h
@@ -86,6 +86,13 @@ struct fom_context
int *n_within;
};
+extern int fom_select_reflections(RefList *list1, RefList *list2,
+ RefList *list1_acc, RefList *list2_acc,
+ UnitCell *cell, SymOpList *sym,
+ int anom, double rmin_fix, double rmax_fix,
+ double sigma_cutoff, int ignore_negs,
+ int zero_negs, int mul_cutoff,
+ double *pmin_I, double *pmax_I);
extern struct fom_context *fom_calculate(RefList *list1, RefList *list2,
UnitCell *cell,