aboutsummaryrefslogtreecommitdiff
path: root/src/process_hkl.c
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2013-05-15 14:34:47 +0200
committerThomas White <taw@physics.org>2013-05-15 14:34:47 +0200
commite8c8ed31440d008b34ed5cbfaa9646a6c1eac671 (patch)
tree55c6d2d183ba197fb093cb8b31ecfa55e37c4e64 /src/process_hkl.c
parent527528e4b4566193ca1bdca4e8b3cfc73e01965b (diff)
process_hkl: Use weighting according to sigmas, and use better variance formula
Diffstat (limited to 'src/process_hkl.c')
-rw-r--r--src/process_hkl.c49
1 files changed, 26 insertions, 23 deletions
diff --git a/src/process_hkl.c b/src/process_hkl.c
index 47ca2feb..ffbc83f0 100644
--- a/src/process_hkl.c
+++ b/src/process_hkl.c
@@ -201,13 +201,14 @@ static int merge_crystal(RefList *model, struct image *image, Crystal *cr,
refl != NULL;
refl = next_refl(refl, iter) )
{
- double intensity;
+ double refl_intensity, refl_sigma;
double xl, yl, zl;
double pol, pa, pb, phi, tt, ool;
signed int h, k, l;
- int cur_redundancy;
- double cur_intensity, cur_sumsq;
+ int model_redundancy;
Reflection *model_version;
+ double w;
+ double temp, delta, R, mean, M2, sumweight;
get_indices(refl, &h, &k, &l);
@@ -219,8 +220,11 @@ static int merge_crystal(RefList *model, struct image *image, Crystal *cr,
model_version = add_refl(model, h, k, l);
}
- intensity = scale * get_intensity(refl);
+ refl_intensity = scale * get_intensity(refl);
+ refl_sigma = scale * get_esd_intensity(refl);
+ w = pow(refl_sigma, -2.0);
+ /* FIXME: Should go before determining the scaling factor */
if ( !config_nopolar ) {
/* Polarisation correction assuming 100% polarisation
@@ -235,23 +239,28 @@ static int merge_crystal(RefList *model, struct image *image, Crystal *cr,
pa = pow(sin(phi)*sin(tt), 2.0);
pb = pow(cos(tt), 2.0);
pol = 1.0 - 2.0*(1.0-pa) + (1.0+pb);
- intensity /= pol;
+ refl_intensity /= pol;
}
- cur_intensity = get_intensity(model_version);
- set_intensity(model_version, cur_intensity + intensity);
+ mean = get_intensity(model_version);
+ sumweight = get_temp1(model_version);
+ M2 = get_temp2(model_version);
- cur_redundancy = get_redundancy(model_version);
- set_redundancy(model_version, cur_redundancy+1);
+ temp = w + sumweight;
+ delta = refl_intensity - mean;
+ R = delta * w / temp;
+ set_intensity(model_version, mean + R);
+ set_temp2(model_version, M2 + sumweight * delta * R);
+ set_temp1(model_version, temp);
- cur_sumsq = get_temp1(model_version);
- set_temp1(model_version, cur_sumsq + pow(intensity, 2.0));
+ model_redundancy = get_redundancy(model_version);
+ set_redundancy(model_version, ++model_redundancy);
if ( hist_vals != NULL ) {
if ( (h==hist_h) && (k==hist_k) && (l==hist_l) ) {
- hist_vals[*hist_n] = intensity;
+ hist_vals[*hist_n] = refl_intensity;
*hist_n += 1;
}
@@ -343,7 +352,7 @@ static int merge_all(Stream *st, RefList *model, RefList *reference,
refl != NULL;
refl = next_refl(refl, iter) )
{
- double intensity, sumsq, esd;
+ double var;
int red;
red = get_redundancy(refl);
@@ -352,13 +361,9 @@ static int merge_all(Stream *st, RefList *model, RefList *reference,
continue;
}
- intensity = get_intensity(refl) / red;
- set_intensity(refl, intensity);
-
- sumsq = get_temp1(refl) / red;
- esd = sqrt(sumsq - pow(intensity, 2.0)) / sqrt(red);
- set_esd_intensity(refl, esd);
-
+ var = get_temp2(refl) / get_temp1(refl);
+ var *= (double)red/(red-1);
+ set_esd_intensity(refl, sqrt(var)/sqrt(red));
}
return 0;
@@ -572,9 +577,7 @@ int main(int argc, char *argv[])
STATUS("Extra pass for scaling...\n");
- reference = copy_reflist(model);
-
- reflist_free(model);
+ reference = model;
model = reflist_new();
r = merge_all(st, model, reference, sym,