diff options
| -rw-r--r-- | .gitignore | 1 | ||||
| -rw-r--r-- | Makefile.am | 9 | ||||
| -rw-r--r-- | doc/man/powder_plot.1 | 159 | ||||
| -rw-r--r-- | src/powder_plot.c | 1204 |
4 files changed, 3 insertions, 1370 deletions
@@ -29,7 +29,6 @@ src/get_hkl src/hdfsee src/indexamajig src/compare_hkl -src/powder_plot src/render_hkl src/partialator src/check_hkl diff --git a/Makefile.am b/Makefile.am index 8feebadc..bf353e2e 100644 --- a/Makefile.am +++ b/Makefile.am @@ -1,10 +1,10 @@ +docdir = ${datadir}/doc/${PACKAGE} EXTRA_DIST = configure m4/gnulib-cache.m4 SUBDIRS = lib doc/reference libcrystfel ACLOCAL_AMFLAGS = -I m4 bin_PROGRAMS = src/pattern_sim src/process_hkl src/get_hkl src/indexamajig \ - src/compare_hkl src/powder_plot \ - src/partialator src/check_hkl src/partial_sim + src/compare_hkl src/partialator src/check_hkl src/partial_sim noinst_PROGRAMS = tests/list_check tests/integration_check \ tests/pr_gradient_check tests/symmetry_check \ @@ -60,8 +60,6 @@ src_compare_hkl_SOURCES = src/compare_hkl.c src_check_hkl_SOURCES = src/check_hkl.c -src_powder_plot_SOURCES = src/powder_plot.c - if HAVE_CAIRO src_render_hkl_SOURCES = src/render_hkl.c endif @@ -108,8 +106,7 @@ EXTRA_DIST += $(crystfel_DATA) man_MANS = doc/man/crystfel.7 doc/man/check_hkl.1 doc/man/compare_hkl.1 \ doc/man/crystfel_geometry.5 doc/man/get_hkl.1 doc/man/hdfsee.1 \ doc/man/indexamajig.1 doc/man/partialator.1 doc/man/partial_sim.1 \ - doc/man/pattern_sim.1 doc/man/powder_plot.1 doc/man/process_hkl.1 \ - doc/man/render_hkl.1 + doc/man/pattern_sim.1 doc/man/process_hkl.1 doc/man/render_hkl.1 EXTRA_DIST += $(man_MANS) diff --git a/doc/man/powder_plot.1 b/doc/man/powder_plot.1 deleted file mode 100644 index 6dac163f..00000000 --- a/doc/man/powder_plot.1 +++ /dev/null @@ -1,159 +0,0 @@ -.\" powder_plot man page -.\" -.\" Copyright © 2012 Andrew Aquila <andrew.aquila@cfel.de> -.\" Copyright © 2012 Thomas White <taw@physics.org> -.\" -.\" Part of CrystFEL - crystallography with a FEL -.\" - -.TH POWDER_PLOT 1 -.SH NAME -powder_plot \- generate 1D powder patterns -.SH SYNOPSIS -.PP -.B powder_plot -\fB-i\fR \fIinput.hkl\fR | \fB-i\fR \fIinput.h5\fR | \fB-i\fR \fIinput.stream\fR -\fB-o\fR \fIoutput.dat\fR -[\fB--min=\fR\fI1/d\fR \fB--max=\fR\fI1/d\fR] -[\fIoptions\fR\] \fB...\fR - -.PP -.BR powder_plot -\fB-i\fR \fIfile.hkl\fR [\fIoptions\fR\] \fB...\fR [\fB--use-redundancy\fR] [\fB--no-d-scaling\fR] - -.PP -.BR powder_plot -\fB-i\fR \fIfile.h5\fR [\fIoptions\fR\] \fB...\fR -\fB-g\fR \fIgeometry.geom\fR -\fB-b\fR \fIbeam.beam\fR [\fB--no-sat-corr\fR] [\fB--ring-corr\fR] - -.PP -.BR powder_plot -\fB-i\fR \fIfile.stream\fR [\fIoptions\fR\] \fB...\fR --data=\fItype\fR -[-g \fIgeometry.geom\fR] [-b \fIbeam.beam\fR] [\fB--no-sat-corr\fR] [\fB--only-indexed\fR] -[\fB--use-redundancy\fR] [\fB--ring-corr\fR] [\fB--no-d-scaling\fR] - -.PP -.BR powder_plot -\fB--help\fR - -.SH DESCRIPTION - -powder_plot calculates one dimensional powder traces by summing many individual intensities, Bragg peaks or pixels. Its input can come from a CrystFEL stream (such as that written by "indexamajig"), an reflection list in CrystFEL format (".hkl" format), or an HDF5 file. - -The output of powder_plot consists of a three line header followed by a tab-delimited list of six values: - - 1/d of the histogram bin, where d is the Bragg law d spacing in meters - the total number of peaks (N) - the total intensity in the N peaks - the mean intensity of the N peaks - the standard deviation of the distribution - the standard deviation of the mean of the data - -The sigma of the mean is not the same as the sigma of the intensities -themselves. The former quantity measures how accurately the mean intensity has -been determined, whereas the latter quantity measures the spread of the -intensities. - -.SH DATA TYPE OPTIONS WHEN READING FROM A STREAM - -When taking input from stream, the d-spacing for a particular intensity can be -generated in a variety of different ways. You can choose which one with -\fB--data=\fR\fItype\fR. Possibilities for \fItype\fR are: -.RS -.IP \fBreflection\fR -.PD -Use peak positions from indexed reflections -.IP \fBhkl\fR -.PD -Use the Miller indices from indexed reflections, combined with a unit cell from PDB file provided with -p. -.IP \fBpeaks\fR -.PD -Use peak positions from peak search -.IP \fBh5\fR -.PD -Use all pixels in the HDF5 file, excluding bad regions -.RE -.PP -The default is \fB--data=hkl\fR. - - -.SH HISTOGRAM OPTIONS - -You can set the mininum and maximum 1/d values, in units of inverse meters, -with the options \fB-min=\fR\fIn\fR and \fB--max=\fR\fIn\fR. -When taking input from peak positions, The default behaviour is to use the entire detector extent from the geometry description file, which you with the \fB-g\fR flag. - -You can also adjust the number of histogram bins with the option --bins=<n>, -where n is an integer. - -Scaling can be set to produce linearly, quadratically or cubically spaced 1/d -values using \fB--spacing=\fR\fItype\fR. Possibilities for \fItype\fR are: -.RS -.IP \fBlinear\fR -.PD -The resolution shells will have equal widths in terms of 1/d. -.IP \fBwilson\fR -.PD -The resolution shells will be quadratically spaced, giving even spacing in a plot against 1/d^2 (a Wilson plot). -.IP \fBvolume\fR -.PD -The resolution shells will all have equal volumes in reciprocal space. -.RE -.PP -The default is \fB--spacing=linear\fR. - -.SH OPTIONS FOR MORE CONTROL - -.B -.IP --no-sat-corr -Don't correct values of saturated peaks using the table included in the HDF5 file. -See the help for indexamajig for more information. - -.B -.IP --only-indexed -Use with --data=peaks or h5 if you want to use the peak list of only indexed patterns. -This is useful for finding differences between patterns which could be indexed and -those which could not. - -.B -.IP --no-d-scaling -Do not scale the intensities in the powder plot by d^2. This should be used when -creating a powder plot from a reflection list. - -.B -.IP --ring-corr -Do not correct for the fractional area sampled of the powder ring. This might be -useful for detectors with gaps. - -.B -.IP --use-redundancy -Divide intensities by the number of measurements (the redundancy column in the -reflection list), and not the number of symmetrical equivalent reflections as the -number of times a reflection occurs in the powder. - -.SH AUTHOR -This page was written by Andrew Aquila and Thomas White. - -.SH REPORTING BUGS -Report bugs to <taw@physics.org>, or visit <http://www.desy.de/~twhite/crystfel>. - -.SH COPYRIGHT AND DISCLAIMER -Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, a research centre of the Helmholtz Association. -.PD -.P -powder_plot, and this manual, are part of CrystFEL. -.P -CrystFEL is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version. -.P -CrystFEL is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. -.P -You should have received a copy of the GNU General Public License along with CrystFEL. If not, see <http://www.gnu.org/licenses/>. - -.SH SEE ALSO -.BR crystfel (7), -.BR indexamajig (1), -.BR process_hkl (1), -.BR check_hkl (1), -.BR crystfel_geometry (5), -.BR render_hkl (1) . diff --git a/src/powder_plot.c b/src/powder_plot.c deleted file mode 100644 index 7917ccd0..00000000 --- a/src/powder_plot.c +++ /dev/null @@ -1,1204 +0,0 @@ -/* - * powder_plot.c - * - * Plot powder patterns - * - * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. - * Copyright © 2012 Richard Kirian - * - * Authors: - * 2010-2012 Thomas White <taw@physics.org> - * 2011-2012 Andrew Aquila <andrew.aquila@cfel.de> - * 2011 Richard Kirian - * - * This file is part of CrystFEL. - * - * CrystFEL is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * CrystFEL is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>. - * - */ - - -#ifdef HAVE_CONFIG_H -#include <config.h> -#endif - -#include <stdarg.h> -#include <stdlib.h> -#include <stdio.h> -#include <string.h> -#include <unistd.h> -#include <getopt.h> - -#include "stream.h" -#include "reflist.h" -#include "utils.h" -#include "image.h" -#include "detector.h" -#include "index.h" -#include "hdf5-file.h" -#include "beam-parameters.h" -#include "reflist-utils.h" -#include "symmetry.h" -#include "cell-utils.h" - - -struct bin_stats { - unsigned int N; - double total; - double mean; - double std_dev; - double q_min; - double q_max; - double q_value; - double fract; -}; - -struct histogram_info { - double q_max; - double q_min; - double q_delta; - unsigned int histsize; - int spacing; //linear, q^2, & equal volume -}; - -enum { - PLOT_PEAKS, - PLOT_HKL, - PLOT_REFL, - PLOT_H5 -}; - -enum { - FILE_STREAM, - FILE_HKL, - FILE_H5 -}; - -enum { - LINEAR, - q2, - VOLUME -}; - - -static int find_q_bin_index(double q, struct histogram_info *info, - struct bin_stats *hist) -{ - /* bisection search alg. find q_bin index of order Log(n) time */ - int mid; - int min = 0; - int max = info->histsize-1; - if (q < hist[min].q_max) return min; - if (q > hist[max].q_min) return max; - do { - mid = (min + max) / 2; - if (q < hist[mid].q_min) { - max = mid; - } else if (q > hist[mid].q_max){ - min = mid ; - } else { - return mid; - } - } while(max - min > 1); - return mid; -} - - -/* Used for HDF5 files, peak list and stream positions */ -static int add_peak_to_histogram(double fs, double ss, double intensity, - struct image *image, - struct histogram_info *info, - struct bin_stats *hist) -{ - struct rvec r; - double q, delta; - int i; - - r = get_q(image, fs, ss, NULL, 1.0/ image->lambda); - q = modulus(r.u, r.v, r.w); - - /* Ignore q value if outside of range */ - if ( (q<info->q_min) || (q>info->q_max) ) return 1; - - i = find_q_bin_index(q, info, hist); - - /* See Knuth TAOCP vol 2, 3rd ed, pg 232 for running variance */ - delta = intensity - hist[i].mean; - hist[i].N++; - hist[i].total += intensity; - hist[i].mean = hist[i].mean + delta /hist[i].N; - hist[i].std_dev = hist[i].std_dev + (delta *(intensity - hist[i].mean)); - - return 0; -} - - -/* Used for d and hkl of stream files where redundancy = 1 */ -static int add_d_to_histogram(double q, double intensity, - struct histogram_info *info, - struct bin_stats *hist) -{ - double delta; - int i; - - /* Ignore q value if outside of range */ - if ( (q<info->q_min) || (q>info->q_max) ) return 1; - - i = find_q_bin_index(q, info, hist); - - delta = intensity - hist[i].mean; - hist[i].N++; - hist[i].total += intensity; - hist[i].mean = hist[i].mean + delta /hist[i].N; - hist[i].std_dev = hist[i].std_dev + (delta *(intensity - hist[i].mean)); - - return 0; -} - - -static int add_hkl_to_histogram(double q, double intensity, int redundancy, - int q_scaling, struct histogram_info *info, - struct bin_stats *hist) -{ - int i = 0; - - /* Ignore q value if outside of range */ - if ( (q<info->q_min) || (q>info->q_max) ) { - return 1; - } - - /* The accounting is the intensity of the reflection times the - * number of occurance of that reflection smeared out over the - * surface area which is 4*pi*q^2 the 4*pi is left out since it is a - * common constant and the total is in arbitrary units. - */ - for ( i=0; i<redundancy; i++ ) { - if ( q_scaling ) { - add_d_to_histogram(q, intensity/(q*q), info, hist); - } else { - add_d_to_histogram(q, intensity, info, hist); - } - } - - return 0; -} - - -static int histogram_setup(struct histogram_info *info, - struct bin_stats *histdata) -{ - int i; - double x; - - if ( info->spacing == LINEAR ) { - x = 1.0; - } else if ( info->spacing == q2 ) { - x = 2.0; - } else { - x = 3.0; - } - - for ( i=0; i<info->histsize; i++ ) { - - double qd, qm; - - histdata[i].N = 0; - histdata[i].total = 0.0; - histdata[i].mean = 0.0; - histdata[i].std_dev = 0.0; - histdata[i].fract = 0.0; - - qd = info->q_delta; - qm = info->q_min; - - histdata[i].q_min = pow( (i*qd) + pow(qm, x), 1.0/x); - - histdata[i].q_max = pow( ((i+1.0)*qd) + pow(qm, x), 1.0/x); - - histdata[i].q_value= pow( ((i+0.5)*qd) + pow(qm, x), 1.0/x); - - } - return 0; -} - - -static int ring_fraction_calc(struct histogram_info *info, - struct bin_stats *hist, struct image *image) -{ - int fs,ss; - int bin; - - /* Check that detector geometry is present and wavelength is valid */ - if ( (image->det == NULL) || (image->lambda < 0.0) ) return 1; - - /* Loop over all pixels */ - for ( ss=0; ss<image->height; ss++ ) { - for ( fs=0; fs<image->width; fs++ ) { - - struct panel *p; - struct rvec r; - int i; - double q, q_fs, q_ss; - - r = get_q(image, fs, ss, NULL, 1.0/image->lambda); - q = modulus(r.u, r.v, r.w); - - /* If pixel is valid (not a bad pixel and not out of range) */ - if ( (q>info->q_min) && (q<info->q_max) && - (in_bad_region(image->det,fs,ss) == 0) ) { - - /* Select the panel, then (sometimes) ask for the q - * of the (corner of the) pixel one step beyond the - * edge, to get the exact size of the required pixel. - */ - p = find_panel(image->det, fs, ss); - - r = get_q_for_panel(p, (fs+1)-(double)p->min_fs, - ss-(double)p->min_ss, - NULL, 1.0/image->lambda); - q_fs = modulus(r.u, r.v, r.w); - - r = get_q_for_panel(p, fs-(double)p->min_fs, - (ss+1) -(double)p->min_ss, - NULL, 1.0/image->lambda); - q_ss = modulus(r.u, r.v, r.w); - - i = find_q_bin_index(q, info, hist); - - hist[i].fract = hist[i].fract + fabs((q_fs-q)*(q_ss-q)); - - } - - } - } - - /* Divide measured area by ring area */ - for ( bin=0; bin<info->histsize; bin++ ) { - - double inner_area, outer_area, ring_area; - - outer_area = pow(hist[bin].q_max, 2.0); - inner_area = pow(hist[bin].q_min, 2.0); - ring_area = M_PI*(outer_area - inner_area); - - hist[bin].fract = hist[bin].fract / ring_area; - - } - - return 0; -} - - -static unsigned int process_h5(struct image *image, struct histogram_info *info, - struct bin_stats *histdata) -{ - int fs, ss; - double intensity; - - for ( ss=0; ss<image->height; ss++ ) { - for ( fs=0; fs<image->width; fs++ ) { - - intensity = image->data[fs + image->width*ss]; - if ( !in_bad_region(image->det,fs,ss) ) { - add_peak_to_histogram(fs, ss, intensity, - image, info, histdata); - } - } - progress_bar(ss, image->height, "Processing"); - } - return 0; -} - - -static unsigned int process_hkl(struct image *image, const SymOpList *sym, - UnitCell *cell, - struct histogram_info *info, - struct bin_stats *histdata, - int q_scaling, int use_redundancy) -{ - Reflection *refl; - RefListIterator *iter; - unsigned int i = 0; - unsigned int n_peaks = 0; - int h, k, l, redundancy; - double q, intensity; - unsigned int nref; - SymOpMask *m; - - m = new_symopmask(sym); - - nref = num_reflections(image->reflections); - - for ( refl = first_refl(image->reflections, &iter); - refl != NULL; - refl = next_refl(refl, iter) ) - { - get_indices(refl, &h, &k, &l); - intensity = get_intensity(refl); - if ( use_redundancy ) { - redundancy = get_redundancy(refl); - } else { - special_position(sym, m, h, k, l); - redundancy = num_equivs(sym, m); - } - - /* Multiply by 2 to get 1/d (in m^-1) */ - q = 2.0 * resolution(cell, h, k, l); - - add_hkl_to_histogram(q, intensity, redundancy, q_scaling, - info, histdata); - - n_peaks += redundancy; - - i++; - progress_bar(i, nref, "Processing"); - - } - - free_symopmask(m); - - return n_peaks; -} - - -static unsigned int process_stream_reflection(FILE *fh, struct image *image, - struct histogram_info *info, - struct bin_stats *histdata, - unsigned int *n_patterns) -{ - int rval; - unsigned int i = 0; - unsigned int n_peaks = 0; - Reflection *refl; - RefListIterator *iter; - double intensity, fs_double, ss_double; - unsigned int processing_total; - - processing_total = count_patterns(fh); - rewind(fh); - - do { - /* Get data from next chunk */ - rval = read_chunk(fh, image); - if ( rval ) continue; - - /* Check if the pattern indexed, if so use those peaks */ - if ( image->reflections != NULL ) { - - (*n_patterns)++; - for ( refl = first_refl(image->reflections, &iter); - refl != NULL; - refl = next_refl(refl, iter) ) { - /* note added fs_double as fs is an int */ - intensity = get_intensity(refl); - get_detector_pos(refl, &fs_double, &ss_double); - - if ( !add_peak_to_histogram(fs_double, - ss_double, - intensity, - image, info, - histdata) ) - { - n_peaks++; - } - - } - - } - - free(image->filename); - reflist_free(image->reflections); - image_feature_list_free(image->features); - cell_free(image->indexed_cell); - - i++; - progress_bar(i, processing_total, "Processing"); - - } while ( rval == 0 ); - - return n_peaks; -} - - -static unsigned int process_stream_hkl(FILE *fh, struct image *image, - struct histogram_info *info, - struct bin_stats *histdata, - UnitCell *cell, unsigned int *n_patterns) -{ - int rval; - unsigned int i = 0; - unsigned int n_peaks = 0; - Reflection *refl; - RefListIterator *iter; - double intensity, q; - unsigned int processing_total; - - processing_total = count_patterns(fh); - rewind(fh); - - do { - - /* Get data from next chunk */ - rval = read_chunk(fh, image); - if ( rval ) continue; - if ( image->reflections != NULL ) { - - (*n_patterns)++; - - for ( refl = first_refl(image->reflections, &iter); - refl != NULL; - refl = next_refl(refl, iter) ) - { - int h, k, l; - - get_indices(refl, &h, &k, &l); - intensity = get_intensity(refl); - q = 2.0 * resolution(cell, h, k, l); - - if ( !add_d_to_histogram(q, intensity, info, - histdata) ) n_peaks++; - } - } - - free(image->filename); - reflist_free(image->reflections); - image_feature_list_free(image->features); - cell_free(image->indexed_cell); - - i++; - progress_bar(i, processing_total, "Processing"); - - } while ( rval == 0 ); - - return n_peaks; -} - - -static int add_features_to_histogram(struct image *image, - struct histogram_info *info, - struct bin_stats *histdata) -{ - int j; - int n_peaks; - - n_peaks = 0; - for ( j=0; j<image_feature_count(image->features); j++) { - - struct imagefeature *f; - - f = image_get_feature(image->features, j); - if ( !f->valid ) continue; - - if ( !in_bad_region(image->det, f->fs,f->ss) ) { - - int r; - - r = add_peak_to_histogram(f->fs, f->ss, f->intensity, - image, info, histdata); - - if ( !r ) n_peaks++; - - } - - } - - return n_peaks; -} - - -static unsigned int process_stream_peaks(FILE *fh, struct image *image, - struct histogram_info *info, - struct bin_stats *histdata, - unsigned int *n_patterns, - int only_indexed) -{ - int rval; - unsigned int i = 0; - unsigned int n_peaks = 0; - unsigned int processing_total; - - processing_total = count_patterns(fh); - rewind(fh); - - do { - - /* Get data from next chunk */ - rval = read_chunk(fh, image); - if ( rval ) continue; - - if ( image->features != NULL ) { - - if ( (!only_indexed) - || ( only_indexed && (image->reflections != NULL)) ) - { - (*n_patterns)++; - n_peaks += add_features_to_histogram(image, - info, - histdata); - } - } - - free(image->filename); - reflist_free(image->reflections); - image_feature_list_free(image->features); - cell_free(image->indexed_cell); - - i++; - progress_bar(i, processing_total, "Processing"); - - } while ( rval == 0 ); - - return n_peaks; -} - - -static unsigned int process_stream_h5(FILE *fh, struct image *image, - struct histogram_info *info, - struct bin_stats *histdata, - int config_satcorr, int only_indexed, - unsigned int *n_patterns, - const char *element) -{ - int fs, ss, rval; - double intensity; - unsigned int i = 0; - unsigned int n_peaks = 0; - struct hdfile *hdfile = NULL; - unsigned int processing_total; - - processing_total = count_patterns(fh); - rewind(fh); - - do { - - /* Get data from next chunk */ - rval = read_chunk(fh, image); - if ( rval ) continue; - if ( !only_indexed || - ( only_indexed && (image->reflections != NULL) ) ) - { - hdfile = hdfile_open(image->filename); - if ( element != NULL ) { - int r; - r = hdfile_set_image(hdfile, element); - if ( r ) { - ERROR("Couldn't select path '%s'\n", - element); - hdfile_close(hdfile); - return 0; - } - } else { - int r; - r = hdfile_set_first_image(hdfile, "/"); - if ( r ) { - ERROR("Couldn't select first path\n"); - hdfile_close(hdfile); - return 0; - } - - } - hdf5_read(hdfile, image, config_satcorr); - hdfile_close(hdfile); - n_patterns++; - for ( ss=0; ss<image->height; ss++ ) { - for ( fs=0; fs<image->width; fs++ ) { - - intensity = image->data[fs + image->width*ss]; - if ( !in_bad_region(image->det,fs,ss) ) { - - add_peak_to_histogram(fs, ss, intensity, - image, info, - histdata); - - } - - } - } - } - - free(image->data); - free(image->filename); - reflist_free(image->reflections); - image_feature_list_free(image->features); - cell_free(image->indexed_cell); - - i++; - progress_bar(i, processing_total, "Processing"); - - } while ( rval == 0 ); - - return n_peaks; -} - - -static void show_help(const char *s) -{ - printf("Syntax: %s [options]\n\n", s); - printf( -"Plot a powder pattern as a 1D graph using the detector geometry.\n" -"\n" -" -h, --help Display this help message.\n" -" -i, --input=<file> Input filename. (*.stream, *.hkl, or *.h5)\n" -" -o, --output=<file> Output filename. Default: stdout.\n" -" -g. --geometry=<file> Get detector geometry from <file>.\n" -" -b, --beam=<file> Get beam parameters (wavelength) from <file>.\n" -" -p, --pdb=<file> Get unit cell from PDB file. (.hkl files only)\n" -" -y, --symmetry=<sym> The symmetry of crystal (.hkl files only)\n" -" -s, --bins=n Makes histogram with n bins (default is 100).\n" -" --spacing=<type> Use 'type' to select the 1/d spacing.\n" -" Choose from:\n" -" linear : linear (default)\n" -" wilson : even spacing in Wilson plots\n" -" volume : constant volume\n" -" --max=n The maximum 1/d to be considered in plot.\n" -" --min=n The minimum 1/d to be considered in plot.\n" -" -d, --data=<type> Use to select the kind of stream data in histogram.\n" -" Choose from:\n" -" reflection : uses peak positons from indexed\n" -" reflection \n" -" hkl : uses the hkl list from indexed\n" -" reflections (requires pdb file)\n" -" peaks : uses all peaks found from peak\n" -" search\n" -" h5 : all points in h5, excluding bad\n" -" regions\n" -" The default is 'hkl'.\n" -" --no-sat-corr Don't correct values of saturated peaks using a\n" -" table included in the HDF5 file.\n" -" --only-indexed Use with -data=peaks or h5 if you want to use the\n" -" peak list of only indexed patterns\n" -" --no-d-scaling Use with .hkl files if you want to not scale the\n" -" powder by d^2\n" -" --ring-corr Use if you want to scale the powder plot to\n" -" correct for the fractional area sampled of the\n" -" powder ring\n" -" --use-redundancy Use with .hkl files if you want to use the number\n" -" of measurements and not the number of symetrical\n" -" equivelent reflections as the number of time a\n" -" reflection occurs in the powder\n" -" -e, --image=<element> Use this image when reading an HDF5 file.\n" -" Example: /data/data0.\n" -" Default: The first one found.\n" -"\n"); -} - - -static void rlim_bailout() -{ - ERROR("Unable to automatically determine the resolution limits.\n"); - ERROR("Try again with --min and --max.\n"); - exit(1); -} - - -int main(int argc, char *argv[]) -{ - FILE *fh = NULL; - UnitCell *cell = NULL; - struct image image; - struct hdfile *hdfile = NULL; - struct bin_stats *histdata = NULL; - struct histogram_info hist_info; - SymOpList *sym; - - /* Default settings */ - hist_info.histsize = 100; - hist_info.q_min = -1.0; - hist_info.q_max = -1.0; - hist_info.spacing = LINEAR; - image.lambda = -1.0; - image.beam = NULL; - image.det = NULL; - char *element = NULL; - - unsigned int n_patterns = 0; - unsigned int n_peaks = 0; - - int c, file_type, data_type; - int config_satcorr = 1; - int need_geometry = 0; - int need_beam = 0; - int need_pdb = 0; - int only_indexed = 0; - int q_scaling = 1; - int ring_corr = 0; - int use_redundancy = 0; - unsigned int i; - - char *filename = NULL; - char *geometry = NULL; - char *beamf = NULL; - char *pdb = NULL; - char *output = NULL; - char *datatype = NULL; - char *sym_str = NULL; - - /* Long options */ - const struct option longopts[] = { - {"help", 0, NULL, 'h'}, - {"input", 1, NULL, 'i'}, - {"output", 1, NULL, 'o'}, - {"geometry", 1, NULL, 'g'}, - {"beam", 1, NULL, 'b'}, - {"pdb", 1, NULL, 'p'}, - {"symmetry", 1, NULL, 'y'}, - {"bins", 1, NULL, 's'}, - {"max", 1, NULL, 1 }, - {"min", 1, NULL, 2 }, - {"spacing", 1, NULL, 3 }, - {"no-sat-corr", 0, &config_satcorr, 0 }, - {"sat-corr", 0, &config_satcorr, 1 }, - {"only-indexed", 0, &only_indexed, 1 }, - {"no-d-scaling", 0, &q_scaling, 0 }, - {"ring-corr", 0, &ring_corr, 1 }, - {"use-redundancy", 0, &use_redundancy, 1 }, - {"data", 1, NULL, 'd'}, - {"image", 1, NULL, 'e'}, - {0, 0, NULL, 0} - }; - - /* Short options */ - while ((c = getopt_long(argc, argv, "hi:o:g:b:p:y:s:e:d:", - longopts, NULL)) != -1) - { - - switch (c) { - - case 'h' : - show_help(argv[0]); - return 0; - - case 'i' : - filename = strdup(optarg); - break; - - case 'o' : - output = strdup(optarg); - break; - - case 'g' : - geometry = strdup(optarg); - break; - - case 'b' : - beamf = strdup(optarg); - break; - - case 'p' : - pdb = strdup(optarg); - break; - - case 'y' : - sym_str = strdup(optarg); - break; - - case 's' : - hist_info.histsize = atoi(optarg); - break; - - case 'e' : - element = strdup(optarg); - break; - - case 1 : - hist_info.q_max = atof(optarg); - break; - - case 2 : - hist_info.q_min = atof(optarg); - break; - - case 3 : - if (strcmp(optarg, "linear") == 0 ) { - hist_info.spacing = LINEAR; - } else if (strcmp(optarg, "wilson") == 0 ) { - hist_info.spacing = q2; - } else if (strcmp(optarg, "volume") == 0) { - hist_info.spacing = VOLUME; - } else { - ERROR("Invalid spacing type: '%s'\n", optarg); - return 1; - } - break; - - case 'd' : - datatype = strdup(optarg); - break; - - case 0 : - break; - - default : - ERROR("Unhandled option '%c'\n", c); - break; - } - - } - - /* Process input file type */ - if ( filename == NULL ) { - - ERROR("You must specify the input filename with -i\n"); - return 1; - - } - - if ( is_stream(filename) == 1 ) { - - file_type = FILE_STREAM; - - } else if ( H5Fis_hdf5(filename) > 0 ) { - - file_type = FILE_H5; - need_geometry = 1; - - } else { - - image.reflections = read_reflections(filename); - - if ( image.reflections != NULL ) { - file_type = FILE_HKL; - need_pdb = 1; - need_geometry = 0; - need_beam = 0; - image.lambda = 0.0; - } else { - ERROR("Couldn't recognise %s as reflection list," - " stream or image.\n", filename); - return 1; - } - - } - - if ( datatype == NULL ) { - data_type = PLOT_HKL; - if ((hist_info.q_min < 0.0) || (hist_info.q_max < 0.0)) { - need_geometry = 1; - } - need_pdb = 1; - - } else if ( strcmp(datatype, "reflection") == 0 ) { - data_type = PLOT_REFL; - need_geometry = 1; - - } else if ( strcmp(datatype, "hkl") == 0 ) { - data_type = PLOT_HKL; - need_pdb = 1; - - } else if ( strcmp(datatype, "peaks") == 0 ) { - data_type = PLOT_PEAKS; - need_geometry = 1; - - } else if ( strcmp(datatype, "h5") == 0 ) { - data_type = PLOT_H5; - need_geometry = 1; - - } else { - - ERROR("Failed to read data plot type: '%s'\n", datatype); - return 1; - } - - /* doubt this is needed, but double check just in case */ - if ( file_type == FILE_HKL ) { - need_geometry = 0; - need_beam = 0; - } - - /* Logic checks */ - if ( need_geometry && (image.lambda < 0.0) ) { - need_beam = 1; - } - if ( hist_info.histsize <= 0 ) { - ERROR("You need to specify a histogram with more then 0 " - "bins\n"); - return 1; - } - - /* Get geometry, beam and pdb files and parameters as needed */ - if ( need_geometry ) { - if ( geometry == NULL ) { - ERROR("You need to specify a geometry file with " - "--geometry\n"); - return 1; - } else { - image.det = get_detector_geometry(geometry); - if ( image.det == NULL ) { - ERROR("Failed to read detector geometry " - "from '%s'\n", geometry); - return 1; - } - image.width = image.det->max_fs; - image.height = image.det->max_ss; - } - } - free(geometry); - - /* Open files to get wavelength if it exists & camera length - if they are not found in the geometry file */ - if (file_type == FILE_STREAM) { - fh = fopen(filename, "r"); - if ( fh == NULL ) { - ERROR("Failed to open input file\n"); - return 1; - } - /* Use wavelength from first chunk */ - read_chunk(fh, &image); - rewind(fh); - } else if (file_type == FILE_H5) { - hdfile = hdfile_open(filename); - if ( element != NULL ) { - int r; - r = hdfile_set_image(hdfile, element); - if ( r ) { - ERROR("Couldn't select path '%s'\n", - element); - hdfile_close(hdfile); - return 0; - } - } else { - int r; - r = hdfile_set_first_image(hdfile, "/"); - if ( r ) { - ERROR("Couldn't select first path\n"); - hdfile_close(hdfile); - return 0; - } - - } - hdf5_read(hdfile, &image, config_satcorr); - hdfile_close(hdfile); - } - free(filename); - - /* Logic checks */ - if ( need_geometry && (image.lambda < 0.0) ) { - need_beam = 1; - } - if ( hist_info.histsize <= 0 ) { - ERROR("You need to specify a histogram with more then 0 " - "bins\n"); - return 1; - } - - if ( need_beam ) { - - if ( beamf == NULL ) { - ERROR("No wavelength in file, so you need to specify " - "a beam parameters file with --beam\n"); - return 1; - } else { - image.beam = get_beam_parameters(beamf); - if ( image.beam == NULL ) { - ERROR("Failed to read beam from '%s'\n", - beamf); - return 1; - } - image.lambda = ph_en_to_lambda(eV_to_J( - image.beam->photon_energy)); - } - - } - free(beamf); - - if ( need_pdb ) { - - if (pdb == NULL) { - ERROR("You need to specify a pdb file with --pdb.\n"); - return 1; - } else { - cell = load_cell_from_pdb(pdb); - if ( cell == NULL ) { - ERROR("Couldn't read unit cell (from %s)\n", - pdb); - return 1; - } - } - } - free(pdb); - - if ( sym_str == NULL ) { - sym_str = strdup("1"); - } - sym = get_pointgroup(sym_str); - free(sym_str); - - /* Set up histogram info*/ - if ( file_type == FILE_HKL || data_type == PLOT_HKL) { - /* get q range from Miller indices in hkl - file. */ - if ((hist_info.q_min < 0.0) && (hist_info.q_max < 0.0)) { - if ( image.reflections == NULL ) rlim_bailout(); - resolution_limits(image.reflections, cell, - &hist_info.q_min, &hist_info.q_max); - } else if (hist_info.q_min < 0.0) { - double dummy; - if ( image.reflections == NULL ) rlim_bailout(); - resolution_limits(image.reflections, cell, - &hist_info.q_min, &dummy); - } else if (hist_info.q_max < 0.0) { - double dummy; - if ( image.reflections == NULL ) rlim_bailout(); - resolution_limits(image.reflections, cell, - &dummy, &hist_info.q_max); - } - } else { - if ( hist_info.q_min < 0.0 ) { - hist_info.q_min = smallest_q(&image); - } - if ( hist_info.q_max < 0.0 ) { - hist_info.q_max = largest_q(&image); - } - } - - if ( hist_info.q_min >= hist_info.q_max ) { - ERROR("the minimum 1/d value (%e) " - "is greater then your max 1/d value (%e).\n", - hist_info.q_min, hist_info.q_max); - return 1; - } - - if ( hist_info.spacing == LINEAR) { - hist_info.q_delta = (hist_info.q_max - hist_info.q_min)/ - hist_info.histsize; - } else if ( hist_info.spacing == q2) { - hist_info.q_delta = (pow(hist_info.q_max, 2.0) - - pow(hist_info.q_min, 2.0)) / - hist_info.histsize; - } else { //by default must be in VOLUME - hist_info.q_delta = (pow(hist_info.q_max, 3.0) - - pow(hist_info.q_min, 3.0)) / - hist_info.histsize; - } - /* Set up histogram data */ - histdata = malloc((hist_info.histsize) * sizeof(struct bin_stats)); - histogram_setup(&hist_info, histdata); - - /* Set up ring scaling */ - if ( ring_corr ) { - - if ( ring_fraction_calc(&hist_info, histdata, &image) ) { - - ERROR("Detector geometry is required to correct for" - " finite ring area.\n"); - return 1; - } - } - - /* Process reflections based on file type and data type */ - switch ( file_type ) { - case FILE_H5 : - n_patterns++; - n_peaks = process_h5(&image, &hist_info, histdata); - free(image.data); - break; - - case FILE_HKL : - n_patterns++; //inc number of patterns used - n_peaks = process_hkl(&image, sym, cell, &hist_info, histdata, - q_scaling, use_redundancy); - break; - - case FILE_STREAM : - switch ( data_type ) { - - case PLOT_REFL : - n_peaks = process_stream_reflection(fh, &image, - &hist_info, histdata, &n_patterns); - break; - - case PLOT_HKL : - n_peaks = process_stream_hkl(fh, &image, &hist_info, - - histdata, cell, &n_patterns); - break; - case PLOT_PEAKS : - n_peaks = process_stream_peaks(fh, &image, &hist_info, - histdata, &n_patterns, only_indexed); - break; - - case PLOT_H5 : - n_peaks = process_stream_h5(fh, &image, &hist_info, - histdata, config_satcorr, only_indexed, - &n_patterns, element); - break; - - default : - break; - } - fclose(fh); - break; - - default : - break; - - } - - /* Sqrt the variance to get the std_dev */ - for( i=0; i<hist_info.histsize; i++ ) { - if (histdata[i].N > 1) { - histdata[i].std_dev = sqrt(histdata[i].std_dev/ - (histdata[i].N-1)); - } - } - if ( ring_corr ) { - for( i=0; i<hist_info.histsize; i++ ) { - histdata[i].N = histdata[i].N / - histdata[i].fract; - histdata[i].total = histdata[i].total / - histdata[i].fract; - histdata[i].mean = histdata[i].mean / - histdata[i].fract; - histdata[i].std_dev = histdata[i].total / - histdata[i].fract; - } - } - - /* Print out the results */ - if ( output != NULL ) { - fh = fopen(output, "w"); - if ( fh == NULL ) { - ERROR("Failed to open output file\n"); - return 1; - } - } else { - fh = stdout; - } - - /* Print header */ - fprintf(fh, "Command line:"); - for ( i=0; i<argc; i++ ) { - fprintf(fh, " %s", argv[i]); - } - fprintf(fh, "\n"); - fprintf(fh, "I read %i patterns with %i peaks\n", n_patterns, n_peaks); - fprintf(fh, "1/d(m^-1)\tN\ttotal\tmean\tstd dev\t std dev of mean\n"); - - for( i=0; i<hist_info.histsize; i++ ) { - fprintf(fh, "%5e\t%i\t%5e\t%5e\t%5e\t%5e\n", - histdata[i].q_min, histdata[i].N, - histdata[i].total, histdata[i].mean, - histdata[i].std_dev, - histdata[i].std_dev/sqrt(histdata[i].N)); - } - - if ( cell != NULL ) cell_free(cell); - if ( image.det != NULL ) free(image.det); - if ( image.beam != NULL ) free(image.beam); - fclose(fh); - free(histdata); - - return 0; -} |