check_hkl.c 9.37 KB
Newer Older
1
2
3
4
5
/*
 * check_hkl.c
 *
 * Characterise reflection lists
 *
Thomas White's avatar
Thomas White committed
6
 * Copyright © 2012 Thomas White <taw@physics.org>
7
 *
Thomas White's avatar
Thomas White committed
8
9
10
11
12
13
14
15
16
17
18
19
20
21
 * This file is part of CrystFEL.
 *
 * CrystFEL is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * CrystFEL is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with CrystFEL.  If not, see <http://www.gnu.org/licenses/>.
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
 *
 */


#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 "symmetry.h"
Thomas White's avatar
Thomas White committed
40
41
#include "reflist.h"
#include "reflist-utils.h"
42
43
44
45
46
47
48
49
50
51
52
53
54


/* 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"
Thomas White's avatar
Thomas White committed
55
"  -y, --symmetry=<sym>       The symmetry of the input file.\n"
56
57
58
"  -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"
Thomas White's avatar
Thomas White committed
59
"      --sigma-cutoff=<n>     Discard reflections with I/sigma(I) < n.\n"
60
61
62
63
"\n");
}


Thomas White's avatar
Thomas White committed
64
static void plot_shells(RefList *list, UnitCell *cell, const SymOpList *sym,
Thomas White's avatar
Thomas White committed
65
                        double rmin_fix, double rmax_fix)
66
67
68
69
70
71
72
73
74
75
{
	double num[NBINS];
	int cts[NBINS];
	int possible[NBINS];
	unsigned int measurements[NBINS];
	unsigned int measured[NBINS];
	double total_vol, vol_per_shell;
	double rmins[NBINS];
	double rmaxs[NBINS];
	double snr[NBINS];
Thomas White's avatar
Thomas White committed
76
77
	double mean[NBINS];
	double var[NBINS];
78
79
80
81
82
83
84
85
	double rmin, rmax;
	signed int h, k, l;
	int i;
	FILE *fh;
	double snr_total = 0;
	int nmeas = 0;
	int nmeastot = 0;
	int nout = 0;
Thomas White's avatar
Thomas White committed
86
87
	Reflection *refl;
	RefListIterator *iter;
88
	RefList *counted;
89
90
91
92
	int hmax, kmax, lmax;
	double asx, asy, asz;
	double bsx, bsy, bsz;
	double csx, csy, csz;
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111

	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;
Thomas White's avatar
Thomas White committed
112
113
		var[i] = 0;
		mean[i] = 0;
114
115
	}

116
	resolution_limits(list, cell, &rmin, &rmax);
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
	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 */
154
	counted = reflist_new();
155
156
157
158
159
160
161
162
163
	cell_get_reciprocal(cell, &asx, &asy, &asz,
	                          &bsx, &bsy, &bsz,
	                          &csx, &csy, &csz);
	hmax = rmax / modulus(asx, asy, asz);
	kmax = rmax / modulus(bsx, bsy, bsz);
	lmax = rmax / modulus(csx, csy, csz);
	for ( h=-hmax; h<hmax; h++ ) {
	for ( k=-kmax; k<kmax; k++ ) {
	for ( l=-lmax; l<lmax; l++ ) {
164
165
166
167
168

		double d;
		signed int hs, ks, ls;
		int bin;

169
		d = 2.0 * resolution(cell, h, k, l);
170
171
172
173
174
175
176
177
178
179

		bin = -1;
		for ( i=0; i<NBINS; i++ ) {
			if ( (d>rmins[i]) && (d<=rmaxs[i]) ) {
				bin = i;
				break;
			}
		}
		if ( bin == -1 ) continue;

Thomas White's avatar
Thomas White committed
180
		get_asymm(sym, h, k, l, &hs, &ks, &ls);
181
182
		if ( find_refl(counted, hs, ks, ls) != NULL ) continue;
		add_refl(counted, hs, ks, ls);
Thomas White's avatar
Thomas White committed
183

184
185
186
187
188
		possible[bin]++;

	}
	}
	}
