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
286
	if ( nsilly ) {
		STATUS("Warning; %i reflections had infinite or invalid values"
		       " of I/sigma(I).\n", nsilly);
	}


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

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

	}

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

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


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

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

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

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

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

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

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

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

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

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

	}

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

	return 0;
}