get_hkl.c 13.2 KB
Newer Older
1
2
3
/*
 * get_hkl.c
 *
4
 * Small program to manipulate reflection lists
5
 *
6
7
8
9
10
 * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY,
 *                  a research centre of the Helmholtz Association.
 *
 * Authors:
 *   2009-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
 *
 */


#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
42
#include "reflist-utils.h"
43
#include "symmetry.h"
44
#include "beam-parameters.h"
45
46
47
48
49
50


static void show_help(const char *s)
{
	printf("Syntax: %s [options]\n\n", s);
	printf(
51
"Manipulate reflection lists.\n"
52
"\n"
Thomas White's avatar
Thomas White committed
53
54
"  -h, --help                 Display this help message.\n"
"\n"
55
56
57
58
"  -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
59
"      --poisson              Simulate Poisson samples.\n"
60
"      --noise                Add 10%% random noise.\n"
61
62
"\n"
"To calculate Poisson samples accurately, you must also give:\n"
63
"      --adu-per-photon=<n>   Number of ADU per photon.\n"
64
"\n"
65
"You can artificially 'twin' the reflections, or expand them out.\n"
66
67
"  -w, --twin=<sym>           Generate twinned data according to the given\n"
"                              point group.\n"
Thomas White's avatar
Thomas White committed
68
"  -e, --expand=<sym>         Expand reflections to this point group.\n"
69
"\n"
70
71
72
73
74
"Use this option with care, and only if you understand why it might sometimes\n"
" be necessary:\n"
"  --trim-centrics            Remove reflections which are duplicated in the\n"
"                              point group specified with the '-y' option.\n"
"\n"
75
76
77
78
"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"
79
80
"      --multiplicity         Multiply intensities by the number of\n"
"                              equivalent reflections.\n"
81
82
83
"\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
84
);
Thomas White's avatar
Thomas White committed
85
86
87
88
}


/* Apply Poisson noise to all reflections */
Thomas White's avatar
Thomas White committed
89
static void poisson_reflections(RefList *list, double adu_per_photon)
Thomas White's avatar
Thomas White committed
90
{
Thomas White's avatar
Thomas White committed
91
92
	Reflection *refl;
	RefListIterator *iter;
Thomas White's avatar
Thomas White committed
93

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

Thomas White's avatar
Thomas White committed
98
		double val, c;
Thomas White's avatar
Thomas White committed
99

Thomas White's avatar
Thomas White committed
100
		val = get_intensity(refl);
101

102
		c = adu_per_photon * poisson_noise(val/adu_per_photon);
103
		set_intensity(refl, c);
Thomas White's avatar
Thomas White committed
104
105

	}
106
107
108
109
}


/* Apply 10% uniform noise to all reflections */
Thomas White's avatar
Thomas White committed
110
static void noise_reflections(RefList *list)
111
{
Thomas White's avatar
Thomas White committed
112
113
	Reflection *refl;
	RefListIterator *iter;
114

Thomas White's avatar
Thomas White committed
115
116
117
	for ( refl = first_refl(list, &iter);
	      refl != NULL;
	      refl = next_refl(refl, iter) ) {
118

Thomas White's avatar
Thomas White committed
119
		double val, r;
120

Thomas White's avatar
Thomas White committed
121
		val = get_intensity(refl);
122
123
124
125

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

126
		set_intensity(refl, val);
127

Thomas White's avatar
Thomas White committed
128
	}
129
130
131
}


Thomas White's avatar
Thomas White committed
132
static RefList *template_reflections(RefList *list, RefList *template)
133
{
Thomas White's avatar
Thomas White committed
134
135
136
	Reflection *refl;
	RefListIterator *iter;
	RefList *out;
137

Thomas White's avatar
Thomas White committed
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
	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
163
                                 const SymOpList *holo, const SymOpList *mero)
