diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/compare_hkl.c | 8 | ||||
| -rw-r--r-- | src/statistics.c | 63 | ||||
| -rw-r--r-- | src/statistics.h | 4 |
3 files changed, 72 insertions, 3 deletions
diff --git a/src/compare_hkl.c b/src/compare_hkl.c index b14a8dab..b2b40da1 100644 --- a/src/compare_hkl.c +++ b/src/compare_hkl.c @@ -50,7 +50,7 @@ int main(int argc, char *argv[]) char *afile = NULL; char *bfile = NULL; signed int h, k, l; - double scale, R; + double scale, R2, Rmerge; unsigned int *c1; unsigned int *c2; int i; @@ -141,8 +141,10 @@ int main(int argc, char *argv[]) ncom += c1[i] && c2[i]; } STATUS("%i,%i reflections: %i in common\n", nc1, nc2, ncom); - R = stat_r2(ref1, c1, ref2, c2, &scale); - STATUS("R2 = %5.4f %% (scale=%5.2f)\n", R*100.0, scale); + R2 = stat_r2(ref1, c1, ref2, c2, &scale); + STATUS("R2 = %5.4f %% (scale=%5.2f)\n", R2*100.0, scale); + Rmerge = stat_rmerge(ref1, c1, ref2, c2, &scale); + STATUS("Rmerge = %5.4f %% (scale=%5.2f)\n", Rmerge*100.0, scale); if ( outfile != NULL ) { write_reflections(outfile, NULL, out, NULL, 1, cell, 1); diff --git a/src/statistics.c b/src/statistics.c index b1589ebe..d78dde9e 100644 --- a/src/statistics.c +++ b/src/statistics.c @@ -73,3 +73,66 @@ double stat_r2(const double *ref1, const unsigned int *c1, return sqrt(top/bot); } + + +double stat_scale_sqrti(const double *ref1, const unsigned int *c1, + const double *ref2, const unsigned int *c2) +{ + double top = 0.0; + double bot = 0.0; + int i; + + /* Start from i=1 -> skip central beam */ + for ( i=1; i<LIST_SIZE; i++ ) { + + if ( c1[i] && c2[i] ) { + + double f1, f2; + + if ( (ref1[i]<0.0) || (ref2[i]<0.0) ) continue; + + f1 = sqrt(ref1[i]) / (double)c1[i]; + f2 = sqrt(ref2[i]) / (double)c2[i]; + + top += f1 * f2; + bot += f2 * f2; + + } /* else reflection not common so don't worry about it */ + + } + + return top/bot; +} + + +double stat_rmerge(const double *ref1, const unsigned int *c1, + const double *ref2, const unsigned int *c2, double *scalep) +{ + double top = 0.0; + double bot = 0.0; + double scale; + int i; + scale = stat_scale_sqrti(ref1, c1, ref2, c2); + *scalep = scale; + + /* Start from i=1 -> skip central beam */ + for ( i=1; i<LIST_SIZE; i++ ) { + + if ( c1[i] && c2[i] ) { + + double f1, f2; + + if ( (ref1[i]<0.0) || (ref2[i]<0.0) ) continue; + + f1 = sqrt(ref1[i]) / (scale*(double)c1[i]); + f2 = sqrt(ref2[i]) / (double)c2[i]; + + top += fabs(f1 - f2); + bot += f1 + f2; + + } /* else reflection not measured so don't worry about it */ + + } + + return 2.0*top/bot; +} diff --git a/src/statistics.h b/src/statistics.h index fe7f57e1..e19d4b29 100644 --- a/src/statistics.h +++ b/src/statistics.h @@ -23,4 +23,8 @@ extern double stat_r2(const double *ref1, const unsigned int *c1, const double *ref2, const unsigned int *c2, double *scalep); +extern double stat_rmerge(const double *ref1, const unsigned int *c1, + const double *ref2, const unsigned int *c2, + double *scalep); + #endif /* STATISTICS_H */ |