diff options
-rw-r--r-- | src/check_hkl.c | 28 |
1 files changed, 21 insertions, 7 deletions
diff --git a/src/check_hkl.c b/src/check_hkl.c index fa5fb0b7..649b86af 100644 --- a/src/check_hkl.c +++ b/src/check_hkl.c @@ -84,9 +84,10 @@ static void plot_shells(RefList *list, UnitCell *cell, const SymOpList *sym, int i; FILE *fh; double snr_total = 0; - int nmeas = 0; + int nrefl = 0; int nmeastot = 0; int nout = 0; + int nsilly = 0; Reflection *refl; RefListIterator *iter; RefList *counted; @@ -193,13 +194,15 @@ static void plot_shells(RefList *list, UnitCell *cell, const SymOpList *sym, refl = next_refl(refl, iter) ) { signed int h, k, l; - double d; + double d, val, esd; int bin; int j; get_indices(refl, &h, &k, &l); d = resolution(cell, h, k, l) * 2.0; + val = get_intensity(refl); + esd = get_esd_intensity(refl); bin = -1; for ( j=0; j<NBINS; j++ ) { @@ -210,6 +213,11 @@ static void plot_shells(RefList *list, UnitCell *cell, const SymOpList *sym, } if ( bin == -1 ) continue; + if ( !isfinite(val/esd) ) { + nsilly++; + continue; + } + measured[bin]++; mean[bin] += get_intensity(refl); } @@ -253,20 +261,26 @@ static void plot_shells(RefList *list, UnitCell *cell, const SymOpList *sym, measurements[bin] += get_redundancy(refl); snr[bin] += val / esd; snr_total += val / esd; - nmeas++; + nrefl++; nmeastot += get_redundancy(refl); var[bin] += pow(val-mean[bin], 2.0); } - STATUS("overall <snr> = %f\n", snr_total/(double)nmeas); + STATUS("overall <snr> = %f\n", snr_total/(double)nrefl); STATUS("%i measurements in total.\n", nmeastot); - STATUS("%i reflections in total.\n", nmeas); + STATUS("%i reflections in total.\n", nrefl); if ( nout ) { STATUS("Warning; %i reflections outside resolution range.\n", nout); } + if ( nsilly ) { + STATUS("Warning; %i reflections had infinite or invalid values" + " of I/sigma(I).\n", nsilly); + } + + fprintf(fh, "1/d centre # refs Possible Compl " "Meas Red SNR Std dev Mean\n"); for ( i=0; i<NBINS; i++ ) { @@ -278,9 +292,9 @@ static void plot_shells(RefList *list, UnitCell *cell, const SymOpList *sym, cen*1.0e-9, measured[i], possible[i], - 100.0*(float)measured[i]/possible[i], + 100.0*(double)measured[i]/possible[i], measurements[i], - (float)measurements[i]/measured[i], + (double)measurements[i]/measured[i], snr[i]/(double)measured[i], sqrt(var[i]/measured[i]), mean[i]); |