cal_tools/cal_tools/agipdlib.py

@@ -215,6 +215,7 @@ class AgipdCorrections:
         self.mg_hard_threshold = 100
         self.hg_hard_threshold = 100
         self.noisy_adc_threshold = 0.25
+        self.ff_gain = 1
 
         # Shared variables for data and constants
         self.shared_dict = []
@@ -955,15 +956,32 @@ class AgipdCorrections:
         if self.corr_bools.get("xray_corr"):
             bpixels |= cons_data["BadPixelsFF"].astype(np.uint32)[..., :bpixels.shape[2], None]  # noqa
             slopesFF = cons_data["SlopesFF"]
+            # This could be used for backward compatibility
+            # for very old SlopesFF constants
             if len(slopesFF.shape) == 4:
                 slopesFF = slopesFF[..., 0]
-            # Memory cell resolved xray_cor correction
-            xray_cor = slopesFF  # (128, 512, mem_cells)
-
-            # relative X-ray correction is normalized by the median
-            # of all pixels
-
-            xray_cor /= np.nanmedian(xray_cor)
+            # This is for backward compatability for old FF constants
+            # (128, 512, mem_cells)
+            if slopesFF.shape[-1] == 2:
+                xray_cor = np.squeeze(slopesFF[...,0])
+                xray_cor_med = np.nanmedian(xray_cor)
+                xray_cor[np.isnan(xray_cor)]= xray_cor_med
+                xray_cor[(xray_cor<0.8) | (xray_cor>1.2)] = xray_cor_med
+                xray_cor = np.dstack([xray_cor]*self.max_cells)
+            else:
+                # Memory cell resolved xray_cor correction
+                xray_cor = slopesFF  # (128, 512, mem_cells)
+                if xray_cor.shape[-1] < self.max_cells:
+                    # In case of having new constant with less memory cells,
+                    # due to lack of enough FF data or during development.
+                    # xray_cor should be expanded by last memory cell.
+                    xray_cor = np.dstack(xray_cor,
+                                         np.dstack([xray_cor[..., -1]]
+                                            *  (self.max_cells - xray_cor.shape[-1])))  # noqa
+                # This is already done for old constants,
+                # but new constant is absolute and we need to have
+                # global ADU output for the moment
+                xray_cor /= self.ff_gain
 
             self.xray_cor[module_idx][...] = xray_cor.transpose()[...]
 
@@ -1078,22 +1096,6 @@ class AgipdCorrections:
                 # cells are used.
                 with h5py.File(mdata["file-path"], "r") as cf:
                     cons_data[cname] = np.copy(cf[f"{dname}/{cname}/0/data"])
-                    shape = cons_data[cname].shape  # (128, 512, mem_cells)
-                    extra_dims = shape[:2] + (self.max_cells-shape[2], )
-
-                    if extra_dims[-1] != 0 and cname == "BadPixelsFF":
-                        extra_temp = np.zeros(extra_dims, dtype=np.int32)
-                        cons_data[cname] = np.concatenate(
-                            (cons_data[cname], extra_temp), axis=2)
-                        print('An extra dimension was added to the constants '
-                              'for the benefit of BadPixelsFF')
-
-                    if extra_dims[-1] != 0 and cname == "SlopesFF":
-                        extra_temp = np.ones(extra_dims, dtype=np.float32)
-                        cons_data[cname] = np.concatenate(
-                            (cons_data[cname], extra_temp), axis=2)
-                        print('An extra dimension was added to the constants '
-                              'for the benefit of SlopesFF')
             else:
                 # Create empty constant using the list elements
                 cons_data[cname] = getattr(np, mdata["file-path"][0])(mdata["file-path"][1])  # noqa

notebooks/AGIPD/AGIPD_Correct_and_Verify.ipynb

 %% Cell type:markdown id: tags:
 
 # AGIPD Offline Correction #
 
 Author: European XFEL Detector Group, Version: 2.0
 
 Offline Calibration for the AGIPD Detector
 
 %% Cell type:code id: tags:
 
 ``` python
 in_folder = "/gpfs/exfel/exp/HED/202031/p900174/raw" # the folder to read data from, required
 out_folder = "/gpfs/exfel/data/scratch/ahmedk/test/hibef_agipd2"  # the folder to output to, required
 sequences = [-1] # sequences to correct, set to -1 for all, range allowed
 modules = [-1] # modules to correct, set to -1 for all, range allowed
 run = 155 # runs to process, required
 
 karabo_id = "HED_DET_AGIPD500K2G" # 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/' # path in the HDF5 file to images
 h5path_idx = 'INDEX/{}/DET/{}:xtdf/' # path in the HDF5 file to images
 h5path_ctrl = '/CONTROL/{}/MDL/FPGA_COMP' # path to control information
 karabo_id_control = "HED_EXP_AGIPD500K2G" # karabo-id for control device
 karabo_da_control = 'AGIPD500K2G00' # karabo DA for control infromation
 
 slopes_ff_from_files = "" # Path to locally stored SlopesFF and BadPixelsFF constants
 
 use_dir_creation_date = True # use the creation data of the input dir for database queries
 cal_db_interface = "tcp://max-exfl016:8015#8045" # the database interface to use
 cal_db_timeout = 30000 # in milli seconds
 creation_date_offset = "00:00:00" # add an offset to creation date, e.g. to get different constants
 
 max_cells = 0 # number of memory cells used, set to 0 to automatically infer
 bias_voltage = 300 # Bias voltage
 acq_rate = 0. # the detector acquisition rate, use 0 to try to auto-determine
 gain_setting = 0.1 # the gain setting, use 0.1 to try to auto-determine
 photon_energy = 9.2 # photon energy in keV
 overwrite = True # set to True if existing data should be overwritten
 max_pulses = [0, 500, 1] # range list [st, end, step] of maximum pulse indices within a train. 3 allowed maximum list input elements.
 mem_cells_db = 0 # set to a value different than 0 to use this value for DB queries
 cell_id_preview = 1 # cell Id used for preview in single-shot plots
 
 # Correction parameters
 blc_noise_threshold = 5000 # above this mean signal intensity now baseline correction via noise is attempted
 cm_dark_fraction = 0.66 # threshold for fraction of  empty pixels to consider module enough dark to perform CM correction
 cm_dark_range = [-50.,30] # range for signal value ADU for pixel to be consider as a dark pixel
 cm_n_itr = 4 # number of iterations for common mode correction
 hg_hard_threshold = 1000 # threshold to force medium gain offset subtracted pixel to high gain
 mg_hard_threshold = 1000 # threshold to force medium gain offset subtracted pixel from low to medium gain
 noisy_adc_threshold = 0.25 # threshold to mask complete adc
+ff_gain = 7.2 # conversion gain for absolute FlatField constants, while applying xray_gain
 
 # Correction Booleans
 only_offset = False # Apply only Offset correction. if False, Offset is applied by Default. if True, Offset is only applied.
 rel_gain = False # do relative gain correction based on PC data
 xray_gain = False # do relative gain correction based on xray data
 blc_noise = False # if set, baseline correction via noise peak location is attempted
 blc_stripes = False # if set, baseline corrected via stripes
 blc_hmatch = False # if set, base line correction via histogram matching is attempted
 match_asics = False # if set, inner ASIC borders are matched to the same signal level
 adjust_mg_baseline = False # adjust medium gain baseline to match highest high gain value
 zero_nans = False # set NaN values in corrected data to 0
 zero_orange = False # set to 0 very negative and very large values in corrected data
 blc_set_min = False # Shift to 0 negative medium gain pixels after offset corr
 corr_asic_diag = False # if set, diagonal drop offs on ASICs are correted
 force_hg_if_below = False # set high gain if mg offset subtracted value is below hg_hard_threshold
 force_mg_if_below = False # set medium gain if mg offset subtracted value is below mg_hard_threshold
 mask_noisy_adc = False # Mask entire ADC if they are noise above a relative threshold
 common_mode = False # Common mode correction
 melt_snow = False # Identify (and optionally interpolate) 'snowy' pixels
 mask_zero_std = False # Mask pixels with zero standard deviation across train
 low_medium_gap = False # 5 sigma separation in thresholding between low and medium gain
 
 # Paralellization parameters
 chunk_size = 1000 # Size of chunk for image-weise correction
 chunk_size_idim = 1  # chunking size of imaging dimension, adjust if user software is sensitive to this.
 n_cores_correct = 16 # Number of chunks to be processed in parallel
 n_cores_files = 4 # Number of files to be processed in parallel
 sequences_per_node = 2 # number of sequence files per cluster node if run as slurm job, set to 0 to not run SLURM parallel
 
 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 copy
 from datetime import timedelta
 from dateutil import parser
 import gc
 import glob
 import itertools
 from IPython.display import HTML, display, Markdown, Latex
 import math
 from multiprocessing import Pool
 import os
 import re
 import sys
 import traceback
 from time import time, sleep, perf_counter
 import tabulate
 import warnings
 warnings.filterwarnings('ignore')
 import yaml
 
 from extra_geom import AGIPD_1MGeometry, AGIPD_500K2GGeometry
 from extra_data import RunDirectory, stack_detector_data
 from iCalibrationDB import Detectors
 from mpl_toolkits.mplot3d import Axes3D
 from matplotlib.ticker import LinearLocator, FormatStrFormatter
 from matplotlib.colors import LogNorm
 from matplotlib import cm as colormap
 import matplotlib.pyplot as plt
 import matplotlib
 matplotlib.use("agg")
 %matplotlib inline
 import numpy as np
 import seaborn as sns
 sns.set()
 sns.set_context("paper", font_scale=1.4)
 sns.set_style("ticks")
 
 from cal_tools.agipdlib import (AgipdCorrections, get_acq_rate,
                                 get_gain_setting, get_num_cells)
 from cal_tools.cython import agipdalgs as calgs
 from cal_tools.ana_tools import get_range
 from cal_tools.enums import BadPixels
 from cal_tools.tools import get_dir_creation_date, map_modules_from_folder
 from cal_tools.step_timing import StepTimer
 
 import seaborn as sns
 sns.set()
 sns.set_context("paper", font_scale=1.4)
 sns.set_style("ticks")
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Evaluated parameters ##
 
 %% Cell type:code id: tags:
 
 ``` python
 # Fill dictionaries comprising bools and arguments for correction and data analysis
 
 # Here the herarichy and dependability for correction booleans are defined
 corr_bools = {}
 
 # offset is at the bottom of AGIPD correction pyramid.
 corr_bools["only_offset"] = only_offset
 
 # Dont apply any corrections if only_offset is requested
 if not only_offset:
     corr_bools["adjust_mg_baseline"] = adjust_mg_baseline
     corr_bools["rel_gain"] = rel_gain
     corr_bools["xray_corr"] = xray_gain
     corr_bools["blc_noise"] = blc_noise
     corr_bools["blc_stripes"] = blc_stripes
     corr_bools["blc_hmatch"] = blc_hmatch
     corr_bools["blc_set_min"] = blc_set_min
     corr_bools["match_asics"] = match_asics
     corr_bools["corr_asic_diag"] = corr_asic_diag
     corr_bools["zero_nans"] = zero_nans
     corr_bools["zero_orange"] = zero_orange
     corr_bools["mask_noisy_adc"] = mask_noisy_adc
     corr_bools["force_hg_if_below"] = force_hg_if_below
     corr_bools["force_mg_if_below"] = force_mg_if_below
     corr_bools["common_mode"] = common_mode
     corr_bools["melt_snow"] = melt_snow
     corr_bools["mask_zero_std"] = mask_zero_std
     corr_bools["low_medium_gap"] = low_medium_gap
 
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 if in_folder[-1] == "/":
     in_folder = in_folder[:-1]
 if sequences[0] == -1:
     sequences = None
 
 control_fname = f'{in_folder}/r{run:04d}/RAW-R{run:04d}-{karabo_da_control}-S00000.h5'
 h5path_ctrl = h5path_ctrl.format(karabo_id_control)
 h5path = h5path.format(karabo_id, receiver_id)
 h5path_idx = h5path_idx.format(karabo_id, receiver_id)
 
 print(f'Path to control file {control_fname}')
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # Create output folder
 os.makedirs(out_folder, exist_ok=overwrite)
 
 # Evaluate detector instance for mapping
 instrument = karabo_id.split("_")[0]
 if instrument == "SPB":
     dinstance = "AGIPD1M1"
     nmods = 16
 elif instrument == "MID":
     dinstance = "AGIPD1M2"
     nmods = 16
 # TODO: Remove DETLAB
 elif instrument == "HED" or instrument == "DETLAB":
     dinstance = "AGIPD500K"
     nmods = 8
 
 # Evaluate requested modules
 if karabo_da[0] == '-1':
     if modules[0] == -1:
         modules = list(range(nmods))
     karabo_da = ["AGIPD{:02d}".format(i) for i in modules]
 else:
     modules = [int(x[-2:]) for x in karabo_da]
 
 def mod_name(modno):
     return f"Q{modno // 4 + 1}M{modno % 4 + 1}"
 
 print("Process modules: ", ', '.join(
     [mod_name(x) for x in modules]))
 print(f"Detector in use is {karabo_id}")
 print(f"Instrument {instrument}")
 print(f"Detector instance {dinstance}")
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # Display Information about the selected pulses indices for correction.
 pulses_lst = list(range(*max_pulses)) if not (len(max_pulses)==1 and max_pulses[0]==0) else max_pulses
 
 try:
     if len(pulses_lst) > 1:
         print("A range of {} pulse indices is selected: from {} to {} with a step of {}"
                .format(len(pulses_lst), pulses_lst[0] , pulses_lst[-1] + (pulses_lst[1] - pulses_lst[0]),
                        pulses_lst[1] - pulses_lst[0]))
     else:
         print("one pulse is selected: a pulse of idx {}".format(pulses_lst[0]))
 except Exception as e:
     raise ValueError('max_pulses input Error: {}'.format(e))
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # set everything up filewise
 mmf = map_modules_from_folder(in_folder, run, path_template,
                               karabo_da, sequences)
 mapped_files, mod_ids, total_sequences, sequences_qm, _ = mmf
 file_list = []
 
 # ToDo: Split table over pages
 print(f"Processing a total of {total_sequences} sequence files in chunks of {n_cores_files}")
 table = []
 ti = 0
 for k, files in mapped_files.items():
     i = 0
     for f in list(files.queue):
         file_list.append(f)
         if i == 0:
             table.append((ti, k, i, f))
         else:
             table.append((ti, "", i,  f))
         i += 1
         ti += 1
 md = display(Latex(tabulate.tabulate(table, tablefmt='latex',
                                      headers=["#", "module", "# module", "file"])))
 file_list = sorted(file_list, key=lambda name: name[-10:])
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 filename = file_list[0]
 channel = int(re.findall(r".*-AGIPD([0-9]+)-.*", filename)[0])
 
 # Evaluate number of memory cells
 mem_cells = get_num_cells(filename, karabo_id, channel)
 if mem_cells is None:
     raise ValueError(f"No raw images found in {filename}")
 
 mem_cells_db = mem_cells if mem_cells_db == 0 else mem_cells_db
 max_cells = mem_cells if max_cells == 0 else max_cells
 
 # Evaluate aquisition rate
 if acq_rate == 0:
     acq_rate = get_acq_rate((filename, karabo_id, channel))
 
 print(f"Maximum memory cells to calibrate: {max_cells}")
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # Evaluate creation time
 creation_time = None
 if use_dir_creation_date:
     creation_time = get_dir_creation_date(in_folder, run)
     offset = parser.parse(creation_date_offset)
     delta = timedelta(hours=offset.hour,
                       minutes=offset.minute, seconds=offset.second)
     creation_time += delta
 
 # Evaluate gain setting
 if gain_setting == 0.1:
     if creation_time.replace(tzinfo=None) < parser.parse('2020-01-31'):
         print("Set gain-setting to None for runs taken before 2020-01-31")
         gain_setting = None
     else:
         try:
             gain_setting = get_gain_setting(control_fname, h5path_ctrl)
         except Exception as e:
             print(f'ERROR: while reading gain setting from: \n{control_fname}')
             print(e)
             print("Set gain setting to 0")
             gain_setting = 0
 
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 print(f"Using {creation_time} as creation time")
 print(f"Operating conditions are:\n• Bias voltage: {bias_voltage}\n• Memory cells: {mem_cells_db}\n"
               f"• Acquisition rate: {acq_rate}\n• Gain setting: {gain_setting}\n• Photon Energy: {photon_energy}\n")
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Data processing ##
 
 %% Cell type:code id: tags:
 
 ``` python
 agipd_corr = AgipdCorrections(max_cells, max_pulses,
                               h5_data_path=h5path,
                               h5_index_path=h5path_idx,
                               corr_bools=corr_bools)
 
 agipd_corr.baseline_corr_noise_threshold = -blc_noise_threshold
 agipd_corr.hg_hard_threshold = hg_hard_threshold
 agipd_corr.mg_hard_threshold = mg_hard_threshold
 
 agipd_corr.cm_dark_min = cm_dark_range[0]
 agipd_corr.cm_dark_max = cm_dark_range[1]
 agipd_corr.cm_dark_fraction = cm_dark_fraction
 agipd_corr.cm_n_itr = cm_n_itr
 agipd_corr.noisy_adc_threshold = noisy_adc_threshold
+agipd_corr.ff_gain = ff_gain
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # Retrieve calibration constants to RAM
 agipd_corr.allocate_constants(modules, (3, mem_cells_db, 512, 128))
 
 const_yaml = None
 if os.path.isfile(f'{out_folder}/retrieved_constants.yml'):
     with open(f'{out_folder}/retrieved_constants.yml', "r") as f:
         const_yaml = yaml.safe_load(f.read())
 
 # retrive constants
 def retrieve_constants(mod):
     """
     Retrieve calibration constants and load them to shared memory
 
     Metadata for constants is taken from yml file or retrieved from the DB
     """
     device = getattr(getattr(Detectors, dinstance), mod_name(mod))
     err = ''
     try:
         # check if there is a yaml file in out_folder that has the device constants.
         if const_yaml and device.device_name in const_yaml:
             when = agipd_corr.initialize_from_yaml(const_yaml, mod, device)
         else:
             when = agipd_corr.initialize_from_db(cal_db_interface, creation_time, mem_cells_db, bias_voltage,
                                                  photon_energy, gain_setting, acq_rate, mod, device, False)
     except Exception as e:
         err = f"Error: {e}\nError traceback: {traceback.format_exc()}"
         when = None
     return err, mod, when, device.device_name
 
 
 ts = perf_counter()
 with Pool(processes=len(modules)) as pool:
     const_out = pool.map(retrieve_constants, modules)
 print(f"Constants were loaded in {perf_counter()-ts:.01f}s")
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # allocate memory for images and hists
 n_images_max = max_cells*256
 data_shape = (n_images_max, 512, 128)
 agipd_corr.allocate_images(data_shape, n_cores_files)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 def batches(l, batch_size):
     """Group a list into batches of (up to) batch_size elements"""
     start = 0
     while start < len(l):
         yield l[start:start + batch_size]
         start += batch_size
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 def imagewise_chunks(img_counts):
     """Break up the loaded data into chunks of up to chunk_size
 
     Yields (file data slot, start index, stop index)
     """
     for i_proc, n_img in enumerate(img_counts):
         n_chunks = math.ceil(n_img / chunk_size)
         for i in range(n_chunks):
             yield i_proc, i * n_img // n_chunks, (i+1) * n_img // n_chunks
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 step_timer = StepTimer()
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 with Pool() as pool:
     for file_batch in batches(file_list, n_cores_files):
         # TODO: Move some printed output to logging or similar
         print(f"Processing next {len(file_batch)} files:")
         for file_name in file_batch:
             print(" ", file_name)
         step_timer.start()
 
         img_counts = pool.starmap(agipd_corr.read_file, enumerate(file_batch))
         step_timer.done_step('Loading data from files')
 
         # Evaluate zero-data-std mask
         pool.starmap(agipd_corr.mask_zero_std, itertools.product(
             range(len(file_batch)), np.array_split(np.arange(agipd_corr.max_cells), n_cores_correct)
         ))
         step_timer.done_step('Mask 0 std')
 
         # Perform offset image-wise correction
         pool.starmap(agipd_corr.offset_correction, imagewise_chunks(img_counts))
         step_timer.done_step("Offset correction")
 
 
         if blc_noise or blc_stripes or blc_hmatch:
             # Perform image-wise correction
             pool.starmap(agipd_corr.baseline_correction, imagewise_chunks(img_counts))
             step_timer.done_step("Base-line shift correction")
 
         if common_mode:
             # Perform cross-file correction parallel over asics
             pool.starmap(agipd_corr.cm_correction, itertools.product(
                 range(len(file_batch)), range(16)  # 16 ASICs per module
             ))
             step_timer.done_step("Common-mode correction")
 
         # Perform image-wise correction
         pool.starmap(agipd_corr.gain_correction, imagewise_chunks(img_counts))
         step_timer.done_step("Image-wise correction")
 
         # Save corrected data
         pool.starmap(agipd_corr.write_file, [
             (i_proc, file_name, os.path.join(out_folder, os.path.basename(file_name).replace("RAW", "CORR")))
             for i_proc, file_name in enumerate(file_batch)
         ])
         step_timer.done_step("Save")
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 print(f"Correction of {len(file_list)} files is finished")
 print(f"Total processing time {step_timer.timespan():.01f} s")
 print(f"Timing summary per batch of {n_cores_files} files:")
 step_timer.print_summary()
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # if there is a yml file that means a leading notebook got processed
 # and the reporting would be generated from it.
 fst_print = True
 
 to_store = []
 line = []
 for i, (error, modno, when, mod_dev) in enumerate(const_out):
     qm = mod_name(modno)
     # expose errors while applying correction
     if error:
         print("Error: {}".format(error) )
 
     if not const_yaml or mod_dev not in const_yaml:
         if fst_print:
             print("Constants are retrieved with creation time: ")
             fst_print = False
 
         line = [qm]
 
         # If correction is crashed
         if not error:
             print(f"{qm}:")
             for key, item in when.items():
                 if hasattr(item, 'strftime'):
                     item = item.strftime('%y-%m-%d %H:%M')
                 when[key] = item
                 print('{:.<12s}'.format(key), item)
 
         # Store few time stamps if exists
         # Add NA to keep array structure
         for key in ['Offset', 'SlopesPC', 'SlopesFF']:
             if when and key in when and when[key]:
                 line.append(when[key])
             else:
                 if error is not None:
                     line.append('Err')
                 else:
                     line.append('NA')
 
         if len(line) > 0:
             to_store.append(line)
 
 seq = sequences[0] if sequences else 0
 
 if len(to_store) > 0:
     with open(f"{out_folder}/retrieved_constants_s{seq}.yml","w") as fyml:
         yaml.safe_dump({"time-summary": {f"S{seq}":to_store}}, fyml)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 def do_3d_plot(data, edges, x_axis, y_axis):
     fig = plt.figure(figsize=(10, 10))
     ax = fig.gca(projection='3d')
 
     # Make data.
     X = edges[0][:-1]
     Y = edges[1][:-1]
     X, Y = np.meshgrid(X, Y)
     Z = data.T
 
     # Plot the surface.
     surf = ax.plot_surface(X, Y, Z, cmap=colormap.coolwarm,
                            linewidth=0, antialiased=False)
     ax.set_xlabel(x_axis)
     ax.set_ylabel(y_axis)
     ax.set_zlabel("Counts")
 
 
 def do_2d_plot(data, edges, y_axis, x_axis):
     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=max(10, np.max(data))))
     ax.set_xlabel(x_axis)
     ax.set_ylabel(y_axis)
     cb = fig.colorbar(im)
     cb.set_label("Counts")
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 def get_trains_data(run_folder, source, include, tid=None, path='*/DET/*', modules=16, fillvalue=np.nan):
     """
     Load single train for all module
 
     :param run_folder: Path to folder with data
     :param source: Data source to be loaded
     :param include: Inset of file name to be considered
     :param tid: Train Id to be loaded. First train is considered if None is given
     :param path: Path to find image data inside h5 file
 
     """
     run_data = RunDirectory(run_folder, include)
     if tid:
         tid, data = run_data.select('*/DET/*', source).train_from_id(tid)
         return tid, stack_detector_data(train=data, data=source, fillvalue=fillvalue, modules=modules)
     else:
         for tid, data in run_data.select('*/DET/*', source).trains(require_all=True):
             return tid, stack_detector_data(train=data, data=source, fillvalue=fillvalue, modules=modules)
     return None, None
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 if dinstance == "AGIPD500K":
     geom = AGIPD_500K2GGeometry.from_origin()
 else:
     geom = AGIPD_1MGeometry.from_quad_positions(quad_pos=[
         (-525, 625),
         (-550, -10),
         (520, -160),
         (542.5, 475),
     ])
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 include = '*S00000*' if sequences is None else f'*S{sequences[0]:05d}*'
 tid, corrected = get_trains_data(f'{out_folder}/', 'image.data', include, modules=nmods)
 _, gains = get_trains_data(f'{out_folder}/', 'image.gain', include, tid, modules=nmods)
 _, mask = get_trains_data(f'{out_folder}/', 'image.mask', include, tid, modules=nmods)
 _, blshift = get_trains_data(f'{out_folder}/', 'image.blShift', include, tid, modules=nmods)
 _, cellId = get_trains_data(f'{out_folder}/', 'image.cellId', include, tid, modules=nmods)
 _, pulseId = get_trains_data(f'{out_folder}/', 'image.pulseId', include, tid,
                              modules=nmods, fillvalue=0)
 _, raw = get_trains_data(f'{in_folder}/r{run:04d}/', 'image.data', include, tid, modules=nmods)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 display(Markdown(f'## Preview and statistics for {gains.shape[0]} images of the train {tid} ##\n'))
 ```
 
 %% Cell type:markdown id: tags:
 
 ### Signal vs. Analogue Gain ###
 
 %% Cell type:code id: tags:
 
 ``` python
 hist, bins_x, bins_y = calgs.histogram2d(raw[:,0,...].flatten().astype(np.float32),
                                          raw[:,1,...].flatten().astype(np.float32),
                                          bins=(100, 100),
                                          range=[[4000, 8192], [4000, 8192]])
 do_2d_plot(hist, (bins_x, bins_y), "Signal (ADU)", "Analogue gain (ADU)")
 do_3d_plot(hist, (bins_x, bins_y), "Signal (ADU)", "Analogue gain (ADU)")
 ```
 
 %% Cell type:markdown id: tags:
 
 ### Signal vs. Digitized Gain ###
 
 The following plot shows plots signal vs. digitized gain
 
 %% Cell type:code id: tags:
 
 ``` python
 hist, bins_x, bins_y = calgs.histogram2d(corrected.flatten().astype(np.float32),
                                          gains.flatten().astype(np.float32), bins=(100, 3),
                                          range=[[-50, 8192], [0, 3]])
 do_2d_plot(hist, (bins_x, bins_y), "Signal (ADU)", "Gain bit value")
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 print(f"Gain statistics in %")
 table = [[f'{gains[gains==0].size/gains.size*100:.02f}',
           f'{gains[gains==1].size/gains.size*100:.03f}',
           f'{gains[gains==2].size/gains.size*100:.03f}']]
 md = display(Latex(tabulate.tabulate(table, tablefmt='latex',
                                      headers=["High", "Medium", "Low"])))
 ```
 
 %% Cell type:markdown id: tags:
 
 ### Intensity per Pulse ###
 
 %% Cell type:code id: tags:
 
 ``` python
 pulse_range = [np.min(pulseId[pulseId>=0]), np.max(pulseId[pulseId>=0])]
 
 mean_data = np.nanmean(corrected, axis=(2, 3))
 hist, bins_x, bins_y = calgs.histogram2d(mean_data.flatten().astype(np.float32),
                                       pulseId.flatten().astype(np.float32),
                                       bins=(100, int(pulse_range[1])),
                                       range=[[-50, 1000], pulse_range])
 
 do_2d_plot(hist, (bins_x, bins_y), "Signal (ADU)", "Pulse id")
 do_3d_plot(hist, (bins_x, bins_y), "Signal (ADU)", "Pulse id")
 
 hist, bins_x, bins_y = calgs.histogram2d(mean_data.flatten().astype(np.float32),
                                       pulseId.flatten().astype(np.float32),
                                       bins=(100,  int(pulse_range[1])),
                                       range=[[-50, 200000], pulse_range])
 
 do_2d_plot(hist, (bins_x, bins_y), "Signal (ADU)", "Pulse id")
 do_3d_plot(hist, (bins_x, bins_y), "Signal (ADU)", "Pulse id")
 ```
 
 %% Cell type:markdown id: tags:
 
 ### Baseline shift ###
 
 Estimated base-line shift with respect to the total ADU counts of corrected image.
 
 %% Cell type:code id: tags:
 
 ``` python
 fig = plt.figure(figsize=(20, 10))
 ax = fig.add_subplot(111)
 h = ax.hist(blshift.flatten(), bins=100, log=True)
 _ = plt.xlabel('Baseline shift [ADU]')
 _ = plt.ylabel('Counts')
 _ = ax.grid()
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 fig = plt.figure(figsize=(10, 10))
 corrected_ave = np.nansum(corrected, axis=(2, 3))
 plt.scatter(corrected_ave.flatten()/10**6, blshift.flatten(), s=0.9)
 plt.xlim(-1, 1000)
 plt.grid()
 plt.xlabel('Illuminated corrected [MADU] ')
 _ = plt.ylabel('Estimated baseline shift [ADU]')
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 display(Markdown('### Raw preview ###\n'))
 display(Markdown(f'Mean over images of the RAW data\n'))
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 fig = plt.figure(figsize=(20, 10))
 ax = fig.add_subplot(111)
 data = np.mean(raw[:, 0, ...], axis=0)
 vmin, vmax = get_range(data, 5)
 ax = geom.plot_data_fast(data, ax=ax, cmap="jet", vmin=vmin, vmax=vmax)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 display(Markdown(f'Single shot of the RAW data from cell {np.max(cellId[cell_id_preview])} \n'))
 fig = plt.figure(figsize=(20, 10))
 ax = fig.add_subplot(111)
 vmin, vmax = get_range(raw[cell_id_preview, 0, ...], 5)
 ax = geom.plot_data_fast(raw[cell_id_preview, 0, ...], ax=ax, cmap="jet", vmin=vmin, vmax=vmax)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 display(Markdown('### Corrected preview ###\n'))
 display(Markdown(f'A single shot image from cell {np.max(cellId[cell_id_preview])} \n'))
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 fig = plt.figure(figsize=(20, 10))
 ax = fig.add_subplot(111)
 vmin, vmax = get_range(corrected[cell_id_preview], 7, -50)
 vmin = - 50
 ax = geom.plot_data_fast(corrected[cell_id_preview], ax=ax, cmap="jet", vmin=vmin, vmax=vmax)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 fig = plt.figure(figsize=(20, 10))
 ax = fig.add_subplot(111)
 vmin, vmax = get_range(corrected[cell_id_preview], 5, -50)
 nbins = np.int((vmax + 50) / 2)
 h = ax.hist(corrected[cell_id_preview].flatten(),
             bins=nbins, range=(-50, vmax),
             histtype='stepfilled', log=True)
 _ = plt.xlabel('[ADU]')
 _ = plt.ylabel('Counts')
 _ = ax.grid()
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 display(Markdown('### Mean CORRECTED Preview ###\n'))
 display(Markdown(f'A mean across one train \n'))
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 fig = plt.figure(figsize=(20, 10))
 ax = fig.add_subplot(111)
 data = np.mean(corrected, axis=0)
 vmin, vmax = get_range(data, 7)
 ax = geom.plot_data_fast(data, ax=ax, cmap="jet", vmin=-50, vmax=vmax)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 fig = plt.figure(figsize=(20, 10))
 ax = fig.add_subplot(111)
 vmin, vmax = get_range(corrected, 10, -100)
 vmax = np.nanmax(corrected)
 if vmax > 50000:
     vmax=50000
 nbins = np.int((vmax + 100) / 5)
 h = ax.hist(corrected.flatten(), bins=nbins,
             range=(-100, vmax), histtype='step', log=True, label = 'All')
 _ = ax.hist(corrected[gains == 0].flatten(), bins=nbins, range=(-100, vmax),
             alpha=0.5, log=True, label='High gain', color='green')
 _ = ax.hist(corrected[gains == 1].flatten(), bins=nbins, range=(-100, vmax),
             alpha=0.5, log=True, label='Medium gain', color='red')
 _ = ax.hist(corrected[gains == 2].flatten(), bins=nbins,
             range=(-100, vmax), alpha=0.5, log=True, label='Low gain', color='yellow')
 _ = ax.legend()
 _ = ax.grid()
 _ = plt.xlabel('[ADU]')
 _ = plt.ylabel('Counts')
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 display(Markdown('### Maximum GAIN Preview ###\n'))
 display(Markdown(f'The per pixel maximum across one train for the digitized gain'))
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 fig = plt.figure(figsize=(20, 10))
 ax = fig.add_subplot(111)
 ax = geom.plot_data_fast(np.max(gains, axis=0), ax=ax,
                          cmap="jet", vmin=-1, vmax=3)
 ```
 
 %% 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, "{:016b}".format(item.value)))
 md = display(Latex(tabulate.tabulate(table, tablefmt='latex',
                                      headers=["Bad pixel type", "Bit mask"])))
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 display(Markdown(f'### Single Shot Bad Pixels ### \n'))
 display(Markdown(f'A single shot bad pixel map from cell {np.max(cellId[cell_id_preview])} \n'))
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 fig = plt.figure(figsize=(20, 10))
 ax = fig.add_subplot(111)
 ax = geom.plot_data_fast(np.log2(mask[cell_id_preview]), ax=ax, vmin=0, vmax=32, cmap="jet")
 ```
 
 %% Cell type:markdown id: tags:
 
 ### Percentage of Bad Pixels across one train  ###
 
 %% Cell type:code id: tags:
 
 ``` python
 fig = plt.figure(figsize=(20, 10))
 ax = fig.add_subplot(111)
 ax = geom.plot_data_fast(np.mean(mask>0, axis=0),
                          vmin=0, ax=ax, vmax=1, cmap="jet")
 ```
 
 %% Cell type:markdown id: tags:
 
 ### Percentage of Bad Pixels across one train. Only Dark Related ###
 
 %% Cell type:code id: tags:
 
 ``` python
 fig = plt.figure(figsize=(20, 10))
 ax = fig.add_subplot(111)
 cm = np.copy(mask)
 cm[cm > BadPixels.NO_DARK_DATA.value] = 0
 ax = geom.plot_data_fast(np.mean(cm>0, axis=0),
                          vmin=0, ax=ax, vmax=1, cmap="jet")
 ```

notebooks/AGIPD/AGIPD_Correct_and_Verify_Summary_NBC.ipynb

 %% Cell type:markdown id: tags:
 
 # Summary of the AGIPD offline correction #
 
 %% Cell type:code id: tags:
 
 ``` python
 run = 11 # runs to process, required
 out_folder = "/gpfs/exfel/data/scratch/ahmedk/test/AGIPD_Corr" # path to output to, required
+karabo_id = "SPB_DET_AGIPD1M-1" # karabo karabo_id
 modules = [-1]
+karabo_da = ['-1']  # a list of data aggregators names, Default [-1] for selecting all data aggregators
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 import os
 
 import dateutil.parser
 import glob
 import numpy as np
 import re
 import yaml
 import warnings
 warnings.filterwarnings('ignore')
 
 import matplotlib.pyplot as plt
 %matplotlib inline
 import tabulate
 from IPython.display import display, Markdown, Latex
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 if os.path.isfile(f'{out_folder}/retrieved_constants.yml'):
     with open(f"{out_folder}/retrieved_constants.yml","r") as fyml:
         main_dict = yaml.load(fyml)
 else:
     main_dict = {"time-summary":dict()}
 
-if modules[0] == -1:
-    modules = list(range(16))
+# Extracting Instrument string
+instrument = karabo_id.split("_")[0]
+# Evaluate detector instance for mapping
+if instrument == "SPB":
+    dinstance = "AGIPD1M1"
+    nmods = 16
+elif instrument == "MID":
+    dinstance = "AGIPD1M2"
+    nmods = 16
+# TODO: Remove DETLAB
+elif instrument == "HED" or instrument == "DETLAB":
+    dinstance = "AGIPD500K"
+    nmods = 8
+
+if karabo_da[0] == '-1':
+    if modules[0] == -1:
+        modules = list(range(nmods))
+    karabo_da = ["AGIPD{:02d}".format(i) for i in modules]
+else:
+    modules = [int(x[-2:]) for x in karabo_da]
 
 # This is needed only if AGIPD Correction notebook had no precorrection notebooks for retrieving constants
 # gather all generated sequence yml files for time summary of retrieved constant in retrieved_constants.yml
 fnames = sorted(glob.glob(f'{out_folder}/retrieved_constants_*yml'))
 for f in fnames:
     with open(f,"r") as fyml:
         fdict = yaml.load(fyml)
     # append different sequences's time summary to the main yaml
     for k, v in fdict["time-summary"].items():
         main_dict["time-summary"][k] = v
     os.remove(f)
 
 with open(f"{out_folder}/retrieved_constants.yml","w") as fyml:
         yaml.safe_dump(main_dict, fyml)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 with open(f"{out_folder}/retrieved_constants.yml","r") as fyml:
     time_summary = yaml.load(fyml)
 # check if pre-notebook has retrieved constants for all modules.
 const_times = []
 seq = []
 for k, v  in sorted(time_summary["time-summary"].items()):
     arr = np.array(v)
     arr = arr.reshape(arr.shape[0]//len(modules), len(modules), arr.shape[1])
     const_times = const_times + list(arr)
     seq.append(k)
 
 const_times = np.array(const_times)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # Function print summary of constant injection time
 # To reduce printouts only unique entries are shown.
 def const_table(const, pos):
     """
     Create a summary table for the creation time differences for
     the retrieved constants (Offset, SlopesPC, SlopesFF).
     """
     print(f"{const} were injected on: ")
 
     # catch timing difference in retrieve constants
     unique, idx, counts = np.unique(const_times[:,:,pos], return_inverse=True, return_counts=True)
     idx = idx.reshape((const_times.shape[0], len(modules)))
 
     table = []
     for i in range(0, counts.shape[0]):
         line = [ const_times[:,:,pos][idx==i][0]  ]
         mods = ''
         for i_s, s in enumerate(seq):
             if(const_times[i_s,:,0][idx[i_s]==i].shape[0] > 0):
                 mods = mods+ '{}: {}, '.format(s, const_times[i_s,:,0][idx[i_s]==i])
         line.append(mods)
         table.append(line)
 
     if counts.shape[0] == 1:
         table[0][1] = 'All modules'
     else:
         table[np.argmax(counts)][1] = 'Rest of the modules'
 
     md = display(Latex(tabulate.tabulate(table, tablefmt='latex',
                                          headers=["Time stamps", "Modules and sequences"])))
 for i_key, key in enumerate(['Offset', 'SlopesPC', 'SlopesFF']):
     if const_times.shape[2] > i_key+1:
         const_table(key, i_key+1)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 ```

notebooks/DSSC/Characterize_DSSC_Darks_NBC.ipynb

 %% Cell type:markdown id: tags:
 
 # DSSC Characterize Dark Images #
 
 Author: S. Hauf, Version: 0.1
 
 The following code analyzes a set of dark images taken with the DSSC detector to deduce detector offsets and noise. Data for the detector is presented in one run and don't acquire multiple gain stages.
 
 The notebook explicitely does what pyDetLib provides in its offset calculation method for streaming data.
 
 %% Cell type:code id: tags:
 
 ``` python
 cluster_profile = "noDB" # The ipcluster profile to use
 in_folder = "/gpfs/exfel/exp/SCS/202030/p900125/raw" # path to input data, required
 out_folder = "/gpfs/exfel/data/scratch/ahmedk/test/DSSC" # path to output to, required
 sequences = [0] # sequence files to evaluate.
 modules = [-1]  # modules to run for
 run = 222 # run numbr in which data was recorded, required
 
 karabo_id = "SCS_DET_DSSC1M-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 dir creation date for determining the creation time
 cal_db_interface = "tcp://max-exfl016:8020" # the database interface to use
 cal_db_timeout = 3000000 # timeout on caldb requests"
 local_output = True # output constants locally
 db_output = False # output constants to database
 
 mem_cells = 0 # number of memory cells used, set to 0 to automatically infer
 bias_voltage = 300 # detector bias voltage
 rawversion = 2 # RAW file format version
 
 thresholds_offset_sigma = 3. # thresholds in terms of n sigma noise for offset deduced bad pixels
 thresholds_offset_hard = [4000, 8500] # thresholds in absolute ADU terms for offset deduced bad pixels
 
 thresholds_noise_sigma = 5. # thresholds in terms of n sigma noise for offset deduced bad pixels
 thresholds_noise_hard = [4, 20] # thresholds in absolute ADU terms for offset deduced bad pixels
 
 instrument = "SCS" # the instrument
 high_res_badpix_3d = False # set this to True if you need high-resolution 3d bad pixel plots. Runtime: ~ 1h
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # imports and things that do not usually need to be changed
 from datetime import datetime
 import os
 import warnings
 warnings.filterwarnings('ignore')
 from collections import OrderedDict
 
 import h5py
 from ipyparallel import Client
 from IPython.display import display, Markdown, Latex
 import matplotlib
 matplotlib.use('agg')
 import matplotlib.pyplot as plt
 %matplotlib inline
 import numpy as np
 import tabulate
 
 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_random_db_interface,
                              map_gain_stages, parse_runs,
                              run_prop_seq_from_path,
                              save_const_to_h5, send_to_db)
 from iCalibrationDB import Constants, Conditions, Detectors, Versions
 view = Client(profile=cluster_profile)[:]
 view.use_dill()
 
 # make sure a cluster is running with ipcluster start --n=32, give it a while to start
 
 h5path = h5path.format(karabo_id, receiver_id)
 h5path_idx = h5path_idx.format(karabo_id, receiver_id)
 gain_names = ['High', 'Medium', 'Low']
 
 if karabo_da[0] == '-1':
     if modules[0] == -1:
         modules = list(range(16))
     karabo_da = ["DSSC{:02d}".format(i) for i in modules]
 else:
     modules = [int(x[-2:]) for x in karabo_da]
 
 max_cells = mem_cells
 
 offset_runs = OrderedDict()
 offset_runs["high"] = run
 
 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 of constant.")
 
 run, prop, seq = run_prop_seq_from_path(in_folder)
 
 dinstance = "DSSC1M1"
 
 print(f"Detector in use is {karabo_id}")
 
 cal_db_interface = get_random_db_interface(cal_db_interface)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 print("Parameters are:")
 print(f"Proposal: {prop}")
 print(f"Memory cells: {mem_cells}/{max_cells}")
 print("Runs: {}".format([ v for v in offset_runs.values()]))
 print(f"Sequences: {sequences}")
 print(f"Using DB: {db_output}")
 print(f"Input: {in_folder}")
 print(f"Output: {out_folder}")
 print(f"Bias voltage: {bias_voltage}V")
 file_loc = f'proposal:{prop} runs:{[ v for v in offset_runs.values()][0]}'
 ```
 
 %% Cell type:markdown id: tags:
 
 The following lines will create a queue of files which will the be executed module-parallel. Distinguishing between different gains.
 
 %% Cell type:code id: tags:
 
 ``` python
 # set everything up filewise
 os.makedirs(out_folder, exist_ok=True)
 gmf = map_gain_stages(in_folder, offset_runs, path_template, karabo_da, sequences)
 gain_mapped_files, total_sequences, total_file_size = gmf
 print(f"Will process a total of {total_sequences} file.")
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Calculate Offsets, Noise and Thresholds ##
 
 The calculation is performed per-pixel and per-memory-cell. Offsets are simply the median value for a set of dark data taken at a given gain, noise the standard deviation, and gain-bit values the medians of the gain array.
 
 %% Cell type:code id: tags:
 
 ``` python
 import copy
 from functools import partial
 def characterize_module(cells, bp_thresh, rawversion, karabo_id, h5path, h5path_idx, inp):
     import numpy as np
     import copy
     import h5py
     from cal_tools.enums import BadPixels
 
     from hashlib import blake2b
     import struct
     import binascii
 
     def get_num_cells(fname, h5path):
         with h5py.File(fname, "r") as f:
 
             cells = f[f"{h5path}/cellId"][()]
             maxcell = np.max(cells)
             options = [100, 200, 400, 500, 600, 700, 800]
             dists = np.array([(o-maxcell) for o in options])
             dists[dists<0] = 10000 # assure to always go higher
             return options[np.argmin(dists)]
 
     filename, channel = inp
 
     h5path = h5path.format(channel)
     h5path_idx = h5path_idx.format(channel)
 
     if cells == 0:
         cells = get_num_cells(filename, h5path)
 
     print(f"Using {cells} memory cells")
 
     pulseid_checksum = None
 
     thresholds_offset_hard, thresholds_offset_sigma, thresholds_noise_hard, thresholds_noise_sigma = bp_thresh
 
     infile = h5py.File(filename, "r", driver="core")
     if rawversion == 2:
         count = np.squeeze(infile[f"{h5path_idx}/count"])
         first = np.squeeze(infile[f"{h5path_idx}/first"])
         last_index = int(first[count != 0][-1]+count[count != 0][-1])
         first_index = int(first[count != 0][0])
         pulseids = infile[f"{h5path}/pulseId"][first_index:int(first[count != 0][1])]
         bveto = blake2b(pulseids.data, digest_size=8)
         pulseid_checksum = struct.unpack('d', binascii.unhexlify(bveto.hexdigest()))[0]
     else:
         status = np.squeeze(infile[f"{h5path_idx}/status"])
         if np.count_nonzero(status != 0) == 0:
             return
         last = np.squeeze(infile[f"{h5path_idx}/last"])
         first = np.squeeze(infile[f"{h5path_idx}/first"])
         last_index = int(last[status != 0][-1]) + 1
         first_index = int(first[status != 0][0])
     im = np.array(infile[f"{h5path}/data"][first_index:last_index,...])
     cellIds = np.squeeze(infile[f"{h5path}/cellId"][first_index:last_index,...])
 
     infile.close()
 
 
     im = im[:, 0, ...].astype(np.float32)
 
     im = np.rollaxis(im, 2)
     im = np.rollaxis(im, 2, 1)
 
     mcells = cells
     offset = np.zeros((im.shape[0], im.shape[1], mcells))
     noise = np.zeros((im.shape[0], im.shape[1], mcells))
 
     for cc in np.unique(cellIds[cellIds < mcells]):
         cellidx = cellIds == cc
         offset[...,cc] = np.median(im[..., cellidx], axis=2)
         noise[...,cc] = np.std(im[..., cellidx], axis=2)
 
 
     # 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
 
 
     return offset, noise, bp, cells, pulseid_checksum
 
 
 offset_g = OrderedDict()
 noise_g = OrderedDict()
 gain_g = OrderedDict()
 badpix_g = OrderedDict()
 gg = 0
 
 start = datetime.now()
 all_cells = []
 checksums = {}
 
 for gain, mapped_files in gain_mapped_files.items():
     inp = []
     dones = []
     for i in modules:
         qm = "Q{}M{}".format(i//4 +1, i % 4 + 1)
         if qm in mapped_files and not mapped_files[qm].empty():
             fname_in = mapped_files[qm].get()
             print("Process file: ", fname_in)
             dones.append(mapped_files[qm].empty())
         else:
             continue
         inp.append((fname_in, i))
 
     p = partial(characterize_module, max_cells,
                (thresholds_offset_hard, thresholds_offset_sigma,
                 thresholds_noise_hard, thresholds_noise_sigma), rawversion, karabo_id, h5path, h5path_idx)
     results = list(map(p, inp))
 
     for ii, r in enumerate(results):
         i = modules[ii]
         offset, noise,  bp, thiscell, pulseid_checksum = r
         all_cells.append(thiscell)
         qm = "Q{}M{}".format(i//4 +1, i % 4 + 1)
         if qm not in offset_g:
             offset_g[qm] = np.zeros((offset.shape[0], offset.shape[1], offset.shape[2]))
             noise_g[qm] = np.zeros_like(offset_g[qm])
 
             badpix_g[qm] = np.zeros_like(offset_g[qm], np.uint32)
             checksums[qm] = pulseid_checksum
 
         offset_g[qm][...] = offset
         noise_g[qm][...] = noise
         badpix_g[qm][...] = bp
     gg +=1
 
 if len(all_cells) > 0:
     max_cells = np.max(all_cells)
     print(f"Using {max_cells} memory cells")
 else:
     raise ValueError("0 processed memory cells. No raw data available.")
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # Retrieve existing constants for comparison
 clist = ["Offset", "Noise"]
 old_const = {}
 old_mdata = {}
 
 print('Retrieve pre-existing constants for comparison.')
 detinst = getattr(Detectors, dinstance)
 for qm in offset_g.keys():
     device = getattr(detinst, qm)
     for const in clist:
         condition = Conditions.Dark.DSSC(memory_cells=max_cells,
                                          bias_voltage=bias_voltage,
                                          pulseid_checksum=checksums[qm])
 
         data, mdata = get_from_db(device,
                                   getattr(Constants.DSSC, const)(),
                                   condition,
                                   None,
                                   cal_db_interface, creation_time=creation_time,
                                   verbosity=2, timeout=cal_db_timeout)
 
         old_const[const] = data
 
         if mdata is not None and data is not None:
             time = mdata.calibration_constant_version.begin_at
             old_mdata[const] = time.isoformat()
             os.makedirs(f'{out_folder}/old/', exist_ok=True)
             save_const_to_h5(device,
                              getattr(Constants.DSSC, const)(),
                              condition, data, file_loc, creation_time,
                              f'{out_folder}/old/')
         else:
             old_mdata[const] = "Not found"
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 res = OrderedDict()
 for i in modules:
     qm = f"Q{i//4+1}M{i%4+1}"
     try:
         res[qm] = {'Offset': offset_g[qm],
                    'Noise': noise_g[qm],
                    #TODO: No badpixelsdark, yet.
                    #'BadPixelsDark': badpix_g[qm]
                    }
     except Exception as e:
         print(f"Error: No constants for {qm}: {e}")
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # Push the same constant two different times.
 # One with the generated pulseID check sum setting for the offline calibration.
 # And another for the online calibration as it doesn't have this pulseID checksum, yet.
 md = None
 for dont_use_pulseIds in [True, False]:
     for qm in res.keys():
         detinst = getattr(Detectors, dinstance)
         device = getattr(detinst, qm)
 
         for const in res[qm].keys():
             dconst = getattr(Constants.DSSC, const)()
             dconst.data = res[qm][const]
 
             pidsum = None if dont_use_pulseIds else checksums[qm]
             # set the operating condition
             condition = Conditions.Dark.DSSC(memory_cells=max_cells,
                                              bias_voltage=bias_voltage,
                                              pulseid_checksum=pidsum)
             if db_output:
                 md = send_to_db(device, dconst, condition, file_loc,
                                 cal_db_interface, creation_time=creation_time, timeout=cal_db_timeout)
 
             if local_output:
                 # Don't save constant localy two times.
                 if dont_use_pulseIds:
                     md = save_const_to_h5(device, dconst, condition,
                                           dconst.data, file_loc,
                                           creation_time, out_folder)
                     print(f"Calibration constant {const} is stored locally.\n")
 
         if not dont_use_pulseIds:
             print("Constants parameter conditions are:\n")
             print(f"• memory_cells: {max_cells}\n• bias_voltage: {bias_voltage}\n"
                   f"• pulseid_checksum: {pidsum}\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
 mnames = []
 for i in modules:
     qm = f"Q{i//4+1}M{i % 4+1}"
-    display(Markdown(f'## Position of the module {mnames} and it\'s ASICs##'))
+    display(Markdown(f'## Position of the module {qm} and its ASICs##'))
     mnames.append(qm)
 
 show_processed_modules(dinstance=dinstance, constants=None, mnames=mnames, mode="position")
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Single-Cell Overviews ##
 
 Single cell overviews allow to identify potential effects on all memory cells, e.g. on 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 = 9
 gain = 0
 out_folder = None
 show_overview(res, cell, gain, out_folder=out_folder, infix="_{}".format(run))
 ```
 
 %% 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 show the results of bad pixel evaluation for all evaluated memory cells.
     Cells are stacked in the Z-dimension, while pixels values in x/y are rebinned 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.
 
     """))
     # set rebin_fac to 1 for avoiding rebining and
     # losing real values of badpixels(High resolution).
     gain = 0
     for mod, data in badpix_g.items():
         plot_badpix_3d(data, cols, title=mod, rebin_fac=2)
         plt.show()
 ```
 
 %% 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.
 
 %% Cell type:code id: tags:
 
 ``` python
 create_constant_overview(offset_g, "Offset (ADU)", max_cells, entries=1)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 create_constant_overview(noise_g, "Noise (ADU)", max_cells, 0, 100, entries=1)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 bad_pixel_aggregate_g = OrderedDict()
 for m, d in badpix_g.items():
     bad_pixel_aggregate_g[m] = d.astype(np.bool).astype(np.float)
 create_constant_overview(bad_pixel_aggregate_g, "Bad pixel fraction", max_cells, entries=1)
 ```
 
 %% 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
 header = ['Parameter',
           "New constant", "Old constant ",
           "New constant", "Old constant ",
           "New constant", "Old constant "]
 
 for const in ['Offset', 'Noise']:
     table = [['','High gain', 'High gain']]
     for qm in res.keys():
 
         data = np.copy(res[qm][const])
 
         if old_const[const] is not None:
             dataold = np.copy(old_const[const])
 
         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]]
             line.append('{:6.1f}'.format(f(data[...,gain])))
             if old_const[const] is not None:
                 line.append('{:6.1f}'.format(f(dataold[...,gain])))
             else:
                 line.append('-')
 
             table.append(line)
 
     display(Markdown('### {} [ADU], good and bad pixels ###'.format(const)))
     md = display(Latex(tabulate.tabulate(table, tablefmt='latex', headers=header)))
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 ```

notebooks/FastCCD/Characterize_Darks_NewDAQ_FastCCD_NBC_New_Common_Mode.ipynb

 %% Cell type:markdown id: tags:
 
 # FastCCD Dark Characterization
 
 Author: I. Klačková, S. Hauf, K. Setoodehnia and M. Cascella
 
 The following notebook provides dark image analysis of the FastCCD detector.
 
 Dark characterization evaluates offset and noise of the FastCCD detector, corrects the noise for Common Mode (CM), and defines bad pixels relative to offset and CM corrected noise. Bad pixels are then excluded and CM corrected noise is recalculated excluding the bad pixels. Resulting offset and CM corrected noise maps, as well as the bad pixel map are sent to the calibration database.
 
 %% Cell type:code id: tags:
 
 ``` python
 cluster_profile = "noDB" # ipcluster profile to use
 in_folder = "/gpfs/exfel/exp/SCS/201930/p900074/raw" # input folder, required
 out_folder = '/gpfs/exfel/data/scratch/ahmedk/test/fastccd' # output folder, required
 sequence = 0 # sequence file to use
 run = 351 # which run to read data from, required
 
 karabo_da = ['DA05'] # data aggregators
 karabo_id = "SCS_CDIDET_FCCD2M" # karabo prefix of PNCCD devices
 receiver_id = "FCCD" # inset for receiver devices
 path_template = 'RAW-R{:04d}-{}-S{{:05d}}.h5' # the template to use to access data
 h5path = '/INSTRUMENT/{}/DAQ/{}:daqOutput/data/image/pixels' # path to the data in the HDF5 file
 h5path_t = '/CONTROL/{}/CTRL/LSLAN/inputA/crdg/value'  # path to find temperature
 h5path_cntrl = '/RUN/{}/DET/FCCD'  # path to find control data
 
 use_dir_creation_date = True # To be used to retrieve calibration constants later on (for database time derivation)
 cal_db_interface = "tcp://max-exfl016:8020" # the calibration database interface to use
 cal_db_timeout = 300000 # timeout on calibration database requests
 db_output = False # Output constants to the calibration database
 local_output = True # output constants locally
 
 number_dark_frames = 0 # number of images to be used, if set to 0 all available images are used
 # The two operation modes for FastCCD have fixed names which cannot be changed:
 operation_mode = "FF" # FS stands for frame-store and FF for full-frame opeartion.
 temp_limits = 5 # to find calibration constants later on, the sensor temperature is allowed to vary by 5 units
 bad_pixel_offset_sigma = 5. # Any pixel whose offset is beyond 5 standard deviations, is a bad pixel
 bad_pixel_noise_sigma = 5. # Any pixel whose noise is beyond 5 standard deviations, is a bad pixel
 sigmaNoise = 5. # Any pixel whose signal exceeds 'sigmaNoise'*noiseCM (common mode corrected noise) will be masked
 fix_temperature = 0. # Fixed operation temperature in Kelvins. If set to 0, mean value of the data file's temperature is
                      # used.
 chunkSize = 100 # Number of images to read per chunk
 cpuCores = 40 # Specifies the number of running cpu cores
 commonModeAxis = 1 # Axis along which common mode will be calculated (0: along rows, 1: along columns)
 run_parallel = True # For parallel computation
 
 
 # According to our gain calibration using 55Fe source, we have the following conversion gains (e.g., 1 ADU = 6.3 e-):
 ADU_to_electron_upper_hg = 6.3 # for upper hemisphere and high gain
 ADU_to_electron_lower_hg = 6.4 # for lower hemisphere and high gain
 ADU_to_electron_upper_mg = 23.4 # for upper hemisphere and medium gain
 ADU_to_electron_lower_mg = 23.4 # for lower hemisphere and medium gain
 ADU_to_electron_upper_lg = 49.3 # for upper hemisphere and low gain
 ADU_to_electron_lower_lg = 47.3 # for lower hemisphere and low gain
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # Required Packages:
 import copy
 import datetime
 import time
 import os
 import warnings
 warnings.filterwarnings('ignore')
 
 import h5py
 from IPython.display import display, Markdown
 import matplotlib.pyplot as plt
 %matplotlib inline
 import numpy as np
 from prettytable import PrettyTable
 
-from cal_tools.tools import get_dir_creation_date, get_random_db_interface, save_const_to_h5
+from cal_tools.tools import (get_dir_creation_date, get_random_db_interface,
+                             save_const_to_h5, send_to_db)
 from cal_tools.enums import BadPixels
 from iCalibrationDB import Constants, Conditions, Detectors, Versions
 from iCalibrationDB.detectors import DetectorTypes
 from XFELDetAna import xfelpyanatools as xana
 from XFELDetAna import xfelpycaltools as xcal
 from XFELDetAna.detectors.fastccd import readerh5 as fastccdreaderh5
 from XFELDetAna.util import env
 env.iprofile = cluster_profile
 import XFELDetAna.xfelprofiler as xprof
 profiler = xprof.Profiler()
 profiler.disable()
 from XFELDetAna.xfelreaders import ChunkReader
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # Output Folder Creation:
 if not os.path.exists(out_folder):
     os.makedirs(out_folder)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # Number of Images:
 def nImagesOrLimit(nImages, limit):
     if limit == 0:
         return nImages
     else:
         return min(nImages, limit)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 proposal = list(filter(None, in_folder.strip('/').split('/')))[-2]
 file_loc = 'proposal:{} runs:{}'.format(proposal, run)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # Detector Operation Mode, Calibration Database Settings, and Some Initial Run Parameters & Paths:
 
 display(Markdown('### Initial Settings'))
 if operation_mode == "FS":
     x = 960 # rows of the FastCCD to analyze in FS mode
     y = 960 # columns of the FastCCD to analyze in FS mode
     print('\nYou are analyzing data in FS mode.')
 else:
     x = 1934 # rows of the FastCCD to analyze in FF mode
     y = 960 # columns of the FastCCD to analyze in FF mode
     print('\nYou are analyzing data in FF mode.')
 
 ped_dir = "{}/r{:04d}".format(in_folder, run)
 fp_name = path_template.format(run, karabo_da[0])
 fp_path = '{}/{}'.format(ped_dir, fp_name)
 filename = fp_path.format(sequence)
 h5path = h5path.format(karabo_id, receiver_id)
 h5path_t = h5path_t.format(karabo_id)
 h5path_cntrl = h5path_cntrl.format(karabo_id)
 
 creation_time = None
 if use_dir_creation_date:
     creation_time = get_dir_creation_date(in_folder, run)
 
 cal_db_interface = get_random_db_interface(cal_db_interface)
 print('Calibration database interface: {}'.format(cal_db_interface))
 print("Sending constants to the calibration database: {}".format(db_output))
 print("HDF5 path to data: {}".format(h5path))
 print("Run number: {}".format(run))
 print("Reading data from: {}".format(filename))
 if creation_time:
     print("Using {} as creation time".format(creation_time.isoformat()))
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # Reading Parameters such as Detector Bias, Gain, etc. from the Data:
 
 memoryCells = 1 # FastCCD has 1 memory cell
 sensorSize = [x, y]
 blockSize = [sensorSize[0]//2, sensorSize[1]] # Sensor area will be analysed according to blocksize
 xcal.defaultBlockSize = blockSize
 nImages = fastccdreaderh5.getDataSize(filename, h5path)[0] # Specifies total number of images to proceed
 nImages = nImagesOrLimit(nImages, number_dark_frames)
 profile = False
 
 with h5py.File(filename, 'r') as f:
     bias_voltage = int(f['{}/biasclock/bias/value'.format(h5path_cntrl)][0])
     det_gain = int(f['{}/exposure/gain/value'.format(h5path_cntrl)][0])
     integration_time = int(f['{}/exposure/exposure_time/value'.format(h5path_cntrl)][0])
     temperature = np.mean(f[h5path_t])
     temperature = round(temperature, 2)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # Printing the Parameters Read from the Data File:
 
 display(Markdown('### Evaluated Parameters'))
 print("Number of dark images to analyze:", nImages)
 
 gain_dict = {
         "high gain" : 8,
         "medium gain" : 2,
         "low gain" : 1,
         "auto gain" : 0
     }
 
 for gain, value in gain_dict.items():
     if det_gain == value:
         gain_setting = gain
 
 print("Bias voltage is {} V.".format(bias_voltage))
 print("Detector gain is set to x{} ({}).".format(det_gain, gain_setting))
 print("Detector integration time is set to {}".format(integration_time), 'ms.')
 
 
 if fix_temperature != 0.:
     print("Using a fixed temperature of {} K".format(fix_temperature))
 else:
     # This is needed while sending the
     # calibration constant to the DB later
     fix_temperature = temperature + 273.15
     print("Temperature is not fixed.")
     print("Mean temperature was {:0.2f} °C / {:0.2f} K".format(temperature, fix_temperature))
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # Reading Files in Chunks:
 
 # Chunk reader returns an iterator to access the data in the file within the ranges:
 
 reader = ChunkReader(filename, fastccdreaderh5.readData, nImages, chunkSize, path = h5path, pixels_x = sensorSize[0],
                      pixels_y = sensorSize[1],)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # Calculator:
 
 # noiseCal is a noise map calculator, which internally also produces a per-pixel mean map, i.e. an offset map:
 
 noiseCal = xcal.NoiseCalculator(sensorSize, memoryCells, cores=cpuCores, blockSize=blockSize, runParallel=run_parallel)
 ```
 
 %% Cell type:markdown id: tags:
 
 ### First Iteration
 
 Characterization of dark images with purpose to create dark maps (offset, noise and bad pixel maps) is an iterative process. Firstly, initial offset and noise maps are produced from raw dark data.
 
 %% Cell type:code id: tags:
 
 ``` python
 counter1 = 1 # To count how many "if data.shape[2] >= chunkSize" instances are there.
 counter2 = 0 # To count how many "if data.shape[2] < chunkSize" instances are there.
 chunkSize_new = 0 # See below
 
 for data in reader.readChunks():
     data = np.bitwise_and(data.astype(np.uint16), 0b0011111111111111).astype(np.float32)
     dx = np.count_nonzero(data, axis=(0, 1))
     data = data[:,:,dx != 0]
     # Some sequences may have less than 500 frames in them. To find out how many images there are, we will temporarily
     # change chunkSize to be the same as whatever number of frames the last chunk of data has:
     if data.shape[2] < chunkSize:
         chunkSize_new = data.shape[2]
         print("Number of images are less than chunkSize. chunkSize is temporarily changed to {}."
               .format(chunkSize_new))
         images = images + chunkSize_new
         counter2 += 1
     else:
         images = counter1*chunkSize + counter2*chunkSize_new
         counter1 += 1
 
     noiseCal.fill(data) # Filling calculators with data
     chunkSize = 100 # resetting the chunkSize to its default value for the next sequence or data-chunk
 
 print('A total number of {} images are processed.'.format(images))
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 offsetMap = noiseCal.getOffset() # Producing offset map
 noiseMap = noiseCal.get() # Producing noise map
 noiseCal.reset() # Resetting noise calculator
 print("Initial maps are created.")
 ```
 
 %% Cell type:markdown id: tags:
 
 ### Offset and Noise Maps prior to Common Mode Correction
 
 In the following, the histogram of the FastCCD offset, FastCCD offset map, as well as the initial uncorrected noise map are plotted:
 
 %% Cell type:code id: tags:
 
 ``` python
 #************** OFFSET MAP HISTOGRAM ***********#
 ho, co = np.histogram(offsetMap.flatten(), bins=700) # ho = offset histogram; co = offset bin centers
 do = {'x': co[:-1],
      'y': ho,
      'y_err': np.sqrt(ho[:]),
      'drawstyle': 'bars',
      'color': 'cornflowerblue',
      'label': 'Raw Signal (ADU)'
      }
 fig = xana.simplePlot(do, figsize='1col', aspect=1, x_label = 'Raw Signal (ADU)', y_label="Counts",
                       x_range = (3400, 4400), title = 'Offset Histogram')
 #fig.savefig('Offset_Hist.svg', format='svg', dpi=1200, bbox_inches='tight')
 
 t0 = PrettyTable()
 t0.title = "Raw Signal"
 t0.field_names = ["Mean", "Median", "Standard Deviation"]
 t0.add_row(["{:0.3f} (ADU)".format(np.mean(data)), "{:0.3f} (ADU)".format(np.median(data)), "{:0.3f} (ADU)"
             .format(np.std(data))])
 print(t0,'\n')
 
 #************** OffsetMAP *******************#
 fig = xana.heatmapPlot(offsetMap[:,:,0], x_label='Column Number', y_label='Row Number',  aspect=1,
                        x_range=(0, y), y_range=(0, x), vmin=3000, vmax=4300, lut_label='Offset (ADU)',
                        panel_x_label='Columns Stat (ADU)', panel_y_label='Rows Stat (ADU)',
                        panel_top_low_lim = 3000, panel_top_high_lim = 4500, panel_side_low_lim = 3000,
                        panel_side_high_lim = 5000, title = 'OffsetMap')
 #fig.savefig('RawOffsetMap.pdf', format='pdf', dpi=400, bbox_inches='tight')
 
 #************** Raw NoiseMAP *******************#
 fig = xana.heatmapPlot(noiseMap[:,:,0], x_label='Column Number', y_label='Row Number', aspect=1,
                        lut_label='Uncorrected Noise (ADU)', x_range=(0, y),
                        y_range=(0, x), vmax=2*np.mean(noiseMap), panel_x_label='Columns Stat (ADU)',
                        panel_y_label='Rows Stat (ADU)', panel_top_low_lim = 0, panel_top_high_lim = 20,
                        panel_side_low_lim = 0, panel_side_high_lim = 50, title = 'Uncorrected NoiseMap')
 #fig.savefig('RawNoiseMap.pdf', format='pdf', dpi=400, bbox_inches='tight')
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # Offset Correction:
 
 offsetCorrection = xcal.OffsetCorrection(sensorSize, offsetMap, nCells = memoryCells, cores=cpuCores, gains=None,
                                          runParallel=run_parallel, blockSize=blockSize)
 
 offsetCorrection.debug()
 
 # Common Mode Correction:
 # This is the new method subtracting the median of all pixels that are read out at the same time along a row:
 cmCorrection = xcal.CommonModeCorrection([data.shape[0], data.shape[1]], [data.shape[0]//2, data.shape[1]],
                                          commonModeAxis, parallel=run_parallel, dType=np.float32, stride=10,
                                          noiseMap=noiseMap.astype(np.float32), minFrac=0)
 
 cmCorrection.debug()
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # Histogram Calculators:
 
 # For offset corrected data:
 histCalCorrected = xcal.HistogramCalculator(sensorSize, bins=400, range=[-200, 200], memoryCells=memoryCells,
                                             cores=cpuCores, gains=None, blockSize=blockSize)
 # For common mode corrected data:
 histCalCMCorrected = xcal.HistogramCalculator(sensorSize, bins=400, range=[-200, 200], memoryCells=memoryCells,
                                               cores=cpuCores, gains=None, blockSize=blockSize)
 ```
 
 %% Cell type:markdown id: tags:
 
 ### Second Iteration
 
 During the second iteration, the data are offset corrected and then common mode corrected to produced a common mode corrected noise map. The common mode correction is calculated by subtracting out the median of all pixels that are read out at the same time along a row.
 
 %% Cell type:code id: tags:
 
 ``` python
 counter1 = 1 # To count how many "if data.shape[2] >= chunkSize" instances are there.
 counter2 = 0 # To count how many "if data.shape[2] < chunkSize" instances are there.
 chunkSize_new = 0 # See below
 
 for data in reader.readChunks():
 
     data = np.bitwise_and(data.astype(np.uint16), 0b0011111111111111).astype(np.float32)
     dx = np.count_nonzero(data, axis=(0, 1))
     data = data[:,:,dx != 0]
     # Some sequences may have less than 500 frames in them. To find out how many images there are, we will temporarily
     # change chunkSize to be the same as whatever number of frames the last chunk of data has:
     if data.shape[2] < chunkSize:
         chunkSize_new = data.shape[2]
         print("Number of images are less than chunkSize. chunkSize is temporarily changed to {}."
               .format(chunkSize_new))
         images = images + chunkSize_new
         counter2 += 1
     else:
         images = counter1*chunkSize + counter2*chunkSize_new
         counter1 += 1
 
     data = offsetCorrection.correct(data) # Offset correction
     offset_corr_data = copy.copy(data) # I am copying this so that I can have access to it in the table below
     histCalCorrected.fill(data)
     cellTable=np.zeros(data.shape[2], np.int32) # Common mode correction
     data = cmCorrection.correct(data.astype(np.float32), cellTable=cellTable) # Common mode correction
     histCalCMCorrected.fill(data)
     noiseCal.fill(data)  # Filling noise calculator with common mode (CM) corrected data
     chunkSize = 100 # resetting the chunkSize to its default value for the next sequence or data-chunk
 
 print('A total number of {} images are processed.'.format(images))
 print("Offset and common mode corrections are applied.")
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 noiseMapCM = noiseCal.get() # Produces CM corrected noise map
 ho, eo, co, so = histCalCorrected.get()
 hCM, eCM, cCM, sCM = histCalCMCorrected.get()
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # We are copying these so that we can replot them later after the calculators are reset:
 
 ho_second_trial = copy.copy(ho)
 co_second_trial = copy.copy(co)
 hCM_second_trial = copy.copy(hCM)
 cCM_second_trial = copy.copy(cCM)
 ```
 
 %% Cell type:markdown id: tags:
 
 ### Signal after Offset and Common Mode Corrections
 
 Here, the offset corrected signal is compared to the common-mode corrected signal (in the form of binned histograms):
 
 %% Cell type:code id: tags:
 
 ``` python
 do = [{'x': co,
      'y': ho,
      'y_err': np.sqrt(ho[:]),
      'drawstyle': 'steps-mid',
      'color': 'cornflowerblue',
      'label': 'Offset Corrected Signal'
      },
      {'x': cCM,
      'y': hCM,
      'y_err': np.sqrt(hCM[:]),
      'drawstyle': 'steps-mid',
      'color': 'red',
      'label': 'Common Mode Corrected Signal'
      }]
 
 fig = xana.simplePlot(do, figsize='2col', aspect=1, x_label = 'Corrected Signal (ADU)', y_label="Counts",
                       x_range = (-20, 20), legend='top-right-frame-1col', title = 'Corrected Signal - 2nd Iteration')
 #fig.savefig('Corrected_Signal_Hist_1.svg', format='svg', dpi=1200, bbox_inches='tight')
 
 t0 = PrettyTable()
 t0.title = "Comparison of the First Round of Corrections - Bad Pixels Not Excluded"
 t0.field_names = ["After Offset Correction", "After Common Mode Correction"]
 t0.add_row(["Mean: {:0.3f} (ADU)".format(np.mean(offset_corr_data)), "Mean: {:0.3f} (ADU)".format(np.mean(data))])
 t0.add_row(["Median: {:0.3f} (ADU)".format(np.median(offset_corr_data)), "Median: {:0.3f} (ADU)"
             .format(np.median(data))])
 t0.add_row(["Standard Deviation: {:0.3f} (ADU)".format(np.std(offset_corr_data)), "Standard Deviation: {:0.3f} (ADU)"
             .format(np.std(data))])
 print(t0,'\n')
 ```
 
 %% Cell type:markdown id: tags:
 
 ### Noise Map after Common Mode Correction
 
 In the following, the effect of common mode correction on the noise is shown. Finally common mode corrected noise map (noiseMapCM) is displayed and compared to the initial uncorrected noise map:
 
 %% Cell type:code id: tags:
 
 ``` python
 #*****NOISE MAP HISTOGRAM FROM THE COMMON MODE CORRECTED DATA*******#
 hn, cn = np.histogram(noiseMap.flatten(), bins=200, range=(0, 40)) # hn: histogram of noise, cn: bin centers for noise
 hn_CM, cn_CM = np.histogram(noiseMapCM.flatten(), bins=200, range=(0, 40))
 
 dn = [{'x': cn[:-1],
      'y': hn,
      #'y_err': np.sqrt(hn[:]),
      'drawstyle': 'steps-mid',#'bars',
      'color': 'blue',#'cornflowerblue',
      'label': 'Uncorrected Noise'
      },
     {'x': cn_CM[:-1],
      'y': hn_CM,
      #'y_err': np.sqrt(hn_CM[:]),
      'drawstyle': 'steps-mid',#'bars',
      'color': 'crimson',#'red',#'cornflowerblue',
      #'ecolor': 'crimson',
      'label': 'Common Mode Corrected Noise'
      }]
 fig = xana.simplePlot(dn, figsize='2col', aspect=1, x_label = 'Noise (ADU)', y_label="Counts",
                       x_range=(0, 40), y_range=(0, 1e6), y_log=True, legend='top-center-frame-1col',
                       title = 'Noise Comparison')
 
 #fig.savefig('Noise_CM_1_Hist.svg', format='svg', dpi=1200, bbox_inches='tight')
 
 fig = xana.heatmapPlot(noiseMapCM[:,:,0], aspect=1, x_label='Column Number', y_label='Row Number',
                        lut_label='Common Mode Corrected Noise (ADU)', x_range=(0,y), y_range=(0,x),
                        vmax=2*np.mean(noiseMapCM), panel_top_low_lim = 0, panel_top_high_lim = 20,
                        panel_side_low_lim = 0, panel_side_high_lim = 50, title = 'Common Mode Corrected Noise',
                        panel_x_label='Columns Stat (ADU)', panel_y_label='Rows Stat (ADU)')
 
 #fig.savefig('NoiseMapCM.pdf', format='pdf', dpi=400, bbox_inches='tight')
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # Resetting the calculators so we can do a third iteration later:
 
 noiseCal.reset()
 histCalCorrected.reset()
 histCalCMCorrected.reset()
 cmCorrection.reset()
 ```
 
 %% Cell type:markdown id: tags:
 
 ### Initial BadPixelMap
 This is generated based on the offset and CM corrected noise maps:
 
 %% Cell type:code id: tags:
 
 ``` python
 bad_pixels = np.zeros(offsetMap.shape, np.uint32)
 mnoffset = np.nanmedian(offsetMap)
 stdoffset = np.nanstd(offsetMap)
 bad_pixels[(offsetMap < mnoffset-bad_pixel_offset_sigma*stdoffset) |
            (offsetMap > mnoffset+bad_pixel_offset_sigma*stdoffset)] = BadPixels.OFFSET_OUT_OF_THRESHOLD.value
 
 mnnoise = np.nanmedian(noiseMapCM)
 stdnoise = np.nanstd(noiseMapCM)
 bad_pixels[(noiseMapCM < mnnoise-bad_pixel_noise_sigma*stdnoise) |
            (noiseMapCM > mnnoise+bad_pixel_noise_sigma*stdnoise)] = BadPixels.NOISE_OUT_OF_THRESHOLD.value
 
 fig = xana.heatmapPlot(np.log2(bad_pixels[:,:,0]),aspect=1, x_label='Column Number', y_label='Row Number',
                        lut_label='2^(Assigned Value to Bad Pixels)', x_range=(0, y), y_range=(0, x),
                        title = 'Bad Pixels Map Excluding Non-Sensitive\n Areas in Middle of CCD',
                        panel_x_label= 'Columns Stat', panel_y_label='Rows Stat')
 ```
 
 %% Cell type:markdown id: tags:
 
 Here, we are adding the pixels in the hole (center of the FastCCD) as well as 4 rows in the center of the detector, which we call overscan region:
 
 %% Cell type:code id: tags:
 
 ``` python
 def create_circular_mask(h, w, center=None, radius=None):
 
     import numpy as np
     import math
 
     if center is None: # use the middle of the image
         center = [int(w/2), int(h/2)]
     if radius is None: # use the smallest distance between the center and image walls
         radius = min(center[0], center[1], w-center[0], h-center[1])
 
     Y, X = np.ogrid[:h, :w]
     dist_from_center = np.sqrt((X - center[0])**2 + (Y-center[1])**2)
 
     mask = dist_from_center < radius
     return mask
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 mask = np.zeros(offsetMap.shape, np.uint32)
 
 
 # Defining a circular mask + a rectangular mask (overscan) for the hole in the middle of the CCD:
 h, w = (x, y)
 hole_mask_bool = create_circular_mask(h-4, w, radius=61.5, center=(w//2, (h-4)//2))
 hole_mask = np.zeros(hole_mask_bool.shape, np.uint32)
 hole_mask[hole_mask_bool] = BadPixels.NON_SENSITIVE.value
 
 overscan_mask = np.full((4, w), BadPixels.OVERSCAN.value)
 
 mask[:,:,0] = np.insert(hole_mask, (h-4)//2, overscan_mask, 0)
 
 # Assigning this masked area as bad pixels:
 bad_pixels = np.bitwise_or(bad_pixels, mask)
 
 fig = xana.heatmapPlot(np.log2(bad_pixels[:,:,0]),aspect=1, x_label='Column Number', y_label='Row Number',
                        lut_label='2^(Assigned Value to Bad Pixels)', x_range=(0, y), y_range=(0, x),
                        panel_top_low_lim = 0, panel_top_high_lim = 20, panel_side_low_lim = 0, panel_side_high_lim = 20,
                        title = 'Bad Pixels Map Including Non-Sensitive\n Areas in Middle of CCD',
                        panel_x_label='Columns Stat', panel_y_label='Rows Stat', vmax=20)
 
 #fig.savefig('BadPixelMap_1.svg', format='svg', dpi=1200, bbox_inches='tight')
 ```
 
 %% Cell type:markdown id: tags:
 
 ### Third Iteration
 
 During the third iteration, the bad pixel map is applied to the data. Bad pixels are masked. Offset and common mode corrections are applied once again to the data, which now have bad pixdels excluded, to produce a common mode corrected noise map:
 
 %% Cell type:code id: tags:
 
 ``` python
 # bad_pixels is an array of (1934, 960, 1) filled with zeros except at indices where we have the actual bad pixels, whose
 # values are set to be: 2 (2^1: BadPixels.OFFSET_OUT_OF_THRESHOLD.value), or
 # 262144 (2^18: BadPixels.OVERSCAN.value), or 524288 (2^19: BadPixels.NON_SENSITIVE.value). These indices can be found
 # using np.argwhere(bad_pixels != 0)
 
 event_threshold = sigmaNoise*np.median(noiseMapCM) # for exclusion of possible cosmic ray events
 noiseCal.setBadPixelMask(bad_pixels != 0)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 counter1 = 1 # To count how many "if data.shape[2] >= chunkSize" instances are there.
 counter2 = 0 # To count how many "if data.shape[2] < chunkSize" instances are there.
 chunkSize_new = 0 # See below
 
 for data in reader.readChunks():
     #data = data.astype(np.float32)
     data = np.bitwise_and(data.astype(np.uint16), 0b0011111111111111).astype(np.float32)
     dx = np.count_nonzero(data, axis=(0, 1))
     data = data[:,:,dx != 0]
     # Some sequences may have less than 500 frames in them. To find out how many images there are, we will temporarily
     # change chunkSize to be the same as whatever number of frames the last chunk of data has:
     if data.shape[2] < chunkSize:
         chunkSize_new = data.shape[2]
         print("Number of images are less than chunkSize. chunkSize is temporarily changed to {}."
               .format(chunkSize_new))
         images = images + chunkSize_new
         counter2 += 1
     else:
         images = counter1*chunkSize + counter2*chunkSize_new
         counter1 += 1
     data_copy = offsetCorrection.correct(copy.copy(data))
     cellTable=np.zeros(data_copy.shape[2], np.int32)
     data_copy = cmCorrection.correct(data_copy.astype(np.float32), cellTable=cellTable)
     data[data_copy > event_threshold] = np.nan # cosmic rays
     data = np.ma.MaskedArray(data, np.isnan(data), fill_value=0) # masking cosmics, the default fill_value is 1e+20
     data = offsetCorrection.correct(data)
     offset_corr_data2 = copy.copy(data) # I am copying this so that I can have access to it in the table below
     histCalCorrected.fill(data)
     cellTable=np.zeros(data.shape[2], np.int32)
     data = cmCorrection.correct(data.astype(np.float32), cellTable=cellTable)
     histCalCMCorrected.fill(data)
     noiseCal.fill(data)
     chunkSize = 100 # resetting the chunkSize to its default value for the next sequence or data-chunk
 
 print('A total number of {} images are processed.'.format(images))
 print("Final iteration is Performed.")
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 noiseMapCM_2nd = noiseCal.get().filled(0) # the masked pixels are filled with zero
 ho2, eo2, co2, so2 = histCalCorrected.get()
 hCM2, eCM2, cCM2, sCM2 = histCalCMCorrected.get()
 ```
 
 %% Cell type:markdown id: tags:
 
 ### Plots of the Final Results
 
 The following plot and table compare the offset and common mode corrected signal with and without the bad pixels:
 
 %% Cell type:code id: tags:
 
 ``` python
 do_Final = [{'x': co_second_trial,
      'y': ho_second_trial,
      'y_err': np.sqrt(ho_second_trial[:]),
      'drawstyle': 'steps-mid',
      'color': 'blue',#'cornflowerblue',
      'errorstyle': 'bars',
      'label': 'Offset Corrected Signal, Bad Pixels Included - 2nd Trial'
      },
     {'x': cCM_second_trial,
      'y': hCM_second_trial,
      'y_err': np.sqrt(hCM_second_trial[:]),
      'drawstyle': 'steps-mid',
      'color': 'red',
      'errorstyle': 'bars',
      'ecolor': 'crimson',
      'label': 'Common Mode Corrected Signal, Bad Pixels Included - 2nd Trial'
      },
     {'x': co2,
      'y': ho2,
      'y_err': np.sqrt(ho2[:]),
      'drawstyle': 'steps-mid',
      'color': 'black', #'cornflowerblue',
      'errorstyle': 'bars',
      'label': 'Offset Corrected Signal, Bad Pixels Excluded - 3rd Trial'
      },
     {'x': cCM2,
      'y': hCM2,
      'y_err': np.sqrt(hCM2[:]),
      'drawstyle': 'steps-mid',
      'color': 'orange', #'cornflowerblue',
      'errorstyle': 'bars',
      'label': 'Common Mode Corrected Signal, Bad Pixels Excluded - 3rd Trial'
      }]
 
 fig = xana.simplePlot(do_Final, figsize='2col', aspect=1, x_label = 'Corrected Signal (ADU)',
                       y_label="Counts (Logarithmic Scale)", y_log=True, x_range=(-40, 40),
                       legend='bottom-left-frame-1col', title = 'Comparison of Corrected Signal')
 #fig.savefig('Corrected_Signal_Hist_2.svg', format='svg', dpi=1200, bbox_inches='tight')
 
 # offset_corr_data2 and data most likely have some nan's => I am going to use nanmean, nanmedian and nanstd functions:
 t0 = PrettyTable()
 t0.title = "Comparison of the Second Round of Corrections - Bad Pixels Excluded"
 t0.field_names = ["After Offset Correction", "After Common Mode Correction"]
 t0.add_row(["Mean: {:0.3f} (ADU)".format(np.nanmean(offset_corr_data2)), "Mean: {:0.3f} (ADU)".format(np.nanmean(data))])
 t0.add_row(["Median: {:0.3f} (ADU)".format(np.nanmedian(offset_corr_data2)), "Median: {:0.3f} (ADU)"
             .format(np.nanmedian(data))])
 t0.add_row(["Standard Deviation: {:0.3f} (ADU)".format(np.nanstd(offset_corr_data2)),
             "Standard Deviation: {:0.3f} (ADU)".format(np.nanstd(data))])
 print(t0,'\n')
 ```
 
 %% Cell type:markdown id: tags:
 
 ### Final NoiseMap
 
 The effect of exclusion of bad pixels on common mode corrected noise is shown below. Finally common mode corrected noise map with bad pixels excluded (noiseMapCM_2nd) is displayed:
 
 %% Cell type:code id: tags:
 
 ``` python
 #*****NOISE MAP HISTOGRAM FROM THE COMMON MODE CORRECTED DATA*******#
 hn_CM2, cn_CM2 = np.histogram(noiseMapCM_2nd.flatten(), bins=200, range=(0, 40))
 
 dn2 = [{'x': cn[:-1],
      'y': hn,
      #'y_err': np.sqrt(hn[:]),
      'drawstyle': 'steps-mid',#'bars',
      'color': 'blue', #'cornflowerblue',
      'label': 'Uncorrected Noise'
      },
     {'x': cn_CM[:-1],
      'y': hn_CM,
      #'y_err': np.sqrt(hn_CM[:]),
      'drawstyle': 'steps-mid',
      'color': 'red',
      #'ecolor': 'crimson',
      'label': 'Common Mode Corrected Noise prior to Bad Pixels Exclusion'
      },
     {'x': cn_CM2[:-1],
      'y': hn_CM2,
      #'y_err': np.sqrt(hn_CM2[:]),
      'drawstyle': 'steps-mid',
      'color': 'black', #'cornflowerblue',
      'label': 'Common Mode Corrected Noise after Bad Pixels Exclusion'
      }]
 
 fig = xana.simplePlot(dn2, figsize='2col', aspect = 1, x_label = 'Noise (ADU)', y_label="Counts", y_log=True,
                       x_range=(0, 40), y_range=(0, 1e6), legend='top-right-frame-1col', title = 'Final Noise Comparison')
 
 #fig.savefig('Noise_Hist_2.svg', format='svg', dpi=1200, bbox_inches='tight')
 
 fig = xana.heatmapPlot(np.log2(noiseMapCM_2nd[:,:,0]), aspect=1, x_label='Column Number', y_label='Row Number',
                        lut_label='Noise (ADU)', x_range=(0, y), y_range=(0, x), vmax=2*np.mean(noiseMapCM_2nd),
                        title = 'Final Common Mode Corrected Noise\n (Bad Pixels Excluded)',
                        panel_x_label='Columns Stat (ADU)', panel_y_label='Rows Stat (ADU)')
 #fig.savefig('NoiseMapCM_2nd.pdf', format='pdf', dpi=400, bbox_inches='tight')
 ```
 
 %% Cell type:markdown id: tags:
 
 ### Final Bad Pixel Map
 
 Lastly, the final bad pixel map is generated based on the OffsetMap and the noiseMapCM_2nd (common mode corrected noise after exclusion of the initial bad pixels):
 
 %% Cell type:code id: tags:
 
 ``` python
 bad_pixels = np.zeros(offsetMap.shape, np.uint32)
 mnoffset = np.nanmedian(offsetMap)
 stdoffset = np.nanstd(offsetMap)
 bad_pixels[(offsetMap < mnoffset-bad_pixel_offset_sigma*stdoffset) |
            (offsetMap > mnoffset+bad_pixel_offset_sigma*stdoffset)] = BadPixels.OFFSET_OUT_OF_THRESHOLD.value
 
 mnnoise = np.nanmedian(noiseMapCM_2nd)
 stdnoise = np.nanstd(noiseMapCM_2nd)
 bad_pixels[(noiseMapCM_2nd < mnnoise-bad_pixel_noise_sigma*stdnoise) |
            (noiseMapCM_2nd > mnnoise+bad_pixel_noise_sigma*stdnoise)] = BadPixels.NOISE_OUT_OF_THRESHOLD.value
 
 bad_pixels = np.bitwise_or(bad_pixels, mask)
 fig = xana.heatmapPlot(np.log2(bad_pixels[:,:,0]),aspect=1, x_label='Column Number', y_label='Row Number',
                        lut_label='2^(Assigned Value to Bad Pixels)', x_range=(0, y), y_range=(0, x),
                        panel_top_low_lim = 0, panel_top_high_lim = 20, panel_side_low_lim = 0, panel_side_high_lim = 20,
                        title = 'Final Bad Pixels Map', panel_x_label='Columns Stat',
                        panel_y_label='Rows Stat', vmax=20)
 #fig.savefig('BadPixelMap_2.svg', format='svg', dpi=1200, bbox_inches='tight')
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 display(Markdown('### Statistics on the Bad Pixels'))
 num_bad_pixels = np.count_nonzero(bad_pixels)
 num_all_pixels = x*y
 percentage_bad_pixels = num_bad_pixels*100/num_all_pixels
 print("Number of bad pixels: {:0.0f}, i.e. {:0.2f}% of all pixels".format(num_bad_pixels, percentage_bad_pixels))
 ```
 
 %% Cell type:markdown id: tags:
 
 ### Electronic Noise
 
 According to Table 6.1 (page 80) of Ivana Klačková's master's thesis: "Conversion gain for the FastCCD high gain is: lower hemisphere = 6.2e-/ADU and upper hemisphere = 6.1e-/ADU." Also, we know that the high gain/medium gain and high gain/low gain ratios are 4 and 8, respectively since high gain = x8, medium gain = x2 and low gain = x1. We do not currently (October - 2019) know the conversion gains for the FastCCD medium and lows gains in electrons. Therefore, we will use those of the high gains (in both hemispheres) together with the gain ratios to convert the noise in ADU to electrons.
 
 The following Tables present the noise along lower hemisphere, upper hemisphere, and the entire FastCCD detector at different stages. Here, the values in the first table (in ADU and e-) are the mean of noise per pixel, where noise is considered to be the initial uncorrected noise, CM corrected noise after second trial (including bad pixels) and CM corrected noise after third trial (excluding bad pixels).
 
 The values of the second table (in electrons) are the standard deviation of noise per pixel.
 
 %% Cell type:code id: tags:
 
 ``` python
 # noiseMap refers to the initial uncorrected noise, noiseMapCM refers to common mode corrected noise with inclusion of
 # bad pixels, and noiseMapCM_2nd refers to common mode corrected noise without inclusion of bad pixels:
 
 ADU_to_electron_hg = (ADU_to_electron_upper_hg + ADU_to_electron_lower_hg)/2 # Average of ADU_to_electron for entire CCD
                                                                              # for high gain
 
 ADU_to_electron_mg = (ADU_to_electron_upper_mg + ADU_to_electron_lower_mg)/2 # Average of ADU_to_electron for entire CCD
                                                                              # for medium gain
 
 
 ADU_to_electron_lg = (ADU_to_electron_upper_lg + ADU_to_electron_lower_lg)/2 # Average of ADU_to_electron for entire CCD
                                                                              # for low gain
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 for gain, value in gain_dict.items():
 
     if det_gain == gain_dict["low gain"]:
         ADU_to_electron = ADU_to_electron_lg
         ADU_to_electron_upper = ADU_to_electron_upper_lg
         ADU_to_electron_lower = ADU_to_electron_lower_lg
 
     elif det_gain == gain_dict["medium gain"]:
         ADU_to_electron = ADU_to_electron_mg
         ADU_to_electron_upper = ADU_to_electron_upper_mg
         ADU_to_electron_lower = ADU_to_electron_lower_mg
 
     else: # Here, we assume the auto gain and high gain conversions from ADU to electrons are the same.
         ADU_to_electron = ADU_to_electron_hg
         ADU_to_electron_upper = ADU_to_electron_upper_hg
         ADU_to_electron_lower = ADU_to_electron_lower_hg
 
 print("Abbreviations:")
 print(" - ED = Entire Detector;\n - LH: Lower Hemisphere;\n - UH: Upper Hemisphere;")
 print(" - CM Noise: Common Mode Corrected Noise;")
 print(" - BP: Bad Pixels\n")
 
 t0 = PrettyTable()
 t0.title = "Averages of Noise per Pixel"
 t0.field_names = ["Uncorrected Noise", "CM Noise, BP Incl.", "CM Noise, BP Excl."]
 t0.add_row(["ED: {:0.2f} ADU = {:0.2f} e-".format(np.mean(noiseMap), np.mean(noiseMap)*ADU_to_electron),
             "ED: {:0.2f} ADU = {:0.2f} e-".format(np.mean(noiseMapCM), np.mean(noiseMapCM)*ADU_to_electron),
             "ED: {:0.2f} ADU = {:0.2f} e-".format(np.mean(noiseMapCM_2nd), np.mean(noiseMapCM_2nd)*ADU_to_electron)])
 t0.add_row(["LH: {:0.2f} ADU = {:0.2f} e-".format(np.mean(noiseMap[:x//2,:]),
                                                   np.mean(noiseMap[:x//2,:])*ADU_to_electron_lower),
             "LH: {:0.2f} ADU = {:0.2f} e-".format(np.mean(noiseMapCM[:x//2,:]),
                                                   np.mean(noiseMapCM[:x//2,:])*ADU_to_electron_lower),
             "LH: {:0.2f} ADU = {:0.2f} e-".format(np.mean(noiseMapCM_2nd[:x//2,:]),
                                                   np.mean(noiseMapCM_2nd[:x//2,:])*ADU_to_electron_lower)])
 t0.add_row(["UH: {:0.2f} ADU = {:0.2f} e-".format(np.mean(noiseMap[x//2:,:]),
                                                   np.mean(noiseMap[x//2:,:])*ADU_to_electron_upper),
             "UH: {:0.2f} ADU = {:0.2f} e-".format(np.mean(noiseMapCM[x//2:,:]),
                                                   np.mean(noiseMapCM[x//2:,:])*ADU_to_electron_upper),
             "UH: {:0.2f} ADU = {:0.2f} e-".format(np.mean(noiseMapCM_2nd[x//2:,:]),
                                                   np.mean(noiseMapCM_2nd[x//2:,:])*ADU_to_electron_upper)])
 print(t0,'\n')
 
 t1 = PrettyTable()
 t1.title = "Standard Deviations of Noise per Pixel"
 t1.field_names = ["Uncorrected Noise", "CM Noise, BP Incl.", "CM Noise, BP Excl."]
 t1.add_row(["ED: {:0.2f} e-".format(np.std(noiseMap)*ADU_to_electron),
             "ED: {:0.2f} e-".format(np.std(noiseMapCM)*ADU_to_electron),
             "ED: {:0.2f} e-".format(np.std(noiseMapCM_2nd)*ADU_to_electron)])
 t1.add_row(["LH: {:0.2f} e-".format(np.std(noiseMap[:x//2,:])*ADU_to_electron_lower),
             "LH: {:0.2f} e-".format(np.std(noiseMapCM[:x//2,:])*ADU_to_electron_lower),
             "LH: {:0.2f} e-".format(np.std(noiseMapCM_2nd[:x//2,:])*ADU_to_electron_lower)])
 t1.add_row(["UH: {:0.2f} e-".format(np.std(noiseMap[x//2:,:])*ADU_to_electron_upper),
             "UH: {:0.2f} e-".format(np.std(noiseMapCM[x//2:,:])*ADU_to_electron_upper),
             "UH: {:0.2f} e-".format(np.std(noiseMapCM_2nd[x//2:,:])*ADU_to_electron_upper)])
 print(t1)
 ```
 
 %% Cell type:markdown id: tags:
 
 ### Calibration Constants
 
 %% Cell type:code id: tags:
 
 ``` python
 dictionary = {}
 dictionary['Offset'] = offsetMap.data
 dictionary['Noise'] = noiseMapCM_2nd.data
 dictionary['BadPixelsDark'] = bad_pixels.data
 md = None
 for const in dictionary:
     dconst = getattr(Constants.CCD(DetectorTypes.fastCCD), const)()
     dconst.data = dictionary[const]
 
     condition = Conditions.Dark.CCD(bias_voltage=bias_voltage,
                                 integration_time=integration_time,
                                 gain_setting=det_gain,
                                 temperature=fix_temperature,
                                 pixels_x=1934,
                                 pixels_y=960)
 
     for parm in condition.parameters:
         if parm.name == "Sensor Temperature":
             parm.lower_deviation = temp_limits
             parm.upper_deviation = temp_limits
 
     device = Detectors.fastCCD1
 
     if db_output:
         md = send_to_db(device, dconst, condition, file_loc,
                         cal_db_interface, creation_time=creation_time, timeout=cal_db_timeout)
 
     if local_output:
         md = save_const_to_h5(device, dconst, condition, dconst.data, file_loc, creation_time, out_folder)
         print(f"Calibration constant {const} is stored locally.")
 
 print("Constants parameter conditions are:\n")
 print(f"• bias_voltage: {bias_voltage}\n• integration_time: {integration_time}\n"
-      f"• temperature: {temperature_k}\n• gain_setting: {gain_setting}\n"
-      f"• in_vacuum: {in_vacuum}\n"
+      f"• temperature: {fix_temperature}\n• gain_setting: {det_gain}\n"
       f"• creation_time: {md.calibration_constant_version.begin_at if md is not None else creation_time}\n")
 ```

notebooks/FastCCD/CorrectionNotebook_NewDAQ_FastCCD_NBC.ipynb

 %% Cell type:markdown id: tags:
 
 # FastCCD Data Correction ##
 
 Authors: I. Klačková, S. Hauf, Version 1.0
 
 The following notebook provides correction of images acquired with the FastCCD.
 
 %% Cell type:code id: tags:
 
 ``` python
 cluster_profile = "noDB" #ipcluster profile to use
 in_folder = "/gpfs/exfel/exp/SCS/201930/p900074/raw" # input folder, required
 out_folder = '/gpfs/exfel/data/scratch/karnem/test/fastccd' # output folder, required
 sequences = [-1] # sequences to correct, set to -1 for all, range allowed
 run = 277 # run number
 
 karabo_da = 'DA05' # data aggregators
 karabo_id = "SCS_CDIDET_FCCD2M" # karabo prefix of PNCCD devices
 receiver_id = "FCCD" # inset for receiver devices
 path_template = 'RAW-R{:04d}-{}-S{{:05d}}.h5' #  path template in hdf5 file
 h5path = '/INSTRUMENT/{}/DAQ/{}:daqOutput/data/image' # path in HDF5 file
 h5path_t = '/CONTROL/{}/CTRL/LSLAN/inputA/crdg/value' # temperature path in HDF5 file
 h5path_cntrl = '/RUN/{}/DET/FCCD'  # path to control data
 
 use_dir_creation_date = True # use dir creation data for calDB queries
 cal_db_interface = "tcp://max-exfl016:8015#8025" # calibration DB interface to use
 cal_db_timeout = 300000000 # timeout on caldb requests
 
 
 cpuCores = 16 #Specifies the number of running cpu cores
 operation_mode = "FF" # FS stands for frame-store and FF for full-frame opeartion
 split_evt_primary_threshold = 7. # primary threshold for split event classification in terms of n sigma noise
 split_evt_secondary_threshold = 4. # secondary threshold for split event classification in terms of n sigma noise
 split_evt_mip_threshold = 1000. # MIP threshold for event classification
 chunk_size_idim = 1 # H5 chunking size of output data
 overwrite = True # overwrite existing files
 sequences_per_node = 1 # sequences to correct per node
 limit_images = 0 # limit images per file
 time_offset_days = 0 # offset in days for calibration parameters
 photon_energy_gain_map = 5.9 # energy in keV
 fix_temperature = 0. # fix temperature to this value, set to 0 to use slow control value
 flipped_between = ["2019-02-01", "2019-04-02"] # detector was flipped during this timespan
 temp_limits = 5 # limits within which temperature is considered the same
 commonModeAxis = 1 # Axis along which common mode will be calculated (0: along rows, 1: along columns)
 
 # Correction Booleans
 only_offset = False # Only apply offset correction
 cti_corr = False # Apply CTI correction
 relgain_corr = True # Apply relative gain correction
 common_mode_corr = False # Apply commonMode correction
 correct_offset_drift = False  # correct for offset drifts
 do_pattern_classification = True # classify split events
 
 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
 # Fill dictionaries comprising bools and arguments for correction and data analysis
 
 # Here the herarichy and dependability for correction booleans are defined
 corr_bools = {}
 
 # offset is at the bottom of AGIPD correction pyramid.
 corr_bools["only_offset"] = only_offset
 
 # Dont apply any corrections if only_offset is requested
 if not only_offset:
 
     # Apply relative gain correction, only if requested
     if relgain_corr:
         corr_bools["relgain"] = relgain_corr
 
         # Apply CTI correction, only if requested
         if cti_corr:
             corr_bools["cti"] = cti_corr
 
     corr_bools["common_mode"] = common_mode_corr
     corr_bools["offset_drift"] = correct_offset_drift
     corr_bools["pattern_class"] = do_pattern_classification
 # Here the herarichy and dependability for data analysis booleans and arguments are defined
 data_analysis_parms = {}
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 import XFELDetAna.xfelprofiler as xprof
 
 profiler = xprof.Profiler()
 profiler.disable()
 from XFELDetAna.util import env
 env.iprofile = cluster_profile
 
 import warnings
 warnings.filterwarnings('ignore')
 
 from XFELDetAna import xfelpycaltools as xcal
 from XFELDetAna import xfelpyanatools as xana
 from XFELDetAna.plotting.util import prettyPlotting
 prettyPlotting=True
 from XFELDetAna.xfelreaders import ChunkReader
 from XFELDetAna.detectors.fastccd import readerh5 as fastccdreaderh5
 
 import numpy as np
 import h5py
 import dateutil.parser
 import matplotlib.pyplot as plt
 from iminuit import Minuit
 
 import time
 import copy
 import os
 
 from prettytable import PrettyTable
 
 from iCalibrationDB import ConstantMetaData, Constants, Conditions, Detectors, Versions
 from iCalibrationDB.detectors import DetectorTypes
 from cal_tools.tools import (get_dir_creation_date,
                              get_random_db_interface,
                              get_constant_from_db_and_time)
 
 from datetime import timedelta
 
 %matplotlib inline
 
 if sequences[0] == -1:
     sequences = None
 
 # select a random port for the data base
 cal_db_interface = get_random_db_interface(cal_db_interface)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 if operation_mode == "FS":
     x = 960 # rows of the FastCCD to analyze in FS mode
     y = 960 # columns of the FastCCD to analyze in FS mode
     print('\nYou are analyzing data in FS mode.')
 else:
     x = 1934 # rows of the FastCCD to analyze in FF mode
     y = 960 # columns of the FastCCD to analyze in FF mode
     print('\nYou are analyzing data in FF mode.')
 
 ped_dir = "{}/r{:04d}".format(in_folder, run)
 fp_name = path_template.format(run, karabo_da)
 fp_path = '{}/{}'.format(ped_dir, fp_name)
 
 h5path = h5path.format(karabo_id, receiver_id)
 h5path_t = h5path_t.format(karabo_id)
 h5path_cntrl = h5path_cntrl.format(karabo_id)
 
 print("Reading data from: {}\n".format(fp_path))
 print("Run is: {}".format(run))
 print("HDF5 path: {}".format(h5path))
 print("Data is output to: {}".format(out_folder))
 
 import datetime
 creation_time = None
 if use_dir_creation_date:
     creation_time = get_dir_creation_date(in_folder, run) + timedelta(days=time_offset_days)
 if creation_time:
     print("Using {} as creation time".format(creation_time.isoformat()))
 
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 
 sensorSize = [x, y]
 chunkSize = 100 #Number of images to read per chunk
 blockSize = [sensorSize[0]//2, sensorSize[1]] #Sensor area will be analysed according to blocksize
 xcal.defaultBlockSize = blockSize
 memoryCells = 1 #FastCCD has 1 memory cell
 #Specifies total number of images to proceed
 
 commonModeBlockSize = blockSize
 # commonModeAxisR = 'row'#Axis along which common mode will be calculated
 run_parallel = True
 profile = False
 
 filename = fp_path.format(sequences[0] if sequences else 0)
 with h5py.File(filename, 'r') as f:
     bias_voltage = int(f['{}/biasclock/bias/value'.format(h5path_cntrl)][0])
     det_gain = int(f['{}/exposure/gain/value'.format(h5path_cntrl)][0])
     integration_time = int(f['{}/exposure/exposure_time/value'.format(h5path_cntrl)][0])
     print("Bias voltage is {} V".format(bias_voltage))
     print("Detector gain is set to x{}".format(det_gain))
     print("Detector integration time is set to {}".format(integration_time))
     temperature = np.mean(f[h5path_t])
     temperature_k = temperature + 273.15
     if fix_temperature != 0.:
         temperature_k = fix_temperature
         print("Using fixed temperature")
     print("Mean temperature was {:0.2f} °C / {:0.2f} K at beginning of run".format(temperature, temperature_k))
 
 
 if not os.path.exists(out_folder):
     os.makedirs(out_folder)
 elif not overwrite:
     # Stop Notebook not only this cell
     raise SystemExit("Output path exists! Exiting")
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 dirlist = sorted(os.listdir(ped_dir))
 file_list = []
 total_sequences = 0
 fsequences = []
 for entry in dirlist:
 
     #only h5 file
     abs_entry = "{}/{}".format(ped_dir, entry)
     if os.path.isfile(abs_entry) and os.path.splitext(abs_entry)[1] == ".h5":
 
         if sequences is None:
             for seq in range(len(dirlist)):
 
                 if path_template.format(run, karabo_da).format(seq) in abs_entry:
                     file_list.append(abs_entry)
                     total_sequences += 1
                     fsequences.append(seq)
         else:
             for seq in sequences:
 
                 if path_template.format(run, karabo_da).format(seq) in abs_entry:
                     file_list.append(os.path.abspath(abs_entry))
                     total_sequences += 1
                     fsequences.append(seq)
 sequences = fsequences
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 import copy
 from IPython.display import HTML, display, Markdown, Latex
 import tabulate
 print("Processing a total of {} sequence files".format(total_sequences))
 table = []
 
 
 for k, f in enumerate(file_list):
     table.append((k, f))
 if len(table):
     md = display(Latex(tabulate.tabulate(table, tablefmt='latex', headers=["#", "file"])))
 ```
 
 %% Cell type:markdown id: tags:
 
 As a first step, dark maps have to be loaded.
 
 %% Cell type:code id: tags:
 
 ``` python
 offsetMap = None
 badPixelMap = None
 noiseMap = None
 
 
 # The following constants are saved in data base as a constant for each gain
 # Hence, a loop over all 3 gains is performed
 
 # This is to be used in messages in reports.
 g_name = {8: "High gain", 2: "Medium gain", 1: "Low gain"}
 
 for i, g in enumerate([8, 2, 1]):
 
     print("Retrieving constants for {}".format(g_name[g]))
 
     # set the operating condition
     condition = Conditions.Dark.CCD(bias_voltage=bias_voltage,
                                     integration_time=integration_time,
                                     gain_setting=g,
                                     temperature=temperature_k,
                                     pixels_x=x,
                                     pixels_y=y)
 
     for parm in condition.parameters:
         if parm.name == "Sensor Temperature":
             parm.lower_deviation = temp_limits
             parm.upper_deviation = temp_limits
 
 
     offset = Constants.CCD(DetectorTypes.fastCCD).Offset()
     noise = Constants.CCD(DetectorTypes.fastCCD).Noise()
     bpix = Constants.CCD(DetectorTypes.fastCCD).BadPixelsDark()
 
     ## retrieve_offset
     offset, offset_time = get_constant_from_db_and_time(device=Detectors.fastCCD1,
                                                         constant=offset,
                                                         condition=condition,
                                                         empty_constant=None,
                                                         cal_db_interface=cal_db_interface,
                                                         creation_time=creation_time,
                                                         timeout=cal_db_timeout, print_once=False)
 
     if offsetMap is None:
         offsetMap = np.zeros(list(offset.shape)+[3], np.float32)
 
     if offset is not None:
         offsetMap[...,i] = offset
     else:
         print("NO OFFSET FOUND IN DB!")
 
     offset_temperature = None
     for parm in condition.parameters:
 
         if parm.name == "Sensor Temperature":
             offset_temperature = parm.value
 
     print("Dark Offset was taken at temperature of {:0.2f}K at {}"
           .format(offset_temperature, offset_time))
 
     ## retrieve_noise
     noise, noise_time = get_constant_from_db_and_time(device=Detectors.fastCCD1,
                                                         constant=noise,
                                                         condition=condition,
                                                         empty_constant=None,
                                                         cal_db_interface=cal_db_interface,
                                                         creation_time=creation_time,
                                                         timeout=cal_db_timeout, print_once=False)
     if noiseMap is None:
         noiseMap = np.zeros(list(noise.shape)+[3], np.float32)
     if noise is not None:
         noiseMap[...,i] = noise
     else:
         print("NO NOISE FOUND IN DB!")
 
     print("Noise at {} was taken at {}"
           .format(g_name[g], noise_time))
 
     ## retrieve_bad pixels
     bpix, bpix_time = get_constant_from_db_and_time(device=Detectors.fastCCD1,
                                                         constant=bpix,
                                                         condition=condition,
                                                         empty_constant=None,
                                                         cal_db_interface=cal_db_interface,
                                                         creation_time=creation_time,
                                                         timeout=cal_db_timeout, print_once=False)
     if badPixelMap is None:
         badPixelMap = np.zeros(list(bpix.shape)+[3], np.uint32)
         print("NO BadPixel FOUND IN DB!")
     if bpix is not None:
         badPixelMap[...,i] = noise
     else:
         print("NO BADPIXEL FOUND IN DB!")
 
     print("BadPixes at {} was taken at {}"
           .format(g_name[g], bpix_time))
 ```
 
 %% Cell type:markdown id: tags:
 
 Loading cti and relative gain values
 
 %% Cell type:code id: tags:
 
 ``` python
 # relative gain
 if corr_bools.get('relgain'):
 
     relgain = Constants.CCD(DetectorTypes.fastCCD).RelativeGain()
 
     # set the operating condition
     condition = Conditions.Illuminated.CCD(bias_voltage=bias_voltage,
                                            integration_time=integration_time,
                                            gain_setting=det_gain,
                                            temperature=temperature_k,
                                            pixels_x=x, pixels_y=y,
                                            photon_energy=photon_energy_gain_map)
 
     relgain, relgain_time = get_constant_from_db_and_time(device=Detectors.fastCCD1,
                                                           constant=relgain,
                                                           condition=condition,
                                                           empty_constant=None,
                                                           cal_db_interface=cal_db_interface,
                                                           creation_time=creation_time,
                                                           timeout=cal_db_timeout, print_once=False)
 
     # TODO: CHECK IF THIS FLIPPING IS CORRECT
     if relgain is None:
         corr_bools["relgain"] = False
         print('Relative Gain was not found. Proceed without Relative Gain correction.')
     else:
         relGain = relgain[::-1, ...]
         print('Relative Gain was taken with creation-date:', relgain_time)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 if corr_bools.get('cti'):
     pass
     ## FASTCCD CTI CURRENTLY IS NOT SUPPORTED BY DET!
 
     # The code is left for tracing past algorithm for later
     # use during the development of the new CTI algorithm.
     # TODO: Proper CTI Retrieval from CTI and the correction
     # in following cells
 
 
 #     relGainCA = copy.copy(relGain)
 #     relGainC = relGainCA[:relGainCA.shape[0]//2,...]
 #     ctiA = np.ones(relGainCA.shape[:2])
 #     cti = np.ones(relGainC.shape[:2])
 #     i = 0
 #     idx = (relGainC[i, :, 0] < 0.9) | (relGainC[i,:,0] > 1.1)
 #     mn1 = np.nanmean(relGainC[i, ~idx, 0])
 
 #     for i in range(1, relGainC.shape[0]):
 #         idx = (relGainC[i, :, 0] < 0.9) | (relGainC[i,:,0] > 1.1)
 #         mn2 = np.nanmean(relGainC[i, ~idx, 0])
 #         cti[i,:] = mn2/mn1
 #     ctiA[:relGainCA.shape[0]//2,...] = cti
 
 #     relGainC = relGainCA[relGainCA.shape[0]//2:,...]
 
 
 #     cti = np.ones(relGainC.shape[:2])
 #     i = -1
 #     idx = (relGainC[i, :, 0] < 0.9) | (relGainC[i,:,0] > 1.1)
 #     mn1 = np.nanmean(relGainC[i, ~idx, 0])
 
 #     for i in range(relGainC.shape[0]-1, 1, -1):
 #         idx = (relGainC[i, :, 0] < 0.9) | (relGainC[i,:,0] > 1.1)
 #         mn2 = np.nanmean(relGainC[i, ~idx, 0])
 #         cti[i,:] = mn2/mn1
 
 #     ctiA[relGainCA.shape[0]//2:,...] = cti
 
 #     relGainCA = copy.copy(relGain)
 #     relGainC = relGainCA[:relGainCA.shape[0]//2,...]
 #     for i in range(relGainC.shape[1]):
 #         idx = (relGainC[:,i, 0] < 0.95) | (relGainC[:,i,0] > 1.05)
 #         relGainC[idx,i,0] = np.nanmean(relGainC[~idx,i,0])
 #         relGainC[idx,i,1] = np.nanmean(relGainC[~idx,i,1])
 #         relGainC[idx,i,2] = np.nanmean(relGainC[~idx,i,2])
 #     relGainCA[:relGainCA.shape[0]//2,...] = relGainC
 #     relGainC = relGainCA[relGainCA.shape[0]//2:,...]
 #     for i in range(relGainC.shape[1]):
 #         idx = (relGainC[:,i, 0] < 0.95) | (relGainC[:,i,0] > 1.05)
 #         relGainC[idx,i,0] = np.nanmean(relGainC[~idx,i,0])
 #         relGainC[idx,i,1] = np.nanmean(relGainC[~idx,i,1])
 #         relGainC[idx,i,2] = np.nanmean(relGainC[~idx,i,2])
 #     relGainCA[relGainCA.shape[0]//2:,...] = relGainC
 #     relGainC = relGainCA*ctiA[...,None]
 
 #     relGain = relGainC
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 if corr_bools.get('relgain'):
     flipped_between = [dateutil.parser.parse(d) for d in flipped_between]
     flip_rgain = creation_time.replace(tzinfo=None) >= flipped_between[0] and creation_time.replace(tzinfo=None) <= flipped_between[1]
     flip_rgain &= (relgain_time.replace(tzinfo=None) >= flipped_between[0]
                 and relgain_time.replace(tzinfo=None) <= flipped_between[1])
     print("Accounting for flipped detector: {}".format(flip_rgain))
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 #************************Calculators************************#
 
 if corr_bools.get('common_mode'):
     cmCorrection = xcal.CommonModeCorrection([x, y],
                                              commonModeBlockSize,
                                              commonModeAxis,
                                              nCells = memoryCells,
                                              stride=10,
                                              runParallel=True,
                                              stats=True, minFrac=0)
 
 patternClassifierLH = xcal.PatternClassifier([x//2, y],
                                               noiseMap[:x//2, :],
                                               split_evt_primary_threshold,
                                              split_evt_secondary_threshold,
                                              split_evt_mip_threshold,
                                               tagFirstSingles = 0,
                                               nCells=memoryCells,
                                               cores=cpuCores,
                                               allowElongated = False,
                                               blockSize=[x//2, y],
                                               runParallel=True)
 
 
 
 patternClassifierUH = xcal.PatternClassifier([x//2, y],
                                               noiseMap[x//2:, :],
                                               split_evt_primary_threshold,
                                               split_evt_secondary_threshold,
                                               split_evt_mip_threshold,
                                               tagFirstSingles = 0,
                                               nCells=memoryCells,
                                               cores=cpuCores,
                                               allowElongated = False,
                                               blockSize=[x//2, y],
                                               runParallel=True)
 
 
 
 
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 #*****************Histogram Calculators******************#
 
 histCalOffsetCor = xcal.HistogramCalculator([x, y],
                                              bins=1050,
                                              range=[-50, 1000],
                                              nCells=memoryCells,
                                              cores=cpuCores,
                                              blockSize=blockSize)
 
 histCalPcorr = xcal.HistogramCalculator([x, y],
                                          bins=1050,
                                          range=[-50, 1000],
                                          nCells=memoryCells,
                                          cores=cpuCores,
                                          blockSize=blockSize)
 
 histCalPcorrS = xcal.HistogramCalculator([x, y],
                                          bins=1050,
                                          range=[-50, 1000],
                                          nCells=memoryCells,
                                          cores=cpuCores,
                                          blockSize=blockSize)
 
 histCalCommonModeCor = xcal.HistogramCalculator([x, y],
                                              bins=1050,
                                              range=[-50, 1000],
                                              nCells=memoryCells,
                                              cores=cpuCores,
                                              blockSize=blockSize)
 ```
 
 %% Cell type:markdown id: tags:
 
 Applying corrections
 
 %% Cell type:code id: tags:
 
 ``` python
 patternClassifierLH._imagesPerChunk = 500
 patternClassifierUH._imagesPerChunk = 500
 patternClassifierLH.debug()
 patternClassifierUH.debug()
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 histCalOffsetCor.debug()
 histCalCommonModeCor.debug()
 histCalPcorr.debug()
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 def copy_and_sanitize_non_cal_data(infile, outfile, h5base):
     """
 
 
     :param infile: Input file
     :param outfile: Otput file
     :param h5base:
     """
 
     if h5base.startswith("/"):
         h5base = h5base[1:]
     dont_copy = ['pixels']
     dont_copy = [h5base+"/{}".format(do)
                 for do in dont_copy]
 
     def visitor(k, item):
         """
 
 
         :param k:
         :param item:
         """
         if k not in dont_copy:
             if isinstance(item, h5py.Group):
                 outfile.create_group(k)
             elif isinstance(item, h5py.Dataset):
                 group = str(k).split("/")
                 group = "/".join(group[:-1])
                 infile.copy(k, outfile[group])
 
     infile.visititems(visitor)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 mean_im = None
 single_im = None
 mean_im_cc = None
 single_im_cc = None
 drift_lh = []
 drift_uh = []
 offsetMap = np.squeeze(offsetMap)
 noiseMap = np.squeeze(noiseMap)
 badPixelMap = np.squeeze(badPixelMap)
 if corr_bools.get('relgain'):
     #TODO: This should be removed after properly injecting gain const for all 3 gains
     if len(relGain.shape) == 3 and relGain.shape[2] == 1:
         # This is a temporary solution for the relGain of one gain in db
         relGain = np.repeat(relGain, 3, axis=2)
     relGain = np.squeeze(relGain)
 
 for k, f in enumerate(file_list):
     with h5py.File(f, 'r', driver='core') as infile:
         out_fileb = "{}/{}".format(out_folder, f.split("/")[-1])
         out_file = out_fileb.replace("RAW", "CORR")
         #out_filed = out_fileb.replace("RAW", "CORR-SC")
 
         data = None
         noise = None
         try:
             with h5py.File(out_file, "w") as ofile:
 
                 copy_and_sanitize_non_cal_data(infile, ofile, h5path)
                 data = infile[h5path+"/pixels"][()]
                 nzidx = np.count_nonzero(data, axis=(1,2))
                 data = data[nzidx != 0, ...]
                 if limit_images > 0:
                     data = data[:limit_images,...]
                 oshape = data.shape
                 data = np.moveaxis(data, 0, 2)
                 ddset = ofile.create_dataset(h5path+"/pixels",
                                              oshape,
                                              chunks=(chunk_size_idim, oshape[1], oshape[2]),
                                              dtype=np.float32)
 
                 ddsetm = ofile.create_dataset(h5path+"/mask",
                                              oshape,
                                              chunks=(chunk_size_idim, oshape[1], oshape[2]),
                                              dtype=np.uint32, compression="gzip")
 
                 ddsetg = ofile.create_dataset(h5path+"/gain",
                                              oshape,
                                              chunks=(chunk_size_idim, oshape[1], oshape[2]),
                                              dtype=np.uint8, compression="gzip")
 
                 # Getting the 14th and 15th bit from the data,
                 # which contains gain before removing them from the data.
                 gain = np.right_shift(data, 14)
 
                 # gains are stored as 00 for High gain,
                 # 01 for Medium gain and 11 for low gain.
                 # Hence, the subtraction from the
                 # gain's int values to have 0, 1 and 2
                 gain[gain != 0] -= 1
 
                 fstride = 1
-                if not flip_rgain:  # rgain was taken during flipped orientation
-                    fstride = -1
+                if corr_bools.get('relgain'):
+                    if not flip_rgain:  # rgain was taken during flipped orientation
+                        fstride = -1
 
                 data = np.bitwise_and(data, 0b0011111111111111).astype(np.float32)
 
                 # creating maps for correction usage
                 omap = np.repeat(offsetMap[...,None,:], data.shape[2], axis=2)
                 nmap = np.repeat(noiseMap[...,None,:], data.shape[2], axis=2)
                 bmap = np.repeat(badPixelMap[...,None,:], data.shape[2], axis=2)
 
                 # selecting element value related to the gain of the pixel.
                 offset = np.choose(gain, (omap[...,0], omap[...,1], omap[...,2]))
                 noise = np.choose(gain, (nmap[...,0], nmap[...,1], nmap[...,2]))
                 bpix = np.choose(gain, (bmap[...,0], bmap[...,1], bmap[...,2]))
 
                 # same for relative gain if correction is available
                 if corr_bools.get('relgain'):
                     rmap = np.repeat(relGain[:,::fstride,None,:], data.shape[2], axis=2)
                     rg = np.choose(gain, (rmap[...,0], rmap[...,1], rmap[...,2]))
 
                 # Apply offset correction
                 data -= offset
 
                 if corr_bools.get('relgain'):
                     # Apply relative gain correction
                     # TODO: check relgain correction in pydetlib
                     # and use it here if the same.
                     data *= rg
 
                 if corr_bools.get("offset_drift"):
                     # Put in consideration the temperature
                     # change of the FastCCD's hole.
                     lhd = np.mean(data[x//2-10:x//2,y//2-5:y//2+5,:], axis=(0,1))
                     data[:x//2, :, :] -= lhd
                     drift_lh.append(lhd)
 
                     uhd = np.mean(data[x//2:x//2+10,y//2-5:y//2+5,:], axis=(0,1))
                     data[x//2:, :, :] -= uhd
                     drift_uh.append(uhd)
 
                 histCalOffsetCor.fill(data)
 
                 ddset[...] = np.moveaxis(data, 2, 0)
                 ddsetm[...] = np.moveaxis(bpix, 2, 0)
                 ddsetg[...] = np.moveaxis(gain, 2, 0).astype(np.uint8)
 
                 if mean_im is None:
                         mean_im = np.nanmean(data, axis=2)
                         single_im = data[...,0]
 
                 if corr_bools.get("pattern_class"):
 
                     ddsetcm = ofile.create_dataset(h5path+"/pixels_cm",
                                          oshape,
                                          chunks=(chunk_size_idim, oshape[1], oshape[2]),
                                          dtype=np.float32)
 
                     ddsetc = ofile.create_dataset(h5path+"/pixels_classified",
                                          oshape,
                                          chunks=(chunk_size_idim, oshape[1], oshape[2]),
                                          dtype=np.float32, compression="gzip")
 
                     ddsetp = ofile.create_dataset(h5path+"/patterns",
                                                  oshape,
                                                  chunks=(chunk_size_idim, oshape[1], oshape[2]),
                                                  dtype=np.int32, compression="gzip")
 
                     # The calculation of the cluster map
                     patternClassifierLH._noisemap = noise[:x//2, :, :]
                     patternClassifierUH._noisemap = noise[x//2:, :, :]
 
                     # common mode correction
                     if corr_bools.get("common_mode"):
                         cellTable = np.zeros(data.shape[2], np.int32) # Common mode correction
                         data = cmCorrection.correct(data.astype(np.float32),
                                                     cellTable=cellTable,
                                                     noiseMap=noise)
 
                         # correct for the row common mode
                         ddsetcm[...] = np.moveaxis(data, 2, 0)
 
                         histCalCommonModeCor.fill(data)
 
                     dataLH = data[:x//2, :, :]
                     dataUH = data[x//2:, :, :]
 
                     dataLH, patternsLH = patternClassifierLH.classify(dataLH)
                     dataUH, patternsUH = patternClassifierUH.classify(dataUH)
 
                     data[:x//2, :, :] = dataLH
                     data[x//2:, :, :] = dataUH
 
                     patterns = np.zeros(data.shape, patternsLH.dtype)
                     patterns[:x//2, :, :] = patternsLH
                     patterns[x//2:, :, :] = patternsUH
 
                     data[data < split_evt_primary_threshold*noise] = 0
                     ddsetc[...] = np.moveaxis(data, 2, 0)
                     ddsetp[...] = np.moveaxis(patterns, 2, 0)
 
                     histCalPcorr.fill(data)
                     data[patterns != 100] = np.nan
                     histCalPcorrS.fill(data)
 
                     if mean_im_cc is None:
                         mean_im_cc = np.nanmean(data, axis=2)
                         single_im_cc = data[...,0]
 
         except Exception as e:
             print("Couldn't calibrate data in {}: {}".format(f, e))
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 if corr_bools.get("offset_drift"):
     lhds = np.concatenate(drift_lh)
     uhds = np.concatenate(drift_uh)
     fig = plt.figure(figsize=(10,5))
     ax = fig.add_subplot(111)
     ax.plot(lhds, label="Lower hem.")
     ax.plot(uhds, label="Upper hem.")
     ax.set_xlabel("Frame #")
     ax.set_xlabel("Offset drift (ADU)")
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 if corr_bools.get("pattern_class"):
     print("******************LOWER HEMISPHERE******************\n")
 
     patternStatsLH = patternClassifierLH.getPatternStats()
     fig = plt.figure(figsize=(15,15))
     ax = fig.add_subplot(4,4,1)
     sfields = ["singles", "first singles", "clusters"]
     mfields = ["doubles", "triples", "quads"]
     relativeOccurances =  []
     labels = []
     for i, f in enumerate(sfields):
         relativeOccurances.append(patternStatsLH[f])
         labels.append(f)
     for i, f in enumerate(mfields):
         for k in range(len(patternStatsLH[f])):
             relativeOccurances.append(patternStatsLH[f][k])
             labels.append(f+"("+str(k)+")")
     relativeOccurances = np.array(relativeOccurances, np.float)
     relativeOccurances/=np.sum(relativeOccurances)
     pie = ax.pie(relativeOccurances, labels=labels, autopct='%1.1f%%', shadow=True)
     ax.set_title("Pattern occurrence")
     # Set aspect ratio to be equal so that pie is drawn as a circle.
     a = ax.axis('equal')
 
     smaps = ["singlemap", "firstsinglemap", "clustermap"]
     for i, m in enumerate(smaps):
 
         ax = fig.add_subplot(4,4,2+i)
 
         pmap = ax.imshow(patternStatsLH[m], interpolation="nearest", vmax=2*np.nanmedian(patternStatsLH[m]))
         ax.set_title(m)
         cb = fig.colorbar(pmap)
 
     mmaps = ["doublemap", "triplemap", "quadmap"]
     k = 0
     for i, m in enumerate(mmaps):
 
         for j in range(4):
             ax = fig.add_subplot(4,4,2+len(smaps)+k)
             pmap = ax.imshow(patternStatsLH[m][j], interpolation="nearest", vmax=2*np.median(patternStatsLH[m][j]))
             ax.set_title(m+"("+str(j)+")")
             cb = fig.colorbar(pmap)
             k+=1
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 if corr_bools.get("pattern_class"):
     print("******************UPPER HEMISPHERE******************\n")
 
     patternStatsUH = patternClassifierUH.getPatternStats()
     fig = plt.figure(figsize=(15,15))
     ax = fig.add_subplot(4,4,1)
     sfields = ["singles", "first singles", "clusters"]
     mfields = ["doubles", "triples", "quads"]
     relativeOccurances =  []
     labels = []
     for i, f in enumerate(sfields):
         relativeOccurances.append(patternStatsUH[f])
         labels.append(f)
     for i, f in enumerate(mfields):
         for k in range(len(patternStatsUH[f])):
             relativeOccurances.append(patternStatsUH[f][k])
             labels.append(f+"("+str(k)+")")
     relativeOccurances = np.array(relativeOccurances, np.float)
     relativeOccurances/=np.sum(relativeOccurances)
     pie = ax.pie(relativeOccurances, labels=labels, autopct='%1.1f%%', shadow=True)
     ax.set_title("Pattern occurrence")
     # Set aspect ratio to be equal so that pie is drawn as a circle.
     a = ax.axis('equal')
 
     smaps = ["singlemap", "firstsinglemap", "clustermap"]
     for i, m in enumerate(smaps):
 
         ax = fig.add_subplot(4,4,2+i)
 
         pmap = ax.imshow(patternStatsUH[m], interpolation="nearest", vmax=2*np.nanmedian(patternStatsUH[m]))
         ax.set_title(m)
         cb = fig.colorbar(pmap)
 
     mmaps = ["doublemap", "triplemap", "quadmap"]
     k = 0
     for i, m in enumerate(mmaps):
 
         for j in range(4):
             ax = fig.add_subplot(4,4,2+len(smaps)+k)
             pmap = ax.imshow(patternStatsUH[m][j], interpolation="nearest", vmax=np.median(patternStatsUH[m][j]))
             ax.set_title(m+"("+str(j)+")")
             cb = fig.colorbar(pmap)
             k+=1
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 if corr_bools.get("pattern_class"):
     t0 = PrettyTable()
     t0.title = "Total number of Counts after all corrections"
     t0.field_names = ["Hemisphere","Singles", "First-singles", "Clusters"]
     t0.add_row(["LH", patternStatsLH['singles'], patternStatsLH['first singles'], patternStatsLH['clusters']])
     t0.add_row(["UH", patternStatsUH['singles'], patternStatsUH['first singles'], patternStatsUH['clusters']])
 
     print(t0)
 
     t1 = PrettyTable()
 
     t1.field_names = ["Index","D-LH", "D-UH", "T-LH", "T-UH", "Q-LH", "Q-UH"]
 
     t1.add_row([0, patternStatsLH['doubles'][0], patternStatsUH['doubles'][0], patternStatsLH['triples'][0], patternStatsUH['triples'][0], patternStatsLH['quads'][0], patternStatsUH['quads'][0]])
     t1.add_row([1, patternStatsLH['doubles'][1], patternStatsUH['doubles'][1], patternStatsLH['triples'][1], patternStatsUH['triples'][1], patternStatsLH['quads'][1], patternStatsUH['quads'][1]])
     t1.add_row([2, patternStatsLH['doubles'][2], patternStatsUH['doubles'][2], patternStatsLH['triples'][2], patternStatsUH['triples'][2], patternStatsLH['quads'][2], patternStatsUH['quads'][2]])
     t1.add_row([3, patternStatsLH['doubles'][3], patternStatsUH['doubles'][3], patternStatsLH['triples'][3], patternStatsUH['triples'][3], patternStatsLH['quads'][3], patternStatsUH['quads'][3]])
 
     print(t1)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 if corr_bools.get("pattern_class"):
     doublesLH = patternStatsLH['doubles'][0] + patternStatsLH['doubles'][1] + patternStatsLH['doubles'][2] + patternStatsLH['doubles'][3]
     triplesLH = patternStatsLH['triples'][0] + patternStatsLH['triples'][1] + patternStatsLH['triples'][2] + patternStatsLH['triples'][3]
     quadsLH = patternStatsLH['quads'][0] + patternStatsLH['quads'][1] + patternStatsLH['quads'][2] + patternStatsLH['quads'][3]
     allsinglesLH = patternStatsLH['singles'] + patternStatsLH['first singles']
     eventsLH = allsinglesLH + doublesLH + triplesLH + quadsLH
 
     doublesUH = patternStatsUH['doubles'][0] + patternStatsUH['doubles'][1] + patternStatsUH['doubles'][2] + patternStatsUH['doubles'][3]
     triplesUH = patternStatsUH['triples'][0] + patternStatsUH['triples'][1] + patternStatsUH['triples'][2] + patternStatsUH['triples'][3]
     quadsUH = patternStatsUH['quads'][0] + patternStatsUH['quads'][1] + patternStatsUH['quads'][2] + patternStatsUH['quads'][3]
     allsinglesUH = patternStatsUH['singles'] + patternStatsUH['first singles']
     eventsUH = allsinglesUH + doublesUH + triplesUH + quadsUH
 
     reloccurLH = np.array([allsinglesLH/eventsLH, doublesLH/eventsLH, triplesLH/eventsLH, quadsLH/eventsLH])
     reloccurUH = np.array([allsinglesUH/eventsUH, doublesUH/eventsUH, triplesUH/eventsUH, quadsUH/eventsUH])
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 if corr_bools.get("pattern_class"):
     fig = plt.figure(figsize=(10,5))
     ax = fig.add_subplot(1,2,1)
     labels = ['singles', 'doubles', 'triples', 'quads']
     pie = ax.pie(reloccurLH, labels=labels, autopct='%1.1f%%', shadow=True)
     ax.set_title("Pattern occurrence LH")
     # Set aspect ratio to be equal so that pie is drawn as a circle.
     a = ax.axis('equal')
     ax = fig.add_subplot(1,2,2)
     pie = ax.pie(reloccurUH, labels=labels, autopct='%1.1f%%', shadow=True)
     ax.set_title("Pattern occurrence UH")
     # Set aspect ratio to be equal so that pie is drawn as a circle.
     a = ax.axis('equal')
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 ho,eo,co,so = histCalOffsetCor.get()
 
 
 d = [{'x': co,
       'y': ho,
       'y_err': np.sqrt(ho[:]),
       'drawstyle': 'steps-mid',
       'errorstyle': 'bars',
       'errorcoarsing': 2,
       'label': 'Offset corr.'
      },
 
      ]
 
 if corr_bools.get("pattern_class") and corr_bools.get("common_mode"):
     hcm,ecm,ccm,scm = histCalCommonModeCor.get()
     d.append({'x': ccm,
               'y': hcm,
               'y_err': np.sqrt(hcm[:]),
               'drawstyle': 'steps-mid',
               'errorstyle': 'bars',
               'errorcoarsing': 2,
               'label': 'CommonMode corr.'
                })
 
 fig = xana.simplePlot(d, aspect=1, x_label='Energy(ADU)',
                       y_label='Number of occurrences', figsize='2col',
                       y_log=True, x_range=(-50,500),
                       legend='top-center-frame-2col')
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 if corr_bools.get("pattern_class"):
     h1,e1L,c1L,s1L = histCalPcorr.get()
     h1s,e1Ls,c1Ls,s1Ls = histCalPcorrS.get()
 
 
     d = [
          {'x': c1L,
           'y': h1,
           'y_err': np.sqrt(h1[:]),
           'drawstyle': 'steps-mid',
           'label': 'Split event corrected'},
          {'x': c1Ls,
           'y': h1s,
           'y_err': np.sqrt(h1s[:]),
           'drawstyle': 'steps-mid',
           'label': 'Single pixel hits'}
     ]
 
 
     fig = xana.simplePlot(d, aspect=1, x_label='Energy(ADU)',
                           y_label='Number of occurrences', figsize='2col',
                           y_log=True, x_range=(0,200),x_log=False,
                           legend='top-center-frame-2col')
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Mean Image of first Sequence ##
 
 %% Cell type:code id: tags:
 
 ``` python
 fig = xana.heatmapPlot(mean_im,
                        x_label='Columns', y_label='Rows',
                        lut_label='Signal (ADU)',
                        x_range=(0,y),
                        y_range=(0,x), vmin=-50, vmax=500)
 
 if corr_bools.get("pattern_class"):
     fig = xana.heatmapPlot(mean_im_cc,
                        x_label='Columns', y_label='Rows',
                        lut_label='Signal (ADU)',
                        x_range=(0,y),
                        y_range=(0,x), vmin=-50, vmax=500)
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Single Shot of first Sequnce ##
 
 %% Cell type:code id: tags:
 
 ``` python
 fig = xana.heatmapPlot(single_im,
                        x_label='Columns', y_label='Rows',
                        lut_label='Signal (ADU)',
                        x_range=(0,y),
                        y_range=(0,x), vmin=-50, vmax=500)
 
 if corr_bools.get("pattern_class"):
     fig = xana.heatmapPlot(single_im_cc,
                        x_label='Columns', y_label='Rows',
                        lut_label='Signal (ADU)',
                        x_range=(0,y),
                        y_range=(0,x), vmin=-50, vmax=500)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 ```

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
 cluster_profile = "noDB" # The ipcluster profile to use
 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
 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
 
 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
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 from collections import OrderedDict
 import copy
 from datetime import datetime
 from functools import partial
 import os
 import warnings
 warnings.filterwarnings('ignore')
 
 import dateutil.parser
 import h5py
 from ipyparallel import Client
 from IPython.display import display, Markdown, Latex
 import matplotlib
 matplotlib.use("agg")
 import matplotlib.patches as patches
 import matplotlib.pyplot as plt
 %matplotlib inline
 import numpy as np
 import tabulate
 
 from cal_tools.enums import BadPixels
 from cal_tools.plotting import (show_overview, plot_badpix_3d,
                                 create_constant_overview,
                                 show_processed_modules)
 from cal_tools.tools import (get_dir_creation_date, get_from_db,
                              get_notebook_name,
                              get_random_db_interface,
                              map_gain_stages, parse_runs,
                              run_prop_seq_from_path,
                              save_const_to_h5, send_to_db)
 from iCalibrationDB import Constants, Conditions, Detectors, Versions
 from XFELDetAna.plotting.heatmap import heatmapPlot
 from XFELDetAna.plotting.simpleplot import simplePlot
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 client = Client(profile=cluster_profile)
 view = client[:]
 view.use_dill()
 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
 # the actual characterization
 def characterize_module(cells, bp_thresh, skip_first_ntrains, ntrains, test_for_normality,
                         h5path, h5path_idx, inp):
     import numpy as np
     import copy
     import h5py
     from cal_tools.enums import BadPixels
     import scipy.stats
 
     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
 
     filename, channel, gg, cap = inp
     thresholds_offset_hard, thresholds_offset_sigma, thresholds_noise_hard, thresholds_noise_sigma = bp_thresh
 
     infile = h5py.File(filename, "r", driver="core")
 
     h5path = h5path.format(channel)
     h5path_idx = h5path_idx.format(channel)
     count = infile[f"{h5path_idx}/count"][()]
     first = infile[f"{h5path_idx}/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(h5path, channel)][first_image:last_image, ...])
     cellid = np.squeeze(np.array(infile["{}/cellId".format(h5path, channel)][first_image:last_image, ...]))
     infile.close()
 
     im, g = splitOffGainLPD(im[:, 0, ...])
     im = im.astype(np.float32)
 
     im = np.rollaxis(im, 2)
     im = np.rollaxis(im, 2, 1)
 
     offset = np.zeros((im.shape[0], im.shape[1], cells))
     noise = np.zeros((im.shape[0], im.shape[1], cells))
     normal_test = np.zeros((im.shape[0], im.shape[1], cells))
     for cc in range(cells):
         idx = cellid == cc
         if np.any(idx):
 
             offset[..., cc] = np.median(im[:, :, idx], axis=2)
             noise[..., cc] = np.std(im[:, :, idx], axis=2)
             if test_for_normality:
                 _, normal_test[..., cc] = scipy.stats.normaltest(
                     im[:, :, idx], axis=2)
 
     # 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
 
 
 offset_g = OrderedDict()
 noise_g = OrderedDict()
 badpix_g = OrderedDict()
 data_g = OrderedDict()
 ntest_g = OrderedDict()
 
 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()
 
     for i in modules:
         qm = "Q{}M{}".format(i//4 + 1, i % 4 + 1)
         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
 
 p = partial(characterize_module, max_cells,
                 (thresholds_offset_hard, thresholds_offset_sigma,
                  thresholds_noise_hard, thresholds_noise_sigma),
                 skip_first_ntrains, ntrains, test_for_normality,
             h5path, h5path_idx)
 
 # Don't remove. Used for Debugging.
 #results = list(map(p, inp))
 results = view.map_sync(p, inp)
 
 for ir, r in enumerate(results):
     offset, noise, i, gg, cap, bp, data, normal = r
     qm = "Q{}M{}".format(i//4 + 1, i % 4 + 1)
     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])
         data_g[cap][qm] = np.full((ntrains, 3), np.nan)
         ntest_g[cap][qm] = np.zeros_like(offset_g[cap][qm])
 
     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
 proposal = list(filter(None, in_folder.strip('/').split('/')))[-2]
 file_loc = 'proposal:{} runs:{} {} {}'.format(proposal, run_low, run_med, run_high)
 ```
 
 %% 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:
     for qm in offset_g[cap].keys():
         for const in clist:
 
             condition = Conditions.Dark.LPD(memory_cells=max_cells, bias_voltage=bias_voltage,
                                             capacitor=cap)
 
             data, mdata = get_from_db(getattr(detinst, qm),
                                       getattr(Constants.LPD, const)(),
                                       condition,
                                       None,
                                       cal_db_interface, creation_time=creation_time,
                                       verbosity=2, timeout=cal_db_timeout)
 
             old_const[const] = data
 
             if mdata is not None and data is not None:
                 time = mdata.calibration_constant_version.begin_at
                 old_mdata[const] = time.isoformat()
                 os.makedirs('{}/old/'.format(out_folder), exist_ok=True)
                 save_const_to_h5(getattr(detinst, qm),
                                  getattr(Constants.LPD, const)(),
                                  condition, data, file_loc, creation_time,
                                  f'{out_folder}/old/')
             else:
                 old_mdata[const] = "Not found"
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 res = OrderedDict()
 for cap in capacitor_settings:
     res[cap] = OrderedDict()
     for i in modules:
         qm = "Q{}M{}".format(i//4+1, i % 4+1)
 
         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]:
         # 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:
             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)
             device = getattr(Detectors.LPD1M1, qm)
 
             if db_output:
                 md = send_to_db(device, dconst, condition, file_loc,
                                 cal_db_interface, creation_time=creation_time, timeout=cal_db_timeout)
 
             if local_output:
                 md = save_const_to_h5(device, dconst, condition,
                                       dconst.data, file_loc, 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"
               f"• capacitor: {cap}\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
 mnames = []
 for i in modules:
     qm = f"Q{modules[0]//4+1}M{modules[0]%4+1}"
     mnames.append(qm)
     show_processed_modules(dinstance=dinstance, constants=None, mnames=mnames, 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 = "Q{}M{}".format(i//4+1, i % 4+1)
         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])]
             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='outside-top', ncol=3,
+                    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 = "Q{}M{}".format(i//4+1, i%4+1)
 
         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
                     data[data == 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(data, 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))
 
                 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
 display(Markdown('The following pre-existing constants are used for comparison: \n'))
 for key in old_mdata:
     display(Markdown('**{}** at {}'.format(key, old_mdata[key])))
 ```
 
 %% 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[const] is not None:
                     ax = fig.add_subplot(122)
 
                     dataold = np.copy(old_const[const][:, :, :, gain])
 
                     label = '$\Delta$ {}'.format(label)
 
                     if const != 'BadPixelsDark':
                         if old_const['BadPixelsDark'] is not None:
                             dataold[old_const['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])
             datau32 = data.astype(np.uint32)
             l_data.append(len(datau32[datau32>0].flatten()))
             for bit in bits:
                 l_data.append(np.count_nonzero(badpix_g[cap][qm][:,:,:,gain].astype(np.uint32) & bit.value))
 
             if old_const['BadPixelsDark'] is not None:
                 dataold = np.copy(old_const['BadPixelsDark'][:, :, :, gain])
                 datau32old = dataold.astype(np.uint32)
                 l_data_old.append(len(datau32old[datau32old>0].flatten()))
                 for bit in bits:
                     l_data_old.append(np.count_nonzero(old_const['BadPixelsDark'][:, :, :, gain].astype(np.uint32) & 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['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[const] is not None and old_const['BadPixelsDark'] is not None :
                 dataold = np.copy(old_const[const])
                 dataold[old_const['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[const] is not None and old_const['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/ePix10K/Characterize_Darks_ePix10K_NBC.ipynb

 %% Cell type:markdown id: tags:
 
 # ePix10K Dark Characterization #
 
 Author: M. Karnevskiy, Version 1.0
 
 The following notebook provides dark image analysis of the ePix10K detector.
 
 Dark characterization evaluates offset and noise of the detector and gives information about bad pixels. Resulting maps are saved as .h5 files for a latter use and injected to the calibration DB.
 
 %% Cell type:code id: tags:
 
 ``` python
 cluster_profile = "noDB" # ipcluster profile to use
 in_folder = '/gpfs/exfel/exp/HED/201922/p002550/raw' # input folder, required
 out_folder = '/gpfs/exfel/data/scratch/ahmedk/test/epix10' # output folder, required
 sequence = 0 # sequence file to use
 run = 55 # which run to read data from, required
 
 karabo_id = "HED_IA1_EPIX10K-1" # karabo karabo_id
 karabo_da = ["EPIX03"]  # data aggregators
 receiver_id = "RECEIVER" # inset for receiver devices
 path_template = 'RAW-R{:04d}-{}-S{{:05d}}.h5' # the template to use to access data
 h5path = '/INSTRUMENT/{}/DET/{}:daqOutput/data/image/pixels' # path in the HDF5 file to images
 h5path_t = '/INSTRUMENT/{}/DET/{}:daqOutput/data/backTemp'  # path to find temperature at
 h5path_cntrl = '/CONTROL/{}/DET'  # path to control data
 
 use_dir_creation_date = True
 cal_db_interface = "tcp://max-exfl016:8020" # calibration DB interface to use
 cal_db_timeout = 300000 # timeout on caldb requests
 db_output = False # Output constants to the calibration database
 local_output = True # output constants locally
 
 temp_limits = 5 # limit for parameter Operational temperature
 number_dark_frames = 0 # number of images to be used, if set to 0 all available images are used
 db_module = 'ePix10K_M43' # detector karabo_id
 bias_voltage = 200 # bias voltage
 in_vacuum = False # detector operated in vacuum
 fix_temperature = 290. # fix temperature to this value
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 import warnings
 warnings.filterwarnings('ignore')
 import os
 
 import h5py
 import matplotlib.pyplot as plt
 %matplotlib inline
 import numpy as np
 
 from cal_tools.tools import (get_dir_creation_date, save_const_to_h5,
                              get_random_db_interface, send_to_db)
 from iCalibrationDB import Constants, Conditions, Detectors, Versions
 from iCalibrationDB.detectors import DetectorTypes
 from XFELDetAna import xfelpyanatools as xana
 from XFELDetAna import xfelpycaltools as xcal
 from XFELDetAna.plotting.util import prettyPlotting
 prettyPlotting = True
 from XFELDetAna.detectors.fastccd import readerh5 as fastccdreaderh5
 from XFELDetAna.util import env
 env.iprofile = cluster_profile
 import XFELDetAna.xfelprofiler as xprof
 profiler = xprof.Profiler()
 profiler.disable()
 from XFELDetAna.xfelreaders import ChunkReader
 
 h5path = h5path.format(karabo_id, receiver_id)
 h5path_t = h5path_t.format(karabo_id, receiver_id)
 h5path_cntrl = h5path_cntrl.format(karabo_id)
 
 def nImagesOrLimit(nImages, limit):
     if limit == 0:
         return nImages
     else:
         return min(nImages, limit)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 proposal = list(filter(None, in_folder.strip('/').split('/')))[-2]
 file_loc = 'proposal:{} runs:{}'.format(proposal, run)
 
 x = 356  # rows of the ePix10K
 y = 384  # columns of the ePix10K
 
 ped_dir = "{}/r{:04d}".format(in_folder, run)
 fp_name = path_template.format(run, karabo_da[0]).format(sequence)
 filename = '{}/{}'.format(ped_dir, fp_name)
 
 print("Reading data from: {}\n".format(filename))
 print("Run number: {}".format(run))
 print("HDF5 path: {}".format(h5path))
 if use_dir_creation_date:
     creation_time = get_dir_creation_date(in_folder, run)
     print("Using {} as creation time".format(creation_time.isoformat()))
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 sensorSize = [x, y]
 chunkSize = 100  #Number of images to read per chunk
 
 #Sensor area will be analysed according to blocksize
 blockSize = [sensorSize[0] // 2, sensorSize[1] // 2]
 xcal.defaultBlockSize = blockSize
 cpuCores = 4  #Specifies the number of running cpu cores
 memoryCells = 1  #No mamery cells
 
 #Specifies total number of images to proceed
 nImages = fastccdreaderh5.getDataSize(filename, h5path)[0]
 nImages = nImagesOrLimit(nImages, number_dark_frames)
 print("\nNumber of dark images to analyze: ", nImages)
 run_parallel = False
 
 with h5py.File(filename, 'r') as f:
     integration_time = int(f['{}/CONTROL/expTime/value'.format(h5path_cntrl)][0])
     gain_setting = int(f['{}/CONTROL/encodedGainSetting/value'.format(h5path_cntrl)][0])
     temperature = np.mean(f[h5path_t])/100.
     temperature_k = temperature + 273.15
     if fix_temperature != 0:
         temperature_k = fix_temperature
         print("Temperature is fixed!")
     print("Bias voltage is {} V".format(bias_voltage))
     print("Gain setting {}".format(gain_setting))
     print("Detector integration time is set to {}".format(integration_time))
     print("Mean temperature was {:0.2f} °C / {:0.2f} K".format(temperature,
                                                                temperature_k))
     print("Operated in vacuum: {} ".format(in_vacuum))
 
 os.makedirs(out_folder, exist_ok=True)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 reader = ChunkReader(filename, fastccdreaderh5.readData,
                      nImages, chunkSize,
                      path=h5path,
                      pixels_x=sensorSize[0],
                      pixels_y=sensorSize[1], )
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 noiseCal = xcal.NoiseCalculator(sensorSize, memoryCells,
                                 cores=cpuCores, blockSize=blockSize,
                                 parallel=run_parallel)
 histCalRaw = xcal.HistogramCalculator(sensorSize, bins=1000,
                                       range=[0, 10000], parallel=False,
                                       memoryCells=memoryCells,
                                       cores=cpuCores, blockSize=blockSize)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 for data in reader.readChunks():
     data = np.bitwise_and(data.astype(np.uint16), 0b0011111111111111).astype(np.float32)
     dx = np.count_nonzero(data, axis=(0, 1))
     data = data[:, :, dx != 0]
     histCalRaw.fill(data)
     noiseCal.fill(data) #Fill calculators with data
 
 constant_maps = {}
 constant_maps['Offset'] = noiseCal.getOffset()  #Produce offset map
 constant_maps['Noise'] = noiseCal.get()  #Produce noise map
 
 noiseCal.reset()  #Reset noise calculator
 print("Initial maps were created")
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 #**************OFFSET MAP HISTOGRAM***********#
 ho, co = np.histogram(constant_maps['Offset'].flatten(), bins=700)
 
 do = {'x': co[:-1],
       'y': ho,
       'y_err': np.sqrt(ho[:]),
       'drawstyle': 'bars',
       'color': 'cornflowerblue',
       }
 
 fig = xana.simplePlot(do, figsize='1col', aspect=2,
                       x_label='Offset (ADU)',
                       y_label="Counts", y_log=True,
                       )
 
 #*****NOISE MAP HISTOGRAM FROM THE OFFSET CORRECTED DATA*******#
 hn, cn = np.histogram(constant_maps['Noise'].flatten(), bins=200)
 
 dn = {'x': cn[:-1],
       'y': hn,
       'y_err': np.sqrt(hn[:]),
       'drawstyle': 'bars',
       'color': 'cornflowerblue',
       }
 
 fig = xana.simplePlot(dn, figsize='1col', aspect=2,
                       x_label='Noise (ADU)',
                       y_label="Counts",
                       y_log=True)
 
 #**************HEAT MAPS*******************#
 fig = xana.heatmapPlot(constant_maps['Offset'][:, :, 0],
                        x_label='Columns', y_label='Rows',
                        lut_label='Offset (ADU)',
                        x_range=(0, y),
                        y_range=(0, x), vmin=1000, vmax=4000)
 
 fig = xana.heatmapPlot(constant_maps['Noise'][:, :, 0],
                        x_label='Columns', y_label='Rows',
                        lut_label='Noise (ADU)',
                        x_range=(0, y),
                        y_range=(0, x), vmax=2 * np.mean(constant_maps['Noise']))
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # Save constants to DB
 dclass="ePix10K"
 md = None
 
 for const_name in constant_maps.keys():
     det = getattr(Constants, dclass)
     const = getattr(det, const_name)()
     const.data = constant_maps[const_name].data
 
     # set the operating condition
     dcond = Conditions.Dark
     condition =  getattr(dcond, dclass)(bias_voltage=bias_voltage,
                                         integration_time=integration_time,
                                         temperature=temperature_k,
                                         in_vacuum=in_vacuum,
                                         gain_setting=gain_setting)
 
     for parm in condition.parameters:
         if parm.name == "Sensor Temperature":
             parm.lower_deviation = temp_limits
             parm.upper_deviation = temp_limits
 
     device = getattr(Detectors, db_module)
 
     if db_output:
-        md = send_to_db(device, const, const.data, condition, file_loc,
+        md = send_to_db(device, const, condition, file_loc,
                         cal_db_interface, creation_time=creation_time, timeout=cal_db_timeout)
 
     if local_output:
         md = save_const_to_h5(device, const, condition,
                               const.data, file_loc, creation_time, out_folder)
         print(f"Calibration constant {const_name} is stored locally.")
 
 print("Constants parameter conditions are:\n")
 print(f"• bias_voltage: {bias_voltage}\n• integration_time: {integration_time}\n"
       f"• temperature: {temperature_k}\n• gain_setting: {gain_setting}\n"
       f"• in_vacuum: {in_vacuum}\n"
       f"• creation_time: {md.calibration_constant_version.begin_at if md is not None else creation_time}\n")
 ```

notebooks/pnCCD/Correct_pnCCD_NBC.ipynb

 %% Cell type:markdown id: tags:
 
 # pnCCD Data Correction #
 
 Authors: DET Group, Modified by Kiana Setoodehnia on March 2020 - Version 2.0
 
 The following notebook provides offset, relative gain, common mode, split events, and pattern classification corrections of images acquired with the pnCCD. This notebook *does not* yet correct for charge transfer inefficiency.
 
 %% Cell type:code id: tags:
 
 ``` python
 in_folder = "/gpfs/exfel/exp/SQS/202022/p002720/raw" # input folder
 out_folder = '/gpfs/exfel/data/scratch/setoodeh' # output folder
 run = 53 # which run to read data from
 sequences = [-1] # sequences to correct, set to -1 for all, range allowed
 
 db_module = "pnCCD_M205_M206"
 karabo_da = 'PNCCD01' # data aggregators
 karabo_da_control = "PNCCD02" # file inset for control data
 karabo_id = "SQS_NQS_PNCCD1MP" # karabo prefix of PNCCD devices
 receiver_id = "PNCCD_FMT-0" # inset for receiver devices
 path_template = 'RAW-R{:04d}-{}-S{{:05d}}.h5' # the template to use to access data
 path_template_seqs = "{}/r{:04d}/*PNCCD01-S*.h5"
 h5path = '/INSTRUMENT/{}/CAL/{}:output/data/' # path to data in the HDF5 file
 h5path_ctrl = '/CONTROL/{}/CTRL/TCTRL'
 
 overwrite = True # keep this as True to not overwrite the output
-use_dir_creation_date = True # required to obtain creation time of the run
+use_dir_creation_date = True # To obtain creation time of the run
 number_dark_frames = 0 # number of images to be used, if set to 0 all available images are used
 chunk_size_idim = 1 # H5 chunking size of output data
 cluster_profile = "noDB" # ipcluster profile to use
 
 cpuCores = 8
 commonModeBlockSize = [512, 512] # size of the detector in pixels for common mode calculations
 commonModeAxis = 0 # axis along which common mode will be calculated, 0 = row, and 1 = column
 split_evt_primary_threshold = 5. # primary threshold for split event classification in terms of n sigma noise
 split_evt_secondary_threshold = 3. # secondary threshold for split event classification in terms of n sigma noise
 split_evt_mip_threshold = 1000. # MIP threshold for event classification
 sequences_per_node = 1
 limit_images = 0
 seq_num = 0  # sequence number for which the last plot at the end of the notebook is plotted
 
 # pnCCD parameters:
 fix_temperature = 0. # fix temperature in K, set to 0. to use value from slow data.
 gain = 0. # the detector's gain setting, It is later read from file and this value is overwritten
 bias_voltage = 0. # the detector's bias voltage. set to 0. to use value from slow data.
 integration_time = 70
 photon_energy = 1.5 # Al fluorescence in keV
 
 cal_db_interface = "tcp://max-exfl016:8015" # calibration DB interface to use
 cal_db_timeout = 300000 # timeout on caldb requests
 creation_time = "" # To overwrite the measured creation_time. Required Format: YYYY-MM-DD HR:MN:SC.ms e.g. 2019-07-04 11:02:41.00
 
 only_offset = False # Only, apply offset.
 common_mode = True # Apply common mode correction
 relgain = True # Apply relative gain correction
 cti = False # Apply charge transfer inefficiency correction (not implemented, yet)
 do_pattern_classification = False # classify split events
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # Here the herarichy and dependability for correction booleans are defined
 corr_bools = {}
 
 corr_bools["only_offset"] = only_offset
 
 # Dont apply any corrections if only_offset is requested.
 if not only_offset:
     corr_bools["relgain"] = relgain
     corr_bools["cti"] = cti
     corr_bools["common_mode"] = common_mode
     corr_bools["pattern_class"] = do_pattern_classification
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 import warnings
 warnings.filterwarnings('ignore')
 import h5py
 import time
 import copy
 import os
 import traceback
 import glob
 import datetime
 from datetime import timedelta
 from prettytable import PrettyTable
 import numpy as np
 import matplotlib.pyplot as plt
 from iminuit import Minuit
 from IPython.display import display, Markdown
 
 import XFELDetAna.xfelprofiler as xprof
 profiler = xprof.Profiler()
 profiler.disable()
 from XFELDetAna.util import env
 env.iprofile = cluster_profile
 from XFELDetAna import xfelpycaltools as xcal
 from XFELDetAna import xfelpyanatools as xana
 from XFELDetAna.plotting.util import prettyPlotting
 prettyPlotting=True
 from XFELDetAna.xfelreaders import ChunkReader
 from XFELDetAna.detectors.fastccd import readerh5 as fastccdreaderh5
 from iCalibrationDB import ConstantMetaData, Constants, Conditions, Detectors, Versions
 from iCalibrationDB.detectors import DetectorTypes
 from cal_tools.tools import get_dir_creation_date, get_random_db_interface, get_constant_from_db_and_time
 
 %matplotlib inline
 
 if sequences[0] == -1:
     sequences = None
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # Calibration Database Settings, and Some Initial Run Parameters & Paths:
 
 display(Markdown('### Initial Settings and Paths'))
 pixels_x = 1024 # rows of pnCCD in pixels
 pixels_y = 1024 # columns of pnCCD in pixels
 print(f"pnCCD size is: {pixels_x}x{pixels_y} pixels.")
 print(f'Calibration database interface selected: {cal_db_interface}')
 
 proposal = list(filter(None, in_folder.strip('/').split('/')))[-2]
 print(f'Proposal: {proposal}, Run: {run}')
 
 # Paths to the data:
 ped_dir = "{}/r{:04d}".format(in_folder, run)
 fp_name = path_template.format(run, karabo_da)
 fp_path = '{}/{}'.format(ped_dir, fp_name)
 h5path = h5path.format(karabo_id, receiver_id)
 print("HDF5 path to data: {}\n".format(h5path))
 
 # Run's creation time:
 if creation_time:
     try:
         creation_time = datetime.datetime.strptime(creation_time, '%Y-%m-%d %H:%M:%S.%f')
     except Exception as e:
         print(f"creation_time value error: {e}."
                "Use same format as YYYY-MM-DD HR:MN:SC.ms e.g. 2019-07-04 11:02:41.00/n")
         creation_time = None
         print("Given creation time wont be used.")
 else:
     creation_time = None
 
 if not creation_time and use_dir_creation_date:
     creation_time = get_dir_creation_date(in_folder, run)
 
 
 print(f"Creation time: {creation_time}")
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # Reading all sequences of the run:
 file_list = []
 total_sequences = 0
 fsequences = []
 
 if sequences is None:
     file_list = glob.glob(fp_path.format(0).replace('00000', '*'))
     file_list = sorted(file_list, key = lambda x: (len (x), x))
     total_sequences = len(file_list)
     fsequences = range(total_sequences)
 else:
     for seq in sequences:
         abs_entry = fp_path.format(seq)
         if os.path.isfile(abs_entry):
             file_list.append(abs_entry)
             total_sequences += 1
             fsequences.append(seq)
 
 sequences = fsequences
 print(f"This run has a total number of {total_sequences} sequences.\n")
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # extract slow data
 if karabo_da_control:
     ctrl_fname = os.path.join(ped_dir, path_template.format(run, karabo_da_control)).format(sequences[0])
     ctrl_path = h5path_ctrl.format(karabo_id)
     mdl_ctrl_path = f"/CONTROL/{karabo_id}/MDL/"
     try:
         with h5py.File(ctrl_fname, "r") as f:
             if bias_voltage == 0.:
                 bias_voltage = abs(f[os.path.join(mdl_ctrl_path, "DAQ_MPOD/u0voltage/value")][0])
             gain = f[os.path.join(mdl_ctrl_path, "DAQ_GAIN/pNCCDGain/value")][0]
             if fix_temperature == 0.:
                 fix_temperature = f[os.path.join(ctrl_path, "inputA/krdg/value")][0]
     except KeyError:
         print("Error !!! during extracting slow data")
         traceback.print_exc(limit=1)
         print("Control file name:", ctrl_fname)
         print("bias voltage control h5path:", os.path.join(mdl_ctrl_path, "DAQ_MPOD/u0voltage/value"))
         print("gain control h5path:", os.path.join(mdl_ctrl_path, "DAQ_GAIN/pNCCDGain/value"))
         print("fix_temperature control h5path:", os.path.join(ctrl_path, "inputA/krdg/value"))
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # Printing the Parameters Read from the Data File:
 
 display(Markdown('### Detector Parameters'))
 print("Bias voltage is {} V.".format(bias_voltage))
 print("Detector gain is set to {}.".format(gain))
 print("Detector integration time is set to {} ms".format(integration_time))
 print(f"Using a fixed temperature of {fix_temperature} K")
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 display(Markdown('### List of Files to be Processed'))
 print("Reading data from the following files:")
 for i in range(len(file_list)):
     print(file_list[i])
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # Sensor size and block size definitions (important for common mode and other corrections):
 
 sensorSize = [pixels_x, pixels_y]
 blockSize = [sensorSize[0]//2, sensorSize[1]//2] # sensor area will be analysed according to blockSize
 xcal.defaultBlockSize = blockSize # for xcal.HistogramCalculators
 memoryCells = 1 # pnCCD has 1 memory cell
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # Output Folder Creation:
 os.makedirs(out_folder, exist_ok=True)
 if not overwrite:
     raise AttributeError("Output path exists! Exiting")
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # Each xcal.HistogramCalculator requires a total number of bins and a binning range. We define these using a
 # dictionary:
 
 # For all xcal histograms:
 if gain == 1:
     Hist_Bin_Dict = {
         "bins": 35000, # number of bins
         "bin_range": [0, 35000]
     }
 
     # For the numpy histograms on the last cell of the notebook:
     Event_Bin_Dict = {
         "event_bins": 1000, # number of bins
         "b_range": [0, 35000] # bin range
     }
 
 #TODO: make it more adaptive for more than only 2 gains [below was for gain==64 only]
 else:
     # For all xcal histograms:
     Hist_Bin_Dict = {
         "bins": 25000, # number of bins
         "bin_range": [0, 25000]
     }
     # For the numpy histograms on the last cell of the notebook:
     Event_Bin_Dict = {
         "event_bins": 1000, # number of bins
         "b_range": [0, 5000] # bin range
     }
 
 bins = Hist_Bin_Dict["bins"]
 bin_range = Hist_Bin_Dict["bin_range"]
 event_bins = Event_Bin_Dict["event_bins"]
 b_range = Event_Bin_Dict["b_range"]
 
 # On the singles spectrum (uploaded in the middle of this notebook), the ADU values correspoding to the boundaries
 # of the first peak region are used as cti_limit_low and cti_limit_high:
 
 if gain == 1:
     cti_limit_low = 1000 # lower limit of cti
     cti_limit_high = 100000 # higher limit of cti
 #TODO: make it more adaptive for more than only 2 gains [below was for gain==64 only
 else:
     cti_limit_low = 50
     cti_limit_high = 2000
 ```
 
 %% Cell type:markdown id: tags:
 
 As a first step, dark constants have to be retrieved from the calibration database
 
 %% Cell type:code id: tags:
 
 ``` python
 def get_dark(db_parms, bias_voltage, integration_time, gain, fix_temperature, creation_time, run):
 # This function is to retrieve the dark constants associated with the run of interest:
 
     cal_db_interface, cal_db_timeout = db_parms
     cal_db_interface = get_random_db_interface(cal_db_interface)
     print(f'Calibration database interface: {cal_db_interface}')
 
     try:
         Offset = None
         Noise = None
         BadPixelsDark = None
 
         constants = {}
         constants['Offset'] = Offset
         constants['Noise'] = Noise
         constants['BadPixelsDark'] = BadPixelsDark
         when = {}
 
         condition = Conditions.Dark.CCD(bias_voltage=bias_voltage,
                             integration_time=integration_time,
                             gain_setting=gain,
                             temperature=fix_temperature,
                             pixels_x=pixels_x,
                             pixels_y=pixels_y)
 
         for const in constants.keys():
             constants[const], when[const] = \
                 get_constant_from_db_and_time(getattr(Detectors, db_module),
                                       getattr(Constants.CCD(DetectorTypes.pnCCD), const)(),
                                       condition,
                                       np.zeros((pixels_x, pixels_y, 1)),
                                       cal_db_interface,
                                       creation_time=creation_time)
         for parm in condition.parameters:
             if parm.name == "Sensor Temperature":
                 print(f"For run {run}: dark images for offset calculation "
                       f"were taken at temperature: {parm.value:0.2f} "
                       f"K @ {when['Offset']}")
     except Exception as e:
         print(f"Error: {e}")
 
 
     return constants
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 display(Markdown('### Dark Data Retrieval'))
 
 db_parms = cal_db_interface, cal_db_timeout
 
 constants = get_dark(db_parms, bias_voltage, integration_time, gain, fix_temperature, creation_time, run)
 
 fig = xana.heatmapPlot(constants["Offset"][:,:,0], x_label='Columns', y_label='Rows', lut_label='Pedestal (ADU)', aspect=1,
                        x_range=(0, pixels_y), y_range=(0, pixels_x), vmax=16000,
                        panel_x_label='Row Stat (ADU)', panel_y_label='Column Stat (ADU)',
                        title = 'Dark Offset Map')
 
 fig = xana.heatmapPlot(constants["Noise"][:,:,0], x_label='Columns', y_label='Rows',
                        lut_label='Common Mode Corrected Noise (ADU)',
                        aspect=1, x_range=(0, pixels_y), y_range=(0, pixels_x), # vmin=-50, vmax=70000,
                        panel_x_label='Row Stat (ADU)', panel_y_label='Column Stat (ADU)',
                        title = 'Dark Noise Map')
 
 fig = xana.heatmapPlot(np.log2(constants["BadPixelsDark"][:,:,0]), x_label='Columns', y_label='Rows',
                        lut_label='Bad Pixel Value (ADU)',
                        aspect=1, x_range=(0, pixels_y), y_range=(0, pixels_x), # vmin=-50, vmax=70000,
                        panel_x_label='Row Stat (ADU)', panel_y_label='Column Stat (ADU)',
                        title = 'Dark Bad Pixels Map')
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 if corr_bools.get('relgain'):
     display(Markdown('We will now retrieve the relative gain map from the calibration database'))
     metadata = ConstantMetaData()
     relgain = Constants.CCD(DetectorTypes.pnCCD).RelativeGain()
     metadata.calibration_constant = relgain
     # set the operating condition
     condition = Conditions.Illuminated.CCD(bias_voltage=bias_voltage,
                                            integration_time=integration_time,
                                            gain_setting=gain,
                                            temperature=fix_temperature,
                                            pixels_x=pixels_x,
                                            pixels_y=pixels_y,
                                            photon_energy=photon_energy)
 
     constants["RelativeGain"], relgain_time = \
         get_constant_from_db_and_time(getattr(Detectors, db_module),
                                       Constants.CCD(DetectorTypes.pnCCD).RelativeGain(),
                                       condition,
                                       np.zeros((pixels_x, pixels_y)),
                                       cal_db_interface,
                                       creation_time=creation_time)
 
     print(f"Retrieved RelativeGain constant with creation time: {relgain_time}")
 
     display(Markdown('### Relative Gain Map Retrieval'))
     fig = xana.heatmapPlot(constants["RelativeGain"], figsize=(8, 8), x_label='Columns', y_label='Rows', lut_label='Relative Gain',
                            aspect=1, x_range=(0, pixels_y), y_range=(0, pixels_x), vmin=0.8, vmax=1.2,
                            panel_x_label='Along Rows', panel_y_label='Along Columns',
                            panel_top_low_lim = 0.5, panel_top_high_lim = 1.5, panel_side_low_lim = 0.5,
                            panel_side_high_lim = 1.5,
                            title = f'1st Relative Gain Map for pnCCD (Gain = {gain})')
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 #************************ Calculators ************************#
 
 if corr_bools.get('common_mode'):
     # Common Mode Correction Calculator:
     cmCorrection = xcal.CommonModeCorrection([pixels_x, pixels_y],
                                              commonModeBlockSize,
                                              commonModeAxis,
                                              parallel=False, dType=np.float32, stride=1,
                                              noiseMap=constants["Noise"].astype(np.float32), minFrac=0)
     cmCorrection.debug()
 
 if corr_bools.get('pattern_class'):
     # Pattern Classifier Calculator:
     # Left Hemisphere:
     patternClassifierLH = xcal.PatternClassifier([pixels_x, pixels_y//2],
                                                  constants["Noise"][:, :pixels_y//2],
                                                  split_evt_primary_threshold,
                                                  split_evt_secondary_threshold,
                                                  split_evt_mip_threshold,
                                                  tagFirstSingles=3, # track along y-axis, left to right (see
                                                  nCells=memoryCells, # split_event.py file in pydetlib/lib/src/
                                                  cores=cpuCores,     # XFELDetAna/algorithms)
                                                  allowElongated=False,
                                                  blockSize=[pixels_x, pixels_y//2],
                                                  runParallel=True)
 
     # Right Hemisphere:
     patternClassifierRH = xcal.PatternClassifier([pixels_x, pixels_y//2],
                                                  constants["Noise"][:, pixels_y//2:],
                                                  split_evt_primary_threshold,
                                                  split_evt_secondary_threshold,
                                                  split_evt_mip_threshold,
                                                  tagFirstSingles=4, # track along y-axis, right to left
                                                  nCells=memoryCells,
                                                  cores=cpuCores,
                                                  allowElongated=False,
                                                  blockSize=[pixels_x, pixels_y//2],
                                                  runParallel=True)
 
     patternClassifierLH._imagesPerChunk = 500
     patternClassifierRH._imagesPerChunk = 500
     patternClassifierLH.debug()
     patternClassifierRH.debug()
 
     # Setting bad pixels:
     patternClassifierLH.setBadPixelMask(constants["BadPixelsDark"][:, :pixels_y//2] != 0)
     patternClassifierRH.setBadPixelMask(constants["BadPixelsDark"][:, pixels_y//2:] != 0)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 #***************** Histogram Calculators ******************#
 # Will contain uncorrected data:
 histCalRaw = xcal.HistogramCalculator(sensorSize,
                                       bins=bins,
                                       range=bin_range,
                                       nCells=memoryCells,
                                       cores=cpuCores,
                                       blockSize=blockSize)
 histCalRaw.debug()
 # Will contain offset corrected data:
 histCalOffsetCor = xcal.HistogramCalculator(sensorSize,
                                             bins=bins,
                                             range=bin_range,
                                             nCells=memoryCells,
                                             cores=cpuCores,
                                             blockSize=blockSize)
 histCalOffsetCor.debug()
 if corr_bools.get('common_mode'):
     # Will contain common mode corrected data:
     histCalCommonModeCor = xcal.HistogramCalculator(sensorSize,
                                                     bins=bins,
                                                     range=bin_range,
                                                     nCells=memoryCells,
                                                     cores=cpuCores,
                                                     blockSize=blockSize)
     histCalCommonModeCor.debug()
 
 if corr_bools.get('pattern_class'):
 # Will contain split events pattern data:
     histCalPcorr = xcal.HistogramCalculator(sensorSize,
                                             bins=bins,
                                             range=bin_range,
                                             nCells=memoryCells,
                                             cores=cpuCores,
                                             blockSize=blockSize)
     histCalPcorr.debug()
 # Will contain singles events data:
     histCalPcorrS = xcal.HistogramCalculator(sensorSize,
                                              bins=bins,
                                              range=bin_range,
                                              nCells=memoryCells,
                                              cores=cpuCores,
                                              blockSize=blockSize)
     histCalPcorrS.debug()
 if corr_bools.get('relgain'):
     # Will contain gain corrected data:
     histCalGainCor = xcal.HistogramCalculator(sensorSize,
                                               bins=bins,
                                               range=bin_range,
                                               nCells=memoryCells,
                                               cores=cpuCores,
                                               blockSize=blockSize)
     histCalGainCor.debug()
 ```
 
 %% Cell type:markdown id: tags:
 
 Applying offset and common mode corrections to the raw data
 
 %% Cell type:code id: tags:
 
 ``` python
 def copy_and_sanitize_non_cal_data(infile, outfile, h5base):
     '''This function reads the .h5 data and writes the corrected .h5 data.'''
     if h5base.startswith("/"):
         h5base = h5base[1:]
     dont_copy = ['image']
     dont_copy = [h5base+"/{}".format(do) for do in dont_copy]
 
     def visitor(k, item):
         if k not in dont_copy:
             if isinstance(item, h5py.Group):
                 outfile.create_group(k)
             elif isinstance(item, h5py.Dataset):
                 group = str(k).split("/")
                 group = "/".join(group[:-1])
                 infile.copy(k, outfile[group])
 
     infile.visititems(visitor)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # Data corrections and event classifications happen here. Also, the corrected data are written to datasets:
 uncor_mean_im = None
 uncor_single_im = None
 offset_mean_im = None
 offset_single_im = None
 cm_mean_im = None
 cm_single_im = None
 gain_mean_im = None
 mean_im_cc = None
 single_im_cc = None
 
 offsetMap = np.squeeze(constants["Offset"])
 noiseMap = np.squeeze(constants["Noise"])
 badPixelMap = np.squeeze(constants["BadPixelsDark"])
 if corr_bools.get('relgain'):
     relGain = constants["RelativeGain"]
 for k, f in enumerate(file_list):
     with h5py.File(f, 'r', driver='core') as infile:
         out_fileb = "{}/{}".format(out_folder, f.split("/")[-1])
         out_file = out_fileb.replace("RAW", "CORR")
 
         data = None
         noise = None
 
         try:
             with h5py.File(out_file, "w") as ofile:
 
                 copy_and_sanitize_non_cal_data(infile, ofile, h5path)
                 data = infile[h5path+"/image"][()]
                 nzidx = np.count_nonzero(data, axis=(1, 2))
                 data = data[nzidx != 0, ...]
 
                 if limit_images > 0:
                     data = data[:limit_images,...]
                 oshape = data.shape
                 data = np.moveaxis(data, 0, 2)
 
                 ddset = ofile.create_dataset(h5path+"/pixels",
                                              oshape,
                                              chunks=(chunk_size_idim, oshape[1], oshape[2]),
                                              dtype=np.float32)
 
                 ddsetm = ofile.create_dataset(h5path+"/mask",
                                              oshape,
                                              chunks=(chunk_size_idim, oshape[1], oshape[2]),
                                              dtype=np.uint32, compression="gzip")
 
                 if corr_bools.get('relgain'):
                     ddsetg = ofile.create_dataset(h5path+"/gain",
                                                   oshape,
                                                   chunks=(chunk_size_idim, oshape[1], oshape[2]),
                                                   dtype=np.float32, compression="gzip")
 
                 data = data.astype(np.float32)
 
                 # Creating maps for correction usage:
                 offset = np.repeat(offsetMap[...,None], data.shape[2], axis=2)
                 noise = np.repeat(noiseMap[...,None], data.shape[2], axis=2)
                 bpix = np.repeat(badPixelMap[...,None], data.shape[2], axis=2)
                 if corr_bools.get('relgain'):
                     rg = np.repeat(relGain[:,:,None], data.shape[2], axis=2) # rg = relative gain
 
                 # uncor_mean_im = averaged over all non-corrected images (in the first sequence only)
                 if uncor_mean_im is None:
                         uncor_mean_im = np.nanmean(data, axis=2)
                         uncor_single_im = data[...,0] # The non-corrected image corresponding to the first frame
 
                 histCalRaw.fill(data) # filling histogram with raw uncorrected data
 
                 # masking data for bad pixels and equating the values to np.nan so that the pattern classifier
                 # ignores them:
                 data[bpix != 0] = np.nan
 
                 data -= offset # offset correction
                 histCalOffsetCor.fill(data) # filling histogram with offset corrected data
 
                 ddset[...] = np.moveaxis(data, 2, 0)
                 ddsetm[...] = np.moveaxis(bpix, 2, 0)
                 ofile.flush()
 
                 # mean_image = averaged over all offset corrected images (in the first sequence only)
                 if offset_mean_im is None:
                         offset_mean_im = np.nanmean(data, axis=2)
                         offset_single_im = data[...,0] # The offset corrected image corresponding to the first frame
 
                 # cm: common mode, c: classifications, p: even patterns
                 if corr_bools.get('common_mode'):
                     ddsetcm = ofile.create_dataset(h5path+"/pixels_cm",
                                                    oshape,
                                                    chunks=(chunk_size_idim, oshape[1], oshape[2]),
                                                    dtype=np.float32)
 
                     data = cmCorrection.correct(data.astype(np.float32),  # common mode correction
                                                     cellTable=np.zeros(data.shape[2], np.int32))
                     histCalCommonModeCor.fill(data) # filling histogram with common mode corrected data
                     # common mode corrected images:
                     if cm_mean_im is None:
                         cm_mean_im = np.nanmean(data, axis=2)
                         cm_single_im = data[...,0] # The common mode corrected image corresponding to the first frame
                     ddsetcm[...] = np.moveaxis(data, 2, 0)
 
                 if corr_bools.get('relgain'):
                     data /= rg  # relative gain correction
                     histCalGainCor.fill(data) # filling histogram with gain corrected data
                     ddsetg[...] = np.moveaxis(rg, 2, 0).astype(np.float32)
 
                     # gain_mean_image = averaged over gain corrected images (in the first sequence only)
                     if gain_mean_im is None:
                             gain_mean_im = np.nanmean(data, axis=2)
                             gain_single_im = data[...,0] # The gain corrected image corresponding to the first frame
 
                 if corr_bools.get('pattern_class'):
 
 
                     ddsetc = ofile.create_dataset(h5path+"/pixels_classified",
                                                   oshape,
                                                   chunks=(chunk_size_idim, oshape[1], oshape[2]),
                                                   dtype=np.float32, compression="gzip")
 
                     ddsetp = ofile.create_dataset(h5path+"/patterns",
                                                  oshape,
                                                  chunks=(chunk_size_idim, oshape[1], oshape[2]),
                                                  dtype=np.int32, compression="gzip")
 
                     # The calculation of the cluster map:
                     patternClassifierLH._noisemap = noise[:, :pixels_x//2,  :]
                     patternClassifierRH._noisemap = noise[:, pixels_x//2:,  :]
 
                     # Dividing the data into left and right hemispheres:
                     dataLH = data[:, :pixels_x//2, :]
                     dataRH = data[:, pixels_x//2:, :]
 
                     dataLH, patternsLH = patternClassifierLH.classify(dataLH) # pattern classification on corrected data
                     dataRH, patternsRH = patternClassifierRH.classify(dataRH)
 
                     data[:, :pixels_x//2, :] = dataLH
                     data[:, pixels_x//2:, :] = dataRH
 
                     patterns = np.zeros(data.shape, patternsLH.dtype)
                     patterns[:, :pixels_x//2, :] = patternsLH
                     patterns[:, pixels_x//2:, :] = patternsRH
 
                     data[data < split_evt_primary_threshold*noise] = 0
                     ddsetc[...] = np.moveaxis(data, 2, 0)
                     ddsetp[...] = np.moveaxis(patterns, 2, 0)
 
                     histCalPcorr.fill(data) # filling histogram with split events corrected data
                     data[patterns != 100] = np.nan # Discard doubles, triples, quadruple, clusters, first singles
                     histCalPcorrS.fill(data) # filling histogram with singles events data
 
                     # mean_im_cc = averaged over all pattern classified corrected images
                     # (in the first sequence only)
                     if mean_im_cc is None:
                         mean_im_cc = np.nanmean(data, axis=2)
                         single_im_cc = data[...,0] # The final corrected image corresponding to the first frame
         except Exception as e:
             print(f"Couldn't calibrate data in {f}: {e}\n")
 
 print("In addition to offset correction, the following corrections were performed:")
 for k, v in corr_bools.items():
     if v:
         print(k)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # Histograming the resulting spectra:
 
 # First _ refers to the bin edges and second _ refers to statistics and we ignore both.
 
 # if you use histCalRaw.get(cumulatative = True) and so on, the cumulatative = True turns the counts array such as
 # RawHistVals and so on into a 1D array instead of keeping the original shape:
 
 RawHistVals, _, RawHistMids, _ = histCalRaw.get()
 off_cor_HistVals, _, off_cor_HistMids, _ = histCalOffsetCor.get()
 if corr_bools.get('common_mode'):
     cm_cor_HistVals, _, cm_HistMids, _ = histCalCommonModeCor.get()
 if corr_bools.get('relgain'):
     gain_cor_HistVals, _, gain_cor_HistMids, _ = histCalGainCor.get()
 if corr_bools.get('pattern_class'):
     split_HistVals, _, split_HistMids, _ = histCalPcorr.get() # split events corrected
     singles_HistVals, _, singles_HistMids, _ = histCalPcorrS.get() # last s in variable names: singles events
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # Saving intermediate data to disk:
 
 np.savez(os.path.join(out_folder, 'Raw_Events.npz'), RawHistMids, RawHistVals)
 np.savez(os.path.join(out_folder, 'Offset_Corrected_Events.npz'), off_cor_HistMids, off_cor_HistVals)
 if corr_bools.get('common_mode'):
     np.savez(os.path.join(out_folder, 'Common_Mode_Corrected_Events.npz'), cm_HistMids, cm_cor_HistVals)
 if corr_bools.get('relgain'):
     np.savez(os.path.join(out_folder, 'Gain_Corrected_Events.npz'), gain_cor_HistMids, gain_cor_HistVals)
 if corr_bools.get('pattern_class'):
     np.savez(os.path.join(out_folder, 'Split_Events.npz'), split_HistMids, split_HistVals)
     np.savez(os.path.join(out_folder, 'Singles_Events.npz'), singles_HistMids, singles_HistVals)
 
 print("Various spectra are saved to disk in the form of histograms. Please check {}".format(out_folder))
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # depending on the gain, the ranges of data are different, so we set a few parameters so that the plots always look
 # good.
 
 if gain == 1:
     x_range = (0, 35000)
 #TODO: make it more adaptive for more than only 2 gains [below was for gain==64 only
 else:
     x_range = (0, 2000)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 display(Markdown('### Raw vs. Corrected Spectra'))
 
 figure = [{'x': RawHistMids,
       'y': RawHistVals,
       'y_err': np.sqrt(RawHistVals[:]),
       'drawstyle': 'steps-post',
       'errorstyle': 'bars',
       'errorcoarsing': 2,
       'label': 'Uncorrected'
      },
      {'x': off_cor_HistMids,
       'y': off_cor_HistVals,
       'y_err': np.sqrt(off_cor_HistVals[:]),
       'drawstyle': 'steps-post',
       'errorstyle': 'bars',
       'errorcoarsing': 2,
       'label': 'Offset Corrected'
      }]
 
 if corr_bools.get('common_mode'):
     figure.append({'x': cm_HistMids,
                    'y': cm_cor_HistVals,
                    'y_err': np.sqrt(cm_cor_HistVals[:]),
                    'drawstyle': 'steps-post',
                    'errorstyle': 'bars',
                    'errorcoarsing': 2,
                    'label': 'Common Mode Corrected'})
 
 if corr_bools.get('relgain'):
     xrange = (cti_limit_low, cti_limit_high)
     figure.append({'x': gain_cor_HistMids,
                    'y': gain_cor_HistVals,
                    'y_err': np.sqrt(gain_cor_HistVals[:]),
                    'drawstyle': 'steps-post',
                    'errorstyle': 'bars',
                    'errorcoarsing': 2,
                    'label': 'Gain Corrected'})
 else:
     xrange = x_range
 
 
 if corr_bools.get('pattern_class'):
     figure.extend([{'x': split_HistMids,
                    'y': split_HistVals,
                    'y_err': np.sqrt(split_HistVals[:]),
                    'drawstyle': 'steps-post',
                    'errorstyle': 'bars',
                    'errorcoarsing': 2,
                    'label': 'Split Events Corrected'
                   },
                   {'x': singles_HistMids,
                    'y': singles_HistVals,
                    'y_err': np.sqrt(singles_HistVals[:]),
                    'drawstyle': 'steps-post',
                    'errorstyle': 'bars',
                    'errorcoarsing': 2,
                    'label': 'Singles Events'
                   }])
 fig = xana.simplePlot(figure, aspect=1, x_label='ADU', y_label='Number of Occurrences', figsize='2col',
                       y_log=True, x_range=x_range, title = '1 ADU per bin is used.',
                       legend='top-right-frame-1col')
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # This function plots pattern statistics:
 
 def classification_plot(patternStats, hemisphere):
 
     print("****************** {} HEMISPHERE ******************\n"
           .format(hemisphere))
     fig = plt.figure(figsize=(15, 15))
     ax = fig.add_subplot(4, 4, 1)
     sfields = ["singles", "first singles", "clusters"]
     mfields = ["doubles", "triples", "quads"]
     relativeOccurances = []
     labels = []
     for i, f in enumerate(sfields):
         relativeOccurances.append(patternStats[f])
         labels.append(f)
     for i, f in enumerate(mfields):
         for k in range(len(patternStats[f])):
             relativeOccurances.append(patternStats[f][k])
             labels.append("{}({})".format(f, k))
     relativeOccurances = np.array(relativeOccurances, np.float)
     relativeOccurances /= np.sum(relativeOccurances)
     pie = ax.pie(relativeOccurances, labels=labels, autopct='%1.1f%%', shadow=True)
     ax.set_title("Pattern Occurrence")
     # Set aspect ratio to be equal so that pie is drawn as a circle.
     a = ax.axis('equal')
 
     smaps = ["singlemap", "firstsinglemap", "clustermap"]
     for i, m in enumerate(smaps):
 
         ax = fig.add_subplot(4, 4, 2+i)
         pmap = ax.imshow(patternStats[m], interpolation="nearest", vmax=2*np.nanmedian(patternStats[m]))
         ax.set_title(m)
         cb = fig.colorbar(pmap)
 
     mmaps = ["doublemap", "triplemap", "quadmap"]
     k = 0
     for i, m in enumerate(mmaps):
 
         for j in range(4):
             ax = fig.add_subplot(4, 4, 2+len(smaps)+k)
             pmap = ax.imshow(patternStats[m][j], interpolation="nearest", vmax=2*np.median(patternStats[m][j]))
             ax.set_title("{}({})".format(m,j))
             cb = fig.colorbar(pmap)
             k+=1
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # The next two cells plot the classification results for left and right hemispheres, respectively:
 display(Markdown('### Classification Results - Plots'))
 if corr_bools.get('pattern_class'):
     patternStatsLH = patternClassifierLH.getPatternStats()
     classification_plot(patternStatsLH, 'Left')
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 if corr_bools.get('pattern_class'):
     patternStatsRH = patternClassifierRH.getPatternStats()
     classification_plot(patternStatsRH, 'Right')
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 display(Markdown('### Classification Results - Tabulated Statistics'))
 
 if corr_bools.get('pattern_class'):
     t0 = PrettyTable()
     t0.title = "Total Number of Counts after All Corrections"
     t0.field_names = ["Hemisphere", "Singles", "First-Singles", "Clusters"]
     t0.add_row(["LH", patternStatsLH['singles'], patternStatsLH['first singles'], patternStatsLH['clusters']])
     t0.add_row(["RH", patternStatsRH['singles'], patternStatsRH['first singles'], patternStatsRH['clusters']])
     print(t0)
 
     print("Abbreviations: D (Doubles), T (Triples), Q (Quadruples), L (Left), R (Right), and H (Hemisphere).")
     t1 = PrettyTable()
     t1.field_names = ["Index", "D-LH", "D-RH", "T-LH", "T-RH", "Q-LH", "Q-RH"]
     t1.add_row([0, patternStatsLH['doubles'][0], patternStatsRH['doubles'][0], patternStatsLH['triples'][0],
                 patternStatsRH['triples'][0], patternStatsLH['quads'][0], patternStatsRH['quads'][0]])
     t1.add_row([1, patternStatsLH['doubles'][1], patternStatsRH['doubles'][1], patternStatsLH['triples'][1],
                 patternStatsRH['triples'][1], patternStatsLH['quads'][1], patternStatsRH['quads'][1]])
     t1.add_row([2, patternStatsLH['doubles'][2], patternStatsRH['doubles'][2], patternStatsLH['triples'][2],
                 patternStatsRH['triples'][2], patternStatsLH['quads'][2], patternStatsRH['quads'][2]])
     t1.add_row([3, patternStatsLH['doubles'][3], patternStatsRH['doubles'][3], patternStatsLH['triples'][3],
                 patternStatsRH['triples'][3], patternStatsLH['quads'][3], patternStatsRH['quads'][3]])
     print(t1)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 if corr_bools.get('pattern_class'):
     doublesLH = patternStatsLH['doubles'][0] + patternStatsLH['doubles'][1] + patternStatsLH['doubles'][2] + \
     patternStatsLH['doubles'][3]
     triplesLH = patternStatsLH['triples'][0] + patternStatsLH['triples'][1] + patternStatsLH['triples'][2] + \
     patternStatsLH['triples'][3]
     quadsLH = patternStatsLH['quads'][0] + patternStatsLH['quads'][1] + patternStatsLH['quads'][2] + \
     patternStatsLH['quads'][3]
     allsinglesLH = patternStatsLH['singles'] + patternStatsLH['first singles']
     eventsLH = allsinglesLH + doublesLH + triplesLH + quadsLH
 
     doublesRH = patternStatsRH['doubles'][0] + patternStatsRH['doubles'][1] + patternStatsRH['doubles'][2] + \
     patternStatsRH['doubles'][3]
     triplesRH = patternStatsRH['triples'][0] + patternStatsRH['triples'][1] + patternStatsRH['triples'][2] + \
     patternStatsRH['triples'][3]
     quadsRH = patternStatsRH['quads'][0] + patternStatsRH['quads'][1] + patternStatsRH['quads'][2] + \
     patternStatsRH['quads'][3]
     allsinglesRH = patternStatsRH['singles'] + patternStatsRH['first singles']
     eventsRH = allsinglesRH + doublesRH + triplesRH + quadsRH
 
     if eventsLH > 0.:
         reloccurLH = np.array([allsinglesLH/eventsLH, doublesLH/eventsLH, triplesLH/eventsLH, quadsLH/eventsLH])
     else:
         reloccurLH = np.array([0]*4)
     if eventsRH > 0.:
         reloccurRH = np.array([allsinglesRH/eventsRH, doublesRH/eventsRH, triplesRH/eventsRH, quadsRH/eventsRH])
     else:
         reloccurRH = np.array([0]*4)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 display(Markdown('### Classification Results - Pie Charts'))
 
 if corr_bools.get('pattern_class'):
     fig = plt.figure(figsize=(12, 7))
     ax = fig.add_subplot(1, 2, 1)
     labels = ['Singles', 'Doubles', 'Triples', 'Quads']
     pie = ax.pie(reloccurLH, labels=labels, autopct='%1.1f%%', shadow=True)
     ax.set_title("Pattern Occurrence in LH")
     # Set aspect ratio to be equal so that pie is drawn as a circle.
     a = ax.axis('equal')
     ax = fig.add_subplot(1, 2, 2)
     pie = ax.pie(reloccurRH, labels=labels, autopct='%1.1f%%', shadow=True)
     ax.set_title("Pattern Occurrence in RH")
     # Set aspect ratio to be equal so that pie is drawn as a circle.
     a = ax.axis('equal')
 ```
 
 %% Cell type:markdown id: tags:
 
 ### Mean Images of the First Sequence ###
 
 %% Cell type:code id: tags:
 
 ``` python
 fig = xana.heatmapPlot(uncor_mean_im, x_label='Columns', y_label='Rows', lut_label='Signal (ADU)', aspect=1,
                        x_range=(0, pixels_y), y_range=(0, pixels_x),
                        panel_x_label='Row Stat (ADU)', panel_y_label='Column Stat (ADU)',
                        title = 'Uncorrected Image Averaged over the First Sequence')
 
 fig = xana.heatmapPlot(offset_mean_im, x_label='Columns', y_label='Rows', lut_label='Signal (ADU)', aspect=1,
                        x_range=(0, pixels_y), y_range=(0, pixels_x),
                        panel_x_label='Row Stat (ADU)', panel_y_label='Column Stat (ADU)',
                        title = 'Offset Corrected Image Averaged over the First Sequence')
 if corr_bools.get('relgain'):
     fig = xana.heatmapPlot(gain_mean_im, x_label='Columns', y_label='Rows',
                            lut_label='Gain Corrected Signal (ADU)', aspect=1,
                            x_range=(0, pixels_y), y_range=(0, pixels_x),
                            panel_x_label='Row Stat (ADU)', panel_y_label='Column Stat (ADU)',
                            title = 'Mean Gain Corrected Image of the First Sequence')
 
 if corr_bools.get('pattern_class') and corr_bools.get('common_mode'):
     fig = xana.heatmapPlot(cm_mean_im, x_label='Columns', y_label='Rows', lut_label='Signal (ADU)', aspect=1,
                            x_range=(0, pixels_y), y_range=(0, pixels_x),
                            panel_x_label='Row Stat (ADU)', panel_y_label='Column Stat (ADU)',
                            title = 'Common Mode Corrected Image Averaged over the First Sequence')
 
     fig = xana.heatmapPlot(mean_im_cc, x_label='Columns', y_label='Rows', lut_label='Signal (ADU)', aspect=1,
                            x_range=(0, pixels_y), y_range=(0, pixels_x),
                            title = 'Image after All Corrections Averaged over the First Sequence')
 ```
 
 %% Cell type:markdown id: tags:
 
 ### Images of the First Frame of the First Sequence ###
 
 %% Cell type:code id: tags:
 
 ``` python
 fig = xana.heatmapPlot(uncor_single_im, x_label='Columns', y_label='Rows', lut_label='Signal (ADU)', aspect=1,
                        x_range=(0, pixels_y), y_range=(0, pixels_x),
                        panel_x_label='Row Stat (ADU)', panel_y_label='Column Stat (ADU)',
                        title = 'Uncorrected Image (First Frame of First Sequence)')
 
 fig = xana.heatmapPlot(offset_single_im, x_label='Columns', y_label='Rows', lut_label='Signal (ADU)', aspect=1,
                        x_range=(0, pixels_y), y_range=(0, pixels_x),
                        panel_x_label='Row Stat (ADU)', panel_y_label='Column Stat (ADU)',
                        title = 'Offset Corrected Image (First Frame of First Sequence)')
 
 if corr_bools.get('pattern_class'):
     if corr_bools.get('common_mode'):
         fig = xana.heatmapPlot(cm_single_im, x_label='Columns', y_label='Rows', lut_label='Signal (ADU)', aspect=1,
                                x_range=(0, pixels_y), y_range=(0, pixels_x),
                                panel_x_label='Row Stat (ADU)', panel_y_label='Column Stat (ADU)',
                                title = 'Common Mode Corrected Image (First Frame of First Sequence)')
 
     fig = xana.heatmapPlot(single_im_cc, x_label='Columns', y_label='Rows', lut_label='Signal (ADU)', aspect=1,
                            x_range=(0, pixels_y), y_range=(0, pixels_x),
                            panel_x_label='Row Stat (ADU)', panel_y_label='Column Stat (ADU)',
                            title = 'Image after All Corrections (First Frame of First Sequence)')
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # Resetting the histogram calculators:
 histCalRaw.reset()
 histCalOffsetCor.reset()
-histCalPcorr.reset()
-histCalPcorrS.reset()
+if corr_bools.get('pattern_class'):
+    histCalPcorr.reset()
+    histCalPcorrS.reset()
 if corr_bools.get('common_mode'):
     histCalCommonModeCor.reset()
 if corr_bools.get('relgain'):
     histCalGainCor.reset()
 ```
 
 %% Cell type:markdown id: tags:
 
 Next, the corrected event patterns are read from the patterns/ dataset created previously and are separated into 4 different categories (singles, doubles, triples and quadruples) using the pattern indices. However, this is done only for one sequence, corresponding to the seq_num variable, as an example.
 
 Note that the number of bins and the bin range for the following histograms may be different from those presented above (depending on gain) to make the counts more noticible and the peaks more defined.
 
 If you are interested in plotting the events from all sequences or the spectra of half of the sensor, execute the spectra_pnCCD_NBC.ipynb notebook.
 
 %% Cell type:code id: tags:
 
 ``` python
 if corr_bools.get('pattern_class'):
     singles = []
     doubles = []
     triples = []
     quads =  []
 
     with h5py.File("{}/CORR-R{:04d}-PNCCD01-S{:05d}.h5".format(out_folder, run, sequences[seq_num]), 'r',
                    driver='core') as infile:
 
             data = infile[h5path+"/pixels_classified"][()].astype(np.float32) # classifications
             patterns = infile[h5path+"/patterns"][()].astype(np.float32) # event patterns
 
             # events' patterns indices are as follows: 100 (singles), 101 (first singles), 200 - 203 (doubles),
             # 300 - 303 (triples), and 400 - 403 (quadruples). Note that for the last three types of patterns,
             # there are left, right, up, and down indices.
 
             # Separating the events:
             # Singles and First Singles:
             for s in range(100, 102):
                 single = copy.copy(data[...])
                 single[patterns != s] = np.nan
                 singles.append(single)
 
 
             for d in range(200, 204):
                 double = copy.copy(data[...])
                 double[patterns != d] = np.nan
                 doubles.append(double)
 
             for t in range(300, 304):
                 triple = copy.copy(data[...])
                 triple[patterns != t] = np.nan
                 triples.append(triple)
 
             for q in range(400, 404):
                 quad = copy.copy(data[...])
                 quad[patterns != q] = np.nan
                 quads.append(quad)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 if corr_bools.get('pattern_class'):
     hA = 0
     h = 0
     for single in singles:
         hs, e = np.histogram(single.flatten(), bins=event_bins, range=b_range) # h: histogram counts, e: bin edges
         h += hs
         hA += hs # hA: counts all events (see below)
 
 
     # bin edges array has one extra element => need to plot from 0 to the one before the last element to have the
     # same size as h-array => in what follows, we use e[:-1] (-1 means one before the last element)
 
     #TODO adapt the title depending on corrections
     display(Markdown('### Histograms of Offset, Common Mode and Gain Corrected Events for One Sequence Only'))
     fig = plt.figure(figsize=(10, 7))
     ax = fig.add_subplot(111)
     ax.step(e[:-1], h, color='blue', label='Events Involving Single Pixels Only')
     ax.semilogy() # y-axis is log, x-axis is linear
     ax.set_xlabel("Energy (ADU) [{} bins per {} ADU.]".format(event_bins, b_range[1]-b_range[0]))
     ax.set_ylabel("Gain Corrected Events (counts) for One Sequence")
     ax.set_xlim(x_range)
 
     h = 0
     for double in doubles:
         hd, e = np.histogram(double.flatten(), bins=event_bins, range=b_range)
         h += hd
         hA += hd
 
     ax.step(e[:-1], h, color='red', label='Events Splitting on Double Pixels')
 
     h = 0
     for triple in triples:
         ht, e = np.histogram(triple.flatten(), bins=event_bins, range=b_range)
         h += ht
         hA += ht
 
     ax.step(e[:-1], h, color='green', label='Events Splitting on Triple Pixels')
 
     h = 0
     for quad in quads:
         hq, e = np.histogram(quad.flatten(), bins=event_bins, range=b_range)
         h += hq
         hA += hq
 
     ax.step(e[:-1], h, color='purple', label='Events Splitting on Quadruple Pixels')
 
     ax.step(e[:-1], hA, color='grey', label='All Valid Events')
     l = ax.legend()
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 ```