Add Lorentz factor

5 files changed, 76 insertions, 13 deletions

diff --git a/doc/reference/libcrystfel/CrystFEL-sections.txt b/doc/reference/libcrystfel/CrystFEL-sections.txt
index ff8001cf..09267538 100644
--- a/doc/reference/libcrystfel/CrystFEL-sections.txt
+++ b/doc/reference/libcrystfel/CrystFEL-sections.txt
@@ -21,10 +21,11 @@ next_found_refl
 get_excitation_error
 get_detector_pos
 get_partiality
+get_lorentz
+get_partial
 get_indices
 get_symmetric_indices
 get_intensity
-get_partial
 get_scalable
 get_refinable
 get_redundancy
@@ -34,6 +35,8 @@ get_temp1
 get_temp2
 <SUBSECTION>
 set_detector_pos
+set_partiality
+set_lorentz
 set_partial
 set_intensity
 set_scalable
diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c
index 450ebfaf..7a553f82 100644
--- a/libcrystfel/src/geometry.c
+++ b/libcrystfel/src/geometry.c
@@ -134,6 +134,10 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
 	Reflection *refl;
 	double cet, cez;
 	double pr;
+	double L;
+
+	/* Don't predict 000 */
+	if ( abs(h)+abs(k)+abs(l) == 0 ) return NULL;
 
 	pr = crystal_get_profile_radius(cryst);
 
@@ -160,6 +164,19 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
 	     && (fabs(rlow) > pr)
 	     && (fabs(rhigh) > pr) ) return NULL;
 
+	if ( rlow < rhigh ) {
+		ERROR("Reflection with rlow < rhigh!\n");
+		ERROR("%3i %3i %3i  rlow = %e, rhigh = %e\n",
+		      h, k, l, rlow, rhigh);
+		ERROR("div = %e\n", image->div);
+		return NULL;
+	}
+
+	/* Lorentz factor is determined direction from the r values, before
+	 * clamping.  The multiplication by the profile radius is to make the
+	 * correction factor vaguely near 1. */
+	L = pr / (rlow - rhigh);
+
 	/* If the "lower" Ewald sphere is a long way away, use the
 	 * position at which the Ewald sphere would just touch the
 	 * reflection.
@@ -185,7 +202,6 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
 		rhigh = +pr;
 		clamp_high = +1;
 	}
-	assert(clamp_low >= clamp_high);
 
 	/* Calculate partiality */
 	part = partiality(rlow, rhigh, pr);
@@ -198,6 +214,7 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
 	refl = reflection_new(h, k, l);
 	set_detector_pos(refl, 0.0, xda, yda);
 	set_partial(refl, rlow, rhigh, part, clamp_low, clamp_high);
+	set_lorentz(refl, L);
 	set_symmetric_indices(refl, h, k, l);
 	set_redundancy(refl, 1);
 
diff --git a/libcrystfel/src/reflist.c b/libcrystfel/src/reflist.c
index 64154774..b3b9f85b 100644
--- a/libcrystfel/src/reflist.c
+++ b/libcrystfel/src/reflist.c
@@ -3,11 +3,11 @@
  *
  * Fast reflection/peak list
  *
- * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY,
- *                  a research centre of the Helmholtz Association.
+ * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY,
+ *                       a research centre of the Helmholtz Association.
  *
  * Authors:
- *   2011-2012 Thomas White <taw@physics.org>
+ *   2011-2013 Thomas White <taw@physics.org>
  *
  * This file is part of CrystFEL.
  *
