notebooks/AGIPD/Characterize_AGIPD_Gain_FlatFields_NBC.ipynb

 %% Cell type:markdown id: tags:
 
 # Gain Characterization #
 
 
 %% Cell type:code id: tags:
 
 ``` python
 in_folder = "/gpfs/exfel/exp/SPB/202030/p900138/scratch/karnem/r0203_r0204_v01/" # the folder to read histograms from, required
 out_folder = ""  # the folder to output to, required
 hist_file_template = "hists_m{:02d}_sum.h5" # the template to use to access histograms
 modules = [10] # modules to correct, set to -1 for all, range allowed
 
 raw_folder = "/gpfs/exfel/exp/MID/202030/p900137/raw" # Path to raw image data used to create histograms
 proc_folder = "" # Path to corrected image data used to create histograms
 
 run = 449 # of the run of image data used to create histograms
 
 karabo_id = "MID_DET_AGIPD1M-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/' # 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
+ctrl_source_template = '{}/MDL/FPGA_COMP' # path to control information
 karabo_id_control = "MID_IRU_AGIPD1M1" # karabo-id for control device
 karabo_da_control = 'AGIPD1MCTRL00' # karabo DA for control infromation
 
 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
 local_output = True # output constants locally
 db_output = False # output constants to database
 
 # Fit parameters
 peak_range = [-30, 30, 35, 70, 95, 135, 145, 220] # where to look for the peaks, [a0, b0, a1, b1, ...] exactly 8 elements
 peak_width_range = [0, 30, 0, 35, 0, 40, 0, 45] # fit limits on the peak widths, [a0, b0, a1, b1, ...] exactly 8 elements
 peak_norm_range = [0.0, -1, 0, -1, 0, -1, 0, -1] #
 
 # Bad-pixel thresholds (gain evaluation error). Contribute to BadPixel bit "Gain_Evaluation_Error"
 peak_lim = [-30, 30] # Limit of position of noise peak
 d0_lim = [10, 80] # hard limits for distance between noise and first peak
 peak_width_lim = [0.9, 1.55, 0.95, 1.65] # hard limits on the peak widths for first and second peak, in units of the noise peak. 4 parameters.
 chi2_lim = [0, 3.0] # Hard limit on chi2/nDOF value
 
 intensity_lim = 15 # Threshold on standard deviation of a histogram in ADU. Contribute to BadPixel bit "No_Entry"
 gain_lim = [0.8, 1.2] # Threshold on gain in relative number. Contribute to BadPixel bit "Gain_deviation"
 
 cell_range = [1, 3] # range of cell to be considered, [0,0] for all
 pixel_range = [0, 0, 32, 32] # range of pixels x1,y1,x2,y2 to consider [0,0,512,128] for all
 max_bins = 0 # Maximum number of bins to consider, 0 for all bins
 batch_size = [1, 8, 8] # batch size: [cell,x,y]
 fit_range = [0, 0] # range of a histogram considered for fitting in ADU. Dynamically evaluated in case [0,0]
 n_peaks_fit = 4 # Number of gaussian peaks to fit including noise peak
 fix_peaks = False # Fix distance between photon peaks
 do_minos = False # This is additional feature of minuit to evaluate errors.
 sigma_limit = 0. # If >0, repeat fit keeping only bins within mu +- sigma_limit*sigma
 
 # Detector conditions
 # NOTE: The below parameters are needed for the summary notebook when running through xfel-calibrate.
 mem_cells = -1  # number of memory cells used, negative values for auto-detection.
 bias_voltage = 300  # Bias voltage.
 acq_rate = 0.  # the detector acquisition rate, use 0 to try to auto-determine.
 gain_setting = -1  # the gain setting, negative values for auto-detection.
 photon_energy = 8.05  # photon energy in keV.
 integration_time = -1  # integration time, negative values for auto-detection.
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 import glob
 import os
 import traceback
 import warnings
 from multiprocessing import Pool
 
 import h5py
 import matplotlib.pyplot as plt
 import numpy as np
 import sharedmem
 import XFELDetAna.xfelpyanatools as xana
 from cal_tools.agipdutils_ff import (
     BadPixelsFF,
     any_in,
     fit_n_peaks,
     gaussian,
     gaussian_sum,
     get_mask,
     get_starting_parameters,
     set_par_limits,
 )
 from cal_tools.ana_tools import get_range, save_dict_to_hdf5
 from iminuit import Minuit
 from XFELDetAna.plotting.heatmap import heatmapPlot
 from XFELDetAna.plotting.simpleplot import simplePlot
 
 # %load_ext autotime
 %matplotlib inline
 warnings.filterwarnings('ignore')
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 peak_range = np.reshape(peak_range,(4,2))
 peak_width_range = np.reshape(peak_width_range,(4,2))
 peak_width_lim = np.reshape(peak_width_lim,(2,2))
 peak_norm_range = [None if x == -1 else x for x in peak_norm_range]
 peak_norm_range = np.reshape(peak_norm_range,(4,2))
 module = modules[0]
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 def idx_gen(batch_start, batch_size):
     """
     This generator iterate across pixels and memory cells starting
     from batch_start until batch_start+batch_size
     """
     for c_idx in range(batch_start[0], batch_start[0]+batch_size[0]):
         for x_idx in range(batch_start[1], batch_start[1]+batch_size[1]):
             for y_idx in range(batch_start[2], batch_start[2]+batch_size[2]):
                 yield(c_idx, x_idx, y_idx)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 n_pixels_x = pixel_range[2]-pixel_range[0]
 n_pixels_y = pixel_range[3]-pixel_range[1]
 
 hist_data = {}
 with h5py.File(f"{in_folder}/{hist_file_template.format(module)}", 'r') as hf:
     hist_data['cellId'] = np.array(hf['cellId'][()])
     hist_data['hRange'] = np.array(hf['hRange'][()])
     hist_data['nBins'] = np.array(hf['nBins'][()])
 
     if cell_range == [0,0]:
         cell_range[1] = hist_data['cellId'].shape[0]
 
     if max_bins == 0:
         max_bins = hist_data['nBins']
 
     hist_data['cellId'] = hist_data['cellId'][cell_range[0]:cell_range[1]]
     hist_data['hist'] = np.array(hf['hist'][cell_range[0]:cell_range[1], :max_bins, :])
 
 n_cells = cell_range[1]-cell_range[0]
 hist_data['hist'] = hist_data['hist'].reshape(n_cells, max_bins, 512, 128)
 hist_data['hist'] = hist_data['hist'][:,:,pixel_range[0]:pixel_range[2],pixel_range[1]:pixel_range[3]]
 
 print(f'Data shape {hist_data["hist"].shape}')
 
 bin_edges = np.linspace(hist_data['hRange'][0], hist_data['hRange'][1], int(hist_data['nBins']+1))
 x = (bin_edges[1:] + bin_edges[:-1])[:max_bins] * 0.5
 
 
 batches = []
 for c_idx in range(0, n_cells, batch_size[0]):
     for x_idx in range(0, n_pixels_x, batch_size[1]):
         for y_idx in range(0, n_pixels_y, batch_size[2]):
             batches.append([c_idx,x_idx,y_idx])
 
 print(f'Number of batches {len(batches)}')
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 def fit_batch(batch_start):
     current_result = {}
     prev = None
     for c_idx, x_idx, y_idx in idx_gen(batch_start, batch_size):
         try:
             y = hist_data['hist'][c_idx, :, x_idx, y_idx]
 
             if prev is None:
                 prev, _ = get_starting_parameters(x, y, peak_range, n_peaks=n_peaks_fit)
 
             if fit_range == [0, 0]:
                 frange = (prev[f'g0mean']-2*prev[f'g0sigma'],
                           prev[f'g{n_peaks_fit-1}mean'] + prev[f'g{n_peaks_fit-1}sigma'])
             else:
                 frange = fit_range
 
             set_par_limits(prev, peak_range, peak_norm_range,
                            peak_width_range, n_peaks_fit)
             minuit = fit_n_peaks(x, y, prev, frange,
                                  do_minos=do_minos, n_peaks=n_peaks_fit,
                                  fix_d01=fix_peaks, sigma_limit=sigma_limit,)
 
             ndof = np.rint(frange[1]-frange[0])-len(minuit.args) ## FIXME: this line is wrong if fix_peaks is True
             current_result['chi2_ndof'] = minuit.fval/ndof
             res = minuit.fitarg
             if fix_peaks : ## set g2 and g3 mean correctly
                 for i in range(2,n_peaks_fit):
                     d = res[f'g1mean'] - res[f'g0mean']
                     res[f'g{i}mean'] = res[f'g0mean'] + d*i
             current_result.update(res)
             current_result.update(minuit.get_fmin())
 
             fit_result['chi2_ndof'][c_idx, x_idx, y_idx] = current_result['chi2_ndof']
 
             for key in res.keys():
                 if key in fit_result:
                     fit_result[key][c_idx, x_idx, y_idx] = res[key]
 
             fit_result['mask'][c_idx, x_idx, y_idx] = get_mask(current_result,
                                                                     peak_lim,
                                                                     d0_lim, chi2_lim,
                                                                     peak_width_lim)
         except Exception as e:
             fit_result['mask'][c_idx, x_idx,
                                     y_idx] = BadPixelsFF.FIT_FAILED.value
             print(c_idx, x_idx, y_idx, e, traceback.format_exc())
 
         if fit_result['mask'][c_idx, x_idx, y_idx] == 0:
             prev = res
         else:
             prev = None
 ```
 
 %% Cell type:markdown id: tags:
 
 # Single fit ##
 
 Left plot shows starting parameters for fitting. Right plot shows result of the fit. Errors are evaluated with minos.
 
 %% Cell type:code id: tags:
 
 ``` python
 hist = hist_data['hist'][1,:,1, 1]
 prev, shapes = get_starting_parameters(x, hist, peak_range, n_peaks=n_peaks_fit)
 
 if fit_range == [0, 0]:
     frange = (prev[f'g0mean']-2*prev[f'g0sigma'],
               prev[f'g3mean'] + prev[f'g3sigma'])
 else:
     frange = fit_range
 
 set_par_limits(prev, peak_range, peak_norm_range,
                peak_width_range, n_peaks=n_peaks_fit)
 minuit = fit_n_peaks(x, hist, prev, frange,
                      do_minos=True, n_peaks=n_peaks_fit,
                      fix_d01=fix_peaks,
                      sigma_limit=sigma_limit,
                     )
 print (minuit.get_fmin())
 minuit.print_matrix()
 print(minuit.get_param_states())
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 res = minuit.fitarg
 if fix_peaks :
     for i in range(2,n_peaks_fit):
         d = res[f'g1mean'] - res[f'g0mean']
         res[f'g{i}mean'] = res[f'g0mean'] + d*i
 err = minuit.errors
 p = minuit.args
 ya = np.arange(0,1e4)
 y = gaussian_sum(x,n_peaks_fit, *p)
 peak_colors = ['g', 'y', 'b', 'orange']
 
 peak_hist = hist.copy()
 d=[]
 if sigma_limit > 0 :
     sel2 = (np.abs(x - res['g0mean']) < sigma_limit*res['g0sigma']) | \
            (np.abs(x - res['g1mean']) < sigma_limit*res['g1sigma']) | \
            (np.abs(x - res['g2mean']) < sigma_limit*res['g2sigma']) | \
            (np.abs(x - res['g3mean']) < sigma_limit*res['g3sigma'])
     peak_hist[~sel2] = 0
     valley_hist = hist.copy()
     valley_hist[sel2] = 0
     d.append({'x': x,
               'y': valley_hist.astype(np.float64),
               'y_err': np.sqrt(valley_hist),
               'drawstyle': 'bars',
               'errorstyle': 'bars',
               'transparency': '95%',
               'errorcoarsing': 3,
               'label': f'X-ray Data)'
              })
     htitle = f'X-ray Data, (μ±{sigma_limit:0.1f}σ)'
 else :
     htitle = 'X-ray Data'
 
 d.append({'x': x,
           'y': peak_hist.astype(np.float64),
           'y_err': np.sqrt(peak_hist),
           'drawstyle': 'bars',
           'errorstyle': 'bars',
           'errorcoarsing': 3,
           'label': htitle,
          }
         )
 d.append({'x': x,
           'y': y,
           'y2': (hist-y)/np.sqrt(hist),
           'drawstyle':'line',
           'drawstyle2': 'steps-mid',
           'label': 'Fit'
          }
         )
 
 for i in range(n_peaks_fit):
     d.append({'x': x,
              'y': gaussian(x, res[f'g{i}n'], res[f'g{i}mean'], res[f'g{i}sigma']),
              'drawstyle':'line',
              'color': peak_colors[i],
              })
     d.append({'x': np.full_like(ya, res[f'g{i}mean']),
               'y': ya,
               'drawstyle': 'line',
               'linestyle': 'dashed',
               'color': peak_colors[i],
               'label': f'peak {i} = {res[f"g{i}mean"]:0.1f} $ \pm $ {err[f"g{i}mean"]:0.2f} ADU' })
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 fig, (ax1, ax2) = plt.subplots(1, 2)
 fig.set_size_inches(16, 7)
 for i, shape in enumerate(shapes):
     idx = shape[3]
     ax1.errorbar(
         x[idx], hist[idx],
         np.sqrt(hist[idx]),
         marker='+', ls='',
     )
     yg = gaussian(x[idx], *shape[:3])
     l = f'Peak {i}: {shape[1]:0.1f} $ \pm $ {shape[2]:0.2f} ADU'
     ax1.plot(x[idx], yg, label=l)
 ax1.grid(True)
 ax1.set_xlabel("Signal [ADU]")
 ax1.set_ylabel("Counts")
 ax1.legend(ncol=2)
 
 _ = xana.simplePlot(
     d,
     use_axis=ax2,
     x_label='Signal [ADU]',
     y_label='Counts',
     secondpanel=True, y_log=False,
     x_range=(frange[0], frange[1]),
     y_range=(1., np.max(hist)*1.6),
     legend='top-left-frame-ncol2',
 )
 
 plt.show()
 ```
 
 %% Cell type:markdown id: tags:
 
 ## All fits ##
 
 %% Cell type:code id: tags:
 
 ``` python
 # Allocate memory for fit results
 fit_result = {}
 keys = list(minuit.fitarg.keys())
 keys = [x for x in keys if 'limit_' not in x and 'fix_' not in x]
 keys += ['chi2_ndof', 'mask', 'gain']
 for key in keys:
     dtype = 'f4'
     if key == 'mask':
         dtype = 'i4'
     fit_result[key] = sharedmem.empty([n_cells, n_pixels_x, n_pixels_y], dtype=dtype)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # Perform fitting
 with Pool() as pool:
     const_out = pool.map(fit_batch, batches)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # Evaluate bad pixels
 fit_result['gain'] = (fit_result['g1mean'] - fit_result['g0mean'])/photon_energy
 
 # Calculate histogram width and evaluate cut
 h_sums = np.sum(hist_data['hist'], axis=1)
 hist_norm = hist_data['hist'] / h_sums[:, None, :, :]
 hist_mean = np.sum(hist_norm[:, :max_bins, ...] *
                    x[None, :, None, None], axis=1)
 hist_sqr = (x[None, :, None, None] - hist_mean[:, None, ...])**2
 hist_std = np.sqrt(np.sum(hist_norm[:, :max_bins, ...] * hist_sqr, axis=1))
 
 fit_result['mask'][hist_std<intensity_lim] |= BadPixelsFF.NO_ENTRY.value
 
 # Bad pixel on gain deviation
 gains = np.copy(fit_result['gain'])
 gains[fit_result['mask']>0] = np.nan
 gain_mean = np.nanmean(gains, axis=(1,2))
 
 fit_result['mask'][fit_result['gain'] > gain_mean[:,None,None]*gain_lim[1] ] |=  BadPixelsFF.GAIN_DEVIATION.value
 fit_result['mask'][fit_result['gain'] < gain_mean[:,None,None]*gain_lim[0] ] |=  BadPixelsFF.GAIN_DEVIATION.value
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # Save fit results
 os.makedirs(out_folder, exist_ok=True)
 out_name = f'{out_folder}/fits_m{module:02d}.h5'
 print(f'Save to file: {out_name}')
 save_dict_to_hdf5({'data': fit_result}, out_name)
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Summary across cells ##
 
 %% Cell type:code id: tags:
 
 ``` python
 labels = [
     "Noise peak [ADU]",
     "First photon peak [ADU]",
     f"gain [ADU/keV] $\gamma$={photon_energy} [keV]",
     "$\chi^2$/nDOF",
     "Fraction of bad pixels",
 ]
 
 for i, key in enumerate(['g0mean', 'g1mean', 'gain', 'chi2_ndof', 'mask']):
     fig = plt.figure(figsize=(20,5))
     ax = fig.add_subplot(121)
     data = fit_result[key]
     if key == 'mask':
         data = data > 0
         vmin, vmax = [0, 1]
     else:
         vmin, vmax = get_range(data, 5)
     _ = heatmapPlot(
         np.mean(data, axis=0).T,
         add_panels=False, cmap='viridis', use_axis=ax,
         vmin=vmin, vmax=vmax, lut_label=labels[i]
     )
 
     if key != 'mask':
         vmin, vmax = get_range(data, 7)
         ax = fig.add_subplot(122)
         _ = xana.histPlot(
             ax, data.flatten(),
             bins=45,range=[vmin, vmax],
             log=True,color='red',histtype='stepfilled'
         )
         ax.set_xlabel(labels[i])
         ax.set_ylabel("Counts")
 ```
 
 %% Cell type:markdown id: tags:
 
 ## histograms of fit parameters ##
 
 %% Cell type:code id: tags:
 
 ``` python
 fig = plt.figure(figsize=(10, 5))
 ax = fig.add_subplot(111)
 a = ax.hist(hist_std.flatten(), bins=100, range=(0,100) )
 ax.plot([intensity_lim, intensity_lim], [0, np.nanmax(a[0])], linewidth=1.5, color='red' )
 ax.set_xlabel('Histogram width [ADU]', fontsize=14)
 ax.set_ylabel('Number of histograms', fontsize=14)
 ax.set_title(f'{hist_std[hist_std<intensity_lim].shape[0]} histograms below threshold in {intensity_lim} ADU',
               fontsize=14, fontweight='bold')
 ax.grid()
 ax.set_yscale('log')
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 def plot_par_distr(par):
     fig = plt.figure(figsize=(16, 5))
     sel = fit_result['mask'] == 0
 
     for i in range(n_peaks_fit) :
         data=fit_result[f"g{i}{par}"]
         plt_range=(-1,50)
         if par =='mean':
             plt_range=[peak_range[i][0] ,peak_range[i][1]]
 
         num_bins = int(plt_range[1] - plt_range[0])
         ax = fig.add_subplot(1,n_peaks_fit,i+1)
         _ = xana.histPlot(ax,data.flatten(),
                           bins= num_bins,range=plt_range,
                           log=True,color='red',
                           label='all fits',)
 
         a = ax.hist(data[sel].flatten(),
                     bins=num_bins, range=plt_range,
                     log=True,color='g',
                     label='good fits only',
                    )
         ax.set_xlabel(f"g{i} {par} [ADU]")
         ax.legend()
 
 plot_par_distr('mean')
 plot_par_distr('sigma')
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 sel = fit_result['mask'] == 0
 
 dsets = {'d01 [ADU]':fit_result[f"g1mean"]-fit_result[f"g0mean"],
          'gain [ADU/keV]':fit_result[f"gain"],
          'gain relative to module mean':fit_result[f"gain"]/np.nanmean(gain_mean),
         }
 fig = plt.figure(figsize=(16,5))
 for i, (par, data) in enumerate(dsets.items()):
     ax = fig.add_subplot(1, 3, i+1)
     plt_range=get_range(data, 10)
     num_bins = 100
     _ = xana.histPlot(ax,data.flatten(),
                       bins= num_bins,range=plt_range,
                       log=True,color='red',
                       label='all fits',)
 
     a = ax.hist(data[sel].flatten(),
                 bins=num_bins, range=plt_range,
                 log=True,color='g',
                 label='good fits only',
                )
     ax.set_xlabel(f"{par}")
     ax.legend()
     if 'd01' in par :
         ax.axvline(d0_lim[0])
         ax.axvline(d0_lim[1])
     if 'rel' in par :
         ax.axvline(gain_lim[0])
         ax.axvline(gain_lim[1])
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Summary across pixels ##
 
 Mean and median values are calculated across all pixels for each memory cell.
 
 %% Cell type:code id: tags:
 
 ``` python
 def plot_error_band(key, x, ax):
 
     cdata = np.copy(fit_result[key])
     cdata[fit_result['mask']>0] = np.nan
 
     mean = np.nanmean(cdata, axis=(1,2))
     median = np.nanmedian(cdata, axis=(1,2))
     std = np.nanstd(cdata, axis=(1,2))
     mad = np.nanmedian(np.abs(cdata - median[:,None,None]), axis=(1,2))
 
     ax.plot(x, mean, 'k', color='#3F7F4C', label=" mean value ")
     ax.plot(x, median, 'o', color='red', label=" median value ")
     ax.fill_between(x, mean-std, mean+std,
                      alpha=0.6, edgecolor='#3F7F4C', facecolor='#7EFF99',
                      linewidth=1, linestyle='dashdot', antialiased=True,
                      label=" mean value $ \pm $ std ")
 
     ax.fill_between(x, median-mad, median+mad,
                      alpha=0.3, edgecolor='red', facecolor='red',
                      linewidth=1, linestyle='dashdot', antialiased=True,
                      label=" median value $ \pm $ mad ")
 
     if f'error_{key}' in fit_result:
         cerr = np.copy(fit_result[f'error_{key}'])
         cerr[fit_result['mask']>0] = np.nan
 
         meanerr = np.nanmean(cerr, axis=(1,2))
         ax.fill_between(x, mean-meanerr, mean+meanerr,
                  alpha=0.6, edgecolor='#089FFF', facecolor='#089FFF',
                  linewidth=1, linestyle='dashdot', antialiased=True,
                  label=" mean fit error ")
 
 
 x = np.linspace(*cell_range, n_cells)
 
 for i, key in enumerate(['g0mean', 'g1mean', 'gain', 'chi2_ndof']):
 
     fig = plt.figure(figsize=(10, 5))
     ax = fig.add_subplot(111)
     plot_error_band(key, x, ax)
 
     ax.set_xlabel('Memory Cell ID', fontsize=14)
     ax.set_ylabel(labels[i], fontsize=14)
     ax.grid()
     ax.legend()
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Cut flow ##
 
 %% Cell type:code id: tags:
 
 ``` python
 fig, ax = plt.subplots()
 fig.set_size_inches(10, 5)
 
 n_bars = 8
 x = np.arange(n_bars)
 width = 0.3
 
 msk = fit_result['mask']
 n_fits = np.prod(msk.shape)
 y = [any_in(msk, BadPixelsFF.FIT_FAILED.value),
      any_in(msk, BadPixelsFF.FIT_FAILED.value | BadPixelsFF.ACCURATE_COVAR.value),
      any_in(msk, BadPixelsFF.FIT_FAILED.value | BadPixelsFF.ACCURATE_COVAR.value |
            BadPixelsFF.CHI2_THRESHOLD.value),
      any_in(msk, BadPixelsFF.FIT_FAILED.value | BadPixelsFF.ACCURATE_COVAR.value |
            BadPixelsFF.CHI2_THRESHOLD.value | BadPixelsFF.GAIN_THRESHOLD.value),
      any_in(msk, BadPixelsFF.FIT_FAILED.value | BadPixelsFF.ACCURATE_COVAR.value |
            BadPixelsFF.CHI2_THRESHOLD.value | BadPixelsFF.GAIN_THRESHOLD.value |
            BadPixelsFF.NOISE_PEAK_THRESHOLD.value),
      any_in(msk, BadPixelsFF.FIT_FAILED.value | BadPixelsFF.ACCURATE_COVAR.value |
            BadPixelsFF.CHI2_THRESHOLD.value | BadPixelsFF.GAIN_THRESHOLD.value |
            BadPixelsFF.NOISE_PEAK_THRESHOLD.value | BadPixelsFF.PEAK_WIDTH_THRESHOLD.value),
      any_in(msk, BadPixelsFF.FIT_FAILED.value | BadPixelsFF.ACCURATE_COVAR.value |
            BadPixelsFF.CHI2_THRESHOLD.value | BadPixelsFF.GAIN_THRESHOLD.value |
            BadPixelsFF.NOISE_PEAK_THRESHOLD.value | BadPixelsFF.PEAK_WIDTH_THRESHOLD.value
            | BadPixelsFF.NO_ENTRY.value),
      any_in(msk, BadPixelsFF.FIT_FAILED.value | BadPixelsFF.ACCURATE_COVAR.value |
            BadPixelsFF.CHI2_THRESHOLD.value | BadPixelsFF.GAIN_THRESHOLD.value |
            BadPixelsFF.NOISE_PEAK_THRESHOLD.value | BadPixelsFF.PEAK_WIDTH_THRESHOLD.value
            | BadPixelsFF.NO_ENTRY.value| BadPixelsFF.GAIN_DEVIATION.value)
     ]
 
 y2 = [any_in(msk, BadPixelsFF.FIT_FAILED.value),
      any_in(msk, BadPixelsFF.ACCURATE_COVAR.value),
      any_in(msk, BadPixelsFF.CHI2_THRESHOLD.value),
      any_in(msk, BadPixelsFF.GAIN_THRESHOLD.value),
      any_in(msk, BadPixelsFF.NOISE_PEAK_THRESHOLD.value),
      any_in(msk, BadPixelsFF.PEAK_WIDTH_THRESHOLD.value),
      any_in(msk, BadPixelsFF.NO_ENTRY.value),
      any_in(msk, BadPixelsFF.GAIN_DEVIATION.value)
     ]
 
 y = (1 - np.sum(y, axis=(1,2,3))/n_fits)*100
 y2 = (1 - np.sum(y2, axis=(1,2,3))/n_fits)*100
 
 labels = ['Fit failes',
          'Accurate covar',
          'Chi2/nDOF',
          'Gain',
          'Noise peak',
          'Peak width',
          'No Entry',
          'Gain deviation']
 
 ax.bar(x, y2, width, label='Only this cut')
 ax.bar(x, y, width, label='Cut flow')
 ax.set_xticks(x)
 ax.set_xticklabels(labels, rotation=90)
 ax.set_ylim(y[5]-0.5, 100)
 ax.grid(True)
 ax.legend()
 plt.show()
 ```

notebooks/AGIPD/Characterize_AGIPD_Gain_FlatFields_Summary.ipynb

 %% Cell type:markdown id: tags:
 
 # Gain Characterization Summary #
 
 
 %% Cell type:code id: tags:
 
 ``` python
 in_folder = "" # in this notebook, in_folder is not used as the data source is in the destination folder
 out_folder = ""  # the folder to output to, required
 metadata_folder = ""  # Directory containing calibration_metadata.yml when run by xfel-calibrate
 hist_file_template = "hists_m{:02d}_sum.h5"
 proc_folder = "" # Path to corrected image data used to create histograms and validation plots
 raw_folder = "/gpfs/exfel/exp/MID/202030/p900137/raw"  # folder of raw data. This is used to save information of source data of generated constants, required
 run = 449 # runs of image data used to create histograms
 
 karabo_id = "MID_DET_AGIPD1M-1" # karabo karabo_id
 ctrl_source_template = '{}/MDL/FPGA_COMP' # path to control information
 karabo_id_control = "MID_EXP_AGIPD1M1" # karabo-id for control device
 
 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
 local_output = True # output constants locally
 db_output = False # output constants to database
 
 # Fit parameters
 peak_range = [-30,30,35,65,80,130,145,200] # where to look for the peaks, [a0, b0, a1, b1, ...] exactly 8 elements
