get_hkl.c 10.7 KB
Newer Older
1
2
3
/*
 * get_hkl.c
 *
4
 * Small program to manipulate reflection lists
5
 *
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
 *
 * 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"
Thomas White's avatar
Thomas White committed
25
#include "reflist-utils.h"
26
#include "symmetry.h"
27
#include "beam-parameters.h"
28
29
30
31
32
33


static void show_help(const char *s)
{
	printf("Syntax: %s [options]\n\n", s);
	printf(
34
"Manipulate reflection lists.\n"
35
"\n"
Thomas White's avatar
Thomas White committed
36
37
"  -h, --help                 Display this help message.\n"
"\n"
38
39
40
41
"  -i, --input=<file>         Read reflections from <file>.\n"
"  -y, --symmetry=<sym>       The symmetry of the input reflection list.\n"
"\n"
"You can add noise to the reflections with either of:\n"
Thomas White's avatar
Thomas White committed
42
"      --poisson              Simulate Poisson samples.\n"
43
"      --noise                Add 10%% random noise.\n"
44
45
46
47
"\n"
"To calculate Poisson samples accurately, you must also give:\n"
"  -b, --beam=<file>          Get beam parameters from file.\n"
"\n"
Thomas White's avatar
Thomas White committed
48
49
"You can artificially 'twin' the reflections, or expand them out.  You can also"
" do both, in which case the 'twinning' will be done first:\n"
50
51
"  -w, --twin=<sym>           Generate twinned data according to the given\n"
"                              point group.\n"
Thomas White's avatar
Thomas White committed
52
"  -e, --expand=<sym>         Expand reflections to this point group.\n"
53
54
55
56
57
"\n"
"You can restrict which reflections are written out:\n"
"  -t, --template=<filename>  Only include reflections mentioned in file.\n"
"\n"
"You might sometimes need to do this:\n"
58
59
"      --multiplicity         Multiply intensities by the number of\n"
"                              equivalent reflections.\n"
60
61
62
"\n"
"Don't forget to specify the output filename:\n"
"  -o, --output=<filename>    Output filename (default: stdout).\n"
Thomas White's avatar
Thomas White committed
63
64
"  -p, --pdb=<filename>       Use unit cell parameters from this PDB file to\n"
"                              generate resolution values in the output file.\n"
Thomas White's avatar
Thomas White committed
65
);
Thomas White's avatar
Thomas White committed
66
67
68
69
}


/* Apply Poisson noise to all reflections */
Thomas White's avatar
Thomas White committed
70
static void poisson_reflections(RefList *list, double adu_per_photon)
Thomas White's avatar
Thomas White committed
71
{
Thomas White's avatar
Thomas White committed
72
73
	Reflection *refl;
	RefListIterator *iter;
Thomas White's avatar
Thomas White committed
74

Thomas White's avatar
Thomas White committed
75
76
77
	for ( refl = first_refl(list, &iter);
	      refl != NULL;
	      refl = next_refl(refl, iter) ) {
Thomas White's avatar
Thomas White committed
78

Thomas White's avatar
Thomas White committed
79
		double val, c;
Thomas White's avatar
Thomas White committed
80

Thomas White's avatar
Thomas White committed
81
		val = get_intensity(refl);
82

83
		c = adu_per_photon * poisson_noise(val/adu_per_photon);
Thomas White's avatar
Thomas White committed
84
		set_int(refl, c);
Thomas White's avatar
Thomas White committed
85
86

	}
87
88
89
90
}