@@ -70,6 +70,7 @@ struct _refldata {
 	double r1;  /* First excitation error */
 	double r2;  /* Second excitation error */
 	double p;   /* Partiality */
+	double L;   /* Lorentz factor */
 	int clamp1; /* Clamp status for r1 */
 	int clamp2; /* Clamp status for r2 */
 
@@ -388,7 +389,7 @@ void get_symmetric_indices(const Reflection *refl,
  * get_partiality:
  * @refl: A %Reflection
  *
- * Returns: The partiality of the reflection.
+ * Returns: The partiality of the reflection.  See get_lorentz().
  **/
 double get_partiality(const Reflection *refl)
 {
@@ -397,6 +398,19 @@ double get_partiality(const Reflection *refl)
 
 
 /**
+ * get_lorentz:
+ * @refl: A %Reflection
+ *
+ * Returns: The Lorentz factor for the reflection.  To "scale up" a partial
+ * reflection, divide by this multiplied by the partiality.
+ **/
+double get_lorentz(const Reflection *refl)
+{
+	return refl->data.L;
+}
+
+
+/**
  * get_intensity:
  * @refl: A %Reflection
  *
@@ -602,10 +616,35 @@ void set_partial(Reflection *refl, double r1, double r2, double p,
 /**
  * set_intensity:
  * @refl: A %Reflection
+ * @p: The partiality for the reflection.
+ *
+ * Set the partiality for the reflection.  See set_lorentz().
+ **/
+void set_partiality(Reflection *refl, double p)
+{
+	refl->data.p = p;
+}
+
+/**
+ * set_lorentz:
+ * @refl: A %Reflection
+ * @L: The Lorentz factor for the reflection.
+ *
+ * Set the Lorentz factor for the reflection.  To "scale up" a partial
+ * reflection, divide by this multiplied by the partiality.
+ **/
+void set_lorentz(Reflection *refl, double L)
+{
+	refl->data.L = L;
+}
+
+
+/**
+ * set_intensity:
+ * @refl: A %Reflection
  * @intensity: The intensity for the reflection.
  *
- * Set the intensity for the reflection.  Note that retrieval is done with
- * get_intensity().
+ * Set the intensity for the reflection.
  **/
 void set_intensity(Reflection *refl, double intensity)
 {
diff --git a/libcrystfel/src/reflist.h b/libcrystfel/src/reflist.h
index 246ef885..7969235c 100644
--- a/libcrystfel/src/reflist.h
+++ b/libcrystfel/src/reflist.h
@@ -3,11 +3,11 @@
  *
  * Fast reflection/peak list
  *
- * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY,
- *                  a research centre of the Helmholtz Association.
+ * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY,
+ *                       a research centre of the Helmholtz Association.
  *
  * Authors:
- *   2011-2012 Thomas White <taw@physics.org>
+ *   2011-2013 Thomas White <taw@physics.org>
  *
  * This file is part of CrystFEL.
  *
@@ -78,6 +78,7 @@ extern Reflection *next_found_refl(Reflection *refl);
 extern double get_excitation_error(const Reflection *refl);
 extern void get_detector_pos(const Reflection *refl, double *fs, double *ss);
 extern double get_partiality(const Reflection *refl);
+extern double get_lorentz(const Reflection *refl);
 extern void get_indices(const Reflection *refl,
                         signed int *h, signed int *k, signed int *l);
 extern void get_symmetric_indices(const Reflection *refl,
@@ -100,6 +101,8 @@ extern void set_detector_pos(Reflection *refl, double exerr,
                              double fs, double ss);
 extern void set_partial(Reflection *refl, double r1, double r2, double p,
                         double clamp_low, double clamp_high);
+extern void set_partiality(Reflection *refl, double p);
+extern void set_lorentz(Reflection *refl, double L);
 extern void set_intensity(Reflection *refl, double intensity);
 extern void set_scalable(Reflection *refl, int scalable);
 extern void set_refinable(Reflection *refl, int refinable);
diff --git a/src/partial_sim.c b/src/partial_sim.c
index 7c93da1b..32e6666a 100644
--- a/src/partial_sim.c
+++ b/src/partial_sim.c
@@ -100,12 +100,13 @@ static void calculate_partials(Crystal *cr,
 	{
 		signed int h, k, l;
 		Reflection *rfull;
-		double p, Ip, If;
+		double L, p, Ip, If;
 		int bin;
 
 		get_indices(refl, &h, &k, &l);
 		get_asymm(sym, h, k, l, &h, &k, &l);
 		p = get_partiality(refl);
+		L = get_lorentz(refl);
 
 		pthread_mutex_lock(full_lock);
 		rfull = find_refl(full, h, k, l);
@@ -139,7 +140,7 @@ static void calculate_partials(Crystal *cr,
 			}
 		}
 
-		Ip = crystal_get_osf(cr) * p * If;
+		Ip = crystal_get_osf(cr) * L * p * If;
 
 		res = resolution(crystal_get_cell(cr), h, k, l);
 		bin = NBINS*2.0*res/max_q;