process_hkl.c 20.9 KB
Newer Older
1
2
3
4
5
/*
 * process_hkl.c
 *
 * Assemble and process FEL Bragg intensities
 *
Thomas White's avatar
Thomas White committed
6
 * (c) 2006-2010 Thomas White <taw@physics.org>
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
 *
 * 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"
#include "statistics.h"
#include "sfac.h"
27
#include "reflections.h"
28
#include "symmetry.h"
29
#include "stream.h"
30
#include "beam-parameters.h"
31
32


33
/* Number of divisions for intensity histograms */
34
#define NBINS (50)
Thomas White's avatar
Thomas White committed
35
36


37
38
39
static void show_help(const char *s)
{
	printf("Syntax: %s [options]\n\n", s);
Thomas White's avatar
Thomas White committed
40
41
42
	printf(
"Assemble and process FEL Bragg intensities.\n"
"\n"
43
44
45
46
"  -h, --help                Display this help message.\n"
"  -i, --input=<filename>    Specify input filename (\"-\" for stdin).\n"
"  -o, --output=<filename>   Specify output filename for merged intensities\n"
"                             (don't specify for no output).\n"
47
"  -p, --pdb=<filename>      PDB file to use (default: molecule.pdb).\n"
48
"  -b, --beam=<file>         Get beam parameters from file (used for sigmas).\n"
Thomas White's avatar
Thomas White committed
49
"\n"
50
51
52
53
54
55
56
57
"      --max-only            Take the integrated intensity to be equal to the\n"
"                             maximum intensity measured for that reflection.\n"
"                             The default is to use the mean value from all\n"
"                             measurements.\n"
"      --sum                 Sum (rather than average) the intensities for the\n"
"                             final output list.  This is useful for comparing\n"
"                             results to radially summed powder patterns, but\n"
"                             will break R-factor analysis.\n"
58
"      --start-after=<n>     Skip n patterns at the start of the stream.\n"
59
60
61
"      --stop-after=<n>      Stop after processing n patterns.  Zero means\n"
"                             keep going until the end of the input, and is\n"
"                             the default.\n"
62
63
"  -g, --histogram=<h,k,l>   Calculate the histogram of measurements for this\n"
"                             reflection.\n"
64
"      --rmerge              Calculate and report Rmerge and Rpim\n"
65
"\n"
66
67
"      --scale               Scale each pattern for best fit with the current\n"
"                             model.\n"
Thomas White's avatar
Thomas White committed
68
"  -y, --symmetry=<sym>      Merge according to point group <sym>.\n"
69
"  -a, --input-symmetry=<a>  Specify the apparent (input) symmetry.\n"
70
71
"      --reference=<file>    Compare against intensities from <file> when\n"
"                             scaling or resolving ambiguities.\n"
72
73
"                             The symmetry of the reference list must be the\n"
"                             same as that given with '-y'.\n"
74
75
"      --outstream=<file>    Write an annotated version of the input stream\n"
"                             to <file>.\n"
76
);
77
78
79
}


80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
static void plot_histogram(double *vals, int n)
{
	int i;
	double max = -INFINITY;
	double min = +INFINITY;
	double step;
	int histo[NBINS];
	FILE *fh;

	fh = fopen("histogram.dat", "w");
	if ( fh == NULL ) {
		ERROR("Couldn't open 'histogram.dat'\n");
		return;
	}

	for ( i=0; i<n; i++ ) {
		if ( vals[i] > max ) max = vals[i];
		if ( vals[i] < min ) min = vals[i];
	}
	STATUS("%f %f\n", min, max);
100
	min--;  max++;
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121

	for ( i=0; i<NBINS; i++ ) {
		histo[i] = 0;
	}

	step = (max-min)/NBINS;

	for ( i=0; i<n; i++ ) {
		int bin;
		bin = (vals[i]-min)/step;
		histo[bin]++;
	}

	for ( i=0; i<NBINS; i++ ) {
		fprintf(fh, "%f %i\n", min+step*i, histo[i]);
	}

	fclose(fh);
}


122
123
124
125
126
/* Note "holo" needn't actually be a holohedral point group, if you want to try
 * something strange like resolving from a low-symmetry group into an even
 * lower symmetry one.
 */
static ReflItemList *get_twin_possibilities(const char *holo, const char *mero)
127
{
128
	ReflItemList *test_items;
129
	ReflItemList *twins;
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
	int np;

	np = num_general_equivs(holo) / num_general_equivs(mero);

	test_items = new_items();

	/* Some arbitrarily chosen reflections which can't be special
	 * reflections in any point group, i.e. lots of odd numbers,
	 * prime numbers and so on.  There's probably an analytical
	 * way of working these out, but this will do. */
	add_item(test_items, 1, 2, 3);
	add_item(test_items, 3, 7, 13);
	add_item(test_items, 5, 2, 1);

	twins = get_twins(test_items, holo, mero);
	delete_items(test_items);

Thomas White's avatar
Thomas White committed
147
148
	/* Idiot check.  Wouldn't be necessary if I could prove that the above
	 * set of arbitrarily chosen reflections were always general. */
149
150
151
152
153
154
155
156
157
	if ( num_items(twins) != np ) {
		ERROR("Whoops! Couldn't find all the twinning possiblities.\n");
		abort();
	}

	return twins;
}


Thomas White's avatar
Thomas White committed
158
159
160
static int resolve_twin(const double *model, ReflItemList *observed,
                        const double *patt, ReflItemList *items,
                        ReflItemList *twins, const char *holo, const char *mero)
161
{
162
163
164
165
166
167
168
169
170
171
172
173
174
	int n, i;
	double best_fom = 0.0;
	int best_op = 0;

	n = num_items(twins);

	for ( i=0; i<n; i++ ) {

		int j;
		int op;
		double *trial_ints = new_list_intensity();
		unsigned int *trial_counts = new_list_count();
		double fom;
Thomas White's avatar
Thomas White committed
175
		ReflItemList *intersection;
176
177
178
179
180
181
182
183
184
185

		op = get_item(twins, i)->op;

		for ( j=0; j<num_items(items); j++ ) {

			signed int h, k, l;
			struct refl_item *r = get_item(items, j);

			get_general_equiv(r->h, r->k, r->l, &h, &k, &l,
			                  holo, op);
186
			get_asymm(h, k, l, &h, &k, &l, mero);
187
188
189
190
191
192
193

			set_intensity(trial_ints, h, k, l,
			              lookup_intensity(patt, r->h, r->k, r->l));
			set_count(trial_counts, h, k, l, 1);

		}

Thomas White's avatar
Thomas White committed
194
		intersection = intersection_items(observed, items);
Thomas White's avatar
Thomas White committed
195
		fom = stat_pearson_i(trial_ints, model, intersection);
Thomas White's avatar
Thomas White committed
196
		delete_items(intersection);
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213

		free(trial_ints);
		free(trial_counts);

		//printf(" %f", fom);
		if ( fom > best_fom ) {
			best_fom = fom;
			best_op = op;
		}

	}
	//printf("\n");

	return best_op;
}


Thomas White's avatar
Thomas White committed
214
215
216
static void merge_pattern(double *model, ReflItemList *observed,
                          const double *new,  ReflItemList *items,
                          unsigned int *model_counts,  int mo,
217
                          ReflItemList *twins,
218
219
220
                          const char *holo, const char *mero,
                          ReflItemList *reference, const double *reference_i,
                          double *hist_vals,
221
                          signed int hist_h, signed int hist_k,
222
                          signed int hist_l, int *hist_n, double *devs,
223
                          double *means, FILE *outfh)
224
{
225
	int i;
226
	int twin;
227
	ReflItemList *sym_items = new_items();
228

229
	if ( twins != NULL ) {
230
231
232
233
234
235
236
		if ( reference != NULL ) {
			twin = resolve_twin(reference_i, reference, new, items,
			                    twins, holo, mero);
		} else {
			twin = resolve_twin(model, observed, new, items,
			                    twins, holo, mero);
		}
237
238
239
	} else {
		twin = 0;
	}
240

241
242
243
244
	if ( outfh != NULL ) {
		fprintf(outfh, "twin=%i\n", twin);
	}

245
	for ( i=0; i<num_items(items); i++ ) {
246
247

		double intensity;
248
		signed int hs, ks, ls;
249
250
		signed int h, k, l;
		struct refl_item *item;
251

252
253
		item = get_item(items, i);

254
255
256
257
		hs = item->h;
		ks = item->k;
		ls = item->l;

258
259
		/* Transform into correct side of the twin law.
		 * "twin" is always zero if no de-twinning is performed. */
260
		get_general_equiv(hs, ks, ls, &h, &k, &l, holo, twin);
261

Thomas White's avatar
Thomas White committed
262
		/* Put into the asymmetric unit for the target group */
263
		get_asymm(h, k, l, &h, &k, &l, mero);
264

265
266
267
		/* Read the intensity from the original location
		 * (i.e. before screwing around with symmetry) */
		intensity = lookup_intensity(new, hs, ks, ls);
268

Thomas White's avatar
Thomas White committed
269
		/* User asked for max only? */
270
271
272
273
274
275
276
277
		if ( !mo ) {
			integrate_intensity(model, h, k, l, intensity);
		} else {
			if ( intensity > lookup_intensity(model, h, k, l) ) {
				set_intensity(model, h, k, l, intensity);
			}
		}

278
		if ( devs != NULL ) {
279
280
			double m;
			m = lookup_intensity(means, h, k, l);
281
282
			integrate_intensity(devs, h, k, l,
			                    pow(intensity-m, 2.0));
283
284
		}

285
		if ( !find_item(sym_items, h, k, l) ) {
Thomas White's avatar
Thomas White committed
286
287
			add_item(sym_items, h, k, l);
		}
288

289
		/* Increase count */
Thomas White's avatar
Thomas White committed
290
291
		integrate_count(model_counts, h, k, l, 1);

292
293
294
295
296
297
298
299
		if ( hist_vals != NULL ) {
			int p = *hist_n;
			if ( (h==hist_h) && (k==hist_k) && (l==hist_l) ) {
				hist_vals[p] = intensity;
				*hist_n = p+1;
			}
		}

300
	}
301
302
303
304
305

	/* Dump the reflections in this pattern into the overall list */
	union_items(observed, sym_items);

	delete_items(sym_items);
306
307
308
}


309
310
311
enum {
	SCALE_NONE,
	SCALE_CONSTINT,
312
313
	SCALE_INTPERBRAGG,
	SCALE_TWOPASS,
314
315
316
};


317
318
static void scale_intensities(const double *model, ReflItemList *model_items,
                              double *new_pattern, ReflItemList *new_items,
319
                              double f0, int f0_valid, const char *sym)
320
321
{
	double s;
322
323
	double top = 0.0;
	double bot = 0.0;
324
	unsigned int i;
325
	const int scaling = SCALE_INTPERBRAGG;
326

327
328
329
330
331
332
333
334
335
	for ( i=0; i<num_items(new_items); i++ ) {

		double i1, i2;
		struct refl_item *it;
		signed int hu, ku, lu;

		/* Get the next item in the list of new reflections */
		it = get_item(new_items, i);

336
337
338
339
340
341
342
343
344
345
346
347
348
		switch ( scaling ) {
		case SCALE_TWOPASS :

			/* Find the (only) partner in the model */
			find_unique_equiv(model_items, it->h, it->k, it->l, sym,
				          &hu, &ku, &lu);

			i1 = lookup_intensity(model, hu, ku, lu);
			i2 = lookup_intensity(new_pattern, it->h, it->k, it->l);

			/* Calculate LSQ estimate of scaling factor */
			top += i1 * i2;
			bot += i2 * i2;
349

350
			break;
351

352
353
354
355
356
357
358
359
		case SCALE_CONSTINT :

			/* Sum up the intensity in the pattern */
			i2 = lookup_intensity(new_pattern, it->h, it->k, it->l);
			top += i2;

			break;

360
361
362
363
364
365
366
367
368
		case SCALE_INTPERBRAGG :

			/* Sum up the intensity in the pattern */
			i2 = lookup_intensity(new_pattern, it->h, it->k, it->l);
			top += i2;
			bot += 1.0;

			break;

369
		}
370
371
372

	}

373
374
375
376
377
378
379
	switch ( scaling ) {
	case SCALE_TWOPASS :
		s = top / bot;
		break;
	case SCALE_CONSTINT :
		s = 1000.0 / top;
		break;
380
381
382
	case SCALE_INTPERBRAGG :
		s = 1000.0 / (top/bot);
		break;
383
384
	}
	//if ( f0_valid ) printf("%f %f\n", s, f0);
385
386
387
388
389
390
391
392

	/* Multiply the new pattern up by "s" */
	for ( i=0; i<LIST_SIZE; i++ ) {
		new_pattern[i] *= s;
	}
}


