diff options
| -rw-r--r-- | src/process_hkl.c | 70 | ||||
| -rw-r--r-- | src/statistics.c | 12 |
2 files changed, 59 insertions, 23 deletions
diff --git a/src/process_hkl.c b/src/process_hkl.c index 097a3eb3..f1325e06 100644 --- a/src/process_hkl.c +++ b/src/process_hkl.c @@ -26,12 +26,23 @@ #include "sfac.h" +/* Number of divisions for R vs |q| graphs */ +#define RVQDV (20) + + static void show_help(const char *s) { printf("Syntax: %s [options]\n\n", s); - printf("Assemble and process FEL Bragg intensities.\n\n"); - printf(" -h Display this help message.\n"); - printf(" -i <filename> Specify input filename (\"-\" for stdin).\n"); + printf( +"Assemble and process FEL Bragg intensities.\n" +"\n" +" -h, --help Display this help message.\n" +" -i, --input=<filename> Specify input filename (\"-\" for stdin).\n" +"\n" +" --max-only Take the integrated intensity to be equal to the\n" +" maximum intensity measured for that reflection.\n" +" The default is to use the mean value from all\n" +" measurements.\n"); } @@ -56,11 +67,14 @@ static void write_RvsQ(const char *name, double *ref, double *trueref, } } - for ( sbracket=0.0; sbracket<smax; sbracket+=smax/10.0 ) { + for ( sbracket=0.0; sbracket<smax; sbracket+=smax/RVQDV ) { double top = 0.0; double bot = 0.0; - int n = 0; + int nhits = 0; + int nrefl = 0; + double R; + double hits_per_refl; for ( h=-INDMAX; h<INDMAX; h++ ) { for ( k=-INDMAX; k<INDMAX; k++ ) { @@ -68,9 +82,12 @@ static void write_RvsQ(const char *name, double *ref, double *trueref, double s; int c; + + if ( (h==0) && (k==0) && (l==0) ) continue; + c = lookup_count(counts, h, k, l); s = resolution(cell, h, k, l); - if ((s>=sbracket) && (s<sbracket+smax/10.0) && (c>0)) { + if ((s>=sbracket) && (s<sbracket+smax/RVQDV) && (c>0)) { double obs, calc, obsi; @@ -78,16 +95,22 @@ static void write_RvsQ(const char *name, double *ref, double *trueref, calc = lookup_intensity(trueref, h, k, l); obsi = obs / (double)c; - top += fabs(obsi/scale - calc); - bot += obsi/scale; - n++; + top += pow(obsi - scale*calc, 2.0); + bot += pow(obsi, 2.0); + nhits += c; + nrefl++; + } } } } - fprintf(fh, "%8.5f %8.5f %i\n", sbracket+smax/20.0, top/bot, n); + R = sqrt(top/bot); + hits_per_refl = nrefl ? (double)nhits/nrefl : 0; + + fprintf(fh, "%8.5f %8.5f %5.2f\n", sbracket+smax/(2.0*RVQDV), + R, hits_per_refl); } fclose(fh); @@ -146,9 +169,9 @@ static void process_reflections(double *ref, double *trueref, double mean_counts; int ctot = 0; int nmeas = 0; - double ff; char name[64]; double R, scale; + double calc_222, obs_222; for ( j=0; j<LIST_SIZE; j++ ) { ctot += counts[j]; @@ -156,13 +179,13 @@ static void process_reflections(double *ref, double *trueref, } mean_counts = (double)ctot/nmeas; - ff = lookup_intensity(ref, 2, 2, 2) / lookup_count(counts, 2, 2, 2); + calc_222 = lookup_intensity(ref, 2, 2, 2) / lookup_count(counts, 2, 2, 2); + obs_222 = lookup_intensity(trueref, 2, 2, 2); R = stat_r2(ref, trueref, counts, LIST_SIZE, &scale); - STATUS("%8u: R=%5.2f%% (sf=%7.4e)" - " mean meas/refl=%5.2f," - " %i reflections measured. %f\n", - n_patterns, R*100.0, scale, mean_counts, nmeas, ff); + STATUS("%8u: R=%5.2f%% (sf=%7.4e) mean meas/refl=%5.2f," + " %i reflections measured, %f\n", + n_patterns, R*100.0, scale, mean_counts, nmeas, calc_222/obs_222); /* Record graph of R against q for this N */ snprintf(name, 63, "results/R_vs_q-%u.dat", n_patterns); @@ -180,11 +203,13 @@ int main(int argc, char *argv[]) unsigned int *counts; char *rval; struct molecule *mol; + int config_maxonly = 0; /* Long options */ const struct option longopts[] = { {"help", 0, NULL, 'h'}, {"input", 1, NULL, 'i'}, + {"max-only", 0, &config_maxonly, 1}, {0, 0, NULL, 0} }; @@ -255,8 +280,17 @@ int main(int argc, char *argv[]) r = sscanf(line, "%i %i %i %i", &h, &k, &l, &intensity); if ( r != 4 ) continue; - integrate_intensity(ref, h, k, l, intensity); - integrate_count(counts, h, k, l, 1); + //if ( (abs(h)>3) || (abs(k)>3) || (abs(l)>3) ) continue; + + if ( !config_maxonly ) { + integrate_intensity(ref, h, k, l, intensity); + integrate_count(counts, h, k, l, 1); + } else { + if ( intensity > lookup_intensity(ref, h, k, l) ) { + set_intensity(ref, h, k, l, intensity); + } + set_count(counts, h, k, l, 1); + } } while ( rval != NULL ); diff --git a/src/statistics.c b/src/statistics.c index 9219b069..50bfd513 100644 --- a/src/statistics.c +++ b/src/statistics.c @@ -29,7 +29,8 @@ static double stat_scale_intensity(double *obs, double *calc, unsigned int *c, double bot = 0.0; int i; - for ( i=0; i<size; i++ ) { + /* Start from i=1 -> skip central beam */ + for ( i=1; i<size; i++ ) { if ( c[i] > 0 ) { double obsi; @@ -54,16 +55,17 @@ double stat_r2(double *obs, double *calc, unsigned int *c, int size, scale = stat_scale_intensity(obs, calc, c, size); *scalep = scale; - for ( i=0; i<size; 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 += fabs(obsi/scale - calc[i]); - bot += obsi/scale; + top += pow(fabs(obsi - scale*calc[i]), 2.0); + bot += pow(obsi, 2.0); } } /* else reflection not measured so don't worry about it */ - return top/bot; + return sqrt(top/bot); } |