Commit e38375cb authored by Thomas White's avatar Thomas White
Browse files

Use ADU per eV in geometry file, rather than ADU per photon in beam file

Conflicts:

	src/get_hkl.c
	tests/morphology_check.c
parent be33bf07
adu_per_eV = 1.0
q0a0/min_fs = 0
q0a0/min_ss = 0
q0a0/max_fs = 193
......
......@@ -14,10 +14,3 @@ beam/bandwidth = 0.001
; Beam divergence in radians
beam/divergence = 0.008
; Detector's quantum efficiency
detector/dqe = 1.0
; Number of detector ADU per photon
detector/adu_per_photon = 1.0
......@@ -14,10 +14,3 @@ beam/bandwidth = 0.001
; Beam divergence in radians
beam/divergence = 0.001
; Detector's quantum efficiency
detector/dqe = 0.9
; Number of detector ADU per photon
detector/adu_per_photon = 167.0
......@@ -5,6 +5,7 @@ mask_bad = 0x00
; These default values will be used unless overridden by the per-panel values
peak_sep = 50.0
integr_radius = 10.0
adu_per_eV = 0.0835
upper/min_fs = 0
upper/max_fs = 1023
......
......@@ -5,6 +5,7 @@ mask_bad = 0x00
; These default values will be used unless overridden by the per-panel values
peak_sep = 50.0
integr_radius = 10.0
adu_per_eV = 0.01045
; Upper panel (nearest the beam)
0/min_fs = 0
......
......@@ -14,10 +14,3 @@ beam/bandwidth = 0.01
; Beam divergence in radians
beam/divergence = 0.0087
; Detector's quantum efficiency
detector/dqe = 0.9
; Number of detector ADU per photon
detector/adu_per_photon = 20.9
......@@ -14,10 +14,3 @@ beam/bandwidth = 0.001
; Beam divergence in radians
beam/divergence = 0.001
; Detector's quantum efficiency
detector/dqe = 1.0
; Number of detector ADU per photon
detector/adu_per_photon = 1.0
......@@ -10,3 +10,4 @@
0/res = 9090 ; 110 micron pixel size
0/badrow_direction = x
0/peak_sep = 4
0/adu_per_eV = 1.0
adu_per_eV = 1.0
; Upper panel
0/min_fs = 0
0/max_fs = 1023
......
......@@ -51,8 +51,6 @@ struct beam_params *get_beam_parameters(const char *filename)
b->photon_energy = -1.0;
b->bandwidth = -1.0;
b->divergence = -1.0;
b->dqe = -1.0;
b->adu_per_photon = -1.0;
do {
......@@ -88,10 +86,6 @@ struct beam_params *get_beam_parameters(const char *filename)
b->bandwidth = atof(bits[2]);
} else if ( strcmp(bits[0], "beam/divergence") == 0 ) {
b->divergence = atof(bits[2]);
} else if ( strcmp(bits[0], "detector/dqe") == 0 ) {
b->dqe = atof(bits[2]);
} else if ( strcmp(bits[0], "detector/adu_per_photon") == 0 ) {
b->adu_per_photon = atof(bits[2]);
} else {
ERROR("Unrecognised field '%s'\n", bits[0]);
}
......@@ -125,15 +119,6 @@ struct beam_params *get_beam_parameters(const char *filename)
ERROR("Invalid or unspecified value for 'beam/divergence'.\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 ( reject ) {
ERROR("Please fix the above problems with the beam"
......
......@@ -40,9 +40,6 @@ struct beam_params
* a rectangular distribution with this as
* its (full) width. */
double divergence; /* divergence (radians) */
double dqe; /* Detector DQE (fraction) */
double adu_per_photon; /* Detector "gain" */
};
......
......@@ -365,15 +365,14 @@ 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 * image->beam->dqe );
counts = poisson_noise(intensity * ph_per_e * sa);
} else {
cf = intensity * ph_per_e * sa * image->beam->dqe;
cf = intensity * ph_per_e * sa;
counts = cf;
}
image->data[x + image->width*y] = counts
* image->beam->adu_per_photon;
image->data[x + image->width*y] = counts * p->adu_per_eV
* ph_lambda_to_eV(image->lambda);
/* Sanity checks */
if ( isinf(image->data[x+image->width*y]) ) n_inf2++;
......@@ -577,6 +576,8 @@ static int parse_field_for_panel(struct panel *panel, const char *key,
panel->cnx = atof(val);
} else if ( strcmp(key, "corner_y") == 0 ) {
panel->cny = atof(val);
} else if ( strcmp(key, "adu_per_eV") == 0 ) {
panel->adu_per_eV = atof(val);
} else if ( strcmp(key, "rigid_group") == 0 ) {
panel->rigid_group = find_or_add_rg(det, val);
} else if ( strcmp(key, "clen") == 0 ) {
......@@ -737,6 +738,7 @@ struct detector *get_detector_geometry(const char *filename)
det->defaults.ssx = 0.0;
det->defaults.ssy = 1.0;
det->defaults.rigid_group = NULL;
det->defaults.adu_per_eV = NAN;
strncpy(det->defaults.name, "", 1023);
do {
......@@ -871,6 +873,11 @@ struct detector *get_detector_geometry(const char *filename)
" panel %s\n", det->panels[i].name);
reject = 1;
}
if ( isnan(det->panels[i].adu_per_eV) ) {
ERROR("Please specify the number of ADU per eV for"
" panel %s\n", det->panels[i].name);
reject = 1;
}
/* It's OK if the badrow direction is '0' */
/* It's not a problem if "no_index" is still zero */
/* The default peak_sep is OK (maybe) */
......@@ -1242,6 +1249,7 @@ int write_detector_geometry(const char *filename, struct detector *det)
fprintf(fh, "%s/ss = %+fx %+fy\n", p->name, p->ssx, p->ssy);
fprintf(fh, "%s/corner_x = %g\n", p->name, p->cnx);
fprintf(fh, "%s/corner_y = %g\n", p->name, p->cny);
fprintf(fh, "%s/adu_per_eV = %g\n", p->name, p->adu_per_eV);
if ( p->no_index ) {
fprintf(fh, "%s/no_index = 1\n", p->name);
......
......@@ -59,6 +59,7 @@ struct panel
double peak_sep; /* Characteristic peak separation */
double integr_radius; /* Peak integration radius */
char *rigid_group; /* Rigid group, or -1 for none */
double adu_per_eV; /* Number of ADU per eV */
double fsx;
double fsy;
......
......@@ -186,6 +186,12 @@ static inline int within_tolerance(double a, double b, double percent)
/* Joules to eV */
#define J_to_eV(a) ((a)/ELECTRON_CHARGE)
/* Photon wavelength (m) to energy (eV) */
#define ph_lambda_to_eV(a) J_to_eV(ph_lambda_to_en(a))
/* Photon energy (eV) to wavelength (m) */
#define ph_eV_to_lambda(a) ph_en_to_lambda(eV_to_J(a))
#define UNUSED __attribute__((unused))
......
......@@ -56,7 +56,7 @@ static void show_help(const char *s)
" --noise Add 10%% random noise.\n"
"\n"
"To calculate Poisson samples accurately, you must also give:\n"
" -b, --beam=<file> Get beam parameters from file.\n"
" --adu-per-photon=<n> Number of ADU per photon.\n"
"\n"
"You can artificially 'twin' the reflections, or expand them out.\n"
" -w, --twin=<sym> Generate twinned data according to the given\n"
......@@ -362,6 +362,8 @@ int main(int argc, char *argv[])
char *beamfile = NULL;
struct beam_params *beam = NULL;
RefList *input;
double adu_per_photon = 0.0;
int have_adu_per_photon = 0;
/* Long options */
const struct option longopts[] = {
......@@ -375,8 +377,8 @@ int main(int argc, char *argv[])
{"expand", 1, NULL, 'e'},
{"intensities", 1, NULL, 'i'},
{"multiplicity", 0, &config_multi, 1},
{"beam", 1, NULL, 'b'},
{"trim-centrics", 0, &config_trimc, 1},
{"adu-per-photon", 1, NULL, 2},
{0, 0, NULL, 0}
};
......@@ -413,8 +415,9 @@ int main(int argc, char *argv[])
expand_str = strdup(optarg);
break;
case 'b' :
beamfile = strdup(optarg);
case 2 :
adu_per_photon = strtof(optarg, NULL);
have_adu_per_photon = 1;
break;
case 0 :
......@@ -476,11 +479,11 @@ int main(int argc, char *argv[])
}
if ( config_poisson ) {
if ( beam != NULL ) {
poisson_reflections(input, beam->adu_per_photon);
if ( have_adu_per_photon ) {
poisson_reflections(input, adu_per_photon);
} else {
ERROR("You must give a beam parameters file in order"
" to calculate Poisson noise.\n");
ERROR("You must give the number of ADU per photon to "
"use --poisson.\n");
return 1;
}
}
......
......@@ -115,6 +115,7 @@ int main(int argc, char *argv[])
det->panels[0].cny = -512.0-sep;
det->panels[0].clen = 100.0e-3;
det->panels[0].res = 9090.91;
det->panels[0].adu_per_eV = 1.0;
det->panels[1].min_fs = 0;
det->panels[1].max_fs = 1023;
......@@ -132,6 +133,7 @@ int main(int argc, char *argv[])
det->panels[1].cny = sep;
det->panels[1].clen = 100.0e-3;
det->panels[1].res = 9090.91;
det->panels[0].adu_per_eV = 1.0;
cpu_image.det = det;
gpu_image.det = det;
......@@ -142,8 +144,6 @@ int main(int argc, char *argv[])
beam->photon_energy = 9000.0;
beam->bandwidth = 0.1 / 100.0;
beam->divergence = 0.0;
beam->dqe = 1.0;
beam->adu_per_photon = 1.0;
cpu_image.beam = beam;
gpu_image.beam = beam;
......
......@@ -291,6 +291,8 @@ int main(int argc, char *argv[])
image.data = malloc(128*128*sizeof(float));
image.flags = NULL;
image.beam = NULL;
image.lambda = ph_eV_to_lambda(1000.0);
image.det = calloc(1, sizeof(struct detector));
image.det->n_panels = 1;
......@@ -313,6 +315,7 @@ int main(int argc, char *argv[])
image.det->panels[0].clen = 1.0;
image.det->panels[0].res = 1.0;
image.det->panels[0].integr_radius = 10.0;
image.det->panels[0].adu_per_eV = 1.0;
image.width = 128;
image.height = 128;
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment