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