diff options
| -rw-r--r-- | .gitignore | 1 | ||||
| -rw-r--r-- | README | 2 | ||||
| -rw-r--r-- | src/Makefile.am | 6 | ||||
| -rw-r--r-- | src/Makefile.in | 37 | ||||
| -rw-r--r-- | src/check_hkl.c | 367 | ||||
| -rw-r--r-- | src/compare_hkl.c | 65 |
6 files changed, 413 insertions, 65 deletions
@@ -19,4 +19,5 @@ src/facetron src/cubeit src/reintegrate src/estimate_background +src/check_hkl *~ @@ -55,6 +55,8 @@ In addition, there is also: - compare_hkl, for working out the differences (e.g. a q-dependent scaling factor) between two lists of reflections. + - check_hkl, for determining things like completeness. + - calibrate_detector, for summing patterns after peak detection to use for calibrating your detector. diff --git a/src/Makefile.am b/src/Makefile.am index dc048526..dbbdc9a2 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -1,6 +1,6 @@ bin_PROGRAMS = pattern_sim process_hkl get_hkl indexamajig compare_hkl \ powder_plot render_hkl calibrate_detector facetron cubeit \ - reintegrate estimate_background + reintegrate estimate_background check_hkl if HAVE_GTK bin_PROGRAMS += hdfsee @@ -46,6 +46,10 @@ compare_hkl_SOURCES = compare_hkl.c sfac.c cell.c utils.c reflections.c \ statistics.c symmetry.c compare_hkl_LDADD = @LIBS@ +check_hkl_SOURCES = check_hkl.c sfac.c cell.c utils.c reflections.c \ + statistics.c symmetry.c +check_hkl_LDADD = @LIBS@ + powder_plot_SOURCES = powder_plot.c cell.c utils.c image.c hdf5-file.c \ detector.c powder_plot_LDADD = @LIBS@ diff --git a/src/Makefile.in b/src/Makefile.in index 5805bef2..00f9af02 100644 --- a/src/Makefile.in +++ b/src/Makefile.in @@ -37,7 +37,7 @@ bin_PROGRAMS = pattern_sim$(EXEEXT) process_hkl$(EXEEXT) \ powder_plot$(EXEEXT) render_hkl$(EXEEXT) \ calibrate_detector$(EXEEXT) facetron$(EXEEXT) cubeit$(EXEEXT) \ reintegrate$(EXEEXT) estimate_background$(EXEEXT) \ - $(am__EXEEXT_1) + check_hkl$(EXEEXT) $(am__EXEEXT_1) @HAVE_GTK_TRUE@am__append_1 = hdfsee @HAVE_OPENCL_TRUE@am__append_2 = diffraction-gpu.c cl-utils.c @HAVE_OPENCL_TRUE@am__append_3 = diffraction-gpu.c cl-utils.c @@ -60,6 +60,11 @@ am_calibrate_detector_OBJECTS = calibrate_detector.$(OBJEXT) \ cell.$(OBJEXT) thread-pool.$(OBJEXT) calibrate_detector_OBJECTS = $(am_calibrate_detector_OBJECTS) calibrate_detector_DEPENDENCIES = +am_check_hkl_OBJECTS = check_hkl.$(OBJEXT) sfac.$(OBJEXT) \ + cell.$(OBJEXT) utils.$(OBJEXT) reflections.$(OBJEXT) \ + statistics.$(OBJEXT) symmetry.$(OBJEXT) +check_hkl_OBJECTS = $(am_check_hkl_OBJECTS) +check_hkl_DEPENDENCIES = am_compare_hkl_OBJECTS = compare_hkl.$(OBJEXT) sfac.$(OBJEXT) \ cell.$(OBJEXT) utils.$(OBJEXT) reflections.$(OBJEXT) \ statistics.$(OBJEXT) symmetry.$(OBJEXT) @@ -165,16 +170,18 @@ am__v_CCLD_0 = @echo " CCLD " $@; AM_V_GEN = $(am__v_GEN_$(V)) am__v_GEN_ = $(am__v_GEN_$(AM_DEFAULT_VERBOSITY)) am__v_GEN_0 = @echo " GEN " $@; -SOURCES = $(calibrate_detector_SOURCES) $(compare_hkl_SOURCES) \ - $(cubeit_SOURCES) $(estimate_background_SOURCES) \ - $(facetron_SOURCES) $(get_hkl_SOURCES) $(hdfsee_SOURCES) \ - $(indexamajig_SOURCES) $(pattern_sim_SOURCES) \ - $(powder_plot_SOURCES) $(process_hkl_SOURCES) \ - $(reintegrate_SOURCES) $(render_hkl_SOURCES) -DIST_SOURCES = $(calibrate_detector_SOURCES) $(compare_hkl_SOURCES) \ - $(cubeit_SOURCES) $(estimate_background_SOURCES) \ - $(facetron_SOURCES) $(get_hkl_SOURCES) \ - $(am__hdfsee_SOURCES_DIST) $(am__indexamajig_SOURCES_DIST) \ +SOURCES = $(calibrate_detector_SOURCES) $(check_hkl_SOURCES) \ + $(compare_hkl_SOURCES) $(cubeit_SOURCES) \ + $(estimate_background_SOURCES) $(facetron_SOURCES) \ + $(get_hkl_SOURCES) $(hdfsee_SOURCES) $(indexamajig_SOURCES) \ + $(pattern_sim_SOURCES) $(powder_plot_SOURCES) \ + $(process_hkl_SOURCES) $(reintegrate_SOURCES) \ + $(render_hkl_SOURCES) +DIST_SOURCES = $(calibrate_detector_SOURCES) $(check_hkl_SOURCES) \ + $(compare_hkl_SOURCES) $(cubeit_SOURCES) \ + $(estimate_background_SOURCES) $(facetron_SOURCES) \ + $(get_hkl_SOURCES) $(am__hdfsee_SOURCES_DIST) \ + $(am__indexamajig_SOURCES_DIST) \ $(am__pattern_sim_SOURCES_DIST) $(powder_plot_SOURCES) \ $(process_hkl_SOURCES) $(reintegrate_SOURCES) \ $(render_hkl_SOURCES) @@ -299,6 +306,10 @@ compare_hkl_SOURCES = compare_hkl.c sfac.c cell.c utils.c reflections.c \ statistics.c symmetry.c compare_hkl_LDADD = @LIBS@ +check_hkl_SOURCES = check_hkl.c sfac.c cell.c utils.c reflections.c \ + statistics.c symmetry.c + +check_hkl_LDADD = @LIBS@ powder_plot_SOURCES = powder_plot.c cell.c utils.c image.c hdf5-file.c \ detector.c @@ -410,6 +421,9 @@ clean-binPROGRAMS: calibrate_detector$(EXEEXT): $(calibrate_detector_OBJECTS) $(calibrate_detector_DEPENDENCIES) @rm -f calibrate_detector$(EXEEXT) $(AM_V_CCLD)$(LINK) $(calibrate_detector_OBJECTS) $(calibrate_detector_LDADD) $(LIBS) +check_hkl$(EXEEXT): $(check_hkl_OBJECTS) $(check_hkl_DEPENDENCIES) + @rm -f check_hkl$(EXEEXT) + $(AM_V_CCLD)$(LINK) $(check_hkl_OBJECTS) $(check_hkl_LDADD) $(LIBS) compare_hkl$(EXEEXT): $(compare_hkl_OBJECTS) $(compare_hkl_DEPENDENCIES) @rm -f compare_hkl$(EXEEXT) $(AM_V_CCLD)$(LINK) $(compare_hkl_OBJECTS) $(compare_hkl_LDADD) $(LIBS) @@ -456,6 +470,7 @@ distclean-compile: @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/beam-parameters.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/calibrate_detector.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/cell.Po@am__quote@ +@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/check_hkl.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/cl-utils.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/compare_hkl.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/cubeit.Po@am__quote@ diff --git a/src/check_hkl.c b/src/check_hkl.c new file mode 100644 index 00000000..449dfa87 --- /dev/null +++ b/src/check_hkl.c @@ -0,0 +1,367 @@ +/* + * check_hkl.c + * + * Characterise reflection lists + * + * (c) 2006-2010 Thomas White <taw@physics.org> + * + * 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 "sfac.h" +#include "reflections.h" +#include "statistics.h" +#include "symmetry.h" + + +/* Number of bins for plot of resolution shells */ +#define NBINS (10) + + +static void show_help(const char *s) +{ + printf("Syntax: %s [options] <file.hkl>\n\n", s); + printf( +"Characterise an intensity list.\n" +"\n" +" -h, --help Display this help message.\n" +" -y, --symmetry=<sym> The symmetry of both the input files.\n" +" -p, --pdb=<filename> PDB file to use (default: molecule.pdb).\n" +" --rmin=<res> Fix lower resolution limit for --shells (m^-1).\n" +" --rmax=<res> Fix upper resolution limit for --shells (m^-1).\n" +"\n"); +} + + +static void plot_shells(const double *ref1, ReflItemList *items, UnitCell *cell, + const char *sym, unsigned int *counts, + const double *sigma, double rmin_fix, double rmax_fix) +{ + double num[NBINS]; + int cts[NBINS]; + int possible[NBINS]; + unsigned int *counted; + unsigned int measurements[NBINS]; + unsigned int measured[NBINS]; + double total_vol, vol_per_shell; + double rmins[NBINS]; + double rmaxs[NBINS]; + double snr[NBINS]; + double den; + double rmin, rmax; + signed int h, k, l; + int i; + int ctot; + FILE *fh; + double snr_total = 0; + int nmeas = 0; + int nmeastot = 0; + int nout = 0; + + if ( cell == NULL ) { + ERROR("Need the unit cell to plot resolution shells.\n"); + return; + } + + fh = fopen("shells.dat", "w"); + if ( fh == NULL ) { + ERROR("Couldn't open 'shells.dat'\n"); + return; + } + + for ( i=0; i<NBINS; i++ ) { + num[i] = 0.0; + cts[i] = 0; + possible[i] = 0; + measured[i] = 0; + measurements[i] = 0; + snr[i] = 0; + } + + rmin = +INFINITY; + rmax = 0.0; + for ( i=0; i<num_items(items); i++ ) { + + struct refl_item *it; + signed int h, k, l; + double d; + + it = get_item(items, i); + h = it->h; k = it->k; l = it->l; + + d = resolution(cell, h, k, l) * 2.0; + if ( d > rmax ) rmax = d; + if ( d < rmin ) rmin = d; + + } + + STATUS("1/d goes from %f to %f nm^-1\n", rmin/1e9, rmax/1e9); + + /* Widen the range just a little bit */ + rmin -= 0.001e9; + rmax += 0.001e9; + + /* Fixed resolution shells if needed */ + if ( rmin_fix > 0.0 ) rmin = rmin_fix; + if ( rmax_fix > 0.0 ) rmax = rmax_fix; + + total_vol = pow(rmax, 3.0) - pow(rmin, 3.0); + vol_per_shell = total_vol / NBINS; + rmins[0] = rmin; + for ( i=1; i<NBINS; i++ ) { + + double r; + + r = vol_per_shell + pow(rmins[i-1], 3.0); + r = pow(r, 1.0/3.0); + + /* Shells of constant volume */ + rmaxs[i-1] = r; + rmins[i] = r; + + /* Shells of constant thickness */ + //rmins[i] = rmins[i-1] + (rmax-rmin)/NBINS; + //rmaxs[i-1] = rmins[i-1] + (rmax-rmin)/NBINS; + + STATUS("Shell %i: %f to %f\n", i-1, + rmins[i-1]/1e9, rmaxs[i-1]/1e9); + + } + rmaxs[NBINS-1] = rmax; + STATUS("Shell %i: %f to %f\n", NBINS-1, + rmins[NBINS-1]/1e9, rmaxs[NBINS-1]/1e9); + + /* Count the number of reflections possible in each shell */ + counted = new_list_count(); + for ( h=-50; h<=+50; h++ ) { + for ( k=-50; k<=+50; k++ ) { + for ( l=-50; l<=+50; l++ ) { + + double d; + signed int hs, ks, ls; + int bin; + + get_asymm(h, k, l, &hs, &ks, &ls, sym); + if ( lookup_count(counted, hs, ks, ls) ) continue; + set_count(counted, hs, ks, ls, 1); + + d = resolution(cell, h, k, l) * 2.0; + + bin = -1; + for ( i=0; i<NBINS; i++ ) { + if ( (d>rmins[i]) && (d<=rmaxs[i]) ) { + bin = i; + break; + } + } + if ( bin == -1 ) continue; + + possible[bin]++; + + } + } + } + free(counted); + + /* Characterise the first data set (only) */ + for ( i=0; i<num_items(items); i++ ) { + + struct refl_item *it; + signed int h, k, l; + double d; + int bin; + int j; + + it = get_item(items, i); + h = it->h; k = it->k; l = it->l; + + d = resolution(cell, h, k, l) * 2.0; + + bin = -1; + for ( j=0; j<NBINS; j++ ) { + if ( (d>rmins[j]) && (d<=rmaxs[j]) ) { + bin = j; + break; + } + } + if ( bin == -1 ) { + nout++; + continue; + } + + measured[bin]++; + measurements[bin] += lookup_count(counts, h, k, l); + snr[bin] += (lookup_intensity(ref1, h, k, l) / + lookup_intensity(sigma, h, k, l)); + snr_total += (lookup_intensity(ref1, h, k, l) / + lookup_intensity(sigma, h, k, l)); + nmeas++; + nmeastot += lookup_count(counts, h, k, l); + + } + STATUS("overall <snr> = %f\n", snr_total/(double)nmeas); + STATUS("%i measurements in total.\n", nmeastot); + STATUS("%i reflections in total.\n", nmeas); + + if ( nout ) { + STATUS("Warning; %i reflections outside resolution range.\n", + nout); + } + + for ( i=0; i<NBINS; i++ ) { + + double r, cen; + cen = rmins[i] + (rmaxs[i] - rmins[i])/2.0; + r = (num[i]/den)*((double)ctot/cts[i]); + fprintf(fh, "%f %i %i %5.2f %i %f %f\n", cen*1.0e-9, measured[i], + possible[i], 100.0*(float)measured[i]/possible[i], + measurements[i], (float)measurements[i]/measured[i], + (snr[i]/(double)measured[i])); + + } + + fclose(fh); +} + + +int main(int argc, char *argv[]) +{ + int c; + double *ref; + UnitCell *cell; + char *file = NULL; + char *sym = NULL; + int i; + ReflItemList *items; + ReflItemList *good_items; + char *pdb = NULL; + double *esd; + int rej = 0; + unsigned int *cts; + float rmin_fix = -1.0; + float rmax_fix = -1.0; + + /* Long options */ + const struct option longopts[] = { + {"help", 0, NULL, 'h'}, + {"symmetry", 1, NULL, 'y'}, + {"pdb", 1, NULL, 'p'}, + {"rmin", 1, NULL, 2}, + {"rmax", 1, NULL, 3}, + {0, 0, NULL, 0} + }; + + /* Short options */ + while ((c = getopt_long(argc, argv, "hy:p:", longopts, NULL)) != -1) { + + switch (c) { + case 'h' : + show_help(argv[0]); + return 0; + + case 'y' : + sym = strdup(optarg); + break; + + case 'p' : + pdb = strdup(optarg); + break; + + case 0 : + break; + + case 2 : + if ( sscanf(optarg, "%e", &rmin_fix) != 1 ) { + ERROR("Invalid value for --rmin\n"); + return 1; + } + break; + + case 3 : + if ( sscanf(optarg, "%e", &rmax_fix) != 1 ) { + ERROR("Invalid value for --rmax\n"); + return 1; + } + break; + + default : + return 1; + } + + } + + if ( argc != (optind+1) ) { + ERROR("Please provide exactly one HKL files to check.\n"); + return 1; + } + + if ( sym == NULL ) { + sym = strdup("1"); + } + + file = strdup(argv[optind++]); + + if ( pdb == NULL ) { + pdb = strdup("molecule.pdb"); + } + cell = load_cell_from_pdb(pdb); + free(pdb); + + ref = new_list_intensity(); + esd = new_list_sigma(); + cts = new_list_count(); + items = read_reflections(file, ref, NULL, cts, esd); + if ( items == NULL ) { + ERROR("Couldn't open file '%s'\n", file); + return 1; + } + + /* Reject reflections */ + good_items = new_items(); + for ( i=0; i<num_items(items); i++ ) { + + struct refl_item *it; + signed int h, k, l; + double val, sig; + int ig = 0; + double d; + + it = get_item(items, i); + h = it->h; k = it->k; l = it->l; + + val = lookup_intensity(ref, h, k, l); + sig = lookup_sigma(esd, h, k, l); + + if ( val < 3.0 * sig ) { + rej++; + ig = 1; + } + + d = 0.5/resolution(cell, h, k, l); + if ( d > 55.0e-10 ) ig = 1; + //if ( d < 15.0e-10 ) ig = 1; + + //if ( ig ) continue; + + add_item(good_items, h, k, l); + + } + + plot_shells(ref, items, cell, sym, cts, esd, rmin_fix, rmax_fix); + + return 0; +} diff --git a/src/compare_hkl.c b/src/compare_hkl.c index a0481559..0b3d5914 100644 --- a/src/compare_hkl.c +++ b/src/compare_hkl.c @@ -51,9 +51,7 @@ static void show_help(const char *s) static void plot_shells(const double *ref1, const double *ref2, ReflItemList *items, double scale, UnitCell *cell, - const char *sym, ReflItemList *characterise, - unsigned int *char_counts, const double *sigma, - double rmin_fix, double rmax_fix) + const char *sym, double rmin_fix, double rmax_fix) { double num[NBINS]; int cts[NBINS]; @@ -71,9 +69,6 @@ static void plot_shells(const double *ref1, const double *ref2, int i; int ctot; FILE *fh; - double snr_total = 0; - int nmeas = 0; - int nmeastot = 0; int nout = 0; if ( cell == NULL ) { @@ -181,53 +176,9 @@ static void plot_shells(const double *ref1, const double *ref2, } free(counted); - /* Characterise the first data set (only) */ - for ( i=0; i<num_items(characterise); i++ ) { - - struct refl_item *it; - signed int h, k, l; - double d; - int bin; - int j; - - it = get_item(characterise, i); - h = it->h; k = it->k; l = it->l; - - d = resolution(cell, h, k, l) * 2.0; - - bin = -1; - for ( j=0; j<NBINS; j++ ) { - if ( (d>rmins[j]) && (d<=rmaxs[j]) ) { - bin = j; - break; - } - } - if ( bin == -1 ) { - nout++; - continue; - } - - measured[bin]++; - measurements[bin] += lookup_count(char_counts, h, k, l); - snr[bin] += (lookup_intensity(ref1, h, k, l) / - lookup_intensity(sigma, h, k, l)); - snr_total += (lookup_intensity(ref1, h, k, l) / - lookup_intensity(sigma, h, k, l)); - nmeas++; - nmeastot += lookup_count(char_counts, h, k, l); - - } - STATUS("overall <snr> = %f\n", snr_total/(double)nmeas); - STATUS("%i measurements in total.\n", nmeastot); - STATUS("%i reflections in total.\n", nmeas); - - if ( nout ) { - STATUS("Warning; %i reflections outside resolution range.\n", - nout); - } - den = 0.0; ctot = 0; + nout = 0; for ( i=0; i<num_items(items); i++ ) { struct refl_item *it; @@ -251,7 +202,10 @@ static void plot_shells(const double *ref1, const double *ref2, } /* Outside resolution range? */ - if ( bin == -1 ) continue; + if ( bin == -1 ) { + nout++; + continue; + } i1 = lookup_intensity(ref1, h, k, l); i2 = lookup_intensity(ref2, h, k, l); @@ -264,6 +218,11 @@ static void plot_shells(const double *ref1, const double *ref2, } + if ( nout ) { + STATUS("Warning; %i reflections outside resolution range.\n", + nout); + } + for ( i=0; i<NBINS; i++ ) { double r, cen; @@ -496,7 +455,7 @@ int main(int argc, char *argv[]) if ( config_shells ) { plot_shells(ref1, ref2_transformed, icommon, scale_r1fi, - cell, sym, i1, cts1, esd1, rmin_fix, rmax_fix); + cell, sym, rmin_fix, rmax_fix); } if ( outfile != NULL ) { |