check_hkl.c 9.47 KB
Newer Older
1
2
3
4
5
/*
 * check_hkl.c
 *
 * Characterise reflection lists
 *
6
7
8
9
10
 * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY,
 *                  a research centre of the Helmholtz Association.
 *
 * Authors:
 *   2009-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
47
48
49
50
51
52
53
54
55
56
57
58


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


Thomas White's avatar
Thomas White committed
68
static void plot_shells(RefList *list, UnitCell *cell, const SymOpList *sym,
Thomas White's avatar
Thomas White committed
69
                        double rmin_fix, double rmax_fix)
70
71
72
73
74
75
76
77
78
79
{
	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
80
81
	double mean[NBINS];
	double var[NBINS];
82
83
84
85
86
87
88
89
	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
90
91
	Reflection *refl;
	RefListIterator *iter;
92
	RefList *counted;
93
94
95
96
	int hmax, kmax, lmax;
	double asx, asy, asz;
	double bsx, bsy, bsz;
	double csx, csy, csz;
97
98
99
100
101
102
103
104
105
106
107
108
109
110

	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
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
196
197
198
199
		signed int h, k, l;
		double d;
		int bin;
		int j;

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

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

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

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

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

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

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

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

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

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

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

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

	}

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

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


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

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

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

	}

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

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

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

	if ( pdb == NULL ) {
Thomas White's avatar
Thomas White committed
384
385
386
		ERROR("You need to provide a PDB file containing"
		       " the unit cell.\n");
		return 1;
387
388
389
390
	}
	cell = load_cell_from_pdb(pdb);
	free(pdb);

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

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

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

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

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

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

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

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

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

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

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

	return 0;
}