Read beam parameters from file (where appropriate)

17 files changed, 288 insertions, 75 deletions

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 <taw@physics.org>
+ *
+ * Part of CrystFEL - crystallography with a FEL
+ *
+ */
+
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdio.h>
+#include <stdlib.h>
+
+#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; i<n1; i++ ) free(bits[i]);
+			free(bits);
+			continue;
+		}
+
+		if ( bits[1][0] != '=' ) {
+			for ( i=0; i<n1; i++ ) free(bits[i]);
+			free(bits);
+			continue;
+		}
+
+		if ( strcmp(bits[0], "beam/fluence") == 0 ) {
+			b->fluence = 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; i<n1; i++ ) free(bits[i]);
+		free(bits);
+
+	} while ( rval != NULL );
+	fclose(fh);
+
+	reject = 0;
+
+	if ( b->fluence < 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 <taw@physics.org>
+ *
+ * Part of CrystFEL - crystallography with a FEL
+ *
+ */
+
+
+#ifndef BEAM_PARAMETERS_H
+#define BEAM_PARAMETERS_H
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#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; x<image->width; x++ ) {
 	for ( y=0; y<image->height; 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=<file>         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=<file>    Get detector geometry from file.\n"
+" -g, --geometry=<file>     Get detector geometry from file.\n"
+" -b, --beam=<file>         Get beam parameters from file.\n"
 "\n"
 "     --near-bragg          Output h,k,l,I near Bragg conditions.\n"
 " -n, --number=<N>          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=<filename>   Specify output filename for merged intensities\n"
 "                             (don't specify for no output).\n"
 "  -p, --pdb=<filename>      PDB file to use (default: molecule.pdb).\n"
+"  -b, --beam=<file>         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,