diff options
| -rw-r--r-- | libcrystfel/src/geometry.c | 5 | ||||
| -rw-r--r-- | src/process_image.c | 73 |
2 files changed, 76 insertions, 2 deletions
diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c index 8952486d..250f889c 100644 --- a/libcrystfel/src/geometry.c +++ b/libcrystfel/src/geometry.c @@ -151,6 +151,7 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst, double pr; double D; double del; + double rlowuc, rhighuc; /* Values before clamping */ /* Don't predict 000 */ if ( abs(h)+abs(k)+abs(l) == 0 ) return NULL; @@ -197,6 +198,8 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst, D = rlow - rhigh; + rlowuc = rlow; rhighuc = rhigh; + /* If the "lower" Ewald sphere is a long way away, use the * position at which the Ewald sphere would just touch the * reflection. @@ -243,7 +246,7 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst, set_detector_pos(refl, 0.0, xda, yda); } - set_partial(refl, rlow, rhigh, part, clamp_low, clamp_high); + set_partial(refl, rlowuc, rhighuc, part, clamp_low, clamp_high); set_lorentz(refl, 1.0); set_symmetric_indices(refl, h, k, l); set_redundancy(refl, 1); diff --git a/src/process_image.c b/src/process_image.c index 1d1dac82..f4392efe 100644 --- a/src/process_image.c +++ b/src/process_image.c @@ -34,6 +34,8 @@ #include <stdlib.h> #include <hdf5.h> #include <gsl/gsl_errno.h> +#include <gsl/gsl_statistics_double.h> +#include <gsl/gsl_sort.h> #include <unistd.h> #include "utils.h" @@ -50,6 +52,63 @@ #include "integration.h" +static int cmpd2(const void *av, const void *bv) +{ + double *ap, *bp; + double a, b; + + ap = (double *)av; + bp = (double *)bv; + + a = ap[1]; + b = bp[1]; + + if ( a > b ) return -1; + return 1; +} + + +static void refine_radius(Crystal *cr) +{ + Reflection *refl; + RefListIterator *iter; + FILE *fh; + double vals[num_reflections(crystal_get_reflections(cr))*2]; + int n = 0; + int i; + + fh = fopen("graph.dat", "w"); + + for ( refl = first_refl(crystal_get_reflections(cr), &iter); + refl != NULL; + refl = next_refl(refl, iter) ) + { + double i = get_intensity(refl); + double rlow, rhigh, p; + int cl, ch; + + get_partial(refl, &rlow, &rhigh, &p, &cl, &ch); + fprintf(fh, "%e %10f %e %e %f %i %i\n", (rhigh+rlow)/2.0, i, + rlow, rhigh, p, cl, ch); + + vals[(2*n)+0] = i; + vals[(2*n)+1] = (rhigh+rlow)/2.0; + n++; + + } + + qsort(vals, n, sizeof(double)*2, cmpd2); + + for ( i=0; i<n-1; i++ ) { + double mean = gsl_stats_mean(vals, 2, i); + double var = gsl_stats_variance_m(vals, 2, i, mean); + printf("%.2f %i\n", mean/var, i); + } + + fclose(fh); +} + + void process_image(const struct index_args *iargs, struct pattern_args *pargs, Stream *st, int cookie, const char *tmpdir, int results_pipe) { @@ -169,7 +228,7 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs, /* Default parameters */ image.div = 0.0; - image.bw = 0.001; + image.bw = 0.00000001; STATUS("Warning: div, bw and pr are hardcoded in this version.\n"); for ( i=0; i<image.n_crystals; i++ ) { crystal_set_profile_radius(image.crystals[i], 0.01e9); @@ -183,6 +242,18 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs, iargs->int_diag, iargs->int_diag_h, iargs->int_diag_k, iargs->int_diag_l, results_pipe); + for ( i=0; i<image.n_crystals; i++ ) { + refine_radius(image.crystals[i]); + reflist_free(crystal_get_reflections(image.crystals[i])); + } + + integrate_all_4(&image, iargs->int_meth, PMODEL_SCSPHERE, + iargs->push_res, + iargs->ir_inn, iargs->ir_mid, iargs->ir_out, + iargs->int_diag, iargs->int_diag_h, + iargs->int_diag_k, iargs->int_diag_l, + results_pipe); + write_chunk(st, &image, hdfile, iargs->stream_peaks, iargs->stream_refls, pargs->filename_p_e->ev); |