A temporar workaround for processing dark runs for HIREX detector using the correct datasources

fc4aeced · Karim Ahmed · c4fef9ce · fc4aeced · fc4aeced
Commit fc4aeced authored 1 year ago by Karim Ahmed
--- a/notebooks/Gotthard2/Characterize_Darks_Gotthard2_NBC.ipynb
+++ b/notebooks/Gotthard2/Characterize_Darks_Gotthard2_NBC.ipynb
@@ -129,7 +129,7 @@
    "run_dc = RunDirectory(in_folder / f\"r{run_high:04d}\")\n",
    "file_loc = f\"proposal:{run_dc.run_metadata()['proposalNumber']} runs:{run_high} {run_med} {run_low}\"  # noqa\n",
    "receiver_names = [f\"*{receiver_template.format(x)}*\" for x in receiver_affixes]\n",
-    "data_sources = list(run_dc.select(receiver_names).all_sources)"
+    "data_sources = sorted(list(run_dc.select(receiver_names).all_sources))"
   ]
  },
  {
@@ -313,7 +313,14 @@
    "    np.uint16\n",
    ")\n",
    "empty_lut = np.stack(1280 * [np.stack([empty_lut] * 2)], axis=0)\n",
-    "for mod, src in zip(karabo_da, data_sources):\n",
+    "for mod in karabo_da:\n",
+    "\n",
+    "    # A temporar workaround for using the correct data_source\n",
+    "    if gh2_detector == \"25um\":\n",
+    "        # Bad assumption for CALCAT DA names e.g. GH200/1 and GH200/2\n",
+    "        src = data_sources[int(mod[-1])-1]\n",
+    "    else:\n",
+    "        src = data_sources[0]\n",
    "\n",
    "    # Retrieve LUT constant\n",
    "    lut, time = get_constant_from_db_and_time(\n",

 %% Cell type:markdown id:49b6577f-96a5-4dd2-bdd9-da661b2c4619 tags:
 # Gotthard2 Dark Image Characterization #
 Author: European XFEL Detector Group, Version: 1.0
 The following is a processing for the dark constants (`Offset`, `Noise`, and `BadPixelsDark`) maps using dark images taken with Gotthard2 detector (GH2 50um or 25um).
 All constants are evaluated per strip, per pulse, and per memory cell. The maps are calculated for each gain stage that is acquired in 3 separate runs.
 The three maps are of shape (stripes, cells, gains): (1280, 2, 3). They can be injected to the database (`db_output`) and/or stored locally (`local_output`).
 %% Cell type:code id:818e24e8 tags:
 ``` python
 in_folder = "/gpfs/exfel/exp/FXE/202231/p900298/raw"  # the folder to read data from, required
 out_folder =  "/gpfs/exfel/data/scratch/ahmedk/test/gotthard2/darks"  # the folder to output to, required
 metadata_folder = ''  # Directory containing calibration_metadata.yml when run by xfel-calibrate
 run_high = 7  # run number for G0 dark run, required
 run_med = 8  # run number for G1 dark run, required
 run_low = 9  # run number for G2 dark run, required
 # Parameters used to access raw data.
 karabo_id = "FXE_XAD_G2XES"  # karabo prefix of Gotthard-II devices
 karabo_da = [""]  # data aggregators
 receiver_template = "RECEIVER{}"  # receiver template used to read INSTRUMENT keys.
 receiver_affixes = [""]  # The affix to format into the receiver template to be able to load the correct receiver name from the data.
 control_template = "CONTROL"  # control template used to read CONTROL keys.
 ctrl_source_template = '{}/DET/{}'  # template for control source name (filled with karabo_id_control)
 karabo_id_control = ""  # Control karabo ID. Set to empty string to use the karabo-id
 # Parameters for the calibration database.
 cal_db_interface = "tcp://max-exfl-cal001:8020"  # calibration DB interface to use
 cal_db_timeout = 300000  # timeout on caldb requests
 creation_time = ""  # To overwrite the measured creation_time. Required Format: YYYY-MM-DD HR:MN:SC e.g. "2022-06-28 13:00:00"
 db_output = False  # Output constants to the calibration database
 local_output = True  # Output constants locally
 # Conditions used for injected calibration constants.
 bias_voltage = -1  # Detector bias voltage, set to -1 to use value in raw file.
 exposure_time = -1.  # Detector exposure time, set to -1 to use value in raw file.
 exposure_period = -1.  # Detector exposure period, set to -1 to use value in raw file.
 acquisition_rate = -1.  # Detector acquisition rate (1.1/4.5), set to -1 to use value in raw file.
 single_photon = -1  # Detector single photon mode (High/Low CDS), set to -1 to use value in raw file.
 # Parameters used during selecting raw data trains.
 min_trains = 1  # Minimum number of trains that should be available to process dark constants. Default 1.
 max_trains = 1000  # Maximum number of trains to use for processing dark constants. Set to 0 to use all available trains.
 badpixel_threshold_sigma = 5.  # bad pixels defined by values outside n times this std from median
 # Don't delete! myMDC sends this by default.
 operation_mode = ''  # Detector dark run acquiring operation mode, optional
 ```
 %% Cell type:code id:8085f9aa tags:
 ``` python
 import numpy as np
 import matplotlib.pyplot as plt
 import pasha as psh
 from IPython.display import Markdown, display
 from extra_data import RunDirectory
 from pathlib import Path
 import yaml
 from cal_tools.calcat_interface import CalCatApi
 from cal_tools.enums import BadPixels
 from cal_tools.gotthard2 import gotthard2algs, gotthard2lib
 from cal_tools.step_timing import StepTimer
 from cal_tools.restful_config import calibration_client
 from cal_tools.tools import (
    calcat_creation_time,
    get_constant_from_db_and_time,
    get_report,
    save_const_to_h5,
    send_to_db,
 )
 from iCalibrationDB import Conditions, Constants
 %matplotlib inline
 ```
 %% Cell type:code id:18fe4379 tags:
 ``` python
 run_nums = [run_high, run_med, run_low]
 in_folder = Path(in_folder)
 out_folder = Path(out_folder)
 out_folder.mkdir(parents=True, exist_ok=True)
 ctrl_src = ctrl_source_template.format(karabo_id_control, control_template)
 # Read report path to associate it later with injected constants.
 report = get_report(metadata_folder)
 # Run's creation time:
 creation_time = calcat_creation_time(in_folder, run_high, creation_time)
 print(f"Creation time: {creation_time}")
 if not karabo_id_control:
    karabo_id_control = karabo_id
 ```
 %% Cell type:code id:c176a86f tags:
 ``` python
 run_dc = RunDirectory(in_folder / f"r{run_high:04d}")
 file_loc = f"proposal:{run_dc.run_metadata()['proposalNumber']} runs:{run_high} {run_med} {run_low}"  # noqa
 receiver_names = [f"*{receiver_template.format(x)}*" for x in receiver_affixes]
-data_sources = list(run_dc.select(receiver_names).all_sources)
+data_sources = sorted(list(run_dc.select(receiver_names).all_sources))
 ```
 %% Cell type:code id:108be688 tags:
 ``` python
 step_timer = StepTimer()
 ```
 %% Cell type:code id:ff9149fc tags:
 ``` python
 # Read parameter conditions and validate the values for the three runs.
 step_timer.start()
 run_dcs_dict = dict()
 ctrl_src = ctrl_source_template.format(karabo_id_control, control_template)
 conditions = {
    "bias_voltage": set(),
    "exposure_time": set(),
    "exposure_period": set(),
    "acquisition_rate": set(),
    "single_photon": set(),
 }
 for gain, run in enumerate(run_nums):
    run_dc = RunDirectory(in_folder / f"r{run:04d}/")
    run_dcs_dict[run] = [gain, run_dc]
    # Read slow data.
    g2ctrl = gotthard2lib.Gotthard2Ctrl(run_dc=run_dc, ctrl_src=ctrl_src)
    conditions["bias_voltage"].add(
        g2ctrl.get_bias_voltage() if bias_voltage == -1 else bias_voltage
    )
    conditions["exposure_time"].add(
        g2ctrl.get_exposure_time() if exposure_time == -1 else exposure_time
    )
    conditions["exposure_period"].add(
        g2ctrl.get_exposure_period() if exposure_period == -1 else exposure_period
    )
    conditions["single_photon"].add(
        g2ctrl.get_single_photon() if single_photon == -1 else single_photon
    )
    conditions["acquisition_rate"].add(
        g2ctrl.get_acquisition_rate() if acquisition_rate == -1 else acquisition_rate
    )
 for c, v in conditions.items():
    assert len(v) == 1, f"{c} value is not the same for the three runs."
 bias_voltage = conditions["bias_voltage"].pop()
 print("Bias voltage: ", bias_voltage)
 exposure_time = conditions["exposure_time"].pop()
 print("Exposure time: ", exposure_time)
 exposure_period = conditions["exposure_period"].pop()
 print("Exposure period: ", exposure_period)
 single_photon = conditions["single_photon"].pop()
 print("Single photon: ", single_photon)
 acquisition_rate = conditions["acquisition_rate"].pop()
 print("Acquisition rate: ", acquisition_rate)
 gh2_detector = g2ctrl.get_det_type()
 print(f"Processing {gh2_detector} Gotthard2.")
 ```
 %% Cell type:code id:f64bc150-cfcd-4f98-83f9-a982fdacedd7 tags:
 ``` python
 calcat = CalCatApi(client=calibration_client())
 detector_id = calcat.detector(karabo_id)['id']
 pdus_by_da = calcat.physical_detector_units(detector_id, pdu_snapshot_at=creation_time)
 da_to_pdu = {da: p['physical_name'] for (da, p) in pdus_by_da.items()}
 if karabo_da != [""]:
    # Filter DA connected to detector in CALCAT
    karabo_da = [da for da in karabo_da if da in da_to_pdu]
    # Exclude non selected DA from processing.
    da_to_pdu = {da: da_to_pdu[da] for da in karabo_da}
 else:
    karabo_da = sorted(da_to_pdu.keys())
 print(f"Processing {karabo_da}")
 ```
 %% Cell type:code id:ac9c5dc3-bc66-4e7e-b6a1-360259be535c tags:
 ``` python
 def specify_trains_to_process(
    img_key_data: "extra_data.KeyData",
    run_num: int,
    src: str,
 ):
    """Specify total number of trains to process.
    Based on given min_trains and max_trains, if given.
    Print number of trains to process and number of empty trains.
    Raise ValueError if specified trains are less than min_trains.
    """
    # Specifies total number of trains to process.
    n_trains = img_key_data.shape[0]
    all_trains = len(img_key_data.train_ids)
    print(
        f"{src} for run {run} has {all_trains - n_trains}"
        f"trains with empty frames out of {all_trains} trains"
    )
    if n_trains < min_trains:
        raise ValueError(
            f"Less than {min_trains} trains are available in RAW data."
            " Not enough data to process darks."
        )
    if max_trains > 0:
        n_trains = min(n_trains, max_trains)
    print(f"Processing {n_trains} trains.")
    return n_trains
 ```
 %% Cell type:code id:3c59c11d tags:
 ``` python
 # set the operating condition
 condition = Conditions.Dark.Gotthard2(
    bias_voltage=bias_voltage,
    exposure_time=exposure_time,
    exposure_period=exposure_period,
    acquisition_rate=acquisition_rate,
    single_photon=single_photon,
 )
 ```
 %% Cell type:code id:e2eb2fc0-df9c-4887-9691-f81474f8c131 tags:
 ``` python
 def convert_train(wid, index, tid, d):
    """Convert a Gotthard2 train from 12bit to 10bit."""
    gotthard2algs.convert_to_10bit(
        d[src]["data.adc"], lut, data_10bit[index, ...]
    )
 ```
 %% Cell type:code id:4e8ffeae tags:
 ``` python
 # Calculate noise and offset per pixel and global average, std and median
 noise_map = dict()
 offset_map = dict()
 badpixels_map = dict()
 context = psh.ProcessContext(num_workers=25)
 empty_lut = (np.arange(2 ** 12).astype(np.float64) * 2 ** 10 / 2 ** 12).astype(
    np.uint16
 )
 empty_lut = np.stack(1280 * [np.stack([empty_lut] * 2)], axis=0)
-for mod, src in zip(karabo_da, data_sources):
+for mod in karabo_da:
+    # A temporar workaround for using the correct data_source
+    if gh2_detector == "25um":
+        # Bad assumption for CALCAT DA names e.g. GH200/1 and GH200/2
+        src = data_sources[int(mod[-1])-1]
+    else:
+        src = data_sources[0]
    # Retrieve LUT constant
    lut, time = get_constant_from_db_and_time(
        constant=Constants.Gotthard2.LUT(),
        condition=condition,
        empty_constant=empty_lut,
        karabo_id=karabo_id,
        karabo_da=mod,
        cal_db_interface=cal_db_interface,
        creation_time=creation_time,
        timeout=cal_db_timeout,
        print_once=False,
    )
    print(f"Retrieved LUT constant with creation-time {time}")
    cshape = (1280, 2, 3)
    offset_map[mod] = context.alloc(shape=cshape, dtype=np.float32)
    noise_map[mod] = context.alloc(like=offset_map[mod])
    badpixels_map[mod] = context.alloc(like=offset_map[mod], dtype=np.uint32)
    for run_num, [gain, run_dc] in run_dcs_dict.items():
        step_timer.start()
        n_trains = specify_trains_to_process(run_dc[src, "data.adc"], run_num, src)
        # Select requested number of trains to process.
        dc = run_dc.select(src, require_all=True).select_trains(
            np.s_[:n_trains]
        )
        step_timer.done_step("preparing raw data")
        step_timer.start()
        # Convert 12bit data to 10bit
        data_10bit = context.alloc(
            shape=dc[src, "data.adc"].shape, dtype=np.float32
        )
        context.map(convert_train, dc)
        step_timer.done_step("convert to 10bit")
        step_timer.start()
        # The first ~20 frames of each train are excluded.
        # The electronics needs some time to reach stable operational conditions
        # at the beginning of each acquisition cycle,
        # hence the first ~20 images have lower quality and should not be used.
        # The rough number of 20 is correct at 4.5 MHz acquisition rate,
        # 5 should be enough at 1.1 MHz. Sticking with 20 excluded frames, as there isn't
        # much of expected difference.
        # Split even and odd data to calculate the two storage cell constants.
        # Detector always operates in burst mode.
        even_data = data_10bit[:, 20::2, :]
        odd_data = data_10bit[:, 21::2, :]
        def offset_noise_cell(wid, index, d):
            offset_map[mod][:, index, gain] = np.mean(d, axis=(0, 1))
            noise_map[mod][:, index, gain] = np.std(d, axis=(0, 1))
        # Calculate Offset and Noise Maps.
        context.map(offset_noise_cell, (even_data, odd_data))
        # Split even and odd gain data.
        data_gain = dc[src, "data.gain"].ndarray()
        even_gain = data_gain[:, 20::2, :]
        odd_gain = data_gain[:, 21::2, :]
        raw_g = 3 if gain == 2 else gain
        def badpixels_cell(wid, index, g):
            """Check if there are wrong bad gain values.
            Indicate pixels with wrong gain value across all trains for each cell."""
            badpixels_map[mod][
                np.mean(g, axis=(0, 1)) != raw_g, index, gain
            ] |= BadPixels.WRONG_GAIN_VALUE.value
        context.map(badpixels_cell, (even_gain, odd_gain))
        step_timer.done_step("Processing darks")
 ```
 %% Cell type:code id:8e410ca2 tags:
 ``` python
 print(f"Total processing time {step_timer.timespan():.01f} s")
 step_timer.print_summary()
 ```
 %% Cell type:code id:3fc17e05-17ab-4ac4-9e79-c95399278bb9 tags:
 ``` python
 def print_bp_entry(bp):
    print(f"{bp.name:<30s} {bp.value:032b} -> {int(bp.value)}")
 print_bp_entry(BadPixels.NOISE_OUT_OF_THRESHOLD)
 print_bp_entry(BadPixels.OFFSET_NOISE_EVAL_ERROR)
 print_bp_entry(BadPixels.WRONG_GAIN_VALUE)
 def eval_bpidx(d):
    mdn = np.nanmedian(d, axis=(0))[None, :, :]
    std = np.nanstd(d, axis=(0))[None, :, :]
    idx = (d > badpixel_threshold_sigma * std + mdn) | (
        d < (-badpixel_threshold_sigma) * std + mdn
    )
    return idx
 ```
 %% Cell type:code id:40c34cc5-fe93-4b83-bf39-f465f37c40b4 tags:
 ``` python
 step_timer.start()
 g_name = ["G0", "G1", "G2"]
 for mod, pdu in da_to_pdu.items():
    display(Markdown(f"### Badpixels for module {mod}:"))
    badpixels_map[mod][
        ~np.isfinite(offset_map[mod])
    ] |= BadPixels.OFFSET_NOISE_EVAL_ERROR.value
    badpixels_map[mod][
        eval_bpidx(noise_map[mod])
    ] |= BadPixels.NOISE_OUT_OF_THRESHOLD.value
    badpixels_map[mod][
        ~np.isfinite(noise_map[mod])
    ] |= BadPixels.OFFSET_NOISE_EVAL_ERROR.value
    if not local_output:
        for cell in [0, 1]:
            fig, ax = plt.subplots(figsize=(10, 5))
            for g_idx in [0, 1, 2]:
                ax.plot(badpixels_map[mod][:, cell, g_idx], label=f"G{g_idx} Bad pixel map")
            ax.set_xticks(np.arange(0, 1281, 80))
            ax.set_xlabel("Stripes #")
            ax.set_xlabel("BadPixels")
            ax.set_title(f"BadPixels map - Cell {cell} - Module {mod} ({pdu})")
            ax.set_ylim([0, 5])
            ax.legend()
            plt.show()
 step_timer.done_step(f"Creating bad pixels constant.")
 ```
 %% Cell type:code id:c8777cfe tags:
 ``` python
 if not local_output:
    for cons, cname in zip([offset_map, noise_map], ["Offset", "Noise"]):
        for mod, pdu in da_to_pdu.items():
            display(Markdown(f"### {cname} for module {mod}:"))
            for cell in [0, 1]:
                fig, ax = plt.subplots(figsize=(10, 5))
                for g_idx in [0, 1, 2]:
                    ax.plot(cons[mod][:, cell, g_idx], label=f"G{g_idx} {cname} map")
                ax.set_xlabel("Stripes #")
                ax.set_xlabel(cname)
                ax.set_title(f"{cname} map - Cell {cell} - Module {mod} ({pdu})")
                ax.legend()
                plt.show()
 ```
 %% Cell type:code id:1c4eddf7-7d6e-49f4-8cbb-12d2bc496a8f tags:
 ``` python
 step_timer.start()
 for mod, db_mod in da_to_pdu.items():
    constants = {
        "Offset": offset_map[mod],
        "Noise": noise_map[mod],
        "BadPixelsDark": badpixels_map[mod],
    }
    md = None
    for key, const_data in constants.items():
        const = getattr(Constants.Gotthard2, key)()
        const.data = const_data
        if db_output:
            md = send_to_db(
                db_module=db_mod,
                karabo_id=karabo_id,
                constant=const,
                condition=condition,
                file_loc=file_loc,
                report_path=report,
                cal_db_interface=cal_db_interface,
                creation_time=creation_time,
                timeout=cal_db_timeout,
            )
        if local_output:
            md = save_const_to_h5(
                db_module=db_mod,
                karabo_id=karabo_id,
                constant=const,
                condition=condition,
                data=const.data,
                file_loc=file_loc,
                report=report,
                creation_time=creation_time,
                out_folder=out_folder,
            )
            print(f"Calibration constant {key} is stored locally at {out_folder}.\n")
 print("Constants parameter conditions are:\n")
 print(
    f"• Bias voltage: {bias_voltage}\n"
    f"• Exposure time: {exposure_time}\n"
    f"• Exposure period: {exposure_period}\n"
    f"• Acquisition rate: {acquisition_rate}\n"
    f"• Single photon: {single_photon}\n"
    f"• Creation time: {md.calibration_constant_version.begin_at if md is not None else creation_time}\n"
 )
 step_timer.done_step("Injecting constants.")
 ```
 %% Cell type:code id:98ca9486 tags:
 ``` python
 # TODO: store old constants for comparison.
 for mod, pdu in da_to_pdu.items():
    mod_file = mod.replace("/", "-")
    with open(f"{metadata_folder or out_folder}/module_metadata_{mod_file}.yml", "w") as fd:
        yaml.safe_dump(
            {
                "module": mod,
                "pdu": pdu,
            },
            fd,
        )
 ```

--- a/notebooks/Gotthard2/Correction_Gotthard2_NBC.ipynb
+++ b/notebooks/Gotthard2/Correction_Gotthard2_NBC.ipynb
@@ -242,7 +242,7 @@
    "\n",
    "# Create the correction receiver name.\n",
    "receiver_names = [f\"*{receiver_template.format(x)}*\" for x in receiver_affixes]\n",
-    "data_sources = list(run_dc.select(receiver_names).all_sources)\n",
+    "data_sources = sorted(list(run_dc.select(receiver_names).all_sources))\n",
    "\n",
    "if not corr_receiver:\n",
    "    # This part assumes this data_source structure: '{karabo_id}/DET/{receiver_name}:{output_channel}'\n",

 %% Cell type:markdown id:bed7bd15-21d9-4735-82c1-c27c1a5e3346 tags:
 # Gotthard2 Offline Correction
 Author: European XFEL Detector Group, Version: 1.0
 Offline Correction for Gotthard2 Detector.
 This notebook is able to correct 25um and 50um GH2 detectors using the same correction steps:
 - Convert 12bit raw data into 10bit, offset subtraction, then multiply with gain constant.
 | Correction | constants   | boolean to enable/disable   |
 |------------|-------------|-----------------------------|
 |   12bit to 10bit  | `LUTGotthard2` |  |
 |   Offset  | `OffsetGotthard2`|`offset_correction`|
 |   Relative gain  | `RelativeGainGotthard2` + `BadPixelsFFGotthard2` |`gain_correction`|
 Beside the corrected data, a mask is stored using the badpixels constant of the same parameter conditions and time.
 - `BadPixelsDarkGotthard2`
 - `BadPixelsFFGotthard2`, if relative gain correction is requested.
 The correction is done per sequence file. If all selected sequence files have no images to correct the notebook will fail.
 The same result would be reached in case the needed dark calibration constants were not retrieved for all modules and `offset_correction` is True.
 In case one of the gain constants were not retrieved `gain_correction` is switched to False and gain correction is disabled.
 The `data` datasets stored in the RECEIVER source along with the corrected image (`adc`) and `mask` are:
  - `gain`
  - `bunchId`
  - `memoryCell`
  - `frameNumber`
  - `timestamp`
  - `trainId`
 %% Cell type:code id:570322ed-f611-4fd1-b2ec-c12c13d55843 tags:
 ``` python
 in_folder = "/gpfs/exfel/exp/DETLAB/202330/p900326/raw"  # the folder to read data from, required
 out_folder = "/gpfs/exfel/data/scratch/ahmedk/test/gotthard2"  # the folder to output to, required
 metadata_folder = ""  # Directory containing calibration_metadata.yml when run by xfel-calibrate
 run = 20  # run to process, required
 sequences = [-1]  # sequences to correct, set to [-1] for all, range allowed
 sequences_per_node = 1  # number of sequence files per node if notebook executed through xfel-calibrate, set to 0 to not run SLURM parallel
 # Parameters used to access raw data.
 karabo_id = "DETLAB_25UM_GH2"  # karabo prefix of Gotthard-II devices
 karabo_da = [""]  # data aggregators
 receiver_template = "RECEIVER{}"  # receiver template used to read INSTRUMENT keys.
 receiver_affixes = [""]  # The affix to format into the receiver template to be able to load the correct receiver name from the data.
 control_template = "CONTROL"  # control template used to read CONTROL keys.
 ctrl_source_template = "{}/DET/{}"  # template for control source name (filled with karabo_id_control)
 karabo_id_control = ""  # Control karabo ID. Set to empty string to use the karabo-id
 corr_source_template = "{}/CORR/{}:daqOutput"  # Correction data source template. filled with karabo_id and correction receiver
 corr_receiver = ""  # The receiver name of the corrected data. Leave empty for using the same receiver name for the 50um GH2 or the first(Master) receiver for the 25um GH2.
 # Parameters for calibration database.
 cal_db_interface = "tcp://max-exfl-cal001:8016#8025"  # the database interface to use.
 cal_db_timeout = 180000  # timeout on caldb requests.
 creation_time = ""  # To overwrite the measured creation_time. Required Format: YYYY-MM-DD HR:MN:SC e.g. "2022-06-28 13:00:00"
 # Parameters affecting corrected data.
 constants_file = ""  # Use constants in given constant file path. /gpfs/exfel/data/scratch/ahmedk/dont_remove/gotthard2/constants/calibration_constants_GH2.h5
 offset_correction = True  # apply offset correction. This can be disabled to only apply LUT or apply LUT and gain correction for non-linear differential results.
 gain_correction = True  # apply gain correction.
 chunks_data = 1  # HDF chunk size for pixel data in number of frames.
 # Parameter conditions.
 bias_voltage = -1  # Detector bias voltage, set to -1 to use value in raw file.
 exposure_time = -1.  # Detector exposure time, set to -1 to use value in raw file.
 exposure_period = -1.  # Detector exposure period, set to -1 to use value in raw file.
 acquisition_rate = -1.  # Detector acquisition rate (1.1/4.5), set to -1 to use value in raw file.
 single_photon = -1  # Detector single photon mode (High/Low CDS), set to -1 to use value in raw file.
 reverse_second_module = -1  # Reverse 25um GH2 second module before interleaving. set to -1 to use value in raw file to reverse based on `CTRL/reverseSlaveReadOutMode`'s value.
 # Parameters for plotting
 skip_plots = False  # exit after writing corrected files
 pulse_idx_preview = 3  # pulse index to preview. The following even/odd pulse index is used for preview. # TODO: update to pulseId preview.
 def balance_sequences(in_folder, run, sequences, sequences_per_node, karabo_da):
    from xfel_calibrate.calibrate import balance_sequences as bs
    return bs(in_folder, run, sequences, sequences_per_node, karabo_da)
 ```
 %% Cell type:code id:6e9730d8-3908-41d7-abe2-d78e046d5de2 tags:
 ``` python
 import warnings
 from logging import warning
 import h5py
 import pasha as psh
 import numpy as np
 import matplotlib.pyplot as plt
 from IPython.display import Markdown, display
 from extra_data import RunDirectory, H5File
 from pathlib import Path
 import cal_tools.restful_config as rest_cfg
 from cal_tools.calcat_interface import CalCatError, GOTTHARD2_CalibrationData
 from cal_tools.files import DataFile
 from cal_tools.gotthard2 import gotthard2algs, gotthard2lib
 from cal_tools.step_timing import StepTimer
 from cal_tools.tools import (
    calcat_creation_time,
    write_constants_fragment,
    map_seq_files,
 )
 from XFELDetAna.plotting.heatmap import heatmapPlot
 warnings.filterwarnings('ignore')
 %matplotlib inline
 ```
 %% Cell type:code id:d7c02c48-4429-42ea-a42e-de45366d7fa3 tags:
 ``` python
 in_folder = Path(in_folder)
 run_folder = in_folder / f"r{run:04d}"
 out_folder = Path(out_folder)
 out_folder.mkdir(parents=True, exist_ok=True)
 if not karabo_id_control:
    karabo_id_control = karabo_id
 ctrl_src = ctrl_source_template.format(karabo_id_control, control_template)
 # Run's creation time:
 creation_time = calcat_creation_time(in_folder, run, creation_time)
 print(f"Creation time: {creation_time}")
 ```
 %% Cell type:code id:f9a8d1eb-ce6a-4ed0-abf4-4a6029734672 tags:
 ``` python
 step_timer = StepTimer()
 ```
 %% Cell type:code id:892172d8 tags:
 ``` python
 run_dc = RunDirectory(run_folder)
 # Read slow data
 g2ctrl = gotthard2lib.Gotthard2Ctrl(run_dc=run_dc, ctrl_src=ctrl_src)
 if bias_voltage == -1:
    bias_voltage = g2ctrl.get_bias_voltage()
 if exposure_time == -1:
    exposure_time = g2ctrl.get_exposure_time()
 if exposure_period == -1:
    exposure_period = g2ctrl.get_exposure_period()
 if acquisition_rate == -1:
    acquisition_rate = g2ctrl.get_acquisition_rate()
 if single_photon == -1:
    single_photon = g2ctrl.get_single_photon()
 gh2_detector = g2ctrl.get_det_type()
 if reverse_second_module == -1 and gh2_detector == "25um":
    reverse_second_module = not g2ctrl.second_module_reversed()
    print(
        "Second module is not reversed. "
        "Reversing second module before interleaving.")
 print("Bias Voltage:", bias_voltage)
 print("Exposure Time:", exposure_time)
 print("Exposure Period:", exposure_period)
 print("Acquisition Rate:", acquisition_rate)
 print("Single Photon:", single_photon)
 print(f"Processing {gh2_detector} Gotthard2.")
 ```
 %% Cell type:code id:21a8953a-8c76-475e-8f4f-b201cc25c159 tags:
 ``` python
 # GH2 calibration data object.
 g2_cal = GOTTHARD2_CalibrationData(
    detector_name=karabo_id,
    sensor_bias_voltage=bias_voltage,
    exposure_time=exposure_time,
    exposure_period=exposure_period,
    acquisition_rate=acquisition_rate,
    single_photon=single_photon,
    event_at=creation_time,
    client=rest_cfg.calibration_client(),
 )
 da_to_pdu = None
 # Keep as long as it is essential to correct
 # RAW data (FXE p003225) before the data mapping was added to CALCAT.
 try:  # in case local constants are used with old RAW data. This can be removed in the future.
    da_to_pdu = g2_cal.mod_to_pdu
 except CalCatError as e:
    print(e)
    db_modules = [None] * len(karabo_da)
 if da_to_pdu:
    if karabo_da == [""]:
        karabo_da = sorted(da_to_pdu.keys())
    else:
        # Exclude non selected DA from processing.
        karabo_da = [da for da in karabo_da if da in da_to_pdu]
    db_modules = [da_to_pdu[da] for da in karabo_da]
 print(f"Process modules: {db_modules} for run {run}")
 # Create the correction receiver name.
 receiver_names = [f"*{receiver_template.format(x)}*" for x in receiver_affixes]
-data_sources = list(run_dc.select(receiver_names).all_sources)
+data_sources = sorted(list(run_dc.select(receiver_names).all_sources))
 if not corr_receiver:
    # This part assumes this data_source structure: '{karabo_id}/DET/{receiver_name}:{output_channel}'
    if gh2_detector == "25um":  # For 25um use virtual karabo_das for CALCAT data mapping.
        corr_receiver = data_sources[0].split("/")[-1].split(":")[0][:-2]
    else:
        corr_receiver = data_sources[0].split("/")[-1].split(":")[0]
    print(f"Using {corr_receiver} as a receiver name for the corrected data.")
 ```
 %% Cell type:code id:2551b923 tags:
 ``` python
 # Check the available trains to correct.
 total_trains = len(RunDirectory(run_folder).select(data_sources, require_all=True).train_ids)
 if total_trains:
    print(f"Correcting {total_trains}.")
 else:
    raise ValueError(f"No trains to correct for run {run}.")
 ```
 %% Cell type:markdown id:8c852392-bb19-4c40-b2ce-3b787538a92d tags:
 ### Retrieving calibration constants
 %% Cell type:code id:5717d722 tags:
 ``` python
 # Used for old FXE (p003225) runs before adding Gotthard2 to CALCAT
 const_data = dict()
 if constants_file:
    for mod in karabo_da:
        const_data[mod] = dict()
        # load constants temporarily using defined local paths.
        with h5py.File(constants_file, "r") as cfile:
            const_data[mod]["LUTGotthard2"] = cfile["LUT"][()]
            const_data[mod]["OffsetGotthard2"] = cfile["offset_map"][()].astype(np.float32)
            const_data[mod]["RelativeGainGotthard2"] = cfile["gain_map"][()].astype(np.float32)
            const_data[mod]["Mask"] = cfile["bpix_ff"][()].astype(np.uint32)
 else:
    constant_names = ["LUTGotthard2", "OffsetGotthard2", "BadPixelsDarkGotthard2"]
    if gain_correction:
        constant_names += ["RelativeGainGotthard2", "BadPixelsFFGotthard2"]
    g2_metadata = g2_cal.metadata(calibrations=constant_names)
    # Display retrieved calibration constants timestamps
    g2_cal.display_markdown_retrieved_constants(metadata=g2_metadata)
    # Validate the constants availability and raise/warn correspondingly.
    for mod, calibrations in g2_metadata.items():
        dark_constants = {"LUTGotthard2"}
        if offset_correction:
            dark_constants |= {"OffsetGotthard2", "BadPixelsDarkGotthard2"}
        missing_dark_constants = dark_constants - set(calibrations)
        if missing_dark_constants:
            karabo_da.remove(mod)
            warning(f"Dark constants {missing_dark_constants} are not available to correct {mod}.")  # noqa
        missing_gain_constants = {
            "BadPixelsFFGotthard2", "RelativeGainGotthard2"} - set(calibrations)
        if gain_correction and missing_gain_constants:
            warning(f"Gain constants {missing_gain_constants} are not retrieved for mod {mod}.")
 if not karabo_da:
    raise ValueError("Dark constants are not available for all modules.")
 ```
 %% Cell type:code id:ac1cdec5 tags:
 ``` python
 # Record constant details in YAML metadata.
 write_constants_fragment(
    out_folder=(metadata_folder or out_folder),
    det_metadata=g2_metadata,
    caldb_root=g2_cal.caldb_root)
 # Load constants data for all constants.
 const_data = g2_cal.ndarray_map(metadata=g2_metadata)
 # Prepare constant arrays.
 if not constants_file:
    for mod in karabo_da:
        # Create the mask array.
        bpix = const_data[mod].get("BadPixelsDarkGotthard2")
        if bpix is None:
            bpix = np.zeros((1280, 2, 3), dtype=np.uint32)
        if const_data[mod].get("BadPixelsFFGotthard2") is not None:
            bpix |= const_data[mod]["BadPixelsFFGotthard2"]
        const_data[mod]["Mask"] = bpix
        # Prepare empty arrays for missing constants.
        if const_data[mod].get("OffsetGotthard2") is None:
            const_data[mod]["OffsetGotthard2"] = np.zeros(
                (1280, 2, 3), dtype=np.float32)
        if const_data[mod].get("RelativeGainGotthard2") is None:
            const_data[mod]["RelativeGainGotthard2"] = np.ones(
                (1280, 2, 3), dtype=np.float32)
        const_data[mod]["RelativeGainGotthard2"] = const_data[mod]["RelativeGainGotthard2"].astype(  # noqa
            np.float32, copy=False)  # Old gain constants are not float32.
 ```
 %% Cell type:code id:2c7dd0bb tags:
 ``` python
 file_da = list({kda.split('/')[0] for kda in karabo_da})
 mapped_files, total_files = map_seq_files(
    run_folder,
    file_da,
    sequences,
 )
 # This notebook doesn't account for processing more
 # than one file data aggregator.
 seq_files = mapped_files[file_da[0]]
 if not len(seq_files):
    raise IndexError(
        "No sequence files available to correct for the selected sequences and karabo_da.")
 print(f"Processing a total of {total_files} sequence files")
 ```
 %% Cell type:code id:23fcf7f4-351a-4df7-8829-d8497d94fecc tags:
 ``` python
 context = psh.ProcessContext(num_workers=23)
 ```
 %% Cell type:code id:daecd662-26d2-4cb8-aa70-383a579cf9f9 tags:
 ``` python
 def correct_train(wid, index, d):
    g = gain[index]
    gotthard2algs.convert_to_10bit(d, const_data[mod]["LUTGotthard2"], data_corr[index, ...])
    gotthard2algs.correct_train(
        data_corr[index, ...],
        mask[index, ...],
        g,
        const_data[mod]["OffsetGotthard2"].astype(np.float32),  # PSI map is in f8
        const_data[mod]["RelativeGainGotthard2"],
        const_data[mod]["Mask"],
        apply_offset=offset_correction,
        apply_gain=gain_correction,
    )
 ```
 %% Cell type:code id:f88c1aa6-a735-4b72-adce-b30162f5daea tags:
 ``` python
 corr_data_source = corr_source_template.format(karabo_id, corr_receiver)
 for raw_file in seq_files:
    out_file = out_folder / raw_file.name.replace("RAW", "CORR")
    # Select module INSTRUMENT sources and deselect empty trains.
    dc = H5File(raw_file).select(data_sources, require_all=True)
    n_trains = len(dc.train_ids)
    if n_trains == 0:
        warning(f"Skipping correction. No trains to correct for this sequence file: {raw_file}.")
        continue
    else:
        print(f"Correcting {n_trains} for {raw_file}.")
    # Initialize GH2 data and gain arrays to store in corrected files.
    if gh2_detector == "25um":
        data_stored = np.zeros((dc[data_sources[0], "data.adc"].shape[:2] + (1280 * 2,)), dtype=np.float32)
        gain_stored = np.zeros((dc[data_sources[0], "data.adc"].shape[:2] + (1280 * 2,)), dtype=np.uint8)
    else:
        data_stored = None
        gain_stored = None
    for i, (src, mod) in enumerate(zip(data_sources, karabo_da)):
        step_timer.start()
        print(f"Correcting {src} for {raw_file}")
        data = dc[src, "data.adc"].ndarray()
        gain = dc[src, "data.gain"].ndarray()
        step_timer.done_step("Preparing raw data")
        dshape = data.shape
        step_timer.start()
        # Allocate shared arrays.
        data_corr = context.alloc(shape=dshape, dtype=np.float32)
        mask = context.alloc(shape=dshape, dtype=np.uint32)
        context.map(correct_train, data)
        step_timer.done_step(f"Correcting one receiver in one sequence file")
        step_timer.start()
        # Provided PSI gain map has 0 values. Set inf values to nan.
        # TODO: This can maybe be removed after creating XFEL gain maps.?
        data_corr[np.isinf(data_corr)] = np.nan
        # Create CORR files and add corrected data sections.
        image_counts = dc[src, "data.adc"].data_counts(labelled=False)
        if reverse_second_module and i == 1:
            data_corr = np.flip(data_corr, axis=-1)
            gain = np.flip(gain, axis=-1)
        if gh2_detector == "25um":
            data_stored[..., i::2] = data_corr.copy()
            gain_stored[..., i::2] = gain.copy()
        else:  # "50um"
            data_stored = data_corr
            gain_stored = gain
    with DataFile(out_file, "w") as ofile:
        # Create INDEX datasets.
        ofile.create_index(dc.train_ids, from_file=dc.files[0])
        ofile.create_metadata(
            like=dc,
            sequence=dc.run_metadata()["sequenceNumber"],
            instrument_channels=(f"{corr_data_source}/data",)
        )
        # Create Instrument section to later add corrected datasets.
        outp_source = ofile.create_instrument_source(corr_data_source)
        # Create count/first datasets at INDEX source.
        outp_source.create_index(data=image_counts)
        # Store uncorrected trainId in the corrected file.
        outp_source.create_key(
                f"data.trainId", data=dc.train_ids,
                chunks=min(50, len(dc.train_ids))
            )
        # Create datasets with the available corrected data
        for field_name, field_data in {
            "adc": data_stored,
            "gain": gain_stored,
        }.items():
            outp_source.create_key(
                f"data.{field_name}", data=field_data,
                chunks=((chunks_data,) + data_corr.shape[1:])
        )
        # For GH2 25um, the data of the second receiver is
        # stored in the corrected file.
        for field in ["bunchId", "memoryCell", "frameNumber", "timestamp"]:
                outp_source.create_key(
                    f"data.{field}", data=dc[src, f"data.{field}"].ndarray(),
                    chunks=(chunks_data, data_corr.shape[1])
            )
        outp_source.create_compressed_key(f"data.mask", data=mask)
        step_timer.done_step("Storing data")
 ```
 %% Cell type:code id:94b8e4d2-9f8c-4c23-a509-39238dd8435c tags:
 ``` python
 print(f"Total processing time {step_timer.timespan():.01f} s")
 step_timer.print_summary()
 ```
 %% Cell type:code id:0ccc7f7e-2a3f-4ac0-b854-7d505410d2fd tags:
 ``` python
 if skip_plots:
    print("Skipping plots")
    import sys
    sys.exit(0)
 ```
 %% Cell type:code id:ff203f77-3811-46f3-bf7d-226d2dcab13f tags:
 ``` python
 mod_dcs = {}
 first_seq_raw = seq_files[0]
 first_seq_corr = out_folder / first_seq_raw.name.replace("RAW", "CORR")
 mod_dcs[corr_data_source] = {}
 with H5File(first_seq_corr) as out_dc:
    tid, mod_dcs[corr_data_source]["train_corr_data"] = next(
        out_dc[corr_data_source, "data.adc"].trains()
    )
 if gh2_detector == "25um":
    mod_dcs[corr_data_source]["train_raw_data"] = np.zeros((data_corr.shape[1], 1280 * 2), dtype=np.float32)
    mod_dcs[corr_data_source]["train_raw_gain"] = np.zeros((data_corr.shape[1], 1280 * 2), dtype=np.uint8)
 for i, src in enumerate(data_sources):
    with H5File(first_seq_raw) as in_dc:
        train_dict = in_dc.train_from_id(tid)[1][src]
        if gh2_detector == "25um":
            if reverse_second_module and i == 1:
                data = np.flip(train_dict["data.adc"], axis=-1)
                gain = np.flip(train_dict["data.gain"], axis=-1)
            else:
                data = train_dict["data.adc"]
                gain = train_dict["data.gain"]
            mod_dcs[corr_data_source]["train_raw_data"][..., i::2] = data
            mod_dcs[corr_data_source]["train_raw_gain"][..., i::2] = gain
        else:
            mod_dcs[corr_data_source]["train_raw_data"] = train_dict["data.adc"]
            mod_dcs[corr_data_source]["train_raw_gain"] = train_dict["data.gain"]
 ```
 %% Cell type:code id:1b379438-eb1d-42b2-ac83-eb8cf88c46db tags:
 ``` python
 display(Markdown("### Mean RAW and CORRECTED across pulses for one train:"))
 display(Markdown(f"Train: {tid}"))
 if gh2_detector == "50um":
    title = f"{{}} data for {karabo_da} ({db_modules})"
 else:
    title = f"Interleaved {{}} data for {karabo_da} ({db_modules})"
 step_timer.start()
 fig, ax = plt.subplots(figsize=(15, 15))
 raw_data = mod_dcs[corr_data_source]["train_raw_data"]
 im = ax.plot(np.mean(raw_data, axis=0))
 ax.set_title(title.format("RAW"), fontsize=20)
 ax.set_xlabel("Strip #", size=20)
 ax.set_ylabel("12-bit ADC output", size=20)
 plt.xticks(fontsize=20)
 plt.yticks(fontsize=20)
 pass
 fig, ax = plt.subplots(figsize=(15, 15))
 corr_data = mod_dcs[corr_data_source]["train_corr_data"]
 im = ax.plot(np.mean(corr_data, axis=0))
 ax.set_title(title.format("CORRECTED"), fontsize=20)
 ax.set_xlabel("Strip #", size=20)
 ax.set_ylabel("10-bit KeV. output", size=20)
 plt.xticks(fontsize=20)
 plt.yticks(fontsize=20)
 pass
 step_timer.done_step("Plotting mean data")
 ```
 %% Cell type:code id:58a6a276 tags:
 ``` python
 display(Markdown(f"### RAW and CORRECTED strips across pulses for train {tid}"))
 step_timer.start()
 for plt_data, dname in zip(
    ["train_raw_data", "train_corr_data"], ["RAW", "CORRECTED"]
 ):
    fig, ax = plt.subplots(figsize=(15, 15))
    plt.rcParams.update({"font.size": 20})
    heatmapPlot(
        mod_dcs[corr_data_source][plt_data],
        y_label="Pulses",
        x_label="Strips",
        title=title.format(dname),
        use_axis=ax,
        cb_pad=0.8,
    )
    pass
 step_timer.done_step("Plotting RAW and CORRECTED data for one train")
 ```
 %% Cell type:code id:cd8f5e08-fcee-4bff-ba63-6452b3d892a2 tags:
 ``` python
 # Validate given "pulse_idx_preview"
 if pulse_idx_preview + 1 > data.shape[1]:
    print(
        f"WARNING: selected pulse_idx_preview {pulse_idx_preview} is not available in data."
        " Previewing 1st pulse."
    )
    pulse_idx_preview = 1
 if data.shape[1] == 1:
    odd_pulse = 1
    even_pulse = None
 else:
    odd_pulse = pulse_idx_preview if pulse_idx_preview % 2 else pulse_idx_preview + 1
    even_pulse = (
        pulse_idx_preview if not (pulse_idx_preview % 2) else pulse_idx_preview + 1
    )
 if pulse_idx_preview + 1 > data.shape[1]:
    pulse_idx_preview = 1
    if data.shape[1] > 1:
        pulse_idx_preview = 2
 ```
 %% Cell type:code id:e5f0d4d8-e32c-4f2c-8469-4ebbfd3f644c tags:
 ``` python
 display(Markdown("### RAW and CORRECTED even/odd pulses for one train:"))
 display(Markdown(f"Train: {tid}"))
 fig, ax = plt.subplots(figsize=(15, 15))
 raw_data = mod_dcs[corr_data_source]["train_raw_data"]
 corr_data = mod_dcs[corr_data_source]["train_corr_data"]
 ax.plot(raw_data[odd_pulse], label=f"Odd Pulse {odd_pulse}")
 if even_pulse:
    ax.plot(raw_data[even_pulse], label=f"Even Pulse {even_pulse}")
 ax.set_title(title.format("RAW"), fontsize=20)
 ax.set_xlabel("Strip #", size=20)
 ax.set_ylabel("12-bit ADC RAW", size=20)
 plt.xticks(fontsize=20)
 plt.yticks(fontsize=20)
 ax.legend()
 pass
 fig, ax = plt.subplots(figsize=(15, 15))
 ax.plot(corr_data[odd_pulse], label=f"Odd Pulse {odd_pulse}")
 if even_pulse:
    ax.plot(corr_data[even_pulse], label=f"Even Pulse {even_pulse}")
 ax.set_title(title.format("CORRECTED"), fontsize=20)
 ax.set_xlabel("Strip #", size=20)
 ax.set_ylabel("10-bit KeV CORRECTED", size=20)
 plt.xticks(fontsize=20)
 plt.yticks(fontsize=20)
 ax.legend()
 pass
 step_timer.done_step("Plotting RAW and CORRECTED odd/even pulses.")
 ```