pattern_sim: Add --background

4 files changed, 32 insertions, 8 deletions

diff --git a/doc/man/pattern_sim.1 b/doc/man/pattern_sim.1
index e34bf11d..2ce1f9f0 100644
--- a/doc/man/pattern_sim.1
+++ b/doc/man/pattern_sim.1
@@ -132,6 +132,11 @@ Include \fIn\fR samples from the spectrum in the calculation.
 .PD
 Use \fItype\fR of spectrum.  \fItype\fR can be one of \fBtophat\fR or \fBsase\fR.
 
+.PD 0
+.IP \fB--background=\fR\fIn\fR
+.PD
+Add \fIn\fR photons of Poisson-distributed background uniformly over the detector surface.
+
 
 .SH REFLECTION LISTS
 
diff --git a/libcrystfel/src/detector.c b/libcrystfel/src/detector.c
index 99205e67..257d1daa 100644
--- a/libcrystfel/src/detector.c
+++ b/libcrystfel/src/detector.c
@@ -312,7 +312,8 @@ double get_tt(struct image *image, double fs, double ss, int *err)
 }
 
 
-void record_image(struct image *image, int do_poisson, gsl_rng *rng)
+void record_image(struct image *image, int do_poisson, int background,
+                  gsl_rng *rng)
 {
 	int x, y;
 	double total_energy, energy_density;
@@ -345,6 +346,7 @@ void record_image(struct image *image, int do_poisson, gsl_rng *rng)
 		double pix_area, Lsq;
 		double xs, ys, rx, ry;
 		double dsq, proj_area;
+		float dval;
 		struct panel *p;
 
 		intensity = (double)image->data[x + image->width*y];
@@ -378,8 +380,13 @@ void record_image(struct image *image, int do_poisson, gsl_rng *rng)
 			counts = cf;
 		}
 
-		image->data[x + image->width*y] = counts * p->adu_per_eV
-		                               * ph_lambda_to_eV(image->lambda);
+		/* Number of photons in pixel */
+		dval = counts + poisson_noise(rng, background);
+
+		/* Convert to ADU */
+		dval *= p->adu_per_eV * ph_lambda_to_eV(image->lambda);
+
+		image->data[x + image->width*y] = dval;
 
 		/* Sanity checks */
 		if ( isinf(image->data[x+image->width*y]) ) n_inf2++;
diff --git a/libcrystfel/src/detector.h b/libcrystfel/src/detector.h
index eaa5618e..e44147eb 100644
--- a/libcrystfel/src/detector.h
+++ b/libcrystfel/src/detector.h
@@ -3,12 +3,12 @@
  *
  * Detector properties
  *
- * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY,
- *                  a research centre of the Helmholtz Association.
+ * Copyright © 2012-2014 Deutsches Elektronen-Synchrotron DESY,
+ *                       a research centre of the Helmholtz Association.
  * Copyright © 2012 Richard Kirian
  *
  * Authors:
- *   2009-2012 Thomas White <taw@physics.org>
+ *   2009-2014 Thomas White <taw@physics.org>
  *   2011-2012 Richard Kirian <rkirian@asu.edu>
  *   2011      Andrew Aquila
  *
@@ -153,7 +153,8 @@ extern double get_tt(struct image *image, double xs, double ys, int *err);
 
 extern int in_bad_region(struct detector *det, double fs, double ss);
 
-extern void record_image(struct image *image, int do_poisson, gsl_rng *rng);
+extern void record_image(struct image *image, int do_poisson, int background,
+                         gsl_rng *rng);
 
 extern struct panel *find_panel(struct detector *det, double fs, double ss);
 
diff --git a/src/pattern_sim.c b/src/pattern_sim.c
index d852d986..08fc75dd 100644
--- a/src/pattern_sim.c
+++ b/src/pattern_sim.c
@@ -89,6 +89,7 @@ static void show_help(const char *s)
 "     --no-noise            Do not calculate Poisson noise.\n"
 " -s, --sample-spectrum=<N> Use N samples from spectrum. Default 3.\n"
 " -x, --spectrum=<type>     Type of spectrum to simulate.\n"
+"     --background=<N>      Add N photons of Poisson background (default 0).\n"
 );
 }
 
@@ -249,6 +250,7 @@ int main(int argc, char *argv[])
 	SymOpList *sym;
 	int nsamples = 3;
 	gsl_rng *rng;
+	int background = 0;
 
 	/* Long options */
 	const struct option longopts[] = {
@@ -273,6 +275,7 @@ int main(int argc, char *argv[])
 		{"gpu-dev",            1, NULL,                2},
 		{"min-size",           1, NULL,                3},
 		{"max-size",           1, NULL,                4},
+		{"background",         1, NULL,                5},
 		{0, 0, NULL, 0}
 	};
 
@@ -364,6 +367,14 @@ int main(int argc, char *argv[])
 			random_size++;
 			break;
 
+			case 5 :
+			background = strtol(optarg, &rval, 10);
+			if ( *rval != '\0' ) {
+				ERROR("Invalid background level.\n");
+				return 1;
+			}
+			break;
+
 			case 0 :
 			break;
 
@@ -622,7 +633,7 @@ int main(int argc, char *argv[])
 			goto skip;
 		}
 
-		record_image(&image, !config_nonoise, rng);
+		record_image(&image, !config_nonoise, background, rng);
 
 		if ( powder_fn != NULL ) {