check_hkl.c 9.05 KB
Newer Older
1
2
3
4
5
/*
 * check_hkl.c
 *
 * Characterise reflection lists
 *
Thomas White's avatar
Thomas White committed
6
 * (c) 2006-2011 Thomas White <taw@physics.org>
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
 *
 * 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 "symmetry.h"
Thomas White's avatar
Thomas White committed
27
28
#include "reflist.h"
#include "reflist-utils.h"
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45


/* 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"
Thomas White's avatar
Thomas White committed
46
"      --sigma-cutoff=<n>     Discard reflections with I/sigma(I) < n.\n"
47
48
49
50
"\n");
}


Thomas White's avatar
Thomas White committed
51
static void plot_shells(RefList *list, UnitCell *cell, const SymOpList *sym,
Thomas White's avatar
Thomas White committed
52
                        double rmin_fix, double rmax_fix)
53
54
55
56
57
58
59
60
61
62
63
{
	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];
Thomas White's avatar
Thomas White committed
64
65
	double mean[NBINS];
	double var[NBINS];
66
67
68
69
70
71
72
73
	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
74
75
	Reflection *refl;
	RefListIterator *iter;
76
77
78
79
	int hmax, kmax, lmax;
	double asx, asy, asz;
	double bsx, bsy, bsz;
	double csx, csy, csz;
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98

	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
99
100
		var[i] = 0;
		mean[i] = 0;
101
102
	}

Thomas White's avatar
Thomas White committed
103
104
105
106
107
108
	/* Iterate over all common reflections and calculate min and max
	 * resolution */
	rmin = +INFINITY;  rmax = 0.0;
	for ( refl = first_refl(list, &iter);
	      refl != NULL;
	      refl = next_refl(refl, iter) ) {
109
110
111
112

		signed int h, k, l;
		double d;

Thomas White's avatar
Thomas White committed
113
		get_indices(refl, &h, &k, &l);
114
115
116
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
154
155
156
157
158

		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();
159
160
161
162
163
164
165
166
167
	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++ ) {
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183

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

		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;

Thomas White's avatar
Thomas White committed
184
		get_asymm(sym, h, k, l, &hs, &ks, &ls);
Thomas White's avatar
Thomas White committed
185
186
187
		if ( lookup_count(counted, hs, ks, ls) ) continue;
		set_count(counted, hs, ks, ls, 1);

188
189
190
191
192
193
194
		possible[bin]++;

	}
	}
	}
	free(counted);

195
	/* Calculate means */
Thomas White's avatar
Thomas White committed
196
197
198
	for ( refl = first_refl(list, &iter);
	      refl != NULL;
	      refl = next_refl(refl, iter) ) {
199
200
201
202
203
204

		signed int h, k, l;
		double d;
		int bin;
		int j;

Thomas White's avatar
Thomas White committed
205
		get_indices(refl, &h, &k, &l);
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221

		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]++;
Thomas White's avatar
Thomas White committed
222
		mean[bin] += get_intensity(refl);
223
224
225
226
227
228
229
230

	}

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

	/* Characterise the data set */
Thomas White's avatar
Thomas White committed
231
232
233
	for ( refl = first_refl(list, &iter);
	      refl != NULL;
	      refl = next_refl(refl, iter) ) {
234
235
236
237
238
239
240

		signed int h, k, l;
		double d;
		int bin;
		int j;
		double val, esd;

Thomas White's avatar
Thomas White committed
241
		get_indices(refl, &h, &k, &l);
242
243

		d = resolution(cell, h, k, l) * 2.0;
Thomas White's avatar
Thomas White committed
244
245
		val = get_intensity(refl);
		esd = get_esd_intensity(refl);
246
247
248
249
250
251
252
253
254
255
256
257
258

		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
259
260
		if ( !isfinite(val/esd) ) continue;

261
		/* measured[bin] was done earlier */
Thomas White's avatar
Thomas White committed
262
		measurements[bin] += get_redundancy(refl);
263
264
		snr[bin] += val / esd;
		snr_total += val / esd;
265
		nmeas++;
Thomas White's avatar
Thomas White committed
266
		nmeastot += get_redundancy(refl);
267

268
269
		var[bin] += pow(val-mean[bin], 2.0);

270
271
272
273
274
275
276
277
278
279
	}
	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);
	}

280
281
	fprintf(fh, "1/d centre   # refs Possible  Compl       "
		    "Meas   Red   SNR    Std dev       Mean\n");
282
283
	for ( i=0; i<NBINS; i++ ) {

284
		double cen;
285
		cen = rmins[i] + (rmaxs[i] - rmins[i])/2.0;
Thomas White's avatar
Thomas White committed
286
287
		fprintf(fh, "%10.3f %8i %8i %6.2f %10i %5.1f"
		            " %5.2f %10.2f %10.2f\n",
288
289
290
291
292
293
294
295
296
		        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]);
297
298
299
300
301
302
303
304
305
306
307
308

	}

	fclose(fh);
}


int main(int argc, char *argv[])
{
	int c;
	UnitCell *cell;
	char *file = NULL;
Thomas White's avatar
Thomas White committed
309
310
	char *sym_str = NULL;
	SymOpList *sym;
Thomas White's avatar
Thomas White committed
311
312
313
314
	RefList *raw_list;
	RefList *list;
	Reflection *refl;
	RefListIterator *iter;
315
316
317
318
	char *pdb = NULL;
	int rej = 0;
	float rmin_fix = -1.0;
	float rmax_fix = -1.0;
Thomas White's avatar
Thomas White committed
319
	float sigma_cutoff = -INFINITY;
320
321
322
323
324
325

	/* 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
326
327
328
		{"rmin",               1, NULL,                2},
		{"rmax",               1, NULL,                3},
		{"sigma-cutoff",       1, NULL,                4},
329
330
331
332
333
334
335
336
337
338
339
340
		{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
341
			sym_str = strdup(optarg);
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
			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
365
366
367
368
369
370
371
		case 4 :
			if ( sscanf(optarg, "%f", &sigma_cutoff) != 1 ) {
				ERROR("Invalid value for --sigma-cutoff\n");
				return 1;
			}
			break;

372
373
374
375
376
377
378
		default :
			return 1;
		}

	}

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

Thomas White's avatar
Thomas White committed
383
384
	if ( sym_str == NULL ) {
		sym_str = strdup("1");
385
	}
Thomas White's avatar
Thomas White committed
386
387
	sym = get_pointgroup(sym_str);
	free(sym_str);
388
389
390
391
392
393
394
395
396

	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
397
398
399
	raw_list = read_reflections(file);
	if ( raw_list == NULL ) {
		ERROR("Couldn't read file '%s'\n", file);
400
401
402
		return 1;
	}

Thomas White's avatar
Thomas White committed
403
	/* Check that the intensities have the correct symmetry */
Thomas White's avatar
Thomas White committed
404
	if ( check_list_symmetry(raw_list, sym) ) {
Thomas White's avatar
Thomas White committed
405
		ERROR("The input reflection list does not appear to"
Thomas White's avatar
Thomas White committed
406
		      " have symmetry %s\n", symmetry_name(sym));
Thomas White's avatar
Thomas White committed
407
408
409
		return 1;
	}

Thomas White's avatar
Thomas White committed
410
411
	/* Reject some reflections */
	list = reflist_new();
Thomas White's avatar
Thomas White committed
412
	for ( refl = first_refl(raw_list, &iter);
Thomas White's avatar
Thomas White committed
413
414
	      refl != NULL;
	      refl = next_refl(refl, iter) ) {
415
416
417
418

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

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

Thomas White's avatar
Thomas White committed
423
424
		val = get_intensity(refl);
		sig = get_esd_intensity(refl);
425

Thomas White's avatar
Thomas White committed
426
		if ( val < sigma_cutoff * sig ) {
427
428
429
430
			rej++;
			ig = 1;
		}

Thomas White's avatar
Thomas White committed
431
		if ( ig ) continue;
432

Thomas White's avatar
Thomas White committed
433
434
		new = add_refl(list, h, k, l);
		copy_data(new, refl);
435
436

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

Thomas White's avatar
Thomas White committed
440
	plot_shells(list, cell, sym, rmin_fix, rmax_fix);
441
442
443

	return 0;
}