diff options
| -rwxr-xr-x | scripts/frequency | 165 | ||||
| -rw-r--r-- | src/process_hkl.c | 104 |
2 files changed, 114 insertions, 155 deletions
diff --git a/scripts/frequency b/scripts/frequency index 26e29bb4..d6caf72c 100755 --- a/scripts/frequency +++ b/scripts/frequency @@ -2,168 +2,37 @@ use strict; -my $bins = 200; my $nplots = 32; +my $sym = "6/mmm"; -open(FH, $ARGV[0]); +my $stream = $ARGV[0]; -my $line; -my %counts; - -# --- Step 1: Count reflections and find the most popular ones - -while ( $line = <FH> ) { - - my $h; - my $k; - my $l; - - chomp($line); - - if ( $line =~ /^\s+(\d+)\s+(\d+)\s+(\d+)/ ) { - - my $key; - - $h = $1; - $k = $2; - $l = $3; - - $key = $h.",".$k.",".$l; - $counts{$key}++; - - } - -} - -my $ref; -my $i; -$i = 0; -foreach $ref ( sort {-($counts{$a} <=> $counts{$b})} keys %counts ) -{ - $i++; - if ( $i > $nplots ) { - delete($counts{$ref}); - } else { - printf("%s %i\n", $ref, $counts{$ref}); - } -} - - -# --- Step 2: Work out histogram bins for each reflection - -seek(FH, 0, 0); - -my %mins; -my %maxs; +#system("~/crystfel/src/process_hkl -i ".$stream." -o test.hkl -y ".$sym); +open(FH, "test.hkl"); +my $max = 0; +my $h; +my $k; +my $l; +my $line; while ( $line = <FH> ) { - my $key; - my $intensity; - - chomp($line); - - if ( $line =~ /^\s+(\d+)\s+(\d+)\s+(\d+)\s+([\d\.]+)/ ) { - - my $h; - my $k; - my $l; + chomp $line; - $h = $1; - $k = $2; - $l = $3; - $intensity = $4; - - $key = $h.",".$k.",".$l; - - next unless exists($counts{$key}); - - unless ( exists($mins{$key}) ) { - $mins{$key} = $intensity; - } - unless ( exists($maxs{$key}) ) { - $maxs{$key} = $intensity; - } - if ( $intensity < $mins{$key} ) { - $mins{$key} = $intensity; - } - if ( $intensity > $maxs{$key} ) { - $maxs{$key} = $intensity; - } - - } - -} - -# --- Step 3: Bin the counts and plot - -open(GP, "| gnuplot"); -printf(GP "set term postscript enhanced font \"Helvetica,20\"\n"); -printf(GP "set output \"histo.ps\"\n"); -printf(GP "set xtics nomirror out rotate by -60\n"); -#printf(GP "set logscale x\n"); -foreach $ref ( sort {-($counts{$a} <=> $counts{$b})} keys %counts ) -{ - my $max = $maxs{$ref}; - my $min = $mins{$ref}; - my $binsize = ($max+1-$min)/$bins; - my @hist; - - printf("%s %f -> %f\n", $ref, $mins{$ref}, $maxs{$ref}); - - for ( $i=0; $i<$bins; $i++ ) { - $hist[$i] = 0; - } - - seek(FH, 0, 0); - while ( $line = <FH> ) { - - my $bin; - my $h; - my $k; - my $l; - my $intensity; - - chomp($line); - - if ( $line =~ /^\s+(\d+)\s+(\d+)\s+(\d+)\s+([\d\.]+)/ ) { - - my $key; + if ( $line =~ /\s+(\d+)$/ ) { + my $n = $1; + if ( $n > $max ) { + $line =~ /^\s+(\d+)\s+(\d+)\s+(\d+)\s+/; $h = $1; $k = $2; $l = $3; - $intensity = $4; - - $key = $h.",".$k.",".$l; - - next unless ( $ref eq $key ); - - $bin = int(($intensity-$min)/$binsize); - $hist[$bin]++; - + $max = $n; } } - - open(OFH, "> ".$ref.".dat"); - for ( $i=0; $i<$bins; $i++ ) { - printf(OFH "%f\t%i\n", $min+(($i+0.5)*$binsize), $hist[$i]); - } - close(OFH); - - my $nref = $ref; - $nref =~ s/\,/\ /g; - my $title = "t \"".$nref." reflection\""; - printf(GP "plot [] [0:20]\"".$ref.".dat\" u 1:2 w histeps ".$title."\n"); - } -close(FH); -close(GP); +printf("%i %i %i = %i\n", $h, $k, $l, $max); -foreach $ref ( sort {-($counts{$a} <=> $counts{$b})} keys %counts ) { - unlink($ref.".dat"); -} -system("ps2pdf histo.ps"); -unlink("histo.ps"); +exit 0; diff --git a/src/process_hkl.c b/src/process_hkl.c index ec23e1cd..05c71985 100644 --- a/src/process_hkl.c +++ b/src/process_hkl.c @@ -29,8 +29,8 @@ #include "symmetry.h" -/* Number of divisions for R vs |q| graphs */ -#define RVQDV (20) +/* Number of divisions for intensity histograms */ +#define NBINS (200) static void show_help(const char *s) @@ -56,6 +56,8 @@ static void show_help(const char *s) " --stop-after=<n> Stop after processing n patterns. Zero means\n" " keep going until the end of the input, and is\n" " the default.\n" +" -g, --histogram=<h,k,l> Calculate the histogram of measurements for this\n" +" reflection.\n" "\n" " --scale Scale each pattern for best fit with the current\n" " model.\n" @@ -64,6 +66,47 @@ static void show_help(const char *s) } +static void plot_histogram(double *vals, int n) +{ + int i; + double max = -INFINITY; + double min = +INFINITY; + double step; + int histo[NBINS]; + FILE *fh; + + fh = fopen("histogram.dat", "w"); + if ( fh == NULL ) { + ERROR("Couldn't open 'histogram.dat'\n"); + return; + } + + for ( i=0; i<n; i++ ) { + if ( vals[i] > max ) max = vals[i]; + if ( vals[i] < min ) min = vals[i]; + } + STATUS("%f %f\n", min, max); + + for ( i=0; i<NBINS; i++ ) { + histo[i] = 0; + } + + step = (max-min)/NBINS; + + for ( i=0; i<n; i++ ) { + int bin; + bin = (vals[i]-min)/step; + histo[bin]++; + } + + for ( i=0; i<NBINS; i++ ) { + fprintf(fh, "%f %i\n", min+step*i, histo[i]); + } + + fclose(fh); +} + + /* Note "holo" needn't actually be a holohedral point group, if you want to try * something strange like resolving from a low-symmetry group into an even * lower symmetry one. @@ -160,7 +203,9 @@ static void merge_pattern(double *model, ReflItemList *observed, const double *new, ReflItemList *items, unsigned int *model_counts, int mo, ReflItemList *twins, - const char *holo, const char *mero) + const char *holo, const char *mero, double *hist_vals, + signed int hist_h, signed int hist_k, + signed int hist_l, int *hist_n) { int i; int twin; @@ -216,6 +261,14 @@ static void merge_pattern(double *model, ReflItemList *observed, /* Increase count count */ integrate_count(model_counts, h, k, l, 1); + if ( hist_vals != NULL ) { + int p = *hist_n; + if ( (h==hist_h) && (k==hist_k) && (l==hist_l) ) { + hist_vals[p] = intensity; + *hist_n = p+1; + } + } + } /* Dump the reflections in this pattern into the overall list */ @@ -230,7 +283,9 @@ static void merge_all(FILE *fh, double **pmodel, ReflItemList **pobserved, int config_maxonly, int config_scale, int config_sum, int config_stopafter, ReflItemList *twins, const char *holo, const char *mero, - int n_total_patterns) + int n_total_patterns, double *hist_vals, + signed int hist_h, signed int hist_k, signed int hist_l, + int *hist_i) { char *rval; float f0; @@ -277,7 +332,9 @@ static void merge_all(FILE *fh, double **pmodel, ReflItemList **pobserved, /* Start of second or later pattern */ merge_pattern(model, observed, new_pattern, items, counts, config_maxonly, - twins, holo, mero); + twins, holo, mero, + hist_vals, hist_h, hist_k, hist_l, + hist_i); if ( n_patterns == config_stopafter ) break; @@ -380,6 +437,10 @@ int main(int argc, char *argv[]) ReflItemList *observed; int i; const char *holo = NULL; + char *histo = NULL; + signed int hist_h, hist_k, hist_l; + double *hist_vals = NULL; + int hist_i; /* Long options */ const struct option longopts[] = { @@ -393,11 +454,12 @@ int main(int argc, char *argv[]) {"scale", 0, &config_scale, 1}, {"symmetry", 1, NULL, 'y'}, {"pdb", 1, NULL, 'p'}, + {"histogram", 1, NULL, 'g'}, {0, 0, NULL, 0} }; /* Short options */ - while ((c = getopt_long(argc, argv, "hi:e:ro:p:y:", + while ((c = getopt_long(argc, argv, "hi:e:ro:p:y:g:", longopts, NULL)) != -1) { switch (c) { @@ -425,6 +487,10 @@ int main(int argc, char *argv[]) sym = strdup(optarg); break; + case 'g' : + histo = strdup(optarg); + break; + case 0 : break; @@ -475,6 +541,22 @@ int main(int argc, char *argv[]) holo = strdup("1"); } + if ( histo != NULL ) { + int r; + r = sscanf(histo, "%i,%i,%i", &hist_h, &hist_k, &hist_l); + if ( r != 3 ) { + ERROR("Invalid indices for '--histogram'\n"); + return 1; + } + hist_vals = malloc(10*1024*sizeof(double)); + free(histo); + STATUS("Histogramming %i %i %i -> ", hist_h, hist_k, hist_l); + /* Put into the asymmetric cell for the target group */ + get_asymm(hist_h, hist_k, hist_l, + &hist_h, &hist_k, &hist_l, sym); + STATUS("%i %i %i\n", hist_h, hist_k, hist_l); + } + /* Open the data stream */ if ( strcmp(filename, "-") == 0 ) { fh = stdin; @@ -492,13 +574,21 @@ int main(int argc, char *argv[]) STATUS("There are %i patterns to process\n", n_total_patterns); rewind(fh); + hist_i = 0; merge_all(fh, &model, &observed, &counts, config_maxonly, config_scale, config_sum, config_stopafter, - twins, holo, sym, n_total_patterns); + twins, holo, sym, n_total_patterns, + hist_vals, hist_h, hist_k, hist_l, &hist_i); rewind(fh); fclose(fh); + 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); + } + if ( output != NULL ) { write_reflections(output, observed, model, NULL, counts, cell); } |