diff options
Diffstat (limited to 'libcrystfel/src')
| -rw-r--r-- | libcrystfel/src/fom.c | 220 | ||||
| -rw-r--r-- | libcrystfel/src/fom.h | 7 |
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, |