diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/partialator.c | 15 | ||||
| -rw-r--r-- | src/rejection.c | 90 | ||||
| -rw-r--r-- | src/rejection.h | 2 |
3 files changed, 96 insertions, 11 deletions
diff --git a/src/partialator.c b/src/partialator.c index dc0d647f..d5e1aef5 100644 --- a/src/partialator.c +++ b/src/partialator.c @@ -301,6 +301,7 @@ static void show_duds(Crystal **crystals, int n_crystals) int n_noref = 0; int n_solve = 0; int n_early = 0; + int n_cc = 0; for ( j=0; j<n_crystals; j++ ) { int flag; @@ -327,6 +328,10 @@ static void show_duds(Crystal **crystals, int n_crystals) n_early++; break; + case 6: + n_cc++; + break; + default: STATUS("Unknown flag %i\n", flag); break; @@ -339,6 +344,7 @@ static void show_duds(Crystal **crystals, int n_crystals) STATUS(" %i not enough reflections.\n", n_noref); STATUS(" %i solve failed.\n", n_solve); STATUS(" %i early rejection.\n", n_early); + STATUS(" %i bad CC.\n", n_cc); } } @@ -741,14 +747,14 @@ int main(int argc, char *argv[]) } /* Make a first pass at cutting out crap */ -// STATUS("Checking patterns.\n"); -// early_rejection(crystals, n_crystals); + STATUS("Checking patterns.\n"); + early_rejection(crystals, n_crystals); /* Make initial estimates */ full = merge_intensities(crystals, n_crystals, nthreads, pmodel, min_measurements, push_res); - check_rejection(crystals, n_crystals); + check_rejection(crystals, n_crystals, full); show_duds(crystals, n_crystals); @@ -773,14 +779,13 @@ int main(int argc, char *argv[]) init_free_dev, final_free_dev); show_duds(crystals, n_crystals); - check_rejection(crystals, n_crystals); + check_rejection(crystals, n_crystals, full); /* Re-estimate all the full intensities */ reflist_free(full); full = merge_intensities(crystals, n_crystals, nthreads, pmodel, min_measurements, push_res); - check_rejection(crystals, n_crystals); write_pgraph(full, crystals, n_crystals, i+1); } diff --git a/src/rejection.c b/src/rejection.c index 0b59e813..28a8d32f 100644 --- a/src/rejection.c +++ b/src/rejection.c @@ -33,10 +33,12 @@ #include <stdlib.h> #include <assert.h> +#include <gsl/gsl_statistics.h> #include "crystal.h" #include "reflist.h" #include "rejection.h" +#include "cell-utils.h" static double mean_intensity(RefList *list) @@ -88,20 +90,98 @@ void early_rejection(Crystal **crystals, int n) STATUS("Mean intensity/peak = %f ADU\n", m/n); STATUS("%i crystals flagged\n", n_flag); - check_rejection(crystals, n); + check_rejection(crystals, n, NULL); } -void check_rejection(Crystal **crystals, int n) +static void check_cc(Crystal *cr, RefList *full) +{ + RefList *list = crystal_get_reflections(cr); + Reflection *refl; + RefListIterator *iter; + double G = crystal_get_osf(cr); + double B = crystal_get_Bfac(cr); + UnitCell *cell = crystal_get_cell(cr); + double *vec1, *vec2; + int n, i; + double cc; + + n = num_reflections(list); + vec1 = malloc(n*sizeof(double)); + vec2 = malloc(n*sizeof(double)); + if ( (vec1 == NULL) || (vec2 == NULL) ) { + ERROR("Not enough memory to check CCs\n"); + return; + } + + i = 0; + for ( refl = first_refl(list, &iter); + refl != NULL; + refl = next_refl(refl, iter) ) + { + signed int h, k, l; + double pobs, pcalc; + double res, corr, Ipart; + Reflection *match; + + if ( !get_flag(refl) ) continue; /* Not free-flagged */ + + get_indices(refl, &h, &k, &l); + res = resolution(cell, h, k, l); + if ( 2.0*res > crystal_get_resolution_limit(cr) ) continue; + + match = find_refl(full, h, k, l); + if ( match == NULL ) continue; + + /* Calculated partiality */ + pcalc = get_partiality(refl); + + /* Observed partiality */ + corr = exp(-G) * exp(B*res*res) * get_lorentz(refl); + Ipart = get_intensity(refl) * corr; + pobs = Ipart / get_intensity(match); + + vec1[i] = pobs; + vec2[i] = pcalc; + i++; + } + + cc = gsl_stats_correlation(vec1, 1, vec2, 1, i); + //printf("%f\n", cc); + if ( cc < 0.5 ) crystal_set_user_flag(cr, 6); + + free(vec1); + free(vec2); +} + + +static void check_ccs(Crystal **crystals, int n_crystals, RefList *full) +{ + int i; + + for ( i=0; i<n_crystals; i++ ) { + check_cc(crystals[i], full); + } +} + + +void check_rejection(Crystal **crystals, int n, RefList *full) { int i; int n_acc = 0; + /* Check according to CCs FIXME: Disabled */ + //if ( full != NULL ) check_ccs(crystals, n, full); + for ( i=0; i<n; i++ ) { - if ( crystal_get_user_flag(crystals[i]) == 0 ) { - n_acc++; - if ( n_acc >= 2 ) break; + + /* Reject if B factor modulus is very large */ + if ( fabs(crystal_get_Bfac(crystals[i])) > 1e18 ) { + crystal_set_user_flag(crystals[i], 1); } + + if ( crystal_get_user_flag(crystals[i]) == 0 ) n_acc++; + } if ( n_acc < 2 ) { diff --git a/src/rejection.h b/src/rejection.h index fb73bd41..ec529941 100644 --- a/src/rejection.h +++ b/src/rejection.h @@ -38,6 +38,6 @@ #include "crystal.h" extern void early_rejection(Crystal **crystals, int n); -extern void check_rejection(Crystal **crystals, int n); +extern void check_rejection(Crystal **crystals, int n, RefList *full); #endif /* REJECTION_H */ |