aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/compare_hkl.c3
-rw-r--r--src/facetron.c3
-rw-r--r--src/get_hkl.c4
-rw-r--r--src/process_hkl.c99
-rw-r--r--src/reflections.c43
-rw-r--r--src/reflections.h8
6 files changed, 87 insertions, 73 deletions
diff --git a/src/compare_hkl.c b/src/compare_hkl.c
index 489e2b43..29dfe2a5 100644
--- a/src/compare_hkl.c
+++ b/src/compare_hkl.c
@@ -456,7 +456,8 @@ int main(int argc, char *argv[])
if ( outfile != NULL ) {
- write_reflections(outfile, icommon, out, NULL, NULL, cell, 1.0);
+ write_reflections(outfile, icommon, out, NULL,
+ NULL, NULL, cell);
STATUS("Sigma(I) values in output file are not meaningful.\n");
}
diff --git a/src/facetron.c b/src/facetron.c
index 53fb00f2..f4148ae1 100644
--- a/src/facetron.c
+++ b/src/facetron.c
@@ -449,8 +449,7 @@ int main(int argc, char *argv[])
}
/* Output results */
- write_reflections(outfile, obs, i_full, NULL, cts, NULL, 1.0);
- STATUS("Sigma(I) values in output file are not (yet) meaningful.\n");
+ write_reflections(outfile, obs, i_full, NULL, NULL, cts, NULL);
/* Clean up */
free(i_full);
diff --git a/src/get_hkl.c b/src/get_hkl.c
index 8a8a71c2..7fb14f7e 100644
--- a/src/get_hkl.c
+++ b/src/get_hkl.c
@@ -429,8 +429,8 @@ int main(int argc, char *argv[])
adu_per_photon = beam->adu_per_photon;
}
- write_reflections(output, write_items, ideal_ref, phases, NULL, cell,
- adu_per_photon);
+ write_reflections(output, write_items, ideal_ref, NULL, phases,
+ NULL, cell);
delete_items(input_items);
delete_items(write_items);
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);
}
diff --git a/src/reflections.c b/src/reflections.c
index 0474d517..98fc9bab 100644
--- a/src/reflections.c
+++ b/src/reflections.c
@@ -22,10 +22,40 @@
#include "beam-parameters.h"
+double *poisson_esds(double *intensities, ReflItemList *items,
+ double adu_per_photon)
+{
+ int i;
+ double *esds = new_list_intensity();
+
+ for ( i=0; i<num_items(items); i++ ) {
+
+ struct refl_item *it;
+ signed int h, k, l;
+ double intensity, sigma;
+
+ it = get_item(items, i);
+ h = it->h; k = it->k; l = it->l;
+
+ intensity = lookup_intensity(intensities, h, k, l);
+
+ if ( intensity > 0.0 ) {
+ sigma = adu_per_photon * sqrt(intensity/adu_per_photon);
+ } else {
+ sigma = 0.0;
+ }
+
+ set_intensity(esds, h, k, l, sigma);
+
+ }
+
+ return esds;
+}
+
+
void write_reflections(const char *filename, ReflItemList *items,
- double *intensities, double *phases,
- unsigned int *counts, UnitCell *cell,
- double adu_per_photon)
+ double *intensities, double *esds, double *phases,
+ unsigned int *counts, UnitCell *cell)
{
FILE *fh;
int i;
@@ -79,8 +109,8 @@ void write_reflections(const char *filename, ReflItemList *items,
s = 0.0;
}
- if ( intensity > 0.0 ) {
- sigma = adu_per_photon * sqrt(intensity/adu_per_photon);
+ if ( esds != NULL ) {
+ sigma = lookup_intensity(esds, h, k, l);
} else {
sigma = 0.0;
}
@@ -90,9 +120,6 @@ void write_reflections(const char *filename, ReflItemList *items,
h, k, l, intensity, ph, sigma, s/1.0e9, N);
}
- STATUS("Warning: Errors have been estimated from Poisson distribution"
- " assuming %5.2f ADU per photon.\n", adu_per_photon);
- STATUS("This is not necessarily a useful estimate.\n");
fclose(fh);
}
diff --git a/src/reflections.h b/src/reflections.h
index 7b5e5e42..20060f21 100644
--- a/src/reflections.h
+++ b/src/reflections.h
@@ -21,10 +21,12 @@
#include "utils.h"
+extern double *poisson_esds(double *intensities, ReflItemList *items,
+ double adu_per_photon);
+
extern void write_reflections(const char *filename, ReflItemList *items,
- double *intensities, double *phases,
- unsigned int *counts, UnitCell *cell,
- double adu_per_photon);
+ double *intensities, double *esds, double *phases,
+ unsigned int *counts, UnitCell *cell);
extern ReflItemList *read_reflections(const char *filename,
double *intensities, double *phases,
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-17 01:41:00 +0000