Multi-spot 5-Samples analysis - Leakage coefficient fit.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Leakage coefficient fit\n",
    "\n",
    "> *This notebook estracts the leakage coefficient from multi-spot smFRET measurements of 4 dsDNA samples.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import software"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from fretbursts import *\n",
    "sns = init_notebook()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import seaborn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import pandas as pd\n",
    "from IPython.display import display\n",
    "from IPython.display import display, Math"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "data_dir = './data/multispot/'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "data_dir = os.path.abspath(data_dir) + '/'\n",
    "assert os.path.exists(data_dir), \"Path '%s' does not exist.\" % data_dir"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from glob import glob\n",
    "file_list = sorted(glob(data_dir + '*.hdf5'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "labels = ['7d', '12d', '17d', '22d', '27d', 'DO']\n",
    "files_dict = {lab: fname for lab, fname in zip(sorted(labels), file_list)}\n",
    "files_dict"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot style"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "PLOT_DIR = './figure/'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import matplotlib as mpl\n",
    "from cycler import cycler\n",
    "\n",
    "bmap = sns.color_palette(\"Set1\", 9)\n",
    "colors = np.array(bmap)[(1,0,2,3,4,8,6,7), :]\n",
    "colors_labels = ['blue', 'red', 'green', 'violet', 'orange', 'gray', 'brown', 'pink', ]\n",
    "for c, cl in zip(colors, colors_labels):\n",
    "    locals()[cl] = tuple(c) # assign variables with color names\n",
    "\n",
    "mpl.rcParams['axes.prop_cycle'] = cycler('color', colors)\n",
    "sns.palplot(colors)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysis Parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "## Background fit\n",
    "bg_kwargs_auto = dict(fun=bg.exp_fit,\n",
    "                 time_s = 30,\n",
    "                 tail_min_us = 'auto',\n",
    "                 F_bg=1.7,\n",
    "                 )\n",
    "\n",
    "## Burst search\n",
    "F = 6\n",
    "ph_sel = Ph_sel(Dex='Dem') \n",
    "#ph_sel = Ph_sel('all') \n",
    "size_min = 80\n",
    "\n",
    "## D-only peak fit with KDE\n",
    "bandwidth = 0.03\n",
    "binwidth = 0.025\n",
    "E_range_do = (-0.05, 0.1)\n",
    "weights = 'size'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "df = pd.DataFrame(index=['7d', '12d', '17d', 'DO'], columns=range(8), dtype=float)\n",
    "df.index.name = 'Sample'\n",
    "df.columns.name = 'Channel'\n",
    "E_do = df.copy()\n",
    "E_do_g = df.copy()\n",
    "nbursts = df.copy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Utility functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def print_fit_report(E_pr, gamma=1, leakage=0, dir_ex_t=0, math=True):\n",
    "    \"\"\"Print fit and standard deviation for both corrected and uncorrected E\n",
    "    Returns d.E_fit.\n",
    "    \"\"\"\n",
    "    E_corr = fretmath.correct_E_gamma_leak_dir(E_pr, gamma=gamma, leakage=leakage, dir_ex_t=dir_ex_t)\n",
    "    \n",
    "    E_pr_mean = E_pr.mean()*100\n",
    "    E_pr_delta = (E_pr.max() - E_pr.min())*100\n",
    "    \n",
    "    E_corr_mean = E_corr.mean()*100\n",
    "    E_corr_delta = (E_corr.max() - E_corr.min())*100\n",
    "    if math:\n",
    "        display(Math(r'\\text{Pre}\\;\\gamma\\quad\\langle{E}_{fit}\\rangle = %.1f\\%% \\qquad'\n",
    "                     '\\Delta E_{fit} = %.2f \\%%' % \\\n",
    "                     (E_pr_mean, E_pr_delta)))\n",
    "        display(Math(r'\\text{Post}\\;\\gamma\\quad\\langle{E}_{fit}\\rangle = %.1f\\%% \\qquad'\n",
    "                     '\\Delta E_{fit} = %.2f \\%%' % \\\n",
    "                     (E_corr_mean, E_corr_delta)))\n",
    "    else:\n",
    "        print('Pre-gamma  E (delta, mean):  %.2f  %.2f' % (E_pr_mean, E_pr_delta))\n",
    "        print('Post-gamma E (delta, mean):  %.2f  %.2f' % (E_corr_mean, E_corr_delta))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 7bp sample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "str(ph_sel)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "data_id = '7d'\n",
    "d7 = loader.photon_hdf5(files_dict[data_id])\n",
    "d7.calc_bg(**bg_kwargs_auto)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "d7.burst_search(m=10, F=F, ph_sel=ph_sel)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ds7 = d7.select_bursts(select_bursts.nd, th1=size_min)\n",
    "dx = ds7"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## KDE Fit\n",
    "E_do.loc[data_id] = bext.fit_bursts_kde_peak(dx, bandwidth=bandwidth, x_range=E_range_do,  \n",
    "                                             weights=weights)\n",
    "\n",
    "## Gaussian fit\n",
    "dx.E_fitter.histogram(binwidth=binwidth, weights=weights)\n",
    "dx.E_fitter.fit_histogram(mfit.factory_gaussian())\n",
    "E_do_g.loc[data_id] = dx.E_fitter.params['center']\n",
    "\n",
    "## D-only selection\n",
    "do_s = dx.select_bursts(select_bursts.E, E2=0.1)\n",
    "nbursts.loc[data_id] = do_s.num_bursts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dplot(dx, hist_fret, binwidth=binwidth, weights=weights, show_kde=True, show_kde_peak=True, show_fit_value=True);\n",
    "plt.xlim(xmax = 0.49)\n",
    "print_fit_report(E_do.loc[data_id])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dplot(dx, hist_fret, binwidth=binwidth, weights=weights, show_model=True, show_fit_value=True, fit_from='center');\n",
    "plt.xlim(xmax = 0.49)\n",
    "print_fit_report(E_do_g.loc[data_id])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Alternative plots"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(4, 2, figsize=(12, 8), sharex=True, sharey=True)\n",
    "fig.subplots_adjust(left=0.08, right=0.96, top=0.93, bottom=0.07,\n",
    "                    wspace=0.06, hspace=0.18)\n",
    "\n",
    "for ich, ax in enumerate(axes.ravel()):\n",
    "    mfit.plot_mfit(dx.E_fitter, ich=ich, ax=ax, plot_model=False, plot_kde=True)\n",
    "plt.xlim(-0.2, 0.3)\n",
    "print_fit_report(E_do.loc[data_id])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(4, 2, figsize=(12, 8), sharex=True, sharey=True)\n",
    "fig.subplots_adjust(left=0.08, right=0.96, top=0.93, bottom=0.07,\n",
    "                    wspace=0.06, hspace=0.18)\n",
    "\n",
    "for ich, ax in enumerate(axes.ravel()):\n",
    "    mfit.plot_mfit(dx.E_fitter, ich=ich, ax=ax)\n",
    "plt.xlim(-0.2, 0.3)\n",
    "print_fit_report(E_do_g.loc[data_id])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 12bp sample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "data_id = '12d'\n",
    "d12 = loader.photon_hdf5(files_dict[data_id])\n",
    "d12.calc_bg(**bg_kwargs_auto)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "d12.burst_search(m=10, F=F, ph_sel=ph_sel)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ds12 = d12.select_bursts(select_bursts.nd, th1=size_min)\n",
    "dx = ds12"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "## KDE Fit\n",
    "E_do.loc[data_id] = bext.fit_bursts_kde_peak(dx, bandwidth=bandwidth, x_range=E_range_do,  \n",
    "                                             weights=weights)\n",
    "\n",
    "## Gaussian fit\n",
    "dx.E_fitter.histogram(binwidth=binwidth, weights=weights)\n",
    "dx.E_fitter.fit_histogram(mfit.factory_gaussian())\n",
    "E_do_g.loc[data_id] = dx.E_fitter.params['center']\n",
    "\n",
    "## D-only selection\n",
    "do_s = dx.select_bursts(select_bursts.E, E2=0.1)\n",
    "nbursts.loc[data_id] = do_s.num_bursts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dplot(dx, hist_fret, binwidth=binwidth, weights=weights, show_kde=True, show_kde_peak=True, show_fit_value=True);\n",
    "plt.xlim(xmax = 0.49)\n",
    "print_fit_report(E_do.loc[data_id])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dplot(dx, hist_fret, binwidth=binwidth, weights=weights, show_model=True, show_fit_value=True, fit_from='center');\n",
    "plt.xlim(xmax = 0.49)\n",
    "print_fit_report(E_do_g.loc[data_id])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 17bp sample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "data_id = '17d'\n",
    "d17 = loader.photon_hdf5(files_dict[data_id])\n",
    "d17.calc_bg(**bg_kwargs_auto)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "d17.burst_search(m=10, F=F, ph_sel=ph_sel)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ds17 = d17.select_bursts(select_bursts.nd, th1=size_min)\n",
    "dx = ds17"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "## KDE Fit\n",
    "E_do.loc[data_id] = bext.fit_bursts_kde_peak(dx, bandwidth=bandwidth, x_range=E_range_do,  \n",
    "                                             weights=weights)\n",
    "\n",
    "## Gaussian fit\n",
    "dx.E_fitter.histogram(binwidth=binwidth, weights=weights)\n",
    "dx.E_fitter.fit_histogram(mfit.factory_two_gaussians(p1_center=0.03, p2_center=0.25))\n",
    "E_do_g.loc[data_id] = dx.E_fitter.params['p1_center']\n",
    "\n",
    "## D-only selection\n",
    "do_s = Sel(dx, select_bursts.E, E2=0.1)\n",
    "nbursts.loc[data_id] = do_s.num_bursts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dplot(dx, hist_fret, binwidth=binwidth, weights=weights, show_kde=True, show_kde_peak=True, show_fit_value=True);\n",
    "plt.xlim(xmax = 0.49)\n",
    "print_fit_report(E_do.loc[data_id])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dplot(dx, hist_fret, binwidth=binwidth, weights=weights, show_model=True, show_fit_value=True, fit_from='p1_center');\n",
    "plt.xlim(xmax = 0.49)\n",
    "print_fit_report(E_do_g.loc[data_id])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# DO sample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "data_id = 'DO'\n",
    "do = loader.photon_hdf5(files_dict[data_id])\n",
    "do.calc_bg(**bg_kwargs_auto)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "do.burst_search(L=10, m=10, F=F, ph_sel=ph_sel)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dos = do.select_bursts(select_bursts.nd, th1=size_min)\n",
    "dx = dos"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "## KDE Fit\n",
    "E_do.loc[data_id] = bext.fit_bursts_kde_peak(dx, bandwidth=bandwidth, x_range=E_range_do,  \n",
    "                                             weights=weights)\n",
    "\n",
    "## Gaussian fit\n",
    "dx.E_fitter.histogram(binwidth=binwidth, weights=weights)\n",
    "dx.E_fitter.fit_histogram(mfit.factory_gaussian())\n",
    "E_do_g.loc[data_id] = dx.E_fitter.params['center']\n",
    "\n",
    "## D-only selection\n",
    "nbursts.loc[data_id] = dx.num_bursts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dplot(dx, hist_fret, binwidth=binwidth, weights=weights, show_kde=True, show_kde_peak=True, show_fit_value=True);\n",
    "plt.xlim(xmax = 0.49)\n",
    "print_fit_report(E_do.loc[data_id])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dplot(dx, hist_fret, binwidth=binwidth, weights=weights, show_model=True, show_fit_value=True, fit_from='center');\n",
    "plt.xlim(xmax = 0.49)\n",
    "print_fit_report(E_do_g.loc[data_id])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ph_sel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "nbursts = nbursts.astype(int)\n",
    "nbursts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "E_do_kde = E_do.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "leakage_kde = (E_do_kde / (1 - E_do_kde)).round(6)\n",
    "leakage_kde"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "leakage_gauss = (E_do_g / (1 - E_do_g)).round(6)\n",
    "leakage_gauss"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "sns.set(style='ticks', font_scale=1.4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "colors = sns.color_palette('Paired', 8)\n",
    "mpl.rcParams['axes.prop_cycle'] = cycler('color', colors)\n",
    "sns.palplot(colors)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax2 = ax.twinx()\n",
    "\n",
    "kws = dict(lw=2, marker='o', ms=8)\n",
    "for i, did in enumerate(('7d', '12d', '17d', 'DO')):\n",
    "    (100*leakage_kde).loc[did].plot(label='%s KDE' % did, ax=ax, color=colors[1+i*2], **kws)\n",
    "    nbursts.loc[did].plot(ax=ax2, ls='--', lw=2.5, color=colors[1+i*2])\n",
    "\n",
    "for i, did in enumerate(('7d', '12d', '17d', 'DO')):    \n",
    "    (100*leakage_gauss).loc[did].plot(label='%s Gauss' % did, ax=ax, color=colors[i*2], **kws)\n",
    "    \n",
    "handles, lab = ax.get_legend_handles_labels()\n",
    "h = handles#[1::2] + handles[::2]\n",
    "l = lab[1::2] + lab[::2]\n",
    "ax.legend(ncol=2, loc=1, bbox_to_anchor=(1, 0.5), borderaxespad=0.)\n",
    "ax.set_ylim(0)\n",
    "\n",
    "ax2.set_ylim(0, 3200)\n",
    "plt.xlim(-0.25, 7.25)\n",
    "plt.xlabel('Channel')\n",
    "ax.set_ylabel('Leakage %')\n",
    "ax2.set_ylabel('# Bursts')\n",
    "sns.despine(offset=10, trim=True, right=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Save results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "leakage_kde.to_csv('results/Multi-spot - leakage coefficient all values KDE %s.csv' % str(ph_sel))\n",
    "leakage_gauss.to_csv('results/Multi-spot - leakage coefficient all values gauss %s.csv' % str(ph_sel))\n",
    "nbursts.to_csv('results/Multi-spot - leakage coefficient all values nbursts %s.csv' % str(ph_sel))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Average leakage"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mean per sample:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "lk_s = pd.DataFrame(index=['mean', 'std'], columns=E_do.index)\n",
    "\n",
    "lk_s.loc['mean'] = leakage_kde.mean(1)*100\n",
    "lk_s.loc['std'] = leakage_kde.std(1)*100\n",
    "\n",
    "lk_s['mean'] = lk_s.mean(1)\n",
    "lk_s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "lk_s.T[['std']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "lk_s.T[['mean']].plot(table=True)\n",
    "plt.gca().xaxis.set_visible(False)   # Hide Ticks"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mean per sample (weighted on the number of bursts):\n",
    "\n",
    "Number of bursts in D-only population:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "nbursts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "lk_sw = pd.DataFrame(index=['mean', 'std'], columns=E_do.index)\n",
    "\n",
    "lk_sw.loc['mean'] = (nbursts*leakage_kde).sum(1)/nbursts.sum(1)*100\n",
    "lk_sw.loc['std'] = (nbursts*leakage_kde).std(1)/nbursts.sum(1)*100\n",
    "\n",
    "lk_sw['mean'] = (nbursts.sum(1)*lk_sw).sum(1)/nbursts.sum(1).sum()\n",
    "lk_sw"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "lk_sw.loc['mean'].plot()\n",
    "ylim(2, 4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mean per channel:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "lk_c = pd.DataFrame(index=['mean', 'std'], columns=E_do.columns)\n",
    "\n",
    "lk_c.loc['mean'] = leakage_kde.mean()*100\n",
    "lk_c.loc['std'] = leakage_kde.std()*100\n",
    "\n",
    "lk_c['mean'] = lk_c.mean(1)\n",
    "lk_c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "lk_c.loc['mean'].plot()\n",
    "ylim(2, 4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mean per channel (weighted on the number of bursts):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "lk_cw = pd.DataFrame(index=['mean', 'std'], columns=E_do.columns)\n",
    "\n",
    "lk_cw.loc['mean'] = (nbursts*leakage_kde).sum()/nbursts.sum()*100\n",
    "lk_cw.loc['std'] = (nbursts*leakage_kde).std()/nbursts.sum()*100\n",
    "\n",
    "lk_cw['mean'] = lk_cw.mean(1)\n",
    "lk_cw"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "lk_cw.loc['mean'].plot()\n",
    "ylim(2, 4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "mch_plot_bg(d7)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> **NOTE:** There is a per-channel trend that cannot be ascribed to the background \n",
    "> because we performend a D-emission burst search and selection and the leakage vs ch\n",
    "> does not resemble the D-background vs channel curve.\n",
    ">\n",
    "> The effect is probably due to slight PDE variations (detectors + optics) that \n",
    "> slightly change $\\gamma$ on a per-spot basis."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Weighted mean of the weighted mean"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "lk = (lk_cw.ix['mean', :-1]*nbursts.mean()).sum()/(nbursts.mean().sum())/100\n",
    "lk_2 = (lk_sw.ix['mean', :-1]*nbursts.mean(1)).sum()/(nbursts.mean(1).sum())/100\n",
    "\n",
    "assert np.allclose(lk, lk_2)\n",
    "\n",
    "print('Mean leakage: %.6f' % lk)"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [default]",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}