Add XTC wavelength calibration using the FEE spectrometer

Adds two parameters, spectrum_eV_per_pixel and spectrum_eV_offset. If present in the XTC locator, these parameters are used to deterimine the energy using a weighted average of the horizontal projection of the 2D FEE spectrometer.

format/FormatXTC.py

@@ -12,6 +12,9 @@
 from dxtbx.format.FormatStill import FormatStill
 import time
 
+import numpy as np
+from cctbx import factor_ev_angstrom
+
 try:
     import psana
 
@@ -48,6 +51,16 @@
   wavelength_offset = None
     .type = float
     .help = Optional constant shift to apply to each wavelength
+  spectrum_eV_per_pixel = None
+    .type = float
+    .help = If not None, use the FEE spectrometer to determine the wavelength. \
+            spectrum_eV_offset should be also specified. A weighted average of \
+            the horizontal projection of the per-shot FEE spectrometer is used.\
+            The equation for each pixel eV is \
+            eV = (spectrum_eV_per_pixel * pixel_number) + spectrum_eV_offset
+  spectrum_eV_offset = None
+    .type = float
+    .help = See spectrum_eV_per_pixel
 """
 locator_scope = parse(locator_str)
 
@@ -87,6 +100,7 @@ def __init__(self, image_file, **kwargs):
         self._beam_index = None
         self._beam_cache = None
         self._initialized = True
+        self._fee = None
 
     @staticmethod
     def understand(image_file):
@@ -251,7 +265,20 @@ def _beam(self, index=None):
         if self._beam_index != index:
             self._beam_index = index
             evt = self._get_event(index)
-            wavelength = cspad_tbx.evt_wavelength(evt)
+            if self.params.spectrum_eV_per_pixel is not None:
+                if self._fee is None:
+                    self._fee = psana.Detector("FEE-SPEC0")
+                fee = self._fee.get(evt)
+                if fee is None:
+                    wavelength = cspad_tbx.evt_wavelength(evt)
+                else:
+                    x = (
+                        self.params.spectrum_eV_per_pixel
+                        * np.array(range(len(fee.hproj())))
+                    ) + self.params.spectrum_eV_offset
+                    wavelength = factor_ev_angstrom / np.average(x, weights=fee.hproj())
+            else:
+                wavelength = cspad_tbx.evt_wavelength(evt)
             if wavelength is None:
                 self._beam_cache = None
             else: