get_hkl.c 12.9 KB
Newer Older
1
2
3
/*
 * get_hkl.c
 *
4
 * Small program to manipulate reflection lists
5
 *
Thomas White's avatar
Thomas White committed
6
 * Copyright © 2012 Thomas White <taw@physics.org>
7
 *
Thomas White's avatar
Thomas White committed
8
9
10
11
12
13
14
15
16
17
18
19
20
21
 * 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/>.
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
 *
 */


#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
38
#include "reflist-utils.h"
39
#include "symmetry.h"
40
#include "beam-parameters.h"
41
42
43
44
45
46


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


/* Apply Poisson noise to all reflections */
Thomas White's avatar
Thomas White committed
85
static void poisson_reflections(RefList *list, double adu_per_photon)
Thomas White's avatar
Thomas White committed
86
{
Thomas White's avatar
Thomas White committed
87
88
	Reflection *refl;
	RefListIterator *iter;
Thomas White's avatar
Thomas White committed
89

Thomas White's avatar
Thomas White committed
90
91
92
	for ( refl = first_refl(list, &iter);
	      refl != NULL;
	      refl = next_refl(refl, iter) ) {
Thomas White's avatar
Thomas White committed
93

Thomas White's avatar
Thomas White committed
94
		double val, c;
Thomas White's avatar
Thomas White committed
95

Thomas White's avatar
Thomas White committed
96
		val = get_intensity(refl);
97

98
		c = adu_per_photon * poisson_noise(val/adu_per_photon);
99
		set_intensity(refl, c);
Thomas White's avatar
Thomas White committed
100
101

	}
102
103
104
105
}


/* Apply 10% uniform noise to all reflections */
Thomas White's avatar
Thomas White committed
106
static void noise_reflections(RefList *list)
107
{
Thomas White's avatar
Thomas White committed
108
109
	Reflection *refl;
	RefListIterator *iter;
110

Thomas White's avatar
Thomas White committed
111
112
113
	for ( refl = first_refl(list, &iter);
	      refl != NULL;
	      refl = next_refl(refl, iter) ) {
114

Thomas White's avatar
Thomas White committed
115
		double val, r;
116

Thomas White's avatar
Thomas White committed
117
		val = get_intensity(refl);
118
119
120
121

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

122
		set_intensity(refl, val);
123

Thomas White's avatar
Thomas White committed
124
	}
125
126
127
}


Thomas White's avatar
Thomas White committed
128
static RefList *template_reflections(RefList *list, RefList *template)
129
{
Thomas White's avatar
Thomas White committed
130
131
132
	Reflection *refl;
	RefListIterator *iter;
	RefList *out;
133

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

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

Thomas White's avatar
Thomas White committed
173
	out = reflist_new();
Thomas White's avatar
Thomas White committed
174
	m = new_symopmask(holo);
Thomas White's avatar
Thomas White committed
175
176
177
178

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

180
		double total, sigma;
181
182
183
184
		signed int h, k, l;
		int n, j;
		int skip;

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

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

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

201
202
203
204
205
		for ( j=0; j<n; j++ ) {

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

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

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

Thomas White's avatar
Thomas White committed
226
227
			total += get_intensity(refl);
			sigma += pow(get_esd_intensity(refl), 2.0);
228
229
230
231
232

		}

		if ( !skip ) {

Thomas White's avatar
Thomas White committed
233
			Reflection *new = add_refl(out, h, k, l);
234
			set_intensity(new, total);
Thomas White's avatar
Thomas White committed
235
			set_esd_intensity(new, sqrt(sigma));
236
237
238
239
240

		}

	}

Thomas White's avatar
Thomas White committed
241
	return out;
242
243
244
}


Thomas White's avatar
Thomas White committed
245
246
static RefList *expand_reflections(RefList *in, const SymOpList *target,
                                                const SymOpList *initial)
Thomas White's avatar
Thomas White committed
247
{
Thomas White's avatar
Thomas White committed
248
249
250
	Reflection *refl;
	RefListIterator *iter;
	RefList *out;
Thomas White's avatar
Thomas White committed
251
	SymOpMask *m;
Thomas White's avatar
Thomas White committed
252

Thomas White's avatar
Thomas White committed
253
	/* FIXME: Check properly */
Thomas White's avatar
Thomas White committed
254
	if ( num_equivs(target, NULL) > num_equivs(initial, NULL) ) {
Thomas White's avatar
Thomas White committed
255
256
		ERROR("%s is not a subgroup of %s!\n", symmetry_name(initial),
		                                       symmetry_name(target));
257
258
259
		return NULL;
	}

Thomas White's avatar
Thomas White committed
260
	out = reflist_new();
Thomas White's avatar
Thomas White committed
261
	m = new_symopmask(initial);
Thomas White's avatar
Thomas White committed
262
263
264
265

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

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

Thomas White's avatar
Thomas White committed
271
272
		get_indices(refl, &h, &k, &l);
		intensity = get_intensity(refl);
Thomas White's avatar
Thomas White committed
273

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

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

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

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

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

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

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

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

	}

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

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


