aboutsummaryrefslogtreecommitdiff
path: root/src/check_hkl.c
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2010-11-16 16:03:49 +0100
committerThomas White <taw@physics.org>2012-02-22 15:27:05 +0100
commit757b17f432be24c47fed65a2dd569fc153535ff7 (patch)
tree110188ae64e76666e9c643c8ff09ba1252db0188 /src/check_hkl.c
parentcca20bccf70791f7ce9bd27c59a6c048716012fb (diff)
check_hkl: Show mean and variance
Diffstat (limited to 'src/check_hkl.c')
-rw-r--r--src/check_hkl.c67
1 files changed, 56 insertions, 11 deletions
diff --git a/src/check_hkl.c b/src/check_hkl.c
index 78e97a3d..650895ec 100644
--- a/src/check_hkl.c
+++ b/src/check_hkl.c
@@ -179,7 +179,7 @@ static void plot_shells(const double *ref, ReflItemList *items, UnitCell *cell,
}
free(counted);
- /* Characterise the data set */
+ /* Calculate means */
for ( i=0; i<num_items(items); i++ ) {
struct refl_item *it;
@@ -206,14 +206,52 @@ static void plot_shells(const double *ref, ReflItemList *items, UnitCell *cell,
}
measured[bin]++;
+ mean[bin] += lookup_intensity(ref, h, k, l);
+
+ }
+
+ for ( i=0; i<NBINS; i++ ) {
+ mean[i] /= (double)measured[i];
+ }
+
+ /* Characterise the data set */
+ for ( i=0; i<num_items(items); i++ ) {
+
+ struct refl_item *it;
+ signed int h, k, l;
+ double d;
+ int bin;
+ int j;
+ double val, esd;
+
+ it = get_item(items, i);
+ h = it->h; k = it->k; l = it->l;
+
+ d = resolution(cell, h, k, l) * 2.0;
+ val = lookup_intensity(ref, h, k, l);
+ esd = lookup_intensity(sigma, h, k, l);
+
+ bin = -1;
+ for ( j=0; j<NBINS; j++ ) {
+ if ( (d>rmins[j]) && (d<=rmaxs[j]) ) {
+ bin = j;
+ break;
+ }
+ }
+ if ( bin == -1 ) {
+ nout++;
+ continue;
+ }
+
+ /* measured[bin] was done earlier */
measurements[bin] += lookup_count(counts, h, k, l);
- snr[bin] += (lookup_intensity(ref1, h, k, l) /
- lookup_intensity(sigma, h, k, l));
- snr_total += (lookup_intensity(ref1, h, k, l) /
- lookup_intensity(sigma, h, k, l));
+ snr[bin] += val / esd;
+ snr_total += val / esd;
nmeas++;
nmeastot += lookup_count(counts, h, k, l);
+ var[bin] += pow(val-mean[bin], 2.0);
+
}
STATUS("overall <snr> = %f\n", snr_total/(double)nmeas);
STATUS("%i measurements in total.\n", nmeastot);
@@ -224,15 +262,22 @@ static void plot_shells(const double *ref, ReflItemList *items, UnitCell *cell,
nout);
}
+ fprintf(fh, "1/d centre # refs Possible Compl "
+ "Meas Red SNR Std dev Mean\n");
for ( i=0; i<NBINS; i++ ) {
- double r, cen;
+ double cen;
cen = rmins[i] + (rmaxs[i] - rmins[i])/2.0;
- r = (num[i]/den)*((double)ctot/cts[i]);
- fprintf(fh, "%f %i %i %5.2f %i %f %f\n", cen*1.0e-9, measured[i],
- possible[i], 100.0*(float)measured[i]/possible[i],
- measurements[i], (float)measurements[i]/measured[i],
- (snr[i]/(double)measured[i]));
+ fprintf(fh, "%10.3f %8i %8i %6.2f %10i %5.1f %5.2f %10.2f %10.2f\n",
+ cen*1.0e-9,
+ measured[i],
+ possible[i],
+ 100.0*(float)measured[i]/possible[i],
+ measurements[i],
+ (float)measurements[i]/measured[i],
+ snr[i]/(double)measured[i],
+ sqrt(var[i]/measured[i]),
+ mean[i]);
}