pnCCD Gain Characterization

cluster_profile = "noDB" # ipcluster profile to use
in_folder = "/gpfs/exfel/exp/SQS/202031/p900166/raw" # input folder for the raw data
out_folder = '/gpfs/exfel/data/scratch/setoodeh/Test' # output folder
run = 347 # which run to read data from
sequences = [-1] # sequences to correct, set to -1 for all, range allowed

db_module = "pnCCD_M205_M206"
karabo_da = 'PNCCD01' # data aggregators
karabo_da_control = "PNCCD02" # file inset for control data
karabo_id = "SQS_NQS_PNCCD1MP" # karabo prefix of PNCCD devices
receiver_id = "PNCCD_FMT-0" # inset for receiver devices
path_template = 'CORR-R{:04d}-PNCCD01-S{{:05d}}.h5' # the template to use to access data
path_template_ctrl = 'RAW-R{:04d}-{}-S{{:05d}}.h5' # the template to use to access data
path_template_seqs = "{}/r{:04d}/*PNCCD01-S*.h5"
h5path = '/INSTRUMENT/{}/CAL/{}:output/data/' # path to data in the HDF5 file 
h5path_ctrl = '/CONTROL/{}/CTRL/TCTRL'

cpuCores = 40 # specifies the number of running cpu cores
use_dir_creation_date = True # this is needed to obtain creation time of the run
sequences_per_node = 1
chunkSize = 100 # number of images to read per chunk
run_parallel = True
db_output = False # if True, the notebook injects dark constants into the calibration database
local_output = True # if True, the notebook saves dark constants locally

cal_db_interface = "tcp://max-exfl016:8015" # calibration DB interface to use
cal_db_timeout = 300000 # timeout on caldb requests
creation_time = "" # To overwrite the measured creation_time. Required Format: YYYY-MM-DD HR:MN:SC.00 e.g. 2019-07-04 11:02:41.00

# pnCCD parameters:
fix_temperature_top = 0.  # fix temperature of top pnCCD sensor in K, set to 0. to use the value from slow data
fix_temperature_bot = 0.  # fix temperature of bottom pnCCD sensor in K, set to 0. to use the value from slow data
gain = 0.1  # the detector's gain setting. Set to 0, to use the value from slow data
bias_voltage = 0. # the detector's bias voltage. set to 0. to use the value from slow data.
integration_time = 70  # detector's integration time
photon_energy = 1.6 # Al fluorescence in keV

def balance_sequences(in_folder, run, sequences, sequences_per_node, karabo_da):
    from xfel_calibrate.calibrate import balance_sequences as bs
    return bs(in_folder, run, sequences, sequences_per_node, karabo_da)

import copy
import datetime
import glob
import os
import traceback
import warnings
from datetime import timedelta

warnings.filterwarnings('ignore')

from functools import partial

import h5py
import iminuit as im
import matplotlib
from iminuit import Minuit
from IPython.display import Markdown, display

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import XFELDetAna.xfelprofiler as xprof
from cal_tools.pnccdlib import VALID_GAINS, extract_slow_data
from cal_tools.tools import (
    get_dir_creation_date,
    get_pdu_from_db,
    get_report,
    save_const_to_h5,
    send_to_db,
)
from iCalibrationDB import Conditions, ConstantMetaData, Constants, Detectors, Versions
from iCalibrationDB.detectors import DetectorTypes
from mpl_toolkits.axes_grid1 import AxesGrid, ImageGrid
from prettytable import PrettyTable

profiler = xprof.Profiler()
profiler.disable()
from XFELDetAna.util import env

env.iprofile = cluster_profile
from XFELDetAna import xfelpyanatools as xana
from XFELDetAna import xfelpycaltools as xcal
from XFELDetAna.plotting.util import prettyPlotting

prettyPlotting=True

if sequences[0] == -1:
    sequences = None

# Calibration Database Settings, and Some Initial Run Parameters & Paths:
display(Markdown('### Initial Settings and Paths'))
pixels_x = 1024 # rows of pnCCD in pixels 
pixels_y = 1024 # columns of pnCCD in pixels
print(f"pnCCD size is: {pixels_x}x{pixels_y} pixels.")
print(f'Calibration database interface to be used: {cal_db_interface}')

proposal = list(filter(None, in_folder.strip('/').split('/')))[-2]
file_loc =f'Proposal: {proposal}, Run: {run}'
print(file_loc)

report = get_report(out_folder)
print("report_path:", report)

# Paths to the data:
ped_dir = "{}/r{:04d}".format(out_folder, run)
ped_dir_ctrl = "{}/r{:04d}".format(in_folder, run)
fp_name = path_template.format(run, karabo_da)
fp_path = '{}/{}'.format(ped_dir, fp_name)
h5path = h5path.format(karabo_id, receiver_id)
print("HDF5 path to data: {}\n".format(h5path))

# Output Folder Creation:
os.makedirs(out_folder, exist_ok=True)

# Run's creation time:
if creation_time:
    try:
        creation_time = datetime.datetime.strptime(creation_time, '%Y-%m-%d %H:%M:%S.%f')
    except Exception as e:
        print(f"creation_time value error: {e}." 
               "Use same format as YYYY-MM-DD HR:MN:SC.ms e.g. 2019-07-04 11:02:41.00/n")
        creation_time = None
        print("Given creation time will not be used.")
else:
    creation_time = None

if not creation_time and use_dir_creation_date:
    creation_time = get_dir_creation_date(in_folder, run)

print(f"Creation time: {creation_time}")

# Reading all sequences of the run:
file_list = []
total_sequences = 0
fsequences = []

if sequences is None:
    file_list = glob.glob(fp_path.format(0).replace('00000', '*'))
    file_list = sorted(file_list, key = lambda x: (len (x), x))
    total_sequences = len(file_list)
    fsequences = range(total_sequences)
else:
    for seq in sequences:
        abs_entry = fp_path.format(seq)
        if os.path.isfile(abs_entry):
            file_list.append(abs_entry)
            total_sequences += 1
            fsequences.append(seq)

sequences = fsequences
print(f"This run has a total number of {total_sequences} sequences.\n")

display(Markdown('### List of Files to be Processed'))
print("Reading data from the following files:")
for i in range(len(file_list)):
    print(file_list[i])

display(Markdown('### Detector Parameters'))

# extract slow data
if karabo_da_control:
    ctrl_fname = os.path.join(ped_dir_ctrl, path_template_ctrl.format(run, karabo_da_control)).format(sequences[0])
    ctrl_path = h5path_ctrl.format(karabo_id)
    mdl_ctrl_path = f"/CONTROL/{karabo_id}/MDL/"

    (bias_voltage, gain,
     fix_temperature_top,
     fix_temperature_bot) = extract_slow_data(karabo_id, karabo_da_control, ctrl_fname, ctrl_path,
                                              mdl_ctrl_path, bias_voltage, gain,
                                              fix_temperature_top, fix_temperature_bot)

print(f"Bias voltage is {bias_voltage:0.2f} V.")
print(f"Detector gain is set to {int(gain)}.")
print(f"Top pnCCD sensor is at temperature of {fix_temperature_top:0.2f} K")
print(f"Bottom pnCCD sensor is at temperature of {fix_temperature_bot:0.2f} K")
print(f"Photon energy is {photon_energy} keV") 

# The ADU values correspoding to the boundaries of the first peak region are used as cti_limit_low and 
# cti_limit_high:
gain_k = [k for k, v in VALID_GAINS.items() if v == gain][0]
if gain_k == 'a':
    cti_limit_low = 300 # lower limit of cti
    cti_limit_high = 8200 # higher limit of cti
    max_points = 100000 # maximum data value
elif gain_k == 'b':
    cti_limit_low = 500 
    cti_limit_high = 2000 
    max_points = 100000 
elif gain_k == 'c':
    cti_limit_low = 100 
    cti_limit_high = 500 
    max_points = 100000 
elif gain_k == 'd':
    if bias_voltage <= 400.:
        cti_limit_low = 50 
        cti_limit_high = 150
        max_points = 20000  # ccoords.shape for each quadrant?
    elif bias_voltage > 400.:
        cti_limit_low = 50 
        cti_limit_high = 120
        max_points = 100000  
elif gain_k == 'e':
    if bias_voltage <= 400.:
        cti_limit_low = 10 
        cti_limit_high = 150 
        max_points = 8000  
    elif bias_voltage > 400.:
        cti_limit_low = 13 
        cti_limit_high = 45 
        max_points = 100000              
else:
    cti_limit_low = 13 
    cti_limit_high = 200 
    max_points = 100000 

# Sensor size and block size definitions (important for common mode and other corrections):

sensorSize = [pixels_x, pixels_y]
blockSize = [sensorSize[0]//2, sensorSize[1]//2] # Sensor area will be analysed according to blocksize
xcal.defaultBlockSize = blockSize
memoryCells = 1 # pnCCD has 1 memory cell

# pattern kernels as defined in pyDetLib:
# possibly at pydetlib/lib/src/XFELDetAna/util/pattern_kernels.py

kernel = {}

kernel["203"] = np.array([[0, 1, 0],
                          [0, 1, 0],
                          [0, 0, 0]], dtype=np.float)[:,:,np.newaxis]

kernel["201"] = np.array([[0, 0, 0],
                          [0, 1, 1],
                          [0, 0, 0]], dtype=np.float)[:,:,np.newaxis]

kernel["202"] = np.array([[0, 0, 0],
                          [0, 1, 0 ],
                          [0, 1, 0]], dtype=np.float)[:,:,np.newaxis]


kernel["200"] = np.array([[0, 0, 0],
                          [1, 1, 0 ],
                          [0, 0, 0]], dtype=np.float)[:,:,np.newaxis]

def separateDouble(pdat, alldat, p):
    """
    Just a helper function to separate doubles back out
    
    :param pdat: the pattern indices
    :param alldat: the corrected data
    :param p: the pattern indice to work on
    """
    
    cdat = np.zeros_like(alldat)
    
    kstack = kernel[str(p)]
    idx = pdat == p
    
    p1 = alldat[idx]
    p2 = None

    if np.roll(kstack, 1, axis=0)[1,1] == 1:
        idx = np.roll(pdat, -1, axis = 0) == p
        p2 = alldat[idx]
    if np.roll(kstack, -1, axis=0)[1,1] == 1:
        idx = np.roll(pdat, 1, axis = 0) == p
        p2 = alldat[idx]
    if np.roll(kstack, 1, axis=1)[1,1] == 1:
        idx = np.roll(pdat, -1, axis = 1) == p
        p2 = alldat[idx]
    if np.roll(kstack, -1, axis=1)[1,1] == 1:
        idx = np.roll(pdat, 1, axis = 1) == p
        p2 = alldat[idx]
        
    r = np.zeros((p1.shape[0], 2), np.float)
    r[:,0] = np.abs(p1)-np.abs(p2)
    r[:,1] = np.abs(p2)
   
    return r

# actual data reading happens here:

# separated doubles will be held in these lists
up_doubles = []  # up major partner
down_doubles = []  # down major partner
left_doubles = []  # left major partner
right_doubles = []  # right major partner

# A struct in the C programming language is a composite data type declaration that defines a physically grouped 
# list of variables under one name in a block of memory, allowing the different variables to be accessed via a 
# single pointer or by the struct declared name which returns the same address. 

# structs for CTI and gain evaluation:
# we work by quadrant:
dstats = {'upper left': {'coords': [0, pixels_x//2, 0, pixels_y//2],  # coordinates of this quadrant
                         'row_coords': [],  # will hold our data's row coordinates
                         'col_coords': [],  # will hold our data's column coordinates
                         'col_dir': 1,  # 1 if columns read out to left, -1 if read out to right
                         'vals': [],  # will hold our data values
                         },
          'lower left': {'coords':[0, pixels_x//2, pixels_y//2, pixels_y],
                         'row_coords': [],
                         'col_coords': [],
                         'col_dir': 1,
                         'vals': [],
                         },
          'upper right': {'coords':[pixels_x//2, pixels_x, 0, pixels_y//2],
                         'row_coords': [],
                         'col_coords': [],
                         'col_dir': -1,
                         'vals': [],
                         },
          'lower right': {'coords': [pixels_x//2, pixels_x, pixels_y//2, pixels_y],
                         'row_coords': [],
                         'col_coords': [],
                         'col_dir': -1,
                         'vals': [],
                         },
    
} 

# structure that identifies and holds pattern statistics
pstats = {'singles': {'p': 100,},
          'first singles': {'p': 101,},
          'clusters': {'p': 1000,},
          'left major doubles': {'p': 200,},
          'right major doubles': {'p': 201,},
          'up major doubles': {'p': 202,},
          'down major doubles': {'p': 203,},
          'left major triples': {'p': 300,},
          'right major triples': {'p': 301,},
          'up major triples': {'p': 302,},
          'down major triples': {'p': 303,},
          'left major quads': {'p': 400,},
          'right major quads': {'p': 401,},
          'up major quads': {'p': 402,},
          'down major quads': {'p': 403,},
}

# total number of patterns
tpats = 0

for k, f in enumerate(file_list):
    with h5py.File(f, 'r') as infile:
        data = infile[h5path+"/pixels_classified"][()]
        patterns = infile[h5path+"/patterns"][()]
        bpix = infile[h5path+"/mask"][()]
        
        # separating out doubles
        up_doubles.append(separateDouble(patterns, data, 202))
        down_doubles.append(separateDouble(patterns, data, 203))
        right_doubles.append(separateDouble(patterns, data, 201))
        left_doubles.append(separateDouble(patterns, data, 200))
        
        # creating pattern statistics
        for pstr, entry in pstats.items():
            cpat = patterns == entry["p"]
            entry["counts"] = np.count_nonzero(cpat)
            hist = entry.get("hist", 0)
            hist += np.sum(cpat, axis=0)
            entry["hist"] = hist 
        
        # for normalization of relative occurences
        tpats += np.count_nonzero((patterns > 0) & (patterns < 1000))
        
        # creating coordinates needed for CTI and relative gain fitting
        xv, yv = np.meshgrid(np.arange(data.shape[0]), np.arange(data.shape[1]))
        xv = np.repeat(xv[...,None], data.shape[2], axis=2)
        yv = np.repeat(yv[...,None], data.shape[2], axis=2)

        # we take all single values
        all_singles = np.zeros_like(data)
        all_singles[(patterns == 100) | (patterns == 101)] = data[(patterns == 100) | (patterns == 101)]
        #all_singles[bpix != 0] = np.nan
        
        # and then save data and coordinates quadrant-wise
        for key, entry in dstats.items():
            co = entry["coords"]
            cd = entry["col_dir"]
            asd = all_singles[:, co[2]:co[3], co[0]:co[1]]
            # taking column direction into account - right hemisphere reads to right, left hemisphere reads to left
            if cd == 1:
                yv, xv = np.meshgrid(np.arange(asd.shape[1]), np.arange(asd.shape[2])+co[2]) # row counting 
                                                                                             # downwards
            else:
                yv, xv = np.meshgrid(np.arange(asd.shape[1], 0, -1), np.arange(asd.shape[2])+co[2]) # row counting 
                                                                                                    # downwards
            xv = np.repeat(xv[None,...], data.shape[0], axis=0)
            yv = np.repeat(yv[None,...], data.shape[0], axis=0)
            entry['row_coords'].append(xv[asd > 0].flatten())
            entry['col_coords'].append(yv[asd > 0].flatten())
            entry['vals'].append(asd[asd > 0].flatten())
        
        #break

up_doubles = np.concatenate(up_doubles, axis=0)
down_doubles = np.concatenate(down_doubles, axis=0)
left_doubles = np.concatenate(left_doubles, axis=0)
right_doubles = np.concatenate(right_doubles, axis=0)

for key, entry in dstats.items():
    entry['row_coords'] = np.concatenate(entry['row_coords'])
    entry['col_coords'] = np.concatenate(entry['col_coords'])
    entry['vals'] = np.concatenate(entry['vals'])

table = PrettyTable()
table.field_names = ["Pattern type", "Absolute Number", "Relative Occurence"]
for key, entry in pstats.items():
    table.add_row([key, entry["counts"], "{:0.2f}".format(float(entry["counts"])/tpats)])
    
print(table)

display(Markdown('### Event Patterns'))

fig = plt.figure(1, (15., 15.))
grid = ImageGrid(fig, 111,  # similar to subplot(111)
                 nrows_ncols=(5, 4),  # creates 2x2 grid of axes
                 axes_pad=0.3,  # pad between axes in inch.
                 )

i = 0
for key, entry in pstats.items():
    d = entry["hist"]
    # The AxesGrid object works as a list of axes.
    grid[i].imshow(d, vmin=0, vmax=2.*np.mean(d[d != 0]), cmap="gray_r") 
    grid[i].set_ylabel("Row")
    grid[i].annotate(key, (20, -30), color="red", annotation_clip=False)
    i += 1
    if i == 3: # For the empty panel
        i += 1
    if i == 16: # I actually don't know how this is working and why only one panel has x-label
        grid[i].set_xlabel("Column")

fig = plt.figure(figsize=(8, 8))
ax = fig.add_subplot(111)
ax.scatter(np.abs(up_doubles[:,0]), np.abs(up_doubles[:,1]), 0.5, alpha=0.5, color='b', label="up doubles")
ax.scatter(np.abs(down_doubles[:,1]), np.abs(down_doubles[:,0]), 0.5, alpha=0.5, color='g', label="down doubles")
ax.set_xlabel('Up Doubles')
ax.set_ylabel('Down Doubles')
ax.legend()
plt.show()

fig = plt.figure(figsize=(8,8))
ax = fig.add_subplot(111)
ax.scatter(np.abs(left_doubles[:,0]), np.abs(left_doubles[:,1]), 0.5, alpha=0.5, color='b', label="left doubles")
ax.scatter(np.abs(right_doubles[:,1]), np.abs(right_doubles[:,0]), 0.5, alpha=0.5, color='g', label="right doubles")
ax.semilogx()
ax.semilogy()  
ax.set_xlabel('Left Doubles')
ax.set_ylabel('Right Doubles')
ax.legend()
plt.show()

display(Markdown('### Display of Signal along the Columns'))

fig = plt.figure(figsize=(15., 15.))
grid = AxesGrid(fig, 111,  # similar to subplot(111)
                 nrows_ncols=(2, 2),  # creates 2x2 grid of axes
                 axes_pad=0.3,  # pad between axes in inch.
                 aspect=False,
                 )

ais = [0, 2, 1, 3]
ai = 0
for key, entry in dstats.items():
    i = ais[ai]
    rcoords = entry["row_coords"][:max_points]
    ccoords = entry["col_coords"][:max_points]
    avals = entry["vals"][:max_points]
    
    idx = (avals > cti_limit_low) & (avals < cti_limit_high)


    xax = ccoords[~idx]
    yax = rcoords[~idx]
    grid[i].scatter(xax, avals[~idx], 1, color="grey", alpha=0.2) # These are the data points which ar outside the
                                                                  # valid CTI limits


    xax = ccoords[idx]
    yax = rcoords[idx]
    pvals = avals[idx]

    unique_coords = np.unique(xax)
    mean_signal = np.asarray([np.mean(pvals[xax==jdx]) for jdx in unique_coords])

    cmap = matplotlib.cm.get_cmap('summer')
    colors = cmap(yax/np.max(yax))
    grid[i].scatter(xax, pvals, 1, color=colors)
    grid[i].scatter(unique_coords, mean_signal, 3, color='tomato')
    grid[i].set_ylabel("Signal Value (ADU)")
    grid[i].set_xlabel("Column Coordinate")
    grid[i].annotate(key, (0.1, 1.02), xycoords="axes fraction", color="red", annotation_clip=False)
    if entry["col_dir"] == -1:  # columns reading out to the right
        grid[i].set_xlim(512, 0)
    l = grid[i].set_ylim(0, 1.25*cti_limit_high)
    ai += 1

display(Markdown('### Display of Signal along the Rows'))
        
fig = plt.figure(figsize=(15., 15.))
grid = AxesGrid(fig, 111,  # similar to subplot(111)
                 nrows_ncols=(2, 2),  # creates 2x2 grid of axes
                 axes_pad=0.3,  # pad between axes in inch.
                 aspect=False,
                 )


ai = 0
for key, entry in dstats.items():
    i = ais[ai]
    rcoords = entry["row_coords"][:max_points]
    ccoords = entry["col_coords"][:max_points]
    avals = entry["vals"][:max_points]
    co = entry["coords"]

    
    yax = ccoords[~idx]
    xax = rcoords[~idx]
    grid[i].scatter(avals[~idx], xax, 1, color="grey", alpha=0.2) # These are the data points which ar outside the
                                                                  # valid CTI limits

    yax = ccoords[idx]
    xax = rcoords[idx]
    pvals = avals[idx]

    unique_coords = np.unique(xax)
    mean_signal = np.asarray([np.mean(pvals[xax==jdx]) for jdx in unique_coords])

    cmap = matplotlib.cm.get_cmap('summer')
    colors = cmap(yax/np.max(yax))
    grid[i].scatter(pvals, xax, 1, color=colors)
    grid[i].scatter(mean_signal, unique_coords, 3, color='tomato')
    grid[i].set_xlabel("Signal Value (ADU)")
    grid[i].set_ylabel("Row Coordinate")
    grid[i].annotate(key, (0.1, 1.02), xycoords="axes fraction", color="red", annotation_clip=False)
    grid[i].set_ylim(co[3], co[2])
    l = grid[i].set_xlim(0, 1.25*cti_limit_high)
    ai += 1

# Fitting is performed in this cell for just a randomly selected sample row as an example:

def show_fit(row_example):
    
    if row_example > 512 or row_example < 0:
        print("Invalid row number. Try again.")
    else:
        max_points = -1
        limit_cols = [0, 512]

        for key, entry in dstats.items():

            rcoords = entry["row_coords"][:max_points]
            ccoords = entry["col_coords"][:max_points]
            avals = entry["vals"][:max_points]
            co = entry["coords"]
            
            idx = (avals > cti_limit_low) & (avals < cti_limit_high) & (ccoords >= limit_cols[0]) & \
                  (ccoords <= limit_cols[1])
            xax = rcoords[idx]-co[2]
            yax = ccoords[idx] # columns
            avalsf = avals[idx]
            min_data = 5

            ms = np.zeros(sensorSize[0]//2, np.float)
            bs = np.zeros(sensorSize[0]//2, np.float)
            mErr = np.zeros(sensorSize[0]//2, np.float)
            bErr = np.zeros(sensorSize[0]//2, np.float)

            axes = None
            fig = None
            fig, ax = plt.subplots(1, sharex=True, figsize=(10, 5))

            def linFunc(x, *p):
                return p[1] * x + p[0]

            for i in [row_example]:
                idx = (xax == i)
                zm = avalsf[idx]
                zx = yax[idx]

                def linWrap(m, b):
                    return np.sum(((linFunc(zx, b, -m) - zm)) ** 2)

                pparm = dict(throw_nan=False, pedantic=False, print_level=0)
                pparm["m"] = 0
                if gain == 1.0:
                    pparm["limit_m"] = (0.9, 1)
                else:
                    pparm["limit_m"] = (0, 1)
                pparm["b"] = np.nanmean(zm)

                mini = Minuit(linWrap, **pparm)
                fmin = mini.migrad()

                if fmin[0]['is_valid'] and fmin[0]['has_covariance'] and (min_data is None or zx.size > min_data):

                    ferr = mini.minos(maxcall = 10000)
                    res = mini.values

                    ms[i] = res["m"]
                    bs[i] = res["b"]

                    mErr[i] = max(-ferr['m']['lower'], ferr['m']['upper'])
                    bErr[i] = max(-ferr['b']['lower'], ferr['b']['upper'])

                    zx, uidx = np.unique(zx, return_index=True)
                    zm = zm[uidx]

                    yt = linFunc(zx, res["b"], res["m"])

                    ax.scatter(zx, zm, 3) # point size = 3
                    ax.plot(zx, -ms[i] * zx + bs[i], color='r',
                               label = 'CTI fit for row {} in {} quadrant'.format(row_example, key))
                    ax.annotate('valid: {}'.format(fmin[0]['is_valid']), 
                                   xy=(0.2, 0.3), xycoords='axes fraction')
                    ax.set_xlabel("Column Coordinate")
                    ax.set_ylabel("Singles Events (between CTI_low and CTI_high)")
                    ax.set_title("Events Are Only Shown for One Row")
                    ax.set_xlim(0, 512)
                    ax.legend()

display(Markdown('### Example of CTI Fit for a Specific Row'))
fig = show_fit(365)

# Fitting is performed in this cell for the entire detector to get the CTI and relative gain:

max_points = -1
limit_cols = [0, 512]
counter = 0 # Total number of bad rows (over all quadrants) for which the fit does not converge

for key, entry in dstats.items():
    
    rcoords = entry["row_coords"][:max_points]
    ccoords = entry["col_coords"][:max_points]
    avals = entry["vals"][:max_points]
    co = entry["coords"]
    idx = (avals > cti_limit_low) & (avals < cti_limit_high) & (ccoords >= limit_cols[0]) & \
          (ccoords <= limit_cols[1])
    xax = rcoords[idx]-co[2]
    yax = ccoords[idx]
    avalsf = avals[idx]
    min_data = 5

    ms = np.zeros(sensorSize[0]//2, np.float)
    bs = np.zeros(sensorSize[0]//2, np.float)
    mErr = np.zeros(sensorSize[0]//2, np.float)
    bErr = np.zeros(sensorSize[0]//2, np.float)

    def linFunc(x, *p):
            return p[1] * x + p[0]
        
    for i in range(sensorSize[0]//2):
        idx = (xax == i)
        zm = avalsf[idx]
        zx = yax[idx]
        
        def linWrap(m, b):
            return np.sum(((linFunc(zx, b, -m) - zm)) ** 2)
        
        pparm = dict(throw_nan=False, pedantic=False, print_level=0)
        pparm["m"] = 0
        if gain == 1.0:
            pparm["limit_m"] = (0.9, 1)
        else:
            pparm["limit_m"] = (0, 1)
        pparm["b"] = np.nanmean(zm)

        mini = Minuit(linWrap, **pparm)
        fmin = mini.migrad()
        
        if fmin[0]['is_valid'] and fmin[0]['has_covariance'] and (min_data is None or zx.size > min_data):
    
            ferr = mini.minos(maxcall = 10000)
            res = mini.values

            ms[i] = res["m"]
            bs[i] = res["b"]

            mErr[i] = max(-ferr['m']['lower'], ferr['m']['upper'])
            bErr[i] = max(-ferr['b']['lower'], ferr['b']['upper'])

            zx, uidx = np.unique(zx, return_index=True)
            zm = zm[uidx]
            
            yt = linFunc(zx, res["b"], res["m"])
        
        else:
            counter += 1 
                
    cti = ms / bs
    ctiErr = np.sqrt((1. / bs * mErr) ** 2 + (ms / bs ** 2 * bErr) ** 2)
    
    relGain = bs/np.nanmean(bs) 
    relGainErr = np.sqrt((1. / np.nanmean(bs) * bErr) ** 2 + (bs / np.nanmean(bs) ** 2 * np.std(bs)) ** 2)

    entry['cti'] = cti
    entry['cti_err'] = ctiErr
    entry['rel_gain'] = relGain
    entry['rel_gain_err'] = relGainErr

print("Fitting is performed for the entire detector, and CTI and relative gain are calculated.")
print(f"Total number of fits which did not converge: {counter}")

display(Markdown('### CTI per Quadrant'))
        
fig = plt.figure(figsize=(15., 15.))
grid = AxesGrid(fig, 111,  # similar to subplot(111)
                 nrows_ncols=(2, 2),  # creates 2x2 grid of axes
                 axes_pad=0.3,  # pad between axes in inch.
                 aspect=False,
                 )

ais = [0, 2, 1, 3]
ai = 0
for key, entry in dstats.items():
    i = ais[ai]
    co = entry["coords"]
    cti = entry["cti"]
    x = np.arange(cti.size)+co[2]
    grid[i].scatter(cti, x, 1., color='k')
    if gain == 1.0:
        grid[i].set_xlim(1e-5, 1e-2)
    else:
        grid[i].set_xlim(1e-8, 1e-2)
    grid[i].semilogx()
    grid[i].set_xlabel("CTI")
    grid[i].set_ylabel("Row Coordinate")
    grid[i].annotate(key, (0.1, 1.02), xycoords="axes fraction", color="red", annotation_clip=False)
    grid[i].set_ylim(co[3], co[2])
    ai += 1

display(Markdown('### Uncertainty in CTI per Quadrant'))
        
fig = plt.figure(figsize=(15., 15.))
grid = AxesGrid(fig, 111,  # similar to subplot(111)
                 nrows_ncols=(2, 2),  # creates 2x2 grid of axes
                 axes_pad=0.3,  # pad between axes in inch.
                 aspect=False,
                 )

ais = [0, 2, 1, 3]
ai = 0
for key, entry in dstats.items():
    i = ais[ai]
    co = entry["coords"]
    cti_err = entry["cti_err"]
    x = np.arange(cti_err.size)+co[2]
    grid[i].scatter(cti_err, x, 1., color='k')
    if gain == 1.0:
        grid[i].set_xlim(1e-12, 1e-8)
    else:
        grid[i].set_xlim(1e-10, 1e-5)
    grid[i].semilogx()
    grid[i].set_xlabel("Uncertainty in CTI")
    grid[i].set_ylabel("Row Coordinate")
    grid[i].annotate(key, (0.1, 1.02), xycoords="axes fraction", color="red", annotation_clip=False)
    grid[i].set_ylim(co[3], co[2])
    ai += 1

display(Markdown('### Relative Gain per Quadrant'))

fig = plt.figure(figsize=(15., 15.))
grid = AxesGrid(fig, 111,  # similar to subplot(111)
                 nrows_ncols=(2, 2),  # creates 2x2 grid of axes
                 axes_pad=0.3,  # pad between axes in inch.
                 aspect=False,
                 )

ais = [0, 2, 1, 3]
ai = 0
for key, entry in dstats.items():
    i = ais[ai]
    co = entry["coords"]
    rgain = entry["rel_gain"]
    x = np.arange(rgain.size)+co[2]
    grid[i].scatter(rgain, x, 1., color='k')
    grid[i].set_xlim(0.5, 1.5)
    grid[i].set_xlabel("Relative Gain")
    grid[i].set_ylabel("Row Coordinate")
    grid[i].annotate(key, (0.1, 1.02), xycoords="axes fraction", color="red", annotation_clip=False)
    grid[i].set_ylim(co[3], co[2])
    ai += 1

display(Markdown('### Uncertainty in Relative Gain per Quadrant'))

fig = plt.figure(figsize=(15., 15.))
grid = AxesGrid(fig, 111,  # similar to subplot(111)
                 nrows_ncols=(2, 2),  # creates 2x2 grid of axes
                 axes_pad=0.3,  # pad between axes in inch.
                 aspect=False,
                 )

ais = [0, 2, 1, 3]
ai = 0
for key, entry in dstats.items():
    i = ais[ai]
    co = entry["coords"]
    rgain_err = entry["rel_gain_err"]
    x = np.arange(rgain_err.size)+co[2]
    grid[i].scatter(rgain_err, x, 1., color='k')
    grid[i].set_xlim(0, 0.2)
    grid[i].set_xlabel("Relative Gain Uncertainty")
    grid[i].set_ylabel("Row Coordinate")
    grid[i].annotate(key, (0.1, 1.02), xycoords="axes fraction", color="red", annotation_clip=False)
    grid[i].set_ylim(co[3], co[2])
    ai += 1

# Sum of the total number of bad columns in each quadrant should be equal to the value of "counter"

didnt_converge = [] # list of columns in each quadrant for which the fit did not converge. Here, we assume the 
                    # number of columns in each quadrant always starts from 0 and ends at 512
    
print("Rows for which fits did not converge:\n")

for key, entry in dstats.items():
    rel_gain = entry["rel_gain"] 
    print("{}:".format(key))
    for index, element in enumerate(rel_gain):
        if element == 0:
            didnt_converge.append(index)
    print("Total number: {}".format(len(didnt_converge)))
    print("These rows are as follows:")
    print(didnt_converge)
    print("\n")
    didnt_converge = []

ais = [0, 2, 1, 3]
ai = 0
for key, entry in dstats.items():
    i = ais[ai]
    cti = entry["cti"]
    idx = np.where(np.isnan(cti))
    for i in idx[0]:
        cti[i] = np.nanmean(cti)
    
    rel_gain = entry["rel_gain"] 
    for index, element in enumerate(rel_gain):
        if element == 0:
            rel_gain[index] = 1

# Calculating CTE and gain maps:

ctemap = np.zeros((sensorSize[0], sensorSize[1]))
ctimap = np.zeros((sensorSize[0], sensorSize[1]))
gainmap = np.zeros((sensorSize[0], sensorSize[1]))
productmap = np.zeros((sensorSize[0], sensorSize[1]))

for key, entry in dstats.items():
    co = entry["coords"]
    rel_gain = entry["rel_gain"]        
    cti = entry["cti"]
    if entry["col_dir"] == 1:
        # plus sign is so that the charge transfer inefficiency is 100% at the center of the CCD, which is the 
        # farthest distance away from the readout. If using a minus sign, you can set CTE = 100% at the readout 
        # => most efficient on the readout and less efficient towards the center of the chip:
        quadcte = (1-cti[:, None])**np.repeat(np.arange(512)[None, :], 512, axis=0)
    else:
        quadcte = (1-cti[:, None])**np.repeat(np.arange(512, 0, -1)[None, :], 512, axis=0)
    
    #quadcti = cti[:, None] 
    quadgain = rel_gain[:, None]
    quadproduct = quadgain*quadcte
    
    ctemap[co[2]:co[3], co[0]:co[1]] = quadcte
    gainmap[co[2]:co[3], co[0]:co[1]] = quadgain 
    productmap[co[2]:co[3], co[0]:co[1]] = quadproduct

display(Markdown('### CTE and Relative Gain Maps'))

fig = xana.heatmapPlot(ctemap, figsize=(8, 8), x_label='Columns', y_label='Rows', 
                       lut_label='Charge Transfer Efficiency', 
                       aspect=1, cmap = 'viridis',x_range=(0, pixels_y), y_range=(0, pixels_x),
                       panel_x_label='Along Rows', panel_y_label='Along Columns', 
                       title = 'CTE Map for pnCCD (Gain = 1/{})'.format(int(gain)))
    
fig = xana.heatmapPlot(gainmap, figsize=(8, 8), x_label='Columns', y_label='Rows', 
                       lut_label='Relative Gain', vmin=0.8, vmax=1.2,
                       aspect=1, x_range=(0, pixels_y), y_range=(0, pixels_x),
                       panel_top_low_lim = 0.9, panel_top_high_lim = 1.1,
                       panel_side_low_lim = 0.9, panel_side_high_lim = 1.1,
                       panel_x_label='Along Rows', panel_y_label='Along Columns', 
                       title = 'Relative Gain Map for pnCCD (Gain = 1/{})'.format(int(gain)))

fig = xana.heatmapPlot(productmap, figsize=(8, 8), x_label='Columns', y_label='Rows', 
                       lut_label='Relative Gain*CTI', 
                       aspect=1, x_range=(0, pixels_y), y_range=(0, pixels_x), vmin=0.8, vmax=1.2, 
                       panel_x_label='Along Rows', panel_y_label='Along Columns', 
                       panel_top_low_lim = 0.9, panel_top_high_lim = 1.1, panel_side_low_lim = 0, 
                       panel_side_high_lim = 1.5, title = 'Relative Gain*CTE Map for pnCCD (Gain = 1/{})'
                       .format(int(gain)))

constant_maps = {
    'RelativeGain' : gainmap,
    'CTE' : ctemap}

md = None
for cname in constant_maps.keys():
    metadata = ConstantMetaData()
    det = Constants.CCD(DetectorTypes.pnCCD)
    const = getattr(det, cname)()
    const.data = constant_maps[cname].data
    metadata.calibration_constant = const
    
    # setting the operating condition
    condition = Conditions.Illuminated.CCD(bias_voltage=bias_voltage,
                                           integration_time=integration_time,
                                           photon_energy=photon_energy,
                                           gain_setting=gain,
                                           temperature=fix_temperature_top,
                                           pixels_x=pixels_x,
                                           pixels_y=pixels_y)
    # This should be used in case of running notebook 
    # by a different method other than myMDC which already
    # sends CalCat info.
    # TODO: Set db_module to "" by default in the first cell
    if not db_module:
        db_module = get_pdu_from_db(karabo_id, karabo_da, const,
                                    condition, cal_db_interface,
                                    snapshot_at=creation_time)[0]
    if db_output:
        md = send_to_db(db_module, karabo_id, const, condition,
                        file_loc, report,
                        cal_db_interface,
                        creation_time=creation_time,
                        timeout=cal_db_timeout)
        
    if local_output:
        md = save_const_to_h5(db_module, karabo_id, 
                              const, condition, const.data,
                              file_loc, report,
                              creation_time, out_folder)
        print(f"Calibration constant {cname} is stored to {out_folder}.")
        
print("\nGenerated constants with conditions:\n")
print(f"• bias_voltage: {bias_voltage}\n• photon_energy: {photon_energy}\n"
      f"• top_temperature: {fix_temperature_top}\n• integration_time: {integration_time}\n"
      f"• gain_setting: {gain}\n• creation_time: {creation_time}\n")