[GH2][CORRECT][DARK] Feat/add support for gh2 25um

Gotthard-II 25um Correction and Dark

Description

This Merge request is related to : https://git.xfel.eu/calibration/planning/-/issues/164

How Has This Been Tested?

• Tested against reference runs for GH2 50um.
• Tested against the production database for GH2 25um.

Relevant Documents (optional)

correct_gh2_25um.pdf

gh2_dark_25um.pdf

Correcting:

Into:

Types of changes

• New feature (non-breaking change which adds functionality)
• Bug fix (non-breaking change which fixes an issue)

Checklist:

• Follow up on learned lessons from LPDMini and propagate it here.
• Update the MR's history.
• Get Feedback on results from @mramilli

Reviewers

@schmidtp @hammerd

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
 
 # Bad-pixel thresholds
 d0_lim = [10, 70] # hard limits for d0 value (distance between noise and first peak)
 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
 n_peaks_fit = 4 # Number of gaussian peaks to fit including noise peak
 
 # 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 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 extra_data import RunDirectory, stack_detector_data
 from extra_geom import AGIPD_1MGeometry, AGIPD_500K2GGeometry
 from iCalibrationDB import Conditions, Constants
 from iminuit import Minuit
 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}, {err}")
 ```
 
 %% 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
 
-if mem_cells == 352:
-    for p in condition.parameters:
-        if p.name == 'Memory cells':
-            p.lower_deviation = 352
+# Add a deviation to enable constant retrieval
+# for runs with operating memory cells less than or
+# equal to the memory cells of this injected SlopesFF constant.
+for p in condition.parameters:
+    if p.name == 'Memory cells':
+        p.lower_deviation = mem_cells
 
 # 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)))
 ```