get_hkl.c 11.2 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
278
279
280
281

		}

	}

Thomas White's avatar
Thomas White committed
282
283
	free_symopmask(m);

Thomas White's avatar
Thomas White committed
284
	return out;
Thomas White's avatar
Thomas White committed
285
286
287
}


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

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

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

		switch (c) {
Thomas White's avatar
Thomas White committed
330
		case 'h' :
331
332
333
			show_help(argv[0]);
			return 0;

Thomas White's avatar
Thomas White committed
334
		case 't' :
Thomas White's avatar
Thomas White committed
335
336
337
			template = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
338
		case 'o' :
Thomas White's avatar
Thomas White committed
339
340
341
			output = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
342
		case 'i' :
Thomas White's avatar
Thomas White committed
343
			input_file = strdup(optarg);
344
345
			break;

346
		case 'y' :
Thomas White's avatar
Thomas White committed
347
			mero_str = strdup(optarg);
348
349
350
			break;

		case 'w' :
Thomas White's avatar
Thomas White committed
351
			holo_str = strdup(optarg);
352
353
			break;

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

358
359
360
361
		case 'b' :
			beamfile = strdup(optarg);
			break;

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

Thomas White's avatar
Thomas White committed
373
		default :
374
375
376
377
378
			return 1;
		}

	}

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

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

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

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

442
	if ( holo != NULL ) {
443

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

Thomas White's avatar
Thomas White committed
457
458
	}

Thomas White's avatar
Thomas White committed
459
	if ( expand != NULL ) {
460

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

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

Thomas White's avatar
Thomas White committed
470
471
	}

472
473
	if ( config_multi ) {

Thomas White's avatar
Thomas White committed
474
475
		Reflection *refl;
		RefListIterator *iter;
Thomas White's avatar
Thomas White committed
476
477
478
		SymOpMask *m;

		m = new_symopmask(mero);
479

Thomas White's avatar
Thomas White committed
480
481
482
		for ( refl = first_refl(input, &iter);
		      refl != NULL;
		      refl = next_refl(refl, iter) ) {
483
484

			double inty;
Thomas White's avatar
Thomas White committed
485
486
487
488
			signed int h, k, l;

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

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

		}
Thomas White's avatar
Thomas White committed
495
496

		free_symopmask(m);
497
498
	}

Thomas White's avatar
Thomas White committed
499
	if ( template ) {
500

Thomas White's avatar
Thomas White committed
501
502
503
504
		RefList *t = read_reflections(template);
		RefList *new = template_reflections(input, t);
		reflist_free(input);
		input = new;
505

Thomas White's avatar
Thomas White committed
506
507
	}

Thomas White's avatar
Thomas White committed
508
	write_reflist(output, input, cell);
Thomas White's avatar
Thomas White committed
509

Thomas White's avatar
Thomas White committed
510
	reflist_free(input);
511
512
513

	return 0;
}