get_hkl.c 11.3 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;
Thomas White's avatar
Thomas White committed
149
	SymOpMask *m;
150

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

Thomas White's avatar
Thomas White committed
158
	out = reflist_new();
Thomas White's avatar
Thomas White committed
159
	m = new_symopmask(holo);
Thomas White's avatar
Thomas White committed
160
161
162
163

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

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

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

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

180
181
		total = 0.0;
		sigma = 0.0;
182
		skip = 0;
Thomas White's avatar
Thomas White committed
183
184
185
		special_position(holo, m, h, k, l);
		n = num_equivs(holo, m);

186
187
188
189
190
		for ( j=0; j<n; j++ ) {

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

Thomas White's avatar
Thomas White committed
191
			get_equiv(holo, m, j, h, k, l, &he, &ke, &le);
192
193
194
195
196
197

			/* 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
198
199
			if ( !find_equiv_in_list(in, he, ke, le, mero,
			                         &hu, &ku, &lu) ) {
200
201
202
203
204
				/* 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
205
206
				      h, k, l, he, ke, le, symmetry_name(mero),
				      h, k, l);
207
208
209
210
				skip = 1;
				break;
			}

Thomas White's avatar
Thomas White committed
211
212
			total += get_intensity(refl);
			sigma += pow(get_esd_intensity(refl), 2.0);
213
214
215
216
217

		}

		if ( !skip ) {

Thomas White's avatar
Thomas White committed
218
219
220
			Reflection *new = add_refl(out, h, k, l);
			set_int(new, total);
			set_esd_intensity(new, sqrt(sigma));
221
222
223
224
225

		}

	}

Thomas White's avatar
Thomas White committed
226
	return out;
227
228
229
}


Thomas White's avatar
Thomas White committed
230
231
static RefList *expand_reflections(RefList *in, const SymOpList *target,
                                                const SymOpList *initial)
Thomas White's avatar
Thomas White committed
232
{
Thomas White's avatar
Thomas White committed
233
234
235
	Reflection *refl;
	RefListIterator *iter;
	RefList *out;
Thomas White's avatar
Thomas White committed
236
	SymOpMask *m;
Thomas White's avatar
Thomas White committed
237

Thomas White's avatar
Thomas White committed
238
	/* FIXME: Check properly */
Thomas White's avatar
Thomas White committed
239
	if ( num_equivs(target, NULL) > num_equivs(initial, NULL) ) {
Thomas White's avatar
Thomas White committed
240
241
		ERROR("%s is not a subgroup of %s!\n", symmetry_name(initial),
		                                       symmetry_name(target));
242
243
244
		return NULL;
	}

Thomas White's avatar
Thomas White committed
245
	out = reflist_new();
Thomas White's avatar
Thomas White committed
246
	m = new_symopmask(initial);
Thomas White's avatar
Thomas White committed
247
248
249
250

	for ( refl = first_refl(in, &iter);
	      refl != NULL;
	      refl = next_refl(refl, iter) ) {
Thomas White's avatar
Thomas White committed
251
252
253
254
255

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

Thomas White's avatar
Thomas White committed
256
257
		get_indices(refl, &h, &k, &l);
		intensity = get_intensity(refl);
Thomas White's avatar
Thomas White committed
258

Thomas White's avatar
Thomas White committed
259
260
		special_position(initial, m, h, k, l);
		n = num_equivs(initial, m);
Thomas White's avatar
Thomas White committed
261

Thomas White's avatar
Thomas White committed
262
		/* For each equivalent in the higher symmetry group */
Thomas White's avatar
Thomas White committed
263
264
265
		for ( j=0; j<n; j++ ) {

			signed int he, ke, le;
Thomas White's avatar
Thomas White committed
266
			Reflection *new;
Thomas White's avatar
Thomas White committed
267
268

			/* Get the equivalent */
Thomas White's avatar
Thomas White committed
269
			get_equiv(initial, m, j, h, k, l, &he, &ke, &le);
Thomas White's avatar
Thomas White committed
270
271

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

			/* Make sure the intensity is in the right place */
Thomas White's avatar
Thomas White committed
275
276
			new = add_refl(out, he, ke, le);
			copy_data(new, refl);
Thomas White's avatar
Thomas White committed
277

Thomas White's avatar
Thomas White committed
278
279
280
			/* FIXME: Make phase negative if the reflection is
			 * separated from the original via an inversion */

Thomas White's avatar
Thomas White committed
281
282
283
284
		}

	}

Thomas White's avatar
Thomas White committed
285
286
	free_symopmask(m);

Thomas White's avatar
Thomas White committed
287
	return out;
Thomas White's avatar
Thomas White committed
288
289
290
}


291
292
293
int main(int argc, char *argv[])
{
	int c;
294
295
	int config_noise = 0;
	int config_poisson = 0;
296
	int config_multi = 0;
Thomas White's avatar
Thomas White committed
297
298
299
300
301
302
	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
303
304
	char *input_file = NULL;
	char *template = NULL;
Thomas White's avatar
Thomas White committed
305
	char *output = NULL;
306
	char *beamfile = NULL;
307
	struct beam_params *beam = NULL;
Thomas White's avatar
Thomas White committed
308
309
	RefList *input;
	UnitCell *cell = NULL;
310
311
312
313

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
Thomas White's avatar
Thomas White committed
314
		{"template",           1, NULL,               't'},
315
316
		{"poisson",            0, &config_poisson,     1},
		{"noise",              0, &config_noise,       1},
Thomas White's avatar
Thomas White committed
317
		{"output",             1, NULL,               'o'},
318
		{"symmetry",           1, NULL,               'y'},
Thomas White's avatar
Thomas White committed
319
320
		{"twin",               1, NULL,               'w'},
		{"expand",             1, NULL,               'e'},
321
		{"intensities",        1, NULL,               'i'},
322
		{"multiplicity",       0, &config_multi,       1},
323
		{"beam",               1, NULL,               'b'},
Thomas White's avatar
Thomas White committed
324
		{"pdb",                1, NULL,               'p'},
325
326
327
328
		{0, 0, NULL, 0}
	};

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

		switch (c) {
Thomas White's avatar
Thomas White committed
333
		case 'h' :
334
335
336
			show_help(argv[0]);
			return 0;

Thomas White's avatar
Thomas White committed
337
		case 't' :
Thomas White's avatar
Thomas White committed
338
339
340
			template = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
341
		case 'o' :
Thomas White's avatar
Thomas White committed
342
343
344
			output = strdup(optarg);
			break;

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

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

		case 'w' :
Thomas White's avatar
Thomas White committed
354
			holo_str = strdup(optarg);
355
356
			break;

Thomas White's avatar
Thomas White committed
357
		case 'e' :
Thomas White's avatar
Thomas White committed
358
			expand_str = strdup(optarg);
Thomas White's avatar
Thomas White committed
359
360
			break;

361
362
363
364
		case 'b' :
			beamfile = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
365
366
367
368
369
370
371
372
		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
373
		case 0 :
374
375
			break;

Thomas White's avatar
Thomas White committed
376
		default :
377
378
379
380
381
			return 1;
		}

	}

Thomas White's avatar
Thomas White committed
382
	if ( (holo_str != NULL) && (expand_str != NULL) ) {
Thomas White's avatar
Thomas White committed
383
384
		ERROR("You cannot 'twin' and 'expand' at the same time.\n");
		ERROR("Decide which one you want to do first.\n");
385
386
387
388
389
390
391
392
393
394
		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
395
396
	}

Thomas White's avatar
Thomas White committed
397
398
399
	if ( cell == NULL ) {
		ERROR("You need to give a PDB file with the unit cell.\n");
		return 1;
400
	}
401

Thomas White's avatar
Thomas White committed
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
	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
421
	input = read_reflections(input_file);
Thomas White's avatar
Thomas White committed
422
423
424
425
	if ( input == NULL ) {
		ERROR("Problem reading input file %s\n", input_file);
		return 1;
	}
Thomas White's avatar
Thomas White committed
426
	free(input_file);
Thomas White's avatar
Thomas White committed
427
	if ( check_list_symmetry(input, mero) ) {
428
		ERROR("The input reflection list does not appear to"
Thomas White's avatar
Thomas White committed
429
		      " have symmetry %s\n", symmetry_name(mero));
430
		return 1;
431
	}
432

433
434
	if ( config_poisson ) {
		if ( beam != NULL ) {
Thomas White's avatar
Thomas White committed
435
			poisson_reflections(input, beam->adu_per_photon);
436
437
438
439
440
441
442
		} 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
443
	if ( config_noise ) noise_reflections(input);
Thomas White's avatar
Thomas White committed
444

445
	if ( holo != NULL ) {
446

Thomas White's avatar
Thomas White committed
447
		RefList *new;
Thomas White's avatar
Thomas White committed
448
449
		STATUS("Twinning from %s into %s\n", symmetry_name(mero),
		                                     symmetry_name(holo));
Thomas White's avatar
Thomas White committed
450
451
452
453
454
455
456
457
458
		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;
459

Thomas White's avatar
Thomas White committed
460
461
	}

Thomas White's avatar
Thomas White committed
462
	if ( expand != NULL ) {
463

Thomas White's avatar
Thomas White committed
464
		RefList *new;
Thomas White's avatar
Thomas White committed
465
466
		STATUS("Expanding from %s into %s\n", symmetry_name(mero),
		                                      symmetry_name(expand));
Thomas White's avatar
Thomas White committed
467
468
469
470
471
		new = expand_reflections(input, expand, mero);

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

Thomas White's avatar
Thomas White committed
473
474
	}

475
476
	if ( config_multi ) {

Thomas White's avatar
Thomas White committed
477
478
		Reflection *refl;
		RefListIterator *iter;
Thomas White's avatar
Thomas White committed
479
480
481
		SymOpMask *m;

		m = new_symopmask(mero);
482

Thomas White's avatar
Thomas White committed
483
484
485
		for ( refl = first_refl(input, &iter);
		      refl != NULL;
		      refl = next_refl(refl, iter) ) {
486
487

			double inty;
Thomas White's avatar
Thomas White committed
488
489
490
491
			signed int h, k, l;

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

Thomas White's avatar
Thomas White committed
493
494
			special_position(mero, m, h, k, l);
			inty *= (double)num_equivs(mero, m);
Thomas White's avatar
Thomas White committed
495
			set_int(refl, inty);
496
497

		}
Thomas White's avatar
Thomas White committed
498
499

		free_symopmask(m);
500
501
	}

Thomas White's avatar
Thomas White committed
502
	if ( template ) {
503

Thomas White's avatar
Thomas White committed
504
505
506
507
		RefList *t = read_reflections(template);
		RefList *new = template_reflections(input, t);
		reflist_free(input);
		input = new;
508

Thomas White's avatar
Thomas White committed
509
510
	}

Thomas White's avatar
Thomas White committed
511
	write_reflist(output, input, cell);
Thomas White's avatar
Thomas White committed
512

Thomas White's avatar
Thomas White committed
513
	reflist_free(input);
514
515
516

	return 0;
}