aboutsummaryrefslogtreecommitdiff
path: root/src/powder_plot.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/powder_plot.c')
-rw-r--r--src/powder_plot.c1204
1 files changed, 0 insertions, 1204 deletions
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;
-}
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-16 13:15:33 +0000