Thomas White's avatar
Thomas White committed
164
165
166
167
{
	Reflection *refl;
	RefListIterator *iter;
	RefList *out;
Thomas White's avatar
Thomas White committed
168
	SymOpMask *m;
169

Thomas White's avatar
Thomas White committed
170
	/* FIXME: Check properly by coset decomposition */
Thomas White's avatar
Thomas White committed
171
	if ( num_equivs(holo, NULL) < num_equivs(mero, NULL) ) {
Thomas White's avatar
Thomas White committed
172
173
		ERROR("%s is not a subgroup of %s!\n", symmetry_name(mero),
		                                       symmetry_name(holo));
174
175
176
		return NULL;
	}

Thomas White's avatar
Thomas White committed
177
	out = reflist_new();
Thomas White's avatar
Thomas White committed
178
	m = new_symopmask(holo);
Thomas White's avatar
Thomas White committed
179
180
181
182

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

184
		double total, sigma;
185
186
187
188
		signed int h, k, l;
		int n, j;
		int skip;

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

191
		/* There is a many-to-one correspondence between reflections
Thomas White's avatar
Thomas White committed
192
193
194
		 * in the merohedral and holohedral groups.  Do the calculation
		 * only once for each reflection in the holohedral group, which
		 * contains fewer reflections.
195
		 */
Thomas White's avatar
Thomas White committed
196
		get_asymm(holo, h, k, l, &h, &k, &l);
Thomas White's avatar
Thomas White committed
197
		if ( find_refl(out, h, k, l) != NULL ) continue;
198

199
200
		total = 0.0;
		sigma = 0.0;
201
		skip = 0;
Thomas White's avatar
Thomas White committed
202
203
204
		special_position(holo, m, h, k, l);
		n = num_equivs(holo, m);

205
206
207
208
209
		for ( j=0; j<n; j++ ) {

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

Thomas White's avatar
Thomas White committed
210
			get_equiv(holo, m, j, h, k, l, &he, &ke, &le);
211
212
213
214
215
216

			/* 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
217
218
			if ( !find_equiv_in_list(in, he, ke, le, mero,
			                         &hu, &ku, &lu) ) {
219
220
221
222
223
				/* 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
224
225
				      h, k, l, he, ke, le, symmetry_name(mero),
				      h, k, l);
226
227
228
229
				skip = 1;
				break;
			}

Thomas White's avatar
Thomas White committed
230
231
			total += get_intensity(refl);
			sigma += pow(get_esd_intensity(refl), 2.0);
232
233
234
235
236

		}

		if ( !skip ) {

Thomas White's avatar
Thomas White committed
237
			Reflection *new = add_refl(out, h, k, l);
238
			set_intensity(new, total);
Thomas White's avatar
Thomas White committed
239
			set_esd_intensity(new, sqrt(sigma));
Thomas White's avatar
Thomas White committed
240
			set_redundancy(new, 1);
241
242
243
244
245

		}

	}

Thomas White's avatar
Thomas White committed
246
	return out;
247
248
249
}


Thomas White's avatar
Thomas White committed
250
251
static RefList *expand_reflections(RefList *in, const SymOpList *target,
                                                const SymOpList *initial)
Thomas White's avatar
Thomas White committed
252
{
Thomas White's avatar
Thomas White committed
253
254
255
	Reflection *refl;
	RefListIterator *iter;
	RefList *out;
Thomas White's avatar
Thomas White committed
256
	SymOpMask *m;
Thomas White's avatar
Thomas White committed
257

Thomas White's avatar
Thomas White committed
258
	/* FIXME: Check properly */
Thomas White's avatar
Thomas White committed
259
	if ( num_equivs(target, NULL) > num_equivs(initial, NULL) ) {
Thomas White's avatar
Thomas White committed
260
261
		ERROR("%s is not a subgroup of %s!\n", symmetry_name(initial),
		                                       symmetry_name(target));
262
263
264
		return NULL;
	}

Thomas White's avatar
Thomas White committed
265
	out = reflist_new();
Thomas White's avatar
Thomas White committed
266
	m = new_symopmask(initial);
Thomas White's avatar
Thomas White committed
267
268
269
270

	for ( refl = first_refl(in, &iter);
	      refl != NULL;
	      refl = next_refl(refl, iter) ) {
Thomas White's avatar
Thomas White committed
271
272
273
274

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

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

Thomas White's avatar
Thomas White committed
277
278
		special_position(initial, m, h, k, l);
		n = num_equivs(initial, m);
Thomas White's avatar
Thomas White committed
279

Thomas White's avatar
Thomas White committed
280
		/* For each equivalent in the higher symmetry group */
Thomas White's avatar
Thomas White committed
281
282
283
		for ( j=0; j<n; j++ ) {

			signed int he, ke, le;
Thomas White's avatar
Thomas White committed
284
			Reflection *new;
285
286
			int have_phase;
			double ph;
Thomas White's avatar
Thomas White committed
287
288

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

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

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

Thomas White's avatar
Thomas White committed
298
299
			/* FIXME: Make phase negative if the reflection is
			 * separated from the original via an inversion */
300
			ph = get_phase(refl, &have_phase);
301
			if ( have_phase ) set_phase(new, -ph);
Thomas White's avatar
Thomas White committed
302

Thomas White's avatar
Thomas White committed
303
304
305
306
		}

	}

Thomas White's avatar
Thomas White committed
307
308
	free_symopmask(m);

Thomas White's avatar
Thomas White committed
309
	return out;
Thomas White's avatar
Thomas White committed
310
311
312
}


313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
static RefList *trim_centrics(RefList *in, const SymOpList *sym)
{
	Reflection *refl;
	RefListIterator *iter;
	RefList *out;

	out = reflist_new();

	for ( refl = first_refl(in, &iter);
	      refl != NULL;
	      refl = next_refl(refl, iter) )
	{
		signed int h, k, l;
		signed int ha, ka, la;
		Reflection *new;

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

		/* Put it into the asymmetric unit */
		get_asymm(sym, h, k, l, &ha, &ka, &la);

		new = find_refl(out, ha, ka, la);
		if ( new != NULL ) {
			STATUS("Trimmed %i %i %i\n", h, k, l);
			continue;
		}

		/* Add new reflection under asymmetric (unique) indices */
		new = add_refl(out, ha, ka, la);
		copy_data(new, refl);
	}

	return out;
}


349
350
351
int main(int argc, char *argv[])
{
	int c;
352
353
	int config_noise = 0;
	int config_poisson = 0;
354
	int config_multi = 0;
355
	int config_trimc = 0;
Thomas White's avatar
Thomas White committed
356
357
358
	char *holo_str = NULL;
	char *mero_str = NULL;
	char *expand_str = NULL;
Thomas White's avatar
Thomas White committed
359
360
361
	SymOpList *holo = NULL;
	SymOpList *mero = NULL;
	SymOpList *expand = NULL;
Thomas White's avatar
Thomas White committed
362
363
	char *input_file = NULL;
	char *template = NULL;
Thomas White's avatar
Thomas White committed
364
	char *output = NULL;
365
	char *beamfile = NULL;
366
	struct beam_params *beam = NULL;
Thomas White's avatar
Thomas White committed
367
	RefList *input;
368
369
	double adu_per_photon = 0.0;
	int have_adu_per_photon = 0;
370
371
372
373

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
Thomas White's avatar
Thomas White committed
374
		{"template",           1, NULL,               't'},
375
376
		{"poisson",            0, &config_poisson,     1},
		{"noise",              0, &config_noise,       1},
Thomas White's avatar
Thomas White committed
377
		{"output",             1, NULL,               'o'},
378
		{"symmetry",           1, NULL,               'y'},
Thomas White's avatar
Thomas White committed
379
380
		{"twin",               1, NULL,               'w'},
		{"expand",             1, NULL,               'e'},
381
		{"intensities",        1, NULL,               'i'},
382
		{"multiplicity",       0, &config_multi,       1},
383
		{"trim-centrics",      0, &config_trimc,       1},
384
		{"adu-per-photon",     1, NULL,                2},
385
386
387
388
		{0, 0, NULL, 0}
	};

	/* Short options */
Thomas White's avatar
Thomas White committed
389
	while ((c = getopt_long(argc, argv, "ht:o:i:w:y:e:b:",
Thomas White's avatar
Thomas White committed
390
	                        longopts, NULL)) != -1) {
391
392

		switch (c) {
Thomas White's avatar
Thomas White committed
393
394

			case 'h' :
395
396
397
			show_help(argv[0]);
			return 0;

Thomas White's avatar
Thomas White committed
398
			case 't' :
Thomas White's avatar
Thomas White committed
399
400
401
			template = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
402
			case 'o' :
Thomas White's avatar
Thomas White committed
403
404
405
			output = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
406
			case 'i' :
Thomas White's avatar
Thomas White committed
407
			input_file = strdup(optarg);
408
409
			break;

Thomas White's avatar
Thomas White committed
410
			case 'y' :
Thomas White's avatar
Thomas White committed
411
			mero_str = strdup(optarg);
412
413
			break;

Thomas White's avatar
Thomas White committed
414
			case 'w' :
Thomas White's avatar
Thomas White committed
415
			holo_str = strdup(optarg);
416
417
			break;

Thomas White's avatar
Thomas White committed
418
			case 'e' :
Thomas White's avatar
Thomas White committed
419
			expand_str = strdup(optarg);
Thomas White's avatar
Thomas White committed
420
421
			break;

Thomas White's avatar
Thomas White committed
422
			case 2 :
423
424
			adu_per_photon = strtof(optarg, NULL);
			have_adu_per_photon = 1;
425
426
			break;

Thomas White's avatar
Thomas White committed
427
			case 0 :
428
429
			break;

Thomas White's avatar
Thomas White committed
430
			default :
431
432
			ERROR("Unhandled option '%c'\n", c);
			break;
Thomas White's avatar
Thomas White committed
433

434
435
436
437
		}

	}

Thomas White's avatar
Thomas White committed
438
	if ( (holo_str != NULL) && (expand_str != NULL) ) {
Thomas White's avatar
Thomas White committed
439
440
		ERROR("You cannot 'twin' and 'expand' at the same time.\n");
		ERROR("Decide which one you want to do first.\n");
441
442
443
444
445
446
447
448
449
450
		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
451
452
	}

Thomas White's avatar
Thomas White committed
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
	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
472
473
474
475
476
477
478
	if ( (expand != NULL) || (holo != NULL) || config_trimc
	  || config_multi ) {
		if ( mero == NULL ) {
			ERROR("You must specify the point group with -y.\n");
		}
	}

Thomas White's avatar
Thomas White committed
479
	input = read_reflections(input_file);
Thomas White's avatar
Thomas White committed
480
481
482
483
	if ( input == NULL ) {
		ERROR("Problem reading input file %s\n", input_file);
		return 1;
	}
Thomas White's avatar
Thomas White committed
484
	free(input_file);
485
486
487

	STATUS("%i reflections in input.\n", num_reflections(input));

Thomas White's avatar
Thomas White committed
488
489
490
	if ( (mero != NULL) && !config_trimc
	  && check_list_symmetry(input, mero) )
	{
491
		ERROR("The input reflection list does not appear to"
Thomas White's avatar
Thomas White committed
492
		      " have symmetry %s\n", symmetry_name(mero));
493
		return 1;
494
	}
495

496
	if ( config_poisson ) {
497
498
		if ( have_adu_per_photon ) {
			poisson_reflections(input, adu_per_photon);
499
		} else {
500
501
			ERROR("You must give the number of ADU per photon to "
			      "use --poisson.\n");
502
503
504
505
			return 1;
		}
	}

Thomas White's avatar
Thomas White committed
506
	if ( config_noise ) noise_reflections(input);
Thomas White's avatar
Thomas White committed
507

508
	if ( holo != NULL ) {
509

Thomas White's avatar
Thomas White committed
510
		RefList *new;
Thomas White's avatar
Thomas White committed
511
512
		STATUS("Twinning from %s into %s\n", symmetry_name(mero),
		                                     symmetry_name(holo));
Thomas White's avatar
Thomas White committed
513
514
515
516
517
518
519
520
521
		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;
522

Thomas White's avatar
Thomas White committed
523
524
	}

Thomas White's avatar
Thomas White committed
525
	if ( expand != NULL ) {
526

Thomas White's avatar
Thomas White committed
527
		RefList *new;
Thomas White's avatar
Thomas White committed
528
529
		STATUS("Expanding from %s into %s\n", symmetry_name(mero),
		                                      symmetry_name(expand));
Thomas White's avatar
Thomas White committed
530
531
532
533
534
		new = expand_reflections(input, expand, mero);

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

Thomas White's avatar
Thomas White committed
536
537
	}

538
539
540
541
542
543
544
545
546
547
548
549
550
	if ( config_trimc ) {

		RefList *new;

		STATUS("Trimming duplicate reflections in %s\n",
		       symmetry_name(mero));
		new = trim_centrics(input, mero);
		reflist_free(input);
		input = new;
		STATUS("%i output reflections\n", num_reflections(input));

	}

551
552
	if ( config_multi ) {

Thomas White's avatar
Thomas White committed
553
554
		Reflection *refl;
		RefListIterator *iter;
Thomas White's avatar
Thomas White committed
555
556
557
		SymOpMask *m;

		m = new_symopmask(mero);
558

Thomas White's avatar
Thomas White committed
559
560
561
		for ( refl = first_refl(input, &iter);
		      refl != NULL;
		      refl = next_refl(refl, iter) ) {
562
563

			double inty;
Thomas White's avatar
Thomas White committed
564
565
566
567
			signed int h, k, l;

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

Thomas White's avatar
Thomas White committed
569
570
			special_position(mero, m, h, k, l);
			inty *= (double)num_equivs(mero, m);
571
			set_intensity(refl, inty);
572
573

		}
Thomas White's avatar
Thomas White committed
574
575

		free_symopmask(m);
576
577
	}

Thomas White's avatar
Thomas White committed
578
	if ( template ) {
579

Thomas White's avatar
Thomas White committed
580
581
582
583
		RefList *t = read_reflections(template);
		RefList *new = template_reflections(input, t);
		reflist_free(input);
		input = new;
584

Thomas White's avatar
Thomas White committed
585
586
	}

587
	write_reflist(output, input);
Thomas White's avatar
Thomas White committed
588

Thomas White's avatar
Thomas White committed
589
	reflist_free(input);
590
591
592

	return 0;
}