diff --git a/acid/fig_flux_dependence.py b/acid/fig_flux_dependence.py
index 82a1561..276ac7a 100644
--- a/acid/fig_flux_dependence.py
+++ b/acid/fig_flux_dependence.py
@@ -7,13 +7,14 @@
 import numpy as np
 from scipy.ndimage import label
 import matplotlib.pyplot as plt
+from matplotlib.gridspec import GridSpec
 import matplotlib
 import random
-
+from utils import load_avg, radial_integral
 from utils import get_potential, get_flux
 
 matplotlib.rcParams.update({'font.size': 8})
-
+plt.ion()
 
 PATH = '/data/visitors/nanomax/20200372/2021062308/raw/sample/%06u.h5'
 CENTER = -1500, 515//2
@@ -75,7 +76,52 @@ def get_hits_for_scans(scans, plot=True):
         hitlist.append(count_hits(arr, plot=plot, cutoff=cutoff))
     return hitlist
 
-if True:
+def integrate_scans(scans, plot=False):
+    integrals = []
+    for scan in scans:
+        print(scan)
+        with h5py.File(PATH%scan, 'r') as fp:
+            arr = fp['entry/measurement/merlin/frames'][:]
+            arr[:] = arr * mask
+        av_image = np.mean(arr, axis=0)
+        m1 = make_mask(arr[0], 181, 20)
+        m2 = make_mask(arr[0], 117, 30)
+        m3 = make_mask(arr[0], 150, 10)
+        if plot:
+            fig, ax = plt.subplots(ncols=4)
+            ax[0].imshow(m1)
+            ax[1].imshow(m2)
+            ax[2].imshow(m3)
+            ax[3].imshow(np.log10(av_image))
+        integrals.append([np.mean(m1 * av_image),
+                          np.mean(m2 * av_image),
+                          np.mean(m3 * av_image)])
+    return np.array(integrals)  # [scan, peak]
+
+### Curve integration
+if 0:
+    fluxes, curves = [], []
+    scans = np.arange(339, 339+15)  # first series, highest flux
+    for i in range(15): # 15 fluxes
+        print('***', i)
+        curves_ = []
+        for scan in (scans + i * 15):
+            print(scan)
+            pix, I = radial_integral(load_avg(scan))
+            curves_.append(I)
+        curves.append(curves_)
+        fluxes.append(get_flux(scans[0] + i * 15).mean())
+    pots = [get_potential(s) for s in scans]
+    curves = np.array(curves)  # indexed [flux, pot, peak]
+    np.savez('fig_flux_dependence_integration.npz', curves=curves, pots=pots, fluxes=fluxes)
+else:
+    dct = np.load('fig_flux_dependence_integration.npz')
+    fluxes = dct['fluxes']
+    pots = dct['pots']
+    curves = dct['curves']
+
+### Hit counting
+if 0:
     fluxes, hits = [], []
     scans = np.arange(339, 339+15)  # first series, highest flux
     for i in range(15): # 15 fluxes
@@ -84,25 +130,57 @@ def get_hits_for_scans(scans, plot=True):
         pots = [get_potential(s) for s in scans]
         fluxes.append(get_flux(scans[0] + i * 15).mean())
     hits = np.array(hits)  # indexed [flux, pot, peak]
-    np.savez('fig_flux_dependence.npz', hits=hits, pots=pots, fluxes=fluxes)
+    np.savez('fig_flux_dependence_hits.npz', hits=hits, pots=pots, fluxes=fluxes)
 else:
-    dct = np.load('fig_flux_dependence.npz')
+    dct = np.load('fig_flux_dependence_hits.npz')
     fluxes = dct['fluxes']
     pots = dct['pots']
     hits = dct['hits']
 
-plt.figure(figsize=(3.33,5))
-plt.subplots_adjust(bottom=.08, left=.16, right=.99, top=.99)
+## Figure
+fig = plt.figure(figsize=(3.33,6))
+gs = GridSpec(2, 1, figure=fig, height_ratios=[1, .25])
+ax = [fig.add_subplot(gs[0, 0]),
+      fig.add_subplot(gs[1, 0])]
+plt.subplots_adjust(bottom=.05, left=.14, right=.99, top=.99)
+
 for iflux in range(15):
     offset = -iflux * 50
-    plt.plot(pots, hits[iflux, :, 1] + offset, 'x-',)# lw=2-iflux/10)
+    ax[0].errorbar(pots, hits[iflux, :, 1] + offset,
+                 yerr=None,  # np.sqrt(hits[iflux, :, 1]),
+                 linestyle='-', marker='o', ms=2.5, lw=.8)
     flux = fluxes[iflux]
     order = int(np.log10(flux))
     prefactor = flux / (10**order)
-    plt.text(.6, offset+10, '%.1f $\\times 10^{%d}$' % (prefactor, order),
+    ax[0].text(.62, offset+10, '%.1f$\\cdot 10^{%d}$' % (prefactor, order),
              ha='right')
-plt.text(.6, 100, 'Flux (s$^{-1}$):', ha='right')
-plt.xlabel('E (V vs. Ag/AgCl sat.)')
-plt.ylabel('$\\beta$-PdH$_x$(111) intensity')
-plt.yticks([0,100,200])
-plt.savefig('fig_flux_dependence.pdf')
+ax[0].text(.62, 100, 'Flux (s$^{-1}$):', ha='right')
+ax[0].set_xlabel('E (V vs. Ag/AgCl sat.)')
+ax[0].set_ylabel('$\\beta$-PdH$_x$(111) diffraction hits', labelpad=-5)
+ax[0].set_yticks([0,100,200])
+
+
+###
+i0, i1 = 158, 183
+p0 = curves[:, :, i0:i1].sum(axis=-1) - (curves[:, :, i0] + curves[:, :, i1]) / 2 * (i1 - i0)
+i0, i1 = 95, 130
+p1 = curves[:, :, i0:i1].sum(axis=-1) - (curves[:, :, i0] + curves[:, :, i1]) / 2 * (i1 - i0)
+# p0 and p1 indexed [flux, pot]
+
+kw = dict(linestyle='-', marker='o', ms=2.5, lw=.8)
+i_pot = 10
+p0 = p0[:, i_pot] / fluxes * 1e10
+p1 = p1[:, i_pot] / fluxes * 1e10
+ax[1].plot(fluxes, p0, **kw, color='k')
+ax[1].plot(fluxes, p1, **kw, color='gray')
+ax[1].set_xscale('log')
+ax[1].set_xlabel('X-ray flux (s$^{-1}$)', labelpad=-5)
+ax[1].set_ylabel('$|I|$ / $I_0$ (arb.)')
+ax[1].text(8e8, p0[-1] + .4, 'Pd(111)')
+ax[1].text(8e8, p1[-1] + .4, '$\\beta$-PdH$_x$(111)', color='gray')
+ax[1].text(1.3e10, -.3, 'E = %.1f V' % pots[i_pot], va='bottom', ha='right')
+ax[1].set_ylim([-.5, 6.3])
+
+fig.text(.01, .98, 'a)')
+fig.text(.01, .23, 'b)')
+plt.savefig('fig_flux_dependence.pdf')
\ No newline at end of file