From 6326b472aca7193b4bc1056e672074c1a44df9f8 Mon Sep 17 00:00:00 2001
From: ahmedk <karim.ahmed@xfel.eu>
Date: Wed, 13 Sep 2023 18:06:43 +0200
Subject: [PATCH] Add rel_gain_mode and Start adding first steps toward using
 CS constants

---
 .../AGIPD/AGIPD_Correct_and_Verify.ipynb      | 32 +++++---
 src/cal_tools/agipdlib.py                     | 81 ++++++-------------
 src/cal_tools/agipdutils.py                   | 76 +++++++++++++++++
 src/cal_tools/calcat_interface.py             |  4 +-
 4 files changed, 126 insertions(+), 67 deletions(-)

diff --git a/notebooks/AGIPD/AGIPD_Correct_and_Verify.ipynb b/notebooks/AGIPD/AGIPD_Correct_and_Verify.ipynb
index cf85b3e52..bca36cae0 100644
--- a/notebooks/AGIPD/AGIPD_Correct_and_Verify.ipynb
+++ b/notebooks/AGIPD/AGIPD_Correct_and_Verify.ipynb
@@ -66,7 +66,8 @@
     "\n",
     "# Correction Booleans\n",
     "only_offset = False # Apply only Offset correction. if False, Offset is applied by Default. if True, Offset is only applied.\n",
-    "rel_gain = False # do relative gain correction based on PC data\n",
+    "rel_gain_mode = \"CS\"  # Select relative gain correct mode. Choices (PC: Pulse Capacitor or CS: Current Source). Default CS.\n",
+    "rel_gain = False # do relative gain correction based on PC or CS data\n",
     "xray_gain = False # do relative gain correction based on xray data\n",
     "blc_noise = False # if set, baseline correction via noise peak location is attempted\n",
     "blc_stripes = False # if set, baseline corrected via stripes\n",
@@ -623,22 +624,35 @@
     "agipd_metadata = agipd_cal.metadata(dark_constants)\n",
     "\n",
     "agipd_cal.gain_mode = None  # gain_mode is not used for gain constants\n",
-    "pc_constants, ff_constants = [], []\n",
-    "if any(agipd_corr.pc_bools):\n",
-    "    pc_constants = [\"SlopesPC\", \"BadPixelsPC\"]\n",
-    "    get_constants_and_update_metadata(\n",
-    "        agipd_cal, agipd_metadata, pc_constants)\n",
+    "pc_constants, ff_constants, cs_constants = [], [], []\n",
     "\n",
     "if agipd_corr.corr_bools.get('xray_corr'):\n",
     "    ff_constants = list(agipd_cal.illuminated_calibrations)\n",
     "    get_constants_and_update_metadata(\n",
     "        agipd_cal, agipd_metadata, ff_constants)\n",
     "\n",
+    "if any(agipd_corr.gain_bools):\n",
+    "    if not rel_gain_mode in [\"CS\", \"PC\"]:\n",
+    "        raise ValueError(\n",
+    "            \"Selected `rel_gain_mode` is unexpected. \"\n",
+    "            \"Please select between CS or PC.\")\n",
+    "\n",
+    "    if rel_gain_mode == \"cs\":\n",
+    "        # Integration time is not used with CS\n",
+    "        agipd_cal.integration_time = None\n",
+    "        cs_constants = [\"SlopesCS\", \"BadPixelsCS\"]\n",
+    "        get_constants_and_update_metadata(\n",
+    "            agipd_cal, agipd_metadata, pc_constants)\n",
+    "    else:  # \"pc\" or \"off\"\n",
+    "        pc_constants = [\"SlopesPC\", \"BadPixelsPC\"]\n",
+    "        get_constants_and_update_metadata(\n",
+    "            agipd_cal, agipd_metadata, pc_constants)\n",
+    "\n",
     "step_timer.done_step(\"Constants were retrieved in\")\n",
     "\n",
     "print(\"Preparing constants (\"\n",
     "      f\"FF: {agipd_corr.corr_bools.get('xray_corr', False)}, \"\n",
-    "      f\"PC: {any(agipd_corr.pc_bools)}, \"\n",
+    "      f\"{rel_gain_mode}: {any(agipd_corr.gain_bools)}, \"\n",
     "      f\"BLC: {any(agipd_corr.blc_bools)})\")\n",
     "# Display retrieved calibration constants timestamps\n",
     "agipd_cal.display_markdown_retrieved_constants(metadata=agipd_metadata)"
@@ -731,7 +745,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# Store timestamps for Offset, SlopesPC, and SlopesFF\n",
+    "# Store timestamps for Offset, SlopesPC/SlopesCS, and SlopesFF\n",
     "# in YAML file for time-summary table.\n",
     "timestamps = {}\n",
     "\n",
@@ -742,7 +756,7 @@
     "\n",
     "    # Store few time stamps if exists\n",
     "    # Add NA to keep array structure\n",
-    "    for key in ['Offset', 'SlopesPC', 'SlopesFF']:\n",
+    "    for key in ['Offset', f'Slopes{rel_gain_mode}', 'SlopesFF']:\n",
     "        if key in mod_mdata:\n",
     "            module_timestamps[key] = mod_mdata[key][\"begin_validity_at\"]\n",
     "        else:\n",
diff --git a/src/cal_tools/agipdlib.py b/src/cal_tools/agipdlib.py
index 74c3b852b..f01ea7221 100644
--- a/src/cal_tools/agipdlib.py
+++ b/src/cal_tools/agipdlib.py
@@ -21,6 +21,7 @@ from cal_tools.agipdutils import (
     cast_array_inplace,
     correct_baseline_via_hist,
     correct_baseline_via_hist_asic,
+    get_gain_pc_slopes,
     make_noisy_adc_mask,
     match_asic_borders,
     melt_snowy_pixels,
@@ -632,8 +633,9 @@ class AgipdCorrections:
             raise Exception('Correction Booleans: '
                             f'{bools} are not available!')
 
+        self.rel_gain_mode = "CS"
         # Flags allowing for pulse capacitor constant retrieval.
-        self.pc_bools = [self.corr_bools.get("rel_gain"),
+        self.gain_bools = [self.corr_bools.get("rel_gain"),
                          self.corr_bools.get("adjust_mg_baseline"),
                          self.corr_bools.get('blc_noise'),
                          self.corr_bools.get('blc_hmatch'),
@@ -1405,16 +1407,22 @@ class AgipdCorrections:
             self.xray_cor[module_idx][...] = xray_cor.transpose()[...]
 
         # add additional bad pixel information
-        if any(self.pc_bools):
-            if "BadPixelsPC" in cons_data:
-                bppc = np.moveaxis(cons_data["BadPixelsPC"].astype(np.uint32),
-                                   0, 2)
-                bpixels |= bppc[..., :bpixels.shape[2], None]
+        if any(self.gain_bools):
+            if f"BadPixels{self.rel_gain_mode}" in cons_data:
+                bp_relgain = np.moveaxis(
+                    cons_data[f"BadPixels{self.rel_gain_mode}"].astype(np.uint32), 0, 2)  # noqa
+                bpixels |= bp_relgain[..., :bpixels.shape[2], None]
 
             # calculate relative gain from the constants
             rel_gain = np.ones((128, 512, self.max_cells, 3), np.float32)
-
-            if "SlopesPC" in cons_data:
+            if "SlopesCS" in cons_data:
+                # No Chance to have both SlopesPC and SlopesCS used at the same time?
+                # What to do for SlopesCS
+                # Switch between SlopesCS and SlopesPC based on rel_gain_mode.
+                # Change pc_bools to gain_bools. or rel_gain_bools??
+                pass
+            
+            elif "SlopesPC" in cons_data:
                 slopesPC = cons_data["SlopesPC"].astype(np.float32, copy=False)
 
                 # This will handle some historical data in a different format
@@ -1423,55 +1431,14 @@ class AgipdCorrections:
                     slopesPC = np.moveaxis(slopesPC, 0, 3)
                     slopesPC = np.moveaxis(slopesPC, 0, 2)
 
-                # high gain slope from pulse capacitor data
-                pc_high_m = slopesPC[..., :self.max_cells, 0]
-                pc_high_l = slopesPC[..., :self.max_cells, 1]
-                # medium gain slope from pulse capacitor data
-                pc_med_m = slopesPC[..., :self.max_cells, 3]
-                pc_med_l = slopesPC[..., :self.max_cells, 4]
-
-                # calculate median for slopes
-                pc_high_med = np.nanmedian(pc_high_m, axis=(0, 1))
-                pc_med_med = np.nanmedian(pc_med_m, axis=(0, 1))
-
-                if variant.get("SlopesPC", 0) == 0:
-                    # calculate median for intercepts:
-                    pc_high_l_med = np.nanmedian(pc_high_l, axis=(0, 1))
-                    pc_med_l_med = np.nanmedian(pc_med_l, axis=(0, 1))
-
-                    # sanitize PC data with CCV variant = 0.
-                    # Sanitization is already done for constants
-                    # with CCV variant = 1
-                    # In the following loop,
-                    # replace `nan`s across memory cells with
-                    # the median value calculated previously.
-                    # Then, values outside of the valid range (0.8 and 1.2)
-                    # are fixed to the median value.
-                    # This is applied for high and medium gain stages
-                    for i in range(self.max_cells):
-                        pc_high_m[
-                            np.isnan(pc_high_m[..., i]), i] = pc_high_med[i]
-                        pc_med_m[
-                            np.isnan(pc_med_m[..., i]), i] = pc_med_med[i]
-
-                        pc_high_l[
-                            np.isnan(pc_high_l[..., i]), i] = pc_high_l_med[i]
-                        pc_med_l[
-                            np.isnan(pc_med_l[..., i]), i] = pc_med_l_med[i]
-
-                        pc_high_m[
-                            (pc_high_m[..., i] < 0.8 * pc_high_med[i]) |
-                            (pc_high_m[..., i] > 1.2 * pc_high_med[i]), i] = pc_high_med[i]  # noqa
-                        pc_med_m[
-                            (pc_med_m[..., i] < 0.8 * pc_med_med[i]) |
-                            (pc_med_m[..., i] > 1.2 * pc_med_med[i]), i] = pc_med_med[i]  # noqa
-
-                        pc_high_l[
-                            (pc_high_l[..., i] < 0.8 * pc_high_l_med[i]) |
-                            (pc_high_l[..., i] > 1.2 * pc_high_l_med[i]), i] = pc_high_l_med[i]  # noqa
-                        pc_med_l[
-                            (pc_med_l[..., i] < 0.8 * pc_med_l_med[i]) |
-                            (pc_med_l[..., i] > 1.2 * pc_med_l_med[i]), i] = pc_med_l_med[i]  # noqa
+                (
+                    pc_high_m, pc_med_m, pc_high_l,
+                    pc_med_l, pc_high_med, pc_med_med
+                ) = get_gain_pc_slopes(
+                    slopes_pc=slopesPC,
+                    mem_cells=self.max_cells,
+                    variant=variant
+                )
 
                 # ration between HG and MG per pixel per mem cell used
                 # for rel gain calculation
diff --git a/src/cal_tools/agipdutils.py b/src/cal_tools/agipdutils.py
index a7ee52de1..e34f5846e 100644
--- a/src/cal_tools/agipdutils.py
+++ b/src/cal_tools/agipdutils.py
@@ -650,3 +650,79 @@ def cast_array_inplace(inp, dtype):
     outp.shape = orig_shape
 
     return outp
+
+def get_gain_pc_slopes(
+    slopes_pc: np.ndarray,
+    mem_cells: int,
+    variant: int = 0,
+    ) -> (
+        np.ndarray, np.ndarray, np.ndarray,
+        np.ndarray, np.ndarray, np.ndarray
+    ):
+    """Calculate gain PC slopes and the high and medium gain medians.
+
+    Args:
+        slopes_pc (np.ndarray): SlopesPC gain constant array.
+        mem_cells (int): Number of operating memory cells.
+        variant (int): Calibration Constant Version variant.
+
+    Returns:
+        np.ndarray: _description_
+        np.ndarray: _description_
+        np.ndarray: _description_
+        np.ndarray: _description_
+        np.ndarray: _description_
+        np.ndarray: _description_
+    """
+    # high gain slope from pulse capacitor data
+    pc_high_m = slopes_pc[..., :mem_cells, 0]
+    pc_high_l = slopes_pc[..., :mem_cells, 1]
+    # medium gain slope from pulse capacitor data
+    pc_med_m = slopes_pc[..., :mem_cells, 3]
+    pc_med_l = slopes_pc[..., :mem_cells, 4]
+
+    # calculate median for slopes
+    pc_high_med = np.nanmedian(pc_high_m, axis=(0, 1))
+    pc_med_med = np.nanmedian(pc_med_m, axis=(0, 1))
+
+    if variant.get("SlopesPC", 0) == 0:
+        # calculate median for intercepts:
+        pc_high_l_med = np.nanmedian(pc_high_l, axis=(0, 1))
+        pc_med_l_med = np.nanmedian(pc_med_l, axis=(0, 1))
+
+        # sanitize PC data with CCV variant = 0.
+        # Sanitization is already done for constants
+        # with CCV variant = 1
+        # In the following loop,
+        # replace `nan`s across memory cells with
+        # the median value calculated previously.
+        # Then, values outside of the valid range (0.8 and 1.2)
+        # are fixed to the median value.
+        # This is applied for high and medium gain stages
+        for i in range(mem_cells):
+
+            pc_high_m[
+                np.isnan(pc_high_m[..., i]), i] = pc_high_med[i]
+            pc_med_m[
+                np.isnan(pc_med_m[..., i]), i] = pc_med_med[i]
+
+            pc_high_l[
+                np.isnan(pc_high_l[..., i]), i] = pc_high_l_med[i]
+            pc_med_l[
+                np.isnan(pc_med_l[..., i]), i] = pc_med_l_med[i]
+
+            pc_high_m[
+                (pc_high_m[..., i] < 0.8 * pc_high_med[i]) |
+                (pc_high_m[..., i] > 1.2 * pc_high_med[i]), i] = pc_high_med[i]  # noqa
+            pc_med_m[
+                (pc_med_m[..., i] < 0.8 * pc_med_med[i]) |
+                (pc_med_m[..., i] > 1.2 * pc_med_med[i]), i] = pc_med_med[i]  # noqa
+
+            pc_high_l[
+                (pc_high_l[..., i] < 0.8 * pc_high_l_med[i]) |
+                (pc_high_l[..., i] > 1.2 * pc_high_l_med[i]), i] = pc_high_l_med[i]  # noqa
+            pc_med_l[
+                (pc_med_l[..., i] < 0.8 * pc_med_l_med[i]) |
+                (pc_med_l[..., i] > 1.2 * pc_med_l_med[i]), i] = pc_med_l_med[i]  # noqa
+
+    return pc_high_m, pc_med_m, pc_high_l, pc_med_l, pc_high_med, pc_med_med
\ No newline at end of file
diff --git a/src/cal_tools/calcat_interface.py b/src/cal_tools/calcat_interface.py
index 164ce0ee0..d5b26c0d3 100644
--- a/src/cal_tools/calcat_interface.py
+++ b/src/cal_tools/calcat_interface.py
@@ -1014,6 +1014,8 @@ class AGIPD_CalibrationData(SplitConditionCalibrationData):
         "BadPixelsDark",
         "BadPixelsPC",
         "SlopesPC",
+        "SlopesCS",
+        "BadPixelsCS",
     }
     illuminated_calibrations = {
         "BadPixelsFF",
@@ -1155,7 +1157,7 @@ class LPD_CalibrationData(SplitConditionCalibrationData):
 
 
 class DSSC_CalibrationData(CalibrationData):
-    """Calibration data for the DSSC detetor."""
+    """Calibration data for the DSSC detector."""
 
     calibrations = {
         "Offset",
-- 
GitLab