diff options
| -rw-r--r-- | src/check_hkl.c | 67 |
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]); } |