aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2010-10-25 16:00:23 +0200
committerThomas White <taw@physics.org>2012-02-22 15:27:03 +0100
commit380ec553c04576fc3dc4c816127078d3c2cf9e32 (patch)
tree7f2e4a384bf1ed2285c15817fbbfe46c1add1d64
parenta1408ca8460ef486730b4e5ce11cfca2504edf3b (diff)
Read beam parameters from file (where appropriate)
-rw-r--r--doc/examples/lcls-dec.beam19
-rw-r--r--doc/examples/lcls-june.beam19
-rw-r--r--src/Makefile.am12
-rw-r--r--src/Makefile.in24
-rw-r--r--src/beam-parameters.c126
-rw-r--r--src/beam-parameters.h37
-rw-r--r--src/compare_hkl.c5
-rw-r--r--src/detector.c18
-rw-r--r--src/diffraction.c5
-rw-r--r--src/facetron.c3
-rw-r--r--src/get_hkl.c25
-rw-r--r--src/image.h1
-rw-r--r--src/parameters-lcls.tmp20
-rw-r--r--src/parameters-xpp.tmp17
-rw-r--r--src/parameters.tmp1
-rw-r--r--src/pattern_sim.c29
-rw-r--r--src/process_hkl.c27
-rw-r--r--src/reflections.c10
-rw-r--r--src/reflections.h3
19 files changed, 326 insertions, 75 deletions
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 <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,