346
347
348
int main(int argc, char *argv[])
{
	int c;
349
350
	int config_noise = 0;
	int config_poisson = 0;
351
	int config_multi = 0;
352
	int config_trimc = 0;
Thomas White's avatar
Thomas White committed
353
354
355
356
357
358
	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
359
360
	char *input_file = NULL;
	char *template = NULL;
Thomas White's avatar
Thomas White committed
361
	char *output = NULL;
362
	char *beamfile = NULL;
363
	struct beam_params *beam = NULL;
Thomas White's avatar
Thomas White committed
364
	RefList *input;
365
366
	double adu_per_photon = 0.0;
	int have_adu_per_photon = 0;
367
368
369
370

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

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

		switch (c) {
Thomas White's avatar
Thomas White committed
390
		case 'h' :
391
392
393
			show_help(argv[0]);
			return 0;

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

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

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

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

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

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

418
419
420
		case 2 :
			adu_per_photon = strtof(optarg, NULL);
			have_adu_per_photon = 1;
421
422
			break;

Thomas White's avatar
Thomas White committed
423
		case 0 :
424
425
			break;

Thomas White's avatar
Thomas White committed
426
		default :
427
428
429
430
431
			return 1;
		}

	}

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

Thomas White's avatar
Thomas White committed
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
	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
466
	input = read_reflections(input_file);
Thomas White's avatar
Thomas White committed
467
468
469
470
	if ( input == NULL ) {
		ERROR("Problem reading input file %s\n", input_file);
		return 1;
	}
Thomas White's avatar
Thomas White committed
471
	free(input_file);
472
473
474
475

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

	if ( !config_trimc && check_list_symmetry(input, mero) ) {
476
		ERROR("The input reflection list does not appear to"
Thomas White's avatar
Thomas White committed
477
		      " have symmetry %s\n", symmetry_name(mero));
478
		return 1;
479
	}
480

481
	if ( config_poisson ) {
482
483
		if ( have_adu_per_photon ) {
			poisson_reflections(input, adu_per_photon);
484
		} else {
485
486
			ERROR("You must give the number of ADU per photon to "
			      "use --poisson.\n");
487
488
489
490
			return 1;
		}
	}

Thomas White's avatar
Thomas White committed
491
	if ( config_noise ) noise_reflections(input);
Thomas White's avatar
Thomas White committed
492

493
	if ( holo != NULL ) {
494

Thomas White's avatar
Thomas White committed
495
		RefList *new;
Thomas White's avatar
Thomas White committed
496
497
		STATUS("Twinning from %s into %s\n", symmetry_name(mero),
		                                     symmetry_name(holo));
Thomas White's avatar
Thomas White committed
498
499
500
501
502
503
504
505
506
		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;
507

Thomas White's avatar
Thomas White committed
508
509
	}

Thomas White's avatar
Thomas White committed
510
	if ( expand != NULL ) {
511

Thomas White's avatar
Thomas White committed
512
		RefList *new;
Thomas White's avatar
Thomas White committed
513
514
		STATUS("Expanding from %s into %s\n", symmetry_name(mero),
		                                      symmetry_name(expand));
Thomas White's avatar
Thomas White committed
515
516
517
518
519
		new = expand_reflections(input, expand, mero);

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

Thomas White's avatar
Thomas White committed
521
522
	}

523
524
525
526
527
528
529
530
531
532
533
534
535
	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));

	}

536
537
	if ( config_multi ) {

Thomas White's avatar
Thomas White committed
538
539
		Reflection *refl;
		RefListIterator *iter;
Thomas White's avatar
Thomas White committed
540
541
542
		SymOpMask *m;

		m = new_symopmask(mero);
543

Thomas White's avatar
Thomas White committed
544
545
546
		for ( refl = first_refl(input, &iter);
		      refl != NULL;
		      refl = next_refl(refl, iter) ) {
547
548

			double inty;
Thomas White's avatar
Thomas White committed
549
550
551
552
			signed int h, k, l;

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

Thomas White's avatar
Thomas White committed
554
555
			special_position(mero, m, h, k, l);
			inty *= (double)num_equivs(mero, m);
556
			set_intensity(refl, inty);
557
558

		}
Thomas White's avatar
Thomas White committed
559
560

		free_symopmask(m);
561
562
	}

Thomas White's avatar
Thomas White committed
563
	if ( template ) {
564

Thomas White's avatar
Thomas White committed
565
566
567
568
		RefList *t = read_reflections(template);
		RefList *new = template_reflections(input, t);
		reflist_free(input);
		input = new;
569

Thomas White's avatar
Thomas White committed
570
571
	}

572
	write_reflist(output, input);
Thomas White's avatar
Thomas White committed
573

Thomas White's avatar
Thomas White committed
574
	reflist_free(input);
575
576
577

	return 0;
}