diff options
Diffstat (limited to 'src/powder_plot.c')
-rw-r--r-- | src/powder_plot.c | 938 |
1 files changed, 549 insertions, 389 deletions
diff --git a/src/powder_plot.c b/src/powder_plot.c index 515c8c4e..13875364 100644 --- a/src/powder_plot.c +++ b/src/powder_plot.c @@ -4,7 +4,7 @@ * Plot powder patterns * * (c) 2011 Andrew Aquila <andrew.aquila@cfel.de> - * (c) 2006-2010 Thomas White <taw@physics.org> + * (c) 2006-2011 Thomas White <taw@physics.org> * * Part of CrystFEL - crystallography with a FEL * @@ -33,8 +33,9 @@ #include "reflist-utils.h" #include "symmetry.h" + struct bin_stats { - unsigned int N; + unsigned int N; double total; double mean; double std_dev; @@ -72,10 +73,11 @@ enum { VOLUME }; -static int find_q_bin_index(double q, struct histogram_info *info, - struct bin_stats *hist) + +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 */ + /* bisection search alg. find q_bin index of order Log(n) time */ int mid; int min = 0; int max = (*info).histsize-1; @@ -94,211 +96,263 @@ static int find_q_bin_index(double q, struct histogram_info *info, return mid; } -static int add_peak_to_histogram(double fs, double ss, double intensity, - struct image *image, struct histogram_info *info, - struct bin_stats *hist) + +/* 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) { - /* used for h5 files, peak/feature list and stream positions*/ - /* See Knuth TAOCP vol 2, 3rd ed, pg 232 for running variance*/ 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)) { + + /* 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].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; + 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) + struct histogram_info *info, + struct bin_stats *hist) { - /* used for d and hkl of stream files where redundancy = 1 */ double delta; int i; - /* ignore q value if outside of range */ - if ((q<(*info).q_min) || (q>(*info).q_max)) { + /* 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].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; + 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) + +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; + int i = 0; - /* ignore q value if outside of range */ - if ((q<(*info).q_min) || (q>(*info).q_max)) { + /* 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 - The number of occurance of that reflection smeered out over the - surface area which is 4 Pi q^2 the 4 pi is left out as it is a - common consant and the total is in arbiturary units. - */ - for (i=0; i<redundancy; i++) { - if (q_scaling) { + + /* 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; + + return 0; } -static int histogram_setup(struct histogram_info *info, +static int histogram_setup(struct histogram_info *info, struct bin_stats *histdata) { int i; double x; - /* setup histogram*/ - 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++) { - histdata[i].N = 0; - histdata[i].total = 0.0; - histdata[i].mean = 0.0; - histdata[i].std_dev = 0.0; + + 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; - histdata[i].q_min = pow(( i *(*info).q_delta) + - pow((*info).q_min, x), 1.0/x); - histdata[i].q_max = pow(((i+1.0)*(*info).q_delta) + - pow((*info).q_min, x), 1.0/x); - histdata[i].q_value = pow(((i+0.5)*(*info).q_delta) + - pow((*info).q_min, x), 1.0/x); + + 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) + struct bin_stats *hist, struct image *image) { - struct panel *p; - struct rvec r; - int i,fs,ss; - double q, q_fs, q_ss; - double pi = 3.14159265; + int fs,ss; + int bin; - /* check if detector geometry is loaded */ - if ( ((*image).det == NULL) || ((*image).lambda < 0.0)) { - return 1; - } + /* Check that detector geometry is present and wavelength is valid */ + if ( (image->det == NULL) || (image->lambda < 0.0) ) return 1; - /* loop over all pixels */ + /* Loop over all pixels */ for ( ss=0; ss<(*image).height; ss++ ) { for ( fs=0; fs<(*image).width; fs++ ) { - r = get_q(image, fs, ss, NULL, 1.0/ (*image).lambda); + + 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) ) { - - /* using the panel to determine the pixel size - takes care of edge pixels correctly */ - p = find_panel((*image).det, fs, ss); + + /* 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); + 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); + (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)); //approx pixel size - } + + hist[i].fract = hist[i].fract + fabs((q_fs-q)*(q_ss-q)); + + } + } } - /* divide measured area by ring area */ - for ( i=0; i<(*info).histsize; i++) { - hist[i].fract = hist[i].fract / - (pi * (pow(hist[i].q_max,2.0) - pow(hist[i].q_min,2.0))); + + /* 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, unsigned int processing_total) + 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]; - progress_bar(fs + (*image).width*ss, processing_total, - "Processing"); - if (in_bad_region((*image).det,fs,ss) == 0) { - add_peak_to_histogram(fs, ss, 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); - } } - } - return 0; -} + } + progress_bar(ss, image->height, "Processing"); + } + return 0; +} + static unsigned int process_hkl(struct image *image, char *sym, UnitCell *cell, - struct histogram_info *info, struct bin_stats *histdata, - int q_scaling,unsigned int processing_total, int use_redundency) + 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; - for ( refl = first_refl((*image).reflections, &iter); - refl != NULL; - refl = next_refl(refl, iter) ) { - i++; - progress_bar(i, processing_total, "Processing"); - get_indices(refl, &h, &k, &l); - intensity = get_intensity(refl); - if (use_redundency == 1) { - redundancy = get_redundancy(refl); - } else { - redundancy = num_equivs(h, k, l, sym); - } - /*note: resolution returns 1/d/2 in m */ - q = 2.0 * resolution(cell, h, k, l); - //sigma = get_esd_intensity(refl); - add_hkl_to_histogram(q, intensity, redundancy, - q_scaling, info, histdata); - n_peaks+=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; + + 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 { + redundancy = num_equivs(h, k, l, sym); } + + /* 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"); + + } + return n_peaks; } -static unsigned int process_stream_reflection(FILE *fh, struct image *image, - struct histogram_info *info, struct bin_stats *histdata, - unsigned int processing_total, unsigned int *n_patterns) + +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; @@ -306,175 +360,263 @@ static unsigned int process_stream_reflection(FILE *fh, struct image *image, 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); + /* Get data from next chunk */ + rval = read_chunk(fh, image); if ( rval ) continue; - i++; - progress_bar(i, processing_total, "Processing"); - // check if the pattern indexed, if so use those peaks - if ((*image).reflections != NULL) { - (*n_patterns)++; //inc number of patterns used - for ( refl = first_refl((*image).reflections, &iter); - refl != NULL; + + /* 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 */ + /* 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) == 0) { + + 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); - } while (rval == 0); + + 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_d(FILE *fh, struct image *image, - struct histogram_info *info, struct bin_stats *histdata, - unsigned int processing_total, unsigned int *n_patterns) + +static unsigned int process_stream_d(FILE *fh, struct image *image, + struct histogram_info *info, + struct bin_stats *histdata, + unsigned int *n_patterns) { - int h, k, l, rval; + int h, k, l, rval; unsigned int i = 0; - unsigned int n_peaks = 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); + + /* Get data from next chunk */ + rval = read_chunk(fh, image); if ( rval ) continue; - i++; - progress_bar(i, processing_total, "Processing"); - // check if the pattern indexed, if so use those peaks - if ((*image).reflections != NULL) { + + if ( image->reflections != NULL ) { + (*n_patterns)++; //inc number of patterns used - for ( refl = first_refl((*image).reflections, &iter); - refl != NULL; - refl = next_refl(refl, iter) ) { + + for ( refl = first_refl(image->reflections, &iter); + refl != NULL; + refl = next_refl(refl, iter) ) + { get_indices(refl, &h, &k, &l); intensity = get_intensity(refl); - /*resolution returns 1/d/2 in nm */ - q = 2.0 * resolution((*image).indexed_cell, h, k, l); - if (add_d_to_histogram(q, intensity, info, - histdata) == 0) { - n_peaks++; - } + q = 2.0 * resolution(image->indexed_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); - } while (rval == 0); + + 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 processing_total, unsigned int *n_patterns) +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 h, k, l, rval; + int rval; unsigned int i = 0; - unsigned int n_peaks = 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); + + /* Get data from next chunk */ + rval = read_chunk(fh, image); if ( rval ) continue; - i++; - progress_bar(i, processing_total, "Processing"); - if ((*image).reflections != NULL) { - n_patterns++; //inc number of patterns used - for ( refl = first_refl((*image).reflections, &iter); + if ( image->reflections != NULL ) { + + (*n_patterns)++; + + for ( refl = first_refl(image->reflections, &iter); refl != NULL; - refl = next_refl(refl, iter) ) { + 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) == 0) { - n_peaks++; - } + + 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); - } while (rval == 0); + } + + 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_peaks(FILE *fh, struct image *image, - struct histogram_info *info, struct bin_stats *histdata, - unsigned int processing_total, unsigned int *n_patterns, int only_indexed) + +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) { - struct imagefeature *f; - int rval; + int rval; unsigned int i = 0; unsigned int n_peaks = 0; - unsigned int j; + unsigned int processing_total; + + processing_total = count_patterns(fh); + rewind(fh); + do { - /* Get data from next chunk */ - rval = read_chunk(fh, image); + + /* Get data from next chunk */ + rval = read_chunk(fh, image); if ( rval ) continue; - i++; - progress_bar(i, processing_total, "Processing"); - if ((*image).features != NULL) { - if ((only_indexed == 0) || - ((only_indexed == 1) && - ((*image).reflections != NULL))) { - (*n_patterns)++; - for (j = 0; - j<image_feature_count((*image).features); - j++) { - f = image_get_feature((*image).features, j); - if (in_bad_region((*image).det, - (*f).fs,(*f).ss) == 0) { - if (add_peak_to_histogram((*f).fs, - (*f).ss, (*f).intensity, image, - info, histdata) == 0) { - n_peaks++; - } - } - } + + 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); - } while (rval == 0); + + 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 processing_total, unsigned int *n_patterns) + +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) { int fs, ss, rval; double intensity; unsigned int i = 0; unsigned int n_peaks = 0; - struct hdfile *hdfile = NULL; + 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); + + /* Get data from next chunk */ + rval = read_chunk(fh, image); if ( rval ) continue; - i++; - progress_bar(i, processing_total, "Processing"); - if ((only_indexed == 0) || - ((only_indexed == 1) && - ((*image).reflections != NULL))) { + if ((only_indexed == 0) || + ((only_indexed == 1) && + ((*image).reflections != NULL))) { hdfile = hdfile_open((*image).filename); hdfile_set_image(hdfile, "/data/data"); hdf5_read(hdfile, image, config_satcorr); @@ -482,23 +624,35 @@ static unsigned int process_stream_h5(FILE *fh, struct image *image, 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) == 0) { - add_peak_to_histogram(fs, ss, intensity, - image, info, histdata); + + 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); - } while (rval == 0); + + 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); @@ -507,19 +661,19 @@ static void show_help(const char *s) "\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" +" -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 wiht n bins (default is 100).\n" +" -s, --bins=n Makes histogram with n bins (default is 100).\n" " --spacing=<type> Use 'type' to select the q spacing.\n" " Choose from:\n" " linear : linear (default)\n" " q2 : even spacing in Wilson plots\n" -" volume : constant volume spacing\n" -" --q-max=n The max q to be considered in plot.\n" -" --q-min=n The min q to be considered in plot.\n" +" volume : constant volume\n" +" --q-max=n The maximum q to be considered in plot.\n" +" --q-min=n The minimum q 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" @@ -534,51 +688,53 @@ static void show_help(const char *s) " regions\n" " --no-sat-corr Don't correct values of saturated peaks using a\n" " table included in the HDF5 file.\n" -" --only-indexed Use wiht -data=peaks or h5 if you want to use the\n" +" --only-indexed Use with -data=peaks or h5 if you want to use the\n" " peak list of only indexed patterns\n" " --no-q-scaling Use with .hkl files if you want to not scale the\n" " powder by 1/q^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-redundency Use with .hkl files if you want to use the number\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" "\n"); } + int main(int argc, char *argv[]) { FILE *fh = NULL; - UnitCell *cell = NULL; - struct image image; + UnitCell *cell = NULL; + struct image image; struct hdfile *hdfile = NULL; struct bin_stats *histdata = NULL; struct histogram_info hist_info; - hist_info.histsize = 100; //default settings + /* 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; - + unsigned int n_patterns = 0; - unsigned int n_peaks = 0; - + unsigned int n_peaks = 0; + int c, rval, file_type, data_type; - int config_satcorr = 1; //true by default - int need_geometry = 0; //false - int need_beam = 0; //false - int need_pdb = 0; //false - int only_indexed = 0; //false - int q_scaling = 1; //true by default - int ring_corr = 0; //false - int use_redundency = 0; //false - unsigned int i, processing_total; - + 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; @@ -605,16 +761,18 @@ int main(int argc, char *argv[]) {"only-indexed", 0, &only_indexed, 1 }, {"no-q-scaling", 0, &q_scaling, 0 }, {"ring-corr", 0, &ring_corr, 1 }, - {"use-redundency", 0, &use_redundency, 1 }, + {"use-redundancy", 0, &use_redundancy, 1 }, {"data", 1, NULL, 'd'}, {0, 0, NULL, 0} }; /* Short options */ - while ((c = getopt_long(argc, argv, "hi:o:g:b:p:s:d:y:", - longopts, NULL)) != -1) { + while ((c = getopt_long(argc, argv, "hi:o:g:b:p:s:d:y:", + longopts, NULL)) != -1) + { switch (c) { + case 'h' : show_help(argv[0]); return 0; @@ -647,15 +805,15 @@ int main(int argc, char *argv[]) hist_info.histsize = atoi(optarg); break; - case 1 : + case 1 : hist_info.q_max = atof(optarg); break; - case 2 : + case 2 : hist_info.q_min = atof(optarg); break; - case 3 : + case 3 : if (strcmp(optarg, "linear") == 0 ) { hist_info.spacing = LINEAR; } else if (strcmp(optarg, "q2") == 0 ) { @@ -663,8 +821,7 @@ int main(int argc, char *argv[]) } else if (strcmp(optarg, "volume") == 0) { hist_info.spacing = VOLUME; } else { - ERROR("Failed to read spacing plot type: '%s'\n", - optarg); + ERROR("Invalid spacing type: '%s'\n", optarg); return 1; } break; @@ -679,79 +836,93 @@ int main(int argc, char *argv[]) default : return 1; } + } /* Process input file type */ if ( filename == NULL ) { + ERROR("You must specify the input filename with -i\n"); return 1; - } else if ( is_stream(filename) == 1 ) { //open .stream file + + } + + if ( is_stream(filename) ) { + file_type = FILE_STREAM; fh = fopen(filename, "r"); if ( fh == NULL ) { ERROR("Failed to open input file\n"); return 1; } - rval = read_chunk(fh, &image); //read chunk to get wavelength + /* Use wavelength from first chunk */ + rval = read_chunk(fh, &image); rewind(fh); - processing_total = count_patterns(fh); - rewind(fh); - - } else if ( H5Fis_hdf5(filename) > 0) { //open .h5 file + + } else if ( H5Fis_hdf5(filename) > 0 ) { + file_type = FILE_H5; - need_geometry = 1; //true + need_geometry = 1; hdfile = hdfile_open(filename); hdfile_set_image(hdfile, "/data/data"); hdf5_read(hdfile, &image, config_satcorr); hdfile_close(hdfile); - processing_total = image.width * image.height; + } else { + image.reflections = read_reflections(filename); - if (image.reflections != NULL) { //open .hkl file + + if ( image.reflections != NULL ) { file_type = FILE_HKL; need_pdb = 1; - processing_total = num_reflections(image.reflections); } else { - ERROR("Input file: %s is an invalid type.\n",filename); + ERROR("Couldn't recognise %s as reflection list," + " stream or image.\n", filename); return 1; } + } free(filename); - /*Process data read type */ if ( datatype == NULL ) { data_type = PLOT_D; if ((hist_info.q_min < 0.0) || (hist_info.q_max < 0.0)) { need_geometry = 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; if ((hist_info.q_min <= 0.0) || (hist_info.q_max <= 0.0)) { need_geometry = 1; } + } else if ( strcmp(datatype, "d") == 0 ) { data_type = PLOT_D; if ((hist_info.q_min <= 0.0) || (hist_info.q_max <= 0.0)) { need_geometry = 1; } + } else if ( strcmp(datatype, "peaks") == 0 ) { data_type = PLOT_PEAKS; need_geometry = 1; - } - else if ( strcmp(datatype, "h5") == 0 ) { + + } 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; } - /* logic checks */ - if ((need_geometry) && (image.lambda < 0.0)) { + /* Logic checks */ + if ( need_geometry && (image.lambda < 0.0) ) { need_beam = 1; } if ( hist_info.histsize <= 0 ) { @@ -759,16 +930,16 @@ int main(int argc, char *argv[]) "bins\n"); return 1; } - if (hist_info.q_min > hist_info.q_max) { + if ( hist_info.q_min > hist_info.q_max ) { ERROR("the minimum q value of: %e " - "is greator then your max q value of: %e\n", + "is greator then your max q value of: %e\n", hist_info.q_min, hist_info.q_max); return 1; } - /*get geometry, beam and pdb files and parameters as needed */ - if (need_geometry) { - if (geometry == NULL) { + /* 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; @@ -782,10 +953,11 @@ int main(int argc, char *argv[]) image.width = image.det->max_fs; image.height = image.det->max_ss; } - } + } free(geometry); - - if (need_beam) { + + if ( need_beam ) { + if ( beamf == NULL ) { ERROR("No wavelength in file, so you need to specify " "a beam parameters file with --beam\n"); @@ -800,10 +972,12 @@ int main(int argc, char *argv[]) image.lambda = ph_en_to_lambda(eV_to_J( image.beam->photon_energy)); } + } free(beamf); - if (need_pdb) { + if ( need_pdb ) { + if (pdb == NULL) { ERROR("You need to specify a pdb file with --pdb.\n"); return 1; @@ -815,6 +989,7 @@ int main(int argc, char *argv[]) return 1; } } + } free(pdb); @@ -822,82 +997,80 @@ int main(int argc, char *argv[]) sym = strdup("1"); } - /* setup histogram info*/ - if (hist_info.q_min <= 0 ) { + /* Set up histogram info*/ + if (hist_info.q_min <= 0.0 ) { hist_info.q_min = smallest_q(&image); } - if (hist_info.q_max <= 0 ) { + if (hist_info.q_max <= 0.0 ) { hist_info.q_max = largest_q(&image); } 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.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)) / + pow(hist_info.q_min, 3.0)) / hist_info.histsize; } - /* setup histogram data */ + /* Set up histogram data */ histdata = malloc((hist_info.histsize) * sizeof(struct bin_stats)); histogram_setup(&hist_info, histdata); - /* setup ring scaling */ - if (ring_corr == 1) { - if (ring_fraction_calc(&hist_info, histdata, &image) == 1) { - ERROR("Detector is not loaded could not correct for" - " finiate ring measurement. Do not use --ring-corr\n"); + /* 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*/ + /* 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, processing_total); + 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, processing_total, use_redundency); + q_scaling, use_redundancy); break; case FILE_STREAM : - switch (data_type) { - case PLOT_REFL : - n_peaks = process_stream_reflection(fh, &image, - &hist_info, histdata, processing_total, - &n_patterns); - break; - case PLOT_D : - n_peaks = process_stream_d(fh, &image, &hist_info, - histdata, processing_total, - &n_patterns); - break; - case PLOT_HKL : - n_peaks = process_stream_hkl(fh, &image, &hist_info, - histdata, cell, processing_total, - &n_patterns); - break; - case PLOT_PEAKS : - n_peaks = process_stream_peaks(fh, &image, &hist_info, - histdata, processing_total, - &n_patterns, only_indexed); - break; - case PLOT_H5 : - n_peaks = process_stream_h5(fh, &image, &hist_info, - histdata, config_satcorr, only_indexed, - processing_total, &n_patterns); - break; - default : - break; - } + switch (data_type) { + case PLOT_REFL : + n_peaks = process_stream_reflection(fh, &image, + &hist_info, histdata, &n_patterns); + break; + case PLOT_D : + n_peaks = process_stream_d(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); + break; + default : + break; + } fclose(fh); break; @@ -906,65 +1079,52 @@ int main(int argc, char *argv[]) } - /* sqrt the variance to get the std_dev */ - - for(i=0; i<hist_info.histsize; i++) { + /* 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 == 1) { - for(i=0; i<hist_info.histsize; i++) { - histdata[i].N = histdata[i].N / + 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].total = histdata[i].total / histdata[i].fract; - histdata[i].mean = histdata[i].mean / + histdata[i].mean = histdata[i].mean / histdata[i].fract; - histdata[i].std_dev = histdata[i].total / + histdata[i].std_dev = histdata[i].total / histdata[i].fract; } } - /* print out the results */ - if (output != NULL) { + /* Print out the results */ + if ( output != NULL ) { fh = fopen(output, "w"); if ( fh == NULL ) { ERROR("Failed to open output file\n"); return 1; } - fprintf(fh,"I read %i patterns with %i peaks\n", - n_patterns, n_peaks); - fprintf(fh,"q\tN\ttotal\tmean\tstd dev\t std dev of mean\n"); - } - else { - printf("I read %i patterns with %i peaks\n", - n_patterns, n_peaks); - printf("q\tN\ttotal\tmean\tstd dev\t std dev of mean\n"); - } - for(i=0; i<hist_info.histsize; i++) { - if (output != NULL) { - 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)); - } - else { - printf("%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)); - } + } else { + fh = stdout; + } + + fprintf(fh, "I read %i patterns with %i peaks\n", n_patterns, n_peaks); + fprintf(fh, "q\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); - if (output != NULL) { - fclose(fh); - free(output); - } + fclose(fh); free(histdata); free(sym); |