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

changing fix_temp from bool to float

parent c1cef4ae
No related branches found
No related tags found
1 merge request!143fix/DARK-fix-temperature
%% Cell type:markdown id: tags:
# FastCCD Dark Characterization
Author: I. Klačková, S. Hauf, K. Setoodehnia and M. Cascella
The following notebook provides dark image analysis of the FastCCD detector.
Dark characterization evaluates offset and noise of the FastCCD detector, corrects the noise for Common Mode (CM), and defines bad pixels relative to offset and CM corrected noise. Bad pixels are then excluded and CM corrected noise is recalculated excluding the bad pixels. Resulting offset and CM corrected noise maps, as well as the bad pixel map are sent to the calibration database.
%% Cell type:code id: tags:
``` python
# Initial Parameters:
in_folder = "/gpfs/exfel/exp/SCS/201930/p900074/raw" # input folder, required
out_folder = '/gpfs/exfel/data/scratch/setoodeh/DarkRuns' # output folder, required
path_template = 'RAW-R{:04d}-DA05-S{{:05d}}.h5' # the template to use to access data
run = 351 # 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
# The two operation modes for FastCCD have fixed names which cannot be changed:
operation_mode = "FF" # FS stands for frame-store and FF for full-frame opeartion.
h5path = '/INSTRUMENT/SCS_CDIDET_FCCD2M/DAQ/FCCD:daqOutput/data/image/pixels' # path to the data in the HDF5 file
h5path_t = '/CONTROL/SCS_CDIDET_FCCD2M/CTRL/LSLAN/inputA/crdg/value' # path to find temperature
h5path_cntrl = '/RUN/SCS_CDIDET_FCCD2M/DET/FCCD' # path to find control data
cal_db_interface = "tcp://max-exfl016:8020" # the calibration database interface to use
cal_db_timeout = 300000 # timeout on calibration database requests
temp_limits = 5 # to find calibration constants later on, the sensor temperature is allowed to vary by 5 units
sequence = 0 # sequallence file to use
use_dir_creation_date = True # To be used to retrieve calibration constants later on (for database time derivation)
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
sigmaNoise = 5. # Any pixel whose signal exceeds 'sigmaNoise'*noiseCM (common mode corrected noise) will be masked
fix_temperature = False # This is required when the temperature is fixed in calibration database
fix_temperature = 0. # This is required when the temperature is fixed in calibration database
temperature_k = 233 # This is only used in the case of a fixed temperature for the calibration database
chunkSize = 100 # Number of images to read per chunk
cpuCores = 40 # Specifies the number of running cpu cores
commonModeAxis = 1 # Axis along which common mode will be calculated (0: along rows, 1: along columns)
ADU_to_electron_upper = 6.1 # According to Table 6.1 of Ivana Klačková's master's thesis, for upper hemisphere: conversion
# gain is 1 ADU = 6.1e-
ADU_to_electron_lower = 6.2 # and for lower hemisphere: conversion gain is 1 ADU = 6.2e-
run_parallel = True # For parallel computation
db_output = True # Output constants to the calibration database
```
%% Cell type:code id: tags:
``` python
# Required Packages:
import copy
import datetime
import os
import time
import warnings
warnings.filterwarnings('ignore')
import h5py
from IPython.display import display, Markdown
import matplotlib.pyplot as plt
%matplotlib inline
import numpy as np
from prettytable import PrettyTable
from iCalibrationDB import ConstantMetaData, Constants, Conditions, Detectors, Versions
from iCalibrationDB.detectors import DetectorTypes
from cal_tools.tools import get_dir_creation_date
from cal_tools.enums import BadPixels
from XFELDetAna import xfelpyanatools as xana
from XFELDetAna import xfelpycaltools as xcal
from XFELDetAna.detectors.fastccd import readerh5 as fastccdreaderh5
from XFELDetAna.util import env
env.iprofile = cluster_profile
import XFELDetAna.xfelprofiler as xprof
profiler = xprof.Profiler()
profiler.disable()
from XFELDetAna.xfelreaders import ChunkReader
```
%% Cell type:code id: tags:
``` python
# Output Folder Creation:
if not os.path.exists(out_folder):
os.makedirs(out_folder)
# Number of Images:
def nImagesOrLimit(nImages, limit):
if limit == 0:
return nImages
else:
return min(nImages, limit)
```
%% Cell type:code id: tags:
``` python
# Detector Operation Mode, Calibration Database Settings, and Some Initial Run Parameters & Paths:
display(Markdown('### Initial Settings'))
if operation_mode == "FS":
x = 960 # rows of the FastCCD to analyze in FS mode
y = 960 # columns of the FastCCD to analyze in FS mode
print('\nYou are analyzing data in FS mode.')
else:
x = 1934 # rows of the FastCCD to analyze in FF mode
y = 960 # columns of the FastCCD to analyze in FF mode
print('\nYou are analyzing data in FF mode.')
ped_dir = "{}/r{:04d}".format(in_folder, run)
fp_name = path_template.format(run)
fp_path = '{}/{}'.format(ped_dir, fp_name)
filename = fp_path.format(sequence)
creation_time = None
if use_dir_creation_date:
creation_time = get_dir_creation_date(in_folder, run)
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("Run number: {}".format(run))
print("Reading data from: {}".format(filename))
if creation_time:
print("Using {} as 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 # FastCCD has 1 memory cell
sensorSize = [x, y]
blockSize = [sensorSize[0]//2, sensorSize[1]] # 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
gain_setting = None
with h5py.File(filename, 'r') as f:
bias_voltage = int(f['{}/biasclock/bias/value'.format(h5path_cntrl)][0])
det_gain = int(f['{}/exposure/gain/value'.format(h5path_cntrl)][0])
integration_time = int(f['{}/acquisitionTime/value'.format(h5path_cntrl)][0])
temperature = np.mean(f[h5path_t])
```
%% Cell type:code id: tags:
``` python
# Printing the Parameters Read from the Data File:
display(Markdown('### Evaluated Parameters'))
print("Number of dark images to analyze:",nImages)
if det_gain == 8:
gain_setting = "high"
elif det_gain == 2:
gain_setting = "medium"
elif det_gain == 1:
gain_setting = "low"
else:
gain_setting = "auto"
print("Bias voltage is {} V".format(bias_voltage))
print("Detector gain is set to x{}".format(det_gain), "({} gain)".format(gain_setting))
print("Detector integration time is set to {}".format(integration_time), 'ms')
if fix_temperature:
temperature_k = fix_temperature + 273.15
if fix_temperature != 0.:
print("Using a fixed temperature of {}".format(temperature_k), "K")
print("Mean temperature was {:0.2f} °C / {:0.2f} K".format(temperature, temperature + 273.15))
else:
# This is needed while sending the
# calibration constant to the DB later
temperature_k = temperature + 273.15
print("Temperature is not fixed.")
print("Mean temperature was {:0.2f} °C / {:0.2f} K".format(temperature, temperature_k))
print("Mean temperature was {:0.2f} °C / {:0.2f} K".format(temperature, temperature + 273.15))
print("Output: {}".format(out_folder))
```
%% 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
# Calculator:
# 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 raw dark data.
%% Cell type:code id: tags:
``` python
for data in reader.readChunks():
data = np.bitwise_and(data.astype(np.uint16), 0b0011111111111111).astype(np.float32)
dx = np.count_nonzero(data, axis=(0, 1))
data = data[:,:,dx != 0]
noiseCal.fill(data) # Filling calculators with data
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 FastCCD offset, FastCCD 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",
x_range = (3400,4000), title = 'Offset Histogram')
#fig.savefig('Offset_Hist.svg', format='svg', dpi=1200, bbox_inches='tight')
t0 = PrettyTable()
t0.title = "Raw Signal"
t0.field_names = ["Mean","Median", "Standard Deviation"]
t0.add_row(["{:0.3f} (ADU)".format(np.mean(data)), "{:0.3f} (ADU)".format(np.median(data)), "{:0.3f} (ADU)".format(np.std(data))])
print(t0,'\n')
#************** OffsetMAP *******************#
fig = xana.heatmapPlot(offsetMap[:,:,0], x_label='Column Number', y_label='Row Number', aspect=1,
x_range=(0,y), y_range=(0,x), vmin=3000, vmax=4300, lut_label='Offset (ADU)',
panel_x_label='Columns Stat (ADU)', panel_y_label='Rows Stat (ADU)',
panel_top_low_lim = 3000, panel_top_high_lim = 4500, panel_side_low_lim = 3000,
panel_side_high_lim = 5000, title = 'OffsetMap')
#fig.savefig('RawOffsetMap.pdf', format='pdf', dpi=400, bbox_inches='tight')
#************** Raw NoiseMAP *******************#
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)', panel_top_low_lim = 0, panel_top_high_lim = 20,
panel_side_low_lim = 0, panel_side_high_lim = 50, title = 'Uncorrected NoiseMap')
#fig.savefig('RawNoiseMap.pdf', format='pdf', dpi=400, bbox_inches='tight')
```
%% Cell type:code id: tags:
``` python
# Offset Correction:
offsetCorrection = xcal.OffsetCorrection(sensorSize, offsetMap, nCells = memoryCells, cores=cpuCores, gains=None,
runParallel=run_parallel, blockSize=blockSize)
offsetCorrection.debug()
# Common Mode Correction:
# This is the new method subtracting the median of all pixels that are read out at the same time along a row:
cmCorrection = xcal.CommonModeCorrection([data.shape[0], data.shape[1]], [data.shape[0]//2, data.shape[1]],
commonModeAxis, parallel=False, dType=np.float32, stride=10,
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=600, range=[-200, 200], memoryCells=memoryCells,
cores=cpuCores, gains=None, blockSize=blockSize)
# For common mode corrected data:
histCalCMCorrected = xcal.HistogramCalculator(sensorSize, bins=600, range=[-200, 200], memoryCells=memoryCells,
cores=cpuCores, gains=None, blockSize=blockSize)
```
%% Cell type:markdown id: tags:
### Second Iteration
During the second iteration, the data are offset corrected and then common mode corrected to produced a common mode corrected noise map. The common mode correction is 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
for data in reader.readChunks():
data = data.astype(np.float32)
dx = np.count_nonzero(data, axis=(0, 1))
data = data[:,:,dx != 0]
data = offsetCorrection.correct(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) # Common mode correction
histCalCMCorrected.fill(data)
noiseCal.fill(data) # Filling noise calculator with common mode (CM) corrected data
noiseMapCM = noiseCal.get() # Produces CM corrected noise map
ho, eo, co , so = histCalCorrected.get()
hCM, eCM, cCM ,sCM = histCalCMCorrected.get()
print("Offset and common mode corrections are applied.")
```
%% Cell type:code id: tags:
``` python
# I am copying these so that I 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 = (-20,20), legend='top-right-frame-1col', title = 'Corrected Signal - 2nd Iteration')
#fig.savefig('Corrected_Signal_Hist_1.svg', format='svg', dpi=1200, bbox_inches='tight')
t0 = PrettyTable()
t0.title = "Comparison of the First Round of Corrections - Bad Pixels Included"
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=200, range=(2,40)) # hn: histogram of noise, cn: bin centers for noise
hn_CM,cn_CM = np.histogram(noiseMapCM.flatten(), bins=200, range=(2,40))
dn = [{'x': cn[:-1],
'y': hn,
#'y_err': np.sqrt(hn[:]),
'drawstyle': 'steps-mid',#'bars',
'color': 'blue',#'cornflowerblue',
'label': 'Uncorrected Noise'
},
{'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'
}]
fig = xana.simplePlot(dn, figsize='2col', aspect=1, x_label = 'Noise (ADU)', y_label="Counts",
x_range=(0,40), y_range=(0,1e6), y_log=True, legend='top-center-frame-1col',
title = 'Noise Comparison')
#fig.savefig('Noise_CM_1_Hist.svg', format='svg', dpi=1200, bbox_inches='tight')
fig = xana.heatmapPlot(noiseMapCM[:,:,0], aspect=1, x_label='Column Number', y_label='Row Number',
lut_label='Common Mode Corrected Noise (ADU)', x_range=(0,y), y_range=(0,x),
vmax=2*np.mean(noiseMapCM), panel_top_low_lim = 0, panel_top_high_lim = 20, panel_side_low_lim = 0,
panel_side_high_lim = 50, title = 'Common Mode Corrected Noise',
panel_x_label='Columns Stat (ADU)', panel_y_label='Rows Stat (ADU)')
#fig.savefig('NoiseMapCM.pdf', format='pdf', dpi=400, bbox_inches='tight')
```
%% Cell type:code id: tags:
``` python
# Resetting the calculators so we can do a third iteration later:
noiseCal.reset()
histCalCorrected.reset()
histCalCMCorrected.reset()
```
%% Cell type:markdown id: tags:
### Initial BadPixelMap
This is generated based on the offset and CM 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]),aspect=1, x_label='Column Number', y_label='Row Number',
lut_label='2^(Assigned Value to Bad Pixels)', x_range=(0,y), y_range=(0,x),
title = 'Bad Pixels Map Excluding Non-Sensitive Areas', panel_x_label= 'Columns Stat',
panel_y_label='Rows Stat')
```
%% Cell type:markdown id: tags:
Here, we are adding the pixels in the hole (center of the FastCCD) as well as 4 rows in the center of the detector, which we call overscan region:
%% Cell type:code id: tags:
``` python
def create_circular_mask(h, w, center=None, radius=None):
import numpy as np
import math
if center is None: # use the middle of the image
center = [int(w/2), int(h/2)]
if radius is None: # use the smallest distance between the center and image walls
radius = min(center[0], center[1], w-center[0], h-center[1])
Y, X = np.ogrid[:h, :w]
dist_from_center = np.sqrt((X - center[0])**2 + (Y-center[1])**2)
mask = dist_from_center < radius
return mask
```
%% Cell type:code id: tags:
``` python
mask = np.zeros(offsetMap.shape, np.uint32)
# Defining a circular mask + a rectangular mask (overscan) for the hole in the middle of the CCD:
h, w = (x,y)
hole_mask_bool = create_circular_mask(h-4, w, radius=61.5, center=(w//2,(h-4)//2))
hole_mask = np.zeros(hole_mask_bool.shape, np.uint32)
hole_mask[hole_mask_bool] = BadPixels.NON_SENSITIVE.value
overscan_mask = np.full((4, w), BadPixels.OVERSCAN.value)
mask[:,:,0] = np.insert(hole_mask, (h-4)//2, overscan_mask, 0)
# Assigning this masked area as bad pixels:
bad_pixels = np.bitwise_or(bad_pixels, mask)
fig = xana.heatmapPlot(np.log2(bad_pixels[:,:,0]),aspect=1, x_label='Column Number', y_label='Row Number',
lut_label='2^(Assigned Value to Bad Pixels)', x_range=(0,y), y_range=(0,x), panel_top_low_lim = 0,
panel_top_high_lim = 20, panel_side_low_lim = 0, panel_side_high_lim = 20,
title = 'Bad Pixels Map Including Non-Sensitive Areas', panel_x_label='Columns Stat',
panel_y_label='Rows Stat', vmax=20)
#fig.savefig('BadPixelMap_1.svg', format='svg', dpi=1200, bbox_inches='tight')
```
%% Cell type:markdown id: tags:
### Third Iteration
During the third iteration, the bad pixel map is applied to the data. Bad pixels are masked. Offset and common mode corrections are applied once again to the data, which now have bad pixdels excluded, to produce a common mode corrected noise map:
%% Cell type:code id: tags:
``` python
# bad_pixels is an array of (1934, 960, 1) filled with zeros except at indices where we have the actual bad pixels, whose
# values are set to be: 2 (2^1: BadPixels.OFFSET_OUT_OF_THRESHOLD.value), or
# 262144 (2^18: BadPixels.OVERSCAN.value), or 524288 (2^19: BadPixels.NON_SENSITIVE.value). These indices can be found
# using np.argwhere(bad_pixels != 0)
event_threshold = sigmaNoise*np.median(noiseMapCM) # for exclusion of possible cosmic ray events
noiseCal.setBadPixelMask(bad_pixels != 0)
```
%% 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]
data_copy = offsetCorrection.correct(copy.copy(data))
cellTable=np.zeros(data_copy.shape[2], np.int32)
data_copy = cmCorrection.correct(data_copy.astype(np.float32), cellTable=cellTable)
data[data_copy > event_threshold] = np.nan # cosmic rays
data = np.ma.MaskedArray(data, np.isnan(data), fill_value=0) # masking cosmics, the default fill_value is 1e+20
data = offsetCorrection.correct(data)
offset_corr_data2 = 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)
data = cmCorrection.correct(data.astype(np.float32), cellTable=cellTable)
histCalCMCorrected.fill(data)
noiseCal.fill(data)
noiseMapCM_2nd = noiseCal.get().filled(0) # the masked pixels are filled with zero
ho2, eo2, co2, so2 = histCalCorrected.get()
hCM2, eCM2, cCM2 ,sCM2 = histCalCMCorrected.get()
print("Final iteration is Performed.")
```
%% Cell type:markdown id: tags:
### Plots of the Final Results
The following plot and table compare the offset and common mode corrected signal with and without the bad pixels:
%% Cell type:code id: tags:
``` python
do_Final = [{'x': co_second_trial,
'y': ho_second_trial,
'y_err': np.sqrt(ho_second_trial[:]),
'drawstyle': 'steps-mid',
'color': 'blue',#'cornflowerblue',
'errorstyle': 'bars',
'label': 'Offset Corrected Signal, Bad Pixels Included - 2nd Trial'
},
{'x': cCM_second_trial,
'y': hCM_second_trial,
'y_err': np.sqrt(hCM_second_trial[:]),
'drawstyle': 'steps-mid',
'color': 'red',
'errorstyle': 'bars',
'ecolor': 'crimson',
'label': 'Common Mode Corrected Signal, Bad Pixels Included - 2nd Trial'
},
{'x': co2,
'y': ho2,
'y_err': np.sqrt(ho2[:]),
'drawstyle': 'steps-mid',
'color': 'black', #'cornflowerblue',
'errorstyle': 'bars',
'label': 'Offset Corrected Signal, Bad Pixels Excluded - 3rd Trial'
},
{'x': cCM2,
'y': hCM2,
'y_err': np.sqrt(hCM2[:]),
'drawstyle': 'steps-mid',
'color': 'orange', #'cornflowerblue',
'errorstyle': 'bars',
'label': 'Common Mode Corrected Signal, Bad Pixels Excluded - 3rd Trial'
}]
fig = xana.simplePlot(do_Final, figsize='2col', aspect=1, x_label = 'Corrected Signal (ADU)',
y_label="Counts (Logarithmic Scale)", y_log=True, x_range=(-40,40), legend='bottom-left-frame-1col',
title = 'Comparison of Corrected Signal')
#fig.savefig('Corrected_Signal_Hist_2.svg', format='svg', dpi=1200, bbox_inches='tight')
# offset_corr_data2 and data most likely have some nan's => I am going to use nanmean, nanmedian and nanstd functions:
t0 = PrettyTable()
t0.title = "Comparison of the Second Round of Corrections - Bad Pixels Excluded"
t0.field_names = ["After Offset Correction","After Common Mode Correction"]
t0.add_row(["Mean: {:0.3f} (ADU)".format(np.nanmean(offset_corr_data2)), "Mean: {:0.3f} (ADU)".format(np.nanmean(data))])
t0.add_row(["Median: {:0.3f} (ADU)".format(np.nanmedian(offset_corr_data2)), "Median: {:0.3f} (ADU)".format(np.nanmedian(data))])
t0.add_row(["Standard Deviation: {:0.3f} (ADU)".format(np.nanstd(offset_corr_data2)), "Standard Deviation: {:0.3f} (ADU)".format(np.nanstd(data))])
print(t0,'\n')
```
%% Cell type:markdown id: tags:
### Final NoiseMap
The effect of exclusion of bad pixels on common mode corrected noise is shown below. Finally common mode corrected noise map with bad pixels excluded (noiseMapCM_2nd) is displayed:
%% Cell type:code id: tags:
``` python
#*****NOISE MAP HISTOGRAM FROM THE COMMON MODE CORRECTED DATA*******#
hn_CM2,cn_CM2 = np.histogram(noiseMapCM_2nd.flatten(), bins=200, range=(2,40))
dn2 = [{'x': cn[:-1],
'y': hn,
#'y_err': np.sqrt(hn[:]),
'drawstyle': 'steps-mid',#'bars',
'color': 'blue', #'cornflowerblue',
'label': 'Uncorrected Noise'
},
{'x': cn_CM[:-1],
'y': hn_CM,
#'y_err': np.sqrt(hn_CM[:]),
'drawstyle': 'steps-mid',
'color': 'red',
#'ecolor': 'crimson',
'label': 'Common Mode Corrected Noise prior to Bad Pixels Exclusion'
},
{'x': cn_CM2[:-1],
'y': hn_CM2,
#'y_err': np.sqrt(hn_CM2[:]),
'drawstyle': 'steps-mid',
'color': 'black', #'cornflowerblue',
'label': 'Common Mode Corrected Noise after Bad Pixels Exclusion'
}]
fig = xana.simplePlot(dn2, figsize='2col', aspect = 1, x_label = 'Noise (ADU)', y_label="Counts", y_log=True,
x_range=(0,40), y_range=(0,1e6), legend='top-right-frame-1col', title = 'Final Noise Comparison')
#fig.savefig('Noise_Hist_2.svg', format='svg', dpi=1200, bbox_inches='tight')
fig = xana.heatmapPlot(np.log2(noiseMapCM_2nd[:,:,0]), aspect=1, x_label='Column Number', y_label='Row Number',
lut_label='Noise (ADU)', x_range=(0,y), y_range=(0,x), vmax=2*np.mean(noiseMapCM_2nd),
title = 'Final Common Mode Corrected Noise (Bad Pixels Excluded)',
panel_x_label='Columns Stat (ADU)', panel_y_label='Rows Stat (ADU)')
#fig.savefig('NoiseMapCM_2nd.pdf', format='pdf', dpi=400, bbox_inches='tight')
```
%% Cell type:markdown id: tags:
### Final Bad Pixel Map
Lastly, the final bad pixel map is generated based on the OffsetMap and the noiseMapCM_2nd (common mode corrected noise after exclusion of the initial bad pixels):
%% 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_2nd)
stdnoise = np.nanstd(noiseMapCM_2nd)
bad_pixels[(noiseMapCM_2nd < mnnoise-bad_pixel_noise_sigma*stdnoise) |
(noiseMapCM_2nd > mnnoise+bad_pixel_noise_sigma*stdnoise)] = BadPixels.NOISE_OUT_OF_THRESHOLD.value
bad_pixels = np.bitwise_or(bad_pixels, mask)
fig = xana.heatmapPlot(np.log2(bad_pixels[:,:,0]),aspect=1, x_label='Column Number', y_label='Row Number',
lut_label='2^(Assigned Value to Bad Pixels)', x_range=(0,y), y_range=(0,x), panel_top_low_lim = 0,
panel_top_high_lim = 20, panel_side_low_lim = 0, panel_side_high_lim = 20,
title = 'Final Bad Pixels Map', panel_x_label='Columns Stat',
panel_y_label='Rows Stat', vmax=20)
#fig.savefig('BadPixelMap_2.svg', format='svg', dpi=1200, bbox_inches='tight')
```
%% Cell type:code id: tags:
``` python
display(Markdown('### Statistics on the Bad Pixels'))
num_bad_pixels = np.count_nonzero(bad_pixels)
num_all_pixels = x*y
percentage_bad_pixels = num_bad_pixels*100/num_all_pixels
print("Number of bad pixels: {:0.0f}, i.e. {:0.2f}% of all pixels".format(num_bad_pixels, percentage_bad_pixels))
```
%% Cell type:markdown id: tags:
### Electronic Noise
According to Table 6.1 (page 80) of Ivana Klačková's master's thesis: "Conversion gain for the FastCCD is: lower hemisphere = 6.2e-/ADU and upper hemisphere = 6.1e-/ADU."
The following Tables present the noise along lower hemisphere, upper hemisphere, and the entire FastCCD detector at different stages. Here, the values in the first table (in ADU and e-) are the mean of noise per pixel, where noise is considered to be the initial uncorrected noise, CM corrected noise after second trial (including bad pixels) and CM corrected noise after third trial (excluding bad pixels).
The values of the second table (in electrons) are the standard deviation of noise per pixel.
%% Cell type:code id: tags:
``` python
# noiseMap refers to the initial uncorrected noise, noiseMapCM refers to common mode corrected noise with inclusion of
# bad pixels, and noiseMapCM_2nd refers to common mode corrected noise without inclusion of bad pixels:
ADU_to_electron = (ADU_to_electron_upper + ADU_to_electron_lower)/2 # Average of ADU_to_electron for the entire detector
print("Abbreviations:")
print(" - ED = Entire Detector; LH: Lower Hemisphere; UH: Upper Hemisphere")
print(" - CM Noise: Common Mode Corrected Noise")
print(" - BP: Bad Pixels\n")
t0 = PrettyTable()
t0.title = "Averages of Noise per Pixel"
t0.field_names = ["Uncorrected Noise","CM Noise, BP Incl.", "CM Noise, BP Excl."]
t0.add_row(["ED: {:0.2f} ADU = {:0.2f} e-".format(np.mean(noiseMap),np.mean(noiseMap)*ADU_to_electron), "ED: {:0.2f} ADU = {:0.2f} e-".format(np.mean(noiseMapCM), np.mean(noiseMapCM)*ADU_to_electron), "ED: {:0.2f} ADU = {:0.2f} e-".format(np.mean(noiseMapCM_2nd), np.mean(noiseMapCM_2nd)*ADU_to_electron)])
t0.add_row(["LH: {:0.2f} ADU = {:0.2f} e-".format(np.mean(noiseMap[:x//2,:]), np.mean(noiseMap[:x//2,:])*ADU_to_electron_lower), "LH: {:0.2f} ADU = {:0.2f} e-".format(np.mean(noiseMapCM[:x//2,:]), np.mean(noiseMapCM[:x//2,:])*ADU_to_electron_lower), "LH: {:0.2f} ADU = {:0.2f} e-".format(np.mean(noiseMapCM_2nd[:x//2,:]), np.mean(noiseMapCM_2nd[:x//2,:])*ADU_to_electron_lower)])
t0.add_row(["UH: {:0.2f} ADU = {:0.2f} e-".format(np.mean(noiseMap[x//2:,:]), np.mean(noiseMap[x//2:,:])*ADU_to_electron_upper), "UH: {:0.2f} ADU = {:0.2f} e-".format(np.mean(noiseMapCM[x//2:,:]), np.mean(noiseMapCM[x//2:,:])*ADU_to_electron_upper), "UH: {:0.2f} ADU = {:0.2f} e-".format(np.mean(noiseMapCM_2nd[x//2:,:]), np.mean(noiseMapCM_2nd[x//2:,:])*ADU_to_electron_upper)])
print(t0,'\n')
t1 = PrettyTable()
t1.title = "Standard Deviations of Noise per Pixel"
t1.field_names = ["Uncorrected Noise","CM Noise, BP Incl.", "CM Noise, BP Excl."]
t1.add_row(["ED: {:0.2f} e-".format(np.std(noiseMap)*ADU_to_electron), "ED: {:0.2f} e-".format(np.std(noiseMapCM)*ADU_to_electron), "ED: {:0.2f} e-".format(np.std(noiseMapCM_2nd)*ADU_to_electron)])
t1.add_row(["LH: {:0.2f} e-".format(np.std(noiseMap[:x//2,:])*ADU_to_electron_lower), "LH: {:0.2f} e-".format(np.std(noiseMapCM[:x//2,:])*ADU_to_electron_lower), "LH: {:0.2f} e-".format(np.std(noiseMapCM_2nd[:x//2,:])*ADU_to_electron_lower)])
t1.add_row(["UH: {:0.2f} e-".format(np.std(noiseMap[x//2:,:])*ADU_to_electron_upper), "UH: {:0.2f} e-".format(np.std(noiseMapCM[x//2:,:])*ADU_to_electron_upper), "UH: {:0.2f} e-".format(np.std(noiseMapCM_2nd[x//2:,:])*ADU_to_electron_upper)])
print(t1)
```
%% Cell type:markdown id: tags:
### Calibration Constants
%% Cell type:code id: tags:
``` python
dictionary = {}
dictionary['Offset'] = offsetMap.data
dictionary['Noise'] = noiseMapCM_2nd.data
dictionary['BadPixelsDark'] = bad_pixels.data
for const in dictionary:
metadata = ConstantMetaData()
dconst = getattr(Constants.CCD(DetectorTypes.fastCCD), const)()
dconst.data = dictionary[const]
metadata.calibration_constant = dconst
condition = Conditions.Dark.CCD(bias_voltage=bias_voltage,
integration_time=integration_time,
gain_setting=det_gain,
temperature=temperature_k,
pixels_x=1934,
pixels_y=960)
for parm in condition.parameters:
if parm.name == "Sensor Temperature":
parm.lower_deviation = temp_limits
parm.upper_deviation = temp_limits
device = Detectors.fastCCD1
metadata.detector_condition = condition
# Specifying 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)
if db_output:
metadata.send(cal_db_interface, timeout=cal_db_timeout)
print("Calibration constants (offsetMap, noiseMapCM_2nd and bad_pixels) are sent to the calibration database.")
print("Creation time is: {}".format(creation_time))
```
%% Cell type:code id: tags:
``` python
```
......
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