Skip to content
GitLab
Commit a4d0acc0 authored by Thomas White's avatar Thomas White
Browse files

process_hkl: Implement Rick's ESD calculation

parent 757b17f4
Hide whitespace changes
Inline Side-by-side
src/compare_hkl.c
View file @ a4d0acc0
......@@ -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");
}
......
src/facetron.c
View file @ a4d0acc0
......@@ -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);
......
src/get_hkl.c
View file @ a4d0acc0
......@@ -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);
......
src/process_hkl.c
View file @ a4d0acc0
......@@ -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);
}
......
src/reflections.c
View file @ a4d0acc0
......@@ -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);
}
......
src/reflections.h
View file @ a4d0acc0
......@@ -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,
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

imprint | privacy | Conditions of Use