aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--src/process_hkl.c84
1 files changed, 59 insertions, 25 deletions
diff --git a/src/process_hkl.c b/src/process_hkl.c
index f94ad4b8..ac2e567e 100644
--- a/src/process_hkl.c
+++ b/src/process_hkl.c
@@ -32,8 +32,7 @@
#include "image.h"
-/* Number of divisions for intensity histograms */
-#define NBINS (50)
+
static void show_help(const char *s)
@@ -63,6 +62,8 @@ static void show_help(const char *s)
" the default.\n"
" -g, --histogram=<h,k,l> Calculate the histogram of measurements for this\n"
" reflection.\n"
+" -z, --hist-parameters Set the range for the histogram and the number of\n"
+" =<min,max,nbins> bins. \n"
"\n"
" --scale Scale each pattern for best fit with the current\n"
" model.\n"
@@ -75,13 +76,13 @@ static void show_help(const char *s)
}
-static void plot_histogram(double *vals, int n)
+static void plot_histogram(double *vals, int n, float hist_min, float hist_max, int nbins)
{
int i;
double max = -INFINITY;
double min = +INFINITY;
double step;
- int histo[NBINS];
+ int histo[nbins];
FILE *fh;
fh = fopen("histogram.dat", "w");
@@ -90,26 +91,33 @@ static void plot_histogram(double *vals, int n)
return;
}
- for ( i=0; i<n; i++ ) {
- if ( vals[i] > max ) max = vals[i];
- if ( vals[i] < min ) min = vals[i];
+ if ( hist_min == hist_max ) {
+ for ( i=0; i<n; i++ ) {
+ if ( vals[i] > max ) max = vals[i];
+ if ( vals[i] < min ) min = vals[i];
+ }
+ } else {
+ min = hist_min;
+ max = hist_max;
}
- STATUS("%f %f\n", min, max);
+ STATUS("min max nbins: \n %f %f %i\n", min, max, nbins);
min--; max++;
- for ( i=0; i<NBINS; i++ ) {
+ for ( i=0; i<nbins; i++ ) {
histo[i] = 0;
}
- step = (max-min)/NBINS;
+ step = (max-min)/nbins;
for ( i=0; i<n; i++ ) {
int bin;
- bin = (vals[i]-min)/step;
- histo[bin]++;
+ if ( (vals[i] > min) && (vals[i] < max) ) {
+ bin = (vals[i]-min)/step;
+ histo[bin]++;
+ }
}
- for ( i=0; i<NBINS; i++ ) {
+ for ( i=0; i<nbins; i++ ) {
fprintf(fh, "%f %i\n", min+step*i, histo[i]);
}
@@ -167,6 +175,7 @@ static void merge_pattern(RefList *model, RefList *new, int max_only,
int cur_redundancy = get_redundancy(model_version);
set_redundancy(model_version, cur_redundancy+1);
+
} else if ( pass == 2 ) {
double dev = get_temp1(model_version);
@@ -420,10 +429,13 @@ int main(int argc, char *argv[])
char *pdb = NULL;
char *histo = NULL;
signed int hist_h, hist_k, hist_l;
+ signed int hist_nbins=50;
+ float hist_min=0.0, hist_max=0.0;
double *hist_vals = NULL;
int hist_i;
struct beam_params *beam = NULL;
int space_for_hist = 0;
+ char *histo_params = NULL;
/* Long options */
const struct option longopts[] = {
@@ -439,12 +451,13 @@ int main(int argc, char *argv[])
{"symmetry", 1, NULL, 'y'},
{"pdb", 1, NULL, 'p'},
{"histogram", 1, NULL, 'g'},
+ {"hist-parameters", 1, NULL, 'z'},
{"beam", 1, NULL, 'b'},
{0, 0, NULL, 0}
};
/* Short options */
- while ((c = getopt_long(argc, argv, "hi:e:o:p:y:g:f:b:",
+ while ((c = getopt_long(argc, argv, "hi:e:o:p:y:g:f:b:z:",
longopts, NULL)) != -1) {
switch (c) {
@@ -480,6 +493,10 @@ int main(int argc, char *argv[])
histo = strdup(optarg);
break;
+ case 'z' :
+ histo_params = strdup(optarg);
+ break;
+
case 'b' :
beam = get_beam_parameters(optarg);
if ( beam == NULL ) {
@@ -562,39 +579,56 @@ int main(int argc, char *argv[])
STATUS("%i %i %i\n", hist_h, hist_k, hist_l);
}
+ if ( histo_params != NULL ) {
+ int rr;
+ rr = sscanf(histo_params, "%f,%f,%i", &hist_min, &hist_max, &hist_nbins);
+ if ( rr != 3 ) {
+ ERROR("Invalid parameters for '--hist-parameters'\n");
+ return 1;
+ }
+ free(histo_params);
+ if ( (hist_max - hist_min) <=0 ) {
+ ERROR("Invalid range for '--hist-parameters' : check if min<max \n");
+ return 1;
+ }
+ }
+
hist_i = 0;
merge_all(fh, model, config_maxonly, config_scale, config_sum,
config_startafter, config_stopafter,
sym, n_total_patterns,
- hist_vals, hist_h, hist_k, hist_l, &hist_i, 1);
+ NULL, 0, 0, 0, NULL, 1);
if ( ferror(fh) ) {
ERROR("Stream read error.\n");
return 1;
}
rewind(fh);
- if ( space_for_hist && (hist_i >= space_for_hist) ) {
- ERROR("Histogram array was too small!\n");
- }
-
- if ( hist_vals != NULL ) {
- STATUS("%i %i %i was seen %i times.\n", hist_h, hist_k, hist_l,
- hist_i);
- plot_histogram(hist_vals, hist_i);
- }
STATUS("Extra pass to calculate ESDs...\n");
rewind(fh);
merge_all(fh, model, config_maxonly, config_scale, 0,
config_startafter, config_stopafter, sym, n_total_patterns,
- NULL, 0, 0, 0, NULL, 2);
+ hist_vals, hist_h, hist_k, hist_l, &hist_i, 2);
if ( ferror(fh) ) {
ERROR("Stream read error.\n");
return 1;
}
+ if ( space_for_hist && (hist_i >= space_for_hist) ) {
+ ERROR("Histogram array was too small!\n");
+ }
+
+ if ( hist_vals != NULL ) {
+ STATUS("%i %i %i was seen %i times.\n", hist_h, hist_k, hist_l,
+ hist_i);
+ plot_histogram(hist_vals, hist_i, hist_min, hist_max,
+ hist_nbins);
+ }
+
write_reflist(output, model, cell);
fclose(fh);
+
free(sym);
reflist_free(model);
free(output);