check_hkl.c 9.86 KB
Newer Older
1
2
3
4
5
/*
 * check_hkl.c
 *
 * Characterise reflection lists
 *
6
7
8
9
 * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY,
 *                  a research centre of the Helmholtz Association.
 *
 * Authors:
10
 *   2010-2012 Thomas White <taw@physics.org>
11
 *
Thomas White's avatar
Thomas White committed
12
13
14
15
16
17
18
19
20
21
22
23
24
25
 * 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/>.
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
 *
 */


#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
44
45
#include "reflist.h"
#include "reflist-utils.h"
46
#include "cell-utils.h"
47
48
49
50
51
52
53
54
55
56
57
58
59


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


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

	fh = fopen("shells.dat", "w");
	if ( fh == NULL ) {
		ERROR("Couldn't open 'shells.dat'\n");
		return;
	}

	for ( i=0; i<NBINS; i++ ) {
		possible[i] = 0;
		measured[i] = 0;
108
		snr_measured[i] = 0;
109
110
		measurements[i] = 0;
		snr[i] = 0;
Thomas White's avatar
Thomas White committed
111
112
		var[i] = 0;
		mean[i] = 0;
113
114
	}

115
	resolution_limits(list, cell, &rmin, &rmax);
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
	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 */
153
	counted = reflist_new();
154
155
156
157
158
159
160
161
162
	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++ ) {
163
164
165
166
167

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

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

		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
179
		get_asymm(sym, h, k, l, &hs, &ks, &ls);
180
181
		if ( find_refl(counted, hs, ks, ls) != NULL ) continue;
		add_refl(counted, hs, ks, ls);
Thomas White's avatar
Thomas White committed
182

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

	}
	}
	}
188
	reflist_free(counted);
189

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

Thomas White's avatar
Thomas White committed
200
		get_indices(refl, &h, &k, &l);
201
202

		d = resolution(cell, h, k, l) * 2.0;
203
204
		val = get_intensity(refl);
		esd = get_esd_intensity(refl);
205
206
207
208
209
210
211
212

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

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

		if ( !isfinite(val/esd) ) nsilly++;

220
221
222
223
224
225
226
	}

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

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

Thomas White's avatar
Thomas White committed
237
		get_indices(refl, &h, &k, &l);
238
239

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

		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] was done earlier */
Thomas White's avatar
Thomas White committed
256
		measurements[bin] += get_redundancy(refl);
257
258
259
260
261

		if ( isfinite(val/esd) ) {
			snr[bin] += val / esd;
			snr_total += val / esd;
			snr_measured[bin]++;
Thomas White's avatar
Thomas White committed
262
263
		} else {
			nsilly++;
264
265
		}

266
		nrefl++;
Thomas White's avatar
Thomas White committed
267
		nmeastot += get_redundancy(refl);
268

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

271
	}
272
	STATUS("overall <snr> = %f\n", snr_total/(double)nrefl);
273
	STATUS("%i measurements in total.\n", nmeastot);
274
	STATUS("%i reflections in total.\n", nrefl);
275
276
277
278
279
280

	if ( nout ) {
		STATUS("Warning; %i reflections outside resolution range.\n",
		       nout);
	}

281
282
283
284
285
	if ( nsilly ) {
		STATUS("Warning; %i reflections had infinite or invalid values"
		       " of I/sigma(I).\n", nsilly);
	}

286
	fprintf(fh, "1/d centre   # refs Possible  Compl       "
287
		    "Meas   Red   SNR    Std dev       Mean     d(A)\n");
288
289
	for ( i=0; i<NBINS; i++ ) {

290
		double cen;
291
		cen = rmins[i] + (rmaxs[i] - rmins[i])/2.0;
Thomas White's avatar
Thomas White committed
292
		fprintf(fh, "%10.3f %8i %8i %6.2f %10i %5.1f"
293
		            " %5.2f %10.2f %10.2f %8.2f\n",
294
295
296
		        cen*1.0e-9,
		        measured[i],
		        possible[i],
297
		        100.0*(double)measured[i]/possible[i],
298
		        measurements[i],
299
		        (double)measurements[i]/measured[i],
300
		        snr[i]/(double)snr_measured[i],
301
		        sqrt(var[i]/measured[i]),
302
		        mean[i], (1.0/cen)*1e10);
303
304
305
306

	}

	fclose(fh);
