diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/process_hkl.c | 74 |
1 files changed, 68 insertions, 6 deletions
diff --git a/src/process_hkl.c b/src/process_hkl.c index 594e4ed7..a34613c7 100644 --- a/src/process_hkl.c +++ b/src/process_hkl.c @@ -58,6 +58,7 @@ static void show_help(const char *s) " the default.\n" " -g, --histogram=<h,k,l> Calculate the histogram of measurements for this\n" " reflection.\n" +" --rmerge Calculate and report Rmerge and Rpim\n" "\n" " --scale Scale each pattern for best fit with the current\n" " model.\n" @@ -205,7 +206,8 @@ static void merge_pattern(double *model, ReflItemList *observed, ReflItemList *twins, const char *holo, const char *mero, double *hist_vals, signed int hist_h, signed int hist_k, - signed int hist_l, int *hist_n) + signed int hist_l, int *hist_n, double *devs, + double *tots, double *means) { int i; int twin; @@ -249,6 +251,13 @@ static void merge_pattern(double *model, ReflItemList *observed, } } + if ( tots != NULL ) { + double m; + m = lookup_intensity(means, h, k, l); + integrate_intensity(tots, h, k, l, intensity); + integrate_intensity(devs, h, k, l, fabs(intensity-m)); + } + /* Already seen this reflection in this pattern? Complain. */ if ( !find_item(sym_items, h, k, l) ) { /* Add the asymmetric version of this reflection to our @@ -285,7 +294,7 @@ static void merge_all(FILE *fh, double **pmodel, ReflItemList **pobserved, ReflItemList *twins, const char *holo, const char *mero, int n_total_patterns, double *hist_vals, signed int hist_h, signed int hist_k, signed int hist_l, - int *hist_i) + int *hist_i, double *devs, double *tots, double *means) { char *rval; float f0; @@ -334,7 +343,7 @@ static void merge_all(FILE *fh, double **pmodel, ReflItemList **pobserved, counts, config_maxonly, twins, holo, mero, hist_vals, hist_h, hist_k, hist_l, - hist_i); + hist_i, devs, tots, means); if ( n_patterns == config_stopafter ) break; @@ -430,6 +439,7 @@ int main(int argc, char *argv[]) int config_stopafter = 0; int config_sum = 0; int config_scale = 0; + int config_rmerge = 0; unsigned int n_total_patterns; char *sym = NULL; char *pdb = NULL; @@ -455,6 +465,7 @@ int main(int argc, char *argv[]) {"symmetry", 1, NULL, 'y'}, {"pdb", 1, NULL, 'p'}, {"histogram", 1, NULL, 'g'}, + {"rmerge", 0, &config_rmerge, 1}, {0, 0, NULL, 0} }; @@ -500,6 +511,12 @@ int main(int argc, char *argv[]) } + if ( config_sum && config_rmerge ) { + ERROR("Options --sum and --rmerge do not make sense" + " together.\n"); + return 1; + } + if ( filename == NULL ) { ERROR("Please specify filename using the -i option\n"); return 1; @@ -578,11 +595,9 @@ int main(int argc, char *argv[]) merge_all(fh, &model, &observed, &counts, config_maxonly, config_scale, config_sum, config_stopafter, twins, holo, sym, n_total_patterns, - hist_vals, hist_h, hist_k, hist_l, &hist_i); + hist_vals, hist_h, hist_k, hist_l, &hist_i, NULL, NULL, NULL); rewind(fh); - fclose(fh); - if ( hist_vals != NULL ) { STATUS("%i %i %i was seen %i times.\n", hist_h, hist_k, hist_l, hist_i); @@ -593,6 +608,53 @@ int main(int argc, char *argv[]) write_reflections(output, observed, model, NULL, counts, cell); } + if ( config_rmerge ) { + + double *devs = new_list_intensity(); + double *tots = new_list_intensity(); + double total_dev = 0.0; + double total_tot = 0.0; + double total_pdev = 0.0; + double total_rdev = 0.0; + int i; + + STATUS("Extra pass to calculate figures of merit...\n"); + + rewind(fh); + merge_all(fh, &model, &observed, &counts, + config_maxonly, config_scale, 0, + config_stopafter, twins, holo, sym, n_total_patterns, + NULL, 0, 0, 0, NULL, devs, tots, model); + + for ( i=0; i<num_items(observed); i++ ) { + + struct refl_item *it; + signed int h, k, l; + double dev; + unsigned int count; + + it = get_item(observed, i); + h = it->h; k = it->k, l = it->l; + + dev = lookup_intensity(devs, h, k, l); + count = lookup_count(counts, h, k, l); + + if ( count < 2 ) continue; + + total_dev += dev; + total_pdev += sqrt(1.0/(count-1.0))*dev; + total_rdev += sqrt(count/(count-1.0))*dev; + total_tot += lookup_intensity(tots, h, k, l); + + } + + STATUS("Rmerge = %f%%\n", 100.0*total_dev/total_tot); + STATUS(" Rpim = %f%%\n", 100.0*total_pdev/total_tot); + STATUS(" Rrim = %f%%\n", 100.0*total_rdev/total_tot); + + } + + fclose(fh); free(sym); free(model); free(counts); |