process_hkl.c 8.37 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
29


Thomas White's avatar
Thomas White committed
30
31
32
33
/* Number of divisions for R vs |q| graphs */
#define RVQDV (20)


34
35
36
static void show_help(const char *s)
{
	printf("Syntax: %s [options]\n\n", s);
Thomas White's avatar
Thomas White committed
37
38
39
40
41
	printf(
"Assemble and process FEL Bragg intensities.\n"
"\n"
"  -h, --help              Display this help message.\n"
"  -i, --input=<filename>  Specify input filename (\"-\" for stdin).\n"
42
43
"  -o, --output=<filename> Specify output filename for merged intensities\n"
"                           (don't specify for no output).\n"
Thomas White's avatar
Thomas White committed
44
45
46
47
"\n"
"      --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"
Thomas White's avatar
Thomas White committed
48
"                           measurements.\n"
Thomas White's avatar
Thomas White committed
49
"  -e, --output-every=<n>  Analyse figures of merit after every n patterns\n"
50
51
"                           Default: 1000.  A value of zero means to do the\n"
"                           analysis only after reading all the patterns.\n"
Thomas White's avatar
Thomas White committed
52
53
"      --no-analyse        Don't perform any kind of analysis, just merge the\n"
"                           intensities.\n"
54
55
56
57
"      --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"
Thomas White's avatar
Thomas White committed
58
"  -r, --rvsq              Output lists of R vs |q| (\"Luzzatti plots\") when\n"
Thomas White's avatar
Thomas White committed
59
"                           analysing figures of merit.\n"
Thomas White's avatar
Thomas White committed
60
61
"      --stop-after=<n>    Stop after processing n patterns.  Zero means\n"
"                           keep going until the end of the input, and is the\n"
Thomas White's avatar
Thomas White committed
62
"                           default.\n"
Thomas White's avatar
Thomas White committed
63
64
"      --zone-axis         Output an [001] zone axis pattern each time the\n"
"                           figures of merit are analysed.\n");
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
}


static void write_RvsQ(const char *name, double *ref, double *trueref,
                       unsigned int *counts, double scale, UnitCell *cell)
{
	FILE *fh;
	double smax, sbracket;
	signed int h, k, l;

	fh = fopen(name, "w");

	smax = 0.0;
	for ( h=-INDMAX; h<INDMAX; h++ ) {
	for ( k=-INDMAX; k<INDMAX; k++ ) {
	for ( l=-INDMAX; l<INDMAX; l++ ) {
Thomas White's avatar
Thomas White committed
81
		double s = 2.0*resolution(cell, h, k, l);
82
83
84
85
86
87
88
		if ( (lookup_count(counts, h, k, l) > 0) && (s > smax) ) {
			smax = s;
		}
	}
	}
	}

Thomas White's avatar
Thomas White committed
89
	for ( sbracket=0.0; sbracket<smax; sbracket+=smax/RVQDV ) {
90
91
92

		double top = 0.0;
		double bot = 0.0;
Thomas White's avatar
Thomas White committed
93
94
95
96
		int nhits = 0;
		int nrefl = 0;
		double R;
		double hits_per_refl;
97
98
99
100
101
102
103

		for ( h=-INDMAX; h<INDMAX; h++ ) {
		for ( k=-INDMAX; k<INDMAX; k++ ) {
		for ( l=-INDMAX; l<INDMAX; l++ ) {

			double s;
			int c;
Thomas White's avatar
Thomas White committed
104
105
106

			if ( (h==0) && (k==0) && (l==0) ) continue;

107
			c = lookup_count(counts, h, k, l);
Thomas White's avatar
Thomas White committed
108
			s = 2.0*resolution(cell, h, k, l);
Thomas White's avatar
Thomas White committed
109
			if ((s>=sbracket) && (s<sbracket+smax/RVQDV) && (c>0)) {
110
111
112
113
114
115
116

				double obs, calc, obsi;

				obs = lookup_intensity(ref, h, k, l);
				calc = lookup_intensity(trueref, h, k, l);

				obsi = obs / (double)c;
Thomas White's avatar
Thomas White committed
117
118
119
120
121
				top += pow(obsi - scale*calc, 2.0);
				bot += pow(obsi, 2.0);
				nhits += c;
				nrefl++;

122
123
124
125
126
127
			}

		}
		}
		}

Thomas White's avatar
Thomas White committed
128
129
130
131
132
		R = sqrt(top/bot);
		hits_per_refl = nrefl ? (double)nhits/nrefl : 0;

		fprintf(fh, "%8.5f %8.5f %5.2f\n", sbracket+smax/(2.0*RVQDV),
		                                R, hits_per_refl);
133
134
135
136
137
138
139
140

	}
	fclose(fh);
}


static void process_reflections(double *ref, double *trueref,
                                unsigned int *counts, unsigned int n_patterns,
Thomas White's avatar
Thomas White committed
141
                                UnitCell *cell, int do_rvsq, int do_zoneaxis)
142
143
144
145
146
147
{
	int j;
	double mean_counts;
	int ctot = 0;
	int nmeas = 0;
	double R, scale;
148
	FILE *fh;
149
150
151
152
153
154
155
156

	for ( j=0; j<LIST_SIZE; j++ ) {
		ctot += counts[j];
		if ( counts[j] > 0 ) nmeas++;
	}
	mean_counts = (double)ctot/nmeas;

	R = stat_r2(ref, trueref, counts, LIST_SIZE, &scale);
Thomas White's avatar
Thomas White committed
157
	STATUS("%8u: R=%5.2f%% (sf=%7.4e)  mean meas/refl=%5.2f,"
158
159
	       " %i reflections measured\n",
	       n_patterns, R*100.0, scale, mean_counts, nmeas);
160

Thomas White's avatar
Thomas White committed
161
162
	if ( do_rvsq ) {
		/* Record graph of R against q for this N */
Thomas White's avatar
Thomas White committed
163
		char name[64];
164
		snprintf(name, 63, "R_vs_q-%u.dat", n_patterns);
Thomas White's avatar
Thomas White committed
165
166
		write_RvsQ(name, ref, trueref, counts, scale, cell);
	}
Thomas White's avatar
Thomas White committed
167
168
169

	if ( do_zoneaxis ) {
		char name[64];
170
		snprintf(name, 63, "ZA-%u.dat", n_patterns);
Thomas White's avatar
Thomas White committed
171
172
		write_reflections(name, counts, ref, 1, cell);
	}
173
174
175
176

	fh = fopen("results/convergence.dat", "a");
	fprintf(fh, "%u %5.2f %5.2f\n", n_patterns, R*100.0, mean_counts);
	fclose(fh);
177
178
179
180
181
182
183
}


int main(int argc, char *argv[])
{
	int c;
	char *filename = NULL;
184
	char *output = NULL;
185
186
	FILE *fh;
	unsigned int n_patterns;
Thomas White's avatar
Thomas White committed
187
	double *ref, *trueref = NULL;
188
189
	unsigned int *counts;
	char *rval;
Thomas White's avatar
Thomas White committed
190
	struct molecule *mol = NULL;
Thomas White's avatar
Thomas White committed
191
	int config_maxonly = 0;
Thomas White's avatar
Thomas White committed
192
	int config_every = 1000;
Thomas White's avatar
Thomas White committed
193
	int config_rvsq = 0;
Thomas White's avatar
Thomas White committed
194
	int config_stopafter = 0;
Thomas White's avatar
Thomas White committed
195
	int config_zoneaxis = 0;
Thomas White's avatar
Thomas White committed
196
	int config_noanalyse = 0;
197
	int config_sum = 0;
198
199
200
201
202

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
		{"input",              1, NULL,               'i'},
Thomas White's avatar
Thomas White committed
203
		{"output",             1, NULL,               'o'},
Thomas White's avatar
Thomas White committed
204
		{"max-only",           0, &config_maxonly,     1},
Thomas White's avatar
Thomas White committed
205
		{"output-every",       1, NULL,               'e'},
Thomas White's avatar
Thomas White committed
206
		{"no-analyse",         0, &config_noanalyse,   1},
Thomas White's avatar
Thomas White committed
207
		{"rvsq",               0, NULL,               'r'},
Thomas White's avatar
Thomas White committed
208
		{"stop-after",         1, NULL,               's'},
Thomas White's avatar
Thomas White committed
209
		{"zone-axis",          0, &config_zoneaxis,    1},
210
		{"sum",                0, &config_sum,    1},
211
212
213
214
		{0, 0, NULL, 0}
	};

	/* Short options */
215
	while ((c = getopt_long(argc, argv, "hi:e:ro:", longopts, NULL)) != -1) {
216
217
218
219
220
221
222
223
224
225
226
227

		switch (c) {
		case 'h' : {
			show_help(argv[0]);
			return 0;
		}

		case 'i' : {
			filename = strdup(optarg);
			break;
		}

228
229
230
231
232
		case 'o' : {
			output = strdup(optarg);
			break;
		}

Thomas White's avatar
Thomas White committed
233
234
235
236
237
		case 'r' : {
			config_rvsq = 1;
			break;
		}

Thomas White's avatar
Thomas White committed
238
239
240
241
242
		case 'e' : {
			config_every = atoi(optarg);
			break;
		}

Thomas White's avatar
Thomas White committed
243
244
245
246
247
		case 's' : {
			config_stopafter = atoi(optarg);
			break;
		}

248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
		case 0 : {
			break;
		}

		default : {
			return 1;
		}
		}

	}

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

	ref = new_list_intensity();
	counts = new_list_count();
Thomas White's avatar
Thomas White committed
266

267
	mol = load_molecule();
Thomas White's avatar
Thomas White committed
268
269
270
271
	if ( !config_noanalyse ) {
		get_reflections_cached(mol, eV_to_J(2.0e3));
		trueref = ideal_intensities(mol->reflections);
	}
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287

	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;
	}

	n_patterns = 0;
	do {

		char line[1024];
Thomas White's avatar
Tidy-up    
Thomas White committed
288
		signed int h, k, l, intensity;
289
290
291
292
		int r;

		rval = fgets(line, 1023, fh);
		if ( strncmp(line, "New pattern", 11) == 0 ) {
Thomas White's avatar
Tidy-up    
Thomas White committed
293

Thomas White's avatar
Thomas White committed
294
295
296
297
			if ( n_patterns == 0 ) {
				n_patterns++;
				continue;
			}
Thomas White's avatar
Tidy-up    
Thomas White committed
298

299
			if (config_every && (n_patterns % config_every == 0)) {
300
				process_reflections(ref, trueref, counts,
Thomas White's avatar
Thomas White committed
301
				                    n_patterns, mol->cell,
Thomas White's avatar
Thomas White committed
302
303
				                    config_rvsq,
				                    config_zoneaxis);
304
			}
Thomas White's avatar
Thomas White committed
305
306
307

			if ( n_patterns == config_stopafter ) break;

Thomas White's avatar
Thomas White committed
308
			n_patterns++;
309
310
311
312
313
		}

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

Thomas White's avatar
Tidy-up    
Thomas White committed
314
315
		if ( (h==0) && (k==0) && (l==0) ) continue;

Thomas White's avatar
Thomas White committed
316
317
		if ( !config_maxonly ) {
			integrate_intensity(ref, h, k, l, intensity);
318
319
320
321
322
			if ( config_sum ) {
				set_count(counts, h, k, l, 1);
			} else {
				integrate_count(counts, h, k, l, 1);
			}
Thomas White's avatar
Thomas White committed
323
324
325
326
327
328
		} else {
			if ( intensity > lookup_intensity(ref, h, k, l) ) {
				set_intensity(ref, h, k, l, intensity);
			}
			set_count(counts, h, k, l, 1);
		}
329
330
331
332
333

	} while ( rval != NULL );

	fclose(fh);

Thomas White's avatar
Thomas White committed
334
335
336
337
	if ( !config_noanalyse ) {
		process_reflections(ref, trueref, counts, n_patterns, mol->cell,
		                    config_rvsq, config_zoneaxis);
	}
338

339
	if ( output != NULL ) {
340
341
342
		UnitCell *cell = NULL;
		if ( mol != NULL ) cell = mol->cell;
		write_reflections(output, counts, ref, 0, cell);
343
	}
Thomas White's avatar
Thomas White committed
344

345
346
347
348
349
350
351
352
	if ( config_zoneaxis ) {
		char name[64];
		UnitCell *cell = NULL;
		if ( mol != NULL ) cell = mol->cell;
		snprintf(name, 63, "ZA-%u.dat", n_patterns);
		write_reflections(name, counts, ref, 1, cell);
	}

Thomas White's avatar
Thomas White committed
353
	STATUS("There were %u patterns.\n", n_patterns);
354
355
356

	return 0;
}