Thomas White's avatar
Thomas White committed
307
308

	STATUS("Resolution shell information written to shells.dat.\n");
309
310
311
312
313
314
315
316
}


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

	/* 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
334
335
336
		{"rmin",               1, NULL,                2},
		{"rmax",               1, NULL,                3},
		{"sigma-cutoff",       1, NULL,                4},
337
338
339
340
341
342
343
		{0, 0, NULL, 0}
	};

	/* Short options */
	while ((c = getopt_long(argc, argv, "hy:p:", longopts, NULL)) != -1) {

		switch (c) {
Thomas White's avatar
Thomas White committed
344
345

			case 'h' :
346
347
348
			show_help(argv[0]);
			return 0;

Thomas White's avatar
Thomas White committed
349
			case 'y' :
Thomas White's avatar
Thomas White committed
350
			sym_str = strdup(optarg);
351
352
			break;

Thomas White's avatar
Thomas White committed
353
			case 'p' :
354
355
356
			pdb = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
357
			case 0 :
358
359
			break;

Thomas White's avatar
Thomas White committed
360
			case 2 :
361
362
363
364
365
366
			if ( sscanf(optarg, "%e", &rmin_fix) != 1 ) {
				ERROR("Invalid value for --rmin\n");
				return 1;
			}
			break;

Thomas White's avatar
Thomas White committed
367
			case 3 :
368
369
370
371
372
373
			if ( sscanf(optarg, "%e", &rmax_fix) != 1 ) {
				ERROR("Invalid value for --rmax\n");
				return 1;
			}
			break;

Thomas White's avatar
Thomas White committed
374
			case 4 :
Thomas White's avatar
Thomas White committed
375
376
377
378
379
380
			if ( sscanf(optarg, "%f", &sigma_cutoff) != 1 ) {
				ERROR("Invalid value for --sigma-cutoff\n");
				return 1;
			}
			break;

Thomas White's avatar
Thomas White committed
381
			default :
382
383
			ERROR("Unhandled option '%c'\n", c);
			break;
384
385
386
387
388
		}

	}

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

Thomas White's avatar
Thomas White committed
393
394
	if ( sym_str == NULL ) {
		sym_str = strdup("1");
395
	}
Thomas White's avatar
Thomas White committed
396
397
	sym = get_pointgroup(sym_str);
	free(sym_str);
398
399
400
401

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

	if ( pdb == NULL ) {
Thomas White's avatar
Thomas White committed
402
403
404
		ERROR("You need to provide a PDB file containing"
		       " the unit cell.\n");
		return 1;
405
406
407
408
	}
	cell = load_cell_from_pdb(pdb);
	free(pdb);

Thomas White's avatar
Thomas White committed
409
410
411
	raw_list = read_reflections(file);
	if ( raw_list == NULL ) {
		ERROR("Couldn't read file '%s'\n", file);
412
413
414
		return 1;
	}

Thomas White's avatar
Thomas White committed
415
	/* Check that the intensities have the correct symmetry */
Thomas White's avatar
Thomas White committed
416
	if ( check_list_symmetry(raw_list, sym) ) {
Thomas White's avatar
Thomas White committed
417
		ERROR("The input reflection list does not appear to"
Thomas White's avatar
Thomas White committed
418
		      " have symmetry %s\n", symmetry_name(sym));
Thomas White's avatar
Thomas White committed
419
420
421
		return 1;
	}

Thomas White's avatar
Thomas White committed
422
423
	/* Reject some reflections */
	list = reflist_new();
Thomas White's avatar
Thomas White committed
424
	for ( refl = first_refl(raw_list, &iter);
Thomas White's avatar
Thomas White committed
425
426
	      refl != NULL;
	      refl = next_refl(refl, iter) ) {
427
428
429
430

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

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

Thomas White's avatar
Thomas White committed
435
436
		val = get_intensity(refl);
		sig = get_esd_intensity(refl);
437

Thomas White's avatar
Thomas White committed
438
		if ( val < sigma_cutoff * sig ) {
439
440
441
442
			rej++;
			ig = 1;
		}

Thomas White's avatar
Thomas White committed
443
		if ( ig ) continue;
444

Thomas White's avatar
Thomas White committed
445
446
		new = add_refl(list, h, k, l);
		copy_data(new, refl);
447
448

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

Thomas White's avatar
Thomas White committed
452
	plot_shells(list, cell, sym, rmin_fix, rmax_fix);
453
454
455

	return 0;
}