03b_bear_necessities.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Stripped-down version of 03"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hi all, this is just a stripped-down version of notebook 3, which took about 3 minutes to execute. That sucks.\n",
    "And yes, I know it should be 'bare' . It's a joke, not a typo."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%run function_definitions.ipynb\n",
    "from scipy.optimize import curve_fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "bg_filenames = [\n",
    "    './PTB_DATA/2016_MAY_31/2016_MAY_31_2_5MeV_NE213_FG_12.mpa',\n",
    "    './PTB_DATA/2016_JUN_02/2016_JUN_02_2_5MeV_NE213_FG_1.mpa']\n",
    "bg_livetimes = [\n",
    "    36001.78,\n",
    "    40002.56\n",
    "    ]\n",
    "ng_filenames = [\n",
    "    './PTB_DATA/2016_MAY_31/2016_MAY_31_2_5MeV_NE213_FG_7.mpa', # 50 kV 2.5 mA 0°\n",
    "    './PTB_DATA/2016_JUN_01/2016_JUN_01_2_5MeV_NE213_FG_1.mpa', # 50 kV 2.5 mA 0° (SAME SETTINGS!)\n",
    "    ]\n",
    "ng_livetimes = [\n",
    "    9806.20,\n",
    "    11950.19,\n",
    "    ]\n",
    "\n",
    "# For the combination\n",
    "fg_livetime = ng_livetimes[0] + ng_livetimes[1]\n",
    "bg_livetime = sum(bg_livetimes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "e_scale = pickle.load(open('./08_ascii_files/e_scale.p', 'rb'))\n",
    "p_scale_to_e_and = [-14.325000000000003, 2.3875]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# NG data\n",
    "spec2d_list = []\n",
    "for ng_filename in ng_filenames:\n",
    "    line_n = find_line_start(ng_filename,'[CDAT2')\n",
    "    spec2d_raw = get_spectrum(ng_filename, line_n)\n",
    "    spec2d = np.reshape(spec2d_raw, (1024,4096))\n",
    "    spec2d_list.append(spec2d)\n",
    "#spec2d = sum(spec2d_list)\n",
    "\n",
    "# NOTE THIS IS A NEW LINE... Only analize the first run, not all combined!\n",
    "# NOTE modified AGAIN August 2016. spec2d = spec2d_list[0] -> spec2d = spec2d_list[0] + spec2d_list[1]\n",
    "spec2d = spec2d_list[0] + spec2d_list[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# BG data\n",
    "spec2d_bg_list = []\n",
    "for bg_filename in bg_filenames:\n",
    "    line_n = find_line_start(bg_filename,'[CDAT2')\n",
    "    spec2d_raw = get_spectrum(bg_filename, line_n)\n",
    "    spec2d_bg = np.reshape(spec2d_raw, (1024,4096))\n",
    "    spec2d_bg_list.append(spec2d_bg)\n",
    "# Sum the background\n",
    "spec2d_bg = sum(spec2d_bg_list)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "class Cut:\n",
    "    '''\n",
    "    Produce a cut function to cut on a percentile list.\n",
    "    '''\n",
    "    \n",
    "    pct = 0.\n",
    "    pct_list = [True for x in range(4096) ] \n",
    "    x_thr = 0.\n",
    "    y_thr = 0.\n",
    "    \n",
    "    def get_percentile(self, spec2d, pct):\n",
    "        self.pct_list = []\n",
    "        self.pct = pct\n",
    "        for x in range(len(spec2d[0])):\n",
    "            self.pct_list.append(self.get_percentile_pos(spec2d[:,x],pct))\n",
    "        \n",
    "    def get_percentile_pos(self, dist, percent):\n",
    "        tot = sum(dist)\n",
    "        count = 0\n",
    "        for i, x in enumerate(dist):\n",
    "            if count >= tot*percent/100.:\n",
    "                return i\n",
    "            count += x\n",
    "        # We should never arrive here\n",
    "        print('WARNING! Error in get_percentile_pos')\n",
    "        return -1\n",
    "    \n",
    "    def f_cut(self, x, y):\n",
    "        '''\n",
    "        FALSE = CUT\n",
    "        TRUE = KEEP\n",
    "        '''\n",
    "        cut_y = self.pct_list[x-1]\n",
    "        if y < cut_y:\n",
    "            return False\n",
    "        if x < self.x_thr:\n",
    "            return False\n",
    "        if y < self.y_thr:\n",
    "            return False\n",
    "        else:\n",
    "            return True\n",
    "        \n",
    "    def plot_bounds(self, **kwargs):\n",
    "        plt.plot(range(len(self.pct_list)), self.pct_list, **kwargs)\n",
    "        plt.axvline(self.x_thr, **kwargs)\n",
    "        plt.axhline(self.y_thr, **kwargs)\n",
    "        \n",
    "    def plot_percentile(self, **kwargs):\n",
    "        plt.plot(range(len(self.pct_list)), self.pct_list, **kwargs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "cut_list = []\n",
    "# 99 percent looks good, add to this list if you disagree\n",
    "for percentile in [99.0]:\n",
    "    cut = Cut()\n",
    "    cut.get_percentile(spec2d_bg, percentile)\n",
    "    cut_list.append(cut)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "cut = cut_list[0]\n",
    "cut.x_thr = 48\n",
    "print('Threshold is at %.1f keVee' % e_scale[cut.x_thr])\n",
    "spec2d_cut = cut_hist2d(spec2d, cut.f_cut)\n",
    "spec2d_bg_cut = cut_hist2d(spec2d_bg, cut.f_cut)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Settings equal for all runs\n",
    "cut = cut_list[0]\n",
    "cut.x_thr = 48\n",
    "\n",
    "# Now, for all runs:\n",
    "# - Take 2d spectrum\n",
    "# - Apply cuts based on BG percentiles\n",
    "# - Project on x axis to obtain 1d (energy) spectrum\n",
    "# - Divide by lifetime or run \n",
    "fg_disc   = 1/fg_livetime*np.sum((cut_hist2d(spec2d_list[0], cut.f_cut) + cut_hist2d(spec2d_list[1], \n",
    "                                                                                                 cut.f_cut)), axis=0)\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# Calculate errors:\n",
    "# - Assume only Poissonian errors\n",
    "# WARNING: - uncertainty on cut is NOT taken into account\n",
    "# WARNING: - zero error for zero count, screws all fits up.\n",
    "fg_disc_err   = 1/fg_livetime*np.sqrt(np.sum((cut_hist2d(spec2d_list[0], cut.f_cut) + cut_hist2d(spec2d_list[1], \n",
    "                                                                                                 cut.f_cut)), axis=0))\n",
    "\n",
    "\n",
    "# Also for BG, because this is how much of the BG leaks into signal!\n",
    "bg_disc = 1/bg_livetime*np.sum(cut_hist2d(spec2d_bg, cut.f_cut), axis=0)\n",
    "bg_disc_err = 1/bg_livetime*np.sqrt(np.sum(cut_hist2d(spec2d_bg, cut.f_cut), axis=0))\n",
    "\n",
    "# Subtract BG leakage\n",
    "fg_disc_sub = fg_disc - bg_disc\n",
    "fg_disc_sub_err = np.sqrt(fg_disc_err**2 + bg_disc_err**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def get_percentile_location(arr, pct):\n",
    "    '''\n",
    "    No idea why fucking numpy doesnt work on 1d arrays gvd\n",
    "    '''\n",
    "    total = sum(arr)\n",
    "    current_sum = 0.\n",
    "    i = 0\n",
    "    while current_sum < total*pct:\n",
    "        current_sum += arr[i]\n",
    "        i += 1\n",
    "    return i"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def rebin(arr, factor=2, mode='avg'):\n",
    "    '''\n",
    "    Rebin an array. Take integer factor. Modes available:\n",
    "      - avg: average. Default.\n",
    "      - sum: Sum.\n",
    "      - quad: adding all in quadrature. Mostly for errors.\n",
    "      - quad_avg: adding all in quadrature, then taking average. For errors on averages.\n",
    "    '''\n",
    "    if mode == 'avg':\n",
    "        rebin = np.array([1/factor*sum(arr[i:i + factor]) for i in np.arange(0, len(arr), factor)])\n",
    "    elif mode == 'sum':\n",
    "        rebin = np.array([sum(arr[i:i + factor]) for i in np.arange(0, len(arr), factor)])\n",
    "    elif mode == 'quad':\n",
    "        rebin = np.array([np.sqrt(sum(arr[i:i + factor]**2))\n",
    "                          for i in np.arange(0, len(arr), factor)])\n",
    "    elif mode == 'quad_avg':\n",
    "        rebin = np.array([1/factor*np.sqrt(sum(arr[i:i + factor]**2))\n",
    "                          for i in np.arange(0, len(arr), factor)])\n",
    "    else:\n",
    "        raise SyntaxError('Did not get the rebinning mode, you specified: %s' % mode)\n",
    "    \n",
    "    return rebin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "pct = 0.99\n",
    "cutoffs = []\n",
    "e_scale_he_bins = []\n",
    "fg_he = []\n",
    "bg_he = []\n",
    "\n",
    "for limit in np.arange(400, 4096, 50):\n",
    "    dist_bg = np.sum(spec2d_bg[:, limit:limit+50], axis=1)\n",
    "    current_cutoff = get_percentile_location(dist_bg, pct)\n",
    "    cutoffs.append(current_cutoff)\n",
    "    e_scale_he_bins.append(limit+25)\n",
    "    fg_he.append(np.sum(spec2d[current_cutoff:, limit:limit+50]))\n",
    "    bg_he.append(np.sum(spec2d_bg[current_cutoff:, limit:limit+50]))\n",
    "    \n",
    "\n",
    "e_scale_he = poly_1(np.array(e_scale_he_bins), *p_scale_to_e_and)\n",
    "\n",
    "# Make arrays to enable numpy magic!\n",
    "# NOTE that these are still NUMBERS and not rate!\n",
    "fg_he = np.array(fg_he)\n",
    "bg_he = np.array(bg_he)\n",
    "\n",
    "# Now calculate Poissonian errors\n",
    "fg_he_err = 1/(ng_livetimes[0] + ng_livetimes[1])*np.sqrt(fg_he)\n",
    "bg_he_err = 1/sum(bg_livetimes)*np.sqrt(bg_he)\n",
    "\n",
    "fg_he = 1/(ng_livetimes[0] + ng_livetimes[1])*fg_he\n",
    "bg_he = 1/sum(bg_livetimes)*bg_he\n",
    "\n",
    "fg_he_sub = fg_he - bg_he\n",
    "fg_he_sub_err = np.sqrt(bg_he_err**2 + fg_he_err**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Need for 'nicer histogram'  instead of imshow\n",
    "# We make event-by-event lists of values that we can histogram, instead of setting bin content.\n",
    "spec2d_counts_x = []\n",
    "spec2d_counts_y = []\n",
    "\n",
    "for i in range(len(spec2d)):\n",
    "    for j in range(len(spec2d[0])):\n",
    "        for k in range(spec2d[i,j]):\n",
    "            spec2d_counts_y.append(i)\n",
    "            spec2d_counts_x.append(j)\n",
    "spec2d_counts_x = np.array(spec2d_counts_x)\n",
    "spec2d_counts_y = np.array(spec2d_counts_y)\n",
    "\n",
    "# Same for bg...\n",
    "spec2d_bg_counts_x = []\n",
    "spec2d_bg_counts_y = []\n",
    "\n",
    "for i in range(len(spec2d_bg)):\n",
    "    for j in range(len(spec2d_bg[0])):\n",
    "        for k in range(spec2d_bg[i,j]):\n",
    "            spec2d_bg_counts_y.append(i)\n",
    "            spec2d_bg_counts_x.append(j)\n",
    "spec2d_bg_counts_x = np.array(spec2d_bg_counts_x)\n",
    "spec2d_bg_counts_y = np.array(spec2d_bg_counts_y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "led_region = np.array(e_scale_he_bins)[np.array(cutoffs) > 200]\n",
    "\n",
    "# For new values first just take the values we found previously\n",
    "fg_he_proper = fg_he.copy()\n",
    "fg_he_err_proper = fg_he_err.copy()\n",
    "bg_he_proper = bg_he.copy()\n",
    "bg_he_err_proper = bg_he_err.copy()\n",
    "fg_he_sub_proper = fg_he_sub.copy()\n",
    "fg_he_sub_err_proper = fg_he_sub_err.copy()\n",
    "\n",
    "for i, e in enumerate(e_scale_he_bins):\n",
    "    if e in led_region:\n",
    "        fg_he_proper[i] = np.sum(spec2d[170:190, e-25:e+25])\n",
    "        bg_he_proper[i] = np.sum(spec2d_bg[170:190, e-25:e+25])\n",
    "        # Now calculate Poissonian errors\n",
    "        fg_he_err_proper[i] = 1/(ng_livetimes[0] + ng_livetimes[1])*np.sqrt(fg_he_proper[i])\n",
    "        bg_he_err_proper[i] = 1/sum(bg_livetimes)*np.sqrt(bg_he_proper[i])\n",
    "\n",
    "        fg_he_proper[i] = 1/(ng_livetimes[0] + ng_livetimes[1])*fg_he_proper[i]\n",
    "        bg_he_proper[i] = 1/sum(bg_livetimes)*bg_he_proper[i]\n",
    "\n",
    "        fg_he_sub_proper[i] = fg_he_proper[i] - bg_he_proper[i]\n",
    "        fg_he_sub_err_proper[i] = np.sqrt(bg_he_err_proper[i]**2 + fg_he_err_proper[i]**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from IPython.utils import io\n",
    "with io.capture_output() as output_14:\n",
    "    %run '14_a_closer_look_at_acceptance.ipynb'\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Rebinning"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "e_scale_r = rebin(e_scale, mode='avg') # r for rebin\n",
    "fg_disc_sub_r = rebin(fg_disc_sub, mode='avg')\n",
    "fg_disc_sub_err_r = rebin(fg_disc_sub_err, mode='quad_avg')\n",
    "# Manual calculation with other method confirms the error bars\n",
    "ehm = np.sqrt(abs(rebin(fg_disc_sub, mode = 'sum') * fg_livetime)) / fg_livetime / 2\n",
    "\n",
    "\n",
    "def correct_for_acc(arr, acceptance_array, cutoff = 1e-6):\n",
    "    assert len(arr) == len(acceptance_array)\n",
    "    arr_c = []\n",
    "    for a, acc in zip(arr, acceptance_array):\n",
    "        if acc < cutoff:\n",
    "            arr_c.append(0)\n",
    "        else:\n",
    "            arr_c.append(a/acc)\n",
    "    return np.array(arr_c)\n",
    "\n",
    "fg_disc_sub_r_a = correct_for_acc(fg_disc_sub_r, acc_e_r)\n",
    "fg_disc_sub_err_r_a = correct_for_acc(fg_disc_sub_err_r, acc_e_r)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import time\n",
    "time.clock()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.4.5"
  },
  "latex_envs": {
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 0
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}