get_hkl.c 11.4 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;
267
268
			int have_phase;
			double ph;
Thomas White's avatar
Thomas White committed
269
270

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

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

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

Thomas White's avatar
Thomas White committed
280
281
			/* FIXME: Make phase negative if the reflection is
			 * separated from the original via an inversion */
282
283
284
285
			get_phase(refl, &have_phase);
			if ( have_phase ) {
				set_ph(new, -ph);
			}
Thomas White's avatar
Thomas White committed
286

Thomas White's avatar
Thomas White committed
287
288
289
290
		}

	}

Thomas White's avatar
Thomas White committed
291
292
	free_symopmask(m);

Thomas White's avatar
Thomas White committed
293
	return out;
Thomas White's avatar
Thomas White committed
294
295
296
}


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

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

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

		switch (c) {
Thomas White's avatar
Thomas White committed
339
		case 'h' :
340
341
342
			show_help(argv[0]);
			return 0;

Thomas White's avatar
Thomas White committed
343
		case 't' :
Thomas White's avatar
Thomas White committed
344
345
346
			template = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
347
		case 'o' :
Thomas White's avatar
Thomas White committed
348
349
350
			output = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
351
		case 'i' :
Thomas White's avatar
Thomas White committed
352
			input_file = strdup(optarg);
353
354
			break;

355
		case 'y' :
Thomas White's avatar
Thomas White committed
356
			mero_str = strdup(optarg);
357
358
359
			break;

		case 'w' :
Thomas White's avatar
Thomas White committed
360
			holo_str = strdup(optarg);
361
362
			break;

Thomas White's avatar
Thomas White committed
363
		case 'e' :
Thomas White's avatar
Thomas White committed
364
			expand_str = strdup(optarg);
Thomas White's avatar
Thomas White committed
365
366
			break;

367
368
369
370
		case 'b' :
			beamfile = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
371
372
373
374
375
376
377
378
		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
379
		case 0 :
380
381
			break;

Thomas White's avatar
Thomas White committed
382
		default :
383
384
385
386
387
			return 1;
		}

	}

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

Thomas White's avatar
Thomas White committed
403
404
405
	if ( cell == NULL ) {
		ERROR("You need to give a PDB file with the unit cell.\n");
		return 1;
406
	}
407

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

439
440
	if ( config_poisson ) {
		if ( beam != NULL ) {
Thomas White's avatar
Thomas White committed
441
			poisson_reflections(input, beam->adu_per_photon);
442
443
444
445
446
447
448
		} 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
449
	if ( config_noise ) noise_reflections(input);
Thomas White's avatar
Thomas White committed
450

451
	if ( holo != NULL ) {
452

Thomas White's avatar
Thomas White committed
453
		RefList *new;
Thomas White's avatar
Thomas White committed
454
455
		STATUS("Twinning from %s into %s\n", symmetry_name(mero),
		                                     symmetry_name(holo));
Thomas White's avatar
Thomas White committed
456
457
458
459
460
461
462
463
464
		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;
465

Thomas White's avatar
Thomas White committed
466
467
	}

Thomas White's avatar
Thomas White committed
468
	if ( expand != NULL ) {
469

Thomas White's avatar
Thomas White committed
470
		RefList *new;
Thomas White's avatar
Thomas White committed
471
472
		STATUS("Expanding from %s into %s\n", symmetry_name(mero),
		                                      symmetry_name(expand));
Thomas White's avatar
Thomas White committed
473
474
475
476
477
		new = expand_reflections(input, expand, mero);

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

Thomas White's avatar
Thomas White committed
479
480
	}

481
482
	if ( config_multi ) {

Thomas White's avatar
Thomas White committed
483
484
		Reflection *refl;
		RefListIterator *iter;
Thomas White's avatar
Thomas White committed
485
486
487
		SymOpMask *m;

		m = new_symopmask(mero);
488

Thomas White's avatar
Thomas White committed
489
490
491
		for ( refl = first_refl(input, &iter);
		      refl != NULL;
		      refl = next_refl(refl, iter) ) {
492
493

			double inty;
Thomas White's avatar
Thomas White committed
494
495
496
497
			signed int h, k, l;

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

Thomas White's avatar
Thomas White committed
499
500
			special_position(mero, m, h, k, l);
			inty *= (double)num_equivs(mero, m);
Thomas White's avatar
Thomas White committed
501
			set_int(refl, inty);
502
503

		}
Thomas White's avatar
Thomas White committed
504
505

		free_symopmask(m);
506
507
	}

Thomas White's avatar
Thomas White committed
508
	if ( template ) {
509

Thomas White's avatar
Thomas White committed
510
511
512
513
		RefList *t = read_reflections(template);
		RefList *new = template_reflections(input, t);
		reflist_free(input);
		input = new;
514

Thomas White's avatar
Thomas White committed
515
516
	}

Thomas White's avatar
Thomas White committed
517
	write_reflist(output, input, cell);
Thomas White's avatar
Thomas White committed
518

Thomas White's avatar
Thomas White committed
519
	reflist_free(input);
520
521
522

	return 0;
}