aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2014-09-18 18:10:11 +0200
committerThomas White <taw@physics.org>2014-09-18 18:10:11 +0200
commit0a533227b1006b4019ea9fccd0957e8be7138f8d (patch)
treed2ff9c81c5ececd501f404388a620eaee9a57f19 /src
parent0b1fa4693e4c4df54623cfeb91b3c887191236fa (diff)
Tidy up scaling/CCs a bit
Diffstat (limited to 'src')
-rw-r--r--src/process_hkl.c94
1 files changed, 53 insertions, 41 deletions
diff --git a/src/process_hkl.c b/src/process_hkl.c
index 5bf44124..9abfc501 100644
--- a/src/process_hkl.c
+++ b/src/process_hkl.c
@@ -85,7 +85,6 @@ static void show_help(const char *s)
" --min-snr=<n> Require individual intensity measurements to\n"
" have I > n * sigma(I). Default: -infinity.\n"
" --min-cc=<n> Reject frames with CC less than n. Default: infinity.\n"
-" have I > n * sigma(I). Default: -infinity.\n"
" --max-adu=<n> Maximum peak value. Default: infinity.\n"
" --min-res=<n> Merge only crystals which diffract above <n> A.\n"
" --push-res=<n> Integrate higher than apparent resolution cutoff.\n"
@@ -182,10 +181,11 @@ static double scale_intensities(RefList *reference, RefList *new,
return s;
}
+
static double cc_intensities(RefList *reference, RefList *new,
- const SymOpList *sym)
+ const SymOpList *sym)
{
- // x is reference
+ /* "x" is "reference" */
float s_xy = 0.0;
float s_x = 0.0;
float s_y = 0.0;
@@ -234,6 +234,7 @@ static double cc_intensities(RefList *reference, RefList *new,
return (s_xy - s_x*s_y/n) / sqrt(t1*t2);
}
+
static double *check_hist_size(int n, double *hist_vals)
{
int ns;
@@ -257,32 +258,38 @@ static int merge_crystal(RefList *model, struct image *image, Crystal *cr,
double **hist_vals, signed int hist_h,
signed int hist_k, signed int hist_l, int *hist_n,
int config_nopolar, double min_snr, double max_adu,
- double push_res, double min_cc, int do_scale, FILE *stat)
+ double push_res, double min_cc, int do_scale,
+ FILE *stat)
{
Reflection *refl;
RefListIterator *iter;
+ RefList *new_refl;
double scale, cc;
+ new_refl = crystal_get_reflections(cr);
+
/* First, correct for polarisation */
if ( !config_nopolar ) {
- polarisation_correction(crystal_get_reflections(cr),
- crystal_get_cell(cr), image);
+ polarisation_correction(new_refl, crystal_get_cell(cr), image);
}
if ( reference != NULL ) {
- scale = scale_intensities(reference,
- crystal_get_reflections(cr), sym);
- cc = cc_intensities(reference, crystal_get_reflections(cr), sym);
- if (stat != NULL) {
- fprintf(stat, "%s,%f,%f\n",image->filename, scale, cc);
+ if ( do_scale ) {
+ scale = scale_intensities(reference, new_refl, sym);
+ } else {
+ scale = 1.0;
}
- if (cc < min_cc || scale <= 0) return 1;
+ cc = cc_intensities(reference, new_refl, sym);
+ if ( cc < min_cc ) return 1;
+ if ( isnan(scale) ) return 1;
+ if ( scale <= 0.0 ) return 1;
+ if ( stat != NULL ) {
+ fprintf(stat, "%s %f %f\n", image->filename, scale, cc);
+ }
+
} else {
scale = 1.0;
}
- if (!do_scale) scale = 1.0;
-
- if ( isnan(scale) ) return 1;
for ( refl = first_refl(crystal_get_reflections(cr), &iter);
refl != NULL;
@@ -391,11 +398,12 @@ static int merge_all(Stream *st, RefList *model, RefList *reference,
int n_crystals_seen = 0;
FILE *stat = NULL;
- if (stat_output != NULL) {
- stat = fopen(stat_output, "w");
- if (stat == NULL) {
- ERROR("Failed to open statistics output file %s\n", stat_output);
- }
+ if ( stat_output != NULL ) {
+ stat = fopen(stat_output, "w");
+ if ( stat == NULL ) {
+ ERROR("Failed to open statistics output file %s\n",
+ stat_output);
+ }
}
do {
@@ -467,8 +475,9 @@ static int merge_all(Stream *st, RefList *model, RefList *reference,
set_esd_intensity(refl, sqrt(var)/sqrt(red));
}
- if (stat != NULL)
- fclose(stat);
+ if ( stat != NULL ) {
+ fclose(stat);
+ }
return 0;
}
@@ -504,7 +513,7 @@ int main(int argc, char *argv[])
double min_res = 0.0;
double push_res = +INFINITY;
double min_cc = -INFINITY;
- int do_scale = 0;
+ int twopass = 0;
/* Long options */
const struct option longopts[] = {
@@ -513,7 +522,7 @@ int main(int argc, char *argv[])
{"output", 1, NULL, 'o'},
{"start-after", 1, NULL, 's'},
{"stop-after", 1, NULL, 'f'},
- {"scale", 0, &do_scale, 1},
+ {"scale", 0, &config_scale, 1},
{"no-polarisation", 0, &config_nopolar, 1},
{"no-polarization", 0, &config_nopolar, 1},
{"symmetry", 1, NULL, 'y'},
@@ -549,10 +558,6 @@ int main(int argc, char *argv[])
output = strdup(optarg);
break;
- case 9 :
- stat_output = strdup(optarg);
- break;
-
case 's' :
errno = 0;
start_after = strtod(optarg, &rval);
@@ -649,7 +654,12 @@ int main(int argc, char *argv[])
ERROR("Invalid value for --min-cc.\n");
return 1;
}
- config_scale = 1;
+ twopass = 1;
+ break;
+
+ case 9 :
+ stat_output = strdup(optarg);
+ twopass = 1;
break;
case 0 :
@@ -663,8 +673,6 @@ int main(int argc, char *argv[])
}
- if (do_scale) config_scale = 1;
-
if ( filename == NULL ) {
ERROR("Please specify filename using the -i option\n");
return 1;
@@ -729,17 +737,21 @@ int main(int argc, char *argv[])
}
+ /* Need to do a second pass if we are scaling */
+ if ( config_scale ) twopass = 1;
+
hist_i = 0;
r = merge_all(st, model, NULL, sym, &hist_vals, hist_h, hist_k, hist_l,
&hist_i, config_nopolar, min_measurements, min_snr,
- max_adu, start_after, stop_after, min_res, push_res, min_cc, do_scale, stat_output);
+ max_adu, start_after, stop_after, min_res, push_res,
+ min_cc, config_scale, stat_output);
fprintf(stderr, "\n");
if ( r ) {
ERROR("Error while reading stream.\n");
return 1;
}
- if ( config_scale ) {
+ if ( twopass ) {
RefList *reference;
@@ -752,7 +764,7 @@ int main(int argc, char *argv[])
int r;
- STATUS("Extra pass for scaling...\n");
+ STATUS("Second pass for scaling and/or CCs...\n");
reference = model;
model = reflist_new();
@@ -763,11 +775,12 @@ int main(int argc, char *argv[])
hist_i = 0;
}
- r = merge_all(st, model, reference, sym,
- &hist_vals, hist_h, hist_k, hist_l, &hist_i,
- config_nopolar, min_measurements, min_snr,
- max_adu, start_after, stop_after, min_res,
- push_res, min_cc, do_scale, stat_output);
+ r = merge_all(st, model, reference, sym, &hist_vals,
+ hist_h, hist_k, hist_l, &hist_i,
+ config_nopolar, min_measurements, min_snr,
+ max_adu, start_after, stop_after, min_res,
+ push_res, min_cc, config_scale,
+ stat_output);
fprintf(stderr, "\n");
if ( r ) {
ERROR("Error while reading stream.\n");
@@ -799,8 +812,7 @@ int main(int argc, char *argv[])
reflist_free(model);
free(output);
free(filename);
- if (stat_output != NULL)
- free(stat_output);
+ free(stat_output);
return 0;
}