Thomas White's avatar
Thomas White committed
393
394
static void merge_all(FILE *fh, double **pmodel, ReflItemList **pobserved,
                      unsigned int **pcounts,
395
                      int config_maxonly, int config_scale, int config_sum,
396
                      int config_startafter, int config_stopafter,
397
                      ReflItemList *twins, const char *holo, const char *mero,
398
399
400
                      int n_total_patterns,
                      ReflItemList *reference, double *reference_i,
                      double *hist_vals,
401
                      signed int hist_h, signed int hist_k, signed int hist_l,
402
                      int *hist_i, double *devs, double *means, FILE *outfh)
403
404
405
406
407
{
	char *rval;
	float f0;
	int n_nof0 = 0;
	int f0_valid = 0;
Thomas White's avatar
Thomas White committed
408
	int n_patterns = 0;
409
410
	double *new_pattern = new_list_intensity();
	ReflItemList *items = new_items();
Thomas White's avatar
Thomas White committed
411
412
413
414
	ReflItemList *observed = new_items();
	double *model = new_list_intensity();
	unsigned int *counts = new_list_count();
	int i;
415

416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
	if ( config_startafter != 0 ) {

		do {

			char line[1024];

			rval = fgets(line, 1023, fh);
			if ( (strncmp(line, "Reflections from indexing", 25)
			                                                   == 0)
			    || (strncmp(line, "New pattern", 11) == 0) ) {
				n_patterns++;
			}

			if ( n_patterns == config_startafter ) break;

		} while ( rval != NULL );

	}

435
436
437
438
439
440
441
442
	do {

		char line[1024];
		signed int h, k, l;
		float intensity;
		int r;

		rval = fgets(line, 1023, fh);
Thomas White's avatar
Thomas White committed
443
444
		if ( ((strncmp(line, "Reflections from indexing", 25) == 0)
		  || (rval == NULL)) && (num_items(items)>0) ) {
445
446
447
448
449
450
451
452
453

			/* Assume a default I0 if we don't have one by now */
			if ( config_scale && !f0_valid ) {
				n_nof0++;
				f0 = 1.0;
			}

			/* Scale if requested */
			if ( config_scale ) {
454
455
456
				if ( reference == NULL ) {
					scale_intensities(model, observed,
					                  new_pattern, items,
457
					                  f0, f0_valid, mero);
458
459
460
461
				} else {
					scale_intensities(reference_i,
					                  reference,
					                  new_pattern, items,
462
					                  f0, f0_valid, mero);
463
				}
464
465
466
			}

			/* Start of second or later pattern */
Thomas White's avatar
Thomas White committed
467
468
			merge_pattern(model, observed, new_pattern, items,
			              counts, config_maxonly,
469
			              twins, holo, mero,
470
			              reference, reference_i,
471
			              hist_vals, hist_h, hist_k, hist_l,
472
			              hist_i, devs, means, outfh);
473

Thomas White's avatar
Thomas White committed
474
			n_patterns++;
475
476
477
478
479
			if ( n_patterns == config_stopafter ) {
				progress_bar(n_total_patterns, n_total_patterns,
				             "Merging");
				break;
			}
Thomas White's avatar
Thomas White committed
480
			progress_bar(n_patterns, n_total_patterns, "Merging");
481

482
			/* Reset for the next pattern */
483
484
485
486
487
488
			clear_items(items);

			f0_valid = 0;

		}

489
490
491
492
		if ( outfh != NULL ) {
			fprintf(outfh, "%s", line);
		}

493
494
495
496
497
498
499
500
501
502
503
504
505
		if ( strncmp(line, "f0 = ", 5) == 0 ) {
			r = sscanf(line, "f0 = %f", &f0);
			if ( r != 1 ) {
				f0 = 1.0;
				f0_valid = 0;
				continue;
			}
			f0_valid = 1;
		}

		r = sscanf(line, "%i %i %i %f", &h, &k, &l, &intensity);
		if ( r != 4 ) continue;

506
		/* Not interested in the central beam */
507
508
		if ( (h==0) && (k==0) && (l==0) ) continue;

509
510
		/* The same raw indices (before mapping into the asymmetric
		 * unit should not turn up twice in one pattern. */
511
512
513
514
515
		if ( find_item(items, h, k, l) != 0 ) {
			ERROR("More than one measurement for %i %i %i in"
			      " pattern number %i\n", h, k, l, n_patterns);
		}
		set_intensity(new_pattern, h, k, l, intensity);
516
517
518

		/* NB: This list contains raw indices, before working out
		 * where they belong in the asymmetric unit. */
519
520
521
522
523
524
525
		add_item(items, h, k, l);

	} while ( rval != NULL );

	delete_items(items);
	free(new_pattern);

526
527
	/* Calculate mean intensity if necessary */
	if ( !config_sum && !config_maxonly ) {
Thomas White's avatar
Thomas White committed
528
529
530
531
532
533
534
535
536
537
538
		for ( i=0; i<IDIM*IDIM*IDIM; i++ ) {
			if ( counts[i] > 0 ) {
				model[i] /= (double)counts[i];
			}
		}
	}

	*pmodel = model;
	*pcounts = counts;
	*pobserved = observed;

539
540
541
542
	STATUS("%i patterns had no f0 valid value.\n", n_nof0);
}


543
544
545
546
int main(int argc, char *argv[])
{
	int c;
	char *filename = NULL;
547
	char *output = NULL;
548
	FILE *fh;
549
	double *model;
Thomas White's avatar
Thomas White committed
550
	unsigned int *counts;
551
	UnitCell *cell = NULL;
Thomas White's avatar
Thomas White committed
552
	int config_maxonly = 0;
553
	int config_startafter = 0;
Thomas White's avatar
Thomas White committed
554
	int config_stopafter = 0;
555
	int config_sum = 0;
556
	int config_scale = 0;
557
	int config_rmerge = 0;
558
	unsigned int n_total_patterns;
559
	char *sym = NULL;
560
	char *in_sym = NULL;
561
	char *pdb = NULL;
562
	ReflItemList *twins;
Thomas White's avatar
Thomas White committed
563
	ReflItemList *observed;
564
	int i;
565
566
567
568
	char *histo = NULL;
	signed int hist_h, hist_k, hist_l;
	double *hist_vals = NULL;
	int hist_i;
569
570
571
572
573
	char *outstream = NULL;
	char *reference = NULL;
	ReflItemList *reference_items;
	double *reference_i;
	FILE *outfh = NULL;
574
	struct beam_params *beam = NULL;
575
	int space_for_hist = 0;
576
577
	double *devs;
	double *esds;
578
579
580
581
582

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
		{"input",              1, NULL,               'i'},
Thomas White's avatar
Thomas White committed
583
		{"output",             1, NULL,               'o'},
Thomas White's avatar
Thomas White committed
584
		{"max-only",           0, &config_maxonly,     1},
Thomas White's avatar
Thomas White committed
585
		{"output-every",       1, NULL,               'e'},
Thomas White's avatar
Thomas White committed
586
		{"stop-after",         1, NULL,               's'},
587
		{"start-after",        1, NULL,               'f'},
588
		{"sum",                0, &config_sum,         1},
589
		{"scale",              0, &config_scale,       1},
590
		{"symmetry",           1, NULL,               'y'},
591
		{"input-symmetry",     1, NULL,               'a'},
592
		{"pdb",                1, NULL,               'p'},
593
		{"histogram",          1, NULL,               'g'},
594
		{"rmerge",             0, &config_rmerge,      1},
595
		{"outstream",          1, NULL,                2},
596
		{"reference",          1, NULL,               'r'},
597
		{"beam",               1, NULL,               'b'},
598
599
600
601
		{0, 0, NULL, 0}
	};

	/* Short options */
