Beam stats

sirepo_bluesky/shadow_handler.py

@@ -1,7 +1,12 @@
+import contextlib
+import os
+
 import numpy as np
 import Shadow.ShadowLibExtensions as sd
 import Shadow.ShadowTools
 
+from . import utils
+
 
 def read_shadow_file_col(filename, parameter=30):
     """Read specified parameter from the Shadow3 output binary file.
@@ -49,7 +54,9 @@ def read_shadow_file_col(filename, parameter=30):
           32   S2-stokes = 2 |Es| |Ep| cos(phase_s-phase_p)
           33   S3-stokes = 2 |Es| |Ep| sin(phase_s-phase_p)
     """
-    data = Shadow.ShadowTools.getshcol(filename, col=parameter)
+    with open(os.devnull, "w") as devnull:
+        with contextlib.redirect_stdout(devnull):
+            data = Shadow.ShadowTools.getshcol(filename, col=parameter)
 
     mean_value = np.mean(data)
 
@@ -73,24 +80,35 @@ def read_shadow_file(filename, histogram_bins=None):
     beam.load(filename)
 
     # 1=X spatial coordinate; 3=Z spatial coordinate
-    data_dict = beam.histo2(1, 3, nolost=1, nbins=histogram_bins)
-    data = data_dict["histogram"]
+    with open(os.devnull, "w") as devnull:
+        with contextlib.redirect_stdout(devnull):
+            data_dict = beam.histo2(1, 3, nolost=1, nbins=histogram_bins)
 
-    # This returns a list of N values (N=number of rays)
-    photon_energy_list = Shadow.ShadowTools.getshcol(filename, col=11)  # 11=Energy [eV]
+            # This returns a list of N values (N=number of rays)
+            photon_energy_list = Shadow.ShadowTools.getshcol(filename, col=11)  # 11=Energy [eV]
+
+    data = data_dict["histogram"]
     photon_energy = np.mean(photon_energy_list)
 
-    return {
+    # convert to um
+    horizontal_extent = 1e3 * np.array(data_dict["xrange"][:2])
+    vertical_extent = 1e3 * np.array(data_dict["yrange"][:2])
+
+    ret = {
         "data": data,
         "shape": data.shape,
-        "mean": np.mean(data),
+        "flux": data.sum(),
+        "mean": data.mean(),
         "photon_energy": photon_energy,
-        "horizontal_extent": data_dict["xrange"][:2],
-        "vertical_extent": data_dict["yrange"][:2],
-        # 'labels': labels,
-        # 'units': units,
+        "horizontal_extent": horizontal_extent,
+        "vertical_extent": vertical_extent,
+        "units": "um",
     }
 
+    ret.update(utils.get_beam_stats(data, horizontal_extent, vertical_extent))
+
+    return ret
+
 
 class ShadowFileHandler:
     specs = {"shadow"}

sirepo_bluesky/sirepo_ophyd.py

@@ -77,7 +77,12 @@ def trigger(self, *args, **kwargs):
 class SirepoWatchpoint(DeviceWithJSONData):
     image = Cpt(ExternalFileReference, kind="normal")
     shape = Cpt(Signal)
-    mean = Cpt(Signal, kind="hinted")
+    flux = Cpt(Signal, kind="hinted")
+    mean = Cpt(Signal, kind="normal")
+    x = Cpt(Signal, kind="normal")
+    y = Cpt(Signal, kind="normal")
+    fwhm_x = Cpt(Signal, kind="hinted")
+    fwhm_y = Cpt(Signal, kind="hinted")
     photon_energy = Cpt(Signal, kind="normal")
     horizontal_extent = Cpt(Signal)
     vertical_extent = Cpt(Signal)
@@ -152,7 +157,12 @@ def trigger(self, *args, **kwargs):
 
         def update_components(_data):
             self.shape.put(_data["shape"])
+            self.flux.put(_data["flux"])
             self.mean.put(_data["mean"])
+            self.x.put(_data["x"])
+            self.y.put(_data["y"])
+            self.fwhm_x.put(_data["fwhm_x"])
+            self.fwhm_y.put(_data["fwhm_y"])
             self.photon_energy.put(_data["photon_energy"])
             self.horizontal_extent.put(_data["horizontal_extent"])
             self.vertical_extent.put(_data["vertical_extent"])
@@ -233,7 +243,12 @@ def trigger(self, *args, **kwargs):
 
         def update_components(_data):
             self.shape.put(_data["shape"])
+            self.flux.put(_data["flux"])
             self.mean.put(_data["mean"])
+            self.x.put(_data["x"])
+            self.y.put(_data["y"])
+            self.fwhm_x.put(_data["fwhm_x"])
+            self.fwhm_y.put(_data["fwhm_y"])
             self.photon_energy.put(_data["photon_energy"])
             self.horizontal_extent.put(_data["horizontal_extent"])
             self.vertical_extent.put(_data["vertical_extent"])
@@ -401,7 +416,7 @@ def create_classes(sirepo_data, connection, create_objects=True, extra_model_fie
 
                 components[k] = Cpt(
                     cpt_class,
-                    value=v,
+                    value=(float(v) if type(v) is int else v),
                     sirepo_dict=sirepo_dict,
                     sirepo_param=k,
                 )

sirepo_bluesky/srw_handler.py

@@ -1,6 +1,8 @@
 import numpy as np
 import srwpy.uti_plot_com as srw_io
 
+from . import utils
+
 
 def read_srw_file(filename, ndim=2):
     data, mode, ranges, labels, units = srw_io.file_load(filename)
@@ -13,18 +15,28 @@ def read_srw_file(filename, ndim=2):
     else:
         raise ValueError(f"The value ndim={ndim} is not supported.")
 
-    return {
+    horizontal_extent = np.array(ranges[3:5])
+    vertical_extent = np.array(ranges[6:8])
+
+    ret = {
         "data": data,
         "shape": data.shape,
-        "mean": np.mean(data),
+        "flux": data.sum(),
+        "mean": data.mean(),
         "photon_energy": photon_energy,
-        "horizontal_extent": ranges[3:5],
-        "vertical_extent": ranges[6:8],
-        # 'mode': mode,
+        "horizontal_extent": horizontal_extent,
+        "vertical_extent": vertical_extent,
         "labels": labels,
-        "units": units,
+        "units": "units",
     }
 
+    if ndim == 1:
+        ret.update({key: np.nan for key in ["x", "y", "fwhm_x", "fwhm_y"]})
+    if ndim == 2:
+        ret.update(utils.get_beam_stats(data, horizontal_extent, vertical_extent))
+
+    return ret
+
 
 class SRWFileHandler:
     specs = {"srw"}

sirepo_bluesky/utils/__init__.py

@@ -0,0 +1,27 @@
+import numpy as np
+
+
+def get_beam_stats(image, x_extent, y_extent):
+    n_y, n_x = image.shape
+
+    if image.sum() > 0:
+        X, Y = np.meshgrid(np.linspace(*x_extent, n_x), np.linspace(*y_extent, n_y))
+
+        mean_x = np.sum(X * image) / np.sum(image)
+        mean_y = np.sum(Y * image) / np.sum(image)
+
+        sigma_x = np.sqrt(np.sum((X - mean_x) ** 2 * image) / np.sum(image))
+        sigma_y = np.sqrt(np.sum((Y - mean_y) ** 2 * image) / np.sum(image))
+
+    else:
+        mean_x, mean_y, sigma_x, sigma_y = np.nan, np.nan, np.nan, np.nan
+
+    return {
+        "shape": (n_y, n_x),
+        "flux": image.sum(),
+        "mean": image.mean(),
+        "x": mean_x,
+        "y": mean_y,
+        "fwhm_x": 2 * np.sqrt(2 * np.log(2)) * sigma_x,
+        "fwhm_y": 2 * np.sqrt(2 * np.log(2)) * sigma_y,
+    }