small refactors

45592638 · Karim Ahmed · 24d22c47 · 45592638
Commit 45592638 authored 3 years ago by Karim Ahmed
--- a/notebooks/AGIPD/Characterize_AGIPD_Gain_FlatFields_Summary.ipynb
+++ b/notebooks/AGIPD/Characterize_AGIPD_Gain_FlatFields_Summary.ipynb
@@ -591,10 +591,10 @@
    "    run_data = RunDirectory(run_folder, include)\n",
    "    if tid:\n",
    "        tid, data = run_data.select('*/DET/*', source).train_from_id(tid)\n",
-    "        return tid, stack_detector_data(data, source,modules=nmods)\n",
+    "        return tid, stack_detector_data(data, source, modules=nmods)\n",
    "    else:\n",
    "        for tid, data in run_data.select('*/DET/*', source).trains(require_all=True):\n",
-    "            return tid, stack_detector_data(data, source,modules=nmods)\n",
+    "            return tid, stack_detector_data(data, source, modules=nmods)\n",
    "    return None, None\n",
    "\n",
    "\n",
@@ -655,8 +655,6 @@
    "fig = plt.figure(figsize=(20,20))\n",
    "ax = fig.add_subplot(111)\n",
    "cdata = np.nanmean(corrected, axis=0)\n",
-    "#Remove me: the scale for comparision should be the same as for orig data\n",
-    "#vmin, vmax = get_range(cdata, 5)\n",
    "ax = geom.plot_data_fast(cdata, ax=ax, cmap=\"jet\", vmin=vmin, vmax=vmax, figsize=(20,20))\n",
    "_ = ax.set_title(\"Corrected data, mean across one train\")"
   ]

 %% 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
 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
 receiver_id = "{}CH0" # inset for receiver devices
 path_template = 'RAW-R{:04d}-{}-S{:05d}.h5' # the template to use to access data
 h5path = 'INSTRUMENT/{}/DET/{}:xtdf/' # path in the HDF5 file to images
 h5path_idx = 'INDEX/{}/DET/{}:xtdf/' # path in the HDF5 file to images
 h5path_ctrl = '/CONTROL/{}/MDL/FPGA_COMP' # path to control information
 karabo_id_control = "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,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
 max_cells = 0 # number of memory cells used, set to 0 to automatically infer
 bias_voltage = 300 # Bias voltage
 acq_rate = 0. # the detector acquisition rate, use 0 to try to auto-determine
 gain_setting = 0.1 # the gain setting, use 0.1 to try to auto-determine
 photon_energy = 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 (
    get_acq_rate,
    get_bias_voltage,
    get_gain_setting,
    get_integration_time,
    get_num_cells,
 )
 from cal_tools.agipdutils_ff import (
    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, BadPixelsFF
 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 RunDirectory, stack_detector_data
 from extra_geom import AGIPD_1MGeometry, AGIPD_500K2GGeometry
 from iCalibrationDB import Conditions, Constants, Detectors
 from iminuit import Minuit
 from IPython.display import HTML, Latex, Markdown, display
 from XFELDetAna.plotting.heatmap import heatmapPlot
 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
 filename = glob.glob(f"{raw_folder}/r{run:04d}/*-AGIPD[0-1][0-9]-*")[0]
 channel = int(re.findall(r".*-AGIPD([0-9]+)-.*", filename)[0])
 control_fname = f'{raw_folder}/r{run:04d}/RAW-R{run:04d}-{karabo_da_control}-S00000.h5'
 h5path_ctrl = h5path_ctrl.format(karabo_id_control)
 # Evaluate number of memory cells
 mem_cells = get_num_cells(filename, karabo_id, channel)
 if mem_cells is None:
    raise ValueError(f"No raw images found in {filename}")
 # Evaluate aquisition rate
 fast_paths = (filename, karabo_id, channel)
 slow_paths = (control_fname, karabo_id_control)
 if acq_rate == 0.:
    acq_rate = get_acq_rate(fast_paths,slow_paths)
 # Evaluate creation time
 creation_time = None
 if use_dir_creation_date:
    creation_time = get_dir_creation_date(raw_folder, run)
 # Evaluate gain setting
 if gain_setting == 0.1:
    if creation_time.replace(tzinfo=None) < parser.parse('2020-01-31'):
        print("Set gain-setting to None for runs taken before 2020-01-31")
        gain_setting = None
    else:
        try:
            gain_setting = get_gain_setting(control_fname, h5path_ctrl)
        except Exception as e:
            print(f'Error while reading gain setting from: \n{control_fname}')
            print(e)
            print("Set gain settion to 0")
            gain_setting = 0
 # Evaluate integration time
 if integration_time < 0:
    integration_time = get_integration_time(control_fname, h5path_ctrl)
 # 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(out_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)
+        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 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)
-#Remove me: the scale for comparision should be the same as for orig data
-#vmin, vmax = get_range(cdata, 5)
 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)))
 ```