diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/rejection.c | 50 |
1 files changed, 38 insertions, 12 deletions
diff --git a/src/rejection.c b/src/rejection.c index 1dccd2da..eca55f2f 100644 --- a/src/rejection.c +++ b/src/rejection.c @@ -99,14 +99,17 @@ static double calculate_cchalf(RefList *template, RefList *full, { Reflection *trefl; RefListIterator *iter; - double all_sum_mean = 0.0; - double all_sumsq_mean = 0.0; - double all_sum_var = 0.0; double hsig2E, sig2Y; int n = 0; signed int oh = 0; signed int ok = 0; signed int ol = 0; + double wSum = 0.0; + double wSum2 = 0.0; + double mean = 0.0; + double S = 0.0; + double all_sum_var = 0.0; + long int total_contribs = 0; /* Iterate over all reflections */ for ( trefl = first_refl(template, &iter); @@ -115,8 +118,8 @@ static double calculate_cchalf(RefList *template, RefList *full, { struct reflection_contributions *c; int j; - double refl_sum = 0.0; - double refl_sumsq = 0.0; + double refl_sum; + double refl_sumsq; double refl_mean, refl_var; signed int h, k, l; int nc = 0; @@ -136,6 +139,8 @@ static double calculate_cchalf(RefList *template, RefList *full, c = get_contributions(refl); assert(c != NULL); + /* Mean of contributions */ + refl_sum = 0.0; for ( j=0; j<c->n_contrib; j++ ) { double Ii; @@ -148,7 +153,6 @@ static double calculate_cchalf(RefList *template, RefList *full, Ii = get_intensity(c->contribs[j]); refl_sum += Ii; - refl_sumsq += Ii * Ii; nc++; } @@ -156,19 +160,41 @@ static double calculate_cchalf(RefList *template, RefList *full, if ( nc < 2 ) continue; refl_mean = refl_sum / nc; - refl_var = (refl_sumsq/nc - refl_mean*refl_mean); - refl_var *= nc/(nc-1.0); - all_sum_mean += refl_mean; - all_sumsq_mean += refl_mean*refl_mean; + /* Variance of contributions */ + refl_sumsq = 0.0; + for ( j=0; j<c->n_contrib; j++ ) { + + double Ii; + + if ( c->contrib_crystals[j] == exclude ) { + continue; + } + + /* FIXME: Apply corrections */ + Ii = get_intensity(c->contribs[j]); + + refl_sumsq += (Ii-refl_mean)*(Ii-refl_mean); + /* (nc already summed above) */ + + } + + refl_var = refl_sumsq/(nc-1.0); all_sum_var += refl_var; n++; + total_contribs += nc; + + double w = 1.0; + wSum += w; + wSum2 += w*w; + double meanOld = mean; + mean = meanOld + (w/wSum) * (refl_mean - meanOld); + S += w * (refl_mean - meanOld) * (refl_mean - mean); } hsig2E = all_sum_var / n; - sig2Y = (all_sumsq_mean/n - all_sum_mean*all_sum_mean/(n*n)); - sig2Y *= n/(n-1.0); + sig2Y = S / (wSum - 1.0); if ( pnref != NULL ) { *pnref = n; |