aboutsummaryrefslogtreecommitdiff
path: root/src/check_hkl.c
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2014-03-30 20:54:23 +0200
committerThomas White <taw@physics.org>2014-03-30 20:54:23 +0200
commitfb6a2d1f8cd57b8e40fd7b33dfebe07057cfc1e9 (patch)
tree4f4d88aadda6e856aca622d208996d859ca3dd4e /src/check_hkl.c
parent42113f1a6f1283e230182b1a6a0d2cb44a6ef3aa (diff)
check_hkl: Add Wilson plot
Diffstat (limited to 'src/check_hkl.c')
-rw-r--r--src/check_hkl.c103
1 files changed, 96 insertions, 7 deletions
diff --git a/src/check_hkl.c b/src/check_hkl.c
index c2ab0788..db41fac5 100644
--- a/src/check_hkl.c
+++ b/src/check_hkl.c
@@ -3,11 +3,11 @@
*
* Characterise reflection lists
*
- * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
*
* Authors:
- * 2010-2013 Thomas White <taw@physics.org>
+ * 2010-2014 Thomas White <taw@physics.org>
*
* This file is part of CrystFEL.
*
@@ -55,15 +55,89 @@ static void show_help(const char *s)
" -h, --help Display this help message.\n"
" -y, --symmetry=<sym> The symmetry of the input file.\n"
" -p, --pdb=<filename> PDB file to use.\n"
-" --rmin=<res> Fix lower resolution limit for resolution shells. (m^-1).\n"
-" --rmax=<res> Fix upper resolution limit for resolution shells. (m^-1).\n"
+" --rmin=<res> Lower resolution limit (1/d in m^-1).\n"
+" --rmax=<res> Upper resolution limit (1/d in m^-1).\n"
" --sigma-cutoff=<n> Discard reflections with I/sigma(I) < n.\n"
-" --nshells=<n> Use <n> resolution shells.\n"
-" --shell-file=<file> Write resolution shells to <file>.\n"
+" --nshells=<n> Use <n> resolution shells or bins.\n"
+" --wilson Calculate a Wilson plot\n"
+" --shell-file=<file> Write results table to <file>.\n"
"\n");
}
+static void wilson_plot(RefList *list, UnitCell *cell, const SymOpList *sym,
+ double rmin_fix, double rmax_fix, int nbins,
+ const char *filename)
+{
+ double rmin, rmax;
+ double s2min, s2max, s2step;
+ Reflection *refl;
+ RefListIterator *iter;
+ FILE *fh;
+ double *plot_i;
+ int *plot_n;
+ int i;
+
+ resolution_limits(list, cell, &rmin, &rmax);
+ 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;
+
+ s2min = pow(rmin/2.0, 2.0);
+ s2max = pow(rmax/2.0, 2.0);
+ s2step = (s2max - s2min)/nbins;
+
+ plot_i = malloc(nbins*sizeof(double));
+ if ( plot_i == NULL ) return;
+ for ( i=0; i<nbins; i++ ) plot_i[i] = 0.0;
+
+ plot_n = calloc(nbins, sizeof(int));
+ if ( plot_n == NULL ) return;
+
+ for ( refl = first_refl(list, &iter);
+ refl != NULL;
+ refl = next_refl(refl, iter) )
+ {
+ signed int h, k, l;
+ double d, intensity;
+ int bin;
+
+ get_indices(refl, &h, &k, &l);
+
+ d = resolution(cell, h, k, l); /* This gives "s" */
+ intensity = get_intensity(refl);
+
+ bin = (pow(d, 2.0) - s2min)/s2step;
+
+ /* FIXME: Divide by epsilon and Sigma */
+ plot_i[bin] += intensity;
+ plot_n[bin]++;
+
+ }
+
+ for ( i=0; i<nbins; i++ ) plot_i[i] = log(plot_i[i] / plot_n[i]);
+
+ fh = fopen(filename, "w");
+ if ( fh == NULL ) return;
+ fprintf(fh, " n s^2/A^-2 d/A <I>\n");
+ for ( i=0; i<nbins; i++ ) {
+ double s2 = s2min + (i+0.5)*s2step;
+ fprintf(fh, "%3i %8.6f %8.4f %10f\n", i,
+ s2/1e20, 0.5e10/sqrt(s2), plot_i[i]);
+ }
+ fclose(fh);
+
+ free(plot_i);
+ free(plot_n);
+}
+
+
static void plot_shells(RefList *list, UnitCell *cell, const SymOpList *sym,
double rmin_fix, double rmax_fix, int nshells,
const char *shell_file)
@@ -373,18 +447,25 @@ int main(int argc, char *argv[])
float rmax_fix = -1.0;
float sigma_cutoff = -INFINITY;
int nshells = 10;
+ int have_nshells = 0;
char *shell_file = NULL;
+ int wilson = 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},
{"sigma-cutoff", 1, NULL, 4},
{"nshells", 1, NULL, 5},
{"shell-file", 1, NULL, 6},
+
+ {"wilson", 0, &wilson, 1},
+
{0, 0, NULL, 0}
};
@@ -434,6 +515,7 @@ int main(int argc, char *argv[])
ERROR("Invalid value for --nshells\n");
return 1;
}
+ have_nshells = 1;
break;
case 6 :
@@ -518,7 +600,14 @@ int main(int argc, char *argv[])
rej, num_reflections(raw_list), sigma_cutoff);
reflist_free(raw_list);
- plot_shells(list, cell, sym, rmin_fix, rmax_fix, nshells, shell_file);
+ if ( wilson ) {
+ if ( !have_nshells ) nshells = 50;
+ wilson_plot(list, cell, sym, rmin_fix, rmax_fix, nshells,
+ shell_file);
+ } else {
+ plot_shells(list, cell, sym, rmin_fix, rmax_fix, nshells,
+ shell_file);
+ }
free_symoplist(sym);
reflist_free(list);