check_hkl.c 9.83 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
262
263

		if ( isfinite(val/esd) ) {
			snr[bin] += val / esd;
			snr_total += val / esd;
			snr_measured[bin]++;
		}

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

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

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

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

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


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

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

	}

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

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


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

	/* 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
333
334
335
		{"rmin",               1, NULL,                2},
		{"rmax",               1, NULL,                3},
		{"sigma-cutoff",       1, NULL,                4},
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) {
Thomas White's avatar
Thomas White committed
343
344

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

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

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

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

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

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

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

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

	}

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

	return 0;
}