Thomas White's avatar
Thomas White committed
602
	while ((c = getopt_long(argc, argv, "hi:e:r:o:p:y:g:f:a:b:",
603
	                        longopts, NULL)) != -1) {
604
605

		switch (c) {
Thomas White's avatar
Thomas White committed
606
		case 'h' :
607
608
609
			show_help(argv[0]);
			return 0;

Thomas White's avatar
Thomas White committed
610
		case 'i' :
611
612
613
			filename = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
614
		case 'o' :
615
616
617
			output = strdup(optarg);
			break;

Thomas White's avatar
Thomas White committed
618
		case 's' :
Thomas White's avatar
Thomas White committed
619
620
621
			config_stopafter = atoi(optarg);
			break;

622
623
624
625
		case 'f' :
			config_startafter = atoi(optarg);
			break;

Thomas White's avatar
Thomas White committed
626
		case 'p' :
627
628
629
			pdb = strdup(optarg);
			break;

630
631
632
633
		case 'y' :
			sym = strdup(optarg);
			break;

634
635
636
637
		case 'g' :
			histo = strdup(optarg);
			break;

638
639
640
641
		case 'r' :
			reference = strdup(optarg);
			break;

642
		case 2 :
643
644
645
			outstream = strdup(optarg);
			break;

646
647
648
649
		case 'a' :
			in_sym = strdup(optarg);
			break;

650
651
652
653
654
655
656
657
658
		case 'b' :
			beam = get_beam_parameters(optarg);
			if ( beam == NULL ) {
				ERROR("Failed to load beam parameters"
				      " from '%s'\n", optarg);
				return 1;
			}
			break;

Thomas White's avatar
Thomas White committed
659
		case 0 :
660
661
			break;

Thomas White's avatar
Thomas White committed
662
		default :
663
664
665
666
667
			return 1;
		}

	}

668
669
670
671
672
673
	if ( config_sum && config_rmerge ) {
		ERROR("Options --sum and --rmerge do not make sense"
		      " together.\n");
		return 1;
	}

674
675
676
677
678
	if ( filename == NULL ) {
		ERROR("Please specify filename using the -i option\n");
		return 1;
	}

679
680
681
682
	if ( pdb == NULL ) {
		pdb = strdup("molecule.pdb");
	}

Thomas White's avatar
Thomas White committed
683
684
685
	cell = load_cell_from_pdb(pdb);
	free(pdb);

686
687
	if ( sym == NULL ) sym = strdup("1");

688
	/* Show useful symmetry information */
689
690
691
	if ( (sym != NULL) && (in_sym != NULL) ) {

		int np = num_general_equivs(in_sym) / num_general_equivs(sym);
692
693
694
695
		if ( np > 1 ) {

			STATUS("Resolving point group %s into %s "
			       "(%i possibilities)\n",
696
			       in_sym, sym, np);
697
			/* Get the list of twin/Bijvoet possibilities */
698
699
			twins = get_twin_possibilities(in_sym, sym);
			STATUS("Twin operation indices from %s are:", in_sym);
700
701
702
703
			for ( i=0; i<num_items(twins); i++ ) {
				STATUS(" %i", get_item(twins, i)->op);
			}
			STATUS("\n");
704

705
		} else {
706
707
			STATUS("No resolution required to get from %s to %s\n",
			       in_sym, sym);
708
709
			twins = NULL;
		}
710

711
	} else {
712
		STATUS("No twin resolution requested (use -a otherwise).\n");
713
		twins = NULL;
714
		in_sym = strdup(sym);
715
	}
716

717
718
719
720
721
722
723
724
725
726
727
	/* Open the data stream */
	if ( strcmp(filename, "-") == 0 ) {
		fh = stdin;
	} else {
		fh = fopen(filename, "r");
	}
	free(filename);
	if ( fh == NULL ) {
		ERROR("Failed to open input file\n");
		return 1;
	}
Thomas White's avatar
Tidy-up    
Thomas White committed
728

729
730
731
732
	/* Read the reference reflections */
	if ( reference != NULL ) {
		reference_i = new_list_intensity();
		reference_items = read_reflections(reference, reference_i,
733
		                                   NULL, NULL, NULL);
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
		if ( reference_items == NULL ) {
			ERROR("Couldn't read '%s'\n", reference);
			return 1;
		}
	} else {
		reference_items = NULL;
		reference_i = NULL;
	}

	if ( outstream != NULL ) {
		outfh = fopen(outstream, "w");
		if ( outfh == NULL ) {
			ERROR("Couldn't open '%s'\n", outstream);
			return 1;
		}
	}

