process_hkl.c 6.99 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
/*
 * process_hkl.c
 *
 * Assemble and process FEL Bragg intensities
 *
 * (c) 2006-2009 Thomas White <thomas.white@desy.de>
 *
 * 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
42
43
44
45
	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"
"\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
46
"                           measurements.\n"
Thomas White's avatar
Thomas White committed
47
"  -e, --output-every=<n>  Analyse figures of merit after every n patterns\n"
48
49
"                           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
50
"  -r, --rvsq              Output lists of R vs |q| (\"Luzzatti plots\") when\n"
Thomas White's avatar
Thomas White committed
51
"                           analysing figures of merit.\n"
Thomas White's avatar
Thomas White committed
52
53
54
"      --stop-after=<n>    Stop after processing n patterns.  Zero means\n"
"                           keep going until the end of the input, and is the\n"
"                           default."
Thomas White's avatar
Thomas White committed
55
56
"      --zone-axis         Output an [001] zone axis pattern each time the\n"
"                           figures of merit are analysed.\n");
57
58
59
60
61
62
63
64
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++ ) {
		double s = resolution(cell, h, k, l);
		if ( (lookup_count(counts, h, k, l) > 0) && (s > smax) ) {
			smax = s;
		}
	}
	}
	}

Thomas White's avatar
Thomas White committed
81
	for ( sbracket=0.0; sbracket<smax; sbracket+=smax/RVQDV ) {
82
83
84

		double top = 0.0;
		double bot = 0.0;
Thomas White's avatar
Thomas White committed
85
86
87
88
		int nhits = 0;
		int nrefl = 0;
		double R;
		double hits_per_refl;
89
90
91
92
93
94
95

		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
96
97
98

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

99
100
			c = lookup_count(counts, h, k, l);
			s = resolution(cell, h, k, l);
Thomas White's avatar
Thomas White committed
101
			if ((s>=sbracket) && (s<sbracket+smax/RVQDV) && (c>0)) {
102
103
104
105
106
107
108

				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
109
110
111
112
113
				top += pow(obsi - scale*calc, 2.0);
				bot += pow(obsi, 2.0);
				nhits += c;
				nrefl++;

114
115
116
117
118
119
			}

		}
		}
		}

Thomas White's avatar
Thomas White committed
120
121
122
123
124
		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);
125
126
127
128
129
130
131
132

	}
	fclose(fh);
}


static void process_reflections(double *ref, double *trueref,
                                unsigned int *counts, unsigned int n_patterns,
Thomas White's avatar
Thomas White committed
133
                                UnitCell *cell, int do_rvsq, int do_zoneaxis)
134
135
136
137
138
139
140
141
142
143
144
145
146
147
{
	int j;
	double mean_counts;
	int ctot = 0;
	int nmeas = 0;
	double R, scale;

	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
148
	STATUS("%8u: R=%5.2f%% (sf=%7.4e)  mean meas/refl=%5.2f,"
149
150
	       " %i reflections measured\n",
	       n_patterns, R*100.0, scale, mean_counts, nmeas);
151

Thomas White's avatar
Thomas White committed
152
153
	if ( do_rvsq ) {
		/* Record graph of R against q for this N */
Thomas White's avatar
Thomas White committed
154
		char name[64];
Thomas White's avatar
Thomas White committed
155
156
157
		snprintf(name, 63, "results/R_vs_q-%u.dat", n_patterns);
		write_RvsQ(name, ref, trueref, counts, scale, cell);
	}
Thomas White's avatar
Thomas White committed
158
159
160
161
162
163

	if ( do_zoneaxis ) {
		char name[64];
		snprintf(name, 63, "results/ZA-%u.dat", n_patterns);
		write_reflections(name, counts, ref, 1, cell);
	}
164
165
166
167
168
169
170
171
172
173
174
175
176
}