/* Apply 10% uniform noise to all reflections */
Thomas White's avatar
Thomas White committed
91
static void noise_reflections(RefList *list)
92
{
Thomas White's avatar
Thomas White committed
93
94
	Reflection *refl;
	RefListIterator *iter;
95

Thomas White's avatar
Thomas White committed
96
97
98
	for ( refl = first_refl(list, &iter);
	      refl != NULL;
	      refl = next_refl(refl, iter) ) {
99

Thomas White's avatar
Thomas White committed
100
		double val, r;
101

Thomas White's avatar
Thomas White committed
102
		val = get_intensity(refl);
103
104
105
106

		r = (double)random()/RAND_MAX;
		val += 0.1 * val * r;

Thomas White's avatar
Thomas White committed
107
		set_int(refl, val);
108

Thomas White's avatar
Thomas White committed
109
	}
110
111
112
}


Thomas White's avatar
Thomas White committed
113
static RefList *template_reflections(RefList *list, RefList *template)
114
{
Thomas White's avatar
Thomas White committed
115
116
117
	Reflection *refl;
	RefListIterator *iter;
	RefList *out;
118

Thomas White's avatar
Thomas White committed
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
	out = reflist_new();

	for ( refl = first_refl(template, &iter);
	      refl != NULL;
	      refl = next_refl(refl, iter) ) {

		signed int h, k, l;
		Reflection *new;
		Reflection *old;

		get_indices(refl, &h, &k, &l);

		old = find_refl(list, h, k, l);
		if ( old == NULL ) continue;

		new = add_refl(out, h, k, l);
		copy_data(new, old);

	}

	return out;
}


static RefList *twin_reflections(RefList *in,
Thomas White's avatar
Thomas White committed
144
                                 const SymOpList *holo, const SymOpList *mero)
Thomas White's avatar
Thomas White committed
145
146
147
148
{
	Reflection *refl;
	RefListIterator *iter;
	RefList *out;
149

Thomas White's avatar
Thomas White committed
150
151
152
153
	/* FIXME: Check properly by coset decomposition */
	if ( num_equivs(holo) < num_equivs(mero) ) {
		ERROR("%s is not a subgroup of %s!\n", symmetry_name(mero),
		                                       symmetry_name(holo));
154
155
156
		return NULL;
	}

Thomas White's avatar
Thomas White committed
157
158
159
160
161
	out = reflist_new();

	for ( refl = first_refl(in, &iter);
	      refl != NULL;
	      refl = next_refl(refl, iter) ) {
162

163
		double total, sigma;
164
165
166
167
		signed int h, k, l;
		int n, j;
		int skip;

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

170
		/* There is a many-to-one correspondence between reflections
Thomas White's avatar
Thomas White committed
171
172
173
		 * in the merohedral and holohedral groups.  Do the calculation
		 * only once for each reflection in the holohedral group, which
		 * contains fewer reflections.
174
		 */
Thomas White's avatar
Thomas White committed
175
		get_asymm(holo, h, k, l, &h, &k, &l);
Thomas White's avatar
Thomas White committed
176
		if ( find_refl(out, h, k, l) != NULL ) continue;
177

178
179
		total = 0.0;
		sigma = 0.0;
180
		skip = 0;
Thomas White's avatar
Thomas White committed
181
		n = num_equivs(holo);
182
183
184
185
186
		for ( j=0; j<n; j++ ) {

			signed int he, ke, le;
			signed int hu, ku, lu;

Thomas White's avatar
Thomas White committed
187
			get_equiv(holo, j, h, k, l, &he, &ke, &le);
188
189
190
191
192
193

			/* Do we have this reflection?
			 * We might not have the particular (merohedral)
			 * equivalent which belongs to our definition of the
			 * asymmetric unit cell, so check them all.
			 */
Thomas White's avatar
Thomas White committed
194
195
			if ( !find_equiv_in_list(in, he, ke, le, mero,
			                         &hu, &ku, &lu) ) {
196
197
198
199
200
				/* Don't have this reflection, so bail out */
				ERROR("Twinning %i %i %i requires the %i %i %i "
				      "reflection (or an equivalent in %s), "
				      "which I don't have. %i %i %i won't "
				      "appear in the output\n",
Thomas White's avatar
Thomas White committed
201
202
				      h, k, l, he, ke, le, symmetry_name(mero),
				      h, k, l);
203
204
205
206
				skip = 1;
				break;
			}

Thomas White's avatar
Thomas White committed
207
208
			total += get_intensity(refl);
			sigma += pow(get_esd_intensity(refl), 2.0);
209
210
211
212
213

		}

		if ( !skip ) {

Thomas White's avatar
Thomas White committed
214
215
216
			Reflection *new = add_refl(out, h, k, l);
			set_int(new, total);
			set_esd_intensity(new, sqrt(sigma));
217
218
219
220
221

		}

	}

Thomas White's avatar
Thomas White committed
222
	return out;
223
224
225
}