751
752
753
754
	/* Count the number of patterns in the file */
	n_total_patterns = count_patterns(fh);
	STATUS("There are %i patterns to process\n", n_total_patterns);
	rewind(fh);
755

756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
	if ( histo != NULL ) {
		int r;
		r = sscanf(histo, "%i,%i,%i", &hist_h, &hist_k, &hist_l);
		if ( r != 3 ) {
			ERROR("Invalid indices for '--histogram'\n");
			return 1;
		}
		space_for_hist = n_total_patterns * num_general_equivs(sym);
		hist_vals = malloc(space_for_hist * sizeof(double));
		free(histo);
		STATUS("Histogramming %i %i %i -> ", hist_h, hist_k, hist_l);
		/* Put into the asymmetric cell for the target group */
		get_asymm(hist_h, hist_k, hist_l,
		          &hist_h, &hist_k, &hist_l, sym);
		STATUS("%i %i %i\n", hist_h, hist_k, hist_l);
	}

773
	hist_i = 0;
Thomas White's avatar
Thomas White committed
774
	merge_all(fh, &model, &observed, &counts,
775
776
	          config_maxonly, config_scale, config_sum,
	          config_startafter, config_stopafter,
777
                  twins, in_sym, sym, n_total_patterns,
778
                  reference_items, reference_i,
779
                  hist_vals, hist_h, hist_k, hist_l, &hist_i, NULL, NULL,
780
                  outfh);
781
	rewind(fh);
782
783
784
	if ( space_for_hist && (hist_i >= space_for_hist) ) {
		ERROR("Histogram array was too small!\n");
	}
785

786
787
788
789
790
791
	if ( hist_vals != NULL ) {
		STATUS("%i %i %i was seen %i times.\n", hist_h, hist_k, hist_l,
		                                        hist_i);
		plot_histogram(hist_vals, hist_i);
	}

792
793
	STATUS("Extra pass to calculate ESDs...\n");
	devs = new_list_intensity();
794
	esds = new_list_sigma();
795
796
797
798
799
800
801
	rewind(fh);
	merge_all(fh, &model, &observed, &counts,
	          config_maxonly, config_scale, 0,
	          config_startafter, config_stopafter,
	          twins, in_sym, sym,
	          n_total_patterns, reference_items, reference_i,
	          NULL, 0, 0, 0, NULL, devs, model, NULL);
802

803
	for ( i=0; i<num_items(observed); i++ ) {
804

805
806
807
808
		struct refl_item *it;
		signed int h, k, l;
		double dev, esd;
		unsigned int count;
809

810
811
		it = get_item(observed, i);
		h = it->h;  k = it->k,  l = it->l;
812

813
814
		dev = lookup_intensity(devs, h, k, l);
		count = lookup_count(counts, h, k, l);
815

Thomas White's avatar
Thomas White committed
816
817
		if ( count < 2 ) {
			/* If we have only one measurement, the error is 100% */
818
			esd = lookup_sigma(model, h, k, l);
Thomas White's avatar
Thomas White committed
819
820
821
		} else {
			esd = sqrt(dev) / (double)count;
		}
822
		set_sigma(esds, h, k, l, esd);
823

824
	}
825

826
	if ( output != NULL ) {
827

828
		double adu_per_photon;
829

830
831
832
833
834
835
		if ( beam == NULL ) {
			adu_per_photon = 167.0;
			STATUS("No beam parameters file provided (use -b), "
			       "so I'm assuming 167.0 ADU per photon.\n");
		} else {
			adu_per_photon = beam->adu_per_photon;
836
837
		}

838
839
		write_reflections(output, observed, model, esds,
		                  NULL, counts, cell);
840
841
842
843

	}

	fclose(fh);
844
	free(sym);
Thomas White's avatar
Thomas White committed
845
	free(model);
Thomas White's avatar
Thomas White committed
846
	free(counts);
Thomas White's avatar
Thomas White committed
847
	free(output);
848
	if ( cell != NULL ) cell_free(cell);
849

850
851
	return 0;
}