189
	reflist_free(counted);
190

191
	/* Calculate means */
Thomas White's avatar
Thomas White committed
192
193
	for ( refl = first_refl(list, &iter);
	      refl != NULL;
194
195
	      refl = next_refl(refl, iter) )
	{
196
197
198
199
200
		signed int h, k, l;
		double d;
		int bin;
		int j;

Thomas White's avatar
Thomas White committed
201
		get_indices(refl, &h, &k, &l);
202
203
204
205
206
207
208
209
210
211

		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;
			}
		}
212
		if ( bin == -1 ) continue;
213
214

		measured[bin]++;
Thomas White's avatar
Thomas White committed
215
		mean[bin] += get_intensity(refl);
216
217
218
219
220
221
222
	}

	for ( i=0; i<NBINS; i++ ) {
		mean[i] /= (double)measured[i];
	}

	/* Characterise the data set */
Thomas White's avatar
Thomas White committed
223
224
	for ( refl = first_refl(list, &iter);
	      refl != NULL;
225
226
	      refl = next_refl(refl, iter) )
	{
227
228
229
230
231
232
		signed int h, k, l;
		double d;
		int bin;
		int j;
		double val, esd;

Thomas White's avatar
Thomas White committed
233
		get_indices(refl, &h, &k, &l);
234
235

		d = resolution(cell, h, k, l) * 2.0;
Thomas White's avatar
Thomas White committed
236
237
		val = get_intensity(refl);
		esd = get_esd_intensity(refl);
238
239
240
241
242
243
244
245
246
247
248
249
250

		bin = -1;
		for ( j=0; j<NBINS; j++ ) {
			if ( (d>rmins[j]) && (d<=rmaxs[j]) ) {
				bin = j;
				break;
			}
		}
		if ( bin == -1 ) {
			nout++;
			continue;
		}

Thomas White's avatar
Thomas White committed
251
252
		if ( !isfinite(val/esd) ) continue;

253
		/* measured[bin] was done earlier */
Thomas White's avatar
Thomas White committed
254
		measurements[bin] += get_redundancy(refl);
255
256
		snr[bin] += val / esd;
		snr_total += val / esd;
257
		nmeas++;
Thomas White's avatar
Thomas White committed
258
		nmeastot += get_redundancy(refl);
259

260
		var[bin] += pow(val-mean[bin], 2.0);
261
262
263
264
265
266
267
268
269
270
	}
	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);
	}

271
272
	fprintf(fh, "1/d centre   # refs Possible  Compl       "
		    "Meas   Red   SNR    Std dev       Mean\n");
273
274
	for ( i=0; i<NBINS; i++ ) {

275
		double cen;
276
		cen = rmins[i] + (rmaxs[i] - rmins[i])/2.0;
Thomas White's avatar
Thomas White committed
277
278
		fprintf(fh, "%10.3f %8i %8i %6.2f %10i %5.1f"
		            " %5.2f %10.2f %10.2f\n",
279
280
281
282
283
284
285
286
287
		        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],
		        sqrt(var[i]/measured[i]),
		        mean[i]);
288
289
290
291

	}

	fclose(fh);
Thomas White's avatar
Thomas White committed
292
293

	STATUS("Resolution shell information written to shells.dat.\n");
294
295
296
297
298
299
300
301
}