Thomas White's avatar
Thomas White committed
226
227
static RefList *expand_reflections(RefList *in, const SymOpList *target,
                                                const SymOpList *initial)
Thomas White's avatar
Thomas White committed
228
{
Thomas White's avatar
Thomas White committed
229
230
231
	Reflection *refl;
	RefListIterator *iter;
	RefList *out;
Thomas White's avatar
Thomas White committed
232

Thomas White's avatar
Thomas White committed
233
234
235
236
	/* FIXME: Check properly */
	if ( num_equivs(target) > num_equivs(initial) ) {
		ERROR("%s is not a subgroup of %s!\n", symmetry_name(initial),
		                                       symmetry_name(target));
237
238
239
		return NULL;
	}

Thomas White's avatar
Thomas White committed
240
241
242
243
244
	out = reflist_new();

	for ( refl = first_refl(in, &iter);
	      refl != NULL;
	      refl = next_refl(refl, iter) ) {
Thomas White's avatar
Thomas White committed
245
246
247
248
249

		signed int h, k, l;
		int n, j;
		double intensity;

Thomas White's avatar
Thomas White committed
250
251
		get_indices(refl, &h, &k, &l);
		intensity = get_intensity(refl);
Thomas White's avatar
Thomas White committed
252

Thomas White's avatar
Thomas White committed
253
		n = num_equivs(initial);
Thomas White's avatar
Thomas White committed
254

Thomas White's avatar
Thomas White committed
255
		/* For each equivalent in the higher symmetry group */
Thomas White's avatar
Thomas White committed
256
257
258
		for ( j=0; j<n; j++ ) {

			signed int he, ke, le;
Thomas White's avatar
Thomas White committed
259
			Reflection *new;
Thomas White's avatar
Thomas White committed
260
261

			/* Get the equivalent */
Thomas White's avatar
Thomas White committed
262
			get_equiv(initial, j, h, k, l, &he, &ke, &le);
Thomas White's avatar
Thomas White committed
263
264

			/* Put it into the asymmetric unit for the target */
Thomas White's avatar
Thomas White committed
265
			get_asymm(target, he, ke, le, &he, &ke, &le);
Thomas White's avatar
Thomas White committed
266
267

			/* Make sure the intensity is in the right place */
Thomas White's avatar
Thomas White committed
268
269
			new = add_refl(out, he, ke, le);
			copy_data(new, refl);
Thomas White's avatar
Thomas White committed
270
271
272
273
274

		}

	}

Thomas White's avatar
Thomas White committed
275
	return out;
Thomas White's avatar
Thomas White committed
276
277
278
}


