aboutsummaryrefslogtreecommitdiff
path: root/src/process_hkl.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/process_hkl.c')
-rw-r--r--src/process_hkl.c99
1 files changed, 42 insertions, 57 deletions
diff --git a/src/process_hkl.c b/src/process_hkl.c
index c05adbe9..2533102a 100644
--- a/src/process_hkl.c
+++ b/src/process_hkl.c
@@ -220,7 +220,7 @@ static void merge_pattern(double *model, ReflItemList *observed,
double *hist_vals,
signed int hist_h, signed int hist_k,
signed int hist_l, int *hist_n, double *devs,
- double *tots, double *means, FILE *outfh)
+ double *means, FILE *outfh)
{
int i;
int twin;
@@ -275,11 +275,11 @@ static void merge_pattern(double *model, ReflItemList *observed,
}
}
- if ( tots != NULL ) {
+ if ( devs != 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));
+ integrate_intensity(devs, h, k, l,
+ pow(intensity-m, 2.0));
}
if ( !find_item(sym_items, h, k, l) ) {
@@ -399,8 +399,7 @@ static void merge_all(FILE *fh, double **pmodel, ReflItemList **pobserved,
ReflItemList *reference, double *reference_i,
double *hist_vals,
signed int hist_h, signed int hist_k, signed int hist_l,
- int *hist_i, double *devs, double *tots, double *means,
- FILE *outfh)
+ int *hist_i, double *devs, double *means, FILE *outfh)
{
char *rval;
float f0;
@@ -470,7 +469,7 @@ static void merge_all(FILE *fh, double **pmodel, ReflItemList **pobserved,
twins, holo, mero,
reference, reference_i,
hist_vals, hist_h, hist_k, hist_l,
- hist_i, devs, tots, means, outfh);
+ hist_i, devs, means, outfh);
n_patterns++;
if ( n_patterns == config_stopafter ) break;
@@ -570,6 +569,8 @@ int main(int argc, char *argv[])
FILE *outfh = NULL;
struct beam_params *beam = NULL;
int space_for_hist = 0;
+ double *devs;
+ double *esds;
/* Long options */
const struct option longopts[] = {
@@ -771,7 +772,7 @@ int main(int argc, char *argv[])
config_startafter, config_stopafter,
twins, in_sym, sym, n_total_patterns,
reference_items, reference_i,
- hist_vals, hist_h, hist_k, hist_l, &hist_i, NULL, NULL, NULL,
+ hist_vals, hist_h, hist_k, hist_l, &hist_i, NULL, NULL,
outfh);
rewind(fh);
if ( space_for_hist && (hist_i >= space_for_hist) ) {
@@ -784,67 +785,51 @@ int main(int argc, char *argv[])
plot_histogram(hist_vals, hist_i);
}
- if ( output != NULL ) {
-
- double adu_per_photon;
-
- if ( beam == NULL ) {
- adu_per_photon = 167.0;
- STATUS("No beam parameters file provided (use -b), "
- "so I'm assuming 167.0 ADU per photon.\n");
- } else {
- adu_per_photon = beam->adu_per_photon;
- }
-
- write_reflections(output, observed, model, NULL, counts, cell,
- adu_per_photon);
- }
-
- if ( config_rmerge ) {
+ STATUS("Extra pass to calculate ESDs...\n");
+ devs = new_list_intensity();
+ esds = new_list_intensity();
+ rewind(fh);
+ merge_all(fh, &model, &observed, &counts,
+ config_maxonly, config_scale, 0,
+ config_startafter, config_stopafter,
+ twins, in_sym, sym,
+ n_total_patterns, reference_items, reference_i,
+ NULL, 0, 0, 0, NULL, devs, model, NULL);
- 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;
+ for ( i=0; i<num_items(observed); i++ ) {
- STATUS("Extra pass to calculate figures of merit...\n");
+ struct refl_item *it;
+ signed int h, k, l;
+ double dev, esd;
+ unsigned int count;
- rewind(fh);
- merge_all(fh, &model, &observed, &counts,
- config_maxonly, config_scale, 0,
- config_startafter, config_stopafter,
- twins, in_sym, sym,
- n_total_patterns, reference_items, reference_i,
- NULL, 0, 0, 0, NULL, devs, tots, model, NULL);
+ it = get_item(observed, i);
+ h = it->h; k = it->k, l = it->l;
- for ( i=0; i<num_items(observed); i++ ) {
+ dev = lookup_intensity(devs, h, k, l);
+ count = lookup_count(counts, h, k, l);
- struct refl_item *it;
- signed int h, k, l;
- double dev;
- unsigned int count;
+ if ( count < 2 ) continue;
- it = get_item(observed, i);
- h = it->h; k = it->k, l = it->l;
+ esd = sqrt(dev) / (double)count;
+ set_intensity(esds, h, k, l, esd);
- dev = lookup_intensity(devs, h, k, l);
- count = lookup_count(counts, h, k, l);
+ }
- if ( count < 2 ) continue;
+ if ( output != NULL ) {
- 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);
+ double adu_per_photon;
+ if ( beam == NULL ) {
+ adu_per_photon = 167.0;
+ STATUS("No beam parameters file provided (use -b), "
+ "so I'm assuming 167.0 ADU per photon.\n");
+ } else {
+ adu_per_photon = beam->adu_per_photon;
}
- STATUS("Rmerge = %f%%\n", 100.0*total_dev/total_tot);
- STATUS(" Rpim = %f%%\n", 100.0*total_pdev/total_tot);
- STATUS(" Rmeas = %f%%\n", 100.0*total_rdev/total_tot);
+ write_reflections(output, observed, model, esds,
+ NULL, counts, cell);
}