-peak_width_range = [0, 30, 0, 35, 0, 40, 0, 45] # fit limits on the peak widths, [a0, b0, a1, b1, ...] exactly 8 elements
 
 # Bad-pixel thresholds
 d0_lim = [10, 70] # hard limits for d0 value (distance between noise and first peak)
-peak_width_lim = [0.97, 1.43, 1.03, 1.57] # hard limits on the peak widths, [a0, b0, a1, b1, ...] in units of the noise peak. 4 parameters.
-chi2_lim = [0,3.0] # Hard limit on chi2/nDOF value
 gain_lim = [0.80, 1.2] # Threshold on gain in relative number. Contribute to BadPixel bit "Gain_deviation"
 
 cell_range = [1,5] # range of cell to be considered, [0,0] for all
 pixel_range = [0,0,512,128] # range of pixels x1,y1,x2,y2 to consider [0,0,512,128] for all
-max_bins = 250 # Maximum number of bins to consider
-batch_size = [1,8,8] # batch size: [cell,x,y]
 n_peaks_fit = 4 # Number of gaussian peaks to fit including noise peak
-fix_peaks = True # Fix distance between photon peaks
-
 
 # Detector conditions
 mem_cells = -1  # number of memory cells used, negative values for auto-detection.
 bias_voltage = 0.  # Bias voltage
 acq_rate = 0.  # the detector acquisition rate, use 0 to try to auto-determine
 gain_setting = -1  # the gain setting, negative values for auto-detection.
 photon_energy = 8.05  # photon energy in keV
 integration_time = -1  # integration time, negative values for auto-detection.
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
-import glob
 import os
 import re
 import traceback
 import warnings
 from multiprocessing import Pool
 
 import h5py
 import matplotlib.pyplot as plt
 import numpy as np
 import tabulate
 from cal_tools.agipdlib import AgipdCtrl
 from cal_tools.agipdutils_ff import (
     BadPixelsFF,
     any_in,
     fit_n_peaks,
     gaussian_sum,
     get_starting_parameters,
 )
 from cal_tools.ana_tools import get_range, save_dict_to_hdf5
 from cal_tools.enums import BadPixels
 from cal_tools.tools import (
     get_dir_creation_date,
     get_pdu_from_db,
     get_report,
     module_index_to_qm,
     send_to_db
 )
