Skip to content
Snippets Groups Projects
Commit 086f33ec authored by Karim Ahmed's avatar Karim Ahmed
Browse files

change False local_output

parent 1acfbe6b
No related branches found
No related tags found
1 merge request!232db_output false by default
%% Cell type:markdown id: tags:
# Jungfrau Dark Image Characterization #
Version: 0.1, Author: M. Ramilli, S. Hauf
Analyzes Jungfrau dark image data to deduce offset, noise and resulting bad pixel maps
%% Cell type:code id: tags:
``` python
in_folder = '/gpfs/exfel/exp/FXE/201931/p900088/raw/' # folder under which runs are located, required
out_folder = '/gpfs/exfel/data/scratch/karnem/test_dark/' # path to place reports at, required
sequences = 1 # number of sequence files in that run
path_template = 'RAW-R{:04d}-{}-S{{:05d}}.h5' # template to use for file name, double escape sequence number
path_inset = "JNGFR02" # file inset for image data
path_inset_control = "JNGFR02" # file inset for control data
cluster_profile = 'noDB' # the ipcluster profile name
cal_db_interface = 'tcp://max-exfl016:8016' # calibrate db interface to connect to
integration_time = 1000 # integration time in us, will be overwritten by value in file
bias_voltage = 90 # sensor bias voltage in V, will be overwritten by value in file
badpixel_threshold_sigma = 20. # bad pixels defined by values outside n times this std from median
offset_abs_threshold_low = [1000, 10000, 10000] # absolute bad pixel threshold in terms of offset, lower values
offset_abs_threshold_high = [8000, 15000, 15000] # absolute bad pixel threshold in terms of offset, upper values
chunkSize = 10 # iteration chunk size, needs to match or be less than number of images in a sequence file
imageRange = [0, 500] # image range in which to evaluate
memoryCells = 16 # number of memory cells
h5path = '/INSTRUMENT/{}/DET/{}:daqOutput/data' # path in H5 file under which images are located
run_high = 144 # run number for G0 dark run, required
run_med = 145 # run number for G1 dark run, required
run_low = 146 # run number for G2 dark run, required
karabo_id = "FXE_XAD_JF500K" # karabo prefix of Jungfrau devices
karabo_id_control = "" # if control is on a different ID, set to empty string for using the same as for data
receiver_id = "RECEIVER" # inset for receiver devices
control_id = "CONTROL" # inset for control devices
db_module = "Jungfrau_M260" # ID of module in calibration database
use_dir_creation_date = True # use dir creation date
manual_slow_data = False # if true, use manually entered bias_voltage and integration_time values
time_limits = 0.025 # to find calibration constants later on, the integration time is allowed to vary by 0.5 us
local_output = False # output constants locally
local_output = True # output constants locally
db_output = False # output constants to database
```
%% Cell type:code id: tags:
``` python
import warnings
warnings.filterwarnings('ignore')
import matplotlib
matplotlib.use('agg')
import numpy as np
import h5py
from h5py import File as h5file
from cal_tools.enums import BadPixels
from cal_tools.tools import get_dir_creation_date
from cal_tools.ana_tools import save_dict_to_hdf5
from iCalibrationDB import ConstantMetaData, Constants, Conditions, Detectors, Versions
from XFELDetAna.util import env
env.iprofile = cluster_profile
from XFELDetAna.detectors.jungfrau import readerPSI as jfreaderPSI
from XFELDetAna.detectors.jungfrau import reader as jfreader
from XFELDetAna.detectors.jungfrau.jf_chunk_reader import JFChunkReader
from XFELDetAna.detectors.jungfrau.util import non_empty_trains, count_n_files, rollout_data, sanitize_data_cellid
import glob
import matplotlib.pyplot as plt
%matplotlib inline
from XFELDetAna.plotting.histogram import histPlot
from XFELDetAna.plotting.heatmap import heatmapPlot
# Constants relevant for the analysis
run_nums = [run_high, run_med, run_low] # run number for G0/HG0, G1, G2
sensorSize = [1024, 512]
blockSize = [1024, 512]
xRange = [0, 0+sensorSize[0]]
yRange = [0, 0+sensorSize[1]]
gains = [0, 1, 2]
h5path = h5path.format(karabo_id, receiver_id)
creation_time = None
if use_dir_creation_date:
creation_time = get_dir_creation_date(in_folder, run_high)
print("Using {} as creation time".format(creation_time))
offset_abs_threshold = [offset_abs_threshold_low, offset_abs_threshold_high]
if karabo_id_control == "":
karabo_id_control = karabo_id
import os
os.makedirs(out_folder, exist_ok=True)
```
%% Cell type:code id: tags:
``` python
proposal = list(filter(None, in_folder.strip('/').split('/')))[-2]
file_loc = 'proposal:{} runs:{} {} {}'.format(proposal, run_high, run_med, run_low)
```
%% Cell type:code id: tags:
``` python
def check_memoryCells(file_name, path):
with h5file(file_name, 'r') as f:
t_stamp = np.array(f[path + '/storageCells/timestamp'])
st_cells = np.array(f[path + '/storageCells/value'])
sc_start = np.array(f[path + '/storageCellStart/value'])
valid_train = t_stamp > 0
n_scs = st_cells[valid_train][0] + 1
sc_s = sc_start[valid_train][0]
return n_scs, sc_s
```
%% Cell type:code id: tags:
``` python
chunkSize = 100
filep_size = 1000
noiseCal = None
noise_map = None
offset_map = None
memoryCells = None
for i, r_n in enumerate(run_nums):
gain = i
print(f"Gain stage {gain}, run {r_n}")
valid_data = []
valid_cellids = []
if r_n is not None:
n_tr = 0
n_empty_trains = 0
n_empty_sc = 0
ped_dir = "{}/r{:04d}/".format(in_folder, r_n)
fp_name = path_template.format(r_n, path_inset_control)
fp_path = '{}/{}'.format(ped_dir, fp_name)
n_files = len(glob.glob("{}/*{}*.h5".format(ped_dir, path_inset)))
myRange = range(0, n_files)
control_path = '/CONTROL/{}/DET/{}'.format(karabo_id_control, control_id)
this_run_mcells, sc_start = check_memoryCells(fp_path.format(0).format(myRange[0]), control_path)
if noise_map is None:
if not manual_slow_data:
with h5py.File(fp_path.format(0), 'r') as f:
integration_time = float(f['/RUN/{}/DET/{}/exposureTime/value'.format(karabo_id_control, control_id)][()]*1e6)
bias_voltage = int(np.squeeze(f['/RUN/{}/DET/{}/vHighVoltage/value'.format(karabo_id_control, control_id)])[0])
print("Integration time is {} us".format(integration_time))
print("Bias voltage is {} V".format(bias_voltage))
if this_run_mcells == 1:
memoryCells = 1
print('Dark runs in single cell mode\n storage cell start: {:02d}'.format(sc_start))
else:
memoryCells = 16
print('Dark runs in burst mode\n storage cell start: {:02d}'.format(sc_start))
noise_map = np.zeros(sensorSize+[memoryCells, 3])
offset_map = np.zeros(sensorSize+[memoryCells, 3])
fp_name = path_template.format(r_n, path_inset)
fp_path = '{}/{}'.format(ped_dir, fp_name)
myRange_P = range(0, sequences)
path = h5path
print("Reading data from {}".format(fp_path))
print("Run is: {}".format(r_n))
print("HDF5 path: {}".format(h5path))
imageRange = [0, filep_size*len(myRange)]
reader = JFChunkReader(filename = fp_path, readFun = jfreader.readData, size = filep_size, chunkSize = chunkSize,
path = h5path, image_range=imageRange, pixels_x = sensorSize[0], pixels_y = sensorSize[1],
x_range = xRange, y_range = yRange, imagesPerChunk=chunkSize, filesRange = myRange,
memoryCells=this_run_mcells, blockSize=blockSize)
for data in reader.readChunks():
images = np.array(data[0], dtype=np.float)
gainmaps = np.array(data[1], dtype=np.uint16)
trainId = np.array(data[2])
fr_num = np.array(data[3])
acelltable = np.array(data[4])
n_tr += acelltable.shape[-1]
this_tr = acelltable.shape[-1]
idxs = non_empty_trains(trainId)
images = images[..., idxs]
gainmaps = gainmaps[..., idxs]
fr_num = fr_num[..., idxs]
acelltable = acelltable[..., idxs]
if memoryCells == 1:
acelltable -= sc_start
n_empty_trains += this_tr - acelltable.shape[-1]
n_empty_sc += len(acelltable[acelltable > 15])
if i > 0 and memoryCells == 16: ## throwing away the second memory cell entry in lower gains
acelltable[1] = 255
images, gainmaps, acelltable = rollout_data([images, gainmaps, acelltable]) # makes 4-dim vecs into 3-dim
# makes 2-dim into 1-dim
# leaves 1-dim and 3-dim vecs
images, gainmaps, acelltable = sanitize_data_cellid([images, gainmaps], acelltable) # removes entries with cellID 255
valid_data.append(images)
valid_cellids.append(acelltable)
valid_data = np.concatenate(valid_data, axis=2)
valid_cellids = np.concatenate(valid_cellids, axis=0)
for cell in range(memoryCells):
thiscell = valid_data[...,valid_cellids == cell]
noise_map[...,cell,gain] = np.std(thiscell, axis=2)
offset_map[...,cell,gain] = np.mean(thiscell, axis=2)
print('G{:01d} dark calibration'.format(i))
print('Missed {:d} out of {:d} trains'.format(n_empty_trains, n_tr))
print('Lost {:d} images out of {:d}'.format(n_empty_sc, this_run_mcells * (n_tr - n_empty_trains)))
else:
print('missing G{:01d}'.format(i))
```
%% Cell type:markdown id: tags:
## Offset and Noise Maps ##
Below offset and noise maps for the high ($g_0$) gain stage are shown, alongside the distribution of these values. One expects block-like structures mapping to the ASICs of the detector
%% Cell type:code id: tags:
``` python
import matplotlib.pyplot as plt
from XFELDetAna.core.util import remove_nans
%matplotlib inline
#%matplotlib notebook
from XFELDetAna.plotting.histogram import histPlot
from XFELDetAna.plotting.heatmap import heatmapPlot
g_name = ['G0', 'G1', 'G2']
g_range = [(0, 8000), (8000, 16000), (8000, 16000)]
n_range = [(0., 50.), (0., 50.), (0., 50.)]
unit = '[ADCu]'
```
%% Cell type:code id: tags:
``` python
for g_idx in gains:
for cell in range(0, memoryCells):
f_o0 = heatmapPlot(np.swapaxes(offset_map[..., cell, g_idx], 0, 1),
y_label="Row",
x_label="Column",
lut_label="Pedestal {} [ADCu]".format(g_name[g_idx]),
aspect=1.,
vmin=g_range[g_idx][0],
vmax=g_range[g_idx][1])
fo0, ax_o0 = plt.subplots()
res_o0 = histPlot(ax_o0, offset_map[..., cell, g_idx],
bins=800,
range=g_range[g_idx],
facecolor='b',
histotype='stepfilled')
ax_o0.tick_params(axis='both',which='major',labelsize=15)
ax_o0.set_title('Module pedestal distribution Cell {:02d}'.format(cell), fontsize=15)
ax_o0.set_xlabel('Pedestal {} [ADCu]'.format(g_name[g_idx]),fontsize=15)
ax_o0.set_yscale('log')
f_n0 = heatmapPlot(np.swapaxes(noise_map[..., cell, g_idx], 0, 1),
y_label="Row",
x_label="Column",
lut_label="RMS noise {} ".format(g_name[g_idx]) + unit,
aspect=1.,
vmin=n_range[g_idx][0],
vmax=n_range[g_idx][1])
fn0, ax_n0 = plt.subplots()
res_n0 = histPlot(ax_n0, noise_map[..., cell, g_idx],
bins=100,
range=n_range[g_idx],
facecolor='b',
histotype='stepfilled')
ax_n0.tick_params(axis='both',which='major',labelsize=15)
ax_n0.set_title('Module noise distribution Cell {:02d}'.format(cell), fontsize=15)
ax_n0.set_xlabel('RMS noise {} '.format(g_name[g_idx]) + unit,fontsize=15)
#ax_n0.set_yscale('log')
plt.show()
```
%% Cell type:markdown id: tags:
## Bad Pixel Map ###
The bad pixel map is deduced by comparing offset and noise of each pixel ($v_i$) and each gain ($g$) against the median value for that gain stage:
$$
v_i > \mathrm{median}(v_{k,g}) + n \sigma_{v_{k,g}}
$$
or
$$
v_i < \mathrm{median}(v_{k,g}) - n \sigma_{v_{k,g}}
$$
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):
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
bad_pixels_map = np.zeros(noise_map.shape, np.uint32)
def eval_bpidx(d):
mdn = np.nanmedian(d, axis=(0, 1, 2))[None, None, None, :]
std = np.nanstd(d, axis=(0, 1, 2))[None, None, None, :]
idx = (d > badpixel_threshold_sigma*std+mdn) | (d < (-badpixel_threshold_sigma)*std+mdn)
return idx
offset_abs_threshold = np.array(offset_abs_threshold)
bad_pixels_map[eval_bpidx(offset_map)] = BadPixels.OFFSET_OUT_OF_THRESHOLD.value
bad_pixels_map[~np.isfinite(offset_map)] |= BadPixels.OFFSET_NOISE_EVAL_ERROR.value
bad_pixels_map[eval_bpidx(noise_map)] |= BadPixels.NOISE_OUT_OF_THRESHOLD.value
bad_pixels_map[~np.isfinite(noise_map)] |= BadPixels.OFFSET_NOISE_EVAL_ERROR.value
bad_pixels_map[(offset_map < offset_abs_threshold[0][None, None, None, :]) | (offset_map > offset_abs_threshold[1][None, None, None, :])] |= BadPixels.OFFSET_OUT_OF_THRESHOLD.value
for g_idx in gains:
for cell in range(0, memoryCells):
bad_pixels = bad_pixels_map[:, :, cell, g_idx]
fn_0 = heatmapPlot(np.swapaxes(bad_pixels, 0, 1),
y_label="Row",
x_label="Column",
lut_label="Badpixels {} [ADCu]".format(g_name[g_idx]),
aspect=1.,
vmin=0)
```
%% Cell type:code id: tags:
``` python
constants = {'Offset': np.moveaxis(offset_map, 0, 1),
'Noise': np.moveaxis(noise_map, 0, 1),
'BadPixelsDark': np.moveaxis(bad_pixels_map, 0, 1)}
for key, const_data in constants.items():
metadata = ConstantMetaData()
const = getattr(Constants.jungfrau, key)()
const.data = const_data
metadata.calibration_constant = const
# set the operating condition
condition = Conditions.Dark.jungfrau(memory_cells=memoryCells, bias_voltage=bias_voltage,
integration_time=integration_time)
for parm in condition.parameters:
if parm.name == "Integration Time":
parm.lower_deviation = time_limits
parm.upper_deviation = time_limits
device = getattr(Detectors, db_module)
metadata.detector_condition = condition
# specify the 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:
metadata.send(cal_db_interface)
print('Constants {} is sent to the data base'.format(key))
# save everything to file.
# one file - one entrt in the DB
if local_output:
dpar = {}
for parm in metadata.detector_condition.parameters:
dpar[parm.name] = {'lower_deviation': parm.lower_deviation,
'upper_deviation': parm.upper_deviation,
'value': parm.value}
data_to_store = {}
data_to_store['condition'] = dpar
data_to_store['db_module'] = db_module
data_to_store['constant'] = key
data_to_store['data'] = const_data
data_to_store['creation_time'] = creation_time
data_to_store['dclass'] = 'jungfrau'
data_to_store['file_loc'] = file_loc
ofile = "{}/const_{}_{}.h5".format(out_folder, key, db_module)
save_dict_to_hdf5(data_to_store, ofile)
```
......
%% Cell type:markdown id: tags:
# pnCCD Dark Characterization
Author: DET Group, Version 1.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. Resulting maps are saved as .h5 files for a latter use.
%% Cell type:code id: tags:
``` python
in_folder = "/gpfs/exfel/exp/SQS/201921/p002430/raw/" # input folder, required
out_folder = 'gpfs/exfel/data/scratch/haufs/test/' # output folder, required
path_template = 'RAW-R{:04d}-PNCCD01-S{{:05d}}.h5' # the template to use to access data
run = 745 # which run to read data from, required
number_dark_frames = 0 # number of images to be used, if set to 0 all available images are used
cluster_profile = "noDB" # ipcluster profile to use
sigma_noise = 10. # Pixel exceeding 'sigmaNoise' * noise value in that pixel will be masked
h5path = '/INSTRUMENT/SQS_NQS_PNCCD1MP/CAL/PNCCD_FMT-0:output/data/image/' # path in the HDF5 file the data is at
cal_db_interface = "tcp://max-exfl016:8021" # calibration DB interface to use
temp_limits = 5 # temperature limits in which to consider calibration parameters equal
sequence = 0 # sequence file to use
multi_iteration = False # use multiple iterations
use_dir_creation_date = True # use dir creation date as data production reference date
bad_pixel_offset_sigma = 5. # deviation in standard deviations from mean offset to consider pixel bad
bad_pixel_noise_sigma = 5. # deviation in standard deviations from mean noise to consider pixel bad
fix_temperature = 233. # fix temperature to this value in K, set to -1 to use value from slow data
gain = 0 # the detector gain setting, only 0 is currently implemented
bias_voltage = 300 # detector bias voltage
integration_time = 70 # detector integration time
db_output = False # Output constants to the calibration database
local_output = False # output constants locally
local_output = True # output constants locally
```
%% Cell type:code id: tags:
``` python
from iCalibrationDB import ConstantMetaData, Constants, Conditions, Detectors, Versions
from iCalibrationDB.detectors import DetectorTypes
```
%% Cell type:code id: tags:
``` python
import XFELDetAna.xfelprofiler as xprof
profiler = xprof.Profiler()
profiler.disable()
from XFELDetAna.util import env
env.iprofile = cluster_profile
import warnings
warnings.filterwarnings('ignore')
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
import numpy as np
import h5py
import matplotlib.pyplot as plt
from iminuit import Minuit
import time
import copy
from prettytable import PrettyTable
%matplotlib inline
def nImagesOrLimit(nImages, limit):
if limit == 0:
return nImages
else:
return min(nImages, limit)
sigmaNoise = sigma_noise
temperature_k = fix_temperature
```
%% Cell type:code id: tags:
``` python
proposal = list(filter(None, in_folder.strip('/').split('/')))[-2]
file_loc = 'proposal:{} runs:{}'.format(proposal, run)
```
%% Cell type:code id: tags:
``` python
x = 1024 # rows of the FastCCD to analyze in FS mode
y = 1024 # columns of the FastCCD to analyze in FS mode
ped_dir = "{}/r{:04d}".format(in_folder, run)
fp_name = path_template.format(run)
import datetime
creation_time = None
if use_dir_creation_date:
creation_time = get_dir_creation_date(in_folder, run)
fp_path = '{}/{}'.format(ped_dir, fp_name)
print("Reading data from: {}\n".format(fp_path))
print("Run is: {}".format(run))
print("HDF5 path: {}".format(h5path))
if creation_time:
print("Using {} as creation time".format(creation_time.isoformat()))
```
%% Cell type:code id: tags:
``` python
filename = fp_path.format(sequence)
sensorSize = [x, y]
chunkSize = 100 #Number of images to read per chunk
#Sensor area will be analysed according to blocksize
blockSize = [sensorSize[0]//2, sensorSize[1]//4]
xcal.defaultBlockSize = blockSize
cpuCores = 8 #Specifies the number of running cpu cores
memoryCells = 1 #FastCCD has 1 memory cell
#Specifies total number of images to proceed
nImages = fastccdreaderh5.getDataSize(filename, h5path)[0]
nImages = nImagesOrLimit(nImages, number_dark_frames)
print("\nNumber of dark images to analyze: ",nImages)
commonModeBlockSize = blockSize
commonModeAxisR = 'row'#Axis along which common mode will be calculated
run_parallel = True
profile = False
```
%% Cell type:code id: tags:
``` python
reader = ChunkReader(filename, fastccdreaderh5.readData,
nImages, chunkSize,
path = h5path,
pixels_x = sensorSize[0],
pixels_y = sensorSize[1],)
```
%% Cell type:code id: tags:
``` python
noiseCal = xcal.NoiseCalculator(sensorSize, memoryCells,
cores=cpuCores, blockSize=blockSize,
runParallel=run_parallel)
histCalRaw = xcal.HistogramCalculator(sensorSize, bins=1000,
range=[0, 10000], parallel=False,
memoryCells=memoryCells,
cores=cpuCores, blockSize=blockSize)
```
%% Cell type:markdown id: tags:
### First Iteration
%% Cell type:markdown id: tags:
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 raw dark data.
%% Cell type:code id: tags:
``` python
for data in reader.readChunks():
data = data.astype(np.float32)
dx = np.count_nonzero(data, axis=(0, 1))
data = data[:,:,dx != 0]
histCalRaw.fill(data)
#Filling calculators with data
noiseCal.fill(data)
offsetMap = noiseCal.getOffset() #Produce offset map
noiseMap = noiseCal.get() #Produce noise map
noiseCal.reset() #Reset noise calculator
print("Initial maps were created")
```
%% Cell type:code id: tags:
``` python
#**************OFFSET MAP HISTOGRAM***********#
ho,co = np.histogram(offsetMap.flatten(), bins=700)
do = {'x': co[:-1],
'y': ho,
'y_err': np.sqrt(ho[:]),
'drawstyle': 'bars',
'color': 'cornflowerblue',
}
fig = xana.simplePlot(do, figsize='1col', aspect=2,
x_label = 'Offset (ADU)',
y_label="Counts", y_log=True,
)
#*****NOISE MAP HISTOGRAM FROM THE OFFSET CORRECTED DATA*******#
hn,cn = np.histogram(noiseMap.flatten(), bins=2000)
dn = {'x': cn[:-1],
'y': hn,
'y_err': np.sqrt(hn[:]),
'drawstyle': 'bars',
'color': 'cornflowerblue',
}
fig = xana.simplePlot(dn, figsize='1col', aspect=2,
x_label = 'Noise (ADU)',
y_label="Counts",
y_log=True,
x_range=(0, 400))
#**************HEAT MAPS*******************#
fig = xana.heatmapPlot(offsetMap[:,:,0],
x_label='Columns', y_label='Rows',
lut_label='Offset (ADU)',
x_range=(0,y),
y_range=(0,x), vmin=6000, vmax=14500)
fig = xana.heatmapPlot(noiseMap[:,:,0],
x_label='Columns', y_label='Rows',
lut_label='Noise (ADU)',
x_range=(0,y),
y_range=(0,x), vmax=2*np.mean(noiseMap))
```
%% Cell type:code id: tags:
``` python
cmCorrection = xcal.CommonModeCorrection(sensorSize,
[512, 512],
'row',
nCells = memoryCells,
dType=np.float32,
noiseMap = noiseMap.astype(np.float32),
runParallel=run_parallel,
stats=True)
cmCorrection.debug()
noiseCalCM = xcal.NoiseCalculator(sensorSize, memoryCells,
cores=cpuCores, blockSize=blockSize,
runParallel=run_parallel)
```
%% Cell type:code id: tags:
``` python
for data in reader.readChunks():
data = data.astype(np.float32)
dx = np.count_nonzero(data, axis=(0, 1))
data = data[:,:,dx != 0]
cellTable=np.zeros(data.shape[2], np.int32) # Common mode correction
data = cmCorrection.correct(data.astype(np.float32), cellTable=cellTable)
#Filling calculators with data
noiseCalCM.fill(data)
offsetMapCM = noiseCalCM.getOffset() #Produce offset map
noiseMapCM = noiseCalCM.get() #Produce noise map
noiseCalCM.reset() #Reset noise calculator
```
%% Cell type:code id: tags:
``` python
#*****NOISE MAP HISTOGRAM FROM THE OFFSET CORRECTED DATA*******#
hn,cn = np.histogram(noiseMapCM.flatten(), bins=2000)
dn = {'x': cn[:-1],
'y': hn,
'y_err': np.sqrt(hn[:]),
'drawstyle': 'bars',
'color': 'cornflowerblue',
}
fig = xana.simplePlot(dn, figsize='1col', aspect=2,
x_label = 'Noise (ADU)',
y_label="Counts",
y_log=True,
x_range=(0, 400))
fig = xana.heatmapPlot(noiseMapCM[:,:,0],
x_label='Columns', y_label='Rows',
lut_label='Noise (ADU)',
x_range=(0,y),
y_range=(0,x), vmax=2*np.mean(noiseMapCM))
```
%% Cell type:code id: tags:
``` python
from cal_tools.enums import BadPixels
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
constant_maps = {
'Offset': offsetMapCM,
'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)
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 locally.".format(const))
```
%% Cell type:code id: tags:
``` python
histCalCorr = xcal.HistogramCalculator(sensorSize, bins=200,
range=[-200, 200], parallel=False,
memoryCells=memoryCells,
cores=cpuCores, blockSize=blockSize)
for data in reader.readChunks():
data = data.astype(np.float32)
data -= offsetMap.data
histCalCorr.fill(data)
```
%% Cell type:code id: tags:
``` python
ho,eo,co,so = histCalCorr.get()
d = [{'x': co,
'y': ho,
'y_err': np.sqrt(ho[:]),
'drawstyle': 'steps-mid',
'errorstyle': 'bars',
'errorcoarsing': 2,
'label': 'Offset corr.'
},
]
fig = xana.simplePlot(d, aspect=1, x_label='Energy(ADU)',
y_label='Number of occurrences', figsize='2col',
y_log=True, x_range=(-50,500),
legend='top-center-frame-2col')
```
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment