check_hkl.c 9.12 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

	}

	fclose(fh);
Thomas White's avatar
Thomas White committed
301
302

	STATUS("Resolution shell information written to shells.dat.\n");
303
304
305
306
307
308
309
310
}


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

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

374
375
376
377
378
379
380
		default :
			return 1;
		}

	}

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

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

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

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

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

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

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

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

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

Thomas White's avatar
Thomas White committed
433
		if ( ig ) continue;
434

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

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

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

	return 0;
}