aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--libcrystfel/src/fom.c220
-rw-r--r--libcrystfel/src/fom.h7
-rw-r--r--src/compare_hkl.c214
3 files changed, 236 insertions, 205 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,
diff --git a/src/compare_hkl.c b/src/compare_hkl.c
index 1d82f01c..1c7865d8 100644
--- a/src/compare_hkl.c
+++ b/src/compare_hkl.c
@@ -342,7 +342,7 @@ int main(int argc, char *argv[])
char *sym_str_fromfile1 = NULL;
char *sym_str_fromfile2 = NULL;
SymOpList *sym;
- int ncom, nrej, nmul, nneg, nres, nbij, ncen;
+ int ncom;
RefList *list1_acc;
RefList *list2_acc;
RefList *list1;
@@ -354,8 +354,6 @@ int main(int argc, char *argv[])
float rmin_fix = -1.0;
float rmax_fix = -1.0;
double rmin, rmax;
- Reflection *refl1;
- RefListIterator *iter;
float sigma_cutoff = -INFINITY;
int config_ignorenegs = 0;
int config_zeronegs = 0;
@@ -367,6 +365,7 @@ int main(int argc, char *argv[])
double max_I = -INFINITY;
float highres, lowres;
int mul_cutoff = 0;
+ int anom;
/* Long options */
const struct option longopts[] = {
@@ -658,215 +657,20 @@ int main(int argc, char *argv[])
reflist_free(list1_raw);
reflist_free(list2_raw);
- /* Select reflections to be used */
- ncom = 0;
- nrej = 0;
- nmul = 0;
- nneg = 0;
- nres = 0;
- nbij = 0;
- ncen = 0;
list1_acc = reflist_new();
list2_acc = reflist_new();
- 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 ( config_ignorenegs && ((val1 < 0.0) || (val2 < 0.0)) ) {
- nneg++;
- continue;
- }
-
- if ( (mul1 < mul_cutoff) || (mul2 < mul_cutoff) ) {
- nmul++;
- continue;
- }
-
- if ( config_zeronegs ) {
- 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++;
-
- }
-
+ anom = ( (fom == FOM_CCANO) || (fom == FOM_CRDANO)
+ || (fom == FOM_RANO) || (fom == FOM_RANORSPLIT) );
+ ncom = fom_select_reflections(list1, list2, list1_acc, list2_acc,
+ cell, sym,
+ anom, rmin_fix, rmax_fix, sigma_cutoff,
+ config_ignorenegs, config_zeronegs,
+ mul_cutoff, &min_I, &max_I);
reflist_free(list1);
reflist_free(list2);
- /* For anomalous figures of merit, we additionally require that we have
- * all the Bijvoet pairs after the above rejection tests */
- if ( (fom == FOM_CCANO) || (fom == FOM_CRDANO)
- || (fom == FOM_RANO) || (fom == FOM_RANORSPLIT) )
- {
- 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++;
- }
- }
-
gsl_set_error_handler_off();
- if ( nrej > 0 ) {
- STATUS("Discarded %i reflection pairs because either or both"
- " versions had I/sigma(I) < %f.\n", nrej, sigma_cutoff);
- }
-
- if ( config_ignorenegs && (nneg > 0) ) {
- STATUS("Discarded %i reflection pairs because either or both"
- " versions had negative intensities.\n", nneg);
- }
-
- if ( config_zeronegs && (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);
- }
-
STATUS("%i reflection pairs accepted.\n", ncom);
resolution_limits(list1_acc, cell, &rmin, &rmax);
generated by cgit v1.2.3 (git 2.25.1) at 2024-10-01 06:38:37 +0000