get_hkl.c 13.3 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
	if ( !is_subgroup(holo, mero) ) {
Thomas White's avatar
Thomas White committed
171
172
		ERROR("%s is not a subgroup of %s!\n", symmetry_name(mero),
		                                       symmetry_name(holo));
173
174
175
		return NULL;
	}

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

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

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

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

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

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

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

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

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

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

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

		}

		if ( !skip ) {

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

		}

	}

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


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

Thomas White's avatar
Thomas White committed
257
	if ( !is_subgroup(target, initial) ) {
Thomas White's avatar
Thomas White committed
258
259
		ERROR("%s is not a subgroup of %s!\n", symmetry_name(initial),
		                                       symmetry_name(target));
260
261
262
		return NULL;
	}

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

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

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

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

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

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

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

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

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

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

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

Thomas White's avatar
Thomas White committed
301
302
303
304
		}

	}

Thomas White's avatar
Thomas White committed
305
306
	free_symopmask(m);

Thomas White's avatar
Thomas White committed
307
	return out;
Thomas White's avatar
Thomas White committed
308
309
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
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;
}


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

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

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

		switch (c) {
Thomas White's avatar
Thomas White committed
391
392

			case 'h' :
393
394
395
			show_help(argv[0]);
			return 0;

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

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

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

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

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

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

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

Thomas White's avatar
Thomas White committed
425
			case 0 :
426
427
			break;

428
429
430
			case '?' :
			break;

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

435
436
437
438
		}

	}

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

Thomas White's avatar
Thomas White committed
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
	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
473
474
475
476
477
478
479
	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
480
	input = read_reflections(input_file);
Thomas White's avatar
Thomas White committed
481
482
483
484
	if ( input == NULL ) {
		ERROR("Problem reading input file %s\n", input_file);
		return 1;
	}
Thomas White's avatar
Thomas White committed
485
	free(input_file);
486
487
488

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

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

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

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

509
	if ( holo != NULL ) {
510

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

Thomas White's avatar
Thomas White committed
524
525
	}

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

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

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

Thomas White's avatar
Thomas White committed
537
538
	}

539
540
541
542
	if ( config_trimc ) {

		RefList *new;

Thomas White's avatar
Thomas White committed
543
544
545
546
547
548
		/* Can't do this if point group is invalid */
		if ( mero == NULL ) {
			ERROR("Need point group to trim centrics.\n");
			return 1;
		}

549
550
551
552
553
554
555
556
557
		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));

	}

558
559
	if ( config_multi ) {

Thomas White's avatar
Thomas White committed
560
561
		Reflection *refl;
		RefListIterator *iter;
Thomas White's avatar
Thomas White committed
562
563
564
		SymOpMask *m;

		m = new_symopmask(mero);
565

Thomas White's avatar
Thomas White committed
566
567
568
		for ( refl = first_refl(input, &iter);
		      refl != NULL;
		      refl = next_refl(refl, iter) ) {
569
570

			double inty;
Thomas White's avatar
Thomas White committed
571
572
573
574
			signed int h, k, l;

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

Thomas White's avatar
Thomas White committed
576
577
			special_position(mero, m, h, k, l);
			inty *= (double)num_equivs(mero, m);
578
			set_intensity(refl, inty);
579
580

		}
Thomas White's avatar
Thomas White committed
581
582

		free_symopmask(m);
583
584
	}

Thomas White's avatar
Thomas White committed
585
	if ( template ) {
586

Thomas White's avatar
Thomas White committed
587
588
589
590
		RefList *t = read_reflections(template);
		RefList *new = template_reflections(input, t);
		reflist_free(input);
		input = new;
591

Thomas White's avatar
Thomas White committed
592
593
	}

594
	write_reflist(output, input);
Thomas White's avatar
Thomas White committed
595

Thomas White's avatar
Thomas White committed
596
	reflist_free(input);
597
598
599

	return 0;
}