279
280
281
int main(int argc, char *argv[])
{
	int c;
282
283
	int config_noise = 0;
	int config_poisson = 0;
284
	int config_multi = 0;
Thomas White's avatar
Thomas White committed
285
286
287
288
289
290
	char *holo_str = NULL;
	char *mero_str = NULL;
	char *expand_str = NULL;
	SymOpList *holo;
	SymOpList *mero;
	SymOpList *expand;
Thomas White's avatar
Thomas White committed
291
292
	char *input_file = NULL;
	char *template = NULL;
Thomas White's avatar
Thomas White committed
293
	char *output = NULL;
294
	char *beamfile = NULL;
295
	struct beam_params *beam = NULL;
Thomas White's avatar
Thomas White committed
296
297
	RefList *input;
	UnitCell *cell = NULL;
298
299
300
301

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
Thomas White's avatar
Thomas White committed
302
		{"template",           1, NULL,               't'},
303
304
		{"poisson",            0, &config_poisson,     1},
		{"noise",              0, &config_noise,       1},
Thomas White's avatar
Thomas White committed
305
		{"output",             1, NULL,               'o'},
306
		{"symmetry",           1, NULL,               'y'},
Thomas White's avatar
Thomas White committed
307
308
		{"twin",               1, NULL,               'w'},
		{"expand",             1, NULL,               'e'},
309
		{"intensities",        1, NULL,               'i'},
310
		{"multiplicity",       0, &config_multi,       1},
311
		{"beam",               1, NULL,               'b'},
Thomas White's avatar
Thomas White committed
312
		{"pdb",                1, NULL,               'p'},
313
314
315
316
		{0, 0, NULL, 0}
	};

	/* Short options */
Thomas White's avatar
Thomas White committed
317
	while ((c = getopt_long(argc, argv, "ht:o:i:w:y:e:b:p:",
Thomas White's avatar
Thomas White committed
318
	                        longopts, NULL)) != -1) {
319
320

		switch (c) {
Thomas White's avatar
Thomas White committed
321
		case 'h' :
322
323
324
			show_help(argv[0]);
			return 0;

Thomas White's avatar
Thomas White committed
325
		case 't' :
Thomas White's avatar
Thomas White committed
326
327
328
			template = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
329
		case 'o' :
Thomas White's avatar
Thomas White committed
330
331
332
			output = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
333
		case 'i' :
Thomas White's avatar
Thomas White committed
334
			input_file = strdup(optarg);
335
336
			break;

337
		case 'y' :
Thomas White's avatar
Thomas White committed
338
			mero_str = strdup(optarg);
339
340
341
			break;

		case 'w' :
Thomas White's avatar
Thomas White committed
342
			holo_str = strdup(optarg);
343
344
			break;

Thomas White's avatar
Thomas White committed
345
		case 'e' :
Thomas White's avatar
Thomas White committed
346
			expand_str = strdup(optarg);
Thomas White's avatar
Thomas White committed
347
348
			break;

349
350
351
352
		case 'b' :
			beamfile = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
353
354
355
356
357
358
359
360
		case 'p' :
			cell = load_cell_from_pdb(optarg);
			if ( cell == NULL ) {
				ERROR("Failed to get cell from '%s'\n", optarg);
				return 1;
			}
			break;

Thomas White's avatar
Thomas White committed
361
		case 0 :
362
363
			break;

Thomas White's avatar
Thomas White committed
364
		default :
365
366
367
368
369
			return 1;
		}

	}

Thomas White's avatar
Thomas White committed
370
	if ( (holo_str != NULL) && (expand_str != NULL) ) {
Thomas White's avatar
Thomas White committed
371
372
		ERROR("You cannot 'twin' and 'expand' at the same time.\n");
		ERROR("Decide which one you want to do first.\n");
373
374
375
376
377
378
379
380
381
382
		return 1;
	}

	if ( beamfile != NULL ) {
		beam = get_beam_parameters(beamfile);
		if ( beam == NULL ) {
			ERROR("Failed to load beam parameters from '%s'\n",
			      beamfile);
			return 1;
		}
Thomas White's avatar
Thomas White committed
383
384
	}

Thomas White's avatar
Thomas White committed
385
386
387
	if ( cell == NULL ) {
		ERROR("You need to give a PDB file with the unit cell.\n");
		return 1;
388
	}
389

