[EPIX100] [DARK]Adding Bad Pixels Map

Description

• Generates the bad pixel map.
• Added badpixel_threshold_sigma input variable: Number of standard deviations considered for bad pixel classification.
• Sensor size is read from hdf5 file instead of being manually input.

How Has This Been Tested?

• Running the notebook
• Terminal command "xfel-calibrate EPIX100 DARK"

Types of changes

• New feature (non-breaking change which adds functionality)

Checklist:

Reviewers

@ahmedk @danilevc

notebooks/ePix100/Characterize_Darks_ePix100_NBC.ipynb

 %% Cell type:markdown id: tags:
 
 # ePix100 Dark Characterization
 
 Author: European XFEL Detector Group, Version: 2.0
 
 The following notebook provides dark image analysis and calibration constants of the ePix100 detector.
 
 Dark characterization evaluates offset and noise of the detector and gives information about bad pixels. Resulting maps are saved as .h5 files for a latter use and injected to calibration DB.
 
 %% Cell type:code id: tags:
 
 ``` python
 in_folder = '/gpfs/exfel/exp/HED/202230/p900247/raw' # input folder, required
 out_folder = '' # output folder, required
 sequence = 0 # sequence file to use
 run = 45 # which run to read data from, required
 
 # Parameters for accessing the raw data.
 karabo_id = "HED_IA1_EPX100-1" # karabo karabo_id
 karabo_da = ["EPIX01"]  # data aggregators
 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 for the calibration database.
 use_dir_creation_date = True
 cal_db_interface = "tcp://max-exfl016:8020" # calibration DB interface to use
 cal_db_timeout = 300000 # timeout on caldb requests
 db_output = False # Output constants to the calibration database
-local_output = True # output constants locally
+local_output = True # Output constants locally
 
 # Conditions used for injected calibration constants.
 bias_voltage = 200 # Bias voltage
 in_vacuum = False # Detector operated in vacuum
-integration_time = -1 # Integration time. Set to -1 to read from .h5 file
+fix_integration_time = -1 # Integration time. Set to -1 to read from .h5 file
 fix_temperature = -1 # Fixed temperature in Kelvin. Set to -1 to read from .h5 file
 temp_limits = 5 # Limit for parameter Operational temperature
 badpixel_threshold_sigma = 5.  # Bad pixels defined by values outside n times this std from median
 
 # Parameters used during selecting raw data trains.
 min_trains = 1 # Minimum number of trains that should be available to process dark constants. Default 1.
 max_trains = 1000  # Maximum number of trains to use for processing dark constants. Set to 0 to use all available trains.
 
 # Don't delete! myMDC sends this by default.
 operation_mode = ''  # Detector operation mode, optional
 
 # TODO: delete after removing from calibration_configurations
 db_module = ''  # ID of module in calibration database, this parameter is ignore in the notebook. TODO: remove from calibration_configurations.
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 import os
 import warnings
 
 import numpy as np
 import pasha as psh
 from extra_data import RunDirectory
 
 import XFELDetAna.xfelprofiler as xprof
 from XFELDetAna import xfelpyanatools as xana
 from XFELDetAna.plotting.util import prettyPlotting
 from cal_tools.enums import BadPixels
 from cal_tools.tools import (
     get_dir_creation_date,
     get_pdu_from_db,
     get_random_db_interface,
     get_report,
     save_const_to_h5,
     send_to_db,
 )
 from iCalibrationDB import Conditions, Constants
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 %matplotlib inline
 
 warnings.filterwarnings('ignore')
 
 prettyPlotting = True
 
 profiler = xprof.Profiler()
 profiler.disable()
 
 instrument_src = instrument_source_template.format(karabo_id, receiver_template)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # Read report path and create file location tuple to add with the injection
 proposal = list(filter(None, in_folder.strip('/').split('/')))[-2]
 file_loc = f'proposal:{proposal} runs:{run}'
 report = get_report(out_folder)
 
 ped_dir = os.path.join(in_folder, f"r{run:04d}")
 fp_name = path_template.format(run, karabo_da[0]).format(sequence)
 run_dir = RunDirectory(ped_dir)
 
 print(f"Run number: {run}")
 print(f"Instrument H5File source: {instrument_src}")
 if use_dir_creation_date:
     creation_time = get_dir_creation_date(in_folder, run)
     print(f"Using {creation_time.isoformat()} as creation time")
 os.makedirs(out_folder, exist_ok=True)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # Read sensor size
 sensor_size = run_dir[instrument_src, 'data.image.dims'].as_single_value(reduce_by='first') # (x=768, y=708) expected
+sensor_size = sensor_size[:2] # data.image.dims for old data is [768, 708, 1]
 assert np.allclose(sensor_size, np.array([768, 708])), 'Unexpected sensor dimensions.'
 
 # Path to pixels ADC values
 pixels_src = (instrument_src, "data.image.pixels")
 
 # Specifies total number of images to proceed
 n_trains = run_dir.get_data_counts(*pixels_src).shape[0]
 
 # Modify n_trains to process based on given maximum
 # and minimun number of trains.
 if max_trains:
     n_trains = min(max_trains, n_trains)
 
 if n_trains < min_trains:
     raise ValueError(
         f"Less than {min_trains} trains are available in RAW data."
          " Not enough data to process darks.")
 
 print(f"Number of dark images to analyze: {n_trains}")
 
 # Read integration time
-if integration_time == -1:
+if fix_integration_time == -1:
     integration_time = run_dir[f"{karabo_id}/DET/CONTROL", 'expTime.value'].as_single_value(reduce_by='first')
     integration_time_str_add = ''
 else:
+    integration_time = fix_integration_time
     integration_time_str_add = '(manual input)'
 
 # Read temperature
 if fix_temperature == -1:
     temperature = np.mean(run_dir.get_array(instrument_src,"data.backTemp").values)/100.
     temperature_k = temperature + 273.15
     temp_str_add = ''
 else:
     temperature_k = fix_temperature
     temperature = fix_temperature - 273.15
     temp_str_add = '(manual input)'
 
 # Print operating conditions
 print(f"Bias voltage: {bias_voltage} V")
 print(f"Detector integration time: {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
 # Calculate noise and offset per pixel and global average, std and median
 data_dc = run_dir.select(
     *pixels_src, require_all=True).select_trains(np.s_[:n_trains])
 print(f"Reading data from: {data_dc.files}\n")
 
 data = data_dc[pixels_src].ndarray()
 
 noise_data = np.std(data, axis=0)
 offset_data = np.mean(data, axis=0)
 
 noise_mean = np.mean(noise_data)
 noise_sigma = np.std(noise_data)
 noise_median = np.median(noise_data)
 noise_min = np.min(noise_data)
 noise_max = np.max(noise_data)
 
 offset_mean = np.mean(offset_data)
 offset_sigma = np.std(offset_data)
 offset_median = np.median(offset_data)
 offset_min = np.min(offset_data)
 offset_max = np.max(offset_data)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # Create and fill offset and noise maps
 constant_maps = {}
 constant_maps['Offset'] = offset_data[..., np.newaxis]
 constant_maps['Noise'] = noise_data[..., np.newaxis]
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Offset ###
 
 %% Cell type:code id: tags:
 
 ``` python
 #************** OFFSET HEAT MAP **************#
 fig = xana.heatmapPlot(
     constant_maps['Offset'][:, :, 0],
     lut_label='[ADU]',
     x_label='Column',
     y_label='Row',
     x_range=(0, sensor_size[0]),
     y_range=(0, sensor_size[1]),
     vmin=max(0, offset_median - badpixel_threshold_sigma*offset_sigma),
     vmax=min(np.power(2,14)-1, offset_median + badpixel_threshold_sigma*offset_sigma)
 )
 fig.suptitle('Offset Map', x=.48, y=.9, fontsize=16)
 fig.set_size_inches(h=15, w=15)
 
 #************** OFFSET HISTOGRAM **************#
 binwidth = offset_sigma/100
 ho, co = np.histogram(
     constant_maps['Offset'].flatten(),
     bins = np.arange(offset_min, offset_max, binwidth)
 )
 do = {'x': co[:-1],
       'y': ho,
       'y_err': np.sqrt(ho[:]),
       'drawstyle': 'bars',
       'color': 'cornflowerblue'}
 
 fig = xana.simplePlot(
     do,
     aspect=1.5,
     x_label='Offset [ADU]',
     y_label='Counts',
     x_range=(0, np.power(2,14)-1),
     y_range=(0, max(ho)*1.1),
     y_log=True
 )
 fig.suptitle('Offset Distribution', x=.5,y =.92, fontsize=16)
 
 stats_str = (
     f'mean: {np.round(offset_mean,2)}\n'
     f'std : {np.round(offset_sigma,2)}\n'
     f'median: {np.round(offset_median,2)}\n'
     f'min: {np.round(offset_min,2)}\n'
     f'max: {np.round(offset_max,2)}'
 )
 fig.text(
     s=stats_str,
     x=.7,
     y=.7,
     fontsize=14,
     bbox=dict(facecolor='yellow', edgecolor='black', alpha=.1));
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Noise ###
 
 %% Cell type:code id: tags:
 
 ``` python
 #************** NOISE HEAT MAP **************#
 fig = xana.heatmapPlot(
     constant_maps['Noise'][:, :, 0],
     lut_label='[ADU]',
     x_label='Column',
     y_label='Row',
     x_range=(0, sensor_size[0]),
     y_range=(0, sensor_size[1]),
     vmin=max(0, noise_median - badpixel_threshold_sigma*noise_sigma),
     vmax=noise_median + badpixel_threshold_sigma*noise_sigma
 )
 fig.suptitle('Noise Map', x=.5, y=.9, fontsize=16)
 fig.set_size_inches(h=15, w=15)
 
 #************** NOISE HISTOGRAM **************#
 binwidth = noise_sigma/100
 hn, cn = np.histogram(constant_maps['Noise'].flatten(),
                       bins=np.arange((min(constant_maps['Noise'].flatten())),
                                  max(constant_maps['Noise'].flatten()) + binwidth,
                                  binwidth))
 dn = {'x': cn[:-1],
       'y': hn,
       'y_err': np.sqrt(hn[:]),
       'drawstyle': 'bars',
       'color': 'cornflowerblue'}
 
 fig = xana.simplePlot(
     dn,
     aspect=1.5,
     x_label='Noise [ADU]',
     y_label='Counts',
     x_range=(max(0, noise_median - badpixel_threshold_sigma*noise_sigma),
              noise_median + badpixel_threshold_sigma*noise_sigma),
     y_range=(0, max(hn)*1.1),
     y_log=True
 )
 fig.suptitle('Noise Distribution',x=.5,y=.92,fontsize=16);
 
 stats_str = (
     f'mean: {np.round(noise_mean,2)}\n'
     f'std: {np.round(noise_sigma,2)}\n'
     f'median: {np.round(noise_median,2)}\n'
     f'min: {np.round(noise_min,2)}\n'
     f'max: {np.round(noise_max,2)}'
 )
 fig.text(
     s=stats_str,
     x=.7,
     y=.7,
     fontsize=14,
     bbox=dict(facecolor='yellow', edgecolor='black', alpha=.1));
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Bad Pixels ###
 
 The bad pixel map is deduced by comparing offset and noise of each pixel ($v_i$) against the median value of the respective maps ($v_k$):
 
 $$
 v_i > \mathrm{median}(v_{k}) + n \sigma_{v_{k}}
 $$
 or
 $$
 v_i < \mathrm{median}(v_{k}) - n \sigma_{v_{k}}
 $$
 
 Values are encoded in a 32 bit mask, where for the dark image deduced bad pixels the following non-zero entries are relevant:
 
 %% Cell type:code id: tags:
 
 ``` python
 def print_bp_entry(bp):
     '''
     Prints bad pixels bit encoding.
 
     Parameters
     ----------
     bp : enum 'BadPixels'
     '''
     print('{:<30s} {:032b}'.format(bp.name, bp.value))
 
 print_bp_entry(BadPixels.OFFSET_OUT_OF_THRESHOLD)
 print_bp_entry(BadPixels.NOISE_OUT_OF_THRESHOLD)
 print_bp_entry(BadPixels.OFFSET_NOISE_EVAL_ERROR)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 def eval_bpidx(const_map, n_sigma):
     '''
     Evaluates bad pixels by comparing offset and noise of each
     pixel against the median value of the respective maps.
 
     Returns boolean array
 
     Parameters
     ----------
     const_map : ndarray
                 Offset or noise constant map to input.
     n_sigma : float
               Standard deviation multiplicity interval outside
               which bad pixels are defined.
     '''
 
     mdn = np.nanmedian(const_map)
     std = np.nanstd(const_map)
     idx = ( (const_map > mdn + n_sigma*std)
           | (const_map < mdn - n_sigma*std) )
 
     return idx
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # Add BadPixels to constant_maps
 constant_maps['BadPixelsDark'] = np.zeros(constant_maps['Offset'].shape, np.uint32)
 
 # Noise related bad pixels
 constant_maps['BadPixelsDark'][eval_bpidx(constant_maps['Noise'], badpixel_threshold_sigma)] = BadPixels.NOISE_OUT_OF_THRESHOLD.value
 constant_maps['BadPixelsDark'][~np.isfinite(constant_maps['Noise'])] = BadPixels.OFFSET_NOISE_EVAL_ERROR.value
 
 # Offset related bad pixels
 constant_maps['BadPixelsDark'][eval_bpidx(constant_maps['Offset'], badpixel_threshold_sigma)] = BadPixels.OFFSET_OUT_OF_THRESHOLD.value
 constant_maps['BadPixelsDark'][~np.isfinite(constant_maps['Offset'])] = BadPixels.OFFSET_NOISE_EVAL_ERROR.value
 
 num_bad_pixels = np.sum(constant_maps['BadPixelsDark']!=0)
 
 print('Number of Bad Pixels: ' + str(num_bad_pixels) + ' (' + str(np.round(100*num_bad_pixels/(sensor_size[0]*sensor_size[1]),2)) + '%)')
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 #************** BAD PIXELS HEAT MAP **************#
 fig = xana.heatmapPlot(
     np.log2(constant_maps['BadPixelsDark'][:, :, 0]),
     lut_label='Bad pixel bit assinged',
     x_label='Column',
     y_label='Row',
     x_range=(0, sensor_size[0]),
     y_range=(0, sensor_size[1])
 )
 fig.suptitle('Bad Pixels Map', x=.5, y=.9, fontsize=16)
 fig.set_size_inches(h=15, w=15)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 # Save constants to DB
 md = None
 
 for const_name in constant_maps.keys():
     const = getattr(Constants.ePix100, const_name)()
     const.data = constant_maps[const_name].data
 
     # Set the operating condition
     condition = Conditions.Dark.ePix100(
         bias_voltage=bias_voltage,
         integration_time=integration_time,
         temperature=temperature_k,
         in_vacuum=in_vacuum)
 
     for parm in condition.parameters:
         if parm.name == "Sensor Temperature":
             parm.lower_deviation = temp_limits
             parm.upper_deviation = temp_limits
 
     # Get physical detector unit
     db_module = get_pdu_from_db(
         karabo_id=karabo_id,
         karabo_da=karabo_da,
         constant=const,
         condition=condition,
         cal_db_interface=cal_db_interface,
         snapshot_at=creation_time)[0]
 
     # Inject or save calibration constants
     if db_output:
         md = send_to_db(
             db_module=db_module,
             karabo_id=karabo_id,
             constant=const,
             condition=condition,
             file_loc=file_loc,
             report_path=report,
             cal_db_interface=cal_db_interface,
             creation_time=creation_time,
             timeout=cal_db_timeout
         )
     if local_output:
         md = save_const_to_h5(
             db_module=db_module,
             karabo_id=karabo_id,
             constant=const,
             condition=condition,
             data=const.data,
             file_loc=file_loc,
             report=report,
             creation_time=creation_time,
             out_folder=out_folder
         )
         print(f"Calibration constant {const_name} is stored locally at {out_folder} \n")
 
 print("Constants parameter conditions are:\n"
     f"• Bias voltage: {bias_voltage}\n"
     f"• Integration time: {integration_time}\n"
     f"• Temperature: {temperature_k}\n"
     f"• In Vacuum: {in_vacuum}\n"
     f"• Creation time: {md.calibration_constant_version.begin_at if md is not None else creation_time}\n")
 ```