[EPIX][CORR] Optimize histograms and plots

Description

Optimization of plots in ePix Correction Notebook, as stated in https://git.xfel.eu/calibration/planning/-/issues/167

How Has This Been Tested?

Run 55 from p003346 (EPIX-1 and EPIX-2).

Relevant Documents (optional)

Comparision between the generated reports before/after this MR, where the histpgrams are optimized and the plots now show meaningfull information:

OLD: MID_EXP_EPIX-1_correct_003346_r55_230330_195859.pdf

NEW:EPIX100CORRECTCalibration_new.pdf

Types of changes

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

Specifically:

• Using np.histogram instead of xcal.HistogramCalculator to calculate histograms
• Using extra_geom instead of xana.heatmapPlot
• Dynamic colorbar limits, replacing the fixed ones
• Added additional plot for cluster corrected data
• Added meaningfull titles to the plots

Reviewers

@ahmedk

notebooks/ePix100/Correction_ePix100_NBC.ipynb

 ``` python
-in_folder = "/gpfs/exfel/exp/HED/202202/p003121/raw" # input folder, required
+in_folder = "/gpfs/exfel/exp/MID/202301/p003346/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 = 156  # which run to read data from, required
+run = 55  # which run to read data from, required
 
 # Parameters for accessing the raw data.
-karabo_id = "HED_IA1_EPX100-1"  # karabo karabo_id
+karabo_id = "MID_EXP_EPIX-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
 import matplotlib.pyplot as plt
 from IPython.display import Latex, Markdown, display
 from extra_data import RunDirectory, H5File
+from extra_geom import Epix100Geometry
+from mpl_toolkits.axes_grid1 import make_axes_locatable
 from pathlib import Path
 
 import cal_tools.restful_config as rest_cfg
-from XFELDetAna import xfelpyanatools as xana
 from XFELDetAna import xfelpycaltools as xcal
 from cal_tools.calcat_interface import EPIX100_CalibrationData
 from cal_tools.epix100 import epix100lib
 from cal_tools.files import DataFile
+from cal_tools.restful_config import restful_config
 from cal_tools.tools import (
     calcat_creation_time,
     write_constants_fragment,
 
 prettyPlotting = True
 
+
 %matplotlib inline
 ``` python
+prop_str = in_folder[in_folder.find('/p')+1:in_folder.find('/p')+8]
+
 in_folder = Path(in_folder)
 out_folder = Path(out_folder)
 
 %% Cell type:markdown id: tags:
 ## Retrieving calibration constants
 
-As a first step, calibration constants have to be loaded.
+As a first step, dark maps have to be loaded.
 %% Cell type:code id: tags:
 # Initializing some parameters.
 hscale = 1
 stats = True
-hrange = np.array([-50, 1000])
-nbins = hrange[1] - hrange[0]
-commonModeBlockSize = [x//2, y//2]
+bins = np.arange(-50,1000)
+hist = {'O': 0} # dictionary to store histograms
 %% 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
-)
+if common_mode:
 
-%% Cell type:code id: tags:
+    commonModeBlockSize = [x//2, y//8]
 
-``` 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,
         stats=stats,
         minFrac=cm_min_frac,
         noiseSigma=cm_noise_sigma,
     )
+
+    hist['CM'] = 0
 %% Cell type:code id: tags:
 ``` python
 if relative_gain:
-    gain_cnst = np.median(const_data["RelativeGainEPix100"])
-    hscale = gain_cnst
-    plot_unit = 'keV'
-    if photon_energy > 0:
-        plot_unit = '$\gamma$'
-        hscale /= photon_energy
+
+    # gain constant is given by the mode of the gain map
+    # because all bad pixels are masked using this value
+    _vals,_counts = np.unique(const_data["RelativeGainEPix100"], return_counts=True)
+    gain_cnst = _vals[np.argmax(_counts)]
 
     gainCorrection = xcal.RelativeGainCorrection(
         sensorSize,
         blockSize=blockSize,
         gains=None,
     )
 
-    histCalRelGainCor = xcal.HistogramCalculator(
-        sensorSize,
-        bins=nbins,
-        range=hrange,
-        parallel=run_parallel,
-        nCells=memoryCells,
-        blockSize=blockSize
-    )
+    hist['RG'] = 0
 
     if absolute_gain:
-        histCalAbsGainCor = xcal.HistogramCalculator(
-            sensorSize,
-            bins=nbins,
-            range=hrange*hscale,
-            parallel=run_parallel,
-            nCells=memoryCells,
-            blockSize=blockSize
-        )
+        hscale = gain_cnst
+        plot_unit = 'keV'
+        if photon_energy > 0:
+            plot_unit = '$\gamma$'
+            hscale /= photon_energy
+        hist['AG'] = 0
 %% Cell type:code id: tags:
         blockSize=[x, y],
         parallel=run_parallel,
     )
-    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
-    )
+    hist['CS'] = 0
+    hist['S'] = 0
 %% Cell type:markdown id: tags:
     d = d[..., np.newaxis].astype(np.float32)
     d = np.compress(
         np.any(d > 0, axis=(0, 1)), d, axis=2)
 
     # Offset correction.
     d -= const_data["OffsetEPix100"]
+    hist['O'] += np.histogram(d,bins=bins)[0]
 
-    histCalOffsetCor.fill(d)
     # Common Mode correction.
     if common_mode:
         # Block CM
         d = cmCorrectionR.correct(d)
         # COL CM
         d = cmCorrectionC.correct(d)
-        histCalCMCor.fill(d)
 
-    # relative gain correction.
+        hist['CM'] += np.histogram(d,bins=bins)[0]
+
+
+    # Relative gain correction.
     if relative_gain:
         d = gainCorrection.correct(d)
-        histCalRelGainCor.fill(d)
+
+        hist['RG'] += np.histogram(d,bins=bins)[0]
 
     """The gain correction is currently applying
     an absolute correction (not a relative correction
         data_clu[index, ...] = np.squeeze(d_clu)
         data_patterns[index, ...] = np.squeeze(patterns)
 
-        histCalCSCor.fill(d_clu)
+        hist['CS'] += np.histogram(d_clu,bins=bins)[0]
 
-    # absolute gain correction
+        d_sing = d_clu[patterns==100] # pattern 100 corresponds to single photons events
+        if len(d_sing):
+            hist['S'] += np.histogram(d_sing,bins=bins)[0]
+
+    # 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)
+
+        hist['AG'] += np.histogram(d,bins=bins)[0]
 
         if pattern_classification:
             # Modify pattern classification.
             d_clu = d_clu * gain_cnst
-
             if photon_energy > 0:
                 d_clu /= photon_energy
 
             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:
     exit(0)
+%% Cell type:markdown id: tags:
+## Plot Histograms
+%% Cell type:code id: tags:
+``` python
+bins_ADU = bins[:-1]+np.diff(bins)[0]/2
+bins_keV = bins_ADU*hscale
 %% Cell type:code id: tags:
 ``` python
-ho, eo, co, so = histCalCor.get()
+# Histogram in ADU
 
-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.'
-})
+plt.figure(figsize=(12,8))
+plt.plot(bins_ADU,hist['O'], 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.'
-    })
-
+    plt.plot(bins_ADU,hist['CM'], label='CM corr')
+if relative_gain:
+    plt.plot(bins_ADU,hist['RG'], label='Relative Gain corr')
 if pattern_classification:
-    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': '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()
+    plt.plot(bins_ADU[bins_ADU>10],hist['CS'][bins_ADU>10], label='Charge Sharing corr')
+    if np.any(hist['S']):
+        plt.plot(bins_ADU,hist['S'], label='Singles')
+
+xtick_step = 50
+plt.xlim(bins[0], bins[-1]+1)
+plt.xticks(np.arange(bins[0],bins[-1]+2,xtick_step))
+plt.xlabel('ADU',fontsize=12)
+
+plt.yscale('log')
+plt.title(f'{karabo_id} | {prop_str}, r{run}', fontsize=14, fontweight='bold')
+plt.legend(fontsize=12)
+plt.grid(ls=':')
 %% Cell type:code id: tags:
 ``` python
-if absolute_gain :
-    d = []
-    ho, eo, co, so = histCalAbsGainCor.get()
-    if co is not None:  # avoid adding None array, if calculator is empty.
-        d.append({
-            'x': co,
-            'y': ho,
-            'y_err': np.sqrt(ho[:]),
-            'drawstyle': 'steps-mid',
-            'errorstyle': 'bars',
-            'errorcoarsing': 2,
-            'label': 'Absolute gain corr.'
-        })
+# Histogram in keV/number of photons
+
+if absolute_gain:
+    colors = plt.rcParams['axes.prop_cycle'].by_key()['color']
+    plt.figure(figsize=(12,8))
+
+    if relative_gain:
+        plt.plot(bins_keV,hist['RG'], label='Absolute Gain corr', c=colors[2])
 
     if pattern_classification:
-        ho, eo, co, so = histCalGainCorClusters.get()
-        if co is not None:  # avoid adding None array, if calculator is empty.
-            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()
-        if co is not None:  # avoid adding None array, if calculator is empty.
-            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}')
+        plt.plot(bins_keV[bins_keV>.5],hist['CS'][bins_keV>.5], label='Charge Sharing corr', c=colors[3])
+        if np.any(hist['S']):
+            plt.plot(bins_keV[bins_keV>.5],hist['S'][bins_keV>.5], label='Singles', c=colors[4])
+
+    if photon_energy==0: # if keV instead of #photons
+        xtick_step = 5
+        plt.xlim(left=-2)
+        plt.xticks(np.arange(0,plt.gca().get_xlim()[1],xtick_step))
+    plt.xlabel(plot_unit,fontsize=12)
+
+    plt.yscale('log')
+    plt.title(f'{karabo_id} | {prop_str}, r{run}', fontsize=14, fontweight='bold')
+    plt.legend(fontsize=12)
+    plt.grid(ls=':')
 %% Cell type:markdown 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.')
+geom = Epix100Geometry.from_relative_positions(top=[386.5, 364.5, 0.], bottom=[386.5, -12.5, 0.])
+
+if pattern_classification:
+    plt.subplots(1,2,figsize=(18,18)) if pattern_classification else plt.subplots(1,1,figsize=(9,9))
+    ax = plt.subplot(1,2,1)
+    ax.set_title(f'Before CS correction',fontsize=12,fontweight='bold');
+else:
+    plt.subplots(1,1,figsize=(9,9))
+    ax = plt.subplot(1,1,1)
+    ax.set_title(f'{karabo_id} | {prop_str}, r{run} | Average of {data.shape[0]} trains',fontsize=12,fontweight='bold');
+
+# Average image before charge sharing corrcetion
+divider = make_axes_locatable(ax)
+cax = divider.append_axes('bottom', size='5%', pad=0.5)
+
+image = data.mean(axis=0)
+vmin = max(image.mean()-2*image.std(),0)
+vmax = image.mean()+3*image.std()
+geom.plot_data(image,
+               ax=ax,
+               colorbar={'cax': cax, 'label': plot_unit, 'orientation': 'horizontal'},
+               origin='upper',
+               vmin=vmin,
+               vmax=vmax)
+
+# Average image after charge sharing corrcetion
+if pattern_classification:
+
+    ax = plt.subplot(1,2,2)
+    divider = make_axes_locatable(ax)
+    cax = divider.append_axes('bottom', size='5%', pad=0.5)
+
+    image = data_clu.mean(axis=0)
+    geom.plot_data(image,
+                   ax=ax,
+                   colorbar={'cax': cax, 'label': plot_unit, 'orientation': 'horizontal'},
+                   origin='upper',
+                   vmin=vmin,
+                   vmax=vmax)
+    ax.set_title(f'After CS correction',fontsize=12,fontweight='bold');
+
+    plt.suptitle(f'{karabo_id} | {prop_str}, r{run} | Average of {data.shape[0]} trains',fontsize=14,fontweight='bold',y=.72);
 %% Cell type:markdown 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.')
+train_idx = -1
+
+if pattern_classification:
+    plt.subplots(1,2,figsize=(18,18)) if pattern_classification else plt.subplots(1,1,figsize=(9,9))
+    ax = plt.subplot(1,2,1)
+    ax.set_title(f'Before CS correction',fontsize=12,fontweight='bold');
+else:
+    plt.subplots(1,1,figsize=(9,9))
+    ax = plt.subplot(1,1,1)
+    ax.set_title(f'{karabo_id} | {prop_str}, r{run} | Single frame',fontsize=12,fontweight='bold');
+
+# Average image before charge sharing corrcetion
+divider = make_axes_locatable(ax)
+cax = divider.append_axes('bottom', size='5%', pad=0.5)
+
+image = data[train_idx]
+vmin = max(image.mean()-2*image.std(),0)
+vmax = image.mean()+3*image.std()
+geom.plot_data(image,
+               ax=ax,
+               colorbar={'cax': cax, 'label': plot_unit, 'orientation': 'horizontal'},
+               origin='upper',
+               vmin=vmin,
+               vmax=vmax)
+
+# Average image after charge sharing corrcetion
+if pattern_classification:
+
+    ax = plt.subplot(1,2,2)
+    divider = make_axes_locatable(ax)
+    cax = divider.append_axes('bottom', size='5%', pad=0.5)
+
+    image = data_clu[train_idx]
+    geom.plot_data(image,
+                   ax=ax,
+                   colorbar={'cax': cax, 'label': plot_unit, 'orientation': 'horizontal'},
+                   origin='upper',
+                   vmin=vmin,
+                   vmax=vmax)
+    ax.set_title(f'After CS correction',fontsize=12,fontweight='bold');
+
+    plt.suptitle(f'{karabo_id} | {prop_str}, r{run} | Single frame',fontsize=14,fontweight='bold',y=.72);