[Gotthard2][Dark] Calculate noise bpixs in multiple iterations

Description

@mramilli mentioned that reiterating on calculating NOISE_OUT_OF_THRESHOLD badpixel type for gotthard2 showed a difference in the number of badpixels identified.

In the below figure are his test results. Each plot is the difference in bpix numbers for each gain for one of the production constants.

Below are my test results showing that sometimes you can get a different number of badpixels with more iterations than 2.

This MR exposes an iteration parameter to control it if needed (default 5).

Now GH2 darks will iterate to calculate the more accurate number of badpixels. each iteration costs =< 1ms, so no obvious performance drawback on dark processing.

How Has This Been Tested?

Relevant Documents (optional)

Types of changes

Checklist:

Reviewers

@mramilli @schmidtp @kluyvert

notebooks/Gotthard2/Characterize_Darks_Gotthard2_NBC.ipynb

 %% 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.
 input_source_template = "{karabo_id}/DET/RECEIVER{input_source_affixes}:daqOutput"  # data source template used to read INSTRUMENT keys.
 ctrl_source_template = "{karabo_id_control}/DET/{control_template}"  # template for control source name (filled with karabo_id_control)
 karabo_id = "FXE_XAD_G2XES"  # karabo prefix of Gotthard-II devices
 karabo_da = [""]  # data aggregators
 input_source_affixes = [""]  # The affix to format into the input source template to be able to load the correct module's data.
 control_template = "CONTROL"  # control template used to read CONTROL keys.
 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 for calculating calibration constants
 bpix_noise_nitrs = 5  # number of maximum iterations to actively try to reach correct number of badpixels of NOISE_OUT_OF_THRESHOLD type. Leave -1 to keep iterating until finding maximum number of badpixels.
 # 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,
     raw_data_location_string,
     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)
 
 if not karabo_id_control:
     karabo_id_control = karabo_id
 
 ctrl_src = ctrl_source_template.format(
     karabo_id_control=karabo_id_control,
     control_template=control_template,
 )
 run_nums = gotthard2lib.sort_dark_runs_by_gain(
     raw_folder=in_folder,
     runs=run_nums,
     ctrl_src=ctrl_src,
 )
 
 # 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_nums[0], creation_time)
 print(f"Creation time: {creation_time}")
 ```
 
 %% Cell type:code id:c176a86f tags:
 
 ``` python
 run_dc = RunDirectory(in_folder / f"r{run_nums[0]:04d}")
 proposal = list(filter(None, str(in_folder).strip('/').split('/')))[-2]
 file_loc = raw_data_location_string(proposal, run_nums)
 ```
 
 %% 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()
 
 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)
 # Module index is stored in the da_to_pdu dict to be used
 # later in selecting the correct data source.
 # TODO: This can be improved later by getting the module index and even data source
 # when we start storing in the CALCAT.
 da_to_pdu = dict()
 for i, (da, p) in enumerate(sorted(pdus_by_da.items())):
     da_to_pdu[da] = p['physical_name']
 
 calcat_das = list(da_to_pdu.keys())
 
 if karabo_da != [""]:
     # Filter DA connected to detector in CALCAT
     karabo_da = [da for da in karabo_da if da in calcat_das]
 else:
     karabo_da = calcat_das
 
 print(f"Processing {karabo_da}")
 
 da_to_source = gotthard2lib.map_da_to_source(
     run_dc,
     calcat_das,
     input_source_template,
     karabo_id,
     input_source_affixes,
 )
 ```
 
 %% 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 in karabo_da:
     # For 50um idx always 0. For 25um pick the right data source
     # from list of 2 sources.
     src = da_to_source[mod]
 
     # 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, :
             ] |= 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 in karabo_da:
     pdu = da_to_pdu[mod]
     display(Markdown(f"### Badpixels for module {mod}:"))
 
     badpixels_map[mod][
         ~np.isfinite(offset_map[mod])
     ] |= BadPixels.OFFSET_NOISE_EVAL_ERROR.value
     badpixels_map[mod][
         ~np.isfinite(noise_map[mod])
     ] |= BadPixels.OFFSET_NOISE_EVAL_ERROR.value
 
     nmap = noise_map[mod].copy()
     last_nbpix = 0
-    # Calculate NOISE_OUT_OF_THRESHOLD bpix in multiple iterations.
-    nitr = -1
-    while nitr != bpix_noise_nitrs:
+    # Iteratively identify noisy pixels (NOISE_OUT_OF_THRESHOLD)
+    # Due to the low number of strips (1280), outliers can significantly impact
+    # the badpixel calculation. We use multiple iterations to gradually refine
+    # the sigma threshold, as a single fixed threshold may not be suitable for
+    # different Gotthard2 detectors.
+    for itr in range(bpix_noise_nitrs):
         badpixels_map[mod][
             eval_bpidx(nmap)
         ] |= BadPixels.NOISE_OUT_OF_THRESHOLD.value
         nmap[badpixels_map[mod] > 0] = np.nan
         recent_nbpix = np.count_nonzero(badpixels_map[mod])
 
         if last_nbpix == recent_nbpix:
-            # skip iterations if no more difference.
-            nitr = bpix_noise_nitrs
-        else:
-            nitr += 1
-        last_nbpix = recent_nbpix
+            # stop iterating if no more badpixels difference
+            break
 
-    if local_output:
+    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 in karabo_da:
             pdu = da_to_pdu[mod]
             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 in karabo_da:
     db_mod = da_to_pdu[mod]
     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 in karabo_da:
     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": da_to_pdu[mod],
             },
             fd,
         )
 ```