indexamajig: Add --wavelength-estimate

This is a more sensible non-indexer-specific and non-Xray-specific
replacement for --pinkIndexer-override-photon-energy

8 files changed, 64 insertions, 35 deletions

diff --git a/doc/man/indexamajig.1 b/doc/man/indexamajig.1
index 1cca53ba..a912e19a 100644
--- a/doc/man/indexamajig.1
+++ b/doc/man/indexamajig.1
@@ -385,6 +385,11 @@ Skip the prediction refinement step.  Usually this will decrease the quality of
 Check that most of the peaks can be accounted for by the indexing solution.  This usually increases the quality of the indexing solutions, but prevents "subtract and retry" multi-lattice indexing from working well.
 
 .PD 0
+.IP \fB--wavelength-estimate=\fIm\fR
+.PD
+Some indexing algorithms need to know the wavelength of the incident radiation in advance, e.g. to prepare an internal look-up table.  However, if the wavelength is taken from image headers, then the wavelength is not available at start-up.  In this case, you will be prompted to add this option to give an approximate wavelength, in metres.  A warning will be generated if the actual wavelength differs from this value by more than 10%.
+
+.PD 0
 .IP \fB--taketwo-member-threshold=\fIn\fR
 .IP \fB--taketwo-len-tolerance=\fIn\fR
 .IP \fB--taketwo-angle-tolerance=\fIn\fR
@@ -446,7 +451,6 @@ These set low-level parameters for the XGANDALF indexing algorithm.
 .IP \fB--pinkIndexer-no-check-indexed\fR
 .IP \fB--pinkIndexer-max-refinement-disbalance=\fIn\fR
 .IP \fB--pinkIndexer-override-bandwidth=\fIn\fR
-.IP \fB--pinkIndexer-override-photon-energy=\fIn\fR
 .IP \fB--pinkIndexer-override-visible-energy-range=\fImin-max\fR
 
 .PD
@@ -474,8 +478,6 @@ These set low-level parameters for the PinkIndexer indexing algorithm.
 .IP
 \fB--pinkIndexer-override-bandwidth=\fIn\fR Overrides the bandwidth in (delta energy)/(mean energy) to use for indexing (which usually is defined in the geometry file). Should be used together with \fB--pinkIndexer-override-photon-energy=\fIn\fR. Note: this option sets the borders for the spectrum, whereas the option in the geometry file sets the standard deviation for a Gaussian that fits the spectrum. Internally, the standard deviation is multiplied by 5 to compute the hard borders of the spectrum. I.e., photon_energy_bandwith = 0.005 in the geometry file equals --pinkIndexer-override-bandwidth=0.025. For monochromatic experiments usually --pinkIndexer-override-bandwidth=0.02 is sufficiently large. If in doubt, use --pinkIndexer-override-visible-energy-range.
 .IP
-\fB--pinkIndexer-override-photon-energy=\fIn\fR Overrides the mean energy in eV to use for indexing (which usually is defined in the geometry file). Should be used together with \fB--pinkIndexer-override-bandwidth=\fIn\fR
-.IP
 \fB--pinkIndexer-override-visible-energy-range=\fImin-max\fR Overrides photon energy and bandwidth according to a range of energies that have high enough intensity to produce visible Bragg spots on the detector. min and max range borders are separated by a minus sign (no whitespace).
 
 .SH INTEGRATION OPTIONS
diff --git a/libcrystfel/src/index.c b/libcrystfel/src/index.c
index 094ed3d7..c06dda26 100644
--- a/libcrystfel/src/index.c
+++ b/libcrystfel/src/index.c
@@ -68,6 +68,7 @@ struct _indexingprivate
 	IndexingFlags flags;
 	UnitCell *target_cell;
 	double tolerance[6];
+	double wavelength_estimate;
 
 	int n_methods;
 	IndexingMethod *methods;
@@ -234,6 +235,7 @@ static char *friendly_indexer_name(IndexingMethod m)
 
 static void *prepare_method(IndexingMethod *m, UnitCell *cell,
                             const DataTemplate *dtempl,
+                            double wavelength_estimate,
                             struct xgandalf_options *xgandalf_opts,
                             struct pinkIndexer_options* pinkIndexer_opts,
                             struct felix_options *felix_opts,
@@ -283,7 +285,7 @@ static void *prepare_method(IndexingMethod *m, UnitCell *cell,
 
 		case INDEXING_PINKINDEXER :
 		priv = pinkIndexer_prepare(m, cell, pinkIndexer_opts,
-		                           dtempl);
+		                           dtempl, wavelength_estimate);
 		break;
 
 		default :
@@ -358,6 +360,7 @@ IndexingMethod *parse_indexing_methods(const char *method_list,
 IndexingPrivate *setup_indexing(const char *method_list, UnitCell *cell,
                                 const DataTemplate *dtempl,
                                 float *tols, IndexingFlags flags,
+                                double wavelength_estimate,
                                 struct taketwo_options *ttopts,
                                 struct xgandalf_options *xgandalf_opts,
                                 struct pinkIndexer_options *pinkIndexer_opts,
@@ -427,6 +430,7 @@ IndexingPrivate *setup_indexing(const char *method_list, UnitCell *cell,
 
 		ipriv->engine_private[i] = prepare_method(&methods[i], cell,
 		                                          dtempl,
+		                                          wavelength_estimate,
 		                                          xgandalf_opts,
 		                                          pinkIndexer_opts,
 		                                          felix_opts,
@@ -469,6 +473,7 @@ IndexingPrivate *setup_indexing(const char *method_list, UnitCell *cell,
 	ipriv->methods = methods;
 	ipriv->n_methods = n;
 	ipriv->flags = flags;
+	ipriv->wavelength_estimate = wavelength_estimate;
 
 	if ( cell != NULL ) {
 		ipriv->target_cell = cell_new_from_cell(cell);
@@ -940,6 +945,18 @@ void index_pattern_3(struct image *image, IndexingPrivate *ipriv, int *ping,
 	image->crystals = NULL;
 	image->n_crystals = 0;
 
+	if ( !isnan(ipriv->wavelength_estimate) ) {
+		if ( !within_tolerance(image->lambda,
+		                          ipriv->wavelength_estimate,
+		                          10.0) )
+		{
+			   ERROR("WARNING: Wavelength for %s %s (%e) differs by "
+			         "more than 10%% from estimated value (%e)\n",
+			         image->filename, image->ev,
+			         image->lambda, ipriv->wavelength_estimate);
+		}
+	}
+
 	orig = image->features;
 
 	for ( n=0; n<ipriv->n_methods; n++ ) {
diff --git a/libcrystfel/src/index.h b/libcrystfel/src/index.h
index 3edfc4cd..708c04a4 100644
--- a/libcrystfel/src/index.h
+++ b/libcrystfel/src/index.h
@@ -9,7 +9,7 @@
  * Copyright © 2012 Lorenzo Galli
  *
  * Authors:
- *   2010-2020 Thomas White <taw@physics.org>
+ *   2010-2021 Thomas White <taw@physics.org>
  *   2010      Richard Kirian
  *   2012      Lorenzo Galli
  *   2015      Kenneth Beyerlein <kenneth.beyerlein@desy.de>
@@ -168,6 +168,7 @@ extern IndexingPrivate *setup_indexing(const char *methods,
                                        const DataTemplate *dtempl,
                                        float *ltl,
                                        IndexingFlags flags,
+                                       double wavelength_estimate,
                                        struct taketwo_options *ttopts,
                                        struct xgandalf_options *xgandalf_opts,
                                        struct pinkIndexer_options *pinkIndexer_opts,
diff --git a/libcrystfel/src/indexers/pinkindexer.c b/libcrystfel/src/indexers/pinkindexer.c
index ac4a42d2..ae4e5642 100644
--- a/libcrystfel/src/indexers/pinkindexer.c
+++ b/libcrystfel/src/indexers/pinkindexer.c
@@ -3,11 +3,12 @@
  *
  * Interface to PinkIndexer
  *
- * Copyright © 2017-2020 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2017-2021 Deutsches Elektronen-Synchrotron DESY,
  *                       a research centre of the Helmholtz Association.
  *
  * Authors:
  *   2017-2019 Yaroslav Gevorkov <yaroslav.gevorkov@desy.de>
+ *   2021 Thomas White <thomas.white@desy.de>
  *
  * This file is part of CrystFEL.
  *
@@ -52,7 +53,6 @@ struct pinkIndexer_options {
 	int min_peaks;
 	int no_check_indexed;
 	float reflectionRadius; /* In m^-1 */
-	float customPhotonEnergy;
 	float customBandwidth;
 	float maxRefinementDisbalance;
 };
@@ -96,6 +96,7 @@ int run_pinkIndexer(struct image *image, void *ipriv)
 	reciprocalPeaks_1_per_A_t* reciprocalPeaks_1_per_A = &(pinkIndexer_private_data->reciprocalPeaks_1_per_A);
 	float *intensities = pinkIndexer_private_data->intensities;
 
+	/* FIXME: Check if wavelength is too far from estimate */
 	int peakCountMax = image_feature_count(image->features);
 	if (peakCountMax < 5) {
 		int goodLatticesCount = 0;
@@ -195,17 +196,17 @@ static int want_center_adjustment(struct pinkIndexer_options *pinkIndexer_opts)
 
 void *pinkIndexer_prepare(IndexingMethod *indm, UnitCell *cell,
                           struct pinkIndexer_options *pinkIndexer_opts,
-                          const DataTemplate *dtempl)
+                          const DataTemplate *dtempl,
+                          double wavelength_estimate)
 {
-	if ( !data_template_has_fixed_wavelength(dtempl)
-	     && (pinkIndexer_opts->customPhotonEnergy <= 0.0) )
-	{
-		ERROR("Geometry file refers to image metadata for "
-		      "wavelength.\n");
-		ERROR("To use PinkIndexer, specify a constant "
-		      "wavelength in the geometry file, or use "
-		      "--pinkIndexer-override-photon-energy.\n");
+	float beamEenergy_eV;
+
+	if ( isnan(wavelength_estimate) ) {
+		ERROR("PinkIndexer requires a wavelength estimate.  "
+		      "Try again with --wavelength-estimate=xx\n");
 		return NULL;
+	} else {
+		beamEenergy_eV = J_to_eV(ph_lambda_to_en(wavelength_estimate));
 	}
 
 	if ( !data_template_has_fixed_geometry(dtempl)
@@ -261,11 +262,7 @@ void *pinkIndexer_prepare(IndexingMethod *indm, UnitCell *cell,
 	}
 
 	/* FIXME: Beam gone */
-	float beamEenergy_eV = beam->photon_energy;
 	float nonMonochromaticity = beam->bandwidth*5;
-	if(pinkIndexer_opts->customPhotonEnergy > 0){
-		beamEenergy_eV = pinkIndexer_opts->customPhotonEnergy;
-	}
 	if(pinkIndexer_opts->customBandwidth >= 0){
 		nonMonochromaticity = pinkIndexer_opts->customBandwidth;
 	}
@@ -462,8 +459,6 @@ static void pinkIndexer_show_help()
 "     --pinkIndexer-max-refinement-disbalance=n\n"
 "                           Maximum disbalance after refinement:\n"
 "                            0 (no disbalance) to 2 (extreme disbalance), default 0.4\n"
-"     --pinkIndexer-override-photon-energy=ev\n"
-"                           Mean energy in eV to use for indexing.\n"
 "     --pinkIndexer-override-bandwidth=n\n"
 "                           Bandwidth in (delta energy)/(mean energy) to use for indexing.\n"
 "     --pinkIndexer-override-visible-energy-range=min-max\n"
@@ -492,7 +487,6 @@ int pinkIndexer_default_options(PinkIndexerOptions **opts_ptr)
 	opts->no_check_indexed = 0;
 	opts->min_peaks = 2;
 	opts->reflectionRadius = -1;
-	opts->customPhotonEnergy = -1;
 	opts->customBandwidth = -1;
 	opts->maxRefinementDisbalance = 0.4;
 
@@ -588,12 +582,9 @@ static error_t pinkindexer_parse_arg(int key, char *arg,
 		break;
 
 		case 11 :
-		if (sscanf(arg, "%f", &(*opts_ptr)->customPhotonEnergy) != 1)
-		{
-			ERROR("Invalid value for --pinkIndexer-override-photon-energy\n");
-			return EINVAL;
-		}
-		break;
+		ERROR("Please use --wavelength-estimate instead of "
+		      "--pinkIndexer-override-photon-energy.\n");
+		return EINVAL;
 
 		case 12 :
 		if (sscanf(arg, "%f", &(*opts_ptr)->customBandwidth) != 1)
diff --git a/libcrystfel/src/indexers/pinkindexer.h b/libcrystfel/src/indexers/pinkindexer.h
index 57df5f62..fc2714c8 100644
--- a/libcrystfel/src/indexers/pinkindexer.h
+++ b/libcrystfel/src/indexers/pinkindexer.h
@@ -3,11 +3,12 @@
  *
  * Interface to PinkIndexer
  *
- * Copyright © 2017-2020 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2017-2021 Deutsches Elektronen-Synchrotron DESY,
  *                       a research centre of the Helmholtz Association.
  *
  * Authors:
  *   2017-2019 Yaroslav Gevorkov <yaroslav.gevorkov@desy.de>
+ *   2021 Thomas White <thomas.white@desy.de>
  *
  * This file is part of CrystFEL.
  *
@@ -40,7 +41,8 @@ extern int run_pinkIndexer(struct image *image, void *ipriv);
 
 extern void *pinkIndexer_prepare(IndexingMethod *indm, UnitCell *cell,
                                  struct pinkIndexer_options *pinkIndexer_opts,
-                                 const DataTemplate *dtempl);
+                                 const DataTemplate *dtempl,
+                                 double wavelength_estimate);
 
 extern void pinkIndexer_cleanup(void *pp);
 
diff --git a/src/gui_index.c b/src/gui_index.c
index 85e36f36..f244f1f5 100644
--- a/src/gui_index.c
+++ b/src/gui_index.c
@@ -516,9 +516,10 @@ static void run_indexing_once(struct crystfelproject *proj)
 	                       &pinkIndexer_opts,
 	                       &felix_opts);
 
-	ipriv = setup_indexing(methods, cell, proj->dtempl,
+	ipriv = setup_indexing(methods, cell,
 	                       proj->indexing_params.tols,
 	                       indexing_flags(&proj->indexing_params),
+	                       proj->cur_image->lambda,
 	                       taketwoopts, xgandalf_opts,
 	                       pinkIndexer_opts, felix_opts);
 	free(methods);
diff --git a/src/indexamajig.c b/src/indexamajig.c
index 3780b33a..f33d3ec8 100644
--- a/src/indexamajig.c
+++ b/src/indexamajig.c
@@ -9,7 +9,7 @@
  * Copyright © 2012 Lorenzo Galli
  *
  * Authors:
- *   2010-2020 Thomas White <taw@physics.org>
+ *   2010-2021 Thomas White <taw@physics.org>
  *   2011      Richard Kirian
  *   2012      Lorenzo Galli
  *   2012      Chunhong Yoon
@@ -455,6 +455,14 @@ static error_t parse_arg(int key, char *arg, struct argp_state *state)
 		 *  better than the user expected. */
 		break;
 
+		case 413 :
+		if (sscanf(arg, "%f", &args->iargs.wavelength_estimate) != 1)
+		{
+			ERROR("Invalid value for --wavelength-estimate\n");
+			return EINVAL;
+		}
+		break;
+
 		/* ---------- Integration ---------- */
 
 		case 501 :
@@ -666,6 +674,7 @@ int main(int argc, char *argv[])
 	args.iargs.fix_divergence = -1.0;
 	args.iargs.no_image_data = 0;
 	args.iargs.no_mask_data = 0;
+	args.iargs.wavelength_estimate = NAN;
 
 	argp_program_version_hook = show_version;
 
@@ -763,6 +772,8 @@ int main(int argc, char *argv[])
 		        "accounted for by the indexing solution"},
 		{"check-peaks", 411, NULL, OPTION_HIDDEN, NULL},
 		{"no-cell-combinations", 412, NULL, OPTION_HIDDEN, NULL},
+		{"wavelength-estimate", 413, "metres", 0,
+		        "Estimate of the incident radiation wavelength, in metres."},
 
 		{NULL, 0, 0, OPTION_DOC, "Integration options:", 5},
 		{"integration", 501, "method", OPTION_NO_USAGE, "Integration method"},
@@ -960,9 +971,12 @@ int main(int argc, char *argv[])
 			flags |= INDEXING_RETRY;
 		}
 
-		args.iargs.ipriv = setup_indexing(args.indm_str, args.iargs.cell,
+		args.iargs.ipriv = setup_indexing(args.indm_str,
+		                                  args.iargs.cell,
 		                                  args.iargs.dtempl,
-		                                  args.iargs.tols, flags,
+		                                  args.iargs.tols,
+		                                  flags,
+		                                  args.iargs.wavelength_estimate,
 		                                  taketwo_opts,
 		                                  xgandalf_opts,
 		                                  pinkindexer_opts,
diff --git a/src/process_image.h b/src/process_image.h
index e64118dd..8174f73d 100644
--- a/src/process_image.h
+++ b/src/process_image.h
@@ -98,6 +98,7 @@ struct index_args
 	signed int wait_for_file; /* -1 means wait forever */
 	int no_image_data;
 	int no_mask_data;
+	float wavelength_estimate;
 };