diff options
| -rw-r--r-- | src/rejection.c | 35 |
1 files changed, 26 insertions, 9 deletions
diff --git a/src/rejection.c b/src/rejection.c index 7b78803e..9f298c0b 100644 --- a/src/rejection.c +++ b/src/rejection.c @@ -124,6 +124,7 @@ static double calculate_cchalf(RefList *template, RefList *full, signed int h, k, l; int nc = 0; Reflection *refl; + double res; get_indices(trefl, &h, &k, &l); @@ -139,19 +140,29 @@ static double calculate_cchalf(RefList *template, RefList *full, c = get_contributions(refl); assert(c != NULL); + /* Calculate the resolution just once, using the cell from the + * first crystal to contribute, otherwise it takes too long */ + res = resolution(crystal_get_cell(c->contrib_crystals[0]), + h, k, l); + /* Mean of contributions */ refl_sum = 0.0; for ( j=0; j<c->n_contrib; j++ ) { - double Ii; + double Ii, G, B; if ( c->contrib_crystals[j] == exclude ) { continue; } - /* FIXME: Apply corrections */ + G = crystal_get_osf(c->contrib_crystals[j]); + B = crystal_get_Bfac(c->contrib_crystals[j]); Ii = get_intensity(c->contribs[j]); + /* Total (multiplicative) correction factor */ + Ii *= 1.0/G * exp(B*res*res) * get_lorentz(c->contribs[j]) + / get_partiality(c->contribs[j]); + refl_sum += Ii; nc++; @@ -164,15 +175,20 @@ static double calculate_cchalf(RefList *template, RefList *full, refl_sumsq = 0.0; for ( j=0; j<c->n_contrib; j++ ) { - double Ii; + double Ii, G, B; if ( c->contrib_crystals[j] == exclude ) { continue; } - /* FIXME: Apply corrections */ + G = crystal_get_osf(c->contrib_crystals[j]); + B = crystal_get_Bfac(c->contrib_crystals[j]); Ii = get_intensity(c->contribs[j]); + /* Total (multiplicative) correction factor */ + Ii *= 1.0/G * exp(B*res*res) * get_lorentz(c->contribs[j]) + / get_partiality(c->contribs[j]); + refl_sumsq += (Ii-refl_mean)*(Ii-refl_mean); /* (nc already summed above) */ @@ -226,12 +242,13 @@ static void check_deltacchalf(Crystal **crystals, int n, RefList *full) RefList *template = crystal_get_reflections(crystals[i]); cchalf = calculate_cchalf(template, full, NULL, NULL); cchalfi = calculate_cchalf(template, full, crystals[i], &nref); - STATUS("Frame %i:", i); - STATUS(" With = %f ", cchalf*100.0); - STATUS("Without = %f", cchalfi*100.0); - STATUS(" Delta = %f ", (cchalf - cchalfi)*100.0); - STATUS("(nref = %i)\n", nref); + //STATUS("Frame %i:", i); + //STATUS(" With = %f ", cchalf*100.0); + //STATUS("Without = %f", cchalfi*100.0); + //STATUS(" Delta = %f ", (cchalf - cchalfi)*100.0); + //STATUS("(nref = %i)\n", nref); vals[i] = cchalf - cchalfi; + progress_bar(i, n-1, "Calculating deltaCChalf"); } mean = gsl_stats_mean(vals, 1, n); |