aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2009-12-01 15:48:00 +0100
committerThomas White <taw@physics.org>2009-12-01 15:48:00 +0100
commit1709a699e5d31da433798b8a312e2a7098c5fa37 (patch)
treec2ac81469d418d55f2166ea4780ebf48bd57cf06
parent7025bd8293fcf72f7287ee9eee661708aad2895c (diff)
Various improvements
-rw-r--r--src/process_hkl.c70
-rw-r--r--src/statistics.c12
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);
}