aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2013-03-07 10:38:18 +0100
committerThomas White <taw@physics.org>2013-03-07 10:38:18 +0100
commit82ac8b4ec60e47e09f33cf925d452810d590b0a5 (patch)
treedf0c362633a3689dabb34389175c0c62b40796c4
parent33260b9a7bf27a87afa1dde789f0b0d9ce28eaab (diff)
check_hkl: Add --nshells
-rw-r--r--src/check_hkl.c123
1 files changed, 92 insertions, 31 deletions
diff --git a/src/check_hkl.c b/src/check_hkl.c
index f82855ab..5aa9874d 100644
--- a/src/check_hkl.c
+++ b/src/check_hkl.c
@@ -3,11 +3,11 @@
*
* Characterise reflection lists
*
- * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY,
- * a research centre of the Helmholtz Association.
+ * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY,
+ * a research centre of the Helmholtz Association.
*
* Authors:
- * 2010-2012 Thomas White <taw@physics.org>
+ * 2010-2013 Thomas White <taw@physics.org>
*
* This file is part of CrystFEL.
*
@@ -46,10 +46,6 @@
#include "cell-utils.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);
@@ -62,23 +58,24 @@ static void show_help(const char *s)
" --rmin=<res> Fix lower resolution limit for resolution shells. (m^-1).\n"
" --rmax=<res> Fix upper resolution limit for resolution shells. (m^-1).\n"
" --sigma-cutoff=<n> Discard reflections with I/sigma(I) < n.\n"
+" --nshells=<n> Use <n> resolution shells.\n"
"\n");
}
static void plot_shells(RefList *list, UnitCell *cell, const SymOpList *sym,
- double rmin_fix, double rmax_fix)
+ double rmin_fix, double rmax_fix, int nshells)
{
- int possible[NBINS];
- unsigned int measurements[NBINS];
- unsigned int measured[NBINS];
- unsigned int snr_measured[NBINS];
+ int *possible;
+ unsigned int *measurements;
+ unsigned int *measured;
+ unsigned int *snr_measured;
double total_vol, vol_per_shell;
- double rmins[NBINS];
- double rmaxs[NBINS];
- double snr[NBINS];
- double mean[NBINS];
- double var[NBINS];
+ double *rmins;
+ double *rmaxs;
+ double *snr;
+ double *mean;
+ double *var;
double rmin, rmax;
signed int h, k, l;
int i;
@@ -96,13 +93,47 @@ static void plot_shells(RefList *list, UnitCell *cell, const SymOpList *sym,
double bsx, bsy, bsz;
double csx, csy, csz;
+ possible = malloc(nshells*sizeof(int));
+ measurements = malloc(nshells*sizeof(unsigned int));
+ measured = malloc(nshells*sizeof(unsigned int));
+ snr_measured = malloc(nshells*sizeof(unsigned int));
+ if ( (possible == NULL) || (measurements == NULL)
+ || (measured == NULL) || (snr_measured == NULL) ) {
+ ERROR("Couldn't allocate memory.\n");
+ free(possible);
+ free(measurements);
+ free(measured);
+ free(snr_measured);
+ return;
+ }
+
+ rmins = malloc(nshells*sizeof(double));
+ rmaxs = malloc(nshells*sizeof(double));
+ snr = malloc(nshells*sizeof(double));
+ mean = malloc(nshells*sizeof(double));
+ var = malloc(nshells*sizeof(double));
+ if ( (rmins == NULL) || (rmaxs == NULL) || (snr == NULL)
+ || (mean == NULL) || (var == NULL) ) {
+ ERROR("Couldn't allocate memory.\n");
+ free(possible);
+ free(measurements);
+ free(measured);
+ free(snr_measured);
+ free(rmins);
+ free(rmaxs);
+ free(snr);
+ free(mean);
+ free(var);
+ return;
+ }
+
fh = fopen("shells.dat", "w");
if ( fh == NULL ) {
ERROR("Couldn't open 'shells.dat'\n");
return;
}
- for ( i=0; i<NBINS; i++ ) {
+ for ( i=0; i<nshells; i++ ) {
possible[i] = 0;
measured[i] = 0;
snr_measured[i] = 0;
@@ -124,9 +155,9 @@ static void plot_shells(RefList *list, UnitCell *cell, const SymOpList *sym,
if ( rmax_fix > 0.0 ) rmax = rmax_fix;
total_vol = pow(rmax, 3.0) - pow(rmin, 3.0);
- vol_per_shell = total_vol / NBINS;
+ vol_per_shell = total_vol / nshells;
rmins[0] = rmin;
- for ( i=1; i<NBINS; i++ ) {
+ for ( i=1; i<nshells; i++ ) {
double r;
@@ -138,16 +169,16 @@ static void plot_shells(RefList *list, UnitCell *cell, const SymOpList *sym,
rmins[i] = r;
/* Shells of constant thickness */
- //rmins[i] = rmins[i-1] + (rmax-rmin)/NBINS;
- //rmaxs[i-1] = rmins[i-1] + (rmax-rmin)/NBINS;
+ //rmins[i] = rmins[i-1] + (rmax-rmin)/nshells;
+ //rmaxs[i-1] = rmins[i-1] + (rmax-rmin)/nshells;
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);
+ rmaxs[nshells-1] = rmax;
+ STATUS("Shell %i: %f to %f\n", nshells-1,
+ rmins[nshells-1]/1e9, rmaxs[nshells-1]/1e9);
/* Count the number of reflections possible in each shell */
counted = reflist_new();
@@ -168,9 +199,12 @@ static void plot_shells(RefList *list, UnitCell *cell, const SymOpList *sym,
d = 2.0 * resolution(cell, h, k, l);
if ( forbidden_reflection(cell, h, k, l) ) continue;
+ if ( h % 2 ) continue;
+ if ( k % 2 ) continue;
+ if ( l % 2 ) continue;
bin = -1;
- for ( i=0; i<NBINS; i++ ) {
+ for ( i=0; i<nshells; i++ ) {
if ( (d>rmins[i]) && (d<=rmaxs[i]) ) {
bin = i;
break;
@@ -200,13 +234,17 @@ static void plot_shells(RefList *list, UnitCell *cell, const SymOpList *sym,
int j;
get_indices(refl, &h, &k, &l);
+ if ( forbidden_reflection(cell, h, k, l) ) continue;
+ if ( h % 2 ) continue;
+ if ( k % 2 ) continue;
+ if ( l % 2 ) continue;
d = resolution(cell, h, k, l) * 2.0;
val = get_intensity(refl);
esd = get_esd_intensity(refl);
bin = -1;
- for ( j=0; j<NBINS; j++ ) {
+ for ( j=0; j<nshells; j++ ) {
if ( (d>rmins[j]) && (d<=rmaxs[j]) ) {
bin = j;
break;
@@ -221,7 +259,7 @@ static void plot_shells(RefList *list, UnitCell *cell, const SymOpList *sym,
}
- for ( i=0; i<NBINS; i++ ) {
+ for ( i=0; i<nshells; i++ ) {
mean[i] /= (double)measured[i];
}
@@ -237,13 +275,17 @@ static void plot_shells(RefList *list, UnitCell *cell, const SymOpList *sym,
double val, esd;
get_indices(refl, &h, &k, &l);
+ if ( forbidden_reflection(cell, h, k, l) ) continue;
+ if ( h % 2 ) continue;
+ if ( k % 2 ) continue;
+ if ( l % 2 ) continue;
d = resolution(cell, h, k, l) * 2.0;
val = get_intensity(refl);
esd = get_esd_intensity(refl);
bin = -1;
- for ( j=0; j<NBINS; j++ ) {
+ for ( j=0; j<nshells; j++ ) {
if ( (d>rmins[j]) && (d<=rmaxs[j]) ) {
bin = j;
break;
@@ -287,7 +329,7 @@ static void plot_shells(RefList *list, UnitCell *cell, const SymOpList *sym,
fprintf(fh, "1/d centre # refs Possible Compl "
"Meas Red SNR Std dev Mean d(A)\n");
- for ( i=0; i<NBINS; i++ ) {
+ for ( i=0; i<nshells; i++ ) {
double cen;
cen = rmins[i] + (rmaxs[i] - rmins[i])/2.0;
@@ -308,6 +350,16 @@ static void plot_shells(RefList *list, UnitCell *cell, const SymOpList *sym,
fclose(fh);
STATUS("Resolution shell information written to shells.dat.\n");
+
+ free(possible);
+ free(measurements);
+ free(measured);
+ free(snr_measured);
+ free(rmins);
+ free(rmaxs);
+ free(snr);
+ free(mean);
+ free(var);
}
@@ -327,6 +379,7 @@ int main(int argc, char *argv[])
float rmin_fix = -1.0;
float rmax_fix = -1.0;
float sigma_cutoff = -INFINITY;
+ int nshells = 10;
/* Long options */
const struct option longopts[] = {
@@ -336,6 +389,7 @@ int main(int argc, char *argv[])
{"rmin", 1, NULL, 2},
{"rmax", 1, NULL, 3},
{"sigma-cutoff", 1, NULL, 4},
+ {"nshells", 1, NULL, 5},
{0, 0, NULL, 0}
};
@@ -380,6 +434,13 @@ int main(int argc, char *argv[])
}
break;
+ case 5 :
+ if ( sscanf(optarg, "%i", &nshells) != 1 ) {
+ ERROR("Invalid value for --nshells\n");
+ return 1;
+ }
+ break;
+
case '?' :
break;
@@ -454,7 +515,7 @@ int main(int argc, char *argv[])
STATUS("Discarded %i reflections (out of %i) with I/sigma(I) < %f\n",
rej, num_reflections(raw_list), sigma_cutoff);
- plot_shells(list, cell, sym, rmin_fix, rmax_fix);
+ plot_shells(list, cell, sym, rmin_fix, rmax_fix, nshells);
return 0;
}