Skip to content
Snippets Groups Projects

Feat/dss cimprove master rebasing

Closed Andrey Samartsev requested to merge feat/DSSCimproveMasterRebasing into feat/DSSCdarksImprove
2 files
+ 55
67
Compare changes
  • Side-by-side
  • Inline
Files
2
%% Cell type:markdown id: tags:
# pnCCD Dark Characterization
Author: Kiana Setoodehnia, DET Group, Version 2.0
The following notebook provides dark image analysis of the pnCCD detector.
Dark characterization evaluates offset and noise of the detector and gives information about bad pixels. The noise map is corrected for common mode. The resulting map together with the offset and bad pixels maps are sent to the database and are also saved as .h5 files to a local path for a later use.
%% Cell type:code id: tags:
``` python
in_folder = "/gpfs/exfel/exp/SQS/201921/p002430/raw" # input folder, required
out_folder = '/gpfs/exfel/exp/SQS/201921/p002430/scratch' # output folder, required
path_template = 'RAW-R{:04d}-PNCCD01-S{{:05d}}.h5' # the template to use to access data
h5path = '/INSTRUMENT/SQS_NQS_PNCCD1MP/CAL/PNCCD_FMT-0:output/data/image/' # path to the data in the HDF5 file
run = 745 # which run to read data from, required
cluster_profile = "noDB" # ipcluster profile to use
in_folder = "/gpfs/exfel/exp/SQS/201921/p002430/raw/" # input folder, required
out_folder = 'gpfs/exfel/data/scratch/karnem/test/' # output folder, required
sequence = 0 # sequence file to use
run = 745 # which run to read data from, required
karabo_da = 'PNCCD01' # data aggregators
karabo_id = "SQS_NQS_PNCCD1MP" # karabo prefix of PNCCD devices
receiver_id = "PNCCD_FMT-0" # inset for receiver devices
path_template = 'RAW-R{:04d}-{}-S{{:05d}}.h5' # the template to use to access data
h5path = '/INSTRUMENT/{}/CAL/{}:output/data/image/' # path in the HDF5 file the data is at
use_dir_creation_date = True # use dir creation date as data production reference date
cal_db_interface = "tcp://max-exfl016:8021" # calibration DB interface to use
cal_db_timeout = 300000 # timeout on caldb requests
db_output = False # if True, the notebook sends dark constants to the calibration database
local_output = True # if True, the notebook saves dark constants locally
number_dark_frames = 0 # number of images to be used, if set to 0 all available images are used
chunkSize = 100 # number of images to read per chunk
fix_temperature = 233. # fix temperature in K, set to -1 to use value from slow data
gain = 0 # the detector's gain setting, only 0 is currently implemented
bias_voltage = 300 # detector's bias voltage
integration_time = 70 # detector's integration time
commonModeAxis = 0 # axis along which common mode will be calculated (0: along rows, 1: along columns)
sigmaNoise = 10. # pixels whose signal value exceeds sigmaNoise*noise value in that pixel will be masked
bad_pixel_offset_sigma = 5. # any pixel whose offset is beyond 5 standard deviations, is a bad pixel
bad_pixel_noise_sigma = 5. # any pixel whose noise is beyond 5 standard deviations, is a bad pixel
cluster_profile = "noDB" # ipcluster profile to use
run_parallel = True # for parallel computation
cpuCores = 40 # specifies the number of running cpu cores
cal_db_interface = "tcp://max-exfl016:8021" # calibration DB interface to use
temp_limits = 5 # temperature limits in which calibration parameters are considered equal
db_output = False # if True, the notebook sends dark constants to the calibration database
local_output = True # if True, the notebook saves dark constants locally
# for database time derivation:
use_dir_creation_date = True # To be used to retrieve calibration constants later on
temp_limits = 5 # temperature limits in which to consider calibration parameters equal
multi_iteration = False # use multiple iterations
```
%% Cell type:code id: tags:
``` python
import os
import copy
import warnings
warnings.filterwarnings('ignore')
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from prettytable import PrettyTable
from IPython.display import display, Markdown
from iCalibrationDB import ConstantMetaData, Constants, Conditions, Detectors, Versions
from iCalibrationDB.detectors import DetectorTypes
import XFELDetAna.xfelprofiler as xprof
profiler = xprof.Profiler()
profiler.disable()
from XFELDetAna.util import env
env.iprofile = cluster_profile
from XFELDetAna import xfelpycaltools as xcal
from XFELDetAna import xfelpyanatools as xana
from XFELDetAna.plotting.util import prettyPlotting
prettyPlotting=True
from XFELDetAna.xfelreaders import ChunkReader
from XFELDetAna.detectors.fastccd import readerh5 as fastccdreaderh5
from cal_tools.tools import get_dir_creation_date, save_const_to_h5, get_random_db_interface
from cal_tools.enums import BadPixels
```
%% Cell type:code id: tags:
``` python
# Output Folder Creation:
if not os.path.exists(out_folder):
os.makedirs(out_folder)
```
%% Cell type:code id: tags:
``` python
def nImagesOrLimit(nImages, limit):
if limit == 0:
return nImages
else:
return min(nImages, limit)
```
%% Cell type:code id: tags:
``` python
proposal = list(filter(None, in_folder.strip('/').split('/')))[-2]
file_loc = 'Proposal: {}, Run: {}'.format(proposal, run)
print(file_loc)
```
%% Cell type:code id: tags:
``` python
# Calibration Database Settings, and Some Initial Run Parameters & Paths:
display(Markdown('### Initial Settings'))
x = 1024 # rows of the FastCCD to analyze in FS mode
y = 1024 # columns of the FastCCD to analyze in FS mode
print("pnCCD size is: {}x{} pixels.".format(x, y))
ped_dir = "{}/r{:04d}".format(in_folder, run)
fp_name = path_template.format(run)
fp_name = path_template.format(run, karabo_da)
fp_path = '{}/{}'.format(ped_dir, fp_name)
filename = fp_path.format(sequence)
h5path = h5path.format(karabo_id, receiver_id)
creation_time = None
if use_dir_creation_date:
creation_time = get_dir_creation_date(in_folder, run)
cal_db_interface = get_random_db_interface(cal_db_interface)
print('Calibration database interface: {}'.format(cal_db_interface))
print("Sending constants to the calibration database: {}".format(db_output))
print("HDF5 path to data: {}".format(h5path))
print("Reading data from: {}".format(filename))
print("Run number: {}".format(run))
if creation_time:
print("Creation time: {}".format(creation_time.isoformat()))
```
%% Cell type:code id: tags:
``` python
# Reading Parameters such as Detector Bias, Gain, etc. from the Data:
memoryCells = 1 # pnCCD has 1 memory cell
sensorSize = [x, y]
blockSize = [sensorSize[0]//2, sensorSize[1]//2]# sensor area will be analysed according to blocksize
xcal.defaultBlockSize = blockSize
nImages = fastccdreaderh5.getDataSize(filename, h5path)[0] # specifies total number of images to proceed
nImages = nImagesOrLimit(nImages, number_dark_frames)
profile = False
```
%% Cell type:code id: tags:
``` python
# Printing the Parameters Read from the Data File:
display(Markdown('### Detector Parameters'))
print("Bias voltage is {} V.".format(bias_voltage))
print("Detector gain is set to {}.".format(gain))
print("Detector integration time is set to {} ms".format(integration_time))
if fix_temperature != 0.:
print("Using a fixed temperature of {} K".format(fix_temperature))
temperature_k = fix_temperature
else:
print("Temperature is not fixed.")
print("Number of dark images to analyze:", nImages)
```
%% Cell type:code id: tags:
``` python
# Reading Files in Chunks:
# Chunk reader returns an iterator to access the data in the file within the ranges:
reader = ChunkReader(filename, fastccdreaderh5.readData, nImages, chunkSize, path = h5path,
pixels_x = sensorSize[0], pixels_y = sensorSize[1],)
```
%% Cell type:code id: tags:
``` python
# Calculators:
# noiseCal is a noise map calculator, which internally also produces a per-pixel mean map, i.e. an offset map:
noiseCal = xcal.NoiseCalculator(sensorSize, memoryCells, cores=cpuCores, blockSize=blockSize,
runParallel=run_parallel)
```
%% Cell type:markdown id: tags:
### First Iteration
Characterization of dark images with purpose to create dark maps (offset, noise and bad pixel maps) is an iterative process. Firstly, initial offset and noise maps are produced from the raw dark data.
%% Cell type:code id: tags:
``` python
counter1 = 1 # To count how many "if data.shape[2] >= chunkSize" instances are there.
counter2 = 0 # To count how many "if data.shape[2] < chunkSize" instances are there.
chunkSize_new = 0 # See below
for data in reader.readChunks():
data = data.astype(np.float32)
dx = np.count_nonzero(data, axis=(0, 1))
data = data[:,:,dx != 0]
# Some sequences may have less than 500 frames in them. To find out how many images there are, we will
# temporarily change chunkSize to be the same as whatever number of frames the last chunk of data has:
if data.shape[2] < chunkSize:
chunkSize_new = data.shape[2]
print("Number of images are less than chunkSize. chunkSize is temporarily changed to {}."
.format(chunkSize_new))
images = images + chunkSize_new
counter2 += 1
else:
images = counter1*chunkSize + counter2*chunkSize_new
counter1 += 1
noiseCal.fill(data) # Filling the histogram calculator with data
chunkSize = 100 # resetting the chunkSize to its default value for the next sequence or data-chunk
print('A total number of {} images are processed.'.format(images))
```
%% Cell type:code id: tags:
``` python
offsetMap = noiseCal.getOffset() # Producing offset map
noiseMap = noiseCal.get() # Producing noise map
noiseCal.reset() # Resetting noise calculator
print("Initial maps are created.")
```
%% Cell type:markdown id: tags:
### Offset and Noise Maps prior to Common Mode Correction
In the following, the histogram of the pnCCD offset, pnCCD offset map, as well as the initial uncorrected noise map are plotted:
%% Cell type:code id: tags:
``` python
#************** OFFSET MAP HISTOGRAM ***********#
ho, co = np.histogram(offsetMap.flatten(), bins=700) # ho = offset histogram; co = offset bin centers
do = {'x': co[:-1],
'y': ho,
'y_err': np.sqrt(ho[:]),
'drawstyle': 'bars',
'color': 'cornflowerblue',
'label': 'Raw Signal (ADU)'
}
fig = xana.simplePlot(do, figsize='1col', aspect=1, x_label = 'Raw Signal (ADU)', y_label="Counts",
title = 'Offset Histogram', x_range=(0, np.nanmax(offsetMap)), y_log=True)
# Calculating mean, median and standard deviation for the above histogram:
mids = 0.5*(co[1:] + co[:-1])
mean = np.average(mids, weights=ho)
variance = np.average((mids - mean)**2, weights=ho)
std = np.sqrt(variance)
# Table of statistics on raw signal:
t0 = PrettyTable()
t0.title = "Statistics on Raw Signal (Offset)"
t0.field_names = ["Mean", "Standard Deviation"]
t0.add_row(["{:0.3f} (ADU)".format(mean), "{:0.3f} (ADU)".format(std)])
print(t0,'\n')
#***** NOISE MAP HISTOGRAM FROM THE OFFSET CORRECTED DATA *******#
hn, cn = np.histogram(noiseMap.flatten(), bins=2000) # hn = noise histogram; cn = noise bin centers
dn = {'x': cn[:-1],
'y': hn,
'y_err': np.sqrt(hn[:]),
'drawstyle': 'bars',
'color': 'cornflowerblue',
'label': 'Noise (ADU)'
}
fig = xana.simplePlot(dn, figsize='1col', aspect=1, x_label = 'Noise (ADU)', y_label="Counts",
title = 'Noise Histogram', y_log=True, x_range=(0, 600))
#************** HEAT MAPS *******************#
fig = xana.heatmapPlot(offsetMap[:,:,0], x_label='Column Number', y_label='Row Number', aspect=1,
x_range=(0, y), y_range=(0, x), vmin=6000, vmax=14500, lut_label='Offset (ADU)',
panel_x_label='Columns Stat (ADU)', panel_y_label='Rows Stat (ADU)', title = 'OffsetMap')
fig = xana.heatmapPlot(noiseMap[:,:,0], x_label='Column Number', y_label='Row Number', aspect=1,
lut_label='Uncorrected Noise (ADU)', x_range=(0, y),
y_range=(0, x), vmax=2*np.mean(noiseMap), panel_x_label='Columns Stat (ADU)',
panel_y_label='Rows Stat (ADU)', title = 'Uncorrected NoiseMap')
```
%% Cell type:markdown id: tags:
### Initial BadPixelMap
This is generated based on the offset and uncorrected noise maps:
%% Cell type:code id: tags:
``` python
bad_pixels = np.zeros(offsetMap.shape, np.uint32)
mnoffset = np.nanmedian(offsetMap)
stdoffset = np.nanstd(offsetMap)
bad_pixels[(offsetMap < mnoffset-bad_pixel_offset_sigma*stdoffset) |
(offsetMap > mnoffset+bad_pixel_offset_sigma*stdoffset)] = BadPixels.OFFSET_OUT_OF_THRESHOLD.value
mnnoise = np.nanmedian(noiseMap)
stdnoise = np.nanstd(noiseMap)
bad_pixels[(noiseMap < mnnoise-bad_pixel_noise_sigma*stdnoise) |
(noiseMap > mnnoise+bad_pixel_noise_sigma*stdnoise)] = BadPixels.NOISE_OUT_OF_THRESHOLD.value
fig = xana.heatmapPlot(np.log2(bad_pixels[:, :, 0]),
x_label='Columns', y_label='Rows',
lut_label='Bad Pixel Value (ADU)',
x_range=(0, y),
y_range=(0, x), vmin=0, vmax=32)
```
%% Cell type:code id: tags:
``` python
# Common Mode Correction:
# In this method, the median of all pixels that are read out at the same time along a row is subtracted from the
# signal in each pixel:
cmCorrection = xcal.CommonModeCorrection(sensorSize, [data.shape[0]//2, data.shape[1]//2],
commonModeAxis, parallel=run_parallel, dType=np.float32, stride=1,
noiseMap=noiseMap.astype(np.float32), minFrac=0)
cmCorrection.debug()
```
%% Cell type:code id: tags:
``` python
# Histogram Calculators:
# For offset corrected data:
histCalCorrected = xcal.HistogramCalculator(sensorSize, bins=400, range=[-200, 200], memoryCells=memoryCells,
cores=cpuCores, gains=None, blockSize=blockSize)
# For common mode corrected data:
histCalCMCorrected = xcal.HistogramCalculator(sensorSize, bins=400, range=[-200, 200], memoryCells=memoryCells,
cores=cpuCores, gains=None, blockSize=blockSize)
noiseCalCM = xcal.NoiseCalculator(sensorSize, memoryCells, cores=cpuCores, blockSize=blockSize,
runParallel=run_parallel)
```
%% Cell type:markdown id: tags:
### Second Iteration
During the second iteration, the data are offset and common mode corrected to produce a common mode corrected noise map. The offset and common mode corrections are calculated by subtracting out the median of all pixels that are read out at the same time along a row.
%% Cell type:code id: tags:
``` python
counter1 = 1 # To count how many "if data.shape[2] >= chunkSize" instances are there.
counter2 = 0 # To count how many "if data.shape[2] < chunkSize" instances are there.
chunkSize_new = 0 # See below
for data in reader.readChunks():
data = data.astype(np.float32)
dx = np.count_nonzero(data, axis=(0, 1))
data = data[:,:,dx != 0]
# Some sequences may have less than 500 frames in them. To find out how many images there are, we will
# temporarily change chunkSize to be the same as whatever number of frames the last chunk of data has:
if data.shape[2] < chunkSize:
chunkSize_new = data.shape[2]
print("Number of images are less than chunkSize. chunkSize is temporarily changed to {}."
.format(chunkSize_new))
images = images + chunkSize_new
counter2 += 1
else:
images = counter1*chunkSize + counter2*chunkSize_new
counter1 += 1
data -= offsetMap.data # Offset correction
offset_corr_data = copy.copy(data) # I am copying this so that I can have access to it in the table below
histCalCorrected.fill(data)
cellTable=np.zeros(data.shape[2], np.int32) # Common mode correction
data = cmCorrection.correct(data.astype(np.float32), cellTable=cellTable)
histCalCMCorrected.fill(data)
noiseCalCM.fill(data) # Filling calculators with data
chunkSize = 100 # resetting the chunkSize to its default value for the next sequence or data-chunk
print('A total number of {} images are processed.'.format(images))
print("Offset and common mode corrections are applied.")
```
%% Cell type:code id: tags:
``` python
noiseMapCM = noiseCalCM.get() # Producing common mode corrected noise map
ho, eo, co, so = histCalCorrected.get()
hCM, eCM, cCM, sCM = histCalCMCorrected.get()
```
%% Cell type:code id: tags:
``` python
# We are copying these so that we can replot them later after the calculators are reset:
ho_second_trial = copy.copy(ho)
co_second_trial = copy.copy(co)
hCM_second_trial = copy.copy(hCM)
cCM_second_trial = copy.copy(cCM)
```
%% Cell type:markdown id: tags:
### Signal after Offset and Common Mode Corrections
Here, the offset corrected signal is compared to the common-mode corrected signal (in the form of binned histograms):
%% Cell type:code id: tags:
``` python
do = [{'x': co,
'y': ho,
'y_err': np.sqrt(ho[:]),
'drawstyle': 'steps-mid',
'color': 'cornflowerblue',
'label': 'Offset Corrected Signal'
},
{'x': cCM,
'y': hCM,
'y_err': np.sqrt(hCM[:]),
'drawstyle': 'steps-mid',
'color': 'red',
'label': 'Common Mode Corrected Signal'
}]
fig = xana.simplePlot(do, figsize='2col', aspect=1, x_label = 'Corrected Signal (ADU)', y_label="Counts",
x_range = (-200, 200), legend='top-right-frame-1col',
title = 'Corrected Signal - 2nd Iteration')
t0 = PrettyTable()
t0.title = "Comparison of the First Round of Corrections - Bad Pixels Not Excluded"
t0.field_names = ["After Offset Correction", "After Common Mode Correction"]
t0.add_row(["Mean: {:0.3f} (ADU)".format(np.mean(offset_corr_data)), "Mean: {:0.3f} (ADU)".format(np.mean(data))])
t0.add_row(["Median: {:0.3f} (ADU)".format(np.median(offset_corr_data)), "Median: {:0.3f} (ADU)"
.format(np.median(data))])
t0.add_row(["Standard Deviation: {:0.3f} (ADU)".format(np.std(offset_corr_data)),
"Standard Deviation: {:0.3f} (ADU)"
.format(np.std(data))])
print(t0,'\n')
```
%% Cell type:markdown id: tags:
### Noise Map after Common Mode Correction
In the following, the effect of common mode correction on the noise is shown. Finally common mode corrected noise map (noiseMapCM) is displayed and compared to the initial uncorrected noise map:
%% Cell type:code id: tags:
``` python
#***** NOISE MAP HISTOGRAM FROM THE COMMON MODE CORRECTED DATA *******#
hn, cn = np.histogram(noiseMap.flatten(), bins=2000, range=(0, 600))
hn_CM, cn_CM = np.histogram(noiseMapCM.flatten(), bins=2000, range=(0, 600))
dn = [{'x': cn[:-1],
'y': hn,
# 'y_err': np.sqrt(hn[:]),
'drawstyle': 'steps-mid', # 'bars',
'color': 'blue', # 'cornflowerblue',
'label': 'Uncorrected Noise (ADU)'
},
{'x': cn_CM[:-1],
'y': hn_CM,
# 'y_err': np.sqrt(hn_CM[:]),
'drawstyle': 'steps-mid', # 'bars',
'color': 'crimson', # 'red',#'cornflowerblue',
# 'ecolor': 'crimson',
'label': 'Common Mode Corrected Noise (ADU)'
}]
fig = xana.simplePlot(dn, figsize='1col', aspect=1, x_label='Noise (ADU)', y_label="Counts",
x_range=(0, 600), y_log=True, legend='top-center-frame-1col',
title='Noise Comparison')
fig = xana.heatmapPlot(noiseMapCM[:, :, 0],
x_label='Columns', y_label='Rows',
lut_label='Common Mode Corrected Noise (ADU)',
x_range=(0, y),
y_range=(0, x), vmax=2*np.mean(noiseMapCM), title='Common Mode Corrected Noise Map')
```
%% Cell type:code id: tags:
``` python
# Resetting the calculators:
noiseCalCM.reset() # resetting noise calculator
histCalCorrected.reset() # resetting histogram calculator
histCalCMCorrected.reset() # resetting histogram calculator
cmCorrection.reset()
```
%% Cell type:markdown id: tags:
### Final BadPixelMap
This is generated based on the offset and common mode corrected noise maps:
%% Cell type:code id: tags:
``` python
bad_pixels = np.zeros(offsetMap.shape, np.uint32)
mnoffset = np.nanmedian(offsetMap)
stdoffset = np.nanstd(offsetMap)
bad_pixels[(offsetMap < mnoffset-bad_pixel_offset_sigma*stdoffset) |
(offsetMap > mnoffset+bad_pixel_offset_sigma*stdoffset)] = BadPixels.OFFSET_OUT_OF_THRESHOLD.value
mnnoise = np.nanmedian(noiseMapCM)
stdnoise = np.nanstd(noiseMapCM)
bad_pixels[(noiseMapCM < mnnoise-bad_pixel_noise_sigma*stdnoise) |
(noiseMapCM > mnnoise+bad_pixel_noise_sigma*stdnoise)] = BadPixels.NOISE_OUT_OF_THRESHOLD.value
fig = xana.heatmapPlot(np.log2(bad_pixels[:, :, 0]),
x_label='Columns', y_label='Rows',
lut_label='Bad Pixel Value (ADU)',
x_range=(0, y),
y_range=(0, x), vmin=0, vmax=32)
```
%% Cell type:markdown id: tags:
### Calibration Constants
%% Cell type:code id: tags:
``` python
constant_maps = {
'Offset': offsetMap,
'Noise': noiseMapCM,
'BadPixelsDark': bad_pixels
}
for const_name in constant_maps.keys():
metadata = ConstantMetaData()
det = Constants.CCD(DetectorTypes.pnCCD)
const = getattr(det, const_name)()
const.data = constant_maps[const_name].data
metadata.calibration_constant = const
# set the operating condition
condition = Conditions.Dark.CCD(bias_voltage=bias_voltage,
integration_time=integration_time,
gain_setting=gain,
temperature=temperature_k,
pixels_x=1024,
pixels_y=1024)
device = Detectors.PnCCD1
metadata.detector_condition = condition
# specify the a version for this constant
if creation_time is None:
metadata.calibration_constant_version = Versions.Now(device=device)
else:
metadata.calibration_constant_version = Versions.Timespan(device=device,
start=creation_time)
metadata.calibration_constant_version.raw_data_location = file_loc
if db_output:
try:
metadata.send(cal_db_interface)
metadata.send(cal_db_interface, timeout=cal_db_timeout)
print("Inject {} constants from {}".format(const_name,
metadata.calibration_constant_version.begin_at))
except Exception as e:
print(e)
if local_output:
save_const_to_h5(metadata, out_folder)
print("Calibration constant {} is stored to {}.".format(const, out_folder))
```
Loading