bin/slurm_calibrate.sh

@@ -43,9 +43,14 @@ then
     sleep 15
 fi
 
-
 echo "Running notebook"
-${python_path} -m princess ${nb_path} --save
+if [ "$caltype" == "CORRECT" ]
+then
+  # calparrot stores and repeats calcat queries
+  ${python_path} -m calparrot -- ${python_path} -m princess ${nb_path} --save
+else
+  ${python_path} -m princess ${nb_path} --save
+fi
 
 # stop the cluster if requested
 if [ "${ipcluster_profile}" != "NO_CLUSTER" ]

notebooks/Gotthard2/Correction_Gotthard2_NBC.ipynb

 %% Cell type:markdown id:bed7bd15-21d9-4735-82c1-c27c1a5e3346 tags:
 
 # Gotthard2 Offline Correction #
 
 Author: European XFEL Detector Group, Version: 1.0
 
 Offline Calibration for the Gothard2 Detector
 
 %% Cell type:code id:570322ed-f611-4fd1-b2ec-c12c13d55843 tags:
 
 ``` python
 in_folder = "/gpfs/exfel/exp/FXE/202221/p003225/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 = 50  # 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 = "FXE_XAD_G2XES"  # karabo prefix of Gotthard-II devices
 karabo_da = ["GH201"]  # data aggregators
 receiver_template = "RECEIVER"  # receiver template used to read INSTRUMENT keys.
 control_template = "CONTROL"  # control template used to read CONTROL keys.
 instrument_source_template = "{}/DET/{}:daqOutput"  # template for source name (filled with karabo_id & receiver_id). e.g. 'SPB_IRDA_JF4M/DET/JNGFR01:daqOutput'
 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 calibration database.
 use_dir_creation_date = True  # use the creation data of the input dir for database queries.
 cal_db_interface = "tcp://max-exfl016:8016#8025"  # the database interface to use.
 cal_db_timeout = 180000  # timeout on caldb requests.
 overwrite_creation_time = ""  # To overwrite the measured creation_time. Required Format: YYYY-MM-DD HR:MN:SC.00 e.g. "2022-06-28 13:00: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.
 
 # 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.
 
 # 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 datetime
 import warnings
 from functools import partial
 
 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
 
 from cal_tools import h5_copy_except
 from cal_tools.gotthard2 import gotthard2algs, gotthard2lib
 from cal_tools.step_timing import StepTimer
 from cal_tools.tools import (
     get_constant_from_db_and_time,
     get_dir_creation_date,
     get_pdu_from_db,
     CalibrationMetadata,
 )
 from iCalibrationDB import Conditions, Constants
 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)
 
 metadata = CalibrationMetadata(metadata_folder or out_folder)
 # NOTE: this notebook will not overwrite calibration metadata file
 const_yaml = metadata.get("retrieved-constants", {})
 
 if not karabo_id_control:
     karabo_id_control = karabo_id
 
 instrument_src = instrument_source_template.format(karabo_id, receiver_template)
 ctrl_src = ctrl_source_template.format(karabo_id_control, control_template)
 
 print(f"Process modules: {karabo_da} for run {run}")
 
 creation_time = None
 if overwrite_creation_time:
     creation_time = datetime.datetime.strptime(
         overwrite_creation_time, "%Y-%m-%d %H:%M:%S.%f"
     )
 elif use_dir_creation_date:
     creation_time = get_dir_creation_date(in_folder, run)
     print(f"Using {creation_time} as creation time")
 ```
 
 %% Cell type:code id:b5eb816e-b5f2-44ce-9907-0273d82341b6 tags:
 
 ``` python
 # Select only sequence files to process for the selected detector.
 if sequences == [-1]:
     possible_patterns = list(f"*{mod}*.h5" for mod in karabo_da)
 else:
     possible_patterns = list(
         f"*{mod}-S{s:05d}.h5" for mod in karabo_da for s in sequences
     )
 
 run_folder = Path(in_folder / f"r{run:04d}")
 seq_files = [
     f for f in run_folder.glob("*.h5") if any(f.match(p) for p in possible_patterns)
 ]
 
 seq_files = sorted(seq_files)
 
 if not seq_files:
     raise IndexError("No sequence files available for the selected sequences.")
 
 print(f"Processing a total of {len(seq_files)} sequence files")
 ```
 
 %% Cell type:code id:f9a8d1eb-ce6a-4ed0-abf4-4a6029734672 tags:
 
 ``` python
 step_timer = StepTimer()
 ```
 
 %% Cell type:code id:892172d8 tags:
 
 ``` python
 # Read slow data
 run_dc = RunDirectory(run_folder)
 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()
 
 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)
 ```
 
 %% 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]["LUT"] = cfile["LUT"][()]
             const_data[mod]["Offset"] = cfile["offset_map"][()].astype(np.float32)
             const_data[mod]["RelativeGain"] = cfile["gain_map"][()].astype(np.float32)
             const_data[mod]["Mask"] = cfile["bpix_ff"][()].astype(np.uint32)
 ```
 
 %% Cell type:code id:1cdbe818 tags:
 
 ``` python
 # Conditions iCalibrationDB object.
 condition = Conditions.Dark.Gotthard2(
     bias_voltage=bias_voltage,
     exposure_time=exposure_time,
     exposure_period=exposure_period,
     single_photon=single_photon,
     acquisition_rate=acquisition_rate,
 )
 
 # TODO: Maybe this condition and previous cell can be removed later after the initial phase.
 if not constants_file:
     # Prepare a dictionary of empty constants to loop on
     # it's keys and initiate non-retrieved constants.
     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)
     empty_constants = {
         "LUT": empty_lut,
         "Offset": np.zeros((1280, 2, 3), dtype=np.float32),
         "BadPixelsDark": np.zeros((1280, 2, 3), dtype=np.uint32),
         "RelativeGain": np.ones((1280, 2, 3), dtype=np.float32),
         "BadPixelsFF": np.zeros((1280, 2, 3), dtype=np.uint32),
     }
 
     for mod in karabo_da:
         const_data[mod] = dict()
         # Only used for printing timestamps within the loop.
         when = dict()
         # Check YAML file for constant metadata of file path and creation-time
         if const_yaml:
             for cname, mdata in const_yaml[mod]["constants"].items():
                 const_data[mod][cname] = dict()
                 when[cname] = mdata["creation-time"]
                 if when[cname]:
                     with h5py.File(mdata["file-path"], "r") as cf:
                         const_data[mod][cname] = np.copy(
                             cf[f"{mdata['dataset-name']}/data"]
                         )
                 else:
                     const_data[mod][cname] = empty_constants[cname]
         else:  # Retrieve constants from CALCAT. Missing YAML file or running notebook interactively.
             for cname, cempty in empty_constants.items():
                 const_data[mod][cname] = dict()
                 const_data[mod][cname], when[cname] = get_constant_from_db_and_time(
                     karabo_id=karabo_id,
                     karabo_da=mod,
                     cal_db_interface=cal_db_interface,
                     creation_time=creation_time,
                     timeout=cal_db_timeout,
                     print_once=False,
                     condition=condition,
                     constant=getattr(Constants.Gotthard2, cname)(),
                     empty_constant=cempty,
                 )
         bpix = const_data[mod]["BadPixelsDark"]
         bpix |= const_data[mod]["BadPixelsFF"]
         const_data[mod]["Mask"] = bpix
 
         # Print timestamps for the retrieved constants.
         print(f"Constants for module {mod}:")
         for cname, ctime in when.items():
             print(f"  {cname} injected at {ctime}")
         del when
 ```
 
 %% 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]["LUT"], data_corr[index, ...])
     gotthard2algs.correct_train(
         data_corr[index, ...],
         mask[index, ...],
         g,
         const_data[mod]["Offset"],
-        const_data[mod]["RelativeGain"],
+        const_data[mod]["RelativeGain"].astype(np.float32, copy=False),
         const_data[mod]["Mask"],
         apply_offset=offset_correction,
         apply_gain=gain_correction,
     )
 ```
 
 %% Cell type:code id:f88c1aa6-a735-4b72-adce-b30162f5daea tags:
 
 ``` python
 for mod in karabo_da:
     # This is used in case receiver template consists of
     # karabo data aggregator index. e.g. detector at DETLAB
     instr_mod_src = instrument_src.format(mod[-2:])
     data_path = "INSTRUMENT/" + instr_mod_src + "/data"
     for raw_file in seq_files:
         step_timer.start()
 
         dc = H5File(raw_file)
         out_file = out_folder / raw_file.name.replace("RAW", "CORR")
 
         # Select module INSTRUMENT source and deselect empty trains.
         dc = dc.select(instr_mod_src, require_all=True)
         data = dc[instr_mod_src, "data.adc"].ndarray()
         gain = dc[instr_mod_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("Correcting 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 sources.
         # Exclude raw data images (data/adc)
         with h5py.File(out_file, "w") as ofile:
             # Copy RAW non-calibrated sources.
             with h5py.File(raw_file, "r") as sfile:
                 h5_copy_except.h5_copy_except_paths(sfile, ofile, [f"{data_path}/adc"])
             # Create datasets with the available corrected data
             ddset = ofile.create_dataset(
                 f"{data_path}/adc",
                 data=data_corr,
                 chunks=((1,) + dshape[1:]),  # 1 chunk == 1 image
                 dtype=np.float32,
             )
             # Create datasets with the available corrected data
             ddset = ofile.create_dataset(
                 f"{data_path}/mask",
                 data=mask,
                 chunks=((1,) + dshape[1:]),  # 1 chunk == 1 image
                 dtype=np.uint32,
                 compression="gzip",
                 compression_opts=1,
                 shuffle=True,
             )
         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")
 for mod in karabo_da:
     mod_dcs[mod] = {}
     with H5File(first_seq_corr) as out_dc:
         tid, mod_dcs[mod]["train_corr_data"] = next(
             out_dc[instr_mod_src, "data.adc"].trains()
         )
     with H5File(first_seq_raw) as in_dc:
         train_dict = in_dc.train_from_id(tid)[1][instr_mod_src]
         mod_dcs[mod]["train_raw_data"] = train_dict["data.adc"]
         mod_dcs[mod]["train_raw_gain"] = train_dict["data.gain"]
 ```
 
 %% Cell type:code id:494b453a tags:
 
 ``` python
 # Keep as long as it is essential to correct
 # RAW data (FXE p003225) before the data mapping was added to CALCAT.
 try:
     db_modules = get_pdu_from_db(
         karabo_id=karabo_id,
         karabo_da=karabo_da,
         constant=Constants.jungfrau.Offset(),
         condition=condition,
         cal_db_interface=cal_db_interface,
         snapshot_at=creation_time,
     )
 except RuntimeError:
     print(
         "No Physical detector units found for this"
         " detector mapping at the RAW data creation time."
     )
     db_modules = [None] * len(karabo_da)
 ```
 
 %% 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}"))
 
 step_timer.start()
 for mod, pdu in zip(karabo_da, db_modules):
 
     fig, ax = plt.subplots(figsize=(20, 10))
     raw_data = mod_dcs[mod]["train_raw_data"]
     im = ax.plot(np.mean(raw_data, axis=0))
     ax.set_title(f"RAW module {mod} ({pdu})")
     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=(20, 10))
     corr_data = mod_dcs[mod]["train_corr_data"]
     im = ax.plot(np.mean(corr_data, axis=0))
     ax.set_title(f"CORRECTED module {mod} ({pdu})")
     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 mod, pdu in zip(karabo_da, db_modules):
     for plt_data, dname in zip(
         ["train_raw_data", "train_corr_data"], ["RAW", "CORRECTED"]
     ):
         fig, ax = plt.subplots(figsize=(15, 20))
         plt.rcParams.update({"font.size": 20})
 
         heatmapPlot(
             mod_dcs[mod][plt_data],
             y_label="Pulses",
             x_label="Strips",
             title=f"{dname} module {mod} ({pdu})",
             use_axis=ax,
         )
         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}"))
 for mod, pdu in zip(karabo_da, db_modules):
     fig, ax = plt.subplots(figsize=(20, 20))
     raw_data = mod_dcs[mod]["train_raw_data"]
     corr_data = mod_dcs[mod]["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(f"RAW module {mod} ({pdu})")
     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=(20, 20))
     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(f"CORRECTED module {mod} ({pdu})")
     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.")
 ```

notebooks/Jungfrau/Jungfrau_Gain_Correct_and_Verify_NBC.ipynb

 %% Cell type:markdown id: tags:
 
 # Jungfrau Offline Correction #
 
 Author: European XFEL Detector Group, Version: 2.0
 
 Offline Calibration for the Jungfrau Detector
 
 %% Cell type:code id: tags:
 
 ``` python
 in_folder = "/gpfs/exfel/exp/SPB/202130/p900204/raw"  # the folder to read data from, required
 out_folder =  "/gpfs/exfel/data/scratch/ahmedk/test/remove"  # the folder to output to, required
 run = 91  # run to process, required
 metadata_folder = ""  # Directory containing calibration_metadata.yml when run by xfel-calibrate
 sequences = [-1]  # sequences to correct, set to [-1] for all, range allowed
 sequences_per_node = 1  # number of sequence files per cluster node if run as slurm job, set to 0 to not run SLURM parallel
 
 # Parameters used to access raw data.
 karabo_id = "SPB_IRDA_JF4M"  # karabo prefix of Jungfrau devices
 karabo_da = ['JNGFR01', 'JNGFR02', 'JNGFR03', 'JNGFR04', 'JNGFR05', 'JNGFR06', 'JNGFR07', 'JNGFR08']  # data aggregators
 receiver_template = "JNGFR{:02d}"  # Detector receiver template for accessing raw data files. e.g. "JNGFR{:02d}"
 instrument_source_template = '{}/DET/{}:daqOutput'  # template for source name (filled with karabo_id & receiver_id). e.g. 'SPB_IRDA_JF4M/DET/JNGFR01:daqOutput'
 ctrl_source_template = '{}/DET/CONTROL'  # template for control source name (filled with karabo_id_control)
 karabo_id_control = ""  # if control is on a different ID, set to empty string if it is the same a karabo-id
 
 # Parameters for calibration database.
 use_dir_creation_date = True  # use the creation data of the input dir for database queries
 cal_db_interface = "tcp://max-exfl016:8017#8025" # the database interface to use
 cal_db_timeout = 180000  # timeout on caldb requests
 
 # Parameters affecting corrected data.
 relative_gain = True  # do relative gain correction.
 strixel_sensor = False  # reordering for strixel detector layout.
 strixel_double_norm = 2.0  # normalization to use for double-size pixels, only applied for strixel sensors.
 limit_trains = 0  # ONLY FOR TESTING. process only first N trains, Use 0 to process all.
 chunks_ids = 32  # HDF chunk size for memoryCell and frameNumber.
 chunks_data = 1  # HDF chunk size for pixel data in number of frames.
 
 # Parameters for retrieving calibration constants
 manual_slow_data = False  # if true, use manually entered bias_voltage, integration_time, gain_setting, and gain_mode values
 integration_time = 4.96  # integration time in us, will be overwritten by value in file
 gain_setting = 0  # 0 for dynamic gain, 1 for dynamic HG0, will be overwritten by value in file
 gain_mode = 0  # 0 for runs with dynamic gain setting, 1 for fixgain. It will be overwritten by value in file, if manual_slow_data is set to True.
 mem_cells = -1  # Set mem_cells to -1 to automatically use the value stored in RAW data.
 bias_voltage = 180  # will be overwritten by value in file
 
 # Parameters for plotting
 skip_plots = False  # exit after writing corrected files
 plot_trains = 500  # Number of trains to plot for RAW and CORRECTED plots. Set to -1 to automatically plot all trains.
 cell_id_preview = 15  # cell Id used for preview in single-shot plots
 
 # Parameters for ROI selection and reduction
 roi_definitions = [-1]  # List with groups of 6 values defining ROIs, e.g. [3, 120, 180, 200, 550, -2] for module 3 (JNGFR03), slice 120:180, 200:550, average along axis -2 (slow scan, or -1 for fast scan)
 
 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: tags:
 
 ``` python
 import multiprocessing
 import sys
 import warnings
 from functools import partial
 from logging import warning
 from pathlib import Path
 
 import h5py
 import matplotlib
 import matplotlib.pyplot as plt
 import numpy as np
 import pasha as psh
 import tabulate
 from IPython.display import Latex, Markdown, display
 from extra_data import H5File, RunDirectory, by_id, components
 from extra_geom import JUNGFRAUGeometry
 from matplotlib.colors import LogNorm
 
 from cal_tools import h5_copy_except
 from cal_tools.jungfraulib import JungfrauCtrl
 from cal_tools.enums import BadPixels
 from cal_tools.files import DataFile
 from cal_tools.step_timing import StepTimer
 from cal_tools.tools import (
     get_constant_from_db_and_time,
     get_dir_creation_date,
     get_pdu_from_db,
     map_seq_files,
     CalibrationMetadata,
 )
 from iCalibrationDB import Conditions, Constants
 
 warnings.filterwarnings('ignore')
 
 matplotlib.use('agg')
 %matplotlib inline
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 in_folder = Path(in_folder)
 out_folder = Path(out_folder)
 run_folder = in_folder / f'r{run:04d}'
 run_dc = RunDirectory(run_folder)
 instrument_src = instrument_source_template.format(karabo_id, receiver_template)
 
 out_folder.mkdir(parents=True, exist_ok=True)
 
 print(f"Run is: {run}")
 print(f"Instrument H5File source: {instrument_src}")
 print(f"Process modules: {karabo_da}")
 
 creation_time = None
 if use_dir_creation_date:
     creation_time = get_dir_creation_date(in_folder, run)
     print(f"Using {creation_time} as creation time")
 
 if karabo_id_control == "":
     karabo_id_control = karabo_id
 
 if any(axis_no not in {-2, -1, 2, 3} for axis_no in roi_definitions[5::6]):
     print("ROI averaging must be on axis 2/3 (or equivalently -2/-1). "
           f"Axis numbers given: {roi_definitions[5::6]}")
     sys.exit(1)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # Read available sequence files to correct.
 mapped_files, num_seq_files = map_seq_files(
     run_folder, karabo_da, sequences)
 
 if not len(mapped_files):
     raise IndexError(
         "No sequence files available to correct for the selected sequences and karabo_da.")
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 print(f"Processing a total of {num_seq_files} sequence files")
 table = []
 fi = 0
 for kda, sfiles in mapped_files.items():
     for k, f in enumerate(sfiles):
         if k == 0:
             table.append((fi, kda, k, f))
         else:
             table.append((fi, "", k,  f))
         fi += 1
 md = display(Latex(tabulate.tabulate(
     table, tablefmt='latex',
     headers=["#", "module", "# module", "file"])))
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 ctrl_src = ctrl_source_template.format(karabo_id_control)
 ctrl_data = JungfrauCtrl(run_dc, ctrl_src)
 
 if mem_cells < 0:
     memory_cells, sc_start = ctrl_data.get_memory_cells()
 
     mem_cells_name = "single cell" if memory_cells == 1 else "burst"
     print(f"Run is in {mem_cells_name} mode.\nStorage cell start: {sc_start:02d}")
 else:
     memory_cells = mem_cells
     mem_cells_name = "single cell" if memory_cells == 1 else "burst"
     print(f"Run is in manually set to {mem_cells_name} mode. With {memory_cells} memory cells")
 
 if not manual_slow_data:
     integration_time = ctrl_data.get_integration_time()
     bias_voltage = ctrl_data.get_bias_voltage()
     gain_setting = ctrl_data.get_gain_setting()
     gain_mode = ctrl_data.get_gain_mode()
 
 print(f"Integration time is {integration_time} us")
 print(f"Gain setting is {gain_setting} (run settings: {ctrl_data.run_settings})")
 print(f"Gain mode is {gain_mode} ({ctrl_data.run_mode})")
 print(f"Bias voltage is {bias_voltage} V")
 print(f"Number of memory cells are {memory_cells}")
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 if strixel_sensor:
     from cal_tools.jfstrixel import STRIXEL_SHAPE as strixel_frame_shape, double_pixel_indices, to_strixel
     Ydouble, Xdouble = double_pixel_indices()
     print('Strixel sensor transformation enabled')
 ```
 
 %% Cell type:markdown id: tags:
 
 ### Retrieving calibration constants ###
 
 %% Cell type:code id: tags:
 
 ``` python
 condition = Conditions.Dark.jungfrau(
     memory_cells=memory_cells,
     bias_voltage=bias_voltage,
     integration_time=integration_time,
     gain_setting=gain_setting,
     gain_mode=gain_mode,
 )
 
 empty_constants = {
     "Offset": np.zeros((512, 1024, memory_cells, 3), dtype=np.float32),
     "BadPixelsDark": np.zeros((512, 1024, memory_cells, 3), dtype=np.uint32),
     "RelativeGain": None,
     "BadPixelsFF": None,
 }
 metadata = CalibrationMetadata(metadata_folder or out_folder)
 # NOTE: this notebook will not overwrite calibration metadata file
 const_yaml = metadata.get("retrieved-constants", {})
 
 def get_constants_for_module(karabo_da: str):
     """ Get calibration constants for given module of Jungfrau
 
     :return:
         offset_map (offset map),
         mask (mask of bad pixels),
         gain_map (map of relative gain factors),
         db_module (name of DB module),
         when (dictionary: constant - creation time)
     """
 
     when = dict()
     const_data = dict()
 
     if const_yaml:
         for cname, mdata in const_yaml[karabo_da]["constants"].items():
             const_data[cname] = dict()
             when[cname] = mdata["creation-time"]
             if when[cname]:
                 with h5py.File(mdata["file-path"], "r") as cf:
                     const_data[cname] = np.copy(
                         cf[f"{mdata['dataset-name']}/data"])
             else:
                 const_data[cname] = empty_constants[cname]
     else:
         retrieval_function = partial(
             get_constant_from_db_and_time,
             karabo_id=karabo_id,
             karabo_da=karabo_da,
             cal_db_interface=cal_db_interface,
             creation_time=creation_time,
             timeout=cal_db_timeout,
             print_once=False,
         )
 
         for cname, cempty in empty_constants.items():
             const_data[cname], when[cname] = retrieval_function(
                 condition=condition,
                 constant=getattr(Constants.jungfrau, cname)(),
                 empty_constant=cempty,
             )
 
     offset_map = const_data["Offset"]
     mask = const_data["BadPixelsDark"]
     gain_map = const_data["RelativeGain"]
     mask_ff = const_data["BadPixelsFF"]
 
     # Combine masks
     if mask_ff is not None:
         mask |= np.moveaxis(mask_ff, 0, 1)
 
     if memory_cells > 1:
         # move from x, y, cell, gain to cell, x, y, gain
         offset_map = np.moveaxis(offset_map, [0, 1], [1, 2])
         mask = np.moveaxis(mask, [0, 1], [1, 2])
     else:
         offset_map = np.squeeze(offset_map)
         mask = np.squeeze(mask)
 
     # masking double size pixels
     mask[..., [255, 256], :, :] |= BadPixels.NON_STANDARD_SIZE
     mask[..., [255, 256, 511, 512, 767, 768], :] |= BadPixels.NON_STANDARD_SIZE
 
     if gain_map is not None:
         if memory_cells > 1:
             gain_map = np.moveaxis(gain_map, [0, 2], [2, 0])
             # add extra empty cell constant
             b = np.ones(((1,)+gain_map.shape[1:]))
             gain_map = np.concatenate((gain_map, b), axis=0)
         else:
             gain_map = np.moveaxis(np.squeeze(gain_map), 1, 0)
 
     return offset_map, mask, gain_map, karabo_da, when
 
 with multiprocessing.Pool() as pool:
     r = pool.map(get_constants_for_module, karabo_da)
 
 # Print timestamps for the retrieved constants.
 constants = {}
 for offset_map, mask, gain_map, k_da, when in r:
     print(f'Constants for module {k_da}:')
     for const in when:
         print(f'  {const} injected at {when[const]}')
 
     if gain_map is None:
         print("No gain map found")
         relative_gain = False
 
     constants[k_da] = (offset_map, mask, gain_map)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # Correct a chunk of images for offset and gain
 def correct_train(wid, index, d):
     d = d.astype(np.float32)  # [cells, x, y]
     g = gain[index]
 
     # Copy gain over first to keep it at the original 3 for low gain.
     if strixel_sensor:
         to_strixel(g, out=gain_corr[index, ...])
     else:
         gain_corr[index, ...] = g
 
     # Jungfrau gains 0[00], 1[01], 3[11]
     # Change low gain to 2 for indexing purposes.
     g[g==3] = 2
 
     # Select memory cells
     if memory_cells > 1:
         """
         Even though it is correct to assume that memory cells pattern
         can be the same across all trains (for one correction run
         taken with one acquisition), it is preferred to not assume
         this to account for exceptions that can happen.
         """
         m = memcells[index].copy()
         # 255 is a cell value pointing to no cell image data (image of 0 pixels).
         # Corresponding image will be corrected with constant of cell 0. To avoid values of 0.
         # This line is depending on not storing the modified memory cells in the corrected data.
         m[m==255] = 0
 
         offset_map_cell = offset_map[m, ...]  # [16 + empty cell, x, y]
         mask_cell = mask[m, ...]
     else:
         offset_map_cell = offset_map
         mask_cell = mask
 
     # Offset correction
     offset = np.choose(g, np.moveaxis(offset_map_cell, -1, 0))
 
     d -= offset
 
     # Gain correction
     if relative_gain:
         if memory_cells > 1:
             gain_map_cell = gain_map[m, ...]
         else:
             gain_map_cell = gain_map
         cal = np.choose(g, np.moveaxis(gain_map_cell, -1, 0))
         d /= cal
 
     msk = np.choose(g, np.moveaxis(mask_cell, -1, 0))
 
     if strixel_sensor:
         to_strixel(d, out=data_corr[index, ...])
         data_corr[index, :, Ydouble, Xdouble] /= strixel_double_norm
         to_strixel(msk, out=mask_corr[index, ...])
     else:
         data_corr[index, ...] = d
         mask_corr[index, ...] = msk
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 step_timer = StepTimer()
 
 n_cpus = multiprocessing.cpu_count()
 context = psh.context.ProcessContext(num_workers=n_cpus)
 print(f"Using {n_cpus} workers for correction.")
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 def save_reduced_rois(ofile, data_corr, mask_corr, karabo_da):
     """If ROIs are defined for this karabo_da, reduce them and save to the output file"""
     rois_defined = 0
     module_no = int(karabo_da[-2:])
     params_source = f'{karabo_id}/ROIPROC/{karabo_da}'
     rois_source = f'{params_source}:output'
     if roi_definitions != [-1]:
         # Create Instrument and Control sections to later add datasets.
         outp_source = ofile.create_instrument_source(rois_source)
         ctrl_source = ofile.create_control_source(params_source)
     for i in range(len(roi_definitions) // 6):
         roi_module, a1, a2, b1, b2, mean_axis = roi_definitions[i*6 : (i+1)*6]
         if roi_module == module_no:
             rois_defined += 1
             # Apply the mask and average remaining pixels to 1D
             roi_data = data_corr[..., a1:a2, b1:b2].mean(
                 axis=mean_axis, where=(mask_corr[..., a1:a2, b1:b2] == 0)
             )
 
             # Add roi corrected datasets
             outp_source.create_key(f'data.roi{rois_defined}.data', data=roi_data)
 
             # Add roi run control datasets.
             ctrl_source.create_run_key(f'roi{rois_defined}.region', np.array([[a1, a2, b1, b2]]))
             ctrl_source.create_run_key(f'roi{rois_defined}.reduce_axis', np.array([mean_axis]))
 
     if rois_defined:
         # Copy the index for the new source
         # Create count/first datasets at INDEX source.
         ofile.copy(f'INDEX/{karabo_id}/DET/{karabo_da}:daqOutput/data',
                    f'INDEX/{rois_source}/data')
         ntrains = ofile['INDEX/trainId'].shape[0]
         ctrl_source.create_index(ntrains)
 ```
 
 %% Cell type:markdown id: tags:
 
 ### Correcting RAW data ###
 
 %% Cell type:code id: tags:
 
 ``` python
 # Loop over modules
 empty_seq = 0
 for local_karabo_da, mapped_files_module in mapped_files.items():
     instrument_src_kda = instrument_src.format(int(local_karabo_da[-2:]))
 
     for sequence_file in mapped_files_module:
         # Save corrected data in an output file with name
         # of corresponding raw sequence file.
         ofile_name = sequence_file.name.replace("RAW", "CORR")
         out_file = out_folder / ofile_name
 
         # Load sequence file data collection, data.adc keydata,
         # the shape for data to later created arrays of the same shape,
         # and number of available trains to correct.
         seq_dc = H5File(sequence_file)
         seq_dc_adc = seq_dc[instrument_src_kda, "data.adc"]
         ishape = seq_dc_adc.shape  # input shape.
         corr_ntrains = ishape[0]  # number of available trains to correct.
         all_train_ids = seq_dc_adc.train_ids
 
         # Raise a WARNING if this sequence has no trains to correct.
         # Otherwise, print number of trains with no data.
         if corr_ntrains == 0:
             warning(f"No trains to correct for {sequence_file.name}: "
                  "Skipping the processing of this file.")
             empty_seq += 1
             continue
         elif len(all_train_ids) != corr_ntrains:
             print(f"{sequence_file.name} has {len(seq_dc_adc.train_ids) - corr_ntrains} "
                   "trains with missing data.")
 
         # For testing, limit corrected trains. i.e. Getting output faster.
         if limit_trains > 0:
             print(f"\nCorrected trains are limited to: {limit_trains} trains")
             corr_ntrains = min(corr_ntrains, limit_trains)
 
         print(f"\nNumber of corrected trains are: {corr_ntrains} for {ofile_name}")
 
         # Load constants from the constants dictionary.
         # These arrays are used by `correct_train()` function
         offset_map, mask, gain_map = constants[local_karabo_da]
 
         # Determine total output shape.
         if strixel_sensor:
             oshape = (*ishape[:-2], *strixel_frame_shape)
         else:
             oshape = ishape
 
         # Allocate shared arrays for corrected data. Used in `correct_train()`
         data_corr = context.alloc(shape=oshape, dtype=np.float32)
         gain_corr = context.alloc(shape=oshape, dtype=np.uint8)
         mask_corr = context.alloc(shape=oshape, dtype=np.uint32)
 
         step_timer.start()
         # Overwrite seq_dc after eliminating empty trains or/and applying limited images.
         seq_dc = seq_dc.select(
             instrument_src_kda, "*", require_all=True).select_trains(np.s_[:corr_ntrains])
 
         # Load raw images(adc), gain, memcells, and frame numbers.
         data = seq_dc[instrument_src_kda, "data.adc"].ndarray()
         gain = seq_dc[instrument_src_kda, "data.gain"].ndarray()
         memcells = seq_dc[instrument_src_kda, "data.memoryCell"].ndarray()
         frame_number = seq_dc[instrument_src_kda, "data.frameNumber"].ndarray()
 
         # Validate that the selected cell id to preview is available in raw data.
         if memory_cells > 1:
             # For plotting, assuming that memory cells are sorted the same for all trains.
             found_cells = memcells[0] == cell_id_preview
             if any(found_cells):
                 cell_idx_preview = np.where(found_cells)[0][0]
             else:
                 print(f"The selected cell_id_preview {cell_id_preview} is not available in burst mode. "
                       f"Previewing cell `{memcells[0]}`.")
                 cell_idx_preview = 0
         else:
             cell_idx_preview = 0
 
         # Correct data per train
         context.map(correct_train, data)
         step_timer.done_step(f"Correction time.")
 
         step_timer.start()
 
         # Create CORR files and add corrected data sections.
         image_counts = seq_dc[instrument_src_kda, "data.adc"].data_counts(labelled=False)
 
         with DataFile(out_file, 'w') as outp_file:
             # Create INDEX datasets.
             outp_file.create_index(seq_dc.train_ids, from_file=seq_dc.files[0])
 
             # Create Instrument section to later add corrected datasets.
             outp_source = outp_file.create_instrument_source(instrument_src_kda)
 
             # Create count/first datasets at INDEX source.
             outp_source.create_index(data=image_counts)
 
             # RAW memoryCell and frameNumber are not corrected. But we are storing only
             # the values for the corrected trains.
             outp_source.create_key(
                 "data.memoryCell", data=memcells,
                 chunks=(min(chunks_ids, memcells.shape[0]), 1))
             outp_source.create_key(
                 "data.frameNumber", data=frame_number,
                 chunks=(min(chunks_ids, frame_number.shape[0]), 1))
             # Add main corrected `data.adc`` dataset and store corrected data.
             outp_source.create_key(
                 "data.adc", data=data_corr,
                 chunks=(min(chunks_data, data_corr.shape[0]), *oshape[1:]))
             outp_source.create_compressed_key(
                 "data.gain", data=gain_corr)
             outp_source.create_compressed_key(
                 "data.mask", data=mask_corr)
 
             # Temporary hotfix for FXE assuming this dataset is in corrected files.
             outp_source.create_key(
                 "data.trainId", data=seq_dc.train_ids,
                 chunks=(min(50, len(seq_dc.train_ids))))
 
             save_reduced_rois(outp_file, data_corr, mask_corr, local_karabo_da)
 
             # Create METDATA datasets
             outp_file.create_metadata(like=seq_dc)
 
         step_timer.done_step(f'Saving data time.')
 if empty_seq == sum([len(i) for i in mapped_files.values()]):
     warning("No valid trains for RAW data to correct.")
     sys.exit(0)
 ```
 
 %% Cell type:markdown id: tags:
 
 ### Processing time summary ###
 
 %% Cell type:code id: tags:
 
 ``` python
 print(f"Total processing time {step_timer.timespan():.01f} s")
 step_timer.print_summary()
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 if skip_plots:
     print('Skipping plots')
     sys.exit(0)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # Positions are given in pixels
 mod_width = (256 * 4) + (2 * 3)  # inc. 2px gaps between tiles
 mod_height = (256 * 2) + 2
 if karabo_id == "SPB_IRDA_JF4M":
     # The first 4 modules are rotated 180 degrees relative to the others.
     # We pass the bottom, beam-right corner of the module regardless of its
     # orientation, requiring a subtraction from the symmetric positions we'd
     # otherwise calculate.
     x_start, y_start = 1125, 1078
     module_pos = [
         (x_start - mod_width, y_start - mod_height - (i * (mod_height + 33)))
         for i in range(4)
     ] + [
         (-x_start, -y_start + (i * (mod_height + 33))) for i in range(4)
     ]
     orientations = [(-1, -1) for _ in range(4)] + [(1, 1) for _ in range(4)]
 elif karabo_id == "FXE_XAD_JF1M":
     module_pos = ((-mod_width//2, 33),(-mod_width//2, -mod_height -33))
     orientations = [(-1,-1), (1,1)]
 else:
     module_pos = ((-mod_width//2,-mod_height//2),)
     orientations = None
 
 geom = JUNGFRAUGeometry.from_module_positions(module_pos, orientations=orientations, asic_gap=0)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 first_seq = 0 if sequences == [-1] else sequences[0]
 
 with RunDirectory(out_folder, f"*{run}*S{first_seq:05d}*") as corr_dc:
     # Reading CORR data for plotting.
     jf_corr = components.JUNGFRAU(
         corr_dc,
         detector_name=karabo_id,
     ).select_trains(np.s_[:plot_trains])
     tid, jf_corr_data = next(iter(jf_corr.trains(require_all=True)))
 
 # Shape = [modules, trains, cells, x, y]
 # TODO: Fix the case if not all modules were requested to be corrected.
 # For example if only one modules was corrected. An assertion error is expected
 # at `geom.plot_data_fast`, while plotting corrected images.
 corrected = jf_corr.get_array("data.adc")[:, :, cell_idx_preview, ...].values
 corrected_train = jf_corr_data["data.adc"][
     :, cell_idx_preview, ...
 ].values  # loose the train axis.
 
 mask = jf_corr.get_array("data.mask")[:, :, cell_idx_preview, ...].values
 mask_train = jf_corr_data["data.mask"][:, cell_idx_preview, ...].values
 
-with RunDirectory(f"{in_folder}/r{run:04d}/", f"*S{first_seq:05d}*") as raw_dc:
+with RunDirectory(f"{in_folder}/r{run:04d}/", f"*S{first_seq:05d}*", _use_voview=False) as raw_dc:
 
     # Reading RAW data for plotting.
     jf_raw = components.JUNGFRAU(raw_dc, detector_name=karabo_id).select_trains(
         np.s_[:plot_trains]
     )
 
 raw = jf_raw.get_array("data.adc")[:, :, cell_idx_preview, ...].values
 raw_train = (
     jf_raw.select_trains(by_id[[tid]])
     .get_array("data.adc")[:, 0, cell_idx_preview, ...]
     .values
 )
 
 gain = jf_raw.get_array("data.gain")[:, :, cell_idx_preview, ...].values
 gain_train_cells = (
     jf_raw.select_trains(by_id[[tid]]).get_array("data.gain")[:, :, :, ...].values
 )
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 db_modules = get_pdu_from_db(
     karabo_id=karabo_id,
     karabo_da=karabo_da,
     constant=Constants.jungfrau.Offset(),
     condition=condition,
     cal_db_interface=cal_db_interface,
     snapshot_at=creation_time,
 )
 ```
 
 %% Cell type:markdown id: tags:
 
 ### Mean RAW Preview
 
 %% Cell type:code id: tags:
 
 ``` python
 print(f"The per pixel mean of the first {raw.shape[1]} trains of the first sequence file")
 
 fig, ax = plt.subplots(figsize=(18, 10))
 raw_mean = np.mean(raw, axis=1)
 geom.plot_data_fast(
     raw_mean,
     ax=ax,
     vmin=min(0.75*np.median(raw_mean[raw_mean > 0]), 2000),
     vmax=max(1.5*np.median(raw_mean[raw_mean > 0]), 16000),
     cmap="jet",
     colorbar={'shrink': 1, 'pad': 0.01},
 )
 ax.set_title(f'{karabo_id} - Mean RAW', size=18)
 plt.show()
 ```
 
 %% Cell type:markdown id: tags:
 
 ### Mean CORRECTED Preview
 
 %% Cell type:code id: tags:
 
 ``` python
 print(f"The per pixel mean of the first {corrected.shape[1]} trains of the first sequence file")
 
 fig, ax = plt.subplots(figsize=(18, 10))
 corrected_mean = np.mean(corrected, axis=1)
 _corrected_vmin = min(0.75*np.median(corrected_mean[corrected_mean > 0]), -0.5)
 _corrected_vmax = max(2.*np.median(corrected_mean[corrected_mean > 0]), 100)
 
 mean_plot_kwargs = dict(
     vmin=_corrected_vmin, vmax=_corrected_vmax, cmap="jet"
 )
 
 if not strixel_sensor:
     geom.plot_data_fast(
         corrected_mean,
         ax=ax,
         colorbar={'shrink': 1, 'pad': 0.01},
         **mean_plot_kwargs
     )
 else:
     ax.imshow(corrected_mean.squeeze(), aspect=10, **mean_plot_kwargs)
 
 ax.set_title(f'{karabo_id} - Mean CORRECTED', size=18)
 
 plt.show()
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 fig, ax = plt.subplots(figsize=(18, 10))
 corrected_masked = corrected.copy()
 corrected_masked[mask != 0] = np.nan
 corrected_masked_mean = np.nanmean(corrected_masked, axis=1)
 del corrected_masked
 
 if not strixel_sensor:
     geom.plot_data_fast(
         corrected_masked_mean,
         ax=ax,
         colorbar={'shrink': 1, 'pad': 0.01},
         **mean_plot_kwargs
     )
 else:
     ax.imshow(corrected_mean.squeeze(), aspect=10, **mean_plot_kwargs)
 
 ax.set_title(f'{karabo_id} - Mean CORRECTED with mask', size=18)
 
 plt.show()
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 display(Markdown((f"#### A single image from train {tid}")))
 
 fig, ax = plt.subplots(figsize=(18, 10))
 
 single_plot_kwargs = dict(
     vmin=min(0.75 * np.median(corrected_train[corrected_train > 0]), -0.5),
     vmax=max(2.0 * np.median(corrected_train[corrected_train > 0]), 100),
     cmap="jet"
 )
 
 if not strixel_sensor:
     geom.plot_data_fast(
         corrected_train,
         ax=ax,
         colorbar={"shrink": 1, "pad": 0.01},
         **single_plot_kwargs
     )
 else:
     ax.imshow(corrected_train.squeeze(), aspect=10, **single_plot_kwargs)
 
 ax.set_title(f"{karabo_id} - CORRECTED train: {tid}", size=18)
 
 plt.show()
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 def do_2d_plot(data, edges, y_axis, x_axis, title):
     fig = plt.figure(figsize=(10, 10))
     ax = fig.add_subplot(111)
     extent = [
         np.min(edges[1]),
         np.max(edges[1]),
         np.min(edges[0]),
         np.max(edges[0]),
     ]
 
     im = ax.imshow(
         data[::-1, :],
         extent=extent,
         aspect="auto",
         norm=LogNorm(vmin=1, vmax=np.max(data))
     )
     ax.set_xlabel(x_axis)
     ax.set_ylabel(y_axis)
     ax.set_title(title)
     cb = fig.colorbar(im)
     cb.set_label("Counts")
 ```
 
 %% Cell type:markdown id: tags:
 
 ### Gain Bit Value
 
 %% Cell type:code id: tags:
 
 ``` python
 for i, (pdu, mod) in enumerate(zip(db_modules, karabo_da)):
     h, ex, ey = np.histogram2d(
         raw[i].flatten(),
         gain[i].flatten(),
         bins=[100, 4],
         range=[[0, 10000], [0, 4]],
     )
     do_2d_plot(
         h,
         (ex, ey),
         "Signal (ADU)",
         "Gain Bit Value (high gain=0[00], medium gain=1[01], low gain=3[11])",
         f"Module {mod} ({pdu})",
     )
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Signal Distribution ##
 
 %% Cell type:code id: tags:
 
 ``` python
 for i, (pdu, mod) in enumerate(zip(db_modules, karabo_da)):
     fig, axs = plt.subplots(nrows=2, ncols=1, figsize=(18, 10))
     corrected_flatten = corrected[i].flatten()
     for ax, hist_range in zip(axs, [(-100, 1000), (-1000, 10000)]):
         h = ax.hist(
             corrected_flatten,
             bins=1000,
             range=hist_range,
             log=True,
         )
         l = ax.set_xlabel("Signal (keV)")
         l = ax.set_ylabel("Counts")
         _ = ax.set_title(f'Module {mod} ({pdu})')
 ```
 
 %% Cell type:markdown id: tags:
 
 ### Maximum GAIN Preview
 
 %% Cell type:code id: tags:
 
 ``` python
 display(Markdown((f"#### The per pixel maximum of train {tid} of the GAIN data")))
 
 fig, ax = plt.subplots(figsize=(18, 10))
 gain_max = np.max(gain_train_cells, axis=(1, 2))
 geom.plot_data_fast(
     gain_max,
     ax=ax,
     cmap="jet",
     colorbar={'shrink': 1, 'pad': 0.01},
 )
 plt.show()
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Bad Pixels ##
 The mask contains dedicated entries for all pixels and memory cells as well as all three gains stages. Each mask entry is encoded in 32 bits as:
 
 %% Cell type:code id: tags:
 
 ``` python
 table = []
 for item in BadPixels:
     table.append(
         (item.name, f"{item.value:016b}"))
 md = display(Latex(tabulate.tabulate(
     table, tablefmt='latex',
     headers=["Bad pixel type", "Bit mask"])))
 ```
 
 %% Cell type:markdown id: tags:
 
 ### Single Image Bad Pixels ###
 
 A single image bad pixel map for the first image of the first train
 
 %% Cell type:code id: tags:
 
 ``` python
 display(Markdown(f"#### Bad pixels image for train {tid}"))
 
 fig, ax = plt.subplots(figsize=(18, 10))
 if not strixel_sensor:
     geom.plot_data_fast(
         np.log2(mask_train),
         ax=ax,
         vmin=0, vmax=32, cmap="jet",
         colorbar={'shrink': 1, 'pad': 0.01},
     )
 else:
     ax.imshow(np.log2(mask_train).squeeze(), vmin=0, vmax=32, cmap='jet', aspect=10)
 
 plt.show()
 ```

notebooks/Jungfrau/Jungfrau_dark_analysis_all_gains_burst_mode_NBC.ipynb

 %% Cell type:markdown id: tags:
 
 # Jungfrau Dark Image Characterization #
 
 Author: European XFEL Detector Group, Version: 2.0
 
 Analyzes Jungfrau dark image data to deduce offset, noise and resulting bad pixel maps
 
 %% Cell type:code id: tags:
 
 ``` python
 in_folder = '/gpfs/exfel/exp/SPB/202130/p900204/raw/'  # folder under which runs are located, required
 out_folder = '/gpfs/exfel/data/scratch/ahmedk/test/remove' # path to place reports at, required
 metadata_folder = ''  # Directory containing calibration_metadata.yml when run by xfel-calibrate
 run_high = 141 # run number for G0 dark run, required
 run_med = 142 # run number for G1 dark run, required
 run_low = 143 # run number for G2 dark run, required
 
 # Parameters used to access raw data.
 karabo_da = ['JNGFR01', 'JNGFR02','JNGFR03','JNGFR04', 'JNGFR05', 'JNGFR06','JNGFR07','JNGFR08'] # list of data aggregators, which corresponds to different JF modules
 karabo_id = 'SPB_IRDA_JF4M'  # karabo_id (detector identifier) prefix of Jungfrau detector to process.
 karabo_id_control = ''  # if control is on a different ID, set to empty string if it is the same a karabo-id
 receiver_template = 'JNGFR{:02}' # inset for receiver devices
 instrument_source_template = '{}/DET/{}:daqOutput'  # template for instrument source name (filled with karabo_id & receiver_id). e.g. 'SPB_IRDA_JF4M/DET/JNGFR01:daqOutput'
 ctrl_source_template = '{}/DET/CONTROL'  # template for control source name (filled with karabo_id_control)
 
 # Parameters for calibration database and storing constants.
 use_dir_creation_date = True  # use dir creation date
 cal_db_interface = 'tcp://max-exfl016:8016#8045'  # calibrate db interface to connect to
 cal_db_timeout = 300000 # timeout on caldb requests
 local_output = True  # output constants locally
 db_output = False  # output constants to database
 
 # Parameters affecting creating dark calibration constants.
 badpixel_threshold_sigma = 5.  # bad pixels defined by values outside n times this std from median
 offset_abs_threshold_low = [1000, 10000, 10000]  # absolute bad pixel threshold in terms of offset, lower values
 offset_abs_threshold_high = [8000, 15000, 15000]  # absolute bad pixel threshold in terms of offset, upper values
-max_trains = 0  # Maximum trains to process darks. Set to 0 to process all available train images.
-min_trains = 1  # Minimum number of trains that should be available to process dark constants. Default 1.
+max_trains = 1000  # Maximum trains to process darks. Set to 0 to process all available train images. 1000 trains is enough resolution to create the dark constants
+min_trains = 100  # Minimum number of trains to process dark constants. Raise a warning if the run has fewer trains.
 manual_slow_data = False  # if true, use manually entered bias_voltage and integration_time values
 time_limits = 0.025  # to find calibration constants later on, the integration time is allowed to vary by 0.5 us
 
 # Parameters to be used for injecting dark calibration constants.
 integration_time = 1000 # integration time in us, will be overwritten by value in file
 gain_setting = 0  # 0 for dynamic, forceswitchg1, forceswitchg2, 1 for dynamichg0, fixgain1, fixgain2. Will be overwritten by value in file
 gain_mode = 0  # 1 if medium and low runs are  fixgain1 and fixgain2, otherwise 0. It will be overwritten by value in file, if manual_slow_data
 bias_voltage = 90  # sensor bias voltage in V, will be overwritten by value in file
 memory_cells = 16  # number of memory cells
 
 # Parameters used for plotting
 detailed_report = False
 
 # TODO: this is used for only Warning check at AGIPD dark.
 # Need to rethink if it makes sense to use it here as well.
 operation_mode = 'ADAPTIVE_GAIN'  # Detector operation mode, optional
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
-import glob
 import os
 import warnings
-from pathlib import Path
+from logging import warning
 warnings.filterwarnings('ignore')
 
 import matplotlib
 import matplotlib.pyplot as plt
 import multiprocessing
 import numpy as np
 import pasha as psh
 import yaml
 from IPython.display import Markdown, display
 from extra_data import RunDirectory
 
 matplotlib.use('agg')
 %matplotlib inline
 
 from XFELDetAna.plotting.heatmap import heatmapPlot
 from XFELDetAna.plotting.histogram import histPlot
 from cal_tools import jungfraulib, step_timing
-from cal_tools.ana_tools import save_dict_to_hdf5
 from cal_tools.enums import BadPixels, JungfrauGainMode
 from cal_tools.tools import (
     get_dir_creation_date,
     get_pdu_from_db,
     get_random_db_interface,
     get_report,
     save_const_to_h5,
     send_to_db,
 )
 from iCalibrationDB import Conditions, Constants
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # Constants relevant for the analysis
 run_nums = [run_high, run_med, run_low]  # run number for G0/HG0, G1, G2
 sensor_size = (1024, 512)
 gains = [0, 1, 2]
 
 fixed_settings = [
     JungfrauGainMode.FIX_GAIN_1.value, JungfrauGainMode.FIX_GAIN_2.value]
 dynamic_settings = [
     JungfrauGainMode.FORCE_SWITCH_HG1.value, JungfrauGainMode.FORCE_SWITCH_HG2.value]
 old_fixed_settings = ["fixgain1", "fixgain2"]
 
 creation_time = None
 if use_dir_creation_date:
     creation_time = get_dir_creation_date(in_folder, run_high)
     print(f"Using {creation_time} as creation time")
 os.makedirs(out_folder, exist_ok=True)
 
 cal_db_interface = get_random_db_interface(cal_db_interface)
 print(f'Calibration database interface: {cal_db_interface}')
 
 if karabo_id_control == "":
     karabo_id_control = karabo_id
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 proposal = list(filter(None, in_folder.strip('/').split('/')))[-2]
 file_loc = f"proposal:{proposal} runs:{run_high} {run_med} {run_low}"
 
 report = get_report(metadata_folder)
 
 step_timer = step_timing.StepTimer()
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Reading control data
 
 %% Cell type:code id: tags:
 
 ``` python
 step_timer.start()
 gain_runs = dict()
 
 med_low_settings = []
 
 ctrl_src = ctrl_source_template.format(karabo_id_control)
 
 for gain, run_n in enumerate(run_nums):
     run_dc = RunDirectory(f"{in_folder}/r{run_n:04d}/")
     gain_runs[run_n] = [gain, run_dc]
     ctrl_data = jungfraulib.JungfrauCtrl(run_dc, ctrl_src)
     # Read control data for the high gain run only.
     if run_n == run_high:
 
         run_mcells, sc_start = ctrl_data.get_memory_cells()
 
         if not manual_slow_data:
             integration_time = ctrl_data.get_integration_time()
             bias_voltage = ctrl_data.get_bias_voltage()
             gain_setting = ctrl_data.get_gain_setting()
             print(f"Gain setting is {gain_setting} ({ctrl_data.run_settings})")
             print(f"Integration time is {integration_time} us")
             print(f"Bias voltage is {bias_voltage} V")
         if run_mcells == 1:
             memory_cells = 1
             print('Dark runs in single cell mode, '
                   f'storage cell start: {sc_start:02d}')
         else:
             memory_cells = 16
             print('Dark runs in burst mode, '
                   f'storage cell start: {sc_start:02d}')
     else:
         gain_mode = ctrl_data.get_gain_mode()
         med_low_settings.append(ctrl_data.run_mode)
 
 # A transperent workaround for old raw data with wrong/missing medium and low settings
 if med_low_settings == [None, None]:
-    print("WARNING: run.settings is not stored in the data to read. "
-          f"Hence assuming gain_mode = {gain_mode} for adaptive old data.")
+    warning("run.settings is not stored in the data to read. "
+            f"Hence assuming gain_mode = {gain_mode} for adaptive old data.")
 elif med_low_settings == ["dynamicgain", "forceswitchg1"]:
-    print(f"WARNING: run.settings for medium and low gain runs are wrong {med_low_settings}. "
-          f"This is an expected bug for old raw data. Setting gain_mode to {gain_mode}.")
+    warning(f"run.settings for medium and low gain runs are wrong {med_low_settings}. "
+            f"This is an expected bug for old raw data. Setting gain_mode to {gain_mode}.")
 # Validate that low_med_settings is not a mix of adaptive and fixed settings.
 elif not (sorted(med_low_settings) in [fixed_settings, dynamic_settings, old_fixed_settings]):  # noqa
     raise ValueError(
         "Medium and low run settings are not as expected. "
         f"Either {dynamic_settings}, {fixed_settings}, or {old_fixed_settings} are expected.\n"
         f"Got {sorted(med_low_settings)} for both runs, respectively.")
 
 print(f"Gain mode is {gain_mode} ({med_low_settings})")
 
 step_timer.done_step(f'Reading control data.')
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # set the operating condition
 condition = Conditions.Dark.jungfrau(
     memory_cells=memory_cells,
     bias_voltage=bias_voltage,
     integration_time=integration_time,
     gain_setting=gain_setting,
     gain_mode=gain_mode,
 )
 
 db_modules = get_pdu_from_db(
     karabo_id=karabo_id,
     karabo_da=karabo_da,
     constant=Constants.jungfrau.Offset(),
     condition=condition,
     cal_db_interface=cal_db_interface,
     snapshot_at=creation_time)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # Start retrieving existing constants for comparison
 mod_x_const = [(mod, const) for const in ["Offset", "Noise", "BadPixelsDark"] for mod in karabo_da]
 
 from cal_tools.tools import get_from_db
 from datetime import timedelta
 
 def retrieve_old_constant(mod, const):
     dconst = getattr(Constants.jungfrau, const)()
 
     data, mdata = get_from_db(
         karabo_id=karabo_id,
         karabo_da=mod,
         constant=dconst,
         condition=condition,
         empty_constant=None,
         cal_db_interface=cal_db_interface,
         creation_time=creation_time-timedelta(seconds=60) if creation_time else None,
         strategy="pdu_prior_in_time",
         verbosity=1,
         timeout=cal_db_timeout
     )
 
     if mdata is None or data is None:
         timestamp = "Not found"
         filepath = None
         h5path = None
     else:
         timestamp = mdata.calibration_constant_version.begin_at.isoformat()
         filepath = os.path.join(
             mdata.calibration_constant_version.hdf5path,
             mdata.calibration_constant_version.filename
         )
         h5path = mdata.calibration_constant_version.h5path
 
     return data, timestamp, filepath, h5path
 
 
 old_retrieval_pool = multiprocessing.Pool()
 old_retrieval_res = old_retrieval_pool.starmap_async(
     retrieve_old_constant, mod_x_const
 )
 old_retrieval_pool.close()
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # Use only high gain threshold for all gains in case of fixed_gain.
 
 if gain_mode:  # fixed_gain
     offset_abs_threshold = [[offset_abs_threshold_low[0]]*3, [offset_abs_threshold_high[0]]*3]
 else:
     offset_abs_threshold = [offset_abs_threshold_low, offset_abs_threshold_high]
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
-context = psh.context.ThreadContext(num_workers=multiprocessing.cpu_count())
+context = psh.context.ThreadContext(num_workers=memory_cells)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 """
 All jungfrau runs are taken through one acquisition, except for the forceswitch runs.
 While taking non-fixed dark runs, a procedure of multiple acquisitions is used to switch the storage cell indices.
 
 This is done for medium and low gain dark dynamic runs, only [forceswitchg1, forceswitchg2]:
 Switching the cell indices in burst mode is a work around for hardware procedure
 deficiency that produces wrong data for dark runs except for the first storage cell.
 This is why multiple acquisitions are taken to switch the used storage cells and
 acquire data through two cells for each of the 16 cells instead of acquiring darks through all 16 cells.
 """
 
 print(f"Maximum trains to process is set to {max_trains}")
 
 noise_map = dict()
 offset_map = dict()
 bad_pixels_map = dict()
 
 for mod in karabo_da:
     step_timer.start()
     instrument_src = instrument_source_template.format(
         karabo_id, receiver_template.format(int(mod[-2:])))
 
     print(f"\n- Instrument data path for {mod} is {instrument_src}.")
 
-    offset_map[mod] = context.alloc(shape=(sensor_size+(memory_cells, 3)), fill=0)
+    offset_map[mod] = context.alloc(
+        shape=(sensor_size+(memory_cells, 3)), fill=0, dtype=np.float32)
     noise_map[mod] = context.alloc(like=offset_map[mod], fill=0)
     bad_pixels_map[mod] = context.alloc(like=offset_map[mod], dtype=np.uint32, fill=0)
 
     for run_n, [gain, run_dc] in gain_runs.items():
 
         def process_cell(worker_id, array_index, cell_number):
             cell_slice_idx = acelltable == cell_number
-            thiscell = images[..., cell_slice_idx]
+            thiscell = images[..., cell_slice_idx]  # [1024, 512, n_trains]
 
             # Identify cells/trains with images of 0 pixels.
             # TODO: An investigation is ongoing by DET to identify reason for these empty images.
-            nonzero_adc = np.any(thiscell != 0 , axis=(0, 1))
+            nonzero_adc = np.any(thiscell != 0 , axis=(0, 1))  # [n_trains]
 
             # Exclude empty images with 0 pixels, before calculating offset and noise
             thiscell = thiscell[..., nonzero_adc]
-            offset_map[mod][..., cell_number, gain] = np.mean(thiscell, axis=2)
-            noise_map[mod][..., cell_number, gain] = np.std(thiscell, axis=2)
-
+            offset_map[mod][..., cell_number, gain] = np.mean(  # [1024, 512]
+                thiscell, axis=2, dtype=np.float32)
+            noise_map[mod][..., cell_number, gain] = np.std(  # [1024, 512]
+                thiscell, axis=2, dtype=np.float32)
+            del thiscell
             # Check if there are wrong bad gain values.
             # 1. Exclude empty images.
             # 2. Indicate pixels with wrong gain value for any train for each cell.
             # TODO: mean is used to use thresholds for accepting gain values, even if not 0 mean value.
-            gain_avg = np.mean(
-                gain_vals[..., cell_slice_idx][..., nonzero_adc], axis=2)
+            gain_avg = np.mean(  # [1024, 512]
+                gain_vals[..., cell_slice_idx][..., nonzero_adc],
+                axis=2, dtype=np.float32
+            )
 
+            # [1024, 512]
             bad_pixels_map[mod][..., cell_number, gain][gain_avg != raw_g] |= BadPixels.WRONG_GAIN_VALUE.value
+
         print(f"Gain stage {gain}, run {run_n}")
 
         # load shape of data for memory cells, and detector size (imgs, cells, x, y)
-        n_imgs = run_dc[instrument_src, "data.adc"].shape[0]
+        n_trains = run_dc[instrument_src, "data.adc"].shape[0]
         # load number of data available, including trains with empty data.
-        n_trains = len(run_dc.train_ids)
+        all_trains = len(run_dc.train_ids)
         instr_dc = run_dc.select(instrument_src, require_all=True)
-        empty_trains = n_trains - n_imgs
+        empty_trains = all_trains - n_trains
         if empty_trains != 0:
-            print(f"\tWARNING: {mod} has {empty_trains} trains with empty data out of {n_trains} trains")  # noqa
+            print(f"{mod} has {empty_trains} empty trains out of {all_trains} trains")
         if max_trains > 0:
-            n_imgs = min(n_imgs, max_trains)
-        print(f"Processing {n_imgs} images.")
-        # Select only requested number of images to process darks.
-        instr_dc = instr_dc.select_trains(np.s_[:n_imgs])
+            n_trains = min(n_trains, max_trains)
+        print(f"Processing {n_trains} images.")
 
-        if n_imgs < min_trains:
-            raise ValueError(
-                f"Less than {min_trains} trains are available in RAW data."
-                " Not enough data to process darks.")
+        if n_trains == 0:
+            raise ValueError(f"{run_n} has no trains to process.")
 
+        if n_trains < min_trains:
+            warning(f"Less than {min_trains} trains are available in RAW data.")
+
+        # Select only requested number of images to process darks.
+        instr_dc = instr_dc.select_trains(np.s_[:n_trains])
         images = np.transpose(
             instr_dc[instrument_src, "data.adc"].ndarray(), (3, 2, 1, 0))
         acelltable = np.transpose(instr_dc[instrument_src, "data.memoryCell"].ndarray())
         gain_vals = np.transpose(
             instr_dc[instrument_src, "data.gain"].ndarray(), (3, 2, 1, 0))
 
         # define gain value as saved in raw gain map
         raw_g = 3 if gain == 2 else gain
 
         if memory_cells == 1:
             acelltable -= sc_start
         # Only for dynamic medium and low gain runs [forceswitchg1, forceswitchg2] in burst mode.
 
         if gain_mode == 0 and gain > 0 and memory_cells == 16:
             # 255 similar to the receiver which uses the 255
             # value to indicate a cell without an image.
             # image shape for forceswitchg1 and forceswitchg2 = (1024, 512, 2, trains)
             # compared to expected shape of (1024, 512, 16, trains) for high gain run.
             acelltable[1:] = 255
 
         # Calculate offset and noise maps
         context.map(process_cell, range(memory_cells))
-
+        del images
+        del acelltable
+        del gain_vals
     step_timer.done_step(f'Creating Offset and noise constants for a module.')
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 if detailed_report:
     display(Markdown("## Offset and Noise Maps:"))
     display(Markdown(
         "Below offset and noise maps for the high ($g_0$) gain stage are shown, "
         "alongside the distribution of these values. One expects block-like "
         "structures mapping to the ASICs of the detector"))
     g_name = ['G0', 'G1', 'G2']
     g_range = [(0, 8000), (8000, 16000), (8000, 16000)]
     n_range = [(0., 50.), (0., 50.), (0., 50.)]
 
     unit = '[ADCu]'
     # TODO: Fix plots arrangment and speed for Jungfrau burst mode.
     step_timer.start()
     for pdu, mod in zip(db_modules, karabo_da):
         for g_idx in gains:
             for cell in range(0, memory_cells):
                 f_o0 = heatmapPlot(
                     np.swapaxes(offset_map[mod][..., cell, g_idx], 0, 1),
                     y_label="Row",
                     x_label="Column",
                     lut_label=unit,
                     aspect=1.,
                     vmin=g_range[g_idx][0],
                     vmax=g_range[g_idx][1],
                     title=f'Pedestal {g_name[g_idx]} - Cell {cell:02d} - Module {mod} ({pdu})')
 
                 fo0, ax_o0 = plt.subplots()
                 res_o0 = histPlot(
                     ax_o0, offset_map[mod][..., cell, g_idx],
                     bins=800,
                     range=g_range[g_idx],
                     facecolor='b',
                     histotype='stepfilled',
                 )
 
                 ax_o0.tick_params(axis='both',which='major',labelsize=15)
                 ax_o0.set_title(
                     f'Module pedestal distribution - Cell {cell:02d} - Module {mod} ({pdu})',
                     fontsize=15)
                 ax_o0.set_xlabel(f'Pedestal {g_name[g_idx]} {unit}',fontsize=15)
                 ax_o0.set_yscale('log')
 
                 f_n0 = heatmapPlot(
                     np.swapaxes(noise_map[mod][..., cell, g_idx], 0, 1),
                     y_label="Row",
                     x_label="Column",
                     lut_label= unit,
                     aspect=1.,
                     vmin=n_range[g_idx][0],
                     vmax=n_range[g_idx][1],
                     title=f"RMS noise {g_name[g_idx]} - Cell {cell:02d} - Module {mod} ({pdu})",
                 )
 
                 fn0, ax_n0 = plt.subplots()
                 res_n0 = histPlot(
                     ax_n0,
                     noise_map[mod][..., cell, g_idx],
                     bins=100,
                     range=n_range[g_idx],
                     facecolor='b',
                     histotype='stepfilled',
                 )
 
                 ax_n0.tick_params(axis='both', which='major', labelsize=15)
                 ax_n0.set_title(
                     f'Module noise distribution - Cell {cell:02d} - Module {mod} ({pdu})',
                     fontsize=15)
                 ax_n0.set_xlabel(
                     f'RMS noise {g_name[g_idx]} ' + unit, fontsize=15)
                 plt.show()
     step_timer.done_step(f'Plotting offset and noise maps.')
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Bad Pixel Map ###
 
 The bad pixel map is deduced by comparing offset and noise of each pixel ($v_i$) and each gain ($g$) against the median value for that gain stage:
 
 $$
 v_i > \mathrm{median}(v_{k,g}) + n \sigma_{v_{k,g}}
 $$
 or
 $$
 v_i < \mathrm{median}(v_{k,g}) - n \sigma_{v_{k,g}}
 $$
 
 Values are encoded in a 32 bit mask, where for the dark image deduced bad pixels the following non-zero entries are relevant:
 
 %% Cell type:code id: tags:
 
 ``` python
 def print_bp_entry(bp):
     print("{:<30s} {:032b} -> {}".format(bp.name, bp.value, int(bp.value)))
 
 print_bp_entry(BadPixels.OFFSET_OUT_OF_THRESHOLD)
 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, 1))[None, None, :, :]
     std = np.nanstd(d, axis=(0, 1))[None, None, :, :]
     idx = (d > badpixel_threshold_sigma*std+mdn) | (d < (-badpixel_threshold_sigma)*std+mdn)
 
     return idx
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 step_timer.start()
 
 for pdu, mod in zip(db_modules, karabo_da):
     display(Markdown(f"### Badpixels for module {mod} ({pdu}):"))
     offset_abs_threshold = np.array(offset_abs_threshold)
 
     bad_pixels_map[mod][eval_bpidx(offset_map[mod])] |= BadPixels.OFFSET_OUT_OF_THRESHOLD.value
 
     bad_pixels_map[mod][~np.isfinite(offset_map[mod])] |= BadPixels.OFFSET_NOISE_EVAL_ERROR.value
 
     bad_pixels_map[mod][eval_bpidx(noise_map[mod])] |= BadPixels.NOISE_OUT_OF_THRESHOLD.value
 
     bad_pixels_map[mod][~np.isfinite(noise_map[mod])] |= BadPixels.OFFSET_NOISE_EVAL_ERROR.value
 
     bad_pixels_map[mod][(offset_map[mod] < offset_abs_threshold[0][None, None, None, :]) | (offset_map[mod] > offset_abs_threshold[1][None, None, None, :])] |= BadPixels.OFFSET_OUT_OF_THRESHOLD.value  # noqa
 
     if detailed_report:
 
         for g_idx in gains:
             for cell in range(memory_cells):
                 bad_pixels = bad_pixels_map[mod][:, :, cell, g_idx]
                 fn_0 = heatmapPlot(
                     np.swapaxes(bad_pixels, 0, 1),
                     y_label="Row",
                     x_label="Column",
                     lut_label=f"Badpixels {g_name[g_idx]} [ADCu]",
                     aspect=1.,
                     vmin=0, vmax=5,
                     title=f'G{g_idx} Bad pixel map - Cell {cell:02d} - Module {mod} ({pdu})')
 step_timer.done_step(f'Creating bad pixels constant')
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Inject and save calibration constants
 
 %% Cell type:code id: tags:
 
 ``` python
 step_timer.start()
 for mod, db_mod in zip(karabo_da, db_modules):
     constants = {
         'Offset': np.moveaxis(offset_map[mod], 0, 1),
         'Noise': np.moveaxis(noise_map[mod], 0, 1),
         'BadPixelsDark': np.moveaxis(bad_pixels_map[mod], 0, 1),
     }
 
     md = None
 
     for key, const_data in constants.items():
 
         const =  getattr(Constants.jungfrau, key)()
         const.data = const_data
 
         for parm in condition.parameters:
             if parm.name == "Integration Time":
                 parm.lower_deviation = time_limits
                 parm.upper_deviation = time_limits
 
         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"• Memory cells: {memory_cells}\n"
     f"• Integration time: {integration_time}\n"
     f"• Gain setting: {gain_setting}\n"
     f"• Gain mode: {gain_mode}\n"
     f"• Creation time: {md.calibration_constant_version.begin_at if md is not None else creation_time}\n")  # noqa
 step_timer.done_step("Injecting constants.")
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 print(f"Total processing time {step_timer.timespan():.01f} s")
 step_timer.print_summary()
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # now we need the old constants
 old_const = {}
 old_mdata = {}
 old_retrieval_res.wait()
 
 for (mod, const), (data, timestamp, filepath, h5path) in zip(
     mod_x_const, old_retrieval_res.get()):
     old_const.setdefault(mod, {})[const] = data
     old_mdata.setdefault(mod, {})[const] = {
         "timestamp": timestamp,
         "filepath": filepath,
         "h5path": h5path,
     }
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 display(Markdown("## The following pre-existing constants are used for comparison:"))
 
 for mod, consts in old_mdata.items():
     pdu = db_modules[karabo_da.index(mod)]
     display(Markdown(f"- {mod} ({pdu})"))
     for const in consts:
         display(Markdown(f"    - {const} at {consts[const]['timestamp']}"))
     # saving locations of old constants for summary notebook
     with open(f"{metadata_folder or out_folder}/module_metadata_{mod}.yml", "w") as fd:
         yaml.safe_dump(
             {
                 "module": mod,
                 "pdu": pdu,
                 "old-constants": old_mdata[mod],
             },
             fd,
         )
 ```

notebooks/LPD/LPDChar_Darks_NBC.ipynb

 %% Cell type:markdown id: tags:
 
 # LPD Offset, Noise and Dead Pixels Characterization #
 
 Author: M. Karnevskiy, S. Hauf
 
 This notebook performs re-characterize of dark images to derive offset, noise and bad-pixel maps. All three types of constants are evaluated per-pixel and per-memory cell.
 
 The notebook will correctly handle veto settings, but note that if you veto cells you will not be able to use these offsets for runs with different veto settings - vetoed cells will have zero offset.
 
 The evaluated calibration constants are stored locally and injected in the calibration data base.
 
 
 %% Cell type:code id: tags:
 
 ``` python
 in_folder = "/gpfs/exfel/exp/FXE/202030/p900121/raw" # path to input data, required
 out_folder = "/gpfs/exfel/data/scratch/ahmedk/test/LPD/" # path to output to, required
 metadata_folder = ""  # Directory containing calibration_metadata.yml when run by xfel-calibrate
 sequence = 0 # sequence files to evaluate
 modules = [-1] # list of modules to evaluate, RANGE ALLOWED
 run_high = 120 # run number in which high gain data was recorded, required
 run_med = 121 # run number in which medium gain data was recorded, required
 run_low = 122 # run number in which low gain data was recorded, required
 
 karabo_id = "FXE_DET_LPD1M-1" # karabo karabo_id
 karabo_da = ['-1']  # a list of data aggregators names, Default [-1] for selecting all data aggregators
 receiver_id = "{}CH0" # inset for receiver devices
 path_template = 'RAW-R{:04d}-{}-S{:05d}.h5' # the template to use to access data
 h5path = '/INSTRUMENT/{}/DET/{}:xtdf/image' # path in the HDF5 file to images
 h5path_idx = '/INDEX/{}/DET/{}:xtdf/image' # path in the HDF5 file to images
 
 use_dir_creation_date = True # use the creation date of the directory for database time derivation
 cal_db_interface = "tcp://max-exfl016:8015#8025" # the database interface to use
 cal_db_timeout = 300000 # timeout on caldb requests"
 local_output = True # output constants locally
 db_output = False # output constants to database
 
 capacitor_setting = 5 # capacitor_setting for which data was taken
 mem_cells = 512 # number of memory cells used
 bias_voltage = 250 # detector bias voltage
 thresholds_offset_sigma = 3. # bad pixel relative threshold in terms of n sigma offset
 thresholds_offset_hard = [400, 1500] # bad pixel hard threshold
 thresholds_noise_sigma = 7. # bad pixel relative threshold in terms of n sigma noise
 thresholds_noise_hard = [1, 35] # bad pixel hard threshold
 skip_first_ntrains = 10 # Number of first trains to skip
 
+# Parameters for plotting
+skip_plots = False  # exit after writing corrected files
+
 instrument = "FXE" # instrument name
 ntrains = 100 # number of trains to use
 high_res_badpix_3d = False # plot bad-pixel summary in high resolution
 test_for_normality = False # permorm normality test
+inject_cell_order = False  # Include memory cell order as part of the detector condition
 operation_mode = ''  # Detector operation mode, optional
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 import copy
 import multiprocessing
 import os
 import warnings
 from collections import OrderedDict
 from datetime import datetime
 from functools import partial
+from logging import warning
 
 warnings.filterwarnings('ignore')
 
 import dateutil.parser
 import h5py
 import matplotlib
 import pasha as psh
 import scipy.stats
 from IPython.display import Latex, Markdown, display
 
 matplotlib.use("agg")
 import matplotlib.patches as patches
 import matplotlib.pyplot as plt
 
 %matplotlib inline
 import numpy as np
 import tabulate
 import yaml
 from iCalibrationDB import Conditions, Constants, Detectors, Versions
 from XFELDetAna.plotting.heatmap import heatmapPlot
 from XFELDetAna.plotting.simpleplot import simplePlot
 
 from cal_tools.enums import BadPixels
 from cal_tools.plotting import (
     create_constant_overview,
     plot_badpix_3d,
     show_overview,
     show_processed_modules,
 )
 from cal_tools.tools import (
     get_dir_creation_date,
     get_from_db,
     get_notebook_name,
     get_pdu_from_db,
     get_random_db_interface,
     get_report,
     map_gain_stages,
     module_index_to_qm,
     parse_runs,
     run_prop_seq_from_path,
     save_const_to_h5,
     send_to_db,
 )
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 gains = np.arange(3)
 max_cells = mem_cells
 cells = np.arange(max_cells)
 gain_names = ['High', 'Medium', 'Low']
 
 if karabo_da[0] == '-1':
     if modules[0] == -1:
         modules = list(range(16))
     karabo_da = ['LPD{:02d}'.format(i) for i in modules]
 else:
     modules = [int(x[-2:]) for x in karabo_da]
 
 gain_runs = OrderedDict()
 if capacitor_setting == 5:
     gain_runs["high_5pf"] = run_high
     gain_runs["med_5pf"] =  run_med
     gain_runs["low_5pf"] =  run_low
 elif capacitor_setting == 50:
     gain_runs["high_50pf"] = run_high
     gain_runs["med_50pf"] =  run_med
     gain_runs["low_50pf"] =  run_low
 
 capacitor_settings = [capacitor_setting]
 capacitor_settings = ['{}pf'.format(c) for c in capacitor_settings]
 
 h5path = h5path.format(karabo_id, receiver_id)
 h5path_idx = h5path_idx.format(karabo_id, receiver_id)
 
 creation_time = None
 if use_dir_creation_date:
     creation_time = get_dir_creation_date(in_folder, run_high)
     print("Using {} as creation time".format(creation_time))
 
 run, prop, seq = run_prop_seq_from_path(in_folder)
 
 cal_db_interface = get_random_db_interface(cal_db_interface)
 
 display(Markdown('## Evaluated parameters'))
 print('CalDB Interface {}'.format(cal_db_interface))
 print("Proposal: {}".format(prop))
 print("Memory cells: {}/{}".format(mem_cells, max_cells))
 print("Runs: {}, {}, {}".format(run_high, run_med, run_low))
 print("Sequence: {}".format(sequence))
 print("Using DB: {}".format(db_output))
 print("Input: {}".format(in_folder))
 print("Output: {}".format(out_folder))
 print("Bias voltage: {}V".format(bias_voltage))
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # set everything up filewise
 gmf = map_gain_stages(in_folder, gain_runs, path_template, karabo_da, [sequence])
 gain_mapped_files, total_sequences, total_file_size = gmf
 print(f"Will process a total of {total_sequences} files.")
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Data processing
 
 %% Cell type:code id: tags:
 
 ``` python
 parallel_num_procs = min(6, len(modules)*3)
 parallel_num_threads = multiprocessing.cpu_count() // parallel_num_procs
 
 # the actual characterization
 def characterize_module(filename, channel, gg, cap):
 
     def splitOffGainLPD(d):
         msk = np.zeros(d.shape, np.uint16)
         msk[...] = 0b0000111111111111
         data = np.bitwise_and(d, msk)
         msk[...] = 0b0011000000000000
         gain = np.bitwise_and(d, msk)//4096
         gain[gain > 2] = 2
         return data, gain
 
     infile = h5py.File(filename, "r")
 
     instrument_src = h5path.format(channel)
     index_src = h5path_idx.format(channel)
     count = infile[f"{index_src}/count"][()]
     first = infile[f"{index_src}/first"][()]
     valid = count != 0
     count, first = count[valid], first[valid]
     first_image = int(first[skip_first_ntrains] if first.shape[0] > skip_first_ntrains else 0)
     last_image = int(first_image + np.sum(count[skip_first_ntrains:skip_first_ntrains+ntrains]))
 
     im = np.array(infile["{}/data".format(instrument_src, channel)][first_image:last_image, ...])
     cellid = np.squeeze(np.array(infile["{}/cellId".format(instrument_src, channel)][first_image:last_image, ...]))
     infile.close()
+    if im.shape[0] == 0:  # No data
+        return None, None, channel, gg, cap, None, None, None, None
+
+    cellid_pattern = cellid[:count[0]]
 
     im, g = splitOffGainLPD(im[:, 0, ...])
     im = im.astype(np.float32)
 
     im = np.rollaxis(im, 2)
     im = np.rollaxis(im, 2, 1)
 
     context = psh.context.ThreadContext(num_workers=parallel_num_threads)
     offset = context.alloc(shape=(im.shape[0], im.shape[1], max_cells), dtype=np.float64)
     noise = context.alloc(like=offset)
     normal_test = context.alloc(like=offset)
     def process_cell(worker_id, array_index, cc):
         idx = cellid == cc
         im_slice = im[..., idx]
         if np.any(idx):
             offset[..., cc] = np.median(im_slice, axis=2)
             noise[..., cc] = np.std(im_slice, axis=2)
             if test_for_normality:
                 _, normal_test[..., cc] = scipy.stats.normaltest(
                     im[:, :, idx], axis=2)
     context.map(process_cell, np.unique(cellid))
 
     # bad pixels
     bp = np.zeros(offset.shape, np.uint32)
     # offset related bad pixels
     offset_mn = np.nanmedian(offset, axis=(0, 1))
     offset_std = np.nanstd(offset, axis=(0, 1))
 
     bp[(offset < offset_mn-thresholds_offset_sigma*offset_std) |
        (offset > offset_mn+thresholds_offset_sigma*offset_std)] |= BadPixels.OFFSET_OUT_OF_THRESHOLD.value
     bp[(offset < thresholds_offset_hard[0]) | (
         offset > thresholds_offset_hard[1])] |= BadPixels.OFFSET_OUT_OF_THRESHOLD.value
     bp[~np.isfinite(offset)] |= BadPixels.OFFSET_NOISE_EVAL_ERROR.value
 
     # noise related bad pixels
     noise_mn = np.nanmedian(noise, axis=(0, 1))
     noise_std = np.nanstd(noise, axis=(0, 1))
 
     bp[(noise < noise_mn-thresholds_noise_sigma*noise_std) |
        (noise > noise_mn+thresholds_noise_sigma*noise_std)] |= BadPixels.NOISE_OUT_OF_THRESHOLD.value
     bp[(noise < thresholds_noise_hard[0]) | (
         noise > thresholds_noise_hard[1])] |= BadPixels.NOISE_OUT_OF_THRESHOLD.value
     bp[~np.isfinite(noise)] |= BadPixels.OFFSET_NOISE_EVAL_ERROR.value
 
-    idx = cellid == 12
-    return offset, noise, channel, gg, cap, bp, im[12, 12, idx], normal_test
+    idx = (cellid == cellid[0])
+    return offset, noise, channel, gg, cap, bp, im[12, 12, idx], normal_test, cellid_pattern
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 offset_g = OrderedDict()
 noise_g = OrderedDict()
 badpix_g = OrderedDict()
 data_g = OrderedDict()
 ntest_g = OrderedDict()
+# Should be the same cell order for all modules & all gain stages
+cellid_patterns_g = {}
 
 gg = 0
 old_cap = None
 start = datetime.now()
 inp = []
 
 for gain, mapped_files in gain_mapped_files.items():
     cap = gain.split("_")[1]
     if cap != old_cap:
         gg = 0
         old_cap = cap
         offset_g[cap] = OrderedDict()
         noise_g[cap] = OrderedDict()
         badpix_g[cap] = OrderedDict()
         data_g[cap] = OrderedDict()
         ntest_g[cap] = OrderedDict()
+        cellid_patterns_g[cap] = {}
 
     for i in modules:
         qm = module_index_to_qm(i)
         if qm in mapped_files and not mapped_files[qm].empty():
             fname_in = mapped_files[qm].get()
             print("Process file: ", fname_in)
             inp.append((fname_in, i, gg, cap))
 
     gg+=1
 
 with multiprocessing.Pool(processes=parallel_num_procs) as pool:
     results = pool.starmap(characterize_module, inp)
 
 for ir, r in enumerate(results):
-    offset, noise, i, gg, cap, bp, data, normal = r
+    offset, noise, i, gg, cap, bp, data, normal, cellid_pattern = r
+    if data is None:
+        warning(f"No data for module {i} of gain {gg}")
+        skip_plots = True
+        continue
     qm = module_index_to_qm(i)
     if qm not in offset_g[cap]:
         offset_g[cap][qm] = np.zeros(
             (offset.shape[0], offset.shape[1], offset.shape[2], 3))
         noise_g[cap][qm] = np.zeros_like(offset_g[cap][qm])
         badpix_g[cap][qm] = np.zeros_like(offset_g[cap][qm], dtype=np.uint32)
         data_g[cap][qm] = np.full((ntrains, 3), np.nan)
         ntest_g[cap][qm] = np.zeros_like(offset_g[cap][qm])
+        cellid_patterns_g[cap][qm] = cellid_pattern
+    else:
+        if not np.array_equal(cellid_pattern, cellid_patterns_g[cap][qm]):
+            raise ValueError("Inconsistent cell ID pattern between gain stages")
+
 
     offset_g[cap][qm][..., gg] = offset
     noise_g[cap][qm][..., gg] = noise
     badpix_g[cap][qm][..., gg] = bp
     data_g[cap][qm][:data.shape[0], gg] = data
     ntest_g[cap][qm][..., gg] = normal
 
     hn, cn = np.histogram(data, bins=20)
     print(f"{gain_names[gg]} gain, Capacitor {cap}, Module: {qm}. "
           f"Number of processed trains per cell: {data.shape[0]}.")
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
+if skip_plots:
+    import sys
+    print('Skipping plots')
+    sys.exit(0)
+```
+
+%% Cell type:code id: tags:
+
+``` python
 # Read report path and create file location tuple to add with the injection
 proposal = list(filter(None, in_folder.strip('/').split('/')))[-2]
 file_loc = 'proposal:{} runs:{} {} {}'.format(proposal, run_low, run_med, run_high)
 
 report = get_report(metadata_folder)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # TODO: add db_module when received from myMDC
 # Create the modules dict of karabo_das and PDUs
 qm_dict = OrderedDict()
 for i, k_da in zip(modules, karabo_da):
     qm = module_index_to_qm(i)
     qm_dict[qm] = {"karabo_da": k_da,
                    "db_module": ""}
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # Retrieve existing constants for comparison
 clist = ["Offset", "Noise", "BadPixelsDark"]
 old_const = {}
 old_mdata = {}
 
 dinstance = "LPD1M1"
 detinst = getattr(Detectors, dinstance)
 print('Retrieve pre-existing constants for comparison.')
 for cap in capacitor_settings:
     old_const[cap] = {}
     old_mdata[cap] = {}
     for qm in offset_g[cap].keys():
         old_const[cap][qm] = {}
         old_mdata[cap][qm] = {}
         qm_db = qm_dict[qm]
         karabo_da = qm_db["karabo_da"]
+        cellid_pattern = cellid_patterns_g[cap][qm]
+        if inject_cell_order:
+            mem_cell_order = ",".join([str(c) for c in cellid_pattern]) + ","
+        else:
+            mem_cell_order = None
 
         condition = Conditions.Dark.LPD(memory_cells=max_cells,
                                         bias_voltage=bias_voltage,
-                                        capacitor=cap)
+                                        capacitor=cap,
+                                        memory_cell_order=mem_cell_order,
+                                       )
         for const in clist:
             constant = getattr(Constants.LPD, const)()
             if not qm_db["db_module"]:
                 # This should be used in case of running notebook
                 # by a different method other than myMDC which already
                 # sends CalCat info.
                 qm_db["db_module"] = get_pdu_from_db(karabo_id, [karabo_da], constant,
                                                      condition, cal_db_interface,
                                                      snapshot_at=creation_time)[0]
 
             data, mdata = get_from_db(karabo_id, karabo_da,
                                       constant,
                                       condition, None,
                                       cal_db_interface,
                                       creation_time=creation_time,
                                       verbosity=2, timeout=cal_db_timeout)
 
             old_const[cap][qm][const] = data
 
             if mdata is None or data is None:
                 old_mdata[cap][qm][const] = {
                     "timestamp": "Not found",
                     "filepath": None,
                     "h5path": None
                 }
             else:
                 timestamp = mdata.calibration_constant_version.begin_at.isoformat()
                 filepath = os.path.join(
                     mdata.calibration_constant_version.hdf5path,
                     mdata.calibration_constant_version.filename
                 )
                 h5path = mdata.calibration_constant_version.h5path
                 old_mdata[cap][qm][const] = {
                     "timestamp": timestamp,
                     "filepath": filepath,
                     "h5path": h5path
                 }
 
         with open(f"{out_folder}/module_metadata_{qm}.yml","w") as fd:
             yaml.safe_dump(
                 {
                     "module": qm,
                     "pdu": qm_db["db_module"],
                     "old-constants": old_mdata[cap][qm]
                 }, fd)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 res = OrderedDict()
 for cap in capacitor_settings:
     res[cap] = OrderedDict()
     for i in modules:
         qm = module_index_to_qm(i)
 
         res[cap][qm] = {'Offset': offset_g[cap][qm],
                         'Noise': noise_g[cap][qm],
                         'BadPixelsDark': badpix_g[cap][qm]
                         }
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # Save constants in the calibration DB
 md = None
 for cap in capacitor_settings:
     for qm in res[cap]:
 
         karabo_da = qm_dict[qm]["karabo_da"]
         db_module = qm_dict[qm]["db_module"]
+        cellid_pattern = cellid_patterns_g[cap][qm]
+        if inject_cell_order:
+            mem_cell_order = ",".join([str(c) for c in cellid_pattern]) + ","
+        else:
+            mem_cell_order = None
 
         # Do not store empty constants
         # In case of 0 trains data_g is initiated with nans and never refilled.
         if np.count_nonzero(~np.isnan(data_g[cap][qm]))==0:
+            print(f"Constant ({cap}, {qm}) would be empty, skipping saving")
             continue
+
         for const in res[cap][qm]:
 
             dconst = getattr(Constants.LPD, const)()
             dconst.data = res[cap][qm][const]
 
             # set the operating condition
+
             condition = Conditions.Dark.LPD(memory_cells=max_cells,
-                                            bias_voltage=bias_voltage,
-                                            capacitor=cap)
+                                        bias_voltage=bias_voltage,
+                                        capacitor=cap,
+                                        memory_cell_order=mem_cell_order,
+                                       )
 
             if db_output:
                 md = send_to_db(db_module, karabo_id, dconst, condition,
                                 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, karabo_id, dconst, condition,
                                       dconst.data, file_loc, report, creation_time, out_folder)
                 print(f"Calibration constant {const} is stored locally.\n")
 
         print("Constants parameter conditions are:\n")
-        print(f"• memory_cells: {max_cells}\n• bias_voltage: {bias_voltage}\n"
+        print(f"• memory_cells: {max_cells}\n"
+              f"• bias_voltage: {bias_voltage}\n"
               f"• capacitor: {cap}\n"
+              f"• memory cell order: {mem_cell_order}\n"
               f"• creation_time: {md.calibration_constant_version.begin_at if md is not None else creation_time}\n")
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 show_processed_modules(
     dinstance=dinstance,
     constants=None,
     mnames=[module_index_to_qm(i) for i in modules],
     mode="position"
 )
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Raw pedestal distribution ##
 
 Distribution of a pedestal (ADUs) over trains for the pixel (12,12), memory cell 12. A median of the distribution is shown in yellow. A standard deviation is shown in red. The green line shows average over all pixels for a given memory cell and gain stage.
 
 %% Cell type:code id: tags:
 
 ``` python
 fig, grid = plt.subplots(3, 1, sharex="col", sharey="row", figsize=(10, 7))
 fig.subplots_adjust(wspace=0, hspace=0)
 
 for cap in capacitor_settings:
     for i in modules:
         qm = module_index_to_qm(i)
         if np.count_nonzero(~np.isnan(data_g[cap][qm])) == 0:
             break
         for gain in range(3):
             data = data_g[cap][qm][:, gain]
             offset = np.nanmedian(data)
             noise = np.nanstd(data)
             xrange = [np.nanmin(data_g[cap][qm]), np.nanmax(data_g[cap][qm])]
             if xrange[1] == xrange[0]:
                 xrange = [0, xrange[0]+xrange[0]//2]
                 nbins = data_g[cap][qm].shape[0]
             else:
                 nbins = int(xrange[1] - xrange[0])
 
             hn, cn = np.histogram(data, bins=nbins, range=xrange)
 
             grid[gain].hist(data, range=xrange, bins=nbins)
             grid[gain].plot([offset-noise, offset-noise], [0, np.nanmax(hn)],
                             linewidth=1.5, color='red',
                             label='1 $\sigma$ deviation')
             grid[gain].plot([offset+noise, offset+noise],
                             [0, np.nanmax(hn)], linewidth=1.5, color='red')
             grid[gain].plot([offset, offset], [0, 0],
                             linewidth=1.5, color='y', label='median')
 
             grid[gain].plot([np.nanmedian(offset_g[cap][qm][:, :, 12, gain]),
                              np.nanmedian(offset_g[cap][qm][:, :, 12, gain])],
                             [0, np.nanmax(hn)], linewidth=1.5, color='green',
                             label='average over pixels')
 
             grid[gain].set_xlim(xrange)
             grid[gain].set_ylim(0, np.nanmax(hn)*1.1)
             grid[gain].set_xlabel("Offset value [ADU]")
             grid[gain].set_ylabel("# of occurance")
 
             if gain == 0:
                 leg = grid[gain].legend(
                     loc='upper center', ncol=3,
                     bbox_to_anchor=(0.1, 0.25, 0.7, 1.0))
 
             grid[gain].text(820, np.nanmax(hn)*0.4,
                             "{} gain".format(gain_names[gain]), fontsize=20)
 
             a = plt.axes([.125, .1, 0.775, .8], frame_on=False)
             a.patch.set_alpha(0.05)
             a.set_xlim(xrange)
             plt.plot([offset, offset], [0, 1], linewidth=1.5, color='y')
             plt.xticks([])
             plt.yticks([])
 
         ypos = 0.9
         x1pos = (np.nanmedian(data_g[cap][qm][:, 0]) +
                  np.nanmedian(data_g[cap][qm][:, 2]))/2.
         x2pos = (np.nanmedian(data_g[cap][qm][:, 2]) +
                  np.nanmedian(data_g[cap][qm][:, 1]))/2.-10
 
         plt.annotate("", xy=(np.nanmedian(data_g[cap][qm][:, 0]), ypos), xycoords='data',
                      xytext=(np.nanmedian(data_g[cap][qm][:, 2]), ypos), textcoords='data',
                      arrowprops=dict(arrowstyle="<->", connectionstyle="arc3"))
 
         plt.annotate('{}'.format(np.nanmedian(data_g[cap][qm][:, 0])-np.nanmedian(data_g[cap][qm][:, 2])),
                      xy=(x1pos, ypos), xycoords='data', xytext=(5, 5), textcoords='offset points')
 
         plt.annotate("", xy=(np.nanmedian(data_g[cap][qm][:, 2]), ypos), xycoords='data',
                      xytext=(np.nanmedian(data_g[cap][qm][:, 1]), ypos), textcoords='data',
                      arrowprops=dict(arrowstyle="<->", connectionstyle="arc3"))
 
         plt.annotate('{}'.format(np.nanmedian(data_g[cap][qm][:, 2])-np.nanmedian(data_g[cap][qm][:, 1])),
                      xy=(x2pos, ypos), xycoords='data', xytext=(5, 5), textcoords='offset points')
 
 plt.show()
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Normality test ##
 
 Distributions of raw pedestal values have been tested if they are normally distributed. A normality test have been performed for each pixel and each memory cell. Plots below show histogram of p-Values and a 2D distribution for the  memory cell 12.
 
 %% Cell type:code id: tags:
 
 ``` python
 # Loop over capacitor settings, modules, constants
 for cap in capacitor_settings:
     if not test_for_normality:
         print('Normality test was not requested. Flag `test_for_normality` False')
         break
     for i in modules:
         qm = module_index_to_qm(i)
 
         data = np.copy(ntest_g[cap][qm][:,:,:,:])
         data[badpix_g[cap][qm][:,:,:,:]>0] = 1.01
 
         hn,cn = np.histogram(data[:,:,:,0], bins=100)
 
         d = [{'x': np.arange(100)*0.01+0.01,
               'y': np.histogram(data[:,:,:,0], bins=100)[0],
               'drawstyle': 'steps-pre',
               'label' : 'High gain',
               },
              {'x': np.arange(100)*0.01+0.01,
               'y': np.histogram(data[:,:,:,1], bins=100)[0],
               'drawstyle': 'steps-pre',
               'label' : 'Medium gain',
               },
              {'x': np.arange(100)*0.01+0.01,
               'y': np.histogram(data[:,:,:,2], bins=100)[0],
               'drawstyle': 'steps-pre',
               'label' : 'Low gain',
               },
             ]
 
 
         fig = plt.figure(figsize=(15,15), tight_layout={'pad': 0.5, 'w_pad': 0.3})
 
         for gain in range(3):
             ax = fig.add_subplot(221+gain)
             heatmapPlot(data[:,:,12,gain], add_panels=False, cmap='viridis', figsize=(10,10),
                 y_label='Rows', x_label='Columns',
                 lut_label='p-Value',
                 use_axis=ax,
                 title='p-Value for cell 12, {} gain'.format(gain_names[gain]) )
 
         ax = fig.add_subplot(224)
         _ = simplePlot(d, #aspect=1.6,
                               x_label = "p-Value".format(gain),
                               y_label="# of occurance",
                               use_axis=ax,
                                y_log=False, legend='outside-top-ncol3-frame', legend_pad=0.05, legend_size='5%')
         ax.ticklabel_format(style='sci', axis='y', scilimits=(4,6))
 
 ```
 
 %% Cell type:raw id: tags:
 
 .. raw:: latex
 
     \newpage
 
 %% Cell type:markdown id: tags:
 
 ## Single-Cell Overviews ##
 
 Single cell overviews allow to identify potential effects on all memory cells, e.g. on a sensor level. Additionally, they should serve as a first sanity check on expected behaviour, e.g. if structuring on the ASIC level is visible in the offsets, but otherwise no immediate artifacts are visible.
 
 %% Cell type:code id: tags:
 
 ``` python
 cell = 12
 for cap in capacitor_settings:
     for gain in range(3):
         display(
             Markdown('### Cell-12 overview - {} gain'.format(gain_names[gain])))
 
         fig = plt.figure(figsize=(18, 22) , tight_layout={'pad': 0.1, 'w_pad': 0.1})
         for qm in res[cap]:
             for iconst, const in enumerate(['Offset', 'Noise', 'BadPixelsDark']):
 
                 ax = fig.add_subplot(321+iconst)
 
                 data = res[cap][qm][const][:, :, 12, gain]
                 vmax = 1.5 * np.nanmedian(res[cap][qm][const][:, :, 12, gain])
                 title = const
                 label = '{} value [ADU]'.format(const)
                 title = '{} value'.format(const)
                 if const == 'BadPixelsDark':
                     vmax = 4
                     bpix_code = data.astype(np.float32)
                     bpix_code[bpix_code == 0] = np.nan
                     title = 'Bad pixel code'
                     label = title
 
                     cb_labels = ['1 {}'.format(BadPixels.NOISE_OUT_OF_THRESHOLD.name),
                                  '2 {}'.format(BadPixels.OFFSET_NOISE_EVAL_ERROR.name),
                                  '3 {}'.format(BadPixels.OFFSET_OUT_OF_THRESHOLD.name),
                                  '4 {}'.format('MIXED')]
 
                     heatmapPlot(bpix_code, add_panels=False, cmap='viridis',
                                 y_label='Rows', x_label='Columns',
                                 lut_label='', vmax=vmax,
                                 use_axis=ax, cb_ticklabels=cb_labels, cb_ticks = np.arange(4)+1,
                                 title='{}'.format(title))
                     del bpix_code
                 else:
 
                     heatmapPlot(data, add_panels=False, cmap='viridis',
                                 y_label='Rows', x_label='Columns',
                                 lut_label=label, vmax=vmax,
                                 use_axis=ax,
                                 title='{}'.format(title))
 
         for qm in res[cap]:
             for iconst, const in enumerate(['Offset', 'Noise']):
                 data = res[cap][qm][const]
                 dataBP = np.copy(data)
                 dataBP[res[cap][qm]['BadPixelsDark'] > 0] = -1
 
                 x_ranges = [[0, 1500], [0, 40]]
                 hn, cn = np.histogram(
                     data[:, :, :, gain], bins=100, range=x_ranges[iconst])
                 hnBP, cnBP = np.histogram(dataBP[:, :, :, gain], bins=cn)
 
                 d = [{'x': cn[:-1],
                       'y': hn,
                       'drawstyle': 'steps-pre',
                       'label': 'All data',
                       },
                      {'x': cnBP[:-1],
                       'y': hnBP,
                       'drawstyle': 'steps-pre',
                       'label': 'Bad pixels masked',
                       },
                      ]
 
                 ax = fig.add_subplot(325+iconst)
                 _ = simplePlot(d, figsize=(5, 7), aspect=1,
                                     x_label="{} value [ADU]".format(const),
                                     y_label="# of occurance",
                                     title='', legend_pad=0.1, legend_size='10%',
                                     use_axis=ax,
                                     y_log=True, legend='outside-top-2col-frame')
 
         plt.show()
 ```
 
 %% Cell type:raw id: tags:
 
 .. raw:: latex
 
     \newpage
 
 %% Cell type:markdown id: tags:
 
 
 %% Cell type:code id: tags:
 
 ``` python
 cols = {BadPixels.NOISE_OUT_OF_THRESHOLD.value: (BadPixels.NOISE_OUT_OF_THRESHOLD.name, '#FF000080'),
         BadPixels.OFFSET_NOISE_EVAL_ERROR.value: (BadPixels.OFFSET_NOISE_EVAL_ERROR.name, '#0000FF80'),
         BadPixels.OFFSET_OUT_OF_THRESHOLD.value: (BadPixels.OFFSET_OUT_OF_THRESHOLD.name, '#00FF0080'),
         BadPixels.OFFSET_OUT_OF_THRESHOLD.value | BadPixels.NOISE_OUT_OF_THRESHOLD.value: ('MIXED', '#DD00DD80')}
 
 if high_res_badpix_3d:
     display(Markdown("""
 
     ## Global Bad Pixel Behaviour ##
 
     The following plots shows the results of a bad pixel evaluation for all evaluated memory cells.
     Cells are stacked in the Z-dimension, while pixels values in x/y are re-binned with a factor of 2.
     This excludes single bad pixels present only in disconnected pixels.
     Hence, any bad pixels spanning at least 4 pixels in the x/y-plane, or across at least two memory cells are indicated.
     Colors encode the bad pixel type, or mixed type.
 
         """))
     # Switch rebin to 1 for full resolution and
     # no interpolation for badpixel values.
     rebin = 2
     for gain in range(3):
         display(Markdown('### Bad pixel behaviour - {} gain ###'.format(gain_names[gain])))
         for cap in capacitor_settings:
             for mod, data in badpix_g[cap].items():
                 plot_badpix_3d(data[...,gain], cols, title='', rebin_fac=rebin)
                 ax = plt.gca()
                 leg = ax.get_legend()
                 leg.set(alpha=0.5)
         plt.show()
 ```
 
 %% Cell type:raw id: tags:
 
 .. raw:: latex
 
     \newpage
 
 %% Cell type:markdown id: tags:
 
 ## Summary across tiles ##
 
 Plots give an overview of calibration constants averaged across tiles. A bad pixel mask is applied. Constants are compared with pre-existing constants retrieved from the calibration database. Differences $\Delta$ between the old and new constants is shown.
 
 %% Cell type:code id: tags:
 
 ``` python
 time_summary = []
 for cap, cap_data in old_mdata.items():
     time_summary.append(f"The following pre-existing constants are used for comparison for capacitor setting **{cap}**:")
     for qm, qm_data in cap_data.items():
         time_summary.append(f"- Module {qm}")
         for const, const_data in qm_data.items():
             time_summary.append(f"    - {const} created at {const_data['timestamp']}")
 display(Markdown("\n".join(time_summary)))
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # Loop over capacitor settings, modules, constants
 for cap in res:
     for qm in res[cap]:
         for gain in range(3):
             display(Markdown('### Summary across tiles - {} gain'.format(gain_names[gain])))
 
             for const in res[cap][qm]:
                 data = np.copy(res[cap][qm][const][:, :, :, gain])
 
                 label = 'Fraction of bad pixels'
 
                 if const != 'BadPixelsDark':
                     data[badpix_g[cap][qm][:, :, :, gain] > 0] = np.nan
                     label = '{} value [ADU]'.format(const)
                 else:
                     data[data>0] = 1.0
 
                 data = data.reshape(
                     int(data.shape[0] / 32),
                     32,
                     int(data.shape[1] / 128),
                     128,
                     data.shape[2])
                 data = np.nanmean(data, axis=(1, 3)).swapaxes(
                     0, 2).reshape(512, 16)
 
                 fig = plt.figure(figsize=(15, 6))
                 ax = fig.add_subplot(121)
 
                 _ = heatmapPlot(data[:510, :], add_panels=True,
                                 y_label='Momery Cell ID', x_label='Tile ID',
                                 lut_label=label, use_axis=ax,
                                 panel_y_label=label, panel_x_label=label,
                                 cmap='viridis',  # cb_loc='right',cb_aspect=15,
                                 x_ticklabels=np.arange(16)+1,
                                 x_ticks=np.arange(16)+0.5)
 
                 if old_const[cap][qm][const] is not None:
                     ax = fig.add_subplot(122)
 
                     dataold = np.copy(old_const[cap][qm][const][:, :, :, gain])
 
                     label = '$\Delta$ {}'.format(label)
 
                     if const != 'BadPixelsDark':
                         if old_const[cap][qm]['BadPixelsDark'] is not None:
                             dataold[old_const[cap][qm]['BadPixelsDark'][:, :, :, gain] > 0] = np.nan
                         else:
                             dataold[:] = np.nan
                     else:
                         dataold[dataold>0]=1.0
 
                     dataold = dataold.reshape(
                         int(dataold.shape[0] / 32),
                         32,
                         int(dataold.shape[1] / 128),
                         128,
                         dataold.shape[2])
                     dataold = np.nanmean(dataold, axis=(
                         1, 3)).swapaxes(0, 2).reshape(512, 16)
                     dataold = dataold - data
 
                     _ = heatmapPlot(dataold[:510, :], add_panels=True,
                                     y_label='Momery Cell ID', x_label='Tile ID',
                                     lut_label=label, use_axis=ax,
                                     panel_y_label=label, panel_x_label=label,
                                     cmap='viridis',  # cb_loc='right',cb_aspect=15,
                                     x_ticklabels=np.arange(16)+1,
                                     x_ticks=np.arange(16)+0.5)
             plt.show()
 ```
 
 %% Cell type:raw id: tags:
 
 .. raw:: latex
 
     \newpage
 
 %% Cell type:markdown id: tags:
 
 ## Variation of offset and noise across Tiles and ASICs ##
 
 The following plots show a standard deviation $\sigma$ of the calibration constant. The plot of standard deviations across tiles show pixels of one tile ($128 \times 32$). Value for each pixel shows a standard deviation across 16 tiles. The standard deviation across ASICs are shown overall tiles. The plot shows pixels of one ASIC ($16 \times 32$), where the value shows a standard deviation across all ACIS of the module.
 
 %% Cell type:code id: tags:
 
 ``` python
 # Loop over capacitor settings, modules, constants
 for cap in res:
     for qm in res[cap]:
         for gain in range(3):
             display(Markdown('### Variation of offset and noise across ASICs - {} gain'.format(gain_names[gain])))
 
             fig = plt.figure(figsize=(15, 6))
             for iconst, const in enumerate(['Offset', 'Noise']):
                 data = np.copy(res[cap][qm][const][:, :, :, gain])
                 data[badpix_g[cap][qm][:, :, :, gain] > 0] = np.nan
                 label = '$\sigma$ {} [ADU]'.format(const)
 
                 dataA = np.nanmean(data, axis=2)  # average over cells
                 dataA = dataA.reshape(8, 32, 16, 16)
                 dataA = np.nanstd(dataA, axis=(0, 2))  # average across ASICs
 
                 ax = fig.add_subplot(121+iconst)
                 _ = heatmapPlot(dataA, add_panels=True,
                                 y_label='rows', x_label='columns',
                                 lut_label=label, use_axis=ax,
                                 panel_y_label=label, panel_x_label=label,
                                 cmap='viridis'
                                 )
 
             plt.show()
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # Loop over capacitor settings, modules, constants
 for cap in res:
     for qm in res[cap]:
         for gain in range(3):
             display(Markdown('### Variation of offset and noise across tiles - {} gain'.format(gain_names[gain])))
 
             fig = plt.figure(figsize=(15, 6))
             for iconst, const in enumerate(['Offset', 'Noise']):
                 data = np.copy(res[cap][qm][const][:, :, :, gain])
                 data[badpix_g[cap][qm][:, :, :, gain] > 0] = np.nan
                 label = '$\sigma$ {} [ADU]'.format(const)
 
                 dataT = data.reshape(
                     int(data.shape[0] / 32),
                     32,
                     int(data.shape[1] / 128),
                     128,
                     data.shape[2])
                 dataT = np.nanstd(dataT, axis=(0, 2))
                 dataT = np.nanmean(dataT, axis=2)
 
                 ax = fig.add_subplot(121+iconst)
                 _ = heatmapPlot(dataT, add_panels=True,
                                 y_label='rows', x_label='columns',
                                 lut_label=label, use_axis=ax,
                                 panel_y_label=label, panel_x_label=label,
                                 cmap='viridis')
             plt.show()
 ```
 
 %% Cell type:raw id: tags:
 
 .. raw:: latex
 
     \newpage
 
 %% Cell type:markdown id: tags:
 
 ## Aggregate values and per cell behaviour ##
 
 The following tables and plots give an overview of statistical aggregates for each constant, as well as per-cell behavior, averaged across pixels.
 
 %% Cell type:code id: tags:
 
 ``` python
 # Loop over capacitor settings, modules, constants
 for cap in res:
     for qm in res[cap]:
         for gain in range(3):
             display(Markdown('### Mean over pixels - {} gain'.format(gain_names[gain])))
 
             fig = plt.figure(figsize=(9,11))
 
             for iconst, const in enumerate(res[cap][qm]):
 
                 ax = fig.add_subplot(311+iconst)
 
                 data = res[cap][qm][const][:,:,:510,gain]
                 if const == 'BadPixelsDark':
                     data[data>0] = 1.0
 
                 dataBP = np.copy(data)
                 dataBP[badpix_g[cap][qm][:,:,:510,gain]>0] = -10
 
                 data = np.nanmean(data, axis=(0,1))
                 dataBP = np.nanmean(dataBP, axis=(0,1))
 
                 d = [{'y': data,
                       'x': np.arange(data.shape[0]),
                       'drawstyle': 'steps-mid',
                       'label' : 'All data'
                      }
                     ]
 
                 if const != 'BadPixelsDark':
                     d.append({'y': dataBP,
                       'x': np.arange(data.shape[0]),
                       'drawstyle': 'steps-mid',
                       'label' : 'good pixels only'
                      })
                     y_title = "{} value [ADU]".format(const)
                     title = "{} value, {} gain".format(const, gain_names[gain])
                 else:
                     y_title = "Fraction of Bad Pixels"
                     title = "Fraction of Bad Pixels, {} gain".format(gain_names[gain])
 
                 data_min = np.min([data, dataBP])if const != 'BadPixelsDark' else np.min([data])
                 data_max = np.max([data[20:], dataBP[20:]])
                 data_dif = data_max - data_min
 
                 local_max = np.max([data[200:300], dataBP[200:300]])
                 frac = 0.35
                 new_max = (local_max - data_min*(1-frac))/frac
                 new_max = np.max([data_max, new_max])
 
                 _ = simplePlot(d, figsize=(10,10), aspect=2, xrange=(-12, 510),
                                   x_label = 'Memory Cell ID',
                                   y_label=y_title, use_axis=ax,
                                   title=title,
                                   title_position=[0.5, 1.15],
                                   inset='xy-coord-right', inset_x_range=(0,20), inset_indicated=True,
                                   inset_labeled=True, inset_coord=[0.2,0.5,0.6,0.95],
                                     inset_lw = 1.0, y_range = [data_min-data_dif*0.05, new_max+data_dif*0.05],
                                   y_log=False, legend='outside-top-ncol2-frame', legend_size='18%',
                                      legend_pad=0.00)
 
                 plt.tight_layout(pad=1.08, h_pad=0.35)
 
             plt.show()
 ```
 
 %% Cell type:raw id: tags:
 
 .. raw:: latex
 
     \newpage
 
 %% Cell type:markdown id: tags:
 
 ## Summary tables ##
 
 The following tables show summary information for the evaluated module. Values for currently evaluated constants are compared with values for pre-existing constants retrieved from the calibration database.
 
 %% Cell type:code id: tags:
 
 ``` python
 table = []
 bits = [BadPixels.NOISE_OUT_OF_THRESHOLD, BadPixels.OFFSET_OUT_OF_THRESHOLD, BadPixels.OFFSET_NOISE_EVAL_ERROR]
 for cap in res:
     for qm in res[cap]:
         for gain in range(3):
 
             l_data = []
             l_data_old = []
 
             data = np.copy(res[cap][qm]['BadPixelsDark'][:,:,:,gain])
             l_data.append(len(data[data>0].flatten()))
             for bit in bits:
                 l_data.append(np.count_nonzero(badpix_g[cap][qm][:,:,:,gain] & bit.value))
 
             if old_const[cap][qm]['BadPixelsDark'] is not None:
                 old_const[cap][qm]['BadPixelsDark'] = old_const[cap][qm]['BadPixelsDark'].astype(np.uint32)
                 dataold = np.copy(old_const[cap][qm]['BadPixelsDark'][:, :, :, gain])
                 l_data_old.append(len(dataold[dataold>0].flatten()))
                 for bit in bits:
                     l_data_old.append(np.count_nonzero(old_const[cap][qm]['BadPixelsDark'][:, :, :, gain] & bit.value))
 
             l_data_name = ['All bad pixels', 'NOISE_OUT_OF_THRESHOLD',
                            'OFFSET_OUT_OF_THRESHOLD', 'OFFSET_NOISE_EVAL_ERROR']
 
             l_threshold = ['', f'{thresholds_noise_sigma}', f'{thresholds_offset_sigma}',
                            f'{thresholds_offset_hard}/{thresholds_noise_hard}']
 
             for i in range(len(l_data)):
                 line = [f'{l_data_name[i]}, gain {gain_names[gain]}', l_threshold[i], l_data[i]]
 
                 if old_const[cap][qm]['BadPixelsDark'] is not None:
                     line += [l_data_old[i]]
                 else:
                     line += ['-']
 
                 table.append(line)
             table.append(['', '', '', ''])
 
 display(Markdown('''
 
 ### Number of bad pixels ###
 
 One pixel can be bad for different reasons, therefore, the sum of all types of bad pixels can be more than the number of all bad pixels.
 
 '''))
 if len(table)>0:
     md = display(Latex(tabulate.tabulate(table, tablefmt='latex',
                                      headers=["Pixel type", "Threshold",
                                               "New constant", "Old constant"])))
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 header = ['Parameter',
           "New constant", "Old constant ",
           "New constant", "Old constant ",
           "New constant", "Old constant "]
 
 for const in ['Offset', 'Noise']:
     table = [['','High gain', 'High gain', 'Medium gain', 'Medium gain', 'Low gain', 'Low gain']]
     for cap in res:
         for qm in res[cap]:
 
             data = np.copy(res[cap][qm][const])
             data[res[cap][qm]['BadPixelsDark']>0] = np.nan
 
             if old_const[cap][qm][const] is not None and old_const[cap][qm]['BadPixelsDark'] is not None :
                 dataold = np.copy(old_const[cap][qm][const])
                 dataold[old_const[cap][qm]['BadPixelsDark']>0] = np.nan
 
             f_list = [np.nanmedian, np.nanmean, np.nanstd, np.nanmin, np.nanmax]
             n_list = ['Median', 'Mean', 'Std', 'Min', 'Max']
 
             for i, f in enumerate(f_list):
                 line = [n_list[i]]
                 for gain in range(3):
                     line.append('{:6.1f}'.format(f(data[...,gain])))
                     if old_const[cap][qm][const] is not None and old_const[cap][qm]['BadPixelsDark'] is not None:
                         line.append('{:6.1f}'.format(f(dataold[...,gain])))
                     else:
                         line.append('-')
 
                 table.append(line)
 
     display(Markdown('### {} [ADU], good pixels only ###'.format(const)))
     md = display(Latex(tabulate.tabulate(table, tablefmt='latex', headers=header)))
 ```

notebooks/LPD/LPD_Correct_Fast.ipynb

 %% Cell type:markdown id: tags:
 
 # LPD Offline Correction #
 
 Author: European XFEL Data Analysis Group
 
 %% Cell type:code id: tags:
 
 ``` python
 # Input parameters
 in_folder = "/gpfs/exfel/exp/FXE/202201/p003073/raw/"  # the folder to read data from, required
 out_folder = "/gpfs/exfel/data/scratch/schmidtp/random/LPD_test"  # the folder to output to, required
 metadata_folder = ''  # Directory containing calibration_metadata.yml when run by xfel-calibrate.
 sequences = [-1]  # Sequences to correct, use [-1] for all
 modules = [-1]  # Modules indices to correct, use [-1] for all, only used when karabo_da is empty
 karabo_da = ['']  # Data aggregators names to correct, use [''] for all
 run = 10  # run to process, required
 
 # Source parameters
 karabo_id = 'FXE_DET_LPD1M-1'  # Karabo domain for detector.
 input_source = '{karabo_id}/DET/{module_index}CH0:xtdf'  # Input fast data source.
 output_source = ''  # Output fast data source, empty to use same as input.
 
 # CalCat parameters
 creation_time = ""  # The timestamp to use with Calibration DB. Required Format: "YYYY-MM-DD hh:mm:ss" e.g. 2019-07-04 11:02:41
 cal_db_interface = ''  # Not needed, compatibility with current webservice.
 cal_db_timeout = 0  # Not needed, compatbility with current webservice.
 cal_db_root = '/gpfs/exfel/d/cal/caldb_store'
 
 # Operating conditions
 mem_cells = 512  # Memory cells, LPD constants are always taken with 512 cells.
 bias_voltage = 250.0  # Detector bias voltage.
 capacitor = '5pF'  # Capacitor setting: 5pF or 50pF
 photon_energy = 9.2  # Photon energy in keV.
 category = 0  # Whom to blame.
+use_cell_order = False  # Whether to use memory cell order as a detector condition (not stored for older constants)
 
 # Correction parameters
 offset_corr = True  # Offset correction.
 rel_gain = True  # Gain correction based on RelativeGain constant.
 ff_map = True  # Gain correction based on FFMap constant.
 gain_amp_map = True  # Gain correction based on GainAmpMap constant.
 
 # Output options
 overwrite = True  # set to True if existing data should be overwritten
 chunks_data = 1  # HDF chunk size for pixel data in number of frames.
 chunks_ids = 32  # HDF chunk size for cellId and pulseId datasets.
 create_virtual_cxi_in = ''  # Folder to create virtual CXI files in (for each sequence).
 
 # Parallelization options
 sequences_per_node = 1  # Sequence files to process per node
 max_nodes = 8  # Maximum number of SLURM jobs to split correction work into
 num_workers = 8  # Worker processes per node, 8 is safe on 768G nodes but won't work on 512G.
 num_threads_per_worker = 32  # Number of threads per worker.
 
 def balance_sequences(in_folder, run, sequences, sequences_per_node, karabo_da, max_nodes):
     from xfel_calibrate.calibrate import balance_sequences as bs
     return bs(in_folder, run, sequences, sequences_per_node, karabo_da, max_nodes=max_nodes)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 from collections import OrderedDict
 from pathlib import Path
 from time import perf_counter
 import gc
 import re
 import warnings
 
 import numpy as np
 import h5py
 
 import matplotlib
 matplotlib.use('agg')
 import matplotlib.pyplot as plt
 %matplotlib inline
 
 from calibration_client import CalibrationClient
 from calibration_client.modules import CalibrationConstantVersion
 import extra_data as xd
 import extra_geom as xg
 import pasha as psh
 
 from extra_data.components import LPD1M
 
 from cal_tools.lpdalgs import correct_lpd_frames
+from cal_tools.lpdlib import get_mem_cell_order
 from cal_tools.tools import CalibrationMetadata, calcat_creation_time
 from cal_tools.files import DataFile
 from cal_tools.restful_config import restful_config
 ```
 
 %% Cell type:markdown id: tags:
 
 # Prepare environment
 
 %% Cell type:code id: tags:
 
 ``` python
 file_re = re.compile(r'^RAW-R(\d{4})-(\w+\d+)-S(\d{5})$')  # This should probably move to cal_tools
 
 run_folder = Path(in_folder) / f'r{run:04d}'
 out_folder = Path(out_folder)
 out_folder.mkdir(exist_ok=True)
 
 output_source = output_source or input_source
 
 cal_db_root = Path(cal_db_root)
 
 metadata = CalibrationMetadata(metadata_folder or out_folder)
 
 creation_time = calcat_creation_time(in_folder, run, creation_time)
 print(f'Using {creation_time.isoformat()} as creation time')
 
 # Pick all modules/aggregators or those selected.
 if not karabo_da or karabo_da == ['']:
     if not modules or modules == [-1]:
         modules = list(range(16))
 
     karabo_da = [f'LPD{i:02d}' for i in modules]
 
 # Pick all sequences or those selected.
 if not sequences or sequences == [-1]:
     do_sequence = lambda seq: True
 else:
     do_sequence = [int(x) for x in sequences].__contains__
 
 # List of detector sources.
 det_inp_sources = [input_source.format(karabo_id=karabo_id, module_index=int(da[-2:])) for da in karabo_da]
 ```
 
 %% Cell type:markdown id: tags:
 
 # Select data to process
 
 %% Cell type:code id: tags:
 
 ``` python
 data_to_process = []
 
 for inp_path in run_folder.glob('RAW-*.h5'):
     match = file_re.match(inp_path.stem)
 
     if match[2] not in karabo_da or not do_sequence(int(match[3])):
         continue
 
     outp_path = out_folder / 'CORR-R{run:04d}-{aggregator}-S{seq:05d}.h5'.format(
         run=int(match[1]), aggregator=match[2], seq=int(match[3]))
 
     data_to_process.append((match[2], inp_path, outp_path))
 
 print('Files to process:')
 for data_descr in sorted(data_to_process, key=lambda x: f'{x[0]}{x[1]}'):
     print(f'{data_descr[0]}\t{data_descr[1]}')
 ```
 
 %% Cell type:markdown id: tags:
 
 # Obtain and prepare calibration constants
 
 %% Cell type:code id: tags:
 
 ``` python
 # Connect to CalCat.
 calcat_config = restful_config['calcat']
 client = CalibrationClient(
     base_api_url=calcat_config['base-api-url'],
     use_oauth2=calcat_config['use-oauth2'],
     client_id=calcat_config['user-id'],
     client_secret=calcat_config['user-secret'],
     user_email=calcat_config['user-email'],
     token_url=calcat_config['token-url'],
     refresh_url=calcat_config['refresh-url'],
     auth_url=calcat_config['auth-url'],
     scope='')
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 metadata = CalibrationMetadata(metadata_folder or out_folder)
 # Constant paths & timestamps are saved under retrieved-constants in calibration_metadata.yml
 const_yaml = metadata.setdefault("retrieved-constants", {})
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 const_data = {}
 const_load_mp = psh.ProcessContext(num_workers=24)
 
 if const_yaml:  # Read constants from YAML file.
     start = perf_counter()
     for da, ccvs in const_yaml.items():
 
         for calibration_name, ccv in ccvs['constants'].items():
+            if ccv['file-path'] is None:
+                warnings.warn(f"Missing {calibration_name} for {da}")
+                continue
 
             dtype = np.uint32 if calibration_name.startswith('BadPixels') else np.float32
 
             const_data[(da, calibration_name)] = dict(
                 path=Path(ccv['file-path']),
                 dataset=ccv['dataset-name'],
                 data=const_load_mp.alloc(shape=(256, 256, mem_cells, 3), dtype=dtype)
             )
 else:  # Retrieve constants from CALCAT.
     dark_calibrations = {
         1: 'Offset',  # np.float32
         14: 'BadPixelsDark'  # should be np.uint32, but is np.float64
     }
 
-    dark_condition = [
+    base_condition = [
         dict(parameter_id=1, value=bias_voltage),  # Sensor bias voltage
         dict(parameter_id=7, value=mem_cells),  # Memory cells
         dict(parameter_id=15, value=capacitor),  # Feedback capacitor
         dict(parameter_id=13, value=256),  # Pixels X
         dict(parameter_id=14, value=256),  # Pixels Y
     ]
+    if use_cell_order:
+        # Read the order of memory cells used
+        raw_data = xd.DataCollection.from_paths([e[1] for e in data_to_process])
+        cell_ids_pattern_s = get_mem_cell_order(raw_data, det_inp_sources)
+        print("Memory cells order:", cell_ids_pattern_s)
+
+        dark_condition = base_condition + [
+            dict(parameter_id=30, value=cell_ids_pattern_s),  # Memory cell order
+        ]
+    else:
+        dark_condition = base_condition.copy()
 
     illuminated_calibrations = {
         20: 'BadPixelsFF',  # np.uint32
         42: 'GainAmpMap',  # np.float32
         43: 'FFMap',  # np.float32
         44: 'RelativeGain'  # np.float32
     }
 
-    illuminated_condition = dark_condition.copy()
-    illuminated_condition += [
+    illuminated_condition = base_condition + [
         dict(parameter_id=3, value=photon_energy),  # Source energy
         dict(parameter_id=25, value=category)  # category
     ]
 
     print('Querying calibration database', end='', flush=True)
     start = perf_counter()
     for calibrations, condition in [
         (dark_calibrations, dark_condition),
         (illuminated_calibrations, illuminated_condition)
     ]:
         resp = CalibrationConstantVersion.get_closest_by_time_by_detector_conditions(
             client, karabo_id, list(calibrations.keys()),
             {'parameters_conditions_attributes': condition},
             karabo_da='', event_at=creation_time.isoformat()
         )
 
         if not resp['success']:
             raise RuntimeError(resp)
 
         for ccv in resp['data']:
             cc = ccv['calibration_constant']
             da = ccv['physical_detector_unit']['karabo_da']
             calibration_name = calibrations[cc['calibration_id']]
 
             dtype = np.uint32 if calibration_name.startswith('BadPixels') else np.float32
 
             const_data[(da, calibration_name)] = dict(
                 path=Path(ccv['path_to_file']) / ccv['file_name'],
                 dataset=ccv['data_set_name'],
                 data=const_load_mp.alloc(shape=(256, 256, mem_cells, 3), dtype=dtype)
             )
         print('.', end='', flush=True)
 
 total_time = perf_counter() - start
 print(f'{total_time:.1f}s')
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 def load_constant_dataset(wid, index, const_descr):
     ccv_entry = const_data[const_descr]
 
     with h5py.File(cal_db_root / ccv_entry['path'], 'r') as fp:
         fp[ccv_entry['dataset'] + '/data'].read_direct(ccv_entry['data'])
 
     print('.', end='', flush=True)
 
 print('Loading calibration data', end='', flush=True)
 start = perf_counter()
 const_load_mp.map(load_constant_dataset, list(const_data.keys()))
 total_time = perf_counter() - start
 
 print(f'{total_time:.1f}s')
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # These are intended in order cell, X, Y, gain
 ccv_offsets = {}
 ccv_gains = {}
 ccv_masks = {}
 
 ccv_shape = (mem_cells, 256, 256, 3)
 
 constant_order = {
     'Offset':        (2, 1, 0, 3),
     'BadPixelsDark': (2, 1, 0, 3),
     'RelativeGain':  (2, 1, 0, 3),
     'FFMap':         (2, 0, 1, 3),
     'BadPixelsFF':   (2, 0, 1, 3),
     'GainAmpMap':    (2, 0, 1, 3),
 }
 
 def prepare_constants(wid, index, aggregator):
     consts = {calibration_name: entry['data']
               for (aggregator_, calibration_name), entry
               in const_data.items()
               if aggregator == aggregator_}
 
     def _prepare_data(calibration_name, dtype):
         return consts[calibration_name] \
             .transpose(constant_order[calibration_name]) \
             .astype(dtype, copy=True)  # Make sure array is contiguous.
 
     if offset_corr and 'Offset' in consts:
         ccv_offsets[aggregator] = _prepare_data('Offset', np.float32)
     else:
         ccv_offsets[aggregator] = np.zeros(ccv_shape, dtype=np.float32)
 
     ccv_gains[aggregator] = np.ones(ccv_shape, dtype=np.float32)
 
     if 'BadPixelsDark' in consts:
         ccv_masks[aggregator] = _prepare_data('BadPixelsDark', np.uint32)
     else:
         ccv_masks[aggregator] = np.zeros(ccv_shape, dtype=np.uint32)
 
     if rel_gain and 'RelativeGain' in consts:
         ccv_gains[aggregator] *= _prepare_data('RelativeGain', np.float32)
 
     if ff_map and 'FFMap' in consts:
         ccv_gains[aggregator] *= _prepare_data('FFMap', np.float32)
 
         if 'BadPixelsFF' in consts:
             np.bitwise_or(ccv_masks[aggregator], _prepare_data('BadPixelsFF', np.uint32),
                           out=ccv_masks[aggregator])
 
     if gain_amp_map and 'GainAmpMap' in consts:
         ccv_gains[aggregator] *= _prepare_data('GainAmpMap', np.float32)
 
     print('.', end='', flush=True)
 
 
 print('Preparing constants', end='', flush=True)
 start = perf_counter()
 psh.ThreadContext(num_workers=len(karabo_da)).map(prepare_constants, karabo_da)
 total_time = perf_counter() - start
 print(f'{total_time:.1f}s')
 
 const_data.clear()  # Clear raw constants data now to save memory.
 gc.collect();
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 def correct_file(wid, index, work):
     aggregator, inp_path, outp_path = work
     module_index = int(aggregator[-2:])
 
     start = perf_counter()
     dc = xd.H5File(inp_path, inc_suspect_trains=False).select('*', 'image.*', require_all=True)
     inp_source = dc[input_source.format(karabo_id=karabo_id, module_index=module_index)]
     open_time = perf_counter() - start
 
     # Load raw data for this file.
     # Reshaping gets rid of the extra 1-len dimensions without
     # mangling the frame axis for an actual frame count of 1.
     start = perf_counter()
     in_raw = inp_source['image.data'].ndarray().reshape(-1, 256, 256)
     in_cell = inp_source['image.cellId'].ndarray().reshape(-1)
     in_pulse = inp_source['image.pulseId'].ndarray().reshape(-1)
     read_time = perf_counter() - start
 
     # Allocate output arrays.
     out_data = np.zeros((in_raw.shape[0], 256, 256), dtype=np.float32)
     out_gain = np.zeros((in_raw.shape[0], 256, 256), dtype=np.uint8)
     out_mask = np.zeros((in_raw.shape[0], 256, 256), dtype=np.uint32)
 
     start = perf_counter()
     correct_lpd_frames(in_raw, in_cell,
                        out_data, out_gain, out_mask,
                        ccv_offsets[aggregator], ccv_gains[aggregator], ccv_masks[aggregator],
                        num_threads=num_threads_per_worker)
     correct_time = perf_counter() - start
 
     image_counts = inp_source['image.data'].data_counts(labelled=False)
 
     start = perf_counter()
     if (not outp_path.exists() or overwrite) and image_counts.sum() > 0:
         outp_source_name = output_source.format(karabo_id=karabo_id, module_index=module_index)
 
         with DataFile(outp_path, 'w') as outp_file:
             outp_file.create_index(dc.train_ids, from_file=dc.files[0])
             outp_file.create_metadata(like=dc, instrument_channels=(f'{outp_source_name}/image',))
 
             outp_source = outp_file.create_instrument_source(outp_source_name)
 
             outp_source.create_index(image=image_counts)
             outp_source.create_key('image.cellId', data=in_cell,
                                    chunks=(min(chunks_ids, in_cell.shape[0]),))
             outp_source.create_key('image.pulseId', data=in_pulse,
                                    chunks=(min(chunks_ids, in_pulse.shape[0]),))
             outp_source.create_key('image.data', data=out_data,
                                    chunks=(min(chunks_data, out_data.shape[0]), 256, 256))
             outp_source.create_compressed_key('image.gain', data=out_gain)
             outp_source.create_compressed_key('image.mask', data=out_mask)
     write_time = perf_counter() - start
 
     total_time = open_time + read_time + correct_time + write_time
     frame_rate = in_raw.shape[0] / total_time
 
     print('{}\t{}\t{:.3f}\t{:.3f}\t{:.3f}\t{:.3f}\t{:.3f}\t{}\t{:.1f}'.format(
         wid, aggregator, open_time, read_time, correct_time, write_time, total_time,
         in_raw.shape[0], frame_rate))
 
     in_raw = None
     in_cell = None
     in_pulse = None
     out_data = None
     out_gain = None
     out_mask = None
     gc.collect()
 
 print('worker\tDA\topen\tread\tcorrect\twrite\ttotal\tframes\trate')
 start = perf_counter()
 psh.ProcessContext(num_workers=num_workers).map(correct_file, data_to_process)
 total_time = perf_counter() - start
 print(f'Total time: {total_time:.1f}s')
 ```
 
 %% Cell type:markdown id: tags:
 
 # Data preview for first train
 
 %% Cell type:code id: tags:
 
 ``` python
 geom = xg.LPD_1MGeometry.from_quad_positions(
     [(11.4, 299), (-11.5, 8), (254.5, -16), (278.5, 275)])
 
 output_paths = [outp_path for _, _, outp_path in data_to_process if outp_path.exists()]
 dc = xd.DataCollection.from_paths(output_paths).select_trains(np.s_[0])
 
 det = LPD1M(dc, detector_name=karabo_id)
 data = det.get_array('image.data')
 ```
 
 %% Cell type:markdown id: tags:
 
 ### Intensity histogram across all cells
 
 %% Cell type:code id: tags:
 
 ``` python
 left_edge_ratio = 0.01
 right_edge_ratio = 0.99
 
 fig, ax = plt.subplots(num=1, clear=True, figsize=(15, 6))
 values, bins, _ = ax.hist(np.ravel(data.data), bins=2000, range=(-1500, 2000))
 
 def find_nearest_index(array, value):
     return (np.abs(array - value)).argmin()
 
 cum_values = np.cumsum(values)
 vmin = bins[find_nearest_index(cum_values, cum_values[-1]*left_edge_ratio)]
 vmax = bins[find_nearest_index(cum_values, cum_values[-1]*right_edge_ratio)]
 
 max_value = values.max()
 ax.vlines([vmin, vmax], 0, max_value, color='red', linewidth=5, alpha=0.2)
 ax.text(vmin, max_value, f'{left_edge_ratio*100:.0f}%',
         color='red', ha='center', va='bottom', size='large')
 ax.text(vmax, max_value, f'{right_edge_ratio*100:.0f}%',
         color='red', ha='center', va='bottom', size='large')
 ax.text(vmax+(vmax-vmin)*0.01, max_value/2, 'Colormap interval',
         color='red', rotation=90, ha='left', va='center', size='x-large')
 
 ax.set_xlim(vmin-(vmax-vmin)*0.1, vmax+(vmax-vmin)*0.1)
 ax.set_ylim(0, max_value*1.1)
 pass
 ```
 
 %% Cell type:markdown id: tags:
 
 ### First memory cell
 
 %% Cell type:code id: tags:
 
 ``` python
 fig, ax = plt.subplots(num=2, figsize=(15, 15), clear=True, nrows=1, ncols=1)
 geom.plot_data_fast(data[:, 0, 0], ax=ax, vmin=vmin, vmax=vmax)
 pass
 ```
 
 %% Cell type:markdown id: tags:
 
 ### Train average
 
 %% Cell type:code id: tags:
 
 ``` python
 fig, ax = plt.subplots(num=3, figsize=(15, 15), clear=True, nrows=1, ncols=1)
 geom.plot_data_fast(data[:, 0].mean(axis=1), ax=ax, vmin=vmin, vmax=vmax)
 pass
 ```
 
 %% Cell type:markdown id: tags:
 
 ### Lowest gain stage per pixel
 
 %% Cell type:code id: tags:
 
 ``` python
 highest_gain_stage = det.get_array('image.gain', pulses=np.s_[:]).max(axis=(1, 2))
 
 fig, ax = plt.subplots(num=4, figsize=(15, 15), clear=True, nrows=1, ncols=1)
 p = geom.plot_data_fast(highest_gain_stage, ax=ax, vmin=0, vmax=2);
 
 cb = ax.images[0].colorbar
 cb.set_ticks([0, 1, 2])
 cb.set_ticklabels(['High gain', 'Medium gain', 'Low gain'])
 ```
 
 %% Cell type:markdown id: tags:
 
 ### Create virtual CXI file
 
 %% Cell type:code id: tags:
 
 ``` python
 if create_virtual_cxi_in:
     vcxi_folder = Path(create_virtual_cxi_in.format(
         run=run, proposal_folder=str(Path(in_folder).parent)))
     vcxi_folder.mkdir(parents=True, exist_ok=True)
 
     def sort_files_by_seq(by_seq, outp_path):
         by_seq.setdefault(int(outp_path.stem[-5:]), []).append(outp_path)
         return by_seq
 
     from functools import reduce
     reduce(sort_files_by_seq, output_paths, output_by_seq := {})
 
     for seq_number, seq_output_paths in output_by_seq.items():
         # Create data collection and detector components only for this sequence.
         try:
             det = LPD1M(xd.DataCollection.from_paths(seq_output_paths), detector_name=karabo_id, min_modules=4)
         except ValueError:  # Couldn't find enough data for min_modules
             continue
         det.write_virtual_cxi(vcxi_folder / f'VCXI-LPD-R{run:04d}-S{seq_number:05d}.cxi')
 ```

notebooks/LPD/LPD_retrieve_constants_precorrection.ipynb

 %% Cell type:markdown id: tags:
 
 # LPD Retrieving Constants Pre-correction #
 
 Author: European XFEL Detector Group, Version: 1.0
 
 The following notebook provides a constants metadata in a YAML file to use while correcting LPD images.
 
 %% Cell type:code id: tags:
 
 ``` python
 # Input parameters
 in_folder = "/gpfs/exfel/exp/FXE/202201/p003073/raw/"  # the folder to read data from, required
 out_folder = "/gpfs/exfel/data/scratch/ahmedk/test/remove/LPD_test"  # the folder to output to, required
 metadata_folder = ''  # Directory containing calibration_metadata.yml when run by xfel-calibrate.
 modules = [-1]  # Modules indices to correct, use [-1] for all, only used when karabo_da is empty
 karabo_da = ['']  # Data aggregators names to correct, use [''] for all
 run = 10  # run to process, required
 
 # Source parameters
 karabo_id = 'FXE_DET_LPD1M-1'  # Karabo domain for detector.
+input_source = '{karabo_id}/DET/{module_index}CH0:xtdf'  # Input fast data source.
 
 # CalCat parameters
 creation_time = ""  # The timestamp to use with Calibration DB. Required Format: "YYYY-MM-DD hh:mm:ss" e.g. 2019-07-04 11:02:41
 
 # Operating conditions
 mem_cells = 512  # Memory cells, LPD constants are always taken with 512 cells.
 bias_voltage = 250.0  # Detector bias voltage.
 capacitor = '5pF'  # Capacitor setting: 5pF or 50pF
 photon_energy = 9.2  # Photon energy in keV.
 category = 0  # Whom to blame.
+use_cell_order = False  # Whether to use memory cell order as a detector condition (not stored for older constants)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 from pathlib import Path
 from time import perf_counter
 
+import numpy as np
+
 from calibration_client import CalibrationClient
 from calibration_client.modules import CalibrationConstantVersion
+import extra_data as xd
 
+from cal_tools.lpdlib import get_mem_cell_order
 from cal_tools.tools import (
     CalibrationMetadata,
     calcat_creation_time,
     save_constant_metadata,
 )
 from cal_tools.restful_config import restful_config
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 out_folder = Path(out_folder)
 out_folder.mkdir(exist_ok=True)
 
 metadata = CalibrationMetadata(metadata_folder or out_folder)
 # Constant paths & timestamps are saved under retrieved-constants in calibration_metadata.yml
 retrieved_constants = metadata.setdefault("retrieved-constants", {})
 
 creation_time = calcat_creation_time(in_folder, run, creation_time)
 print(f'Using {creation_time.isoformat()} as creation time')
 
 # Pick all modules/aggregators or those selected.
 if not karabo_da or karabo_da == ['']:
     if not modules or modules == [-1]:
         modules = list(range(16))
 
     karabo_da = [f'LPD{i:02d}' for i in modules]
+
+# List of detector sources.
+det_inp_sources = [input_source.format(karabo_id=karabo_id, module_index=int(da[-2:])) for da in karabo_da]
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # Connect to CalCat.
 calcat_config = restful_config['calcat']
 client = CalibrationClient(
     base_api_url=calcat_config['base-api-url'],
     use_oauth2=calcat_config['use-oauth2'],
     client_id=calcat_config['user-id'],
     client_secret=calcat_config['user-secret'],
     user_email=calcat_config['user-email'],
     token_url=calcat_config['token-url'],
     refresh_url=calcat_config['refresh-url'],
     auth_url=calcat_config['auth-url'],
     scope='')
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 dark_calibrations = {
     1: 'Offset',
     14: 'BadPixelsDark',
 }
 
-dark_condition = [
+base_condition = [
     dict(parameter_id=1, value=bias_voltage),  # Sensor bias voltage
     dict(parameter_id=7, value=mem_cells),  # Memory cells
     dict(parameter_id=15, value=capacitor),  # Feedback capacitor
     dict(parameter_id=13, value=256),  # Pixels X
     dict(parameter_id=14, value=256),  # Pixels Y
 ]
+if use_cell_order:
+    # Read the order of memory cells used
+    raw_data = xd.RunDirectory(Path(in_folder, f'r{run:04d}'))
+    cell_ids_pattern_s = get_mem_cell_order(raw_data, det_inp_sources)
+    print("Memory cells order:", cell_ids_pattern_s)
+
+    dark_condition = base_condition + [
+        dict(parameter_id=30, value=cell_ids_pattern_s),  # Memory cell order
+    ]
+else:
+    dark_condition = base_condition.copy()
 
 illuminated_calibrations = {
     20: 'BadPixelsFF',
     42: 'GainAmpMap',
     43: 'FFMap',
     44: 'RelativeGain',
 }
 
-illuminated_condition = dark_condition.copy()
-illuminated_condition += [
+illuminated_condition = base_condition + [
     dict(parameter_id=3, value=photon_energy),  # Source energy
     dict(parameter_id=25, value=category)  # category
 ]
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 const_data = {}
 
 print('Querying calibration database', end='', flush=True)
 start = perf_counter()
 for k_da in karabo_da:
     pdu = None
-    if k_da in retrieved_constants:
-        print(f"Constant for {k_da} already in {metadata.filename}, won't query again.")  # noqa
-        continue
     retrieved_constants[k_da] = dict()
     const_mdata = retrieved_constants[k_da]["constants"] = dict()
     for calibrations, condition in [
         (dark_calibrations, dark_condition),
         (illuminated_calibrations, illuminated_condition)
     ]:
         resp = CalibrationConstantVersion.get_closest_by_time_by_detector_conditions(
             client, karabo_id, list(calibrations.keys()),
             {'parameters_conditions_attributes': condition},
             karabo_da=k_da, event_at=creation_time.isoformat())
 
         if not resp["success"]:
             print(f"ERROR: Constants {list(calibrations.values())} "
             f"were not retrieved, {resp['app_info']}")
             for cname in calibrations.values():
                 const_mdata[cname] = dict()
                 const_mdata[cname]["file-path"] = None
                 const_mdata[cname]["dataset-name"] = None
                 const_mdata[cname]["creation-time"] = None
             continue
 
         for ccv in resp["data"]:
             cc = ccv['calibration_constant']
             cname = calibrations[cc['calibration_id']]
             const_mdata[cname] = dict()
             const_mdata[cname]["file-path"] = str(Path(ccv['path_to_file']) / ccv['file_name'])
             const_mdata[cname]["dataset-name"] = ccv['data_set_name']
             const_mdata[cname]["creation-time"] = ccv['begin_at']
             pdu = ccv['physical_detector_unit']['physical_name']
 
         print('.', end='', flush=True)
     retrieved_constants[k_da]["physical-detector-unit"] = pdu
 metadata.save()
 
 total_time = perf_counter() - start
 print(f'{total_time:.1f}s')
 print(f"Stored retrieved constants in {metadata.filename}")
 ```

setup.py

@@ -108,8 +108,9 @@ install_requires = [
 
 if "readthedocs.org" not in sys.executable:
     install_requires += [
-        "iCalibrationDB @ git+ssh://git@git.xfel.eu:10022/detectors/cal_db_interactive.git@2.3.0",  # noqa
+        "iCalibrationDB @ git+ssh://git@git.xfel.eu:10022/detectors/cal_db_interactive.git@2.4.0",  # noqa
         "XFELDetectorAnalysis @ git+ssh://git@git.xfel.eu:10022/karaboDevices/pyDetLib.git@2.7.0",  # noqa
+        "CalParrot @ git+ssh://git@git.xfel.eu:10022/calibration/calparrot.git@0.1",  # noqa
     ]
 
 setup(

src/cal_tools/h5_copy_except.py

@@ -25,7 +25,7 @@ def copy_except_tree(src_group: h5py.Group, dest_group: h5py.Group, except_tree)
         if except_tree_part is True:  # Totally excluded
             pass
         elif except_tree_part is None:  # Not excluded
-            src_group.copy(name, dest_group, name)
+            src_group.copy(name, dest_group, name, without_attrs=True)
         else:  # Partially excluded
             src_subgroup = src_group[name]
             assert isinstance(src_subgroup, h5py.Group)

src/cal_tools/lpdlib.py

 import copy
 from typing import List, Optional, Tuple
+from warnings import warn
 
 import h5py
 import numpy as np
@@ -769,3 +770,23 @@ class LpdCorrections:
 
         self.initialize(offsets, rel_gains, rel_gains_b, bpixels, noises,
                         flat_fields)
+
+
+def get_mem_cell_order(run, sources) -> str:
+    """Load the memory cell order to use as a condition to find constants"""
+    res = set()
+    for source in sources:
+        cell_id_data = run[source, 'image.cellId'].drop_empty_trains()
+        if len(cell_id_data.train_ids) == 0:
+            continue  # No data for this module
+        cell_ids = cell_id_data[0].ndarray()
+        # Trailing comma required so e.g. "...,1" doesn't match "...,10"
+        res.add(",".join([str(c) for c in cell_ids.flatten()]) + ",")
+
+    if len(res) > 1:
+        warn("Memory cell order varies between detector modules: "
+             "; ".join([f"{s[:10]}...{s[-10:]}" for s in res]))
+    elif not res:
+        raise ValueError("Couldn't find memory cell order for any modules")
+
+    return res.pop()

src/cal_tools/restful_config.py

 from pathlib import Path
 
-from calibration_client import CalibrationClient
 from dynaconf import Dynaconf
 
 config_dir = Path(__file__).parent.resolve()
@@ -17,6 +16,8 @@ restful_config = Dynaconf(
 
 
 def calibration_client():
+    from calibration_client import CalibrationClient
+
     # Create client for CalCat.
     calcat_config = restful_config.get('calcat')
     return CalibrationClient(

src/xfel_calibrate/calibrate.py

@@ -385,11 +385,14 @@ class JobArgs:
         """Run this job in a local process, return exit status"""
         return call(self.format_cmd(python), cwd=work_dir)
 
-    def submit_job(self, work_dir, python, slurm_opts, after_ok=(), after_any=()):
+    def submit_job(
+            self, work_dir, python, slurm_opts, after_ok=(), after_any=(), env=None
+    ):
         """Submit this job to Slurm, return its job ID"""
         cmd = slurm_opts.get_launcher_command(work_dir, after_ok, after_any)
         cmd += self.format_cmd(python)
-        output = check_output(cmd, cwd=work_dir).decode('utf-8')
+        # sbatch propagates environment variables into the job by default
+        output = check_output(cmd, cwd=work_dir, env=env).decode('utf-8')
         return output.partition(';')[0].strip()  # job ID
 
 
@@ -440,7 +443,7 @@ class JobChain:
                 'steps': [step.to_dict() for step in self.steps]
             }, f, indent=2)
 
-    def submit_jobs(self, slurm_opts: SlurmOptions):
+    def submit_jobs(self, slurm_opts: SlurmOptions, env=None):
         """Submit these jobs to Slurm, return a list of job IDs
 
         Slurm dependencies are used to manage the sequence of jobs.
@@ -453,7 +456,9 @@ class JobChain:
             step_job_ids = []
             kw = {('after_any' if step.after_error else 'after_ok'): dep_job_ids}
             for job_desc in step.jobs:
-                jid = job_desc.submit_job(self.work_dir, self.python, slurm_opts, **kw)
+                jid = job_desc.submit_job(
+                    self.work_dir, self.python, slurm_opts, env=env, **kw
+                )
                 step_job_ids.append(jid)
             dep_job_ids = step_job_ids
             all_job_ids.extend(step_job_ids)

src/xfel_calibrate/repeat.py

@@ -135,12 +135,12 @@ def main(argv=None):
     out_folder = parameters['out-folder']
     params_to_set = {'metadata_folder': "."}
     if args.out_folder:
-        out_folder = parameters['out-folder'] = args.out_folder
+        out_folder = parameters['out-folder'] = os.path.abspath(args.out_folder)
         params_to_set['out_folder'] = out_folder
     update_notebooks_params(cal_work_dir, params_to_set)
 
     if args.report_to:
-        report_to = args.report_to
+        report_to = os.path.abspath(args.report_to)
     else:  # Default to saving report in output folder
         report_to = str(Path(out_folder, f'xfel-calibrate-repeat-{run_uuid}'))
     cal_metadata['report-path'] = f'{report_to}.pdf'
@@ -159,12 +159,17 @@ def main(argv=None):
         job_chain.run_direct()
         joblist = []
     else:
+        # The queries to look up constants should all be the same as those
+        # from the previous calibration - tell CalParrot to warn if not.
+        env = os.environ.copy()
+        env['CALPARROT_NEW_QUERY'] = 'warn'
+
         joblist = job_chain.submit_jobs(SlurmOptions(
             partition=args.slurm_partition,
             mem=args.slurm_mem,
-        ))
+        ), env=env)
 
-    fmt_args = {'run_path': cal_work_dir,
+    fmt_args = {'cal_work_dir': cal_work_dir,
                 'out_path': out_folder,
                 'version': get_pycalib_version(),
                 'report_to': report_to,

tests/test_cal_tools.py

@@ -253,7 +253,7 @@ def test_raise_get_from_db(_agipd_const_cond):
     assert str(excinfo.value) == "Resource temporarily unavailable"
 
     # Wrong type for creation_time.
-    with pytest.raises(ValueError):
+    with pytest.raises(TypeError):
         _call_get_from_db(
             constant=constant, condition=condition,
             karabo_id=AGIPD_KARABO_ID, karabo_da="AGIPD00",

webservice/job_monitor.py

@@ -279,17 +279,16 @@ class JobsMonitor:
         log.debug("Update MDC for %s, %s: %s", r['action'], r['mymdc_id'], msg)
 
         if r['action'] == 'CORRECT':
-            status = 'A' if success else 'NA'  # Not-/Available
+            status = 'A' if success else 'E'  # Available/Error
             self.mymdc_update_run(r['mymdc_id'], msg, status)
         else:  # r['action'] == 'DARK'
             status = 'F' if success else 'E'  # Finished/Error
             self.mymdc_update_dark(r['mymdc_id'], msg, status)
 
-    def mymdc_update_run(self, run_id, msg, status='R'):
+    def mymdc_update_run(self, run_id, msg, status='IP'):
         """Update correction status in MyMdC"""
-        data = {'flg_cal_data_status': status,
-                'cal_pipeline_reply': msg}
-        if status != 'R':
+        data = {'cal_pipeline_reply': msg, 'flg_cal_data_status': status}
+        if status != 'IP':
             data['cal_last_end_at'] = datetime.now(tz=timezone.utc).isoformat()
         response = self.mdc.update_run_api(run_id, data)
         if response.status_code != 200:

webservice/update_mdc.py

@@ -27,7 +27,8 @@ parser.add_argument('--runs', type=parse_int_range,
 parser.add_argument('--rid', type=int, help='Run id from MDC')
 parser.add_argument('--msg', type=str, help='Message string to MDC',
                     default='Error while job submission')
-parser.add_argument('--really', help="Actually make changes (otherwise dry-run)")
+parser.add_argument('--really', action='store_true',
+                    help="Actually make changes (otherwise dry-run)")
 args = parser.parse_args()
 
 if args.conf_file is not None:

webservice/webservice.py

@@ -742,17 +742,20 @@ async def update_mdc_status(mdc: MetadataClient, action: str,
     https://git.xfel.eu/gitlab/detectors/pycalibration/wikis/MyMDC-Communication
     """
     if message.split(':')[0] in ('FAILED', 'WARN'):  # Errors
-        flag = 'NA' if action == 'correct' else 'E'
+        flag = 'E'
     elif message.split(':')[0] == 'SUCCESS':  # Success
-        flag = 'R' if action == 'correct' else 'IP'
+        flag = 'IP'
         if 'Uploaded' in message or 'Finished' in message:
             flag = 'A' if action == 'correct' else 'F'
     else:  # MDC Timeout
-        flag = 'NA' if action == 'correct' else 'T'
+        flag = 'E' if action == 'correct' else 'T'
 
     if action == 'correct':
         func = mdc.update_run_api
-        data = {'flg_cal_data_status': flag, 'cal_pipeline_reply': message}
+        data = {'cal_pipeline_reply': message}
+        # Don't send In Progress; job_monitor will send this when jobs start.
+        if flag != 'IP':
+            data['flg_cal_data_status'] = flag
 
     elif action == 'dark_request':
         func = mdc.update_dark_run_api