int main(int argc, char *argv[])
{
	int c;
	char *filename = NULL;
	FILE *fh;
	unsigned int n_patterns;
	double *ref, *trueref;
	unsigned int *counts;
	char *rval;
	struct molecule *mol;
Thomas White's avatar
Thomas White committed
177
	int config_maxonly = 0;
Thomas White's avatar
Thomas White committed
178
	int config_every = 1000;
Thomas White's avatar
Thomas White committed
179
	int config_rvsq = 0;
Thomas White's avatar
Thomas White committed
180
	int config_stopafter = 0;
Thomas White's avatar
Thomas White committed
181
	int config_zoneaxis = 0;
182
183
184
185
186

	/* Long options */
	const struct option longopts[] = {
		{"help",               0, NULL,               'h'},
		{"input",              1, NULL,               'i'},
Thomas White's avatar
Thomas White committed
187
		{"max-only",           0, &config_maxonly,     1},
Thomas White's avatar
Thomas White committed
188
		{"output-every",       1, NULL,               'e'},
Thomas White's avatar
Thomas White committed
189
		{"rvsq",               0, NULL,               'r'},
Thomas White's avatar
Thomas White committed
190
		{"stop-after",         1, NULL,               's'},
Thomas White's avatar
Thomas White committed
191
		{"zone-axis",          0, &config_zoneaxis,    1},
192
193
194
195
		{0, 0, NULL, 0}
	};

	/* Short options */
Thomas White's avatar
Thomas White committed
196
	while ((c = getopt_long(argc, argv, "hi:e:r", longopts, NULL)) != -1) {
197
198
199
200
201
202
203
204
205
206
207
208

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

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

Thomas White's avatar
Thomas White committed
209
210
211
212
213
		case 'r' : {
			config_rvsq = 1;
			break;
		}

Thomas White's avatar
Thomas White committed
214
215
216
217
218
		case 'e' : {
			config_every = atoi(optarg);
			break;
		}

Thomas White's avatar
Thomas White committed
219
220
221
222
223
		case 's' : {
			config_stopafter = atoi(optarg);
			break;
		}

224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
		case 0 : {
			break;
		}

		default : {
			return 1;
		}
		}

	}

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

	mol = load_molecule();
	get_reflections_cached(mol, eV_to_J(2.0e3));

	ref = new_list_intensity();
	counts = new_list_count();
	trueref = ideal_intensities(mol->reflections);

	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
262
		signed int h, k, l, intensity;
263
264
265
266
		int r;

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

Thomas White's avatar
Thomas White committed
268
269
270
271
			if ( n_patterns == 0 ) {
				n_patterns++;
				continue;
			}
Thomas White's avatar
Tidy-up    
Thomas White committed
272

273
			if (config_every && (n_patterns % config_every == 0)) {
274
				process_reflections(ref, trueref, counts,
Thomas White's avatar
Thomas White committed
275
				                    n_patterns, mol->cell,
Thomas White's avatar
Thomas White committed
276
277
				                    config_rvsq,
				                    config_zoneaxis);
278
			}
Thomas White's avatar
Thomas White committed
279
280
281

			if ( n_patterns == config_stopafter ) break;

Thomas White's avatar
Thomas White committed
282
			n_patterns++;
283
284
285
286
287
		}

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

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

Thomas White's avatar
Thomas White committed
290
291
292
293
294
295
296
297
298
		if ( !config_maxonly ) {
			integrate_intensity(ref, h, k, l, intensity);
			integrate_count(counts, h, k, l, 1);
		} else {
			if ( intensity > lookup_intensity(ref, h, k, l) ) {
				set_intensity(ref, h, k, l, intensity);
			}
			set_count(counts, h, k, l, 1);
		}
299
300
301
302
303

	} while ( rval != NULL );

	fclose(fh);

304
305
306
	process_reflections(ref, trueref, counts, n_patterns, mol->cell,
	                    config_rvsq, config_zoneaxis);

Thomas White's avatar
Thomas White committed
307
308
309
	write_reflections("results/reflections.hkl", counts, ref, 0, NULL);

	STATUS("There were %u patterns.\n", n_patterns);
310
311
312

	return 0;
}