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
	printf(
"Assemble and process FEL Bragg intensities.\n"
"\n"
40
41
42
43
"  -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"
Thomas White's avatar
Thomas White committed
44
"\n"
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
"      --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"
"      --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"
"  -c, --compare-with=<file> Compare with reflection intensities in this file\n"
"\n"
"  -e, --output-every=<n>    Analyse figures of merit after every n patterns\n"
"                             Default: 1000.  A value of zero means to do the\n"
"                             analysis only after reading all the patterns.\n"
"  -r, --rvsq                Output lists of R vs |q| (\"Luzzatti plots\")\n"
"                             when analysing figures of merit.\n"
"      --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);
171
		write_reflections(name, counts, ref, 1, cell, 1);
Thomas White's avatar
Thomas White committed
172
	}
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;
190
	UnitCell *cell;
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;
196
	int config_sum = 0;
197
	char *intfile = NULL;
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
		{"rvsq",               0, NULL,               'r'},
Thomas White's avatar
Thomas White committed
207
		{"stop-after",         1, NULL,               's'},
Thomas White's avatar
Thomas White committed
208
		{"zone-axis",          0, &config_zoneaxis,    1},
209
		{"compare-with",       0, NULL,               'c'},
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
		case 'c' : {
			intfile = strdup(optarg);
			break;
		}

253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
		case 0 : {
			break;
		}

		default : {
			return 1;
		}
		}

	}

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

269
270
	if ( intfile != NULL ) {
		STATUS("Comparing against '%s'\n", intfile);
271
		trueref = read_reflections(intfile, NULL);
272
273
274
275
276
		free(intfile);
	} else {
		trueref = NULL;
	}

277
278
	ref = new_list_intensity();
	counts = new_list_count();
279
	cell = load_cell_from_pdb("molecule.pdb");
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295

	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];
296
297
		signed int h, k, l;
		float intensity;
298
299
300
		int r;

		rval = fgets(line, 1023, fh);
301
		if ( strncmp(line, "Reflections from indexing", 25) == 0 ) {
Thomas White's avatar
Tidy-up    
Thomas White committed
302

Thomas White's avatar
Thomas White committed
303
304
305
306
			if ( n_patterns == 0 ) {
				n_patterns++;
				continue;
			}
Thomas White's avatar
Tidy-up    
Thomas White committed
307

308
			if (config_every && (n_patterns % config_every == 0)) {
309
				process_reflections(ref, trueref, counts,
310
				                    n_patterns, cell,
Thomas White's avatar
Thomas White committed
311
312
				                    config_rvsq,
				                    config_zoneaxis);
313
			}
Thomas White's avatar
Thomas White committed
314
315
316

			if ( n_patterns == config_stopafter ) break;

Thomas White's avatar
Thomas White committed
317
			n_patterns++;
318
319
		}

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

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

Thomas White's avatar
Thomas White committed
325
326
		if ( !config_maxonly ) {
			integrate_intensity(ref, h, k, l, intensity);
327
328
329
330
331
			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
332
333
334
335
336
337
		} else {
			if ( intensity > lookup_intensity(ref, h, k, l) ) {
				set_intensity(ref, h, k, l, intensity);
			}
			set_count(counts, h, k, l, 1);
		}
338
339
340
341
342

	} while ( rval != NULL );

	fclose(fh);

343
	if ( trueref != NULL ) {
344
		process_reflections(ref, trueref, counts, n_patterns, cell,
Thomas White's avatar
Thomas White committed
345
346
		                    config_rvsq, config_zoneaxis);
	}
347

348
	if ( output != NULL ) {
349
		write_reflections(output, counts, ref, 0, cell, 1);
350
	}
Thomas White's avatar
Thomas White committed
351

352
353
354
	if ( config_zoneaxis ) {
		char name[64];
		snprintf(name, 63, "ZA-%u.dat", n_patterns);
355
		write_reflections(name, counts, ref, 1, cell, 10);
356
357
	}

Thomas White's avatar
Thomas White committed
358
	STATUS("There were %u patterns.\n", n_patterns);
359
360
361

	return 0;
}