Merge branch 'fix/path_inset_fastccd_correct' into 'master'

add path_inset to balance sequences See merge request detectors/pycalibration!229

Merge branch 'fix/path_inset_fastccd_correct' into 'master'
add path_inset to balance sequences See merge request detectors/pycalibration!229
40a0ea42 · Karim Ahmed · c603ba08 · 1892b687 · 40a0ea42
Commit 40a0ea42 authored 5 years ago by Karim Ahmed
--- a/notebooks/FastCCD/CorrectionNotebook_NewDAQ_FastCCD_NBC.ipynb
+++ b/notebooks/FastCCD/CorrectionNotebook_NewDAQ_FastCCD_NBC.ipynb
@@ -53,7 +53,7 @@
    "flipped_between = [\"2019-02-01\", \"2019-04-02\"] # detector was flipped during this timespan\n",
    "temp_limits = 5 # limits within which temperature is considered the same\n",
    "\n",
-    "def balance_sequences(in_folder, run, sequences, sequences_per_node):\n",
+    "def balance_sequences(in_folder, run, sequences, sequences_per_node, path_inset):\n",
    "    import glob\n",
    "    import re\n",
    "    import numpy as np\n",

 %% Cell type:markdown id: tags:

 # FastCCD Data Correction ##

 Authors: I. Klačková, S. Hauf, Version 1.0

 The following notebook provides correction of images acquired with the FastCCD.

 %% Cell type:code id: tags:

 ``` python
 in_folder = "/gpfs/exfel/exp/SCS/201802/p002170/raw/" # input folder, required
 out_folder = '/gpfs/exfel/data/scratch/xcal/test/' # output folder, required
 path_template = 'RAW-R{:04d}-{}-S{{:05d}}.h5' #  path template in hdf5 file
 path_inset = 'DA05'
 run = 277 # run number
 h5path = '/INSTRUMENT/SCS_CDIDET_FCCD2M/DAQ/FCCD:daqOutput/data/image' # path in HDF5 file
 h5path_t = '/CONTROL/SCS_CDIDET_FCCD2M/CTRL/LSLAN/inputA/crdg/value' # temperature path in HDF5 file
 h5path_cntrl = '/RUN/SCS_CDIDET_FCCD2M/DET/FCCD'  # path to control data
 cluster_profile = "noDB" #ipcluster profile to use
 cpuCores = 16 #Specifies the number of running cpu cores
 operation_mode = "FF" # FS stands for frame-store and FF for full-frame opeartion
 split_evt_primary_threshold = 7. # primary threshold for split event classification in terms of n sigma noise
 split_evt_secondary_threshold = 4. # secondary threshold for split event classification in terms of n sigma noise
 split_evt_mip_threshold = 1000. # MIP threshold for event classification
 cal_db_interface = "tcp://max-exfl016:8015#8025" # calibration DB interface to use
 cal_db_timeout = 300000000 # timeout on caldb requests
 sequences = [-1] # sequences to correct, set to -1 for all, range allowed
 chunk_size_idim = 1 # H5 chunking size of output data
 overwrite = True # overwrite existing files
 do_pattern_classification = True # classify split events
 sequences_per_node = 1 # sequences to correct per node
 limit_images = 0 # limit images per file
 correct_offset_drift = False  # correct for offset drifts
 use_dir_creation_date = True # use dir creation data for calDB queries
 time_offset_days = 0 # offset in days for calibration parameters
 photon_energy_gain_map = 2. # energy in keV
 fix_temperature = 0. # fix temperature to this value, set to 0 to use slow control value
 flipped_between = ["2019-02-01", "2019-04-02"] # detector was flipped during this timespan
 temp_limits = 5 # limits within which temperature is considered the same

-def balance_sequences(in_folder, run, sequences, sequences_per_node):
+def balance_sequences(in_folder, run, sequences, sequences_per_node, path_inset):
    import glob
    import re
    import numpy as np
    if sequences[0] == -1:
        sequence_files = glob.glob("{}/r{:04d}/*{}-S*.h5".format(in_folder, run, path_inset))
        seq_nums = set()
        for sf in sequence_files:
            seqnum = re.findall(r".*-S([0-9]*).h5", sf)[0]
            seq_nums.add(int(seqnum))
        seq_nums -= set(sequences)
        nsplits = len(seq_nums)//sequences_per_node+1
        while nsplits > 8:
            sequences_per_node += 1
            nsplits = len(seq_nums)//sequences_per_node+1
            print("Changed to {} sequences per node to have a maximum of 8 concurrent jobs".format(sequences_per_node))
        return [l.tolist() for l in np.array_split(list(seq_nums), nsplits)]
    else:
        return sequences
 ```

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

 import numpy as np
 import h5py
 import matplotlib.pyplot as plt
 from iminuit import Minuit

 import time
 import copy
 import os

 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 datetime import timedelta

 %matplotlib inline

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

 offset_correction_args = (0.2459991787617141, 243.21639920846485)
 t_base = 247.82

 if "#" in cal_db_interface:
    prot, serv, ran = cal_db_interface.split(":")
    r1, r2 = ran.split("#")
    cal_db_interface = ":".join(
        [prot, serv, str(np.random.randint(int(r1), int(r2)))])
 ```

 %% Cell type:code id: tags:

 ``` python
 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, path_inset)
 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))
 print("Data is output to: {}".format(out_folder))

 import datetime
 creation_time = None
 if use_dir_creation_date:
    creation_time = get_dir_creation_date(in_folder, run) + timedelta(days=time_offset_days)
 if creation_time:
    print("Using {} as creation time".format(creation_time.isoformat()))

 ```

 %% Cell type:code id: tags:

 ``` python

 sensorSize = [x, y]
 chunkSize = 100 #Number of images to read per chunk
 blockSize = [sensorSize[0]//2, sensorSize[1]//4] #Sensor area will be analysed according to blocksize
 xcal.defaultBlockSize = blockSize
 memoryCells = 1 #FastCCD has 1 memory cell
 #Specifies total number of images to proceed

 commonModeBlockSize = blockSize
 commonModeAxisR = 'row'#Axis along which common mode will be calculated
 run_parallel = True
 profile = False

 temperature_k = 291
 filename = fp_path.format(sequences[0] if sequences else 0)
 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])
    print("Bias voltage is {} V".format(bias_voltage))
    print("Detector gain is set to x{}".format(det_gain))
    print("Detector integration time is set to {}".format(integration_time))
    temperature = np.mean(f[h5path_t])
    temperature_k = temperature + 273.15
    if fix_temperature != 0.:
        temperature_k = fix_temperature
        print("Using fixed temperature")
    print("Mean temperature was {:0.2f} °C / {:0.2f} K at beginning of run".format(temperature, temperature_k))


 if not os.path.exists(out_folder):
    os.makedirs(out_folder)
 elif not overwrite:
    raise AttributeError("Output path exists! Exiting")
 ```

 %% Cell type:code id: tags:

 ``` python
 dirlist = sorted(os.listdir(ped_dir))
 file_list = []
 total_sequences = 0
 fsequences = []
 for entry in dirlist:

    #only h5 file
    abs_entry = "{}/{}".format(ped_dir, entry)
    if os.path.isfile(abs_entry) and os.path.splitext(abs_entry)[1] == ".h5":

        if sequences is None:
            for seq in range(len(dirlist)):

                if path_template.format(run, path_inset).format(seq) in abs_entry:
                    file_list.append(abs_entry)
                    total_sequences += 1
                    fsequences.append(seq)
        else:
            for seq in sequences:

                if path_template.format(run, path_inset).format(seq) in abs_entry:
                    file_list.append(os.path.abspath(abs_entry))
                    total_sequences += 1
                    fsequences.append(seq)
 sequences = fsequences
 ```

 %% Cell type:code id: tags:

 ``` python
 import copy
 from IPython.display import HTML, display, Markdown, Latex
 import tabulate
 print("Processing a total of {} sequence files".format(total_sequences))
 table = []


 for k, f in enumerate(file_list):
    table.append((k, f))
 if len(table):
    md = display(Latex(tabulate.tabulate(table, tablefmt='latex', headers=["#", "file"])))
 ```

 %% Cell type:markdown id: tags:

 As a first step, dark maps have to be loaded.

 %% Cell type:code id: tags:

 ``` python
 offsetMap = None
 badPixelMap = None
 noiseMap = None
 for i, g in enumerate([8, 2, 1]):
    ## offset
    metadata = ConstantMetaData()
    offset = Constants.CCD(DetectorTypes.fastCCD).Offset()
    metadata.calibration_constant = offset

    # set the operating condition
    condition = Conditions.Dark.CCD(bias_voltage=bias_voltage,
                                    integration_time=integration_time,
                                    gain_setting=g,
                                    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



    # specify the version for this constant
    if creation_time is None:
        metadata.calibration_constant_version = Versions.Now(device=device)
        metadata.retrieve(cal_db_interface)
    else:
        metadata.calibration_constant_version = Versions.Timespan(device=device,
                                                                  start=creation_time)
        metadata.retrieve(cal_db_interface, when=creation_time.isoformat(), timeout=3000000)


    if offsetMap is None:
        offsetMap = np.zeros(list(offset.data.shape)+[3], np.float32)
    offsetMap[...,i] = offset.data

    offset_temperature = None
    for parm in condition.parameters:

        if parm.name == "Sensor Temperature":
            offset_temperature = parm.value

    print("Temperature of detector when dark images (gain {}) for offset calculation ".format(g) +
              "were taken at: {:0.2f} K @ {}".format(offset_temperature,
                                                     metadata.calibration_constant_version.begin_at))

    ## noise
    metadata = ConstantMetaData()
    noise = Constants.CCD(DetectorTypes.fastCCD).Noise()
    metadata.calibration_constant = noise

    # set the operating condition
    condition = Conditions.Dark.CCD(bias_voltage=bias_voltage,
                                    integration_time=integration_time,
                                    gain_setting=g,
                                    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

    # specify the version for this constant
    if creation_time is None:
        metadata.calibration_constant_version = Versions.Now(device=device)
        metadata.retrieve(cal_db_interface)
    else:
        metadata.calibration_constant_version = Versions.Timespan(device=device,
                                                                  start=creation_time)
        metadata.retrieve(cal_db_interface, when=creation_time.isoformat(), timeout=3000000)

    if noiseMap is None:
        noiseMap = np.zeros(list(noise.data.shape)+[3], np.float32)
    noiseMap[...,i] = noise.data


    ## bad pixels

    metadata = ConstantMetaData()
    bpix = Constants.CCD(DetectorTypes.fastCCD).BadPixelsDark()
    metadata.calibration_constant = bpix

    # set the operating condition
    condition = Conditions.Dark.CCD(bias_voltage=bias_voltage,
                                    integration_time=integration_time,
                                    gain_setting=g,
                                    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

    # specify the version for this constant
    if creation_time is None:
        metadata.calibration_constant_version = Versions.Now(device=device)
        metadata.retrieve(cal_db_interface)
    else:
        metadata.calibration_constant_version = Versions.Timespan(device=device,
                                                                  start=creation_time)
        metadata.retrieve(cal_db_interface, when=creation_time.isoformat(), timeout=3000000)

    if badPixelMap is None:
        badPixelMap = np.zeros(list(bpix.data.shape)+[3], np.uint32)
    badPixelMap[...,i] = bpix.data

 ```

 %% Cell type:markdown id: tags:

 Loading cti and relative gain values

 %% Cell type:code id: tags:

 ``` python
 ## relative gain

 metadata = ConstantMetaData()
 relgain = Constants.CCD(DetectorTypes.fastCCD).RelativeGain()
 metadata.calibration_constant = relgain

 # set the operating condition
 condition = Conditions.Illuminated.CCD(bias_voltage=bias_voltage,
                                       integration_time=integration_time,
                                       gain_setting=0,
                                       temperature=temperature_k,
                                       pixels_x=1934,
                                       pixels_y=960, photon_energy=photon_energy_gain_map)
 device = Detectors.fastCCD1


 metadata.detector_condition = condition

 # specify the a version for this constant
 metadata.calibration_constant_version = Versions.Now(device=device)
 metadata.retrieve(cal_db_interface)

 relGain = relgain.data[::-1,...]
 ```

 %% Cell type:code id: tags:

 ``` python
 relGainCA = copy.copy(relGain)
 relGainC = relGainCA[:relGainCA.shape[0]//2,...]
 ctiA = np.ones(relGainCA.shape[:2])
 cti = np.ones(relGainC.shape[:2])
 i = 0
 idx = (relGainC[i, :, 0] < 0.9) | (relGainC[i,:,0] > 1.1)
 mn1 = np.nanmean(relGainC[i, ~idx, 0])

 for i in range(1, relGainC.shape[0]):
    idx = (relGainC[i, :, 0] < 0.9) | (relGainC[i,:,0] > 1.1)
    mn2 = np.nanmean(relGainC[i, ~idx, 0])
    cti[i,:] = mn2/mn1
 ctiA[:relGainCA.shape[0]//2,...] = cti

 relGainC = relGainCA[relGainCA.shape[0]//2:,...]


 cti = np.ones(relGainC.shape[:2])
 i = -1
 idx = (relGainC[i, :, 0] < 0.9) | (relGainC[i,:,0] > 1.1)
 mn1 = np.nanmean(relGainC[i, ~idx, 0])

 for i in range(relGainC.shape[0]-1, 1, -1):
    idx = (relGainC[i, :, 0] < 0.9) | (relGainC[i,:,0] > 1.1)
    mn2 = np.nanmean(relGainC[i, ~idx, 0])
    cti[i,:] = mn2/mn1

 ctiA[relGainCA.shape[0]//2:,...] = cti

 relGainCA = copy.copy(relGain)
 relGainC = relGainCA[:relGainCA.shape[0]//2,...]
 for i in range(relGainC.shape[1]):
    idx = (relGainC[:,i, 0] < 0.95) | (relGainC[:,i,0] > 1.05)
    relGainC[idx,i,0] = np.nanmean(relGainC[~idx,i,0])
    relGainC[idx,i,1] = np.nanmean(relGainC[~idx,i,1])
    relGainC[idx,i,2] = np.nanmean(relGainC[~idx,i,2])
 relGainCA[:relGainCA.shape[0]//2,...] = relGainC
 relGainC = relGainCA[relGainCA.shape[0]//2:,...]
 for i in range(relGainC.shape[1]):
    idx = (relGainC[:,i, 0] < 0.95) | (relGainC[:,i,0] > 1.05)
    relGainC[idx,i,0] = np.nanmean(relGainC[~idx,i,0])
    relGainC[idx,i,1] = np.nanmean(relGainC[~idx,i,1])
    relGainC[idx,i,2] = np.nanmean(relGainC[~idx,i,2])
 relGainCA[relGainCA.shape[0]//2:,...] = relGainC
 relGainC = relGainCA*ctiA[...,None]

 relGain = relGainC
 ```

 %% Cell type:code id: tags:

 ``` python
 import dateutil.parser
 flipped_between = [dateutil.parser.parse(d) for d in flipped_between]
 flip_rgain = creation_time >= flipped_between[0] and creation_time <= flipped_between[1]
 flip_rgain &= (metadata.calibration_constant_version.begin_at.replace(tzinfo=None) >= flipped_between[0]
               and metadata.calibration_constant_version.begin_at.replace(tzinfo=None) <= flipped_between[1])
 print("Accounting for flipped detector: {}".format(flip_rgain))
 ```

 %% Cell type:code id: tags:

 ``` python
 #************************Calculators************************#


 cmCorrection = xcal.CommonModeCorrection([x, y],
                                         commonModeBlockSize,
                                         commonModeAxisR,
                                         nCells = memoryCells,
                                         noiseMap = noiseMap,
                                        runParallel=True,
                                        stats=True)

 patternClassifierLH = xcal.PatternClassifier([x//2, y],
                                              noiseMap[:x//2, :],
                                              split_evt_primary_threshold,
                                             split_evt_secondary_threshold,
                                             split_evt_mip_threshold,
                                              tagFirstSingles = 0,
                                              nCells=memoryCells,
                                              cores=cpuCores,
                                              allowElongated = False,
                                              blockSize=[x//2, y],
                                              runParallel=True)



 patternClassifierUH = xcal.PatternClassifier([x//2, y],
                                              noiseMap[x//2:, :],
                                              split_evt_primary_threshold,
                                              split_evt_secondary_threshold,
                                              split_evt_mip_threshold,
                                              tagFirstSingles = 0,
                                              nCells=memoryCells,
                                              cores=cpuCores,
                                              allowElongated = False,
                                              blockSize=[x//2, y],
                                              runParallel=True)




 ```

 %% Cell type:code id: tags:

 ``` python
 #*****************Histogram Calculators******************#

 histCalOffsetCor = xcal.HistogramCalculator([x, y],
                                             bins=500,
                                             range=[-50, 1000],
                                             nCells=memoryCells,
                                             cores=cpuCores,
                                             blockSize=blockSize)

 histCalPcorr = xcal.HistogramCalculator([x, y],
                                         bins=500,
                                         range=[-50, 1000],
                                         nCells=memoryCells,
                                         cores=cpuCores,
                                         blockSize=blockSize)

 histCalPcorrS = xcal.HistogramCalculator([x, y],
                                         bins=500,
                                         range=[-50, 1000],
                                         nCells=memoryCells,
                                         cores=cpuCores,
                                         blockSize=blockSize)
 ```

 %% Cell type:markdown id: tags:

 Applying corrections

 %% Cell type:code id: tags:

 ``` python
 patternClassifierLH._imagesPerChunk = 500
 patternClassifierUH._imagesPerChunk = 500
 patternClassifierLH.debug()
 patternClassifierUH.debug()
 ```

 %% Cell type:code id: tags:

 ``` python
 histCalOffsetCor.debug()
 histCalPcorr.debug()
 ```

 %% Cell type:code id: tags:

 ``` python
 def copy_and_sanitize_non_cal_data(infile, outfile, h5base):

    if h5base.startswith("/"):
        h5base = h5base[1:]
    dont_copy = ['pixels']
    dont_copy = [h5base+"/{}".format(do)
                for do in dont_copy]

    def visitor(k, item):
        if k not in dont_copy:
            if isinstance(item, h5py.Group):
                outfile.create_group(k)
            elif isinstance(item, h5py.Dataset):
                group = str(k).split("/")
                group = "/".join(group[:-1])
                infile.copy(k, outfile[group])

    infile.visititems(visitor)
 ```

 %% Cell type:code id: tags:

 ``` python
 mean_im = None
 single_im = None
 mean_im_cc = None
 single_im_cc = None
 drift_lh = []
 drift_uh = []
 offsetMap = np.squeeze(offsetMap)
 noiseMap = np.squeeze(noiseMap)
 badPixelMap = np.squeeze(badPixelMap)
 relGain = np.squeeze(relGain)
 for k, f in enumerate(file_list):
    with h5py.File(f, 'r', driver='core') as infile:
        out_fileb = "{}/{}".format(out_folder, f.split("/")[-1])
        out_file = out_fileb.replace("RAW", "CORR")
        #out_filed = out_fileb.replace("RAW", "CORR-SC")

        data = None
        noise = None
        try:
            with h5py.File(out_file, "w") as ofile:

                copy_and_sanitize_non_cal_data(infile, ofile, h5path)
                data = infile[h5path+"/pixels"][()]
                nzidx = np.count_nonzero(data, axis=(1,2))
                data = data[nzidx != 0, ...]
                if limit_images > 0:
                    data = data[:limit_images,...]
                oshape = data.shape
                data = np.moveaxis(data, 0, 2)
                ddset = ofile.create_dataset(h5path+"/pixels",
                                             oshape,
                                             chunks=(chunk_size_idim, oshape[1], oshape[2]),
                                             dtype=np.float32)

                ddsetm = ofile.create_dataset(h5path+"/mask",
                                             oshape,
                                             chunks=(chunk_size_idim, oshape[1], oshape[2]),
                                             dtype=np.uint32, compression="gzip")

                ddsetg = ofile.create_dataset(h5path+"/gain",
                                             oshape,
                                             chunks=(chunk_size_idim, oshape[1], oshape[2]),
                                             dtype=np.uint8, compression="gzip")

                gain = np.right_shift(data, 14)

                gain[gain != 0] -= 1

                fstride = 1
                if not flip_rgain:  # rgain was taken during flipped orientation
                    fstride = -1

                data = np.bitwise_and(data, 0b0011111111111111).astype(np.float32)
                omap = np.repeat(offsetMap[...,None,:], data.shape[2], axis=2)
                rmap = np.repeat(relGain[:,::fstride,None,:], data.shape[2], axis=2)
                nmap = np.repeat(noiseMap[...,None,:], data.shape[2], axis=2)
                bmap = np.repeat(badPixelMap[...,None,:], data.shape[2], axis=2)
                offset = np.choose(gain, (omap[...,0], omap[...,1], omap[...,2]))
                rg = np.choose(gain, (rmap[...,0], rmap[...,1], rmap[...,2]))
                noise = np.choose(gain, (nmap[...,0], nmap[...,1], nmap[...,2]))
                bpix = np.choose(gain, (bmap[...,0], bmap[...,1], bmap[...,2]))

                data -= offset
                data *= rg

                if correct_offset_drift:
                    lhd = np.mean(data[x//2-10:x//2,y//2-5:y//2+5,:], axis=(0,1))
                    data[:x//2, :, :] -= lhd
                    drift_lh.append(lhd)

                    uhd = np.mean(data[x//2:x//2+10,y//2-5:y//2+5,:], axis=(0,1))
                    data[x//2:, :, :] -= uhd
                    drift_uh.append(lhd)

                histCalOffsetCor.fill(data)


                ddset[...] = np.moveaxis(data, 2, 0)
                ddsetm[...] = np.moveaxis(bpix, 2, 0)
                ddsetg[...] = np.moveaxis(gain, 2, 0).astype(np.uint8)

                if mean_im is None:
                        mean_im = np.nanmean(data, axis=2)
                        single_im = data[...,0]

                if do_pattern_classification:

                    ddsetcm = ofile.create_dataset(h5path+"/pixels_cm",
                                         oshape,
                                         chunks=(chunk_size_idim, oshape[1], oshape[2]),
                                         dtype=np.float32)

                    ddsetc = ofile.create_dataset(h5path+"/pixels_classified",
                                         oshape,
                                         chunks=(chunk_size_idim, oshape[1], oshape[2]),
                                         dtype=np.float32, compression="gzip")

                    ddsetp = ofile.create_dataset(h5path+"/patterns",
                                                 oshape,
                                                 chunks=(chunk_size_idim, oshape[1], oshape[2]),
                                                 dtype=np.int32, compression="gzip")


                    patternClassifierLH._noisemap = noise[:x//2, :, :]
                    patternClassifierUH._noisemap = noise[x//2:, :, :]

                    data = cmCorrection.correct(data)  # correct for the row common mode
                    ddsetcm[...] = np.moveaxis(data, 2, 0)

                    dataLH = data[:x//2, :, :]
                    dataUH = data[x//2:, :, :]

                    dataLH, patternsLH = patternClassifierLH.classify(dataLH)
                    dataUH, patternsUH = patternClassifierUH.classify(dataUH)

                    data[:x//2, :, :] = dataLH
                    data[x//2:, :, :] = dataUH

                    patterns = np.zeros(data.shape, patternsLH.dtype)
                    patterns[:x//2, :, :] = patternsLH
                    patterns[x//2:, :, :] = patternsUH

                    data[data < split_evt_primary_threshold*noise] = 0
                    ddsetc[...] = np.moveaxis(data, 2, 0)
                    ddsetp[...] = np.moveaxis(patterns, 2, 0)

                    histCalPcorr.fill(data)
                    data[patterns != 100] = np.nan
                    histCalPcorrS.fill(data)

                    if mean_im_cc is None:
                        mean_im_cc = np.nanmean(data, axis=2)
                        single_im_cc = data[...,0]

        except Exception as e:
            print("Couldn't calibrate data in {}: {}".format(f, e))
 ```

 %% Cell type:code id: tags:

 ``` python
 if correct_offset_drift:
    lhds = np.concatenate(drift_lh)
    uhds = np.concatenate(drift_uh)
    fig = plt.figure(figsize=(10,5))
    ax = fig.add_subplot(111)
    ax.plot(lhds, label="Lower hem.")
    ax.plot(uhds, label="Upper hem.")
    ax.set_xlabel("Frame #")
    ax.set_xlabel("Offset drift (ADU)")
 ```

 %% Cell type:code id: tags:

 ``` python
 if do_pattern_classification:
    print("******************LOWER HEMISPHERE******************\n")

    patternStatsLH = patternClassifierLH.getPatternStats()
    fig = plt.figure(figsize=(15,15))
    ax = fig.add_subplot(4,4,1)
    sfields = ["singles", "first singles", "clusters"]
    mfields = ["doubles", "triples", "quads"]
    relativeOccurances =  []
    labels = []
    for i, f in enumerate(sfields):
        relativeOccurances.append(patternStatsLH[f])
        labels.append(f)
    for i, f in enumerate(mfields):
        for k in range(len(patternStatsLH[f])):
            relativeOccurances.append(patternStatsLH[f][k])
            labels.append(f+"("+str(k)+")")
    relativeOccurances = np.array(relativeOccurances, np.float)
    relativeOccurances/=np.sum(relativeOccurances)
    pie = ax.pie(relativeOccurances, labels=labels, autopct='%1.1f%%', shadow=True)
    ax.set_title("Pattern occurrence")
    # Set aspect ratio to be equal so that pie is drawn as a circle.
    a = ax.axis('equal')

    smaps = ["singlemap", "firstsinglemap", "clustermap"]
    for i, m in enumerate(smaps):

        ax = fig.add_subplot(4,4,2+i)

        pmap = ax.imshow(patternStatsLH[m], interpolation="nearest", vmax=2*np.nanmedian(patternStatsLH[m]))
        ax.set_title(m)
        cb = fig.colorbar(pmap)

    mmaps = ["doublemap", "triplemap", "quadmap"]
    k = 0
    for i, m in enumerate(mmaps):

        for j in range(4):
            ax = fig.add_subplot(4,4,2+len(smaps)+k)
            pmap = ax.imshow(patternStatsLH[m][j], interpolation="nearest", vmax=2*np.median(patternStatsLH[m][j]))
            ax.set_title(m+"("+str(j)+")")
            cb = fig.colorbar(pmap)
            k+=1
 ```

 %% Cell type:code id: tags:

 ``` python
 if do_pattern_classification:
    patternStatsUH = patternClassifierUH.getPatternStats()
    fig = plt.figure(figsize=(15,15))
    ax = fig.add_subplot(4,4,1)
    sfields = ["singles", "first singles", "clusters"]
    mfields = ["doubles", "triples", "quads"]
    relativeOccurances =  []
    labels = []
    for i, f in enumerate(sfields):
        relativeOccurances.append(patternStatsUH[f])
        labels.append(f)
    for i, f in enumerate(mfields):
        for k in range(len(patternStatsUH[f])):
            relativeOccurances.append(patternStatsUH[f][k])
            labels.append(f+"("+str(k)+")")
    relativeOccurances = np.array(relativeOccurances, np.float)
    relativeOccurances/=np.sum(relativeOccurances)
    pie = ax.pie(relativeOccurances, labels=labels, autopct='%1.1f%%', shadow=True)
    ax.set_title("Pattern occurrence")
    # Set aspect ratio to be equal so that pie is drawn as a circle.
    a = ax.axis('equal')

    smaps = ["singlemap", "firstsinglemap", "clustermap"]
    for i, m in enumerate(smaps):

        ax = fig.add_subplot(4,4,2+i)

        pmap = ax.imshow(patternStatsUH[m], interpolation="nearest", vmax=2*np.nanmedian(patternStatsUH[m]))
        ax.set_title(m)
        cb = fig.colorbar(pmap)

    mmaps = ["doublemap", "triplemap", "quadmap"]
    k = 0
    for i, m in enumerate(mmaps):

        for j in range(4):
            ax = fig.add_subplot(4,4,2+len(smaps)+k)
            pmap = ax.imshow(patternStatsUH[m][j], interpolation="nearest", vmax=np.median(patternStatsUH[m][j]))
            ax.set_title(m+"("+str(j)+")")
            cb = fig.colorbar(pmap)
            k+=1

    print("******************UPPER HEMISPHERE******************\n")
 ```

 %% Cell type:code id: tags:

 ``` python
 if do_pattern_classification:
    t0 = PrettyTable()
    t0.title = "Total number of Counts after all corrections"
    t0.field_names = ["Hemisphere","Singles", "First-singles", "Clusters"]
    t0.add_row(["LH", patternStatsLH['singles'], patternStatsLH['first singles'], patternStatsLH['clusters']])
    t0.add_row(["UH", patternStatsUH['singles'], patternStatsUH['first singles'], patternStatsUH['clusters']])

    print(t0)

    t1 = PrettyTable()

    t1.field_names = ["Index","D-LH", "D-UH", "T-LH", "T-UH", "Q-LH", "Q-UH"]

    t1.add_row([0, patternStatsLH['doubles'][0], patternStatsUH['doubles'][0], patternStatsLH['triples'][0], patternStatsUH['triples'][0], patternStatsLH['quads'][0], patternStatsUH['quads'][0]])
    t1.add_row([1, patternStatsLH['doubles'][1], patternStatsUH['doubles'][1], patternStatsLH['triples'][1], patternStatsUH['triples'][1], patternStatsLH['quads'][1], patternStatsUH['quads'][1]])
    t1.add_row([2, patternStatsLH['doubles'][2], patternStatsUH['doubles'][2], patternStatsLH['triples'][2], patternStatsUH['triples'][2], patternStatsLH['quads'][2], patternStatsUH['quads'][2]])
    t1.add_row([3, patternStatsLH['doubles'][3], patternStatsUH['doubles'][3], patternStatsLH['triples'][3], patternStatsUH['triples'][3], patternStatsLH['quads'][3], patternStatsUH['quads'][3]])

    print(t1)
 ```

 %% Cell type:code id: tags:

 ``` python
 if do_pattern_classification:
    doublesLH = patternStatsLH['doubles'][0] + patternStatsLH['doubles'][1] + patternStatsLH['doubles'][2] + patternStatsLH['doubles'][3]
    triplesLH = patternStatsLH['triples'][0] + patternStatsLH['triples'][1] + patternStatsLH['triples'][2] + patternStatsLH['triples'][3]
    quadsLH = patternStatsLH['quads'][0] + patternStatsLH['quads'][1] + patternStatsLH['quads'][2] + patternStatsLH['quads'][3]
    allsinglesLH = patternStatsLH['singles'] + patternStatsLH['first singles']
    eventsLH = allsinglesLH + doublesLH + triplesLH + quadsLH

    doublesUH = patternStatsUH['doubles'][0] + patternStatsUH['doubles'][1] + patternStatsUH['doubles'][2] + patternStatsUH['doubles'][3]
    triplesUH = patternStatsUH['triples'][0] + patternStatsUH['triples'][1] + patternStatsUH['triples'][2] + patternStatsUH['triples'][3]
    quadsUH = patternStatsUH['quads'][0] + patternStatsUH['quads'][1] + patternStatsUH['quads'][2] + patternStatsUH['quads'][3]
    allsinglesUH = patternStatsUH['singles'] + patternStatsUH['first singles']
    eventsUH = allsinglesUH + doublesUH + triplesUH + quadsUH

    reloccurLH = np.array([allsinglesLH/eventsLH, doublesLH/eventsLH, triplesLH/eventsLH, quadsLH/eventsLH])
    reloccurUH = np.array([allsinglesUH/eventsUH, doublesUH/eventsUH, triplesUH/eventsUH, quadsUH/eventsUH])
 ```

 %% Cell type:code id: tags:

 ``` python
 if do_pattern_classification:
    fig = plt.figure(figsize=(10,5))
    ax = fig.add_subplot(1,2,1)
    labels = ['singles', 'doubles', 'triples', 'quads']
    pie = ax.pie(reloccurLH, labels=labels, autopct='%1.1f%%', shadow=True)
    ax.set_title("Pattern occurrence LH")
    # Set aspect ratio to be equal so that pie is drawn as a circle.
    a = ax.axis('equal')
    ax = fig.add_subplot(1,2,2)
    pie = ax.pie(reloccurUH, labels=labels, autopct='%1.1f%%', shadow=True)
    ax.set_title("Pattern occurrence UH")
    # Set aspect ratio to be equal so that pie is drawn as a circle.
    a = ax.axis('equal')
 ```

 %% Cell type:code id: tags:

 ``` python
 ho,eo,co,so = histCalOffsetCor.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')
 ```

 %% Cell type:code id: tags:

 ``` python
 if do_pattern_classification:
    h1,e1L,c1L,s1L = histCalPcorr.get()
    h1s,e1Ls,c1Ls,s1Ls = histCalPcorrS.get()


    d = [
         {'x': c1L,
          'y': h1,
          'y_err': np.sqrt(h1[:]),
          'drawstyle': 'steps-mid',
          'label': 'Split event corrected'},
         {'x': c1Ls,
          'y': h1s,
          'y_err': np.sqrt(h1s[:]),
          'drawstyle': 'steps-mid',
          'label': 'Single pixel hits'}
    ]


    fig = xana.simplePlot(d, aspect=1, x_label='Energy(ADU)',
                          y_label='Number of occurrences', figsize='2col',
                          y_log=True, x_range=(0,200),x_log=False,
                          legend='top-center-frame-2col')
 ```

 %% Cell type:markdown id: tags:

 ## Mean Image of first Sequence ##

 %% Cell type:code id: tags:

 ``` python
 fig = xana.heatmapPlot(mean_im,
                       x_label='Columns', y_label='Rows',
                       lut_label='Signal (ADU)',
                       x_range=(0,y),
                       y_range=(0,x), vmin=-50, vmax=500)

 if do_pattern_classification:
    fig = xana.heatmapPlot(mean_im_cc,
                       x_label='Columns', y_label='Rows',
                       lut_label='Signal (ADU)',
                       x_range=(0,y),
                       y_range=(0,x), vmin=-50, vmax=500)
 ```

 %% Cell type:markdown id: tags:

 ## Single Shot of first Sequnce ##

 %% Cell type:code id: tags:

 ``` python
 fig = xana.heatmapPlot(single_im,
                       x_label='Columns', y_label='Rows',
                       lut_label='Signal (ADU)',
                       x_range=(0,y),
                       y_range=(0,x), vmin=-50, vmax=500)

 if do_pattern_classification:
    fig = xana.heatmapPlot(single_im_cc,
                       x_label='Columns', y_label='Rows',
                       lut_label='Signal (ADU)',
                       x_range=(0,y),
                       y_range=(0,x), vmin=-50, vmax=500)
 ```

 %% Cell type:code id: tags:

 ``` python
 ```