aboutsummaryrefslogtreecommitdiff
path: root/src/statistics.c
blob: e9aaff98a349eb459404987a4846f4caab1da2df (plain)
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
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
/*
 * statistics.c
 *
 * Structure-factor statistics
 *
 * (c) 2007-2009 Thomas White <taw@physics.org>
 *
 * Part of CrystFEL - crystallography with a FEL
 *
 */


#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <math.h>
#include <stdlib.h>

#include "statistics.h"


/* By what (best-fitted) factor must the list "list" be multiplied by,
 * if it were to be merged with "target"? */
static double stat_scale_intensity(double *obs, double *calc, unsigned int *c,
                                   int size)
{
	double top = 0.0;
	double bot = 0.0;
	int i;

	/* Start from i=1 -> skip central beam */
	for ( i=1; i<size; i++ ) {

		if ( c[i] > 0 ) {
			double obsi;
			obsi = obs[i] / (double)c[i];
			top += obsi * calc[i];
			bot += calc[i] * calc[i];
		} /* else reflection not measured so don't worry about it */

	}

	return top/bot;
}


double stat_r2(double *obs, double *calc, unsigned int *c, int size,
               double *scalep)
{
	double top = 0.0;
	double bot = 0.0;
	double scale;
	int i;
	scale = stat_scale_intensity(obs, calc, c, size);
	*scalep = scale;

	/* Start from i=1 -> skip central beam */
	for ( i=1; i<size; i++ ) {

		if ( c[i] > 0 ) {

			double obsi;

			obsi = obs[i] / (double)c[i];
			obsi = obsi / scale;

			top += pow(fabs(obsi - calc[i]), 2.0);
			bot += pow(obsi, 2.0);

		} /* else reflection not measured so don't worry about it */

	}

	return sqrt(top/bot);
}