Thomas White's avatar
Thomas White committed
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
	if ( holo_str != NULL ) {
		holo = get_pointgroup(holo_str);
		free(holo_str);
	} else {
		holo = NULL;
	}
	if ( mero_str != NULL ) {
		mero = get_pointgroup(mero_str);
		free(mero_str);
	} else {
		mero = NULL;
	}
	if ( expand_str != NULL ) {
		expand = get_pointgroup(expand_str);
		free(expand_str);
	} else {
		expand = NULL;
	}

Thomas White's avatar
Thomas White committed
409
	input = read_reflections(input_file);
Thomas White's avatar
Thomas White committed
410
	free(input_file);
Thomas White's avatar
Thomas White committed
411
	if ( check_list_symmetry(input, mero) ) {
412
		ERROR("The input reflection list does not appear to"
Thomas White's avatar
Thomas White committed
413
		      " have symmetry %s\n", symmetry_name(mero));
414
		return 1;
415
	}
416

417
418
	if ( config_poisson ) {
		if ( beam != NULL ) {
Thomas White's avatar
Thomas White committed
419
			poisson_reflections(input, beam->adu_per_photon);
420
421
422
423
424
425
426
		} else {
			ERROR("You must give a beam parameters file in order"
			      " to calculate Poisson noise.\n");
			return 1;
		}
	}

Thomas White's avatar
Thomas White committed
427
	if ( config_noise ) noise_reflections(input);
Thomas White's avatar
Thomas White committed
428

429
	if ( holo != NULL ) {
430

Thomas White's avatar
Thomas White committed
431
		RefList *new;
Thomas White's avatar
Thomas White committed
432
433
		STATUS("Twinning from %s into %s\n", symmetry_name(mero),
		                                     symmetry_name(holo));
Thomas White's avatar
Thomas White committed
434
435
436
437
438
439
440
441
442
		new = twin_reflections(input, holo, mero);

		/* Replace old with new */
		reflist_free(input);
		input = new;

		/* The symmetry of the list has changed */
		free(mero);
		mero = holo;
443

Thomas White's avatar
Thomas White committed
444
445
	}

Thomas White's avatar
Thomas White committed
446
	if ( expand != NULL ) {
447

Thomas White's avatar
Thomas White committed
448
		RefList *new;
Thomas White's avatar
Thomas White committed
449
450
		STATUS("Expanding from %s into %s\n", symmetry_name(mero),
		                                      symmetry_name(expand));
Thomas White's avatar
Thomas White committed
451
452
453
454
455
		new = expand_reflections(input, expand, mero);

		/* Replace old with new */
		reflist_free(input);
		input = new;
456

Thomas White's avatar
Thomas White committed
457
458
	}

459
460
	if ( config_multi ) {

Thomas White's avatar
Thomas White committed
461
462
		Reflection *refl;
		RefListIterator *iter;
463

Thomas White's avatar
Thomas White committed
464
465
466
		for ( refl = first_refl(input, &iter);
		      refl != NULL;
		      refl = next_refl(refl, iter) ) {
467
468

			double inty;
Thomas White's avatar
Thomas White committed
469
470
471
472
			signed int h, k, l;

			get_indices(refl, &h, &k, &l);
			inty = get_intensity(refl);
473

Thomas White's avatar
Thomas White committed
474
			inty *= (double)num_equivs(mero);
Thomas White's avatar
Thomas White committed
475
			set_int(refl, inty);
476
477
478
479

		}
	}

Thomas White's avatar
Thomas White committed
480
	if ( template ) {
481

Thomas White's avatar
Thomas White committed
482
483
484
485
		RefList *t = read_reflections(template);
		RefList *new = template_reflections(input, t);
		reflist_free(input);
		input = new;
486

Thomas White's avatar
Thomas White committed
487
488
	}

Thomas White's avatar
Thomas White committed
489
	write_reflist(output, input, cell);
Thomas White's avatar
Thomas White committed
490

Thomas White's avatar
Thomas White committed
491
	reflist_free(input);
492
493
494

	return 0;
}