int main(int argc, char *argv[])
{
	int c;
	UnitCell *cell;
	char *file = NULL;
Thomas White's avatar
Thomas White committed
302
303
	char *sym_str = NULL;
	SymOpList *sym;
Thomas White's avatar
Thomas White committed
304
305
306
307
	RefList *raw_list;
	RefList *list;
	Reflection *refl;
	RefListIterator *iter;
308
309
310
311
	char *pdb = NULL;
	int rej = 0;
	float rmin_fix = -1.0;
	float rmax_fix = -1.0;
Thomas White's avatar
Thomas White committed
312
	float sigma_cutoff = -INFINITY;
313
314
315
316
317
318

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
		{"symmetry",           1, NULL,               'y'},
		{"pdb",                1, NULL,               'p'},
Thomas White's avatar
Thomas White committed
319
320
321
		{"rmin",               1, NULL,                2},
		{"rmax",               1, NULL,                3},
		{"sigma-cutoff",       1, NULL,                4},
322
323
324
325
326
327
328
329
330
331
332
333
		{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' :
Thomas White's avatar
Thomas White committed
334
			sym_str = strdup(optarg);
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
			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;

Thomas White's avatar
Thomas White committed
358
359
360
361
362
363
364
		case 4 :
			if ( sscanf(optarg, "%f", &sigma_cutoff) != 1 ) {
				ERROR("Invalid value for --sigma-cutoff\n");
				return 1;
			}
			break;

365
366
367
368
369
370
371
		default :
			return 1;
		}

	}

	if ( argc != (optind+1) ) {
Thomas White's avatar
Thomas White committed
372
		ERROR("Please provide exactly one HKL file to check.\n");
373
374
375
		return 1;
	}

Thomas White's avatar
Thomas White committed
376
377
	if ( sym_str == NULL ) {
		sym_str = strdup("1");
378
	}
Thomas White's avatar
Thomas White committed
379
380
	sym = get_pointgroup(sym_str);
	free(sym_str);
381
382
383
384
385
386
387
388
389

	file = strdup(argv[optind++]);

	if ( pdb == NULL ) {
		pdb = strdup("molecule.pdb");
	}
	cell = load_cell_from_pdb(pdb);
	free(pdb);

Thomas White's avatar
Thomas White committed
390
391
392
	raw_list = read_reflections(file);
	if ( raw_list == NULL ) {
		ERROR("Couldn't read file '%s'\n", file);
393
394
395
		return 1;
	}

Thomas White's avatar
Thomas White committed
396
	/* Check that the intensities have the correct symmetry */
Thomas White's avatar
Thomas White committed
397
	if ( check_list_symmetry(raw_list, sym) ) {
Thomas White's avatar
Thomas White committed
398
		ERROR("The input reflection list does not appear to"
Thomas White's avatar
Thomas White committed
399
		      " have symmetry %s\n", symmetry_name(sym));
Thomas White's avatar
Thomas White committed
400
401
402
		return 1;
	}

Thomas White's avatar
Thomas White committed
403
404
	/* Reject some reflections */
	list = reflist_new();
Thomas White's avatar
Thomas White committed
405
	for ( refl = first_refl(raw_list, &iter);
Thomas White's avatar
Thomas White committed
406
407
	      refl != NULL;
	      refl = next_refl(refl, iter) ) {
408
409
410
411

		signed int h, k, l;
		double val, sig;
		int ig = 0;
Thomas White's avatar
Thomas White committed
412
		Reflection *new;
413

Thomas White's avatar
Thomas White committed
414
		get_indices(refl, &h, &k, &l);
415

Thomas White's avatar
Thomas White committed
416
417
		val = get_intensity(refl);
		sig = get_esd_intensity(refl);
418

Thomas White's avatar
Thomas White committed
419
		if ( val < sigma_cutoff * sig ) {
420
421
422
423
			rej++;
			ig = 1;
		}

Thomas White's avatar
Thomas White committed
424
		if ( ig ) continue;
425

Thomas White's avatar
Thomas White committed
426
427
		new = add_refl(list, h, k, l);
		copy_data(new, refl);
428
429

	}
Thomas White's avatar
Thomas White committed
430
431
	STATUS("Discarded %i reflections (out of %i) with I/sigma(I) < %f\n",
	       rej, num_reflections(raw_list), sigma_cutoff);
432

Thomas White's avatar
Thomas White committed
433
	plot_shells(list, cell, sym, rmin_fix, rmax_fix);
434
435
436

	return 0;
}