From c68ee0a5c446023b7f550d723c8f3e8b109b9bc5 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 22 Sep 2010 17:35:06 +0200 Subject: compare_hkl: Reject reflections with low SNR --- src/compare_hkl.c | 29 +++++++++++++++++++++++++++-- src/get_hkl.c | 6 ++++-- src/indexamajig.c | 3 ++- src/pattern_sim.c | 3 ++- src/process_hkl.c | 2 +- src/reflections.c | 5 ++++- src/reflections.h | 2 +- src/render_hkl.c | 2 +- src/utils.h | 5 +++++ 9 files changed, 47 insertions(+), 10 deletions(-) (limited to 'src') diff --git a/src/compare_hkl.c b/src/compare_hkl.c index e10e9567..f2e9a184 100644 --- a/src/compare_hkl.c +++ b/src/compare_hkl.c @@ -142,6 +142,10 @@ int main(int argc, char *argv[]) ReflItemList *i1, *i2, *icommon; int config_luzzati = 0; char *pdb = NULL; + double *esd1; + double *esd2; + int rej1 = 0; + int rej2 = 0; /* Long options */ const struct option longopts[] = { @@ -202,13 +206,15 @@ int main(int argc, char *argv[]) free(pdb); ref1 = new_list_intensity(); - i1 = read_reflections(afile, ref1, NULL, NULL); + esd1 = new_list_sigma(); + i1 = read_reflections(afile, ref1, NULL, NULL, esd1); if ( i1 == NULL ) { ERROR("Couldn't open file '%s'\n", afile); return 1; } ref2 = new_list_intensity(); - i2 = read_reflections(bfile, ref2, NULL, NULL); + esd2 = new_list_sigma(); + i2 = read_reflections(bfile, ref2, NULL, NULL, esd2); if ( i2 == NULL ) { ERROR("Couldn't open file '%s'\n", bfile); return 1; @@ -226,6 +232,8 @@ int main(int argc, char *argv[]) signed int h, k, l; signed int he, ke, le; double val1, val2; + double sig1, sig2; + int ig = 0; it = get_item(i1, i); h = it->h; k = it->k; l = it->l; @@ -238,12 +246,29 @@ int main(int argc, char *argv[]) val1 = lookup_intensity(ref1, h, k, l); val2 = lookup_intensity(ref2, he, ke, le); + sig1 = lookup_sigma(esd1, h, k, l); + sig2 = lookup_sigma(esd2, h, k, l); + + if ( val1 < 3.0 * sig1 ) { + rej1++; + ig = 1; + } + if ( val2 < 3.0 * sig2 ) { + rej2++; + ig = 1; + } + if ( ig ) continue; + set_intensity(ref2_transformed, h, k, l, val2); set_intensity(out, h, k, l, val1/val2); add_item(icommon, h, k, l); } ncom = num_items(icommon); + STATUS("%i reflections with I < 3.0*sigma(I) rejected from '%s'\n", + rej1, afile); + STATUS("%i reflections with I < 3.0*sigma(I) rejected from '%s'\n", + rej2, bfile); STATUS("%i,%i reflections: %i in common\n", num_items(i1), num_items(i2), ncom); diff --git a/src/get_hkl.c b/src/get_hkl.c index 8c7db771..a86b1450 100644 --- a/src/get_hkl.c +++ b/src/get_hkl.c @@ -238,7 +238,8 @@ int main(int argc, char *argv[]) phases, input_items); } else { ideal_ref = new_list_intensity(); - input_items = read_reflections(input, ideal_ref, phases, NULL); + input_items = read_reflections(input, ideal_ref, phases, + NULL, NULL); free(input); } @@ -275,7 +276,8 @@ int main(int argc, char *argv[]) /* Write out only reflections which are in the template * (and which we have in the input) */ ReflItemList *template_items; - template_items = read_reflections(template, NULL, NULL, NULL); + template_items = read_reflections(template, + NULL, NULL, NULL, NULL); write_items = intersection_items(input_items, template_items); delete_items(template_items); } else { diff --git a/src/indexamajig.c b/src/indexamajig.c index 281294f3..8163465b 100644 --- a/src/indexamajig.c +++ b/src/indexamajig.c @@ -597,7 +597,8 @@ int main(int argc, char *argv[]) if ( intfile != NULL ) { ReflItemList *items; - items = read_reflections(intfile, intensities, NULL, NULL); + items = read_reflections(intfile, intensities, + NULL, NULL, NULL); delete_items(items); } else { intensities = NULL; diff --git a/src/pattern_sim.c b/src/pattern_sim.c index 0d414653..01383161 100644 --- a/src/pattern_sim.c +++ b/src/pattern_sim.c @@ -332,7 +332,8 @@ int main(int argc, char *argv[]) } intensities = new_list_intensity(); phases = new_list_phase(); - items = read_reflections(intfile, intensities, phases, NULL); + items = read_reflections(intfile, intensities, phases, + NULL, NULL); free(intfile); delete_items(items); } diff --git a/src/process_hkl.c b/src/process_hkl.c index a00d1e99..82bbbb74 100644 --- a/src/process_hkl.c +++ b/src/process_hkl.c @@ -696,7 +696,7 @@ int main(int argc, char *argv[]) if ( reference != NULL ) { reference_i = new_list_intensity(); reference_items = read_reflections(reference, reference_i, - NULL, NULL); + NULL, NULL, NULL); if ( reference_items == NULL ) { ERROR("Couldn't read '%s'\n", reference); return 1; diff --git a/src/reflections.c b/src/reflections.c index 9a1a2910..4c3cf527 100644 --- a/src/reflections.c +++ b/src/reflections.c @@ -102,7 +102,7 @@ void write_reflections(const char *filename, ReflItemList *items, */ ReflItemList *read_reflections(const char *filename, double *intensities, double *phases, - unsigned int *counts) + unsigned int *counts, double *esds) { FILE *fh; char *rval; @@ -163,6 +163,9 @@ ReflItemList *read_reflections(const char *filename, if ( counts != NULL ) { set_count(counts, h, k, l, cts); } + if ( esds != NULL ) { + set_sigma(esds, h, k, l, sigma); + } } while ( rval != NULL ); diff --git a/src/reflections.h b/src/reflections.h index de5feb72..84345260 100644 --- a/src/reflections.h +++ b/src/reflections.h @@ -27,7 +27,7 @@ extern void write_reflections(const char *filename, ReflItemList *items, extern ReflItemList *read_reflections(const char *filename, double *intensities, double *phases, - unsigned int *counts); + unsigned int *counts, double *esds); extern double *ideal_intensities(double complex *sfac); diff --git a/src/render_hkl.c b/src/render_hkl.c index b6deac2b..5b27e4d5 100644 --- a/src/render_hkl.c +++ b/src/render_hkl.c @@ -609,7 +609,7 @@ int main(int argc, char *argv[]) } ref = new_list_intensity(); cts = new_list_count(); - ReflItemList *items = read_reflections(infile, ref, NULL, cts); + ReflItemList *items = read_reflections(infile, ref, NULL, cts, NULL); if ( ref == NULL ) { ERROR("Couldn't open file '%s'\n", infile); return 1; diff --git a/src/utils.h b/src/utils.h index dfd1913f..a8461d6a 100644 --- a/src/utils.h +++ b/src/utils.h @@ -175,6 +175,11 @@ static inline double angle_between(double x1, double y1, double z1, #define TYPE unsigned int #include "list_tmp.h" +/* As above, but for sigmas */ +#define LABEL(x) x##_sigma +#define TYPE double +#include "list_tmp.h" + /* ----------- Reflection lists indexed by sequence (not indices) ----------- */ -- cgit v1.2.3