-from dateutil import parser
-from extra_data import H5File, RunDirectory, stack_detector_data
+from extra_data import RunDirectory, stack_detector_data
 from extra_geom import AGIPD_1MGeometry, AGIPD_500K2GGeometry
-from iCalibrationDB import Conditions, Constants, Detectors
+from iCalibrationDB import Conditions, Constants
 from iminuit import Minuit
-from IPython.display import HTML, Latex, Markdown, display
-from XFELDetAna.plotting.heatmap import heatmapPlot
+from IPython.display import Latex, display
 from XFELDetAna.plotting.simpleplot import simplePlot
 
 %matplotlib inline
 warnings.filterwarnings('ignore')
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 peak_range = np.reshape(peak_range,(4,2))
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # Get operation conditions
 ctrl_source = ctrl_source_template.format(karabo_id_control)
 
 run_folder = f'{raw_folder}/r{run:04d}/'
 
 raw_dc = RunDirectory(run_folder)
 
 # Read operating conditions from AGIPD00 files
 instrument_src_mod = [
     s for s in list(raw_dc.all_sources) if "0CH" in s][0]
 ctrl_src = [
     s for s in list(raw_dc.all_sources) if ctrl_source in s][0]
 # Evaluate creation time
 creation_time = None
 if use_dir_creation_date:
     creation_time = get_dir_creation_date(raw_folder, run)
 
 agipd_cond = AgipdCtrl(
     run_dc=raw_dc,
     image_src=instrument_src_mod,
     ctrl_src=ctrl_src,
     raise_error=False,  # to be able to process very old data without mosetting value
 )
 if mem_cells < 0:
     mem_cells = agipd_cond.get_num_cells()
 if mem_cells is None:
     raise ValueError(f"No raw images found in {run_folder}")
 if acq_rate == 0.:
     acq_rate = agipd_cond.get_acq_rate()
 if gain_setting < 0:
     gain_setting = agipd_cond.get_gain_setting(creation_time)
 if bias_voltage == 0.:
     bias_voltage = agipd_cond.get_bias_voltage(karabo_id_control)
 if integration_time < 0:
     integration_time = agipd_cond.get_integration_time()
 
 # Evaluate detector instance for mapping
 instrument = karabo_id.split("_")[0]
 if instrument == "HED":
     nmods = 8
 else:
     nmods = 16
 
 print(f"Using {creation_time} as creation time")
 print(f"Operating conditions are:\n• Bias voltage: {bias_voltage}\n• Memory cells: {mem_cells}\n"
       f"• Acquisition rate: {acq_rate}\n• Gain setting: {gain_setting}\n• Integration time: {integration_time}\n"
       f"• Photon Energy: {photon_energy}\n")
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # Load constants for all modules
 keys = ['g0mean', 'g1mean', 'gain', 'chi2_ndof', 'mask']
 all_keys = set(keys)
 for i in range(n_peaks_fit) :
     all_keys.add(f'g{i}mean')
     all_keys.add(f'g{i}sigma')
 
 fit_data = {}
 labels = {'g0mean': 'Noise peak position [ADU]',
           'g1mean': 'First photon peak [ADU]',
           'gain': f"Gain [ADU/keV], $\gamma$={photon_energy} [keV]",
           'chi2_ndof': '$\chi^2$/nDOF',
           'mask': 'Fraction of bad pixels over cells' }
 
 modules = []
 karabo_da = []
 for mod in range(nmods):
     qm = module_index_to_qm(mod)
     fit_data[mod] = {}
     try:
         hf = h5py.File(f'{out_folder}/fits_m{mod:02d}.h5', 'r')
         shape = hf['data/g0mean'].shape
         for key in keys:
             fit_data[mod][key] = hf[f'data/{key}'][()]
 
         print(f"{in_folder}/{hist_file_template.format(mod)}")
         modules.append(mod)
         karabo_da.append(f"AGIPD{mod:02d}")
     except Exception as e:
         err = f"Error: {e}\nError traceback: {traceback.format_exc()}"
         print(f"No fit data available for module {qm}")
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # Calculate SlopesFF and BadPixels to be send to DB
 bpmask = {}
 slopesFF = {}
 
 for mod in modules:
     bpmask[mod] = np.zeros(fit_data[mod]['mask'].shape).astype(np.int32)
     bpmask[mod][ any_in(fit_data[mod]['mask'], BadPixelsFF.NO_ENTRY.value) ] = BadPixels.FF_NO_ENTRIES.value
     bpmask[mod][ any_in(fit_data[mod]['mask'],
                         BadPixelsFF.GAIN_DEVIATION.value) ] |= BadPixels.FF_GAIN_DEVIATION.value
     bpmask[mod][ any_in(fit_data[mod]['mask'],
                         BadPixelsFF.FIT_FAILED.value | BadPixelsFF.ACCURATE_COVAR.value |
                         BadPixelsFF.CHI2_THRESHOLD.value | BadPixelsFF.GAIN_THRESHOLD.value |
                         BadPixelsFF.NOISE_PEAK_THRESHOLD.value | BadPixelsFF.PEAK_WIDTH_THRESHOLD.value) ] |= BadPixels.FF_GAIN_EVAL_ERROR.value
 
     # Set value for bad pixel to average across pixels for a given module
     slopesFF[mod] = np.copy(fit_data[mod]['gain'])
     slopesFF[mod][fit_data[mod]['mask']>0] = np.nan
     gain_mean = np.nanmean(slopesFF[mod], axis=(1,2))
 
     for i in range(slopesFF[mod].shape[0]):
         slopesFF[mod][i][ fit_data[mod]['mask'][i] > 0 ] = gain_mean[i]
 
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # Read report path and create file location tuple to add with the injection
 proposal = list(filter(None, raw_folder.strip('/').split('/')))[-2]
 file_loc = f'Proposal: {proposal}, Run: {run}'
 
 report = get_report(metadata_folder)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # set the operating condition
 condition = Conditions.Illuminated.AGIPD(mem_cells, bias_voltage, 9.2,
                                          pixels_x=512, pixels_y=128, beam_energy=None,
                                          acquisition_rate=acq_rate, gain_setting=gain_setting,
                                          integration_time=integration_time)
 
 # Modify acceptable deviations for integration time condition if and only if
 # the integration time is not using the standard value (12).
 if integration_time != 12:
     for p in condition.parameters:
         if p.name == 'Integration Time':
             p.lower_deviation = 5
             p.upper_deviation = 5
 
 # Retrieve a list of all modules corresponding to processed karabo_das
 db_modules = get_pdu_from_db(karabo_id, karabo_da, Constants.AGIPD.SlopesFF(),
                              condition, cal_db_interface,
                              snapshot_at=creation_time)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # Send constants to DB
 def send_const(mod, pdu):
     try:
         # gain
         constant = Constants.AGIPD.SlopesFF()
         constant.data = np.moveaxis(np.moveaxis(slopesFF[mod], 0, 2), 0, 1)
         send_to_db(
             pdu, karabo_id, constant, condition, file_loc,
             report, cal_db_interface, creation_time,
             timeout=cal_db_timeout,
         )
 
         # bad pixels
         constant_bp = Constants.AGIPD.BadPixelsFF()
         constant_bp.data = np.moveaxis(np.moveaxis(bpmask[mod], 0, 2), 0, 1)
         send_to_db(
             pdu, karabo_id, constant_bp, condition, file_loc,
             report, cal_db_interface, creation_time,
             timeout=cal_db_timeout,
         )
 
     except Exception as e:
         err = f"Error: {e}\nError traceback: {traceback.format_exc()}"
         when = None
 
 # Check, if we have a shape we expect
 if db_output:
     if slopesFF[modules[0]].shape == (mem_cells, 512, 128):
         with Pool(processes=len(modules)) as pool:
             const_out = pool.starmap(send_const, zip(modules, db_modules))
     else:
         print(f"Constants are not sent to the DB because of the shape mismatsh")
         print(f"Expected {(mem_cells, 512, 128)}, observed {slopesFF[modules[0]].shape}")
 
 
 condition_dict ={}
 
 for entry in condition.to_dict()['parameters']:
     key = entry.pop('parameter_name')
     del entry['description']
     del entry['flg_available']
     condition_dict[key] = entry
 
 # Create the same file structure as database constants files, in which
 # each constant type has its corresponding condition and data.
 if local_output:
     for mod, pdu in zip(modules, db_modules):
         qm = module_index_to_qm(mod)
         file = f"{out_folder}/slopesff_bpmask_module_{qm}.h5"
         dic = {
             pdu:{
                'SlopesFF': {
                    0:{
                        'condition': condition_dict,
                        'data': np.moveaxis(np.moveaxis(slopesFF[mod],0,2),0,1)}
                },
                'BadPixelsFF':{
                    0:{
                        'condition': condition_dict,
                        'data': np.moveaxis(np.moveaxis(bpmask[mod],0,2),0,1)}
                },
            }
         }
         save_dict_to_hdf5(dic, file)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 #Define AGIPD geometry
 #To do: find the better way to do it?
 if instrument == "HED":
     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
 # Create the arrays that will be used for figures.
 # A dictionary contains all the data for each of the processing stages (gains, mean, slopesFF...).
 # Each array correponds to the data for all processed modules.
 # These are updated with their fit/slopes data in the following loops.
 if cell_range==[0,0]:
     cell_range[1] = shape[0]
 
 const_data = {}
 for key in keys:
     const_data[key] = np.full((nmods, shape[0],512,128), np.nan)
     for i in range(nmods):
         if key in fit_data[i]:
             const_data[key][i,:,pixel_range[0]:pixel_range[2],
                                pixel_range[1]:pixel_range[3]] = fit_data[i][key]
 
 const_data['slopesFF'] = np.full((nmods, shape[0],512,128), np.nan)
 labels['slopesFF'] = f'slopesFF [ADU/keV], $\gamma$={photon_energy} [keV]'
 for i in range(nmods):
     if i in slopesFF:
         const_data['slopesFF'][i,:,pixel_range[0]:pixel_range[2],
                                pixel_range[1]:pixel_range[3]] = slopesFF[i]
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Summary across pixels ##
 
 %% Cell type:code id: tags:
 
 ``` python
 for key in const_data.keys():
     fig = plt.figure(figsize=(20,20))
     ax = fig.add_subplot(111)
     if key=='mask':
         data = np.nanmean(const_data[key]>0, axis=1)
         vmin, vmax = (0,1)
     else:
         data = np.nanmean(const_data[key], axis=1)
         vmin, vmax = get_range(data, 5)
     ax = geom.plot_data_fast(data, ax=ax, cmap="jet", vmin=vmin, vmax=vmax, figsize=(20,20))
     _ = ax.set_title(labels[key])
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Summary histograms ##
 
 %% Cell type:code id: tags:
 
 ``` python
 sel = (const_data['mask'] == 0)
 
 module_mean = np.nanmean(const_data[f"gain"],axis=(1,2,3))
 module_mean = module_mean[:,np.newaxis,np.newaxis,np.newaxis]
 
 
 dsets = {'d01 [ADU]':const_data[f"g1mean"]-const_data[f"g0mean"],
          'gain [ADU/keV]':const_data[f"gain"],
          'gain relative to module mean':const_data[f"gain"]/module_mean,
         }
 fig = plt.figure(figsize=(16,5))
 for i, (par, data) in enumerate(dsets.items()):
     ax = fig.add_subplot(1, 3, i+1)
     plt_range= np.nanmin(data), np.nanmax(data)
     if 'd01' in par :
         ax.axvline(d0_lim[0])
         ax.axvline(d0_lim[1])
     elif 'rel' in par :
         ax.axvline(gain_lim[0])
         ax.axvline(gain_lim[1])
     num_bins = 100
     _ = ax.hist(data.flatten(),
                   bins= num_bins,range=plt_range,
                   log=True,color='red',
                   label='all fits',)
 
     a = ax.hist(data[sel].flatten(),
                 bins=num_bins, range=plt_range,
                 log=True,color='g',
                 label='good fits only',
                )
     ax.set_xlabel(f"{par}")
     ax.legend()
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Summary across cells ##
 
 Good pixels only.
 
 %% Cell type:code id: tags:
 
 ``` python
 for key in const_data.keys():
     data = np.copy(const_data[key])
     if key=='mask':
         data = data>0
     else:
         data[const_data['mask']>0] = np.nan
 
     d = []
     for i in range(nmods):
         d.append({'x': np.arange(data[i].shape[0]),
                   'y': np.nanmean(data[i], axis=(1,2)),
                   'drawstyle': 'steps-pre',
                   'label': f'{i}',
                   'linewidth': 2,
                   'linestyle': '--' if i>7 else '-'
                   })
 
     fig = plt.figure(figsize=(15, 6))
     ax = fig.add_subplot(111)
 
     _ = simplePlot(d, xrange=(-12, 510),
                         x_label='Memory Cell ID',
                         y_label=labels[key],
                         use_axis=ax,
                         legend='top-left-frame-ncol8',)
     ylim = ax.get_ylim()
     ax.set_ylim(ylim[0], ylim[1] + np.abs(ylim[1]-ylim[0])*0.2)
     ax.grid()
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Summary table ##
 
 %% Cell type:code id: tags:
 
 ``` python
 table = []
 for i in modules:
     table.append((i,
                   f"{np.nanmean(slopesFF[i]):0.1f} +- {np.nanstd(slopesFF[i]):0.2f}",
                   f"{np.nanmean(bpmask[i]>0)*100:0.1f} ({np.nansum(bpmask[i]>0)})"
                         ))
 
 all_SFF = np.array([list(sff) for sff in slopesFF.values()])
 all_MSK = np.array([list(msk) for msk in bpmask.values()])
 
 table.append(('overall',
               f"{np.nanmean(all_SFF):0.1f} +- {np.nanstd(all_SFF):0.2f}",
               f"{np.nanmean(all_MSK>0)*100:0.1f} ({np.nansum(all_MSK>0)})"
                     ))
 
 md = display(Latex(tabulate.tabulate(table, tablefmt='latex',
                                      headers=["Module", "Gain [ADU/keV]", "Bad pixels [%(Count)]"])))
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Performance plots
 
 %% Cell type:code id: tags:
 
 ``` python
 def get_trains_data(run_folder, source, include, tid=None):
     """
     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
 
     """
     run_data = RunDirectory(run_folder, include)
     if tid:
         tid, data = run_data.select('*/DET/*', source).train_from_id(tid)
         return tid, stack_detector_data(data, source, modules=nmods)
     else:
         for tid, data in run_data.select('*/DET/*', source).trains(require_all=True):
             return tid, stack_detector_data(data, source, modules=nmods)
     return None, None
 
 
 include = '*S00000*'
 tid, orig = get_trains_data(f'{proc_folder}/r{run:04d}/', 'image.data', include)
 orig = orig[cell_range[0]:cell_range[1], ...]
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # FIXME: mask bad pixels from median
 # mask = const_data['BadPixelsFF']
 
 corrections = const_data['slopesFF'] # (16,shape[0],512,128) shape[0]= cell_range[1]-cell_range[0] /
 corrections = np.moveaxis(corrections, 1, 0) # (shape[0],16,512,128)
 rel_corr = corrections/np.nanmedian(corrections)
 corrected = orig / rel_corr
 ```
 
 %% Cell type:markdown id: tags:
 
 ### Mean value not corrected (train 0)
 
 %% Cell type:code id: tags:
 
 ``` python
 fig = plt.figure(figsize=(20,20))
 ax = fig.add_subplot(111)
 odata = np.nanmean(orig, axis=0)
 vmin, vmax = get_range(odata, 5)
 ax = geom.plot_data_fast(odata, ax=ax, cmap="jet", vmin=vmin, vmax=vmax, figsize=(20,20))
 _ = ax.set_title("Original data, mean across one train")
 ```
 
 %% Cell type:markdown id: tags:
 
 ### Mean value corrected (train 0)
 
 %% Cell type:code id: tags:
 
 ``` python
 fig = plt.figure(figsize=(20,20))
 ax = fig.add_subplot(111)
 cdata = np.nanmean(corrected, axis=0)
 ax = geom.plot_data_fast(cdata, ax=ax, cmap="jet", vmin=vmin, vmax=vmax, figsize=(20,20))
 _ = ax.set_title("Corrected data, mean across one train")
 ```
 
 %% Cell type:markdown id: tags:
 
 ### Laplace transform of mean image
 
 %% Cell type:code id: tags:
 
 ``` python
 from scipy.ndimage import laplace
 
 cmax = np.max(cdata)
 omax = np.max(odata)
 clap = np.zeros_like(cdata)
 olap = np.zeros_like(odata)
 
 for i in range(nmods) :
     clap[i] = np.abs(laplace(cdata[i].astype(float)/cmax))
     olap[i] = np.abs(laplace(odata[i].astype(float)/omax))
 
 fig = plt.figure(figsize=(20,10))
 vmin, vmax = get_range(olap, 2)
 
 ax = fig.add_subplot(121)
 ax = geom.plot_data_fast(olap, ax=ax, cmap="jet", vmin=vmin, vmax=vmax, )
 _ = ax.set_title("Laplace (original data)")
 
 ax = fig.add_subplot(122)
 ax = geom.plot_data_fast(clap, ax=ax, cmap="jet", vmin=vmin, vmax=vmax, )
 _ = ax.set_title("Laplace (gain corrected data)")
 ```
 
 %% Cell type:markdown id: tags:
 
 ### Histogram of corrected and uncorrected spectrum (train 0)
 
 
 %% Cell type:code id: tags:
 
 ``` python
 ######################################
 #            FIT PEAKS
 ######################################
 
 x_range = [peak_range[0][0], peak_range[-1][-1]]
 nb = x_range[1] - x_range[0]+1
 
 sel = ~np.isnan(corrected)
 
 fig = plt.figure(figsize=(10, 5))
 ax = fig.add_subplot(111)
 y,xe, _ = ax.hist(corrected[sel].flatten(), bins=nb, range=x_range, label='corrected', alpha=0.5)
 
 # get the bin centers from the bin edges
 xc=xe[:-1]+(xe[1]-xe[0])/2
 pars, _ = get_starting_parameters(xc, y, peak_range,4)
 minuit = fit_n_peaks(xc, y, pars, x_range,fix_d01=False,sigma_limit=1)
 
 pc = minuit.args
 resc=minuit.fitarg
 
 yfc = gaussian_sum(xc,4, *pc)
 plt.plot(xc, yfc, label='corrected fit')
 
 y,_, _ = ax.hist(orig[sel].flatten(), bins=nb, range=x_range, label='original',alpha=0.5)
 pars, _ = get_starting_parameters(xc, y, peak_range,4)
 minuit = fit_n_peaks(xc, y, pars, x_range,fix_d01=False,sigma_limit=1)
 
 po = minuit.args
 reso=minuit.fitarg
 
 yfo = gaussian_sum(xc,4, *po)
 plt.plot(xc, yfo, label='original fit')
 
 plt.title(f"Signal spectrum, first train")
 plt.xlabel('[ADU]')
 plt.legend()
 plt.show()
 ```
 
 %% Cell type:markdown id: tags:
 
 ### Summary table ##
 
 %% Cell type:code id: tags:
 
 ``` python
 from scipy.stats import median_absolute_deviation as mad
 
 table = []
 headers = ["Parameter",
            "Value (original data)",
            "Value (gain corrected data)",
            "Relative difference"]
 
 for i in range(4):
     table.append((f"Sigma{i} (ADU)",
                   f"{reso[f'g{i}sigma']:0.2f} ",
                   f"{resc[f'g{i}sigma']:0.2f} ",
                   f"{(reso[f'g{i}sigma']-resc[f'g{i}sigma'])/reso[f'g{i}sigma']:0.2f} ",
                  ))
 
 ovar = np.std(odata)
 cvar = np.std(cdata)
 table.append((f"RMS of mean image",
               f"{ovar:0.3f} ",
               f"{cvar:0.3f} ",
               f"{(ovar-cvar)/ovar:0.3f} ",
              ))
 
 omin, omax = get_range(odata, 5)
 cmin, cmax = get_range(cdata, 5)
 ovar = np.std(odata[(odata > omin) & (odata<omax)])
 cvar = np.std(cdata[(cdata > cmin) & (cdata<cmax)])
 table.append((f"RMS of mean image (mu+-5sigma)",
               f"{ovar:0.3f} ",
               f"{cvar:0.3f} ",
               f"{(ovar-cvar)/ovar:0.3f} ",
              ))
 
 ovar = mad(odata.flatten())
 cvar = mad(cdata.flatten())
 table.append((f"MAD of mean image",
               f"{ovar:0.3f} ",
               f"{cvar:0.3f} ",
               f"{(ovar-cvar)/ovar:0.3f} ",
              ))
 
 ovar = np.median(olap)
 cvar = np.median(clap)
 table.append((f"Median Laplace",
               f"{ovar:0.3f} ",
               f"{cvar:0.3f} ",
               f"{(ovar-cvar)/ovar:0.3f} ",
              ))
 
 md = display(Latex(tabulate.tabulate(table,
                                      tablefmt='latex',
                                      headers=headers)))
 ```

notebooks/Jungfrau/Jungfrau_Gain_Correct_and_Verify_NBC.ipynb

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

notebooks/LPD/LPD_Correct_Fast.ipynb

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

notebooks/LPD/LPD_retrieve_constants_precorrection.ipynb

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

notebooks/ePix100/Correction_ePix100_NBC.ipynb

 %% Cell type:markdown id: tags:
 
 # ePix100 Data Correction
 
 Author: European XFEL Detector Group, Version: 2.0
 
 The following notebook provides data correction of images acquired with the ePix100 detector.
 
+The sequence of correction applied are:
+Offset --> Common Mode Noise --> Relative Gain --> Charge Sharing --> Absolute Gain.
+
+Offset, common mode and gain corrected data is saved to /data/image/pixels in the CORR files.
+
+If pattern classification is applied (charge sharing correction), this data will be saved to /data/image/pixels_classified, while the corresponding patterns will be saved to /data/image/patterns in the CORR files.
+
 %% Cell type:code id: tags:
 
 ``` python
-in_folder = "/gpfs/exfel/exp/CALLAB/202031/p900113/raw"  # input folder, required
-out_folder = "/gpfs/exfel/data/scratch/ahmedk/test/remove/epix_correct"  # output folder, required
+in_folder = "/gpfs/exfel/exp/HED/202202/p003121/raw" # input folder, required
+out_folder = ""  # output folder, required
 metadata_folder = ""  # Directory containing calibration_metadata.yml when run by xfel-calibrate
 sequences = [-1]  # sequences to correct, set to -1 for all, range allowed
 sequences_per_node = 1  # number of sequence files per cluster node if run as slurm job, set to 0 to not run SLURM parallel
-run = 9988  # which run to read data from, required
+run = 156  # which run to read data from, required
 
 # Parameters for accessing the raw data.
-karabo_id = "MID_EXP_EPIX-1"  # karabo karabo_id
+karabo_id = "HED_IA1_EPX100-1"  # karabo karabo_id
 karabo_da = "EPIX01"  # data aggregators
 db_module = ""  # module id in the database
 receiver_template = "RECEIVER"  # detector receiver template for accessing raw data files
 path_template = 'RAW-R{:04d}-{}-S{{:05d}}.h5'  # the template to use to access data
 instrument_source_template = '{}/DET/{}:daqOutput'  # instrument detector data source in h5files
 
 # Parameters affecting writing corrected data.
 chunk_size_idim = 1  # H5 chunking size of output data
 
 # Only for testing
 limit_images = 0  # ONLY FOR TESTING. process only first N images, 0 - process all.
 
 # Parameters for the calibration database.
 cal_db_interface = "tcp://max-exfl016:8015#8025"  # calibration DB interface to use
 cal_db_timeout = 300000  # timeout on caldb requests
 creation_time = ""  # The timestamp to use with Calibration DBe. Required Format: "YYYY-MM-DD hh:mm:ss" e.g. 2019-07-04 11:02:41
 
 # Conditions for retrieving calibration constants.
 bias_voltage = 200  # bias voltage
 in_vacuum = False  # detector operated in vacuum
 integration_time = -1  # Detector integration time, Default value -1 to use the value from the slow data.
 fix_temperature = -1  # fixed temperature value in Kelvin, Default value -1 to use the value from files.
-gain_photon_energy = 9.0  # Photon energy used for gain calibration
+gain_photon_energy = 8.048  # Photon energy used for gain calibration
 photon_energy = 0.  # Photon energy to calibrate in number of photons, 0 for calibration in keV
 
 # Flags to select type of applied corrections.
 pattern_classification = True  # do clustering.
 relative_gain = True  # Apply relative gain correction.
 absolute_gain = True  # Apply absolute gain correction (implies relative gain).
 common_mode = True  # Apply common mode correction.
 
 # Parameters affecting applied correction.
 cm_min_frac = 0.25  # No CM correction is performed if after masking the ratio of good pixels falls below this
 cm_noise_sigma = 5.  # CM correction noise standard deviation
 split_evt_primary_threshold = 7.  # primary threshold for split event correction
 split_evt_secondary_threshold = 5.  # secondary threshold for split event correction
 split_evt_mip_threshold = 1000.  # minimum ionizing particle threshold
 
 
 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 tabulate
 import warnings
 
 import h5py
 import pasha as psh
 import numpy as np
 import matplotlib.pyplot as plt
 from IPython.display import Latex, display
 from extra_data import RunDirectory, H5File
 from pathlib import Path
 
 from XFELDetAna import xfelpyanatools as xana
 from XFELDetAna import xfelpycaltools as xcal
 from cal_tools import h5_copy_except
 from cal_tools.epix100 import epix100lib
 from cal_tools.tools import (
     calcat_creation_time,
     get_dir_creation_date,
     get_constant_from_db,
     load_specified_constants,
     CalibrationMetadata,
 )
 from cal_tools.step_timing import StepTimer
 from iCalibrationDB import (
     Conditions,
     Constants,
 )
 
 warnings.filterwarnings('ignore')
 
 prettyPlotting = True
 
 %matplotlib inline
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 x = 708  # rows of the ePix100
 y = 768  # columns of the ePix100
 
 if absolute_gain:
     relative_gain = True
 
 plot_unit = 'ADU'
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 in_folder = Path(in_folder)
 out_folder = Path(out_folder)
 
 out_folder.mkdir(parents=True, exist_ok=True)
 
 run_folder = in_folder / f"r{run:04d}"
 
 instrument_src = instrument_source_template.format(
     karabo_id, receiver_template)
 
 print(f"Correcting run: {run_folder}")
 print(f"Instrument H5File source: {instrument_src}")
 print(f"Data corrected files are stored at: {out_folder}")
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 creation_time = calcat_creation_time(in_folder, run, creation_time)
 print(f"Using {creation_time.isoformat()} as creation time")
 
 metadata = CalibrationMetadata(metadata_folder or out_folder)
 # Constant paths are saved under retrieved-constants in calibration_metadata.yml.
 # NOTE: this notebook shouldn't overwrite calibration metadata file.
 const_yaml = metadata.get("retrieved-constants", {})
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 run_dc = RunDirectory(run_folder, _use_voview=False)
 
 seq_files = [Path(f.filename) for f in run_dc.select(f"*{karabo_id}*").files]
 
 # If a set of sequences requested to correct,
 # adapt seq_files list.
 if sequences != [-1]:
     seq_files = [f for f in seq_files if any(f.match(f"*-S{s:05d}.h5") for s in sequences)]
 
 if not len(seq_files):
     raise IndexError("No sequence files available for the selected sequences.")
 
 print(f"Processing a total of {len(seq_files)} sequence files")
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 step_timer = StepTimer()
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 step_timer.start()
 
 sensorSize = [x, y]
 # Sensor area will be analysed according to blocksize
 blockSize = [sensorSize[0]//2, sensorSize[1]//2]
 xcal.defaultBlockSize = blockSize
 memoryCells = 1  # ePIX has no memory cells
 run_parallel = False
 
 # Read control data.
 ctrl_data = epix100lib.epix100Ctrl(
     run_dc=run_dc,
     instrument_src=f"{karabo_id}/DET/{receiver_template}:daqOutput",
     ctrl_src=f"{karabo_id}/DET/CONTROL",
     )
 
 if integration_time < 0:
     integration_time = ctrl_data.get_integration_time()
     integration_time_str_add = ""
 else:
     integration_time_str_add = "(manual input)"
 
 if fix_temperature < 0:
     temperature = ctrl_data.get_temprature()
     temperature_k = temperature + 273.15
     temp_str_add = ""
 else:
     temperature_k = fix_temperature
     temperature = fix_temperature - 273.15
     temp_str_add = "(manual input)"
 
 print(f"Bias voltage is {bias_voltage} V")
 print(f"Detector integration time is set to {integration_time} \u03BCs {integration_time_str_add}")
 print(f"Mean temperature: {temperature:0.2f}°C / {temperature_k:0.2f} K {temp_str_add}")
 print(f"Operated in vacuum: {in_vacuum}")
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # Table of sequence files to process
 table = [(k, f) for k, f in enumerate(seq_files)]
 
 if len(table):
     md = display(Latex(tabulate.tabulate(
         table,
         tablefmt='latex',
         headers=["#", "file"]
     )))
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Retrieving calibration constants
 
 As a first step, dark maps have to be loaded.
 
 %% Cell type:code id: tags:
 
 ``` python
 cond_dict = {
     "bias_voltage": bias_voltage,
     "integration_time": integration_time,
     "temperature": temperature_k,
     "in_vacuum": in_vacuum,
 }
 
 dark_condition = Conditions.Dark.ePix100(**cond_dict)
 
 # update conditions with illuminated conditins.
 cond_dict.update({
         "photon_energy": gain_photon_energy
     })
 
 illum_condition = Conditions.Illuminated.ePix100(**cond_dict)
 
 const_cond = {
     "Offset": dark_condition,
     "Noise": dark_condition,
     "RelativeGain": illum_condition,
 }
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 empty_constant = np.zeros((708, 768, 1), dtype=np.float32)
 if const_yaml:  #  Used while reproducing corrected data.
     print(f"Using stored constants in {metadata.filename}")
     const_data, _ = load_specified_constants(const_yaml[karabo_da]["constants"])
     for cname, cval in const_data.items():
         if cval is None and cname != "RelativeGain":
             const_data[cname] = empty_constant
 else:  # First correction attempt.
     const_data = dict()
     for cname, condition in const_cond.items():
         # Avoid retrieving RelativeGain, if not needed for correction.
         if cname == "RelativeGain" and not relative_gain:
             const_data[cname] = None
         else:
             const_data[cname] = get_constant_from_db(
                 karabo_id=karabo_id,
                 karabo_da=karabo_da,
                 constant=getattr(Constants.ePix100, cname)(),
                 condition=condition,
                 empty_constant=None if cname == "RelativeGain" else empty_constant,
                 cal_db_interface=cal_db_interface,
                 creation_time=creation_time,
                 print_once=2,
                 timeout=cal_db_timeout
     )
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 if relative_gain and const_data.get("RelativeGain", None) is None:
     print(
         "WARNING: RelativeGain map is requested, but not found.\n"
         "No gain correction will be applied"
     )
     relative_gain = False
     absolute_gain = False
 
 # Initializing some parameters.
 hscale = 1
 stats = True
 hrange = np.array([-50, 1000])
 nbins = hrange[1] - hrange[0]
 commonModeBlockSize = [x//2, y//2]
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 histCalOffsetCor = xcal.HistogramCalculator(
     sensorSize,
     bins=nbins,
     range=hrange,
     parallel=run_parallel,
     nCells=memoryCells,
     blockSize=blockSize
 )
 
-
 # *****************Histogram Calculators****************** #
 histCalCor = xcal.HistogramCalculator(
     sensorSize,
     bins=1050,
     range=[-50, 1000],
     parallel=run_parallel,
     nCells=memoryCells,
     blockSize=blockSize
 )
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 if common_mode:
     histCalCMCor = xcal.HistogramCalculator(
         sensorSize,
         bins=nbins,
         range=hrange,
         parallel=run_parallel,
         nCells=memoryCells,
         blockSize=blockSize,
     )
     cmCorrectionB = xcal.CommonModeCorrection(
         shape=sensorSize,
         blockSize=commonModeBlockSize,
         orientation='block',
         nCells=memoryCells,
         noiseMap=const_data['Noise'],
         runParallel=run_parallel,
         parallel=run_parallel,
         stats=stats,
         minFrac=cm_min_frac,
         noiseSigma=cm_noise_sigma,
     )
     cmCorrectionR = xcal.CommonModeCorrection(
         shape=sensorSize,
         blockSize=commonModeBlockSize,
         orientation='row',
         nCells=memoryCells,
         noiseMap=const_data['Noise'],
         runParallel=run_parallel,
         parallel=run_parallel,
         stats=stats,
         minFrac=cm_min_frac,
         noiseSigma=cm_noise_sigma,
     )
     cmCorrectionC = xcal.CommonModeCorrection(
         shape=sensorSize,
         blockSize=commonModeBlockSize,
         orientation='col',
         nCells=memoryCells,
         noiseMap=const_data['Noise'],
         runParallel=run_parallel,
         parallel=run_parallel,
         stats=stats,
         minFrac=cm_min_frac,
         noiseSigma=cm_noise_sigma,
     )
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 if relative_gain:
     gain_cnst = np.median(const_data["RelativeGain"])
     hscale = gain_cnst
     plot_unit = 'keV'
     if photon_energy > 0:
         plot_unit = '$\gamma$'
         hscale /= photon_energy
 
     gainCorrection = xcal.RelativeGainCorrection(
         sensorSize,
         gain_cnst/const_data["RelativeGain"][..., None],
         nCells=memoryCells,
         parallel=run_parallel,
         blockSize=blockSize,
         gains=None,
     )
 
     histCalRelGainCor = xcal.HistogramCalculator(
         sensorSize,
         bins=nbins,
         range=hrange,
         parallel=run_parallel,
         nCells=memoryCells,
         blockSize=blockSize
     )
 
     if absolute_gain:
         histCalAbsGainCor = xcal.HistogramCalculator(
             sensorSize,
             bins=nbins,
             range=hrange*hscale,
             parallel=run_parallel,
             nCells=memoryCells,
             blockSize=blockSize
         )
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 if pattern_classification :
     patternClassifier = xcal.PatternClassifier(
         [x, y],
         const_data["Noise"],
         split_evt_primary_threshold,
         split_evt_secondary_threshold,
         split_evt_mip_threshold,
         tagFirstSingles=0,
         nCells=memoryCells,
         allowElongated=False,
         blockSize=[x, y],
         parallel=run_parallel,
     )
-    histCalSECor = xcal.HistogramCalculator(
+    histCalCSCor = xcal.HistogramCalculator(
         sensorSize,
         bins=nbins,
         range=hrange,
         parallel=run_parallel,
         nCells=memoryCells,
         blockSize=blockSize,
     )
+    histCalGainCorClusters = xcal.HistogramCalculator(
+        sensorSize,
+        bins=nbins,
+        range=hrange*hscale,
+        parallel=run_parallel,
+        nCells=memoryCells,
+        blockSize=blockSize
+    )
     histCalGainCorSingles = xcal.HistogramCalculator(
         sensorSize,
         bins=nbins,
         range=hrange*hscale,
         parallel=run_parallel,
         nCells=memoryCells,
         blockSize=blockSize
     )
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Applying corrections
 
 %% Cell type:code id: tags:
 
 ``` python
 def correct_train(wid, index, tid, d):
 
     d = d[pixel_data[0]][pixel_data[1]][..., np.newaxis].astype(np.float32)
     d = np.compress(
         np.any(d > 0, axis=(0, 1)), d, axis=2)
 
     # Offset correction.
     d -= const_data["Offset"]
 
     histCalOffsetCor.fill(d)
     # Common Mode correction.
     if common_mode:
         # Block CM
         d = cmCorrectionB.correct(d)
         # Row CM
         d = cmCorrectionR.correct(d)
         # COL CM
         d = cmCorrectionC.correct(d)
         histCalCMCor.fill(d)
 
     # relative gain correction.
     if relative_gain:
         d = gainCorrection.correct(d)
         histCalRelGainCor.fill(d)
 
-    data[index, ...] = np.squeeze(d)
-
     """The gain correction is currently applying
     an absolute correction (not a relative correction
     as the implied by the name);
     it changes the scale (the unit of measurement)
     of the data from ADU to either keV or n_of_photons.
     But the pattern classification relies on comparing
     data with the noise map, which is still in ADU.
 
     The best solution is to do a relative gain
     correction first and apply the global absolute
     gain to the data at the end, after clustering.
     """
 
     if pattern_classification:
 
         d_clu, patterns = patternClassifier.classify(d)
-
         d_clu[d_clu < (split_evt_primary_threshold*const_data["Noise"])] = 0
 
-        data_patterns[index, ...] = np.squeeze(patterns)
-
         data_clu[index, ...] = np.squeeze(d_clu)
+        data_patterns[index, ...] = np.squeeze(patterns)
 
-        d_clu[patterns != 100] = np.nan
-
-        histCalSECor.fill(d_clu)
+        histCalCSCor.fill(d_clu)
 
     # absolute gain correction
     # changes data from ADU to keV (or n. of photons)
     if absolute_gain:
 
         d = d * gain_cnst
         if photon_energy > 0:
             d /= photon_energy
         histCalAbsGainCor.fill(d)
 
         if pattern_classification:
             # Modify pattern classification.
             d_clu = d_clu * gain_cnst
+
             if photon_energy > 0:
                 d_clu /= photon_energy
 
-            histCalGainCorSingles.fill(d_clu)
-
             data_clu[index, ...] = np.squeeze(d_clu)
 
+            histCalGainCorClusters.fill(d_clu)
+
+            d_sing = d_clu[patterns==100] # pattern 100 corresponds to single photons events
+            if len(d_sing):
+                histCalGainCorSingles.fill(d_sing)
+
+    data[index, ...] = np.squeeze(d)
     histCalCor.fill(d)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 pixel_data = (instrument_src, "data.image.pixels")
 
 # 10 is a number chosen after testing 1 ... 71 parallel threads
 context = psh.context.ThreadContext(num_workers=10)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 for f in seq_files:
 
     seq_dc = H5File(f)
 
     n_imgs = seq_dc.get_data_counts(*pixel_data).shape[0]
 
     # Data shape in seq_dc excluding trains with empty images.
     dshape = seq_dc[pixel_data].shape
     dataset_chunk = ((chunk_size_idim,) + dshape[1:])  # e.g. (1, pixels_x, pixels_y)
 
     if n_imgs - dshape[0] != 0:
         print(f"- WARNING: {f} has {n_imgs - dshape[0]} trains with empty data.")
 
     # This parameter is only used for testing.
     if limit_images > 0:
         n_imgs = min(n_imgs, limit_images)
 
     data = context.alloc(shape=dshape, dtype=np.float32)
 
     if pattern_classification:
         data_clu = context.alloc(shape=dshape, dtype=np.float32)
         data_patterns = context.alloc(shape=dshape, dtype=np.int32)
 
     step_timer.start()
 
     context.map(
         correct_train, seq_dc.select(
             *pixel_data, require_all=True).select_trains(np.s_[:n_imgs])
     )
     step_timer.done_step(f'Correcting {n_imgs} trains.')
 
     # Store detector h5 information in the corrected file
     # and deselect data to correct and store later.
     step_timer.start()
 
     out_file = out_folder / f.name.replace("RAW", "CORR")
     data_path = "INSTRUMENT/"+instrument_src+"/data/image"
     pixels_path = f"{data_path}/pixels"
 
     # First copy all raw data source to the corrected file,
     # while excluding the raw data image /data/image/pixels.
     with h5py.File(out_file, 'w') as ofile:
         # Copy RAW non-calibrated sources.
         with h5py.File(f, 'r') as sfile:
             h5_copy_except.h5_copy_except_paths(
                 sfile, ofile,
                 [pixels_path])
 
         # Create dataset in CORR h5 file and add corrected images.
         dataset = ofile.create_dataset(
             pixels_path,
             data=data,
             chunks=dataset_chunk,
             dtype=np.float32)
 
         if pattern_classification:
-            # Save /data/image//pixels_classified in corrected file.
+            # Save /data/image/pixels_classified in corrected file.
             datasetc = ofile.create_dataset(
                 f"{data_path}/pixels_classified",
                 data=data_clu,
                 chunks=dataset_chunk,
                 dtype=np.float32)
 
-            # Save /data/image//patterns in corrected file.
+            # Save /data/image/patterns in corrected file.
             datasetp = ofile.create_dataset(
                 f"{data_path}/patterns",
                 data=data_patterns,
                 chunks=dataset_chunk,
                 dtype=np.int32)
 
         step_timer.done_step('Storing data.')
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 ho, eo, co, so = histCalCor.get()
 
 d = [{
     'x': co,
     'y': ho,
     'y_err': np.sqrt(ho[:]),
     'drawstyle': 'steps-mid',
     'errorstyle': 'bars',
     'errorcoarsing': 2,
     'label': 'Total corr.'
 }]
 
 ho, eo, co, so = histCalOffsetCor.get()
 
 d.append({
     'x': co,
     'y': ho,
     'y_err': np.sqrt(ho[:]),
     'drawstyle': 'steps-mid',
     'errorstyle': 'bars',
     'errorcoarsing': 2,
     'label': 'Offset corr.'
 })
 
 if common_mode:
     ho, eo, co, so = histCalCMCor.get()
     d.append({
         'x': co,
         'y': ho,
         'y_err': np.sqrt(ho[:]),
         'drawstyle': 'steps-mid',
         'errorstyle': 'bars',
         'errorcoarsing': 2,
         'label': 'CM corr.'
     })
 
 if relative_gain :
     ho, eo, co, so = histCalRelGainCor.get()
     d.append({
         'x': co,
         'y': ho,
         'y_err': np.sqrt(ho[:]),
         'drawstyle': 'steps-mid',
         'errorstyle': 'bars',
         'errorcoarsing': 2,
         'label': 'Relative gain corr.'
     })
 
-
 if pattern_classification:
-    ho, eo, co, so = histCalSECor.get()
+    ho, eo, co, so = histCalCSCor.get()
     d.append({
         'x': co,
         'y': ho,
         'y_err': np.sqrt(ho[:]),
         'drawstyle': 'steps-mid',
         'errorstyle': 'bars',
         'errorcoarsing': 2,
-        'label': 'Isolated photons (singles)'
+        'label': 'Charge sharing corr.'
     })
 
 fig = xana.simplePlot(
     d, aspect=1, x_label=f'Energy (ADU)',
     y_label='Number of occurrences', figsize='2col',
     y_log=True, x_range=(-50, 500),
     legend='top-center-frame-2col',
 )
 plt.title(f'run {run} - {karabo_da}')
 plt.grid()
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 if absolute_gain :
     d=[]
     ho, eo, co, so = histCalAbsGainCor.get()
     d.append({
         'x': co,
         'y': ho,
         'y_err': np.sqrt(ho[:]),
         'drawstyle': 'steps-mid',
         'errorstyle': 'bars',
         'errorcoarsing': 2,
         'label': 'Absolute gain corr.'
     })
 
     if pattern_classification:
+        ho, eo, co, so = histCalGainCorClusters.get()
+        d.append({
+            'x': co,
+            'y': ho,
+            'y_err': np.sqrt(ho[:]),
+            'drawstyle': 'steps-mid',
+            'errorstyle': 'bars',
+            'errorcoarsing': 2,
+            'label': 'Charge sharing corr.'
+        })
+
         ho, eo, co, so = histCalGainCorSingles.get()
         d.append({
             'x': co,
             'y': ho,
             'y_err': np.sqrt(ho[:]),
             'drawstyle': 'steps-mid',
             'errorstyle': 'bars',
             'errorcoarsing': 2,
             'label': 'Isolated photons (singles)'
         })
 
     fig = xana.simplePlot(
         d, aspect=1, x_label=f'Energy ({plot_unit})',
         y_label='Number of occurrences', figsize='2col',
         y_log=True,
         x_range=np.array((-50, 500))*hscale,
         legend='top-center-frame-2col',
     )
     plt.grid()
     plt.title(f'run {run} - {karabo_da}')
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Mean Image of the corrected data
 
 %% Cell type:code id: tags:
 
 ``` python
 step_timer.start()
 fig = xana.heatmapPlot(
     np.nanmedian(data, axis=0),
     x_label='Columns', y_label='Rows',
     lut_label=f'Signal ({plot_unit})',
     x_range=(0, y),
     y_range=(0, x),
     vmin=-50, vmax=50)
 step_timer.done_step(f'Plotting mean image of {data.shape[0]} trains.')
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Single Shot of the corrected data
 
 %% Cell type:code id: tags:
 
 ``` python
 step_timer.start()
 fig = xana.heatmapPlot(
     data[0, ...],
     x_label='Columns', y_label='Rows',
     lut_label=f'Signal ({plot_unit})',
     x_range=(0, y),
     y_range=(0, x),
     vmin=-50, vmax=50)
 step_timer.done_step(f'Plotting single shot of corrected data.')
 ```

setup.py

@@ -24,12 +24,6 @@ ext_modules = [
                             '-ftree-vectorize', '-frename-registers'],
         extra_link_args=['-fopenmp'],
     ),
-    Extension(
-        "cal_tools.jfalgs",
-        ['src/cal_tools/jfalgs.pyx'],
-        extra_compile_args=['-O3', '-march=native', '-ftree-vectorize',
-                            '-frename-registers']
-    ),
     Extension(
         "cal_tools.gotthard2.gotthard2algs",
         ["src/cal_tools/gotthard2/gotthard2algs.pyx"],
@@ -65,7 +59,7 @@ install_requires = [
         "astcheck==0.2.5",
         "astsearch==0.2.0",
         "cfelpyutils==2.0.6",
-        "calibration_client==10.0.0",
+        "calibration_client==11.2.0",
         "dill==0.3.0",
         "docutils==0.17.1",
         "dynaconf==3.1.4",

src/cal_tools/jfalgs.pyx

-
-from cython cimport boundscheck, wraparound, cdivision
-from cython.view cimport contiguous
-
-import numpy as np
-
-
-ctypedef fused jf_data_t:
-    unsigned short  # raw pixel data
-    float  # corrected pixel data
-    unsigned int  # mask data
-    unsigned char  # gain data
-
-
-DEF STRIXEL_Y = 86
-DEF STRIXEL_X = 1024 * 3 + 18
-
-
-def strixel_shape():
-    return STRIXEL_Y, STRIXEL_X
-
-
-def strixel_double_pixels():
-    """Build index arrays for double-size pixels.
-
-    In raw data, the entire columns 255, 256, 511, 512, 767 and 768
-    are double-size pixels. After strixelation, these end up in columns
-    765-776, 1539-1550 and 2313-2324 on rows 0-85 or 0-83, with a set
-    of four columns with 86 rows followed by a set of 84 and 86 again.
-
-    This function builds the index arrays after strixelation.
-    """
-
-    Ydouble = []
-    Xdouble = []
-
-    for double_col in [765, 1539, 2313]:
-        for col in range(double_col, double_col+12):
-            for row in range(84 if ((col-double_col) // 4) == 1 else 86):
-                Ydouble.append(row)
-                Xdouble.append(col)
-
-    return np.array(Ydouble), np.array(Xdouble)
-
-
-@boundscheck(False)
-@wraparound(False)
-@cdivision(True)
-def strixel_transform(jf_data_t[:, :, ::contiguous] data,
-                      jf_data_t[:, :, ::contiguous] out = None):
-    """Reorder raw data to physical strixel sensor layout. """
-
-    if data.shape[1] < 256 or data.shape[2] < 256:
-        raise ValueError('Pixel shape of data may not be below (256, 256)')
-
-    if out is None:
-        import numpy as np
-        out = np.zeros((data.shape[0], STRIXEL_Y, STRIXEL_X), dtype=np.float32)
-    elif data.shape[0] > out.shape[0]:
-        raise ValueError('Cell shape of data exceeds out')
-    elif out.shape[1] < STRIXEL_Y or out.shape[2] < STRIXEL_X:
-        raise ValueError(f'Pixel shape of out may not be below '
-                         f'({STRIXEL_Y}, {STRIXEL_X})')
-
-    cdef int cell, yin, xin, xout, yout, igap
-
-    for cell in range(data.shape[0]):
-        # Normal pixels.
-        for yin in range(256):
-            yout = yin // 3
-
-            for xin in range(1024) :
-                xout = 774 * (xin // 256) + 3 * (xin % 256) + yin % 3
-                out[cell, yout, xout] = data[cell, yin, xin]
-
-        # Gap pixels.
-        for yin in range(256):
-            yout = 2 * (yin // 6)
-
-            for igap in range(3) :
-                # Left side of the gap.
-                xin = igap * 256 + 255
-                xout = igap * 774 + 765 + yin % 6
-                out[cell, yout, xout] = data[cell, yin, xin]
-                out[cell, yout+1, xout] = data[cell, yin, xin]
-
-                # Right side of the gap.
-                xin = igap * 256 + 255 + 1
-                xout = igap * 774 + 765 + 11 - yin % 6
-                out[cell, yout, xout] = data[cell, yin, xin]
-                out[cell, yout+1, xout] = data[cell, yin, xin]
-
-    return out

src/cal_tools/jfstrixel.py

+
+import numpy as np
+
+
+REGULAR_SHAPE = (512, 1024)
+STRIXEL_SHAPE = (86, 3090)
+
+
+def _normal_indices():
+    """Build normal size pixel indices."""
+    
+    # Normal pixels
+    yin = np.arange(256)
+    xin = np.arange(1024)
+
+    Yin, Xin = np.meshgrid(yin, xin)
+    Yout, Xout = np.meshgrid(yin // 3, (xin // 256 * 774) + (xin % 256) * 3)
+    Xout += (yin % 3).astype(int)[None, :]
+    
+    return Yout, Xout, Yin, Xin
+
+
+def _gap_indices(in_gap_offset=0, out_gap_offset=0,
+                 xout_factor=+1, yout_offset=0):
+    """Build one half of double size gap pixel indices."""
+    
+    igap = np.arange(3)
+    yin = np.arange(256)
+
+    Yin, Xin = np.meshgrid(yin, igap * 256 + 255 + in_gap_offset)
+    Yout, Xout = np.meshgrid(yin // 6 * 2, igap * 774 + 765 + out_gap_offset)
+    Xout += xout_factor * (yin % 6).astype(int)[None, :]
+    Yout += yout_offset
+    
+    return Yout, Xout, Yin, Xin
+
+
+def transformation_indices2d():
+    """Build 2D strixel transformation index arrays."""
+    
+    # Each of this index sets contains four 2D index arrays
+    # Yout, Xout, Yin, Xin from different parts constituting the full
+    # strixel frame. They are each concatenated across these parts into
+    # four final index arrays to be used for translating between the
+    # regular frame and the strixel frame.
+    index_sets = [
+        _normal_indices(),
+        
+        # Left gap
+        _gap_indices(0, 0, +1, 0), _gap_indices(0, 0, +1, 1),
+        
+        # Right gap
+        _gap_indices(1, 11, -1, 0), _gap_indices(1, 11, -1, 1)
+    ]
+    
+    # Yout, Xout, Yin, Xin
+    # Casting to int64 improves indexing performance by up to 30%.
+    return [np.concatenate(index_set).astype(np.int64)
+            for index_set in zip(*index_sets)]
+
+
+def transformation_indices1d():
+    """Build 1D strixel transformation index arrays.
+    
+    Internally this function reduces the 2D index arrays to a single
+    dimension to operate on raveled data arrays. This improves the
+    transformation performance substantially by up to 3x.
+    """
+
+    Yout, Xout, Yin, Xin = transformation_indices2d()
+     
+    regular_pixel_idx = np.arange(np.prod(REGULAR_SHAPE), dtype=np.uint32) \
+        .reshape(REGULAR_SHAPE)
+    strixel_pixel_idx = np.empty(STRIXEL_SHAPE, dtype=np.int64)
+    strixel_pixel_idx.fill(-1)
+    strixel_pixel_idx[Yout, Xout] = regular_pixel_idx[Yin, Xin]
+
+    Iout = np.where(strixel_pixel_idx.ravel() != -1)[0].astype(np.int64)
+    Iin = strixel_pixel_idx.ravel()[Iout].astype(np.int64)
+    
+    return Iout, Iin
+
+
+def double_pixel_indices():
+    """Build index arrays for double-size pixels.
+
+    In raw data, the entire columns 255, 256, 511, 512, 767 and 768
+    are double-size pixels. After strixelation, these end up in columns
+    765-776, 1539-1550 and 2313-2324 on rows 0-85 or 0-83, with a set
+    of four columns with 86 rows followed by a set of 84 and 86 again.
+
+    This function builds the index arrays for double pixels after
+    strixelation.
+
+    Returns: 
+        (ndarray, ndarray) 2D index arrays for double pixel Y and X.
+    """
+
+    Ydouble = []
+    Xdouble = []
+
+    for double_col in [765, 1539, 2313]:
+        for col in range(double_col, double_col+12):
+            for row in range(84 if ((col-double_col) // 4) == 1 else 86):
+                Ydouble.append(row)
+                Xdouble.append(col)
+
+    return np.array(Ydouble), np.array(Xdouble)
+
+
+def to_strixel(data, out=None):
+    """Transform from regular to strixel geometry.
+
+    Only the last two axes are considered for transformation, input data
+    may have any number of additional axes in front.
+    
+    Args:
+        data (array_like): Data in regular geometry.
+        out (array_like, optional): Buffer for transformed output, a new
+            one is allocated if omitted. Must match all non-frame axes
+            of input data and able to hold strixel frame.
+
+    Returns:
+        (array_like) Data in strixel geometry.
+    """
+
+    if out is None:
+        out = np.zeros((*data.shape[:-2], *STRIXEL_SHAPE), dtype=data.dtype)
+
+    out.reshape(*out.shape[:-2], -1)[..., Iout] = data.reshape(
+        *data.shape[:-2], -1)[..., Iin]
+
+    return out
+
+
+def from_strixel(data, out=None):
+    """Transform from strixel to regular geometry.
+
+    Only the last two axes are considered for transformation, input data
+    may have any number of additional axes in front.
+
+    Args:
+        data (array_like): Data in strixel geometry.
+        out (array_like, optional): Buffer for transformed output, a new
+            one is allocated if omitted. Must match all non-frame axes
+            of input data and able to hold regular frame.
+
+    Returns:
+        (array_like): Data in regular geometry.
+    """
+
+    if out is None:
+        out = np.zeros((*data.shape[:-2], *REGULAR_SHAPE), dtype=data.dtype)
+
+    out.reshape(*out.shape[:-2], -1)[..., Iin] = data.reshape(
+        *data.shape[:-2], -1)[..., Iout]
+
+    return out
+
+
+Iout, Iin = transformation_indices1d()