reviewed version

5f95bdad · Nuno Duarte · Nuno Duarte · 4e11910c · 5f95bdad
Commit 5f95bdad authored 3 years ago by Nuno Duarte Committed by Nuno Duarte 2 years ago
--- a/notebooks/ePix100/Characterize_Darks_ePix100_NBC.ipynb
+++ b/notebooks/ePix100/Characterize_Darks_ePix100_NBC.ipynb
@@ -202,7 +202,6 @@
    "noise_min = np.min(noise_data)\n",
    "noise_max = np.max(noise_data)\n",
    "\n",
-    "\n",
    "offset_mean = np.mean(offset_data)\n",
    "offset_sigma = np.std(offset_data)\n",
    "offset_median = np.median(offset_data)\n",

 %% 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

 # Conditions used for injected calibration constants.
 bias_voltage = 200 # bias voltage
 in_vacuum = False # detector operated in vacuum
 fix_temperature = 0. # Fixed temperature in Kelvin. Set to 0 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.detectors.fastccd import readerh5 as fastccdreaderh5
 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)
 filename = os.path.join(ped_dir, fp_name)
 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 = np.array(run_dir.get_array(
    instrument_src,
    "data.image.dims")[0,:2]) # (x,y)

 # Path to pixel ADC values
 pixel_data_dir = (instrument_src, "data.image.pixels")

 # Specifies total number of images to proceed
 n_trains = run_dir.get_data_counts(*pixel_data_dir).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}")

 integration_time = int(run_dir.get_array(
    f"{karabo_id}/DET/CONTROL",
    "expTime.value")[0])

 if fix_temperature:
    temperature_k = fix_temperature
    temperature = fix_temperature - 273.15
    print("Temperature is fixed!")
 else:
    temperature = np.mean(run_dir.get_array(
        instrument_src,
        "data.backTemp").values) / 100.
    temperature_k = temperature + 273.15

 print(f"Bias voltage is {bias_voltage} V")
 print(f"Detector integration time is set to {integration_time} \u03BCs")
 print(f"Mean temperature was {temperature:0.1f}°C / {temperature_k:0.1f} K")
 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(
    *pixel_data_dir, require_all=True).select_trains(np.s_[:n_trains])
 print(f"Reading data from: {data_dc.files}\n")

 data = data_dc[pixel_data_dir].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 Map ###

 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['BadPixels'] = np.zeros(constant_maps['Offset'].shape, np.uint32)

 # Noise related bad pixels
 constant_maps['BadPixels'][eval_bpidx(constant_maps['Noise'], badpixel_threshold_sigma)] = BadPixels.NOISE_OUT_OF_THRESHOLD.value
 constant_maps['BadPixels'][~np.isfinite(constant_maps['Noise'])] = BadPixels.OFFSET_NOISE_EVAL_ERROR.value

 # Offset related bad pixels
 constant_maps['BadPixels'][eval_bpidx(constant_maps['Offset'], badpixel_threshold_sigma)] = BadPixels.OFFSET_OUT_OF_THRESHOLD.value
 constant_maps['BadPixels'][~np.isfinite(constant_maps['Offset'])] = BadPixels.OFFSET_NOISE_EVAL_ERROR.value

 num_bad_pixels = np.sum(constant_maps['BadPixels']!=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['BadPixels'][:, :, 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")
 ```