Fix for the BOZ analysis

• in inspecting the flat field, use np.roll to avoid issues with rois of different size
• make find roi a bit more robust by dividing the used sensor in quarters within which the rois boundaries will be found
• fitting non-linearity of UP2589 did not work. The issue could be with not ideal dark, affecting certain data point at low intensity. Here we use weighted mean and standard deviation of both p/0 and n/0 orders for the fit instead of deviation from a transmission of 1.

src/toolbox_scs/routines/boz.py

@@ -342,11 +342,11 @@ def find_rois(data_mean, threshold):
     pY = data_mean.sum(axis=1)
 
     # along X
-    lowX = int(np.argmax(pX[:128] > threshold))  # 1st occurrence returned
-    highX = int(np.argmax(pX[128:] < threshold) + 128)  # 1st occ. returned
-    midX = (lowX + highX)//2
-    leftX = int(np.argmin(pX[(lowX+20):midX]) + lowX + 20)
-    rightX = int(np.argmin(pX[midX:highX-20]) + midX)
+    lowX = int(np.argmax(pX[:64] > threshold))  # 1st occurrence returned
+    highX = int(np.argmax(pX[192:] <= threshold) + 192)  # 1st occ. returned
+
+    leftX = int(np.argmin(pX[64:128]) + 64)
+    rightX = int(np.argmin(pX[128:192]) + 128)
 
     # along Y
     lowY = int(np.argmax(pY[:64] > threshold))  # 1st occurrence returned
@@ -795,26 +795,24 @@ def inspect_plane_fitting(avg, rois, vmin=None, vmax=None):
     fig, axs = plt.subplots(2, 3, sharex=True, figsize=(6, 6))
 
     img_rois = {}
+    centers = {}
+
+    for k, r in enumerate(['n', '0', 'p']):
+        roi = rois[r]
+        centers[r] = np.array([(roi['yl'] + roi['yh'])//2,
+                      (roi['xl'] + roi['xh'])//2])
 
+    d = '0'
+    roi = rois[d]
     for k, r in enumerate(['n', '0', 'p']):
-        img_rois[r] = avg[rois[r]['yl']:rois[r]['yh'],
-                          rois[r]['xl']:rois[r]['xh']]
+        img_rois[r] = np.roll(avg, tuple(centers[r] - centers[d]))[
+        roi['yl']:roi['yh'], roi['xl']:roi['xh']]
         im = axs[0, k].imshow(img_rois[r],
                               vmin=vmin,
                               vmax=vmax)
 
     for k, r in enumerate(['n', '0', 'p']):
-        if img_rois[r].shape[1] != img_rois['0'].shape[1]:
-            if k == 0:
-                n1 = img_rois[r].shape[1]
-                n = img_rois['0'].shape[1]
-                v = img_rois[r][:, (n1-n):]/img_rois['0']
-            else:
-                n1 = img_rois[r].shape[1]
-                n = img_rois['0'].shape[1]
-                v = img_rois[r][:, :-(n1-n)]/img_rois['0']
-        else:
-            v = img_rois[r]/img_rois['0']
+        v = img_rois[r]/img_rois['0']
         im2 = axs[1, k].imshow(v, vmin=0.2, vmax=1.1, cmap='RdBu_r')
 
     cbar = fig.colorbar(im, ax=axs[0, :], orientation="horizontal")
@@ -1073,16 +1071,17 @@ def nl_crit(p, domain, alpha, arr_dark, arr, tid, rois, mask, flat_field,
     # drop saturated shots
     d = data.where(data['sat_sat'] == False, drop=True)
 
-    # calculated error from transmission of 1.0
-    v = d['n'].values.flatten()/d['0'].values.flatten()
-    err1 = 1e8*np.nanmean((v - 1.0)**2)
+    v_1 = snr(d['n'].values.flatten(), d['0'].values.flatten(),
+            methods=['weighted'])
+    err_1 = 1e8*v_1['weighted']['s']**2
 
-    err2 = np.sum((Fmodel-np.arange(2**9))**2)
+    v_2 = snr(d['p'].values.flatten(), d['0'].values.flatten(),
+            methods=['weighted'])
+    err_2 = 1e8*v_2['weighted']['s']**2
 
-    # print(f'{err}: {p}')
-    # logging.info(f'{err}: {p}')
+    err_a = np.sum((Fmodel-np.arange(2**9))**2)
 
-    return (1.0 - alpha)*err1 + alpha*err2
+    return (1.0 - alpha)*0.5*(err_1 + err_2) + alpha*err_a
 
 
 def nl_fit(params, domain):
@@ -1175,44 +1174,68 @@ def inspect_nl_fit(res_fit):
     return f
 
 
-def snr(sig, ref, verbose=False):
-    """ Compute mean, std and SNR with and without weight from transmitted signal sig
-        and I0 signal ref
+def snr(sig, ref, methods=None, verbose=False):
+    """ Compute mean, std and SNR from transmitted signal sig and I0 signal ref.
+
+    Inputs
+    ------
+    sig: 1D signal samples
+    ref: 1D reference samples
+    methods: None by default or list of strings to select which methods to use.
+        Possible values are 'direct', 'weighted', 'diff'. In case of None, all
+        methods will be calculated.
+    verbose: booleand, if True prints calculated values
+
+    Returns
+    -------
+    dictionnary of [methods][value] where value is 'mu' for mean and 's' for
+    standard deviation.
+
     """
+    if methods is None:
+        methods = ['direct', 'weighted', 'diff']
+
     w = ref
     x = sig/ref
 
     mask = np.isfinite(x) & np.isfinite(sig) & np.isfinite(ref)
-    
+
     w = w[mask]
     sig = sig[mask]
     ref = ref[mask]
     x = x[mask]
 
-    # direct mean and std
-    mu = np.mean(x)
-    s = np.std(x)
-    if verbose:
-        print(f'mu: {mu}, s: {s}, snr: {mu/s}')
     res = {}
-    res['direct'] = {'mu': mu, 's':s}
+
+    # direct mean and std
+    if 'direct' in methods:
+        mu = np.mean(x)
+        s = np.std(x)
+        if verbose:
+            print(f'mu: {mu}, s: {s}, snr: {mu/s}')
+
+        res['direct'] = {'mu': mu, 's':s}
 
     # weighted mean and std
-    wmu = np.sum(sig)/np.sum(ref)
-    v1 = np.sum(w)
-    v2 = np.sum(w**2)
-    ws = np.sqrt(np.sum(w*(x - wmu)**2)/(v1 - v2/v1))
+    if 'weighted' in methods:
+        wmu = np.sum(sig)/np.sum(ref)
+        v1 = np.sum(w)
+        v2 = np.sum(w**2)
+        ws = np.sqrt(np.sum(w*(x - wmu)**2)/(v1 - v2/v1))
+
+        if verbose:
+            print(f'weighted mu: {wmu}, s: {ws}, snr: {wmu/ws}')
 
-    if verbose:
-        print(f'weighted mu: {wmu}, s: {ws}, snr: {wmu/ws}')
-    res['weighted'] = {'mu': wmu, 's':ws}
+        res['weighted'] = {'mu': wmu, 's':ws}
 
     # noise from diff
-    dmu = np.mean(x)
-    ds = np.std(np.diff(x))/np.sqrt(2)
-    if verbose:
-        print(f'diff mu: {dmu}, s: {ds}, snr: {dmu/ds}')
-    res['diff'] = {'mu': dmu, 's':ds}
+    if 'diff' in methods:
+        dmu = np.mean(x)
+        ds = np.std(np.diff(x))/np.sqrt(2)
+        if verbose:
+            print(f'diff mu: {dmu}, s: {ds}, snr: {dmu/ds}')
+
+        res['diff'] = {'mu': dmu, 's':ds}
 
     return res
 
@@ -1256,7 +1279,7 @@ def inspect_correction(params, gain=None):
 
     scale = 1e-6
 
-    f, axs = plt.subplots(3, 3, figsize=(6, 6), sharex=True, sharey=True)
+    f, axs = plt.subplots(3, 3, figsize=(8, 6), sharex=True)
 
     # nbins = np.linspace(0.01, 1.0, 100)
 
@@ -1278,16 +1301,11 @@ def inspect_correction(params, gain=None):
             snr_v = snr(good_d[n].values.flatten(),
                         good_d[r].values.flatten(), verbose=True)
 
-            if k == 0:
-                m = np.nanmean(good_d[n].values.flatten()
-                    /good_d[r].values.flatten())
-            else:
-                m = 1
-
+            m = snr_v['direct']['mu']
             h, xedges, yedges, img = axs[l, k].hist2d(
                 g*scale*good_d[r].values.flatten(),
-                good_d[n].values.flatten()/good_d[r].values.flatten()/m,
-                [photon_scale, np.linspace(0.95, 1.05, 150)],
+                good_d[n].values.flatten()/good_d[r].values.flatten(),
+                [photon_scale, np.linspace(0.95, 1.05, 150)*m],
                 cmap='Blues',
                 vmax=200,
                 norm=LogNorm(),
@@ -1295,13 +1313,14 @@ def inspect_correction(params, gain=None):
                 )
             h, xedges, yedges, img2 = axs[l, k].hist2d(
                 g*scale*sat_d[r].values.flatten(),
-                sat_d[n].values.flatten()/sat_d[r].values.flatten()/m,
-                [photon_scale, np.linspace(0.95, 1.05, 150)],
+                sat_d[n].values.flatten()/sat_d[r].values.flatten(),
+                [photon_scale, np.linspace(0.95, 1.05, 150)*m],
                 cmap='Reds',
                 vmax=200,
                 norm=LogNorm(),
                 # alpha=0.5 # make  the plot looks ugly with lots of white lines
                 )
+
             v = snr_v['direct']['mu']/snr_v['direct']['s']
             axs[l, k].text(0.4, 0.15, f'SNR: {v:.0f}',
                             transform = axs[l, k].transAxes)
@@ -1312,9 +1331,11 @@ def inspect_correction(params, gain=None):
             # axs[l, k].plot(3*nbins, 1+np.sqrt(2/(1e6*nbins)), c='C1', ls='--')
             # axs[l, k].plot(3*nbins, 1-np.sqrt(2/(1e6*nbins)), c='C1', ls='--')
 
+            axs[l, k].set_ylim([0.95*m, 1.05*m])
+
     for k in range(3):
-        for l in range(3):
-            axs[l, k].set_ylim([0.95, 1.05])
+        #for l in range(3):
+        #    axs[l, k].set_ylim([0.95, 1.05])
         if gain:
             axs[2, k].set_xlabel('#ph (10$^6$)')
         else: