Commit 380ec553 authored by Thomas White's avatar Thomas White
Browse files

Read beam parameters from file (where appropriate)

parent a1408ca8
; 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
; 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
......@@ -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
......@@ -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@
......
/*
* 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;
}
/*
* 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 */
......@@ -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;
......
......@@ -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);
}
......
......@@ -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;
}
......
......@@ -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);
......
......@@ -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);
......
......@@ -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;
......
/* 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)
/* 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)
#include "parameters-lcls.tmp"
......@@ -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);
......
......@@ -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;