diff --git a/src/toolbox_scs/routines/boz.py b/src/toolbox_scs/routines/boz.py
index d280b309a12a37ab1d2bb089d674ec6f9e7acf4e..cd0d512df7b7b1cdcfd476d29f94773062c4f055 100644
--- a/src/toolbox_scs/routines/boz.py
+++ b/src/toolbox_scs/routines/boz.py
@@ -1459,7 +1459,7 @@ def ff_refine_crit_sk(p, alpha, params, arr_dark, arr, tid, rois,
     # drop saturated shots
     d = data.where(data['sat_sat'] == False, drop=True)
 
-    r = xas(d, 40, Iokey='np_mean', Itkey='0', nrjkey='0', fluorescence=True)
+    r = xas(d, 40, Iokey='np_mean_sk', Itkey='0', nrjkey='0', fluorescence=True)
 
     err_sigma = np.nansum(r['sigmaA'])
     err_mean = (1.0 - np.nanmean(r['muA']))**2
@@ -1650,7 +1650,7 @@ def nl_crit_sk(p, domain, alpha, arr_dark, arr, tid, rois, mask, flat_field,
     # drop saturated shots
     d = data.where(data['sat_sat'] == False, drop=True)
 
-    v = snr(d['np_mean'].values.flatten(), d['0'].values.flatten(),
+    v = snr(d['np_mean_sk'].values.flatten(), d['0'].values.flatten(),
             methods=['weighted'])
     err = 1e8*v['weighted']['s']**2
 
@@ -2014,7 +2014,7 @@ def inspect_correction_sk(params, ff, gain=None):
 
     scale = 1e-6
 
-    f, axs = plt.subplots(1, 3, figsize=(8, 2), sharex=True)
+    f, axs = plt.subplots(1, 3, figsize=(8, 2.5), sharex=True)
 
     # nbins = np.linspace(0.01, 1.0, 100)
 
@@ -2029,16 +2029,13 @@ def inspect_correction_sk(params, ff, gain=None):
         good_d = d.where(d['sat_sat'] == False, drop=True)
         sat_d = d.where(d['sat_sat'], drop=True)
 
-        good_d['np_mean'] = 0.5*(good_d['n']+good_d['p'])
-        sat_d['np_mean'] = 0.5*(sat_d['n']+sat_d['p'])
-
-        snr_v = snr(good_d['np_mean'].values.flatten(),
+        snr_v = snr(good_d['np_mean_sk'].values.flatten(),
                     good_d['0'].values.flatten(), verbose=True)
 
         m = snr_v['direct']['mu']
         h, xedges, yedges, img = axs[k].hist2d(
             g*scale*good_d['0'].values.flatten(),
-            good_d['np_mean'].values.flatten()/good_d['0'].values.flatten(),
+            good_d['np_mean_sk'].values.flatten()/good_d['0'].values.flatten(),
             [photon_scale, np.linspace(0.95, 1.05, 150)*m],
             cmap='Blues',
             norm=LogNorm(vmin=0.2, vmax=200),
@@ -2046,7 +2043,7 @@ def inspect_correction_sk(params, ff, gain=None):
             )
         h, xedges, yedges, img2 = axs[k].hist2d(
             g*scale*sat_d['0'].values.flatten(),
-            sat_d['np_mean'].values.flatten()/sat_d['0'].values.flatten(),
+            sat_d['np_mean_sk'].values.flatten()/sat_d['0'].values.flatten(),
             [photon_scale, np.linspace(0.95, 1.05, 150)*m],
             cmap='Reds',
             norm=LogNorm(vmin=0.2, vmax=200),
@@ -2288,7 +2285,7 @@ def process_module(arr, tid, dark, rois, mask=None, sat_level=511,
     
     # np_mean roi where we normalize the sum of flat_field
     np_mean  = (r['n'] + r['p'])/(ff['n'] + ff['p'])
-    v['np_mean'] = np_mean.sum(axis=(2,3))
+    v['np_mean_sk'] = np_mean.sum(axis=(2,3))
 
     res = xr.Dataset()
 
@@ -2298,9 +2295,9 @@ def process_module(arr, tid, dark, rois, mask=None, sat_level=511,
         res[n] = xr.DataArray(ensure_on_host(v[n]), coords=r_coords, dims=dims)
         res[n + '_sat'] = xr.DataArray(ensure_on_host(r_sat[n][:, :]),
                                        coords=r_coords, dims=dims)
-    res['np_mean'] = xr.DataArray(ensure_on_host(v['np_mean']),
+    res['np_mean_sk'] = xr.DataArray(ensure_on_host(v['np_mean_sk']),
                                   coords=r_coords, dims=dims)
-    res['np_mean_sat'] = res['n_sat'] + res['p_sat']
+    res['np_mean_sk_sat'] = res['n_sat'] + res['p_sat']
 
     for n in rois.keys():
         roi = rois[n]