diff options
Diffstat (limited to 'src/merge.c')
| -rw-r--r-- | src/merge.c | 128 |
1 files changed, 121 insertions, 7 deletions
diff --git a/src/merge.c b/src/merge.c index 17e57c90..14460953 100644 --- a/src/merge.c +++ b/src/merge.c @@ -50,7 +50,7 @@ /* Minimum partiality of a reflection for it to be merged */ -#define MIN_PART_MERGE (0.05) +#define MIN_PART_MERGE (0.3) struct merge_queue_args @@ -59,7 +59,6 @@ struct merge_queue_args pthread_rwlock_t full_lock; Crystal **crystals; int n_started; - PartialityModel pmodel; double push_res; int use_weak; long long int n_reflections; @@ -159,7 +158,6 @@ static void run_merge_job(void *vwargs, int cookie) signed int h, k, l; double mean, sumweight, M2, temp, delta, R; double corr, res, w; - //double esd; if ( get_partiality(refl) < MIN_PART_MERGE ) continue; @@ -189,7 +187,7 @@ static void run_merge_job(void *vwargs, int cookie) } /* Total (multiplicative) correction factor */ - corr = exp(-G) * exp(B*res*res) * get_lorentz(refl) + corr = 1.0/G * exp(B*res*res) * get_lorentz(refl) / get_partiality(refl); if ( isnan(corr) ) { ERROR("NaN corr:\n"); @@ -200,7 +198,7 @@ static void run_merge_job(void *vwargs, int cookie) continue; } - //esd = get_esd_intensity(refl) * corr; + /* Reflections count less the more they have to be scaled up */ w = 1.0; /* Running mean and variance calculation */ @@ -229,7 +227,7 @@ static void finalise_merge_job(void *vqargs, void *vwargs) RefList *merge_intensities(Crystal **crystals, int n, int n_threads, - PartialityModel pmodel, int min_meas, + int min_meas, double push_res, int use_weak) { RefList *full; @@ -245,7 +243,6 @@ RefList *merge_intensities(Crystal **crystals, int n, int n_threads, qargs.full = full; qargs.n_started = 0; qargs.crystals = crystals; - qargs.pmodel = pmodel; qargs.push_res = push_res; qargs.use_weak = use_weak; qargs.n_reflections = 0; @@ -285,3 +282,120 @@ RefList *merge_intensities(Crystal **crystals, int n, int n_threads, reflist_free(full); return full2; } + + +double residual(Crystal *cr, const RefList *full, int free, + int *pn_used, const char *filename) +{ + Reflection *refl; + RefListIterator *iter; + int n_used = 0; + double num = 0.0; + double den = 0.0; + double G = crystal_get_osf(cr); + double B = crystal_get_Bfac(cr); + UnitCell *cell = crystal_get_cell(cr); + + for ( refl = first_refl(crystal_get_reflections(cr), &iter); + refl != NULL; + refl = next_refl(refl, iter) ) + { + double p, w, corr, res; + signed int h, k, l; + Reflection *match; + double I_full; + double int1, int2; + + if ( free && !get_flag(refl) ) continue; + + get_indices(refl, &h, &k, &l); + res = resolution(cell, h, k, l); + match = find_refl(full, h, k, l); + if ( match == NULL ) continue; + I_full = get_intensity(match); + + if ( get_redundancy(match) < 2 ) continue; + + p = get_partiality(refl); + //if ( p < 0.2 ) continue; + + corr = get_lorentz(refl) / (G * exp(-B*res*res)); + int1 = get_intensity(refl) * corr; + int2 = p*I_full; + w = p; + + num += fabs(int1 - int2) * w; + den += fabs(int1 + int2) * w/2.0; + + n_used++; + + } + + if ( pn_used != NULL ) *pn_used = n_used; + return num/(den*sqrt(2)); +} + + +double log_residual(Crystal *cr, const RefList *full, int free, + int *pn_used, const char *filename) +{ + double dev = 0.0; + double G, B; + Reflection *refl; + RefListIterator *iter; + int n_used = 0; + FILE *fh = NULL; + + G = crystal_get_osf(cr); + B = crystal_get_Bfac(cr); + if ( filename != NULL ) { + fh = fopen(filename, "a"); + if ( fh == NULL ) { + ERROR("Failed to open '%s'\n", filename); + } + } + + for ( refl = first_refl(crystal_get_reflections(cr), &iter); + refl != NULL; + refl = next_refl(refl, iter) ) + { + double p, L, s, w; + signed int h, k, l; + Reflection *match; + double esd, I_full, I_partial; + double fx; + + if ( free && !get_flag(refl) ) continue; + + get_indices(refl, &h, &k, &l); + match = find_refl(full, h, k, l); + if ( match == NULL ) continue; + + p = get_partiality(refl); + L = get_lorentz(refl); + I_partial = get_intensity(refl); + I_full = get_intensity(match); + esd = get_esd_intensity(refl); + s = resolution(crystal_get_cell(cr), h, k, l); + + if ( I_partial <= 3.0*esd ) continue; /* Also because of log */ + if ( get_redundancy(match) < 2 ) continue; + if ( I_full <= 0 ) continue; /* Because log */ + if ( p <= 0.0 ) continue; /* Because of log */ + + fx = log(G) - B*s*s + log(p) + log(I_full) - log(I_partial) - log(L); + w = 1.0; + dev += w*fx*fx; + + if ( fh != NULL ) { + fprintf(fh, "%4i %4i %4i %e %e\n", + h, k, l, s, dev); + } + + } + + if ( fh != NULL ) fclose(fh); + + if ( pn_used != NULL ) *pn_used = n_used; + return dev; +} |