From 380ec553c04576fc3dc4c816127078d3c2cf9e32 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 25 Oct 2010 16:00:23 +0200 Subject: Read beam parameters from file (where appropriate) --- doc/examples/lcls-dec.beam | 19 +++++++ doc/examples/lcls-june.beam | 19 +++++++ src/Makefile.am | 12 +++-- src/Makefile.in | 24 +++++---- src/beam-parameters.c | 126 ++++++++++++++++++++++++++++++++++++++++++++ src/beam-parameters.h | 37 +++++++++++++ src/compare_hkl.c | 5 +- src/detector.c | 18 ++++--- src/diffraction.c | 5 +- src/facetron.c | 3 +- src/get_hkl.c | 25 ++++++++- src/image.h | 1 + src/parameters-lcls.tmp | 20 ------- src/parameters-xpp.tmp | 17 ------ src/parameters.tmp | 1 - src/pattern_sim.c | 29 ++++++++-- src/process_hkl.c | 27 +++++++++- src/reflections.c | 10 ++-- src/reflections.h | 3 +- 19 files changed, 326 insertions(+), 75 deletions(-) create mode 100644 doc/examples/lcls-dec.beam create mode 100644 doc/examples/lcls-june.beam create mode 100644 src/beam-parameters.c create mode 100644 src/beam-parameters.h delete mode 100644 src/parameters-lcls.tmp delete mode 100644 src/parameters-xpp.tmp delete mode 100644 src/parameters.tmp diff --git a/doc/examples/lcls-dec.beam b/doc/examples/lcls-dec.beam new file mode 100644 index 00000000..e19709a3 --- /dev/null +++ b/doc/examples/lcls-dec.beam @@ -0,0 +1,19 @@ +; Number of photons per pulse +beam/fluence = 2.0e11 + +; Radius of X-ray beam +beam/radius = 3.0e-6 / 2.0 + +; Photon energy in eV +beam/photon_energy = 2000.0 + + +; Detector's quantum efficiency +detector/dqe = 0.9 + +; Number of detector ADU per photon +detector/adu_per_photon = 167.0 + + +; Radius in metres +jet/radius = 3.0e-6 / 2.0 diff --git a/doc/examples/lcls-june.beam b/doc/examples/lcls-june.beam new file mode 100644 index 00000000..cc6cf28f --- /dev/null +++ b/doc/examples/lcls-june.beam @@ -0,0 +1,19 @@ +; Number of photons per pulse +beam/fluence = 2.0e11 + +; Radius of X-ray beam +beam/radius = 3.0e-6 / 2.0 + +; Photon energy in eV +beam/photon_energy = 2000.0 + + +; Detector's quantum efficiency +detector/dqe = 0.9 + +; Number of detector ADU per photon +detector/adu_per_photon = 20.9 + + +; Radius in metres +jet/radius = 3.0e-6 / 2.0 diff --git a/src/Makefile.am b/src/Makefile.am index 9f9e2517..3dd67727 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -10,7 +10,8 @@ AM_CFLAGS = -Wall AM_CPPFLAGS = -DDATADIR=\""$(datadir)"\" pattern_sim_SOURCES = pattern_sim.c diffraction.c utils.c image.c cell.c \ - hdf5-file.c detector.c sfac.c peaks.c reflections.c + hdf5-file.c detector.c sfac.c peaks.c reflections.c \ + beam-parameters.c if HAVE_OPENCL pattern_sim_SOURCES += diffraction-gpu.c cl-utils.c endif @@ -18,7 +19,7 @@ endif pattern_sim_LDADD = @LIBS@ process_hkl_SOURCES = process_hkl.c sfac.c statistics.c cell.c utils.c \ - reflections.c symmetry.c stream.c + reflections.c symmetry.c stream.c beam-parameters.c process_hkl_LDADD = @LIBS@ indexamajig_SOURCES = indexamajig.c hdf5-file.c utils.c cell.c image.c \ @@ -37,7 +38,8 @@ hdfsee_SOURCES = hdfsee.c displaywindow.c render.c hdf5-file.c utils.c image.c \ hdfsee_LDADD = @LIBS@ endif -get_hkl_SOURCES = get_hkl.c sfac.c cell.c utils.c reflections.c symmetry.c +get_hkl_SOURCES = get_hkl.c sfac.c cell.c utils.c reflections.c symmetry.c \ + beam-parameters.c get_hkl_LDADD = @LIBS@ compare_hkl_SOURCES = compare_hkl.c sfac.c cell.c utils.c reflections.c \ @@ -78,6 +80,6 @@ INCLUDES = "-I$(top_srcdir)/data" EXTRA_DIST = cell.h hdf5-file.h image.h utils.h diffraction.h detector.h \ sfac.h reflections.h list_tmp.h statistics.h displaywindow.h \ render.h hdfsee.h dirax.h peaks.h index.h filters.h \ - diffraction-gpu.h cl-utils.h parameters-lcls.tmp symmetry.h \ + diffraction-gpu.h cl-utils.h symmetry.h \ povray.h index-priv.h geometry.h templates.h render_hkl.h \ - stream.h thread-pool.h parameters.tmp + stream.h thread-pool.h beam-parameters.h diff --git a/src/Makefile.in b/src/Makefile.in index 7e34ee08..c5a6c346 100644 --- a/src/Makefile.in +++ b/src/Makefile.in @@ -82,7 +82,8 @@ am_facetron_OBJECTS = facetron.$(OBJEXT) cell.$(OBJEXT) \ facetron_OBJECTS = $(am_facetron_OBJECTS) facetron_DEPENDENCIES = am_get_hkl_OBJECTS = get_hkl.$(OBJEXT) sfac.$(OBJEXT) cell.$(OBJEXT) \ - utils.$(OBJEXT) reflections.$(OBJEXT) symmetry.$(OBJEXT) + utils.$(OBJEXT) reflections.$(OBJEXT) symmetry.$(OBJEXT) \ + beam-parameters.$(OBJEXT) get_hkl_OBJECTS = $(am_get_hkl_OBJECTS) get_hkl_DEPENDENCIES = am__hdfsee_SOURCES_DIST = hdfsee.c displaywindow.c render.c \ @@ -109,11 +110,12 @@ indexamajig_OBJECTS = $(am_indexamajig_OBJECTS) indexamajig_DEPENDENCIES = am__pattern_sim_SOURCES_DIST = pattern_sim.c diffraction.c utils.c \ image.c cell.c hdf5-file.c detector.c sfac.c peaks.c \ - reflections.c diffraction-gpu.c cl-utils.c + reflections.c beam-parameters.c diffraction-gpu.c cl-utils.c am_pattern_sim_OBJECTS = pattern_sim.$(OBJEXT) diffraction.$(OBJEXT) \ utils.$(OBJEXT) image.$(OBJEXT) cell.$(OBJEXT) \ hdf5-file.$(OBJEXT) detector.$(OBJEXT) sfac.$(OBJEXT) \ - peaks.$(OBJEXT) reflections.$(OBJEXT) $(am__objects_1) + peaks.$(OBJEXT) reflections.$(OBJEXT) \ + beam-parameters.$(OBJEXT) $(am__objects_1) pattern_sim_OBJECTS = $(am_pattern_sim_OBJECTS) pattern_sim_DEPENDENCIES = am_powder_plot_OBJECTS = powder_plot.$(OBJEXT) cell.$(OBJEXT) \ @@ -123,7 +125,8 @@ powder_plot_OBJECTS = $(am_powder_plot_OBJECTS) powder_plot_DEPENDENCIES = am_process_hkl_OBJECTS = process_hkl.$(OBJEXT) sfac.$(OBJEXT) \ statistics.$(OBJEXT) cell.$(OBJEXT) utils.$(OBJEXT) \ - reflections.$(OBJEXT) symmetry.$(OBJEXT) stream.$(OBJEXT) + reflections.$(OBJEXT) symmetry.$(OBJEXT) stream.$(OBJEXT) \ + beam-parameters.$(OBJEXT) process_hkl_OBJECTS = $(am_process_hkl_OBJECTS) process_hkl_DEPENDENCIES = am_reintegrate_OBJECTS = reintegrate.$(OBJEXT) cell.$(OBJEXT) \ @@ -270,10 +273,10 @@ AM_CFLAGS = -Wall AM_CPPFLAGS = -DDATADIR=\""$(datadir)"\" pattern_sim_SOURCES = pattern_sim.c diffraction.c utils.c image.c \ cell.c hdf5-file.c detector.c sfac.c peaks.c reflections.c \ - $(am__append_2) + beam-parameters.c $(am__append_2) pattern_sim_LDADD = @LIBS@ process_hkl_SOURCES = process_hkl.c sfac.c statistics.c cell.c utils.c \ - reflections.c symmetry.c stream.c + reflections.c symmetry.c stream.c beam-parameters.c process_hkl_LDADD = @LIBS@ indexamajig_SOURCES = indexamajig.c hdf5-file.c utils.c cell.c image.c \ @@ -285,7 +288,9 @@ indexamajig_LDADD = @LIBS@ @HAVE_GTK_TRUE@ filters.c @HAVE_GTK_TRUE@hdfsee_LDADD = @LIBS@ -get_hkl_SOURCES = get_hkl.c sfac.c cell.c utils.c reflections.c symmetry.c +get_hkl_SOURCES = get_hkl.c sfac.c cell.c utils.c reflections.c symmetry.c \ + beam-parameters.c + get_hkl_LDADD = @LIBS@ compare_hkl_SOURCES = compare_hkl.c sfac.c cell.c utils.c reflections.c \ statistics.c symmetry.c @@ -323,9 +328,9 @@ INCLUDES = "-I$(top_srcdir)/data" EXTRA_DIST = cell.h hdf5-file.h image.h utils.h diffraction.h detector.h \ sfac.h reflections.h list_tmp.h statistics.h displaywindow.h \ render.h hdfsee.h dirax.h peaks.h index.h filters.h \ - diffraction-gpu.h cl-utils.h parameters-lcls.tmp symmetry.h \ + diffraction-gpu.h cl-utils.h symmetry.h \ povray.h index-priv.h geometry.h templates.h render_hkl.h \ - stream.h thread-pool.h parameters.tmp + stream.h thread-pool.h beam-parameters.h all: all-am @@ -444,6 +449,7 @@ mostlyclean-compile: distclean-compile: -rm -f *.tab.c +@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/beam-parameters.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/calibrate_detector.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/cell.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/cl-utils.Po@am__quote@ diff --git a/src/beam-parameters.c b/src/beam-parameters.c new file mode 100644 index 00000000..73dc2897 --- /dev/null +++ b/src/beam-parameters.c @@ -0,0 +1,126 @@ +/* + * beam-parameters.c + * + * Beam parameters + * + * (c) 2006-2010 Thomas White + * + * Part of CrystFEL - crystallography with a FEL + * + */ + + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include + +#include "beam-parameters.h" +#include "utils.h" + + +struct beam_params *get_beam_parameters(const char *filename) +{ + FILE *fh; + struct beam_params *b; + char *rval; + int reject; + + fh = fopen(filename, "r"); + if ( fh == NULL ) return NULL; + + b = calloc(1, sizeof(struct beam_params)); + if ( b == NULL ) return NULL; + + b->fluence = -1.0; + b->beam_radius = -1.0; + b->photon_energy = -1.0; + b->dqe = -1.0; + b->adu_per_photon = -1.0; + b->water_radius = -1.0; + + do { + + int n1; + char line[1024]; + int i; + char **bits; + + rval = fgets(line, 1023, fh); + if ( rval == NULL ) break; + chomp(line); + + n1 = assplode(line, " \t", &bits, ASSPLODE_NONE); + if ( n1 < 3 ) { + for ( i=0; ifluence = atof(bits[0]); + } else if ( strcmp(bits[0], "beam/radius") == 0 ) { + b->beam_radius = atof(bits[0]); + } else if ( strcmp(bits[0], "beam/photon_energy") == 0 ) { + b->photon_energy = atof(bits[0]); + } else if ( strcmp(bits[0], "detector/dqe") == 0 ) { + b->dqe = atof(bits[0]); + } else if ( strcmp(bits[0], "detector/adu_per_photon") == 0 ) { + b->adu_per_photon = atof(bits[0]); + } else if ( strcmp(bits[0], "jet/radius") == 0 ) { + b->water_radius = atof(bits[0]); + } else { + ERROR("Unrecognised field '%s'\n", bits[0]); + } + + for ( i=0; ifluence < 0.0 ) { + ERROR("Invalid or unspecified value for 'beam/fluence'.\n"); + reject = 1; + } + if ( b->beam_radius < 0.0 ) { + ERROR("Invalid or unspecified value for 'beam/radius'.\n"); + reject = 1; + } + if ( b->photon_energy < 0.0 ) { + ERROR("Invalid or unspecified value for" + " 'beam/photon_energy'.\n"); + reject = 1; + } + if ( b->dqe < 0.0 ) { + ERROR("Invalid or unspecified value for 'detector/dqe'.\n"); + reject = 1; + } + if ( b->adu_per_photon < 0.0 ) { + ERROR("Invalid or unspecified value for" + " 'detector/adu_per_photon'.\n"); + reject = 1; + } + if ( b->water_radius < 0.0 ) { + ERROR("Invalid or unspecified value for 'jet/radius'.\n"); + reject = 1; + } + + if ( reject ) { + ERROR("Please fix the above problems with the beam" + " parameters file and try again.\n"); + return NULL; + } + + return b; +} diff --git a/src/beam-parameters.h b/src/beam-parameters.h new file mode 100644 index 00000000..84d95b49 --- /dev/null +++ b/src/beam-parameters.h @@ -0,0 +1,37 @@ +/* + * beam-parameters.h + * + * Beam parameters + * + * (c) 2006-2010 Thomas White + * + * Part of CrystFEL - crystallography with a FEL + * + */ + + +#ifndef BEAM_PARAMETERS_H +#define BEAM_PARAMETERS_H + +#ifdef HAVE_CONFIG_H +#include +#endif + + +struct beam_params +{ + double fluence; /* photons per pulse */ + double beam_radius; /* metres */ + double photon_energy; /* eV per photon */ + + double dqe; /* Detector DQE (fraction) */ + double adu_per_photon; /* Detector "gain" */ + + double water_radius; /* metres */ +}; + + +extern struct beam_params *get_beam_parameters(const char *filename); + + +#endif /* BEAM_PARAMETERS_H */ diff --git a/src/compare_hkl.c b/src/compare_hkl.c index 1100a05e..2bdbd5a4 100644 --- a/src/compare_hkl.c +++ b/src/compare_hkl.c @@ -505,7 +505,10 @@ int main(int argc, char *argv[]) } if ( outfile != NULL ) { - write_reflections(outfile, icommon, out, NULL, NULL, cell); + + write_reflections(outfile, icommon, out, NULL, NULL, cell, 1.0); + STATUS("Sigma(I) values in output file are not meaningful.\n"); + } return 0; diff --git a/src/detector.c b/src/detector.c index cc5663af..2a4e3579 100644 --- a/src/detector.c +++ b/src/detector.c @@ -20,7 +20,7 @@ #include "utils.h" #include "diffraction.h" #include "detector.h" -#include "parameters.tmp" +#include "beam-parameters.h" int atob(const char *a) @@ -88,14 +88,14 @@ void record_image(struct image *image, int do_poisson) double max_tt = 0.0; /* How many photons are scattered per electron? */ - area = M_PI*pow(BEAM_RADIUS, 2.0); - total_energy = FLUENCE * ph_lambda_to_en(image->lambda); + area = M_PI*pow(image->beam->beam_radius, 2.0); + total_energy = image->beam->fluence * ph_lambda_to_en(image->lambda); energy_density = total_energy / area; - ph_per_e = (FLUENCE/area) * pow(THOMSON_LENGTH, 2.0); + ph_per_e = (image->beam->fluence /area) * pow(THOMSON_LENGTH, 2.0); STATUS("Fluence = %8.2e photons, " "Energy density = %5.3f kJ/cm^2, " "Total energy = %5.3f microJ\n", - FLUENCE, energy_density/1e7, total_energy*1e6); + image->beam->fluence, energy_density/1e7, total_energy*1e6); for ( x=0; xwidth; x++ ) { for ( y=0; yheight; y++ ) { @@ -135,13 +135,15 @@ void record_image(struct image *image, int do_poisson) sa = proj_area / (dsq + Lsq); if ( do_poisson ) { - counts = poisson_noise(intensity * ph_per_e * sa * DQE); + counts = poisson_noise(intensity * ph_per_e + * sa * image->beam->dqe ); } else { - cf = intensity * ph_per_e * sa * DQE; + cf = intensity * ph_per_e * sa * image->beam->dqe; counts = cf; } - image->data[x + image->width*y] = counts * DETECTOR_GAIN; + image->data[x + image->width*y] = counts + * image->beam->adu_per_photon; if ( isinf(image->data[x+image->width*y]) ) { ERROR("Processed infinity at %i,%i\n", x, y); } diff --git a/src/diffraction.c b/src/diffraction.c index aea47035..00a73b94 100644 --- a/src/diffraction.c +++ b/src/diffraction.c @@ -22,7 +22,7 @@ #include "cell.h" #include "diffraction.h" #include "sfac.h" -#include "parameters.tmp" +#include "beam-parameters.h" #define SAMPLING (4) @@ -374,7 +374,8 @@ void get_diffraction(struct image *image, int na, int nb, int nc, /* Add intensity contribution from water */ image->data[x + image->width*y] += water_diffraction(q, ph_lambda_to_en(image->lambda), - BEAM_RADIUS, WATER_RADIUS) * sw; + image->beam->beam_radius, + image->beam->water_radius) * sw; } diff --git a/src/facetron.c b/src/facetron.c index 9bef1682..8d66d549 100644 --- a/src/facetron.c +++ b/src/facetron.c @@ -428,7 +428,8 @@ int main(int argc, char *argv[]) } /* Output results */ - write_reflections(outfile, obs, i_full, NULL, NULL, NULL); + write_reflections(outfile, obs, i_full, NULL, NULL, NULL, 1.0); + STATUS("Sigma(I) values in output file are not (yet) meaningful.\n"); /* Clean up */ free(i_full); diff --git a/src/get_hkl.c b/src/get_hkl.c index dba341c6..8a8a71c2 100644 --- a/src/get_hkl.c +++ b/src/get_hkl.c @@ -25,6 +25,7 @@ #include "sfac.h" #include "reflections.h" #include "symmetry.h" +#include "beam-parameters.h" static void show_help(const char *s) @@ -50,6 +51,7 @@ static void show_help(const char *s) " --no-phases Do not try to use phases in the input file.\n" " --multiplicity Multiply intensities by the number of\n" " equivalent reflections.\n" +" -b, --beam= Get beam parameters from file (used for sigmas).\n" ); } @@ -262,6 +264,8 @@ int main(int argc, char *argv[]) ReflItemList *input_items; ReflItemList *write_items; UnitCell *cell = NULL; + double adu_per_photon; + struct beam_params *beam = NULL; /* Long options */ const struct option longopts[] = { @@ -277,6 +281,7 @@ int main(int argc, char *argv[]) {"pdb", 1, NULL, 'p'}, {"no-phases", 0, &config_nophase, 1}, {"multiplicity", 0, &config_multi, 1}, + {"beam", 1, NULL, 'b'}, {0, 0, NULL, 0} }; @@ -317,6 +322,15 @@ int main(int argc, char *argv[]) expand = strdup(optarg); break; + case 'b' : + beam = get_beam_parameters(optarg); + if ( beam == NULL ) { + ERROR("Failed to load beam parameters" + " from '%s'\n", optarg); + return 1; + } + break; + case 0 : break; @@ -407,7 +421,16 @@ int main(int argc, char *argv[]) union_items(write_items, input_items); } - write_reflections(output, write_items, ideal_ref, phases, NULL, cell); + if ( beam == NULL ) { + adu_per_photon = 167.0; + STATUS("No beam parameters file provided (use -b), " + "so I'm assuming 167.0 ADU per photon.\n"); + } else { + adu_per_photon = beam->adu_per_photon; + } + + write_reflections(output, write_items, ideal_ref, phases, NULL, cell, + adu_per_photon); delete_items(input_items); delete_items(write_items); diff --git a/src/image.h b/src/image.h index c39c19a0..0f6db034 100644 --- a/src/image.h +++ b/src/image.h @@ -73,6 +73,7 @@ struct image { UnitCell *candidate_cells[MAX_CELL_CANDIDATES]; int ncells; struct detector *det; + struct beam_params *beam; char *filename; struct cpeak *cpeaks; int n_cpeaks; diff --git a/src/parameters-lcls.tmp b/src/parameters-lcls.tmp deleted file mode 100644 index a35c712f..00000000 --- a/src/parameters-lcls.tmp +++ /dev/null @@ -1,20 +0,0 @@ -/* Number of photons in pulse */ -#define FLUENCE (2.0e11) - -/* Detector's quantum efficiency */ -#define DQE (0.9) - -/* ADU per photon (front detector) */ -/* Front detector, December */ -#define DETECTOR_GAIN (167.0) -/* Front detector, June */ -//#define DETECTOR_GAIN (20.9) - -/* Radius of the water column */ -#define WATER_RADIUS (3.0e-6 / 2.0) - -/* Radius of X-ray beam */ -#define BEAM_RADIUS (3.0e-6 / 2.0) - -/* Photon energy in eV */ -#define PHOTON_ENERGY (1780.0) diff --git a/src/parameters-xpp.tmp b/src/parameters-xpp.tmp deleted file mode 100644 index 29173ed9..00000000 --- a/src/parameters-xpp.tmp +++ /dev/null @@ -1,17 +0,0 @@ -/* Number of photons in pulse */ -#define FLUENCE (2.0e11) - -/* Detector's quantum efficiency */ -#define DQE (0.9) - -/* ADU per photon (front detector) */ -#define DETECTOR_GAIN (1.0) - -/* Radius of the water column */ -#define WATER_RADIUS (0.0e-6 / 2.0) - -/* Radius of X-ray beam */ -#define BEAM_RADIUS (1.0e-6 / 2.0) - -/* Photon energy in eV */ -#define PHOTON_ENERGY (8000.0) diff --git a/src/parameters.tmp b/src/parameters.tmp deleted file mode 100644 index f550cd1a..00000000 --- a/src/parameters.tmp +++ /dev/null @@ -1 +0,0 @@ -#include "parameters-lcls.tmp" diff --git a/src/pattern_sim.c b/src/pattern_sim.c index 5fd7e7ea..a1c72ef4 100644 --- a/src/pattern_sim.c +++ b/src/pattern_sim.c @@ -31,7 +31,7 @@ #include "peaks.h" #include "sfac.h" #include "reflections.h" -#include "parameters.tmp" +#include "beam-parameters.h" static void show_help(const char *s) @@ -47,7 +47,8 @@ static void show_help(const char *s) " intensities file)\n" " --simulation-details Show technical details of the simulation.\n" " --gpu Use the GPU to speed up the calculation.\n" -" -g, --geometry= Get detector geometry from file.\n" +" -g, --geometry= Get detector geometry from file.\n" +" -b, --beam= Get beam parameters from file.\n" "\n" " --near-bragg Output h,k,l,I near Bragg conditions.\n" " -n, --number= Generate N images. Default 1.\n" @@ -189,6 +190,7 @@ int main(int argc, char *argv[]) char *grad_str = NULL; char *outfile = NULL; char *geometry = NULL; + char *beamfile = NULL; GradientMethod grad; int ndone = 0; /* Number of simulations done (images or not) */ int number = 1; /* Number used for filename of image */ @@ -213,11 +215,12 @@ int main(int argc, char *argv[]) {"pdb", 1, NULL, 'p'}, {"output", 1, NULL, 'o'}, {"geometry", 1, NULL, 'g'}, + {"beam", 1, NULL, 'b'}, {0, 0, NULL, 0} }; /* Short options */ - while ((c = getopt_long(argc, argv, "hrn:i:t:p:o:g:", + while ((c = getopt_long(argc, argv, "hrn:i:t:p:o:g:b:", longopts, NULL)) != -1) { switch (c) { @@ -257,6 +260,10 @@ int main(int argc, char *argv[]) geometry = strdup(optarg); break; + case 'b' : + beamfile = strdup(optarg); + break; + case 0 : break; @@ -316,6 +323,12 @@ int main(int argc, char *argv[]) return 1; } + if ( beamfile == NULL ) { + ERROR("You need to specify a beam parameter file" + " with --beam\n"); + return 1; + } + if ( intfile == NULL ) { /* Gentle reminder */ STATUS("You didn't specify the file containing the "); @@ -345,10 +358,17 @@ int main(int argc, char *argv[]) } free(geometry); + image.beam = get_beam_parameters(beamfile); + if ( image.beam == NULL ) { + ERROR("Failed to read beam parameters from '%s'\n", beamfile); + return 1; + } + free(beamfile); + /* Define image parameters */ image.width = image.det->max_x + 1; image.height = image.det->max_y + 1; - image.lambda = ph_en_to_lambda(eV_to_J(PHOTON_ENERGY)); /* Wavelength */ + image.lambda = ph_en_to_lambda(eV_to_J(image.beam->photon_energy)); cell = load_cell_from_pdb(filename); if ( cell == NULL ) { exit(1); @@ -478,6 +498,7 @@ skip: free(image.det->panels); free(image.det); + free(image.beam); free(powder); cell_free(cell); free(intensities); diff --git a/src/process_hkl.c b/src/process_hkl.c index 14acd6b9..22cb8016 100644 --- a/src/process_hkl.c +++ b/src/process_hkl.c @@ -27,6 +27,7 @@ #include "reflections.h" #include "symmetry.h" #include "stream.h" +#include "beam-parameters.h" /* Number of divisions for intensity histograms */ @@ -44,6 +45,7 @@ static void show_help(const char *s) " -o, --output= Specify output filename for merged intensities\n" " (don't specify for no output).\n" " -p, --pdb= PDB file to use (default: molecule.pdb).\n" +" -b, --beam= Get beam parameters from file (used for sigmas).\n" "\n" " --max-only Take the integrated intensity to be equal to the\n" " maximum intensity measured for that reflection.\n" @@ -564,6 +566,7 @@ int main(int argc, char *argv[]) ReflItemList *reference_items; double *reference_i; FILE *outfh = NULL; + struct beam_params *beam = NULL; /* Long options */ const struct option longopts[] = { @@ -582,6 +585,7 @@ int main(int argc, char *argv[]) {"rmerge", 0, &config_rmerge, 1}, {"outstream", 1, NULL, 'a'}, {"reference", 1, NULL, 'r'}, + {"beam", 1, NULL, 'b'}, {0, 0, NULL, 0} }; @@ -630,6 +634,15 @@ int main(int argc, char *argv[]) outstream = strdup(optarg); break; + case 'b' : + beam = get_beam_parameters(optarg); + if ( beam == NULL ) { + ERROR("Failed to load beam parameters" + " from '%s'\n", optarg); + return 1; + } + break; + case 0 : break; @@ -758,7 +771,19 @@ int main(int argc, char *argv[]) } if ( output != NULL ) { - write_reflections(output, observed, model, NULL, counts, cell); + + double adu_per_photon; + + if ( beam == NULL ) { + adu_per_photon = 167.0; + STATUS("No beam parameters file provided (use -b), " + "so I'm assuming 167.0 ADU per photon.\n"); + } else { + adu_per_photon = beam->adu_per_photon; + } + + write_reflections(output, observed, model, NULL, counts, cell, + adu_per_photon); } if ( config_rmerge ) { diff --git a/src/reflections.c b/src/reflections.c index b26d9528..0474d517 100644 --- a/src/reflections.c +++ b/src/reflections.c @@ -19,12 +19,13 @@ #include "utils.h" #include "cell.h" #include "reflections.h" -#include "parameters.tmp" +#include "beam-parameters.h" void write_reflections(const char *filename, ReflItemList *items, double *intensities, double *phases, - unsigned int *counts, UnitCell *cell) + unsigned int *counts, UnitCell *cell, + double adu_per_photon) { FILE *fh; int i; @@ -79,7 +80,7 @@ void write_reflections(const char *filename, ReflItemList *items, } if ( intensity > 0.0 ) { - sigma = DETECTOR_GAIN * sqrt(intensity/DETECTOR_GAIN); + sigma = adu_per_photon * sqrt(intensity/adu_per_photon); } else { sigma = 0.0; } @@ -90,7 +91,8 @@ void write_reflections(const char *filename, ReflItemList *items, } STATUS("Warning: Errors have been estimated from Poisson distribution" - " assuming %5.2f ADU per photon.\n", DETECTOR_GAIN); + " assuming %5.2f ADU per photon.\n", adu_per_photon); + STATUS("This is not necessarily a useful estimate.\n"); fclose(fh); } diff --git a/src/reflections.h b/src/reflections.h index 84345260..7b5e5e42 100644 --- a/src/reflections.h +++ b/src/reflections.h @@ -23,7 +23,8 @@ extern void write_reflections(const char *filename, ReflItemList *items, double *intensities, double *phases, - unsigned int *counts, UnitCell *cell); + unsigned int *counts, UnitCell *cell, + double adu_per_photon); extern ReflItemList *read_reflections(const char *filename, double *intensities, double *phases, -- cgit v1.2.3