diff options
| -rw-r--r-- | .gitignore | 2 | ||||
| -rw-r--r-- | src/Makefile.am | 7 | ||||
| -rw-r--r-- | src/integr_sim.c | 232 | ||||
| -rw-r--r-- | src/process_hkl.c | 268 | ||||
| -rw-r--r-- | src/utils.h | 2 |
5 files changed, 274 insertions, 237 deletions
@@ -10,5 +10,5 @@ stamp-h1 src/*.o src/.deps src/pattern_sim -src/integr_sim +src/process_hkl *~ diff --git a/src/Makefile.am b/src/Makefile.am index 47f6f002..efaf09a8 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -1,4 +1,4 @@ -bin_PROGRAMS = pattern_sim integr_sim +bin_PROGRAMS = pattern_sim process_hkl AM_CFLAGS = -Wall @@ -6,6 +6,5 @@ pattern_sim_SOURCES = pattern_sim.c diffraction.c utils.c image.c cell.c \ hdf5-file.c ewald.c detector.c sfac.c intensities.c pattern_sim_LDADD = @LIBS@ -integr_sim_SOURCES = integr_sim.c diffraction.c utils.c ewald.c sfac.c \ - statistics.c cell.c -integr_sim_LDADD = @LIBS@ +process_hkl_SOURCES = process_hkl.c sfac.c statistics.c cell.c utils.c +process_hkl_LDADD = @LIBS@ diff --git a/src/integr_sim.c b/src/integr_sim.c deleted file mode 100644 index 08788277..00000000 --- a/src/integr_sim.c +++ /dev/null @@ -1,232 +0,0 @@ -/* - * integr_sim.c - * - * Test integration of intensities - * - * (c) 2006-2009 Thomas White <thomas.white@desy.de> - * - * Part of CrystFEL - crystallography with a FEL - * - */ - - -#ifdef HAVE_CONFIG_H -#include <config.h> -#endif - -#include <stdarg.h> -#include <stdlib.h> -#include <stdio.h> -#include <string.h> -#include <unistd.h> - -#include "cell.h" -#include "image.h" -#include "utils.h" -#include "statistics.h" -#include "sfac.h" - - -static void main_show_help(const char *s) -{ - printf("Syntax: %s\n\n", s); - printf("Test relrod integration\n\n"); - printf(" -h Display this help message\n"); -} - - -static void write_RvsQ(const char *name, double *ref, double *trueref, - unsigned int *counts, double scale, UnitCell *cell) -{ - FILE *fh; - double smax, sbracket; - signed int h, k, l; - - fh = fopen(name, "w"); - - smax = 0.0; - for ( h=-INDMAX; h<INDMAX; h++ ) { - for ( k=-INDMAX; k<INDMAX; k++ ) { - for ( l=-INDMAX; l<INDMAX; l++ ) { - double s = resolution(cell, h, k, l); - if ( (lookup_count(counts, h, k, l) > 0) && (s > smax) ) { - smax = s; - } - } - } - } - - for ( sbracket=0.0; sbracket<smax; sbracket+=smax/10.0 ) { - - double top = 0.0; - double bot = 0.0; - int n = 0; - - for ( h=-INDMAX; h<INDMAX; h++ ) { - for ( k=-INDMAX; k<INDMAX; k++ ) { - for ( l=-INDMAX; l<INDMAX; l++ ) { - - double s; - int c; - c = lookup_count(counts, h, k, l); - s = resolution(cell, h, k, l); - if ((s>=sbracket) && (s<sbracket+smax/10.0) && (c>0)) { - - double obs, calc, obsi; - - obs = lookup_intensity(ref, h, k, l); - calc = lookup_intensity(trueref, h, k, l); - - obsi = obs / (double)c; - top += fabs(obsi/scale - calc); - bot += obsi/scale; - n++; - } - - } - } - } - - fprintf(fh, "%8.5f %8.5f %i\n", sbracket+smax/20.0, top/bot, n); - - } - fclose(fh); -} - - -int main(int argc, char *argv[]) -{ - int c; - int i; - struct image image; - double *ref, *trueref; - unsigned int *counts; - size_t ref_size; - signed int h, k, l; - double R, scale; - FILE *fh; - - while ((c = getopt(argc, argv, "h")) != -1) { - - switch ( c ) { - - case 'h' : { - main_show_help(argv[0]); - return 0; - } - - } - - } - - /* Define image parameters */ - image.width = 1024; - image.height = 1024; - image.fmode = FORMULATION_CLEN; - image.x_centre = 512.5; - image.y_centre = 512.5; - image.camera_len = 0.05; /* 5 cm (front CCD can move from 5cm-20cm) */ - image.resolution = 13333.3; /* 75 micron pixel size */ - image.xray_energy = eV_to_J(2.0e3); /* 2 keV energy */ - image.lambda = ph_en_to_lambda(image.xray_energy); /* Wavelength */ - image.molecule = load_molecule(); - - /* Prepare array for integrated intensities */ - ref = new_list_intensity(); - - /* Array for sample counts */ - counts = new_list_count(); - - /* Calculate true intensities */ - get_reflections_cached(image.molecule, image.xray_energy); - /* Complex structure factors now in image.molecule->reflections */ - - for ( i=1; i<=10e3; i++ ) { - - #if 0 - image.orientation = random_quaternion(); - - /* Calculate reflections using large smax - * (rather than the actual value) */ - //get_reflections(&image, cell, 1.0e9); - - nrefl = image_feature_count(image.rflist); - for ( j=0; j<nrefl; j++ ) { - - struct imagefeature *f; - double t, s, intensity, F; - - f = image_get_feature(image.rflist, j); - - t = 100.0e-9; /* Thickness 100 nm */ - s = f->s; /* Get excitation error */ - F = structure_factor(f->h, f->k, f->l); - - /* Calculate intensity from this reflection */ - intensity = pow( F * SINC(M_PI*t*s), 2.0); - - if ( intensity < 0.1 ) continue; - - if ( (f->h == 2) && (f->k == 2) && (f->l == 2) ) { - fprintf(fh1, "%f %f\n", s, intensity); - } - - if ( (f->h == 15) && (f->k == 15) && (f->l == 15) ) { - fprintf(fh2, "%f %f\n", s, intensity); - } - - integrate_reflection(ref, f->h, f->k, f->l, intensity); - add_count(counts, f->h, f->k, f->l, 1); - } - #endif - - if ( i % 1000 == 0 ) { - - int j; - double mean_counts; - int ctot = 0; - int nmeas = 0; - double ff; - char name[64]; - - for ( j=0; j<ref_size; j++ ) { - ctot += counts[j]; - if ( counts[j] > 0 ) nmeas++; - } - mean_counts = (double)ctot/nmeas; - - ff = lookup_intensity(ref, 2, 2, 2) - / lookup_count(counts, 2, 2, 2); - - R = stat_r2(ref, trueref, counts, ref_size, &scale); - printf("%8i: R=%5.2f%% (sf=%7.4f)" - " mean meas/refl=%5.2f," - " %i reflections measured. %f\n", - i, R*100.0, scale, mean_counts, nmeas, - ff); - - /* Record graph of R against q for this N */ - snprintf(name, 63, "R_vs_q-%i.dat", i); - write_RvsQ(name, ref, trueref, counts, - scale, image.molecule->cell); - - } - - } - - fh = fopen("reflections.dat", "w"); - for ( h=-INDMAX; h<INDMAX; h++ ) { - for ( k=-INDMAX; k<INDMAX; k++ ) { - for ( l=-INDMAX; l<INDMAX; l++ ) { - int N; - N = lookup_count(counts, h, k, l); - if ( N == 0 ) continue; - double F = lookup_intensity(ref, h, k, l) / N; - fprintf(fh, "%3i %3i %3i %f\n", h, k, l, F); - } - } - } - fclose(fh); - - return 0; -} diff --git a/src/process_hkl.c b/src/process_hkl.c new file mode 100644 index 00000000..097a3eb3 --- /dev/null +++ b/src/process_hkl.c @@ -0,0 +1,268 @@ +/* + * process_hkl.c + * + * Assemble and process FEL Bragg intensities + * + * (c) 2006-2009 Thomas White <thomas.white@desy.de> + * + * Part of CrystFEL - crystallography with a FEL + * + */ + + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include <stdarg.h> +#include <stdlib.h> +#include <stdio.h> +#include <string.h> +#include <unistd.h> +#include <getopt.h> + +#include "utils.h" +#include "statistics.h" +#include "sfac.h" + + +static void show_help(const char *s) +{ + printf("Syntax: %s [options]\n\n", s); + printf("Assemble and process FEL Bragg intensities.\n\n"); + printf(" -h Display this help message.\n"); + printf(" -i <filename> Specify input filename (\"-\" for stdin).\n"); +} + + +static void write_RvsQ(const char *name, double *ref, double *trueref, + unsigned int *counts, double scale, UnitCell *cell) +{ + FILE *fh; + double smax, sbracket; + signed int h, k, l; + + fh = fopen(name, "w"); + + smax = 0.0; + for ( h=-INDMAX; h<INDMAX; h++ ) { + for ( k=-INDMAX; k<INDMAX; k++ ) { + for ( l=-INDMAX; l<INDMAX; l++ ) { + double s = resolution(cell, h, k, l); + if ( (lookup_count(counts, h, k, l) > 0) && (s > smax) ) { + smax = s; + } + } + } + } + + for ( sbracket=0.0; sbracket<smax; sbracket+=smax/10.0 ) { + + double top = 0.0; + double bot = 0.0; + int n = 0; + + for ( h=-INDMAX; h<INDMAX; h++ ) { + for ( k=-INDMAX; k<INDMAX; k++ ) { + for ( l=-INDMAX; l<INDMAX; l++ ) { + + double s; + int c; + c = lookup_count(counts, h, k, l); + s = resolution(cell, h, k, l); + if ((s>=sbracket) && (s<sbracket+smax/10.0) && (c>0)) { + + double obs, calc, obsi; + + obs = lookup_intensity(ref, h, k, l); + calc = lookup_intensity(trueref, h, k, l); + + obsi = obs / (double)c; + top += fabs(obsi/scale - calc); + bot += obsi/scale; + n++; + } + + } + } + } + + fprintf(fh, "%8.5f %8.5f %i\n", sbracket+smax/20.0, top/bot, n); + + } + fclose(fh); +} + + +static void write_reflections(const char *filename, unsigned int *counts, + double *ref) +{ + FILE *fh; + signed int h, k, l; + + fh = fopen(filename, "w"); + for ( h=-INDMAX; h<INDMAX; h++ ) { + for ( k=-INDMAX; k<INDMAX; k++ ) { + for ( l=-INDMAX; l<INDMAX; l++ ) { + int N; + N = lookup_count(counts, h, k, l); + if ( N == 0 ) continue; + double F = lookup_intensity(ref, h, k, l) / N; + fprintf(fh, "%3i %3i %3i %f\n", h, k, l, F); + } + } + } + fclose(fh); +} + + +static double *ideal_intensities(double complex *sfac) +{ + double *ref; + signed int h, k, l; + + ref = new_list_intensity(); + + /* Generate ideal reflections from complex structure factors */ + for ( h=-INDMAX; h<=INDMAX; h++ ) { + for ( k=-INDMAX; k<=INDMAX; k++ ) { + for ( l=-INDMAX; l<=INDMAX; l++ ) { + double complex F = lookup_sfac(sfac, h, k, l); + double intensity = pow(cabs(F), 2.0); + set_intensity(ref, h, k, l, intensity); + } + } + } + + return ref; +} + + +static void process_reflections(double *ref, double *trueref, + unsigned int *counts, unsigned int n_patterns, + UnitCell *cell) +{ + int j; + double mean_counts; + int ctot = 0; + int nmeas = 0; + double ff; + char name[64]; + double R, scale; + + for ( j=0; j<LIST_SIZE; j++ ) { + ctot += counts[j]; + if ( counts[j] > 0 ) nmeas++; + } + mean_counts = (double)ctot/nmeas; + + ff = lookup_intensity(ref, 2, 2, 2) / lookup_count(counts, 2, 2, 2); + + R = stat_r2(ref, trueref, counts, LIST_SIZE, &scale); + STATUS("%8u: R=%5.2f%% (sf=%7.4e)" + " mean meas/refl=%5.2f," + " %i reflections measured. %f\n", + n_patterns, R*100.0, scale, mean_counts, nmeas, ff); + + /* Record graph of R against q for this N */ + snprintf(name, 63, "results/R_vs_q-%u.dat", n_patterns); + write_RvsQ(name, ref, trueref, counts, scale, cell); +} + + +int main(int argc, char *argv[]) +{ + int c; + char *filename = NULL; + FILE *fh; + unsigned int n_patterns; + double *ref, *trueref; + unsigned int *counts; + char *rval; + struct molecule *mol; + + /* Long options */ + const struct option longopts[] = { + {"help", 0, NULL, 'h'}, + {"input", 1, NULL, 'i'}, + {0, 0, NULL, 0} + }; + + /* Short options */ + while ((c = getopt_long(argc, argv, "hi:", longopts, NULL)) != -1) { + + switch (c) { + case 'h' : { + show_help(argv[0]); + return 0; + } + + case 'i' : { + filename = strdup(optarg); + break; + } + + case 0 : { + break; + } + + default : { + return 1; + } + } + + } + + if ( filename == NULL ) { + ERROR("Please specify filename using the -i option\n"); + return 1; + } + + mol = load_molecule(); + get_reflections_cached(mol, eV_to_J(2.0e3)); + + ref = new_list_intensity(); + counts = new_list_count(); + trueref = ideal_intensities(mol->reflections); + + if ( strcmp(filename, "-") == 0 ) { + fh = stdin; + } else { + fh = fopen(filename, "r"); + } + free(filename); + if ( fh == NULL ) { + ERROR("Failed to open input file\n"); + return 1; + } + + n_patterns = 0; + do { + + char line[1024]; + int h, k, l, intensity; + int r; + + rval = fgets(line, 1023, fh); + if ( strncmp(line, "New pattern", 11) == 0 ) { + n_patterns++; + if ( n_patterns % 1000 == 0 ) { + process_reflections(ref, trueref, counts, + n_patterns, mol->cell); + } + } + + r = sscanf(line, "%i %i %i %i", &h, &k, &l, &intensity); + if ( r != 4 ) continue; + + integrate_intensity(ref, h, k, l, intensity); + integrate_count(counts, h, k, l, 1); + + } while ( rval != NULL ); + + fclose(fh); + + write_reflections("results/reflections.hkl", counts, ref); + + return 0; +} diff --git a/src/utils.h b/src/utils.h index 1de206bf..2f31db57 100644 --- a/src/utils.h +++ b/src/utils.h @@ -131,6 +131,8 @@ static inline double distance3d(double x1, double y1, double z1, /* Array size */ #define IDIM (INDMAX*2 +1) +#define LIST_SIZE (IDIM*IDIM*IDIM) + /* Create functions for storing reflection intensities indexed as h,k,l */ #define LABEL(x) x##_intensity #define TYPE double |