Initial commit

AGIPD/correct_agipd_batch.py

+
+# coding: utf-8
+
+# In[29]:
+
+import sys
+#in_folder, out_folder, base_store, offset_store, mem_cells, sequences
+in_folder = sys.argv[1] # "/gpfs/exfel/exp/SPB/201701/p002012/raw/r0100"
+out_folder = sys.argv[2]# "./corrected_test"
+base_store = sys.argv[3]
+offset_store = sys.argv[4]
+sequences = sys.argv[6] #[0]
+mem_cells = int(sys.argv[5]) # 30
+max_cells = mem_cells
+overwrite = True if sys.argv[7] == "True" else False
+if sequences.upper()  != "all":
+    sequences = [int(s) for s in  sequences.split(",")]
+else:
+    sequences = None
+
+# make sure a cluster is running with ipcluster start --n=32, give it a while to start
+import os
+import h5py
+import numpy as np
+import matplotlib.pyplot as plt
+from ipyparallel import Client
+view = Client()[:]
+view.use_dill()
+gains = np.arange(3)
+cells = np.arange(max_cells)
+
+
+QUADRANTS = 4
+MODULES_PER_QUAD = 4
+DET_FILE_INSET = "AGIPD"
+max_cells = mem_cells
+
+if in_folder[-1] == "/":
+    in_folder = in_folder[:-1]
+out_folder = "{}/{}".format(out_folder, os.path.split(in_folder)[-1])
+print("Outputting to {}".format(out_folder))
+
+if not os.path.exists(out_folder):
+    os.makedirs(out_folder)
+elif not overwrite:
+    raise AttributeError("Output path exists! Exiting")
+
+# In[2]:
+
+def combine_stack(d, sdim):
+    combined = np.zeros((sdim, 2048,2048))
+    combined[...] = np.nan
+    
+    dy = 0
+    for i in range(16):
+        
+        
+        if i < 8:
+            dx = -512
+            mx = 1
+            my = i % 8
+            combined[:, my*128+dy:(my+1)*128+dy,
+                     mx*512-dx:(mx+1)*512-dx] = np.rollaxis(d[i],2,1)[:,:,::-1]
+            dy += 30
+            if i == 3:
+                dy += 30
+        elif i < 12:
+            dx = 100
+            if i == 8:
+                dy = 4*30 + 30 +50
+            
+            mx = 1
+            my = i % 8 +4
+            combined[:, my*128+dy:(my+1)*128+dy,
+                     mx*512-dx:(mx+1)*512-dx] = np.rollaxis(d[i],2,1)[:,::-1,:]
+            dy += 30
+        else:
+            dx = 100
+            if i == 11:
+                dy = 50
+            
+            mx = 1
+            my = i - 14
+            
+            combined[:, my*128+dy:(my+1)*128+dy,
+                     mx*512-dx:(mx+1)*512-dx] = np.rollaxis(d[i],2,1)[:,::-1,:]
+            dy += 30
+            
+            
+        
+    return combined
+
+
+# In[3]:
+
+rel_gains = []
+offsets = []
+base_offsets = []
+thresholds = []
+masks = []
+saveFile = h5py.File(base_store, "r")
+for i in range(16):
+    qm = "Q{}M{}".format(i//4+1, i%4+1)
+    data = np.array(saveFile["{}/RelativeGain/0/data".format(qm)])
+    mask = np.array(saveFile["{}/RelativeGain/0/mask".format(qm)])
+    rel_gains.append(data)
+    data = np.array(saveFile["{}/BaseOffset/0/data".format(qm)])
+    mask = np.maximum(np.array(saveFile["{}/BaseOffset/0/mask".format(qm)]), np.max(mask, axis=3)[...,None])
+    base_offsets.append(data)
+    data = np.array(saveFile["{}/Threshold/0/data".format(qm)])
+    mask = np.maximum(np.array(saveFile["{}/BaseOffset/0/mask".format(qm)]), mask)
+    masks.append(mask)
+    thresholds.append(data)
+saveFile.close()
+
+
+saveFile = h5py.File(offset_store, "r")
+for i in range(16):
+    qm = "Q{}M{}".format(i//4+1, i%4+1)
+    data = np.array(saveFile["{}/Offset/0/data".format(qm)])
+    offsets.append(data)
+saveFile.close()
+
+
+# In[4]:
+
+# set everything up filewise
+from queue import Queue
+#if not os.path.exists(out_folder):
+#    os.makedirs(out_folder)
+#elif not overwrite:
+#    raise AttributeError("Output path exists! Exiting")
+ 
+def map_modules_from_files(filelist):
+    module_files = {}
+    mod_ids = {}
+    for quadrant in range(0, QUADRANTS):
+        for module in range(0, MODULES_PER_QUAD):
+            name = "Q{}M{}".format(quadrant + 1, module + 1)
+            module_files[name] = Queue()
+            num = quadrant * 4 + module
+            mod_ids[name] = num
+            file_infix = "{}{:02d}".format(DET_FILE_INSET, num)
+            for file in filelist:
+                if file_infix in file:
+                    module_files[name].put(file)
+    return module_files, mod_ids
+
+dirlist = os.listdir(in_folder)
+file_list = []
+for entry in dirlist:
+    #only h5 file
+    abs_entry = "{}/{}".format(in_folder, entry)
+    if os.path.isfile(abs_entry) and os.path.splitext(abs_entry)[1] == ".h5":
+        
+        if sequences is None:
+            file_list.append(abs_entry)
+        else:
+            for seq in sequences:
+                if "{:05d}.h5".format(seq) in abs_entry:
+                    file_list.append(os.path.abspath(abs_entry))
+                    
+mapped_files, mod_ids = map_modules_from_files(file_list)
+
+
+# In[5]:
+
+import copy
+from functools import partial
+def correct_module(cells, inp):
+    import numpy as np
+    import copy
+    import h5py
+    if True:
+    
+        filename, filename_out, channel, offset, base_offset, rel_gain, threshold, mask = inp
+        print(filename, filename_out)
+
+        infile = h5py.File(filename, "r", driver="core")
+        im = np.array(infile["/INSTRUMENT/SPB_DET_AGIPD1M-1/DET/{}CH0:xtdf/image/data".format(channel)])
+
+        cellid = np.array(infile["/INSTRUMENT/SPB_DET_AGIPD1M-1/DET/{}CH0:xtdf/image/cellId".format(channel)][::2])
+        pulses = np.array(infile["/INSTRUMENT/SPB_DET_AGIPD1M-1/DET/{}CH0:xtdf/image/pulseId".format(channel)][::2])
+        trains = np.array(infile["/INSTRUMENT/SPB_DET_AGIPD1M-1/DET/{}CH0:xtdf/image/trainId".format(channel)][::2])
+        statii = np.array(infile["/INSTRUMENT/SPB_DET_AGIPD1M-1/DET/{}CH0:xtdf/image/status".format(channel)][::2])
+
+
+        dont_copy = ["data", "pulseId", "cellId", "trainId", "status"]
+        dont_copy = ["INSTRUMENT/SPB_DET_AGIPD1M-1/DET/{}CH0:xtdf/image/{}".format(channel, do)
+                    for do in dont_copy]
+
+        outfile = h5py.File(filename_out, "w", driver="core")
+        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)
+        outfile.flush()
+
+        outfile["/INSTRUMENT/SPB_DET_AGIPD1M-1/DET/{}CH0:xtdf/image/cellId".format(channel)] = cellid
+        outfile["/INSTRUMENT/SPB_DET_AGIPD1M-1/DET/{}CH0:xtdf/image/pulseId".format(channel)] = pulses
+        outfile["/INSTRUMENT/SPB_DET_AGIPD1M-1/DET/{}CH0:xtdf/image/trainId".format(channel)] = trains
+        outfile["/INSTRUMENT/SPB_DET_AGIPD1M-1/DET/{}CH0:xtdf/image/status".format(channel)] = statii
+
+        outfile.flush()
+        infile.close()
+
+        ga = im[1::2, 0, ...]
+        im = im[0::2, 0, ...].astype(np.float32)
+
+
+        im = np.rollaxis(im, 2)
+        im = np.rollaxis(im, 2, 1)
+
+        ga = np.rollaxis(ga, 2)
+        ga = np.rollaxis(ga, 2, 1)
+
+        tmap = threshold
+
+        gain = np.zeros(ga.shape, np.uint8)
+
+        for cc in range(im.shape[2]//cells):
+            tg = gain[...,cc*cells:(cc+1)*cells]
+            tga = ga[...,cc*cells:(cc+1)*cells]
+            tg[...] = 255
+            tg[tga < tmap[..., 0]] = 0
+            tg[tga > tmap[..., 1]] = 2
+            tg[tg == 255] = 1
+            gain[...,cc*cells:(cc+1)*cells] = tg
+
+        om = offset
+        rc = rel_gain
+
+        offsetb = om
+
+        do = np.median(base_offset-om, axis=(0,1))
+        offsetb[...,1:] -= do[...,1:]
+        msk = np.zeros(list(im.shape)+[2], np.uint8)
+        for cc in range(im.shape[2]//cells):
+            tg = gain[...,cc*cells:(cc+1)*cells]
+            offset = np.choose(tg, (offsetb[...,0], offsetb[...,1], offsetb[...,2]))
+            rel_cor = np.choose(tg, (rc[...,0], rc[...,1], rc[...,2]))
+            tim = im[...,cc*cells:(cc+1)*cells]
+            tim = tim - offset
+
+            tim /= rel_cor
+            im[...,cc*cells:(cc+1)*cells] = tim
+            msk[...,cc*cells:(cc+1)*cells,:] = mask
+
+        outfile["INSTRUMENT/SPB_DET_AGIPD1M-1/DET/{}CH0:xtdf/image/data".format(channel)] = np.rollaxis(np.rollaxis(im,1), 2)
+        outfile["INSTRUMENT/SPB_DET_AGIPD1M-1/DET/{}CH0:xtdf/image/gain".format(channel)] = np.rollaxis(np.rollaxis(gain,1), 2)
+        outfile["INSTRUMENT/SPB_DET_AGIPD1M-1/DET/{}CH0:xtdf/image/mask".format(channel)] = np.rollaxis(np.rollaxis(msk,1), 2)
+        outfile.close()
+    
+
+done = False
+first_files = []
+while not done:
+    inp = []
+    dones = []
+    first = True
+    for i in range(16):
+        qm = "Q{}M{}".format(i//4 +1, i % 4 + 1)    
+        if qm in mapped_files and not mapped_files[qm].empty():
+            fname_in = mapped_files[qm].get()            
+            dones.append(mapped_files[qm].empty())
+        else:
+            first_files.append((None, None))
+            continue
+        fout = os.path.abspath("{}/{}".format(out_folder, (os.path.split(fname_in)[-1]).replace("RAW", "CORR")))
+        if first:
+            first_files.append((fname_in, fout))
+        inp.append((fname_in, fout, i, offsets[i][...,:max_cells,:].astype(np.float32), base_offsets[i][...,:max_cells,:].astype(np.float32),
+                    rel_gains[i][...,:max_cells,:], thresholds[i][...,:max_cells,:], masks[i][...,:max_cells,:]))
+    first = False
+    p = partial(correct_module, max_cells)
+    
+    r = view.map_sync(p, inp)
+    done = all(dones)
+
+
+
+
+# In[6]:
+
+corrected = []
+raw = []
+gains = []
+for i, ff in enumerate(first_files):
+    try:
+        rf, cf = ff
+        if rf is None:
+            raise Exception("File not present")
+        infile = h5py.File(rf, "r")
+        raw.append(np.array(infile["/INSTRUMENT/SPB_DET_AGIPD1M-1/DET/{}CH0:xtdf/image/data".format(i)][:max_cells:2,0,...]))
+        infile.close()
+        
+        infile = h5py.File(cf, "r")
+        corrected.append(np.array(infile["/INSTRUMENT/SPB_DET_AGIPD1M-1/DET/{}CH0:xtdf/image/data".format(i)][:max_cells,...]))
+        gains.append(np.array(infile["/INSTRUMENT/SPB_DET_AGIPD1M-1/DET/{}CH0:xtdf/image/gain".format(i)][:max_cells,...]))
+        infile.close()
+        
+    except Exception as e:
+        print(e)
+        corrected.append(np.zeros((max_cells, 512, 128)))
+
+
+# In[16]:
+
+combined = combine_stack(corrected, corrected[0].shape[0])
+combined_raw = combine_stack(raw, raw[0].shape[0])
+combined_g = combine_stack(gains, gains[0].shape[0])
+
+
+# In[24]:
+
+fig = plt.figure(figsize=(10,5))
+ax = fig.add_subplot(111)
+im = ax.imshow(np.mean(combined_raw[...][:1300,400:1600],axis=0), vmin=0, vmax=8000, cmap="jet")
+fig.colorbar(im, ax=ax)
+fig.savefig("{}/mean_first_train_RAW.png".format(out_folder))
+
+
+fig = plt.figure(figsize=(10,5))
+ax = fig.add_subplot(111)
+im = ax.imshow(np.mean(combined[...][:1300,400:1600], axis=0), vmin=0, vmax=1000, cmap="jet")
+fig.colorbar(im, ax=ax)
+fig.savefig("{}/mean_first_train_CORR.png".format(out_folder))
+
+fig = plt.figure(figsize=(10,5))
+ax = fig.add_subplot(111)
+im = ax.imshow(np.mean(combined_g[...][:1300,400:1600], axis=0), vmin=0, vmax=3, cmap="jet")
+fig.colorbar(im, ax=ax)
+fig.savefig("{}/mean_first_train_GAIN.png".format(out_folder))
+
+
+# In[26]:
+
+fig = plt.figure(figsize=(10,5))
+ax = fig.add_subplot(111)
+im = ax.imshow(np.std(combined_raw[...][:1300,400:1600],axis=0), vmin=0, vmax=8000, cmap="jet")
+plt.colorbar(im, ax=ax)
+fig.savefig("{}/std_first_train_RAW.png".format(out_folder))
+
+
+fig = plt.figure(figsize=(10,5))
+ax = fig.add_subplot(111)
+im = ax.imshow(np.std(combined[...][:1300,400:1600], axis=0), vmin=0, vmax=20000, cmap="jet")
+plt.colorbar(im, ax=ax)
+fig.savefig("{}/std_first_train_CORR.png".format(out_folder))
+
+fig = plt.figure(figsize=(10,5))
+ax = fig.add_subplot(111)
+im = ax.imshow(np.std(combined_g[...][:1300,400:1600], axis=0), vmin=0, vmax=3, cmap="jet")
+plt.colorbar(im, ax=ax)
+fig.savefig("{}/std_first_train_GAIN.png".format(out_folder))
+
+
+# In[25]:
+
+fig = plt.figure(figsize=(10,5))
+ax = fig.add_subplot(111)
+im = ax.imshow(np.max(combined_raw[...][:1300,400:1600],axis=0), vmin=0, vmax=8000, cmap="jet")
+plt.colorbar(im)
+fig.savefig("{}/max_first_train_RAW.png".format(out_folder))
+
+
+fig = plt.figure(figsize=(10,5))
+ax = fig.add_subplot(111)
+im = ax.imshow(np.max(combined[...][:1300,400:1600], axis=0), vmin=0, vmax=20000, cmap="jet")
+plt.colorbar(im)
+fig.savefig("{}/max_first_train_CORR.png".format(out_folder))
+
+fig = plt.figure(figsize=(10,5))
+ax = fig.add_subplot(111)
+im = ax.imshow(np.max(combined_g[...][:1300,400:1600], axis=0), vmin=0, vmax=3, cmap="jet")
+plt.colorbar(im)
+fig.savefig("{}/max_first_train_GAIN.png".format(out_folder))
+
+
+

README.rst

+Python Scripted Calibration
+===========================
+
+First: do not run this on the Maxell gateway. Rather, `salloc`
+a node for yourself first::
+
+   salloc -p exfel/upex -t 01:00:00
+
+where `-p` gives the partition to use: exfel or upex and `-t`
+the duration the node should be allocated. Then `ssh` onto 
+that node.
+
+(optionally) Set up the environment::
+
+   source setup_env.sh
+
+On Maxwell execution rights might be stripped. If you get
+an error for the above command, first run::
+
+   chmod +x setup_env.sh
+
+If running headless, be sure to set `MPLBACKEND=Agg`, via::
+
+   export MPLBACKEND=Agg
+
+This is automatically done if you source the environment script.
+
+
+Run the script::
+
+    python calibrate.py --input /gpfs/exfel/exp/SPB/201701/p002012/raw/r0100 \
+       --output ../../test_out --mem-cells 30 --detector AGIPD --sequences 0,1   
+
+Here `--input` should point to a directory of `RAW` files for the detector you
+are calibrating. They will be output into the folder specified by `--output`, 
+which will have the run number or the last folder in the hiearchy of the input
+appended. Additionally, you need to specify the number of `--mem-cells` used
+for the run, as well as the `--detector`. Finally, you can optionally 
+specify to only process certain `--sequences` of files, matching the sequence
+numbers of the `RAW` input. These should be given as a comma-separated list.
+
+You'll get a series of plots in the output directory as well.
+
+
+
+
+
+
+
+

calibrate.py

+import argparse
+from subprocess import Popen
+parser = argparse.ArgumentParser(description="Main entry point "
+                                             "for offline calibration")
+parser.add_argument("--input", type=str)
+parser.add_argument("--output", type=str)
+parser.add_argument("--base-cal-store", default="base", type=str)
+parser.add_argument("--offset-cal-store", default="base", type=str)
+parser.add_argument("--mem-cells", type=int)
+parser.add_argument("--detector", type=str)
+parser.add_argument("--sequences", type=str, default="All")
+parser.add_argument("--overwrite", action="store_true")
+
+args = vars(parser.parse_args())
+
+in_folder = args["input"]
+out_folder = args["output"]
+base_store = args["base_cal_store"]
+offset_store = args["offset_cal_store"]
+mem_cells = args["mem_cells"]
+detector = args["detector"]
+sequences = args["sequences"]
+overwrite = str(args["overwrite"])
+if detector.upper() not in ["AGIPD"]:
+    raise AttributeError("Unknown detector: {}".format(detector))
+
+if base_store == "base":
+    if detector.upper() == "AGIPD":
+        base_store = "/gpfs/exfel/data/user/haufs/scripted_offline_cal/calstore/agipd_base_cal.h5"
+if offset_store == "base":
+    if detector.upper() == "AGIPD":
+        offset_store = "/gpfs/exfel/data/user/haufs/scripted_offline_cal/calstore/agipd_offset_cal.h5"
+
+if detector.upper() == "AGIPD":
+    script = "AGIPD/correct_agipd_batch.py"
+cmd = ["python", script, in_folder, out_folder, base_store, offset_store, str(mem_cells), sequences, overwrite]
+Popen(cmd).wait()

setup_env.sh

+#/bin/bash
+export MPLBACKEND=Agg
+source /gpfs/exfel/data/user/haufs/karabo/activate
+ipcluster start --n=16 --daemon