Get Planck 2015 Parameter Fisher.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Get Planck MCMC parameter Covariance Matrices and Process into Fisher matrices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "59\n",
      "57\n",
      "82\n",
      "line  0 :  omegabh2\t\\Omega_b h^2\n",
      "        :  omegabh2\t\\Omega_b h^2\n",
      "        :  omegabh2\t\\Omega_b h^2\n",
      "\n",
      "line  1 :  omegach2\t\\Omega_c h^2\n",
      "        :  omegach2\t\\Omega_c h^2\n",
      "        :  omegach2\t\\Omega_c h^2\n",
      "\n",
      "line  2 :  theta\t100\\theta_{MC}\n",
      "        :  theta\t100\\theta_{MC}\n",
      "        :  theta\t100\\theta_{MC}\n",
      "\n",
      "line  3 :  tau\t\\tau\n",
      "        :  tau\t\\tau\n",
      "        :  tau\t\\tau\n",
      "\n",
      "line  4 :  logA\t{\\rm{ln}}(10^{10} A_s)\n",
      "        :  logA\t{\\rm{ln}}(10^{10} A_s)\n",
      "        :  logA\t{\\rm{ln}}(10^{10} A_s)\n",
      "\n",
      "line  5 :  ns\tn_s\n",
      "        :  ns\tn_s\n",
      "        :  ns\tn_s\n",
      "\n",
      "line  6 :  calPlanck\ty_{\\rm cal}\n",
      "        :  galfEE100\tA^{{\\rm dust}EE}_{100}\n",
      "        :  galfEE100\tA^{{\\rm dust}EE}_{100}\n",
      "\n",
      "line  7 :  galfEE100\tA^{{\\rm dust}EE}_{100}\n",
      "        :  galfEE100143\tA^{{\\rm dust}EE}_{100\\times 143}\n",
      "        :  galfEE100143\tA^{{\\rm dust}EE}_{100\\times 143}\n",
      "\n",
      "line  8 :  galfEE100143\tA^{{\\rm dust}EE}_{100\\times 143}\n",
      "        :  galfEE100217\tA^{{\\rm dust}EE}_{100\\times 217}\n",
      "        :  galfEE100217\tA^{{\\rm dust}EE}_{100\\times 217}\n",
      "\n",
      "line  9 :  galfEE100217\tA^{{\\rm dust}EE}_{100\\times 217}\n",
      "        :  galfEE143\tA^{{\\rm dust}EE}_{143}\n",
      "        :  galfEE143\tA^{{\\rm dust}EE}_{143}\n",
      "\n",
      "line  10 :  galfEE143\tA^{{\\rm dust}EE}_{143}\n",
      "        :  galfEE143217\tA^{{\\rm dust}EE}_{143\\times 217}\n",
      "        :  galfEE143217\tA^{{\\rm dust}EE}_{143\\times 217}\n",
      "\n",
      "line  11 :  galfEE143217\tA^{{\\rm dust}EE}_{143\\times 217}\n",
      "        :  galfEE217\tA^{{\\rm dust}EE}_{217}\n",
      "        :  galfEE217\tA^{{\\rm dust}EE}_{217}\n",
      "\n",
      "line  12 :  galfEE217\tA^{{\\rm dust}EE}_{217}\n",
      "        :  calPlanck\ty_{\\rm cal}\n",
      "        :  calPlanck\ty_{\\rm cal}\n",
      "\n",
      "line  13 :  H0*\tH_0\n",
      "        :  H0*\tH_0\n",
      "        :  H0*\tH_0\n",
      "\n",
      "line  14 :  omegal*\t\\Omega_\\Lambda\n",
      "        :  omegal*\t\\Omega_\\Lambda\n",
      "        :  omegal*\t\\Omega_\\Lambda\n",
      "\n",
      "line  15 :  omegam*\t\\Omega_m\n",
      "        :  omegam*\t\\Omega_m\n",
      "        :  omegam*\t\\Omega_m\n",
      "\n",
      "line  16 :  omegamh2*\t\\Omega_m h^2\n",
      "        :  omegamh2*\t\\Omega_m h^2\n",
      "        :  omegamh2*\t\\Omega_m h^2\n",
      "\n",
      "line  17 :  omeganuh2*\t\\Omega_\\nu h^2\n",
      "        :  omeganuh2*\t\\Omega_\\nu h^2\n",
      "        :  omeganuh2*\t\\Omega_\\nu h^2\n",
      "\n",
      "line  18 :  omegamh3*\t\\Omega_m h^3\n",
      "        :  omegamh3*\t\\Omega_m h^3\n",
      "        :  omegamh3*\t\\Omega_m h^3\n",
      "\n",
      "line  19 :  sigma8*\t\\sigma_8\n",
      "        :  sigma8*\t\\sigma_8\n",
      "        :  sigma8*\t\\sigma_8\n",
      "\n",
      "line  20 :  s8omegamp5*\t\\sigma_8 \\Omega_m^{0.5}\n",
      "        :  s8omegamp5*\t\\sigma_8 \\Omega_m^{0.5}\n",
      "        :  s8omegamp5*\t\\sigma_8 \\Omega_m^{0.5}\n",
      "\n",
      "line  21 :  s8omegamp25*\t\\sigma_8 \\Omega_m^{0.25}\n",
      "        :  s8omegamp25*\t\\sigma_8 \\Omega_m^{0.25}\n",
      "        :  s8omegamp25*\t\\sigma_8 \\Omega_m^{0.25}\n",
      "\n",
      "line  22 :  s8h5*\t\\sigma_8/h^{0.5}\n",
      "        :  s8h5*\t\\sigma_8/h^{0.5}\n",
      "        :  s8h5*\t\\sigma_8/h^{0.5}\n",
      "\n",
      "line  23 :  rmsdeflect*\t\\langle d^2\\rangle^{1/2}\n",
      "        :  rmsdeflect*\t\\langle d^2\\rangle^{1/2}\n",
      "        :  rmsdeflect*\t\\langle d^2\\rangle^{1/2}\n",
      "\n",
      "line  24 :  zrei*\tz_{\\rm re}\n",
      "        :  zrei*\tz_{\\rm re}\n",
      "        :  zrei*\tz_{\\rm re}\n",
      "\n",
      "line  25 :  A*\t10^9 A_s\n",
      "        :  A*\t10^9 A_s\n",
      "        :  A*\t10^9 A_s\n",
      "\n",
      "line  26 :  clamp*\t10^9 A_s e^{-2\\tau}\n",
      "        :  clamp*\t10^9 A_s e^{-2\\tau}\n",
      "        :  clamp*\t10^9 A_s e^{-2\\tau}\n",
      "\n",
      "line  27 :  DL40*\tD_{40}\n",
      "        :  DL40*\tD_{40}\n",
      "        :  DL40*\tD_{40}\n",
      "\n",
      "line  28 :  DL220*\tD_{220}\n",
      "        :  DL220*\tD_{220}\n",
      "        :  DL220*\tD_{220}\n",
      "\n",
      "line  29 :  DL810*\tD_{810}\n",
      "        :  DL810*\tD_{810}\n",
      "        :  DL810*\tD_{810}\n",
      "\n",
      "line  30 :  DL1420*\tD_{1420}\n",
      "        :  DL1420*\tD_{1420}\n",
      "        :  DL1420*\tD_{1420}\n",
      "\n",
      "line  31 :  DL2000*\tD_{2000}\n",
      "        :  DL2000*\tD_{2000}\n",
      "        :  DL2000*\tD_{2000}\n",
      "\n",
      "line  32 :  ns02*\tn_{s,0.002}\n",
      "        :  ns02*\tn_{s,0.002}\n",
      "        :  ns02*\tn_{s,0.002}\n",
      "\n",
      "line  33 :  yheused*\tY_P\n",
      "        :  yheused*\tY_P\n",
      "        :  yheused*\tY_P\n",
      "\n",
      "line  34 :  YpBBN*\tY_P^{\\rm{BBN}}\n",
      "        :  YpBBN*\tY_P^{\\rm{BBN}}\n",
      "        :  YpBBN*\tY_P^{\\rm{BBN}}\n",
      "\n",
      "line  35 :  DHBBN*\t10^5D/H\n",
      "        :  DHBBN*\t10^5D/H\n",
      "        :  DHBBN*\t10^5D/H\n",
      "\n",
      "line  36 :  age*\t{\\rm{Age}}/{\\rm{Gyr}}\n",
      "        :  age*\t{\\rm{Age}}/{\\rm{Gyr}}\n",
      "        :  age*\t{\\rm{Age}}/{\\rm{Gyr}}\n",
      "\n",
      "line  37 :  zstar*\tz_*\n",
      "        :  zstar*\tz_*\n",
      "        :  zstar*\tz_*\n",
      "\n",
      "line  38 :  rstar*\tr_*\n",
      "        :  rstar*\tr_*\n",
      "        :  rstar*\tr_*\n",
      "\n",
      "line  39 :  thetastar*\t100\\theta_*\n",
      "        :  thetastar*\t100\\theta_*\n",
      "        :  thetastar*\t100\\theta_*\n",
      "\n",
      "line  40 :  DAstar*\tD_{\\rm{A}}/{\\rm{Gpc}}\n",
      "        :  DAstar*\tD_{\\rm{A}}/{\\rm{Gpc}}\n",
      "        :  DAstar*\tD_{\\rm{A}}/{\\rm{Gpc}}\n",
      "\n",
      "line  41 :  zdrag*\tz_{\\rm{drag}}\n",
      "        :  zdrag*\tz_{\\rm{drag}}\n",
      "        :  zdrag*\tz_{\\rm{drag}}\n",
      "\n",
      "line  42 :  rdrag*\tr_{\\rm{drag}}\n",
      "        :  rdrag*\tr_{\\rm{drag}}\n",
      "        :  rdrag*\tr_{\\rm{drag}}\n",
      "\n",
      "line  43 :  kd*\tk_{\\rm D}\n",
      "        :  kd*\tk_{\\rm D}\n",
      "        :  kd*\tk_{\\rm D}\n",
      "\n",
      "line  44 :  thetad*\t100\\theta_{\\rm{D}}\n",
      "        :  thetad*\t100\\theta_{\\rm{D}}\n",
      "        :  thetad*\t100\\theta_{\\rm{D}}\n",
      "\n",
      "line  45 :  zeq*\tz_{\\rm{eq}}\n",
      "        :  zeq*\tz_{\\rm{eq}}\n",
      "        :  zeq*\tz_{\\rm{eq}}\n",
      "\n",
      "line  46 :  keq*\tk_{\\rm{eq}}\n",
      "        :  keq*\tk_{\\rm{eq}}\n",
      "        :  keq*\tk_{\\rm{eq}}\n",
      "\n",
      "line  47 :  thetaeq*\t100\\theta_{\\rm{eq}}\n",
      "        :  thetaeq*\t100\\theta_{\\rm{eq}}\n",
      "        :  thetaeq*\t100\\theta_{\\rm{eq}}\n",
      "\n",
      "line  48 :  thetarseq*\t100\\theta_{\\rm{s,eq}}\n",
      "        :  thetarseq*\t100\\theta_{\\rm{s,eq}}\n",
      "        :  thetarseq*\t100\\theta_{\\rm{s,eq}}\n",
      "\n",
      "line  49 :  rsDv057*\tr_{\\rm{drag}}/D_V(0.57)\n",
      "        :  rsDv057*\tr_{\\rm{drag}}/D_V(0.57)\n",
      "        :  rsDv057*\tr_{\\rm{drag}}/D_V(0.57)\n",
      "\n",
      "line  50 :  Hubble057*\tH(0.57)\n",
      "        :  Hubble057*\tH(0.57)\n",
      "        :  Hubble057*\tH(0.57)\n",
      "\n",
      "line  51 :  DA057*\tD_A(0.57)\n",
      "        :  DA057*\tD_A(0.57)\n",
      "        :  DA057*\tD_A(0.57)\n",
      "\n",
      "line  52 :  FAP057*\tF_{\\rm AP}(0.57)\n",
      "        :  FAP057*\tF_{\\rm AP}(0.57)\n",
      "        :  FAP057*\tF_{\\rm AP}(0.57)\n",
      "\n",
      "line  53 :  fsigma8z057*\tf\\sigma_8(0.57)\n",
      "        :  fsigma8z057*\tf\\sigma_8(0.57)\n",
      "        :  fsigma8z057*\tf\\sigma_8(0.57)\n",
      "\n",
      "line  54 :  sigma8z057*\t\\sigma_8(0.57)\n",
      "        :  sigma8z057*\t\\sigma_8(0.57)\n",
      "        :  sigma8z057*\t\\sigma_8(0.57)\n",
      "\n",
      "line  55 :  chi2_lowTEB*\t\\chi^2_{\\rm lowTEB}\n",
      "        :  chi2_plikEE*\t\\chi^2_{\\rm plikEE}\n",
      "        :  chi2_plikEE*\t\\chi^2_{\\rm plikEE}\n",
      "\n",
      "line  56 :  chi2_plikEE*\t\\chi^2_{\\rm plikEE}\n",
      "        :  chi2_prior*\t\\chi^2_{\\rm prior}\n",
      "        :  chi2_prior*\t\\chi^2_{\\rm prior}\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# load the data\n",
    "planckDir = 'Planck_data/'\n",
    "\n",
    "covmatFile1 = 'base_plikHM_EE_lowTEB.covmat'\n",
    "paramNamesFile1 = 'base_plikHM_EE_lowTEB.paramnames'\n",
    "covmatFile2 = 'base_plikHM_EE.covmat'\n",
    "paramNamesFile2 = 'base_plikHM_EE.paramnames'\n",
    "covmatFile3 = 'base_plikHM_TTTEEE_lowTEB.covmat'\n",
    "paramNamesFile3 = 'base_plikHM_TTTEEE_lowTEB.paramnames'\n",
    "\n",
    "# data from pages:\n",
    "# https://wiki.cosmos.esa.int/planckpla2015/index.php/Cosmological_Parameters\n",
    "# http://pla.esac.esa.int/pla/#cosmology\n",
    "\n",
    "covmat1 = np.loadtxt(planckDir+covmatFile1)\n",
    "text_file = open(planckDir+paramNamesFile1, \"r\")\n",
    "paramNames1 = text_file.readlines()\n",
    "#print paramNames1\n",
    "print len(paramNames1)\n",
    "text_file.close()\n",
    "\n",
    "covmat2 = np.loadtxt(planckDir+covmatFile2)\n",
    "text_file = open(planckDir+paramNamesFile2, \"r\")\n",
    "paramNames2 = text_file.readlines()\n",
    "#print paramNames2\n",
    "print len(paramNames2)\n",
    "text_file.close()\n",
    "\n",
    "covmat3 = np.loadtxt(planckDir+covmatFile3)\n",
    "text_file = open(planckDir+paramNamesFile3, \"r\")\n",
    "paramNames3 = text_file.readlines()\n",
    "#print paramNames3\n",
    "print len(paramNames3)\n",
    "text_file.close()\n",
    "\n",
    "for lineNum in range(len(paramNames2)): #13):\n",
    "    print 'line ',lineNum,': ',paramNames1[lineNum], \\\n",
    "                  '        : ',paramNames2[lineNum], \\\n",
    "                  '        : ',paramNames2[lineNum]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "covmat1 shape:  (13, 13)\n",
      "covmat2 shape:  (13, 13)\n",
      "covmat3 shape:  (33, 33)\n",
      "covmat1:  ['omegabh2\\t\\\\Omega_b h^2\\n', 'omegach2\\t\\\\Omega_c h^2\\n', 'theta\\t100\\\\theta_{MC}\\n', 'logA\\t{\\\\rm{ln}}(10^{10} A_s)\\n', 'ns\\tn_s\\n', 'tau\\t\\\\tau\\n']\n",
      "covmat2:  ['omegabh2\\t\\\\Omega_b h^2\\n', 'omegach2\\t\\\\Omega_c h^2\\n', 'theta\\t100\\\\theta_{MC}\\n', 'logA\\t{\\\\rm{ln}}(10^{10} A_s)\\n', 'ns\\tn_s\\n', 'tau\\t\\\\tau\\n']\n",
      "covmat3:  ['omegabh2\\t\\\\Omega_b h^2\\n', 'omegach2\\t\\\\Omega_c h^2\\n', 'theta\\t100\\\\theta_{MC}\\n', 'logA\\t{\\\\rm{ln}}(10^{10} A_s)\\n', 'ns\\tn_s\\n', 'tau\\t\\\\tau\\n']\n"
     ]
    }
   ],
   "source": [
    "# get indices of interest from the above output\n",
    "print 'covmat1 shape: ', covmat1.shape\n",
    "print 'covmat2 shape: ', covmat2.shape\n",
    "print 'covmat3 shape: ', covmat3.shape\n",
    "\n",
    "# covariance matrix is only 20x20.  If I use this then I'll have to use theta_mc instead of H0.\n",
    "selectIndices = [0,1,2,4,5,3] # used later to put them in order compatible with my Fij\n",
    "print 'covmat1: ',[paramNames1[si] for si in selectIndices]\n",
    "print 'covmat2: ',[paramNames2[si] for si in selectIndices]\n",
    "print 'covmat3: ',[paramNames3[si] for si in selectIndices]\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(13, 13)\n",
      "(13, 13)\n",
      "(33, 33)\n",
      "[14392681.1929545   1049528.32399962  1437035.41942155    46769.37259216\n",
      "    15165.26890736   189447.056113  ]\n",
      "[1.43899404e+07 1.08746731e+06 1.38833734e+06 4.72009326e+04\n",
      " 1.22552384e+04 1.90637457e+05]\n",
      "[93644638.67707175  4268555.25616327 14958494.88509028   804119.32282487\n",
      "   369396.33712423  3204935.51385855]\n",
      "[14392681.1929545   1049528.32399962  1437035.41942155    46769.37259216\n",
      "    15165.26890736   189447.056113  ]\n",
      "[1.43899404e+07 1.08746731e+06 1.38833734e+06 4.72009326e+04\n",
      " 1.22552384e+04 1.90637457e+05]\n",
      "[93644638.67707175  4268555.25616327 14958494.88509028   804119.32282487\n",
      "   369396.33712423  3204935.51385855]\n"
     ]
    }
   ],
   "source": [
    "# fix the nuissance parameters\n",
    "\n",
    "# invert covariance matrices\n",
    "invCov1 = np.linalg.inv(covmat1)\n",
    "invCov2 = np.linalg.inv(covmat2)\n",
    "invCov3 = np.linalg.inv(covmat3)\n",
    "\n",
    "# print shapes and diagonals\n",
    "print invCov1.shape\n",
    "print invCov2.shape\n",
    "print invCov3.shape\n",
    "print np.array([invCov1[si,si] for si in selectIndices])\n",
    "print np.array([invCov2[si,si] for si in selectIndices])\n",
    "print np.array([invCov3[si,si] for si in selectIndices])\n",
    "\n",
    "#subInvCov = [invCov[si1][si2] for si1 in selectIndices for si2 in selectIndices]\n",
    "#subInvCov = [[invCov[:,si1] for si1 in selectIndices[si2]] for si2 in selectIndices]\n",
    "#subInvCov1 = [invCov[si1] for si1 in selectIndices]\n",
    "#print subInvCov1\n",
    "#subInvCov = [subInvCov1[:,si2] for si2 in selectIndices]\n",
    "\n",
    "nParams = selectIndices.__len__()\n",
    "subInvCov1 = np.empty((nParams,nParams))\n",
    "subInvCov2 = np.empty((nParams,nParams))\n",
    "subInvCov3 = np.empty((nParams,nParams))\n",
    "# pull out the parameters indicated in selectIndices in order indicated there\n",
    "for si1 in range(6):\n",
    "    for si2 in range(6):\n",
    "        subInvCov1[si1,si2] = invCov1[selectIndices[si1],selectIndices[si2]]\n",
    "        subInvCov2[si1,si2] = invCov2[selectIndices[si1],selectIndices[si2]]\n",
    "        subInvCov3[si1,si2] = invCov3[selectIndices[si1],selectIndices[si2]]\n",
    "#print subInvCov\n",
    "print np.diag(subInvCov1)\n",
    "print np.diag(subInvCov2)\n",
    "print np.diag(subInvCov3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get desired combinations before finding sigmas\n",
    "\"\"\"\n",
    "covmatFile1 = 'base_plikHM_EE_lowTEB.covmat'\n",
    "paramNamesFile1 = 'base_plikHM_EE_lowTEB.paramnames'\n",
    "covmatFile2 = 'base_plikHM_EE.covmat'\n",
    "paramNamesFile2 = 'base_plikHM_EE.paramnames'\n",
    "covmatFile3 = 'base_plikHM_TTTEEE_lowTEB.covmat'\n",
    "paramNamesFile3 = 'base_plikHM_TTTEEE_lowTEB.paramnames'\n",
    "\"\"\"\n",
    "Fij_plikHM_EE_lowTEB = subInvCov1\n",
    "Fij_plikHM_EE = subInvCov2\n",
    "Fij_plikHM_TTTEEE_lowTEB = subInvCov3\n",
    "\n",
    "Fij_plikHM_lowTEB = Fij_plikHM_EE_lowTEB - Fij_plikHM_EE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.00013805 0.00139624 0.00032041 0.03085309 0.00389213 0.01601033]\n",
      "[0.00315051 0.00444641 0.00440098 0.01898294 0.01669401 0.02325072]\n"
     ]
    }
   ],
   "source": [
    "# re-invert and get sigmas\n",
    "getSigmas = lambda Fij: np.sqrt(np.diag(np.linalg.inv(Fij)))\n",
    "sigmas_TTTEEE_lowTEB = getSigmas(Fij_plikHM_TTTEEE_lowTEB)\n",
    "sigmas_lowTEB = getSigmas(Fij_plikHM_lowTEB)\n",
    "print sigmas_TTTEEE_lowTEB\n",
    "print sigmas_lowTEB"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# save that sub-Fisher matrix to use in another notebook\n",
    "# I will want it in \"Rebuild covar with noise\"\n",
    "#np.savetxt('PlanckInvCovFromMCMC.txt',subInvCov)\n",
    "np.savetxt('Planck_TTTEEE_lowTEB_FromMCMC.txt',Fij_plikHM_TTTEEE_lowTEB)\n",
    "np.savetxt('Planck_lowTEB_FromMCMC.txt',Fij_plikHM_lowTEB)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Transform to H0 parameter via Jacobian"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# can I get a jacobian to translate theta_mc to H0?\n",
    "# yes.  Use discrete derivative with CAMB.  \n",
    "# dtheta/dparams vector below from notebook \"Rebuild covar with noise\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# first need to convert 100*theta_MC to theta_MC\n",
    "mcIndex = 2\n",
    "mcFac = 100\n",
    "Fij_plikHM_TTTEEE_lowTEB[:,mcIndex] *= mcFac\n",
    "Fij_plikHM_TTTEEE_lowTEB[mcIndex,:] *= mcFac\n",
    "Fij_plikHM_lowTEB[:,mcIndex] *= mcFac\n",
    "Fij_plikHM_lowTEB[mcIndex,:] *= mcFac"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00\n",
      "   0.00000000e+00  0.00000000e+00]\n",
      " [ 0.00000000e+00  1.00000000e+00  0.00000000e+00  0.00000000e+00\n",
      "   0.00000000e+00  0.00000000e+00]\n",
      " [-2.48333991e-02  1.04570657e-02  2.96219550e-05  0.00000000e+00\n",
      "   0.00000000e+00  0.00000000e+00]\n",
      " [ 0.00000000e+00  0.00000000e+00  0.00000000e+00  1.00000000e+00\n",
      "   0.00000000e+00  0.00000000e+00]\n",
      " [ 0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00\n",
      "   1.00000000e+00  0.00000000e+00]\n",
      " [ 0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00\n",
      "   0.00000000e+00  1.00000000e+00]]\n"
     ]
    }
   ],
   "source": [
    "#dtheta/dparams\n",
    "nCosParams = 6\n",
    "# paramList after transforming to (H0, 10^9 As) from (cosmomc_theta, ln(10^10*As) )\n",
    "paramList = ['ombh2','omch2','H0','10^9 As','ns','tau'] \n",
    "dthetas = [-0.024833399132992862, 0.010457065703734872, 2.9621955038288544e-05, 0.0, 0.0, 0.0]\n",
    "\n",
    "# The only derivatives in the Jacobian are on the dtheta/dparam row.\n",
    "mcIndex = 2\n",
    "jacobian = np.diag(np.ones(nCosParams))\n",
    "jacobian[mcIndex,:] = dthetas#*100\n",
    "#jacobian[:,mcIndex] = dthetas#*100\n",
    "\n",
    "# multiply these all by d{100thetamc}/d{thetamc} = 100\n",
    "#jacobian[mcIndex,:] *= 100\n",
    "#jacobian[:,mcIndex] *= 100\n",
    "\n",
    "print jacobian"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# jacobian-it-up\n",
    "Fij_plikHM_TTTEEE_lowTEB = np.dot(jacobian.T,np.dot(Fij_plikHM_TTTEEE_lowTEB,jacobian))\n",
    "Fij_plikHM_lowTEB = np.dot(jacobian.T,np.dot(Fij_plikHM_lowTEB,jacobian))\n",
    "\n",
    "#print Fij_plikHM_TTTEEE_lowTEB\n",
    "#print Fij_plikHM_lowTEB"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Transform to $10^9 As$ via Jacobian"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# only one partial derivative needed\n",
    "# d( ln(10^10*As)) / d(10^9*As) = 1/(10^9*As)*d(10*10^9*As)/d(10^9*As) = 1/(10^9*As)\n",
    "TenToTheNineAs_fid = 2.2\n",
    "AsIndex = 3\n",
    "Fij_plikHM_TTTEEE_lowTEB[:,AsIndex] *= 1/TenToTheNineAs_fid\n",
    "Fij_plikHM_TTTEEE_lowTEB[AsIndex,:] *= 1/TenToTheNineAs_fid\n",
    "Fij_plikHM_lowTEB[:,AsIndex] *= 1/TenToTheNineAs_fid\n",
    "Fij_plikHM_lowTEB[AsIndex,:] *= 1/TenToTheNineAs_fid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Load  FisherMatrix objects that have the CMB-S4 fiducial parameters for combination with this result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "IMPORT ERROR: /opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/quicklens/mpi.pyc (No module named pypar). Could not load pbs or pypar. MPI will not be used.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['ombh2', 'omch2', 'H0', 'As', 'ns', 'tau', 'mnu', 'w', 'bin1']\n",
      "{'cosmomc_theta': None, 'r': 0, 'kPivot': 0.05, 'omk': 0, 'H0': 69.0, 'mnu': 0.06, 'nnu': 3.046, 'tau': 0.06, 'neutrino_hierarchy': 'normal', 'standard_neutrino_neff': 3.046, 'As': 2.2e-09, 'omch2': 0.12, 'w': -1.0, 'num_massive_neutrinos': 1, 'ns': 0.966, 'ombh2': 0.0222}\n"
     ]
    }
   ],
   "source": [
    "# load saved Fobj\n",
    "import FisherCl as fcl\n",
    "import noiseCl as ncl\n",
    "import crosspower as cp\n",
    "import pickle\n",
    "\n",
    "#filehandler = open('Fisher_7_6X4_w2_wan.obj','r') # AccuracyBoost=3, lmax=3000, lmaxP=5000, nz=3e5; nonlinear=False,\n",
    "#filehandler = open('Fisher_7_6X4_w2n.obj','r') # AccuracyBoost=3, lmax=3000, lmaxP=5000, nz=3e5; nonlinear=False,\n",
    "#filehandler = open('Fisher_7_6X4n.obj','r') # AccuracyBoost=3, lmax=3000, lmaxP=5000, nz=3e5; nonlinear=False,\n",
    "#neutrino_hierarchy = 'normal', usePrimaryCMB = True; from FisherCl branch quickCl\n",
    "# kludge for this one: told cp.getDNDZinterp to use modelDNDZ3; set cp.tophat beesBins=True\n",
    "# biasByBin = False (uses sigma(b)/b), dw = 0.05; X2n version has noise included\n",
    "\n",
    "#filehandler = open('Fisher_7_6X2_w2_wan.obj','r') # AccuracyBoost=3, lmax=3000, lmaxP=5000, nz=3e5; nonlinear=True,\n",
    "#filehandler = open('Fisher_7_6X2_w2n.obj','r') # AccuracyBoost=3, lmax=3000, lmaxP=5000, nz=3e5; nonlinear=True,\n",
    "#filehandler = open('Fisher_7_6X2n.obj','r') # AccuracyBoost=3, lmax=3000, lmaxP=5000, nz=3e5; nonlinear=True,\n",
    "#neutrino_hierarchy = 'normal', usePrimaryCMB = True; from FisherCl branch quickCl\n",
    "# kludge for this one: told cp.getDNDZinterp to use modelDNDZ3; set cp.tophat beesBins=True\n",
    "# biasByBin = False (uses sigma(b)/b), dw = 0.05; X2n version has noise included\n",
    "\n",
    "# this one has TE polarazation to l=5000, but only one galaxy bin\n",
    "filehandler = open('Fisher_4_1Yn.obj','r') # AccuracyBoost=3, lmax=3000, lmaxP=5000; nz=3e5; nonlinear=True,\n",
    "#neutrino_hierarchy = 'normal', usePrimaryCMB = True; from FisherCl branch quickCl\n",
    "# this one used fiducial params and values to match CMB-S4 Science book verification table 8.1\n",
    "\n",
    "# this on like the Yn version but uses 'total' rather than 'unlensed' power spectra\n",
    "#filehandler = open('Fisher_4_1Zn.obj','r') # AccuracyBoost=3, lmax=3000, lmaxP=5000; nz=3e5; nonlinear=True,\n",
    "#neutrino_hierarchy='normal', usePrimaryCMB=True, nMapsP=3, useNoise=True; myClName='total'; \n",
    "#from FisherCl branch quickCl\n",
    "# this one used fiducial params and values to match CMB-S4 Science book verification table 8.1\n",
    "\n",
    "Fobj = pickle.load(filehandler) \n",
    "\n",
    "print Fobj.paramList\n",
    "print Fobj.cosParams"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# break from the main point of this notebook for a moment to check the primary CMB Cls"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(2, 2, 4999)\n",
      "['TT', 'TE', 'EE']\n"
     ]
    }
   ],
   "source": [
    "crossClsP = Fobj.crossClsP\n",
    "#crossClsP = Fobj.crossClsPPlus[:,:,0,:]\n",
    "print crossClsP.shape\n",
    "#print Fobj.fieldNamesP\n",
    "print Fobj.obsNamesP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAEDCAYAAAAVyO4LAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzt3Xd4XFeZ+PHvO12jmVG3XBRbtuOS\n5pKIVAghIR0CDxtIoyxkf4altx8kyz5LYJ/9LW0puywsWWDZBZIQICGQhUAKIZQQY8cldty7XNS7\nNP38/rh3ZHWNpBlp7vj9PI8yM/feuXOOIr969d5zzxFjDEoppZzDNdcNUEopNTUauJVSymE0cCul\nlMNo4FZKKYfRwK2UUg6jgVsppRwmb4FbRL4rIs0isiNH53tCRDpF5PER25eKyAsisk9EfiQivlx8\nnlJKFap8ZtzfA27I4fm+CLxtjO2fB75ijFkBdAB35/AzlVKq4OQtcBtjngPah24TkeV25rxZRH4v\nIquncL6ngZ4R5xPgauAn9qb/Bt44s5YrpVRh88zy590PvMcYs09ELgG+gRV4p6sK6DTGJO3XjcCi\nGbZRKaUK2qwFbhEJAZcDP7YSZQD89r43AZ8d423HjTHXT3TaMbbpPfxKqaI2mxm3Cys7XjdyhzHm\nEeCRaZyzFSgXEY+dddcBJ2bWTKWUKmyzNhzQGNMNHBKRN4NVnxaRtTM8pwF+C9xqb3oH8NiMGqqU\nUgVO8jU7oIg8CFwFVANNwKeBZ4BvAgsAL/CQMWasEslY5/s9sBoIAW3A3caYX4vIMuAhoBLYArzV\nGBPLbW+UUqpw5C1wK6WUyg+9c1IppRwmLxcnq6urTX19fT5OrZRSRWnz5s2txpiabI7NS+Cur69n\n06ZN+Ti1UkoVJRE5ku2xWipRSimH0cCtlFIOo4FbKaUcZrbnKlFKqRlLJBI0NjYSjUbnuilTFggE\nqKurw+v1TvscGriVUo7T2NhIOBymvr6eIXMfFTxjDG1tbTQ2NrJ06dJpnyerUomIfEREdorIDhF5\nUEQC0/5EpZSaoWg0SlVVlaOCNoCIUFVVNeO/FCYN3CKyCPgg0GCMOR9wA7fP6FOVUmqGnBa0M3LR\n7mwvTnqAEhHxAEEKaAa+rcc6eWZ301w3QymlZs2kNW5jzHER+RJwFBgAfmOM+c3I40RkA7ABYPHi\nxblu55i6owne9I0/kjbw8/dfwZq68ln5XKXUma2trY1rrrkGgFOnTuF2u6msrMTtdhOPx2lvb8ft\ndlNTY90IuXHjRny+3C2HO2ngFpEK4A3AUqATayGEtxpjfjD0OGPM/Vgr3NDQ0DArM1f9+UAbafuT\nnthxSgO3UmpWVFVVsXXrVgDuu+8+QqEQH//4xwf3j7Utl7IplbwWOGSMaTHGJLAWPLg8L62Zos1H\nOvC5XayYF2Lrsc65bo5SSs2KbIYDHgUuFZEgVqnkGqAgJiLZdaqHlfNDnL+wjCdf1jq3Umeiz/xi\nJy+f6M7pOc9dGOHTrz8vp+fMpUkzbmPMC1irqL8IvGS/5/48tysrR9v6WFJVSn11KW19cbqjiblu\nklJK5V1WN+AYYz6NtYJNwUim0jR2DHDjBQuorwoCcKS1nwvqyua4ZUqp2VTImXG+OHaukpNdUZJp\nw5LKIEuqSgE40t43x61SSqn8c2zgPtbeD8DiyiALy0oAaOrWpSaVUsXPsXOVNPdYQbq2LECkxIPf\n46K523kTziilnO2+++7LalsuOTbjbrEDd03Yj4hQGwnQpIFbKXUGcG7g7o3h97gI+60/Gmojfi2V\nKKXOCM4N3D2xwWwbYF4kQFOPZtxKqeLn2MDd2hujOuQffD0v7KdZM26l1BnAsYE7k3FnVJX66I0l\niSVTc9gqpZTKv6IJ3BWl1sxbnf1696RSqrg5cjhgIpWmvT9OzZBSSUXQCtwd/XFqI7pAj1Iqf8aa\n1jUzheu2bdtYu3bt4LG3334799xzT04/35GBu7M/gTFQFTo9v20mcLf3xeeqWUqpM8RE07qGQqHB\nffniyFJJ14AVnMtKTq+SXGmXSjr6tFSilCpujsy4uwas4FweHJJxl1pBvL1fM26lzii/ugdOvZTb\nc86/AG783LTeOjAwwLp16wZf33vvvdx22225ahng0MCduQBZPiTjLi/JZNwauJVSc6ekpCTvpRJn\nB+7g6cDts++i7NCMW6kzyzQzYyebtMYtIqtEZOuQr24R+fBsNG48mVLJ0Bo3WEMCNeNWShW7bFZ5\n3wOsAxARN3AceDTP7ZpQ50ACEQgHRgfuNg3cSqk5NLLGfcMNN/C5z+X2r4KplkquAQ4YY47ktBVT\n1NUfJxLw4nbJsO1lJd7BbFwppWbDyClcU6n837091eGAtwMPjrVDRDaIyCYR2dTS0jLzlk2gcyAx\nqkwCEAl46NHArZQqclkHbhHxAbcAPx5rvzHmfmNMgzGmIXMHUb50DSSGXZjMiJR46Y4m8/rZSik1\n16aScd8IvGiMacpXY7LV2T92xh0OeHSld6VU0ZtK4L6Dccoks83KuH2jtkcCXuLJNNGEzhColCpe\nWQVuEQkC1wKP5Lc52ekaSFBWMvq6asTOwjXrVkoVs6wCtzGm3xhTZYzpyneDJpNOGzr744N3Sg4V\nCVjBvHtA69xKqeLluDsne+NJ0oZxL06CZtxKqfxzu91ccMEFg68z07deddVVnDx5kpKSEgDOPvts\nfvKTn+T0sx0XuLvt4X7hwBilEvuGnB4dWaKUyrOJ5iT54Q9/SENDQ94+23HTuvbGrKA88q5JYLDu\n3a1juZVSRcxxGXevnU2H/KObngnmWipR6szx+Y2fZ3f77pyec3Xlaj558ScnPGai6VvvuuuuwVLJ\ntddeyxe/+MWcts9xgbvHzrhDE5RK9OKkUirf5rJU4rjAncm4w2Nk3AGvC69bNONW6gwyWWZcjIqq\nxi0iRAJerXErpYqaYzPusUolYA0J1FElSql8m2j61qE17urqap566qmcfrbjAndPLIkIBL3uMfdH\ndL4SpdQsGG/61meffTbvn+28Ukk0ScjnwTViLu6MsJZKlFJFznmBO5YYt0wCECnx6NSuSqmi5sDA\nnRxzDHdGyO8ZrIMrpYqXMWaumzAtuWi34wJ3TzQ5YcYd8nvpi2ngVqqYBQIB2traHBe8jTG0tbUR\nCARmdB7HXZycNOMOeKyJqNJm3Dq4UsrZ6urqaGxsJN/LJOZDIBCgrq5uRudwXuCOJllQNv5vq7Df\ngzHQn0hNGOCVUs7l9XpZunTpXDdjzmS7kEK5iPxERHaLyC4RuSzfDRvPZBl3qb1P69xKqWKVbUr6\nNeAJY8yt9qLBwTy2aUK90SQh/+i7JjMy9e/eWAKYWR1JKaUK0aSBW0QiwJXAXwMYY+JAPL/NGls6\nbeiNT3xxMjOHid49qZQqVtmUSpYBLcB/icgWEfm2iJSOPEhENojIJhHZlK8LBn3xJMaMPcFURiao\n98V0wWClVHHKJnB7gAuBbxpj1gN9wD0jDzLG3G+MaTDGNNTU1OS4mZbeCaZ0zcjUv61SiVJKFZ9s\nAncj0GiMecF+/ROsQD7rJlpEISOkpRKlVJGbNHAbY04Bx0Rklb3pGuDlvLZqHBMtopBxOuPWwK2U\nKk7Zjir5APBDe0TJQeCd+WvS+DIZd2SCwK3DAZVSxS6rwG2M2Qrkbx2eLA3WuCcYDujzuPB7XPTG\nNXArpYqTo+YqmWwRhQydaEopVcwcFbgHa9yT3MoeCni0xq2UKlqOCtzZjCrJ7NeMWylVrJwVuGMJ\ngj437klm/Qv5PYPZuVJKFRuHBe6JJ5jKCAc8Oie3UqpoOSpwT7aIQkapX2vcSqni5ajA3RtLTjhP\nSYbWuJVSxcxZgTvLjDsU0Bq3Uqp4OStwZ1vj9nuIJ9PEkjpDoFKq+DgqcPdMsohCRia469SuSqli\n5KjA3RtLEs7y4iSgI0uUUkXJMYHbGDOl4YCgU7sqpYqTYwL3QCJFKm2yuzhpl1N0SKBSqhg5JnBn\ne7s7jFwwWCmliotjAndmeF82NW5dBUcpVcyymo9bRA4DPUAKSBpjZn1u7kzGPZXAraNKlFLFKNsV\ncABeY4xpzVtLJpHNIgoZWipRShUz55RKplDjDnrdiOjyZUqp4pRt4DbAb0Rks4hsGOsAEdkgIptE\nZFNLS0vuWmjrzXIRBQCXSwj59LZ3pVRxyjZwX2GMuRC4EXifiFw58gBjzP3GmAZjTENNTU1OGwnQ\nG7XKHtnUuMGeIVAzbqVUEcoqcBtjTtiPzcCjwMX5bNRYMhl3aRYZN+jyZUqp4jVp4BaRUhEJZ54D\n1wE78t2wkXpiSXweFz5Pdn8khHRObqVUkcomfa0FHhWRzPEPGGOeyGurxtAbzW4u7oywZtxKqSI1\naSQ0xhwE1s5CWybUF8tuLu6MkN/Dqa5oHluklFJzwzHDAbOdYCpDly9TShUrxwRuay7uqWXcOqpE\nKVWMHBO4s52LOyMc8NAbT5JOmzy2SimlZp+jAne2QwHBCtzGQF9cs26lVHFxTuCeYqkkHLDmNOnW\ncolSqsg4JnD3THFUScQO3D1RnWhKKVVcHBG4Y8kU8WR6yuO4QefkVkoVH0cE7sy82lMrlWQCt2bc\nSqni4ojAPbhsWWDyubgzwoOlEs24lVLFxRmBewpTumZE7IxbL04qpYqNowL31MZx68VJpVRxckjg\ntoLvVMZxB7wuPC7RUolSqug4InBPZdmyDBEhUuLVjFspVXQcEbinUyrJHK8Zt1Kq2DgjcE8j4wYr\ncHcPaMatlCouWQduEXGLyBYReTyfDRpLbyyJCAR97im9L+z3asatlCo6U8m4PwTsyldDJpKZ0tVe\nhSdrWipRShWjrAK3iNQBNwPfzm9zxtYbm9qyZRnhgF6cVEoVn2wz7q8CnwDSeWzLuHqjU5vSNUMz\nbqVUMcpmlffXAc3GmM2THLdBRDaJyKaWlpacNRDsZcumOKIErLsndTEFpVSxySbjvgK4RUQOAw8B\nV4vID0YeZIy53xjTYIxpqKmpyWkjp7reZEakxIsx0KuLKSilisikgdsYc68xps4YUw/cDjxjjHlr\n3ls2xFSXLcvQqV2VUsXIMeO4p5Nx63wlSqliNKVoaIx5Fng2Ly2ZgFUqyX5K14xMxt09oBm3Uqp4\nFHzGnU4bO3BP7eYb0IxbKVWcCj5w99jzlERKpp9xa41bKVVMCj5wZ+YamU7gnuqCwU/vauLTj+3g\nROfAlD9LKaVmy9Sv+M2ybjvoRqawbFlGmR3sO/snD9x7m3p4zw82k0gZthzr5GfvvQKXa2q32Btj\n+ObvDrDzeDefvGE1i6uCU26zUkpNxgEZd6ZUMvXfMT6Pi1Kfm44sAvd//fEwHpeLv7tpNdsbu3j+\nYNuUP++nLx7nC0/s4X9fOsl7H9isN/4opfKi8AP3DDJugPKgj87++ITHJFJpHt92ghsvmM/bL6sn\nEvDw0xcbp/Q56bThm8/u5/xFEb5w6xp2HO/mz9MI/kopNZnCD9x2jbtsGjVugIpSLx2TBO7tjV30\nxJK89pxaAl43r1k9j+f2tkwpY95yrJMDLX2847J6blm7kEjAw6Nbjk+rzUopNZHCD9z2iJBpZ9wl\nvklLJc8faAXg0mVVAFy5oobW3jgvn+zO+nOe3tWE2yVcd+58Al43r1pRw3P7WjBGyyVKqdwq/MBt\nZ9zTmWQKoDzopWuSVXD+criDVbVhKkt9AFx+thXANx5qz/pzntndzCvqKygLWr9grlxZTVN3jL1N\nvdNqt1JKjafwA3c0QdjvwT3FER4ZFUHfpKWSPad6OG9hZPD1grISFpQF2HKsM6vP6OiLs/tUD69a\ncXpyrUz2/uLRjmm0Wimlxlf4gXsgOa0x3BkVdsadGqde3dWf4FR3lJXzw8O2r19czpYsg+7Wxs7B\n92QsrgxSHvSyLcvgr5RS2Sr8wB1NTGtmwIzyoA9jGHfR4D1NPQCsGhm4z6qgsWOA5p7opJ+x9Wgn\nLoE1dacDt4iwpq6cbY1d0267UkqNpfAD90BiZhl3qfXe8colg4G7dnTGDVZQnsyWY52srA2PmsFw\nbV0Ze5t6iCZSU263UkqNp/ADdzQ57RElYGXcwLgjS/ac6iYc8LCgLDBs+3kLy3C7hO2TZMzptGHb\nsc5hZZKM1fMjpNKG/c16gVIplTuFH7gHEtO6azKjwg7c492Es/dUL6tqw6NWkC/xuVlZG2Zb48QZ\n96G2ProGEqw7a3TgXlkbAtDArZTKqcIP3NHEjDLuimCmVDI64zbGsPtU96j6dsbaujK2N3ZNOBY7\nU0pZv7hi1L766lI8LmGvXY5RSqlcyGax4ICIbBSRbSKyU0Q+MxsNg9Nzcc+kxp0Zm93aGxu1r6k7\nRnc0OX7gPqucroEER9r6xz3/lmMdhP0ezq4JjdrndbtYWl2qY7mVUjmVTcYdA642xqwF1gE3iMil\n+W2WpSeWxBhrtfbpCvk9BH1umrtHB+7dp6w7I0demMxYU1cGMGG5ZOuxTtacVTbuTIIra8Psb9aM\nWymVO9ksFmyMMZmU0Wt/zcp93DOZiztDRJgX9o85rC9Twlg5TuBeWRvG73GNe4FyIJ5i18ke1p81\nukySsaI2xJH2fh1ZopTKmaxq3CLiFpGtQDPwpDHmhTGO2SAim0RkU0tLS04al7lVfSYZN8C8cIDm\nnrEy7h7mhf1U2OWUkbxuF+cvKmP7OBn3jhNdpNJmzAuTGSvmhTEGDrRouUQplRtZBW5jTMoYsw6o\nAy4WkfPHOOZ+Y0yDMaahpqZm9EmmoWtwZsCxA2u2aiJ+WsYI3Hubesatb2esqSvjpeNdJFPpUfsy\nd1auG2MoYMbS6lIADrX2TaXJSik1rimNKjHGdGKt8n5DXlozQuammcxNNNM1L+ynuXt4qSSVNuxr\n6h23vp2xtq6caCLNvjGG9L14pJMlVUGqQ/5x319fba2CM9EFTqWUmopsRpXUiEi5/bwEeC2wO98N\ng9ND+DJjsadrXjhAXzxFX+z0osFH2vqIJdNZZdzAqHKJMYYXj3awfoIyCUDQ56E24teMWymVM9lk\n3AuA34rIduAvWDXux/PbLEtnn5VxlwdnnnEDw+rce06NPUfJSPVVpUQCnlFzjpzsitLcExtz/PZY\n5zisgVsplSOTXvUzxmwH1s9CW0bp6E8Q9Lnxe9wzOs+Ccut29hOdA4M1512nenDJ+CNKMlwua7Ko\nkRn35iN2fXuSjBusOvdTu5qm03SllBqloO+c7OyPz7hMAtYUqzC8zrzrZDdLq0sJeCf/pbCmrozd\nJ4dPFvWnA62EA55h83iPZ0lVKa29cXqiky9arJRSkynowN3RH5/xhUmwFkbwuV0caT9drth1sptz\nFkwedMHKqpP2ZFJg1bef29vKZcuq8Lgn/xYutS9QHm7VC5RKqZkr8MCdyEnG7XYJdZUlHLUz7q6B\nBI0dA1kH7suWV+Fzu3jyZavcsb+5l+OdA7xqRXVW76/PDAls0zq3UmrmCjpwd/bHB6dlnakllcHB\nUsnLJ6xb3c/NoswBEA54eeWKan614xTptOGxrSdwCVx/3vwsP9sK3HqBUimVCwUduK2Me+alEoAV\ntWH2t/QST6bZfMRaBPjCCW5VH+mN6xdxvHOAH206xgMbj/LqlTXMiwQmfyPWFLELygIauJVSOVGw\ngTuVNnRHEznLuC9YVEY8mWZvUw8b7VXdy6bwS+HmCxawen6Yex95ie6BBB+5duWUPr++qpTDWipR\nSuXAzCYByaOugQTGkLOMOzNs73d7W3jhYBt3XLx4Su93u4T/ftfFfO9Ph7l8edWw9SWzUV8d5Nc7\ndUigUmrmCjZwD97unqOMu66ihFW1Yb746z0A3LxmwZTPURsJ8MkbVk/r8+urSmnvi9M1kKBsBrMd\nKqVUwZZKMkuNzfSuyQwR4f9cuQyAK86uomFJ9vXtXMiMLNE6t1Jqpgo34+7LzTwlQ916UR2XLK2k\nNhIYtcZkvg2dJXBtFndbKqXUeAo3cOe4VJJxln0X5WxbUhXEJXBQ5+VWSs1QwZZK2vtyM6VrofB7\n3CyuDHKgRUslSqmZKdjA3dobw+9xEfIX7B8FU7a8JqQr4SilZqyAA3ec6pB/1mvR+bR8XoiDrX2k\n0rOyZKdSqkgVcOCOUR0ef2UZJ1peU0o8meZ4x8BcN0Up5WDZrIBzloj8VkR2ichOEfnQbDSspSdG\nTSi3Fybn2rKaEKALByulZiabjDsJfMwYcw5wKfA+ETk3v806XSopJss1cCulcmDSwG2MOWmMedF+\n3gPsAhbls1GptKG9L1Z0gbuy1EdF0KuBWyk1I1OqcYtIPdYyZi+MsW+DiGwSkU0tLS0zalRHf5y0\ngeoiK5VAZmSJDglUSk1f1oFbRELAT4EPG2O6R+43xtxvjGkwxjTU1NTMqFGtvdaivsV2cRKswL2v\nqQdjdGSJUmp6sgrcIuLFCto/NMY8kt8mQWuPdfNNsZVKwFq8oaM/wanu6Fw3RSnlUNmMKhHgO8Au\nY8yX89+kIRl3EQbu8xdZq+7sOD7qjxallMpKNhn3FcDbgKtFZKv9dVM+G5UJ3DVFGLjPWRBBBHYc\n75r02K7+BJsOtw9bXV4ppSa9n9wY8wdgVm9fbOmN4XO7iJQUz+3uGUGfh+U1IXaemDjjfurlJj74\n0Bb64ynqKkr4/t2XDM4wqJQ6sxXknZOtPXGqQr6iut19qPMXRnjpeOe4Fyhbe2N85OGtLKsp5au3\nraM/nuI9399MLKmZt1KqQAN3U3c064V4nejCJRU0dcc41j72re9feGI3A/EUX71tPW9cv4gvvXkN\ne5p6+P7zR2a5pUqpQlSQgftUd5QFRRy4L19eBcDzB1tH7dtytIOHNzVy9yuXcvY8607Lq1fXcvny\nKv7jdwe13q2UKtDA3RVlflnxBu7lNSFqwn7+dKBt2PZU2vAPj+1kXtjPB65ZMWzf+19zNq29MX6x\n7cRsNlUpVYAKLnD3RBP0xpJFHbhFhFetqObZPS3D6tYPbzrGS8e7+NTN54yah/yy5VUsrynlgY1H\nZ7u5SqkCU3CBu8m+MWVBEQdugDesW0TXQILf7m4GrBV/vvDEbi6ur+SWtQtHHS8i3HnJErYc7WTX\nSR0DrtSZrKAC9z/978uDF+Bqi7jGDXDF8ioWlZfwjWcPkEyl+fufvURvLMln33jeuKNp3rR+ET63\nix9vapzl1iqlCklBBe4HNx7jZ1utGm6xZ9wet4tP3LCK7Y1dXP65Z/jlS6f46LWrWD0/Mu57Kkp9\nXHtuLT/bepx4Mj2LrVVKFZKCCtwlPjddAwmg+DNusMol//jG81lWU8p9rz+X97x62aTvubWhjva+\nOM/sbpqFFiqlClFB3ZoY9LkBa97qgNc9x62ZHW+7dAlvu3RJ1sdfuaKG2oifhzc1csP5C/LYMqVU\noSqojDvos36PLK4MznFLCpfbJfzVhXU8u6eZZp1hUKkzUoEFbivLXlKlgXsit15UR9rAI1uOz3VT\nlFJzoKACt8dljaZYohn3hJbVhHhFfQUPbzqW9YIMqbQu3KBUsSiowJ0Zw720RmfBm8ybLzqLgy19\nvHi0c8LjfrPzFNd95Xcs/7tfcv1XnuOP+0ffZq+UcpaCCtwN9ZUAXLS4co5bUvhuWrOAoM/Nw385\nNu4x3/vjITZ8fzOC8N6rlpNIpXn7dzfy5Ms6IkUpJ8tmBZzvikiziOzId2M+/fpzeeS9l7NYa9yT\nCvk9vHH9Ih7Z0six9v5R+7/zh0Pc94uXuf68Wn7xgVfyiRtW89j7r+C8hRE+/NCWMd+jlHKGbDLu\n7wE35LkdAIQDXi5cXDEbH1UUPnj1Ctwu4bOPvzys1v2t3x3gHx9/mRvPn8/X77wQn8f63xwOePnm\nWy/CAP/w2A5dsFgph5o0cBtjngPaZ6EtaormlwX42LWrePLlJv7fL3dxqLWPTz+2g3/+1W5et2YB\n/3rHerzu4f+LF5WX8NFrV/LbPS08tat5jlqulJoJySbrEpF64HFjzPkTHLMB2ACwePHii44c0Un/\nZ4Mxhk/9bAcPvHB61sB3XbGUT918Dm7X2HOeJFNprv3Kc/g9Ln75wVfhGuc4pdTsEZHNxpiGrI7N\nVeAeqqGhwWzatCmbQ1UOGGN48WgHe0710lBfwcra8KTveWzrcT700Fa+fud6Xrdm9GyESqnZNZXA\nXVCjStT0iAgXLankzksWZxW0AV63ZiEra0N8+cm9JFM6YZVSTqKB+wzldgkfvXYlB1v6eGyrrqqj\nlJNkMxzwQeB5YJWINIrI3flvlpoN1583n/MXRfjq03t1mlilHCSbUSV3GGMWGGO8xpg6Y8x3ZqNh\nKv9EhI9dt4pj7QM8vGn8G3mUUoVFSyVnuKtW1nDRkgr+7Zl9uoK8Ug6hgfsMJyJ8/LpVNHXH+Pbv\nD851c5RSWdDArbhseRU3XTCff316P/ube+a6OUqpSWjgVgB85pbzCfrdvOcHLw4uHzeWtt4Y+5p6\n6OiLz2LrlFJDFdTSZWru1IT9fOOuC3nHdzdy27ee51/espbzFpYRT6bZ1tjJs3uaeWZ3C7tOdg++\n59wFEd71yqW8af0ivftSqVmU1Z2TU6V3TjrXc3tb+MCDW+gaSFAR9NIXTxFPpnG7hIYlFVy1ah6L\nKko43jHAz7edYNfJbi5aUsHXbl9HXYXO6qjUdOX8lvep0sDtbF39CR7d0sje5l5Cfg8XLi7nsuXV\nlJV4hx2XThse2XKcz/x8Jy6X8OW3rOWac2rnqNVKOZsGbjWrDrf28b4HXmTniW7e/epl/N/rVuFx\n6+UTpaZC5ypRs6q+upSf/u3l3HnJYr71u4Pc+Z8vcKJzYK6bpVTR0oxb5dSjWxr5u0d2YDC87dIl\nvGHdIlbND+N1u+iPJzneMcCJrigtPTE6+uIEfG4WlQdYW1dOVcg/181Xas5oqUTNqWPt/XzpN3t4\nfPtJUmmDS8DjchGfZBbCS5dV8uaLzuL1axcOrtqj1JlCA7cqCK29MZ7b28Kh1j4SKUM44KGuooSF\n5SXMC/spD/qIJVIcau3j+YNtPLb1BIda+1hYFmDDlcu47RWLKfG557obSs0KDdzKkYwxPLu3hX9/\nZj+bjnRQEfRy1yVLePtlS5iz7Z7mAAANxklEQVQXCcx185TKqzM7cHefhG0PQPtBqFwGa++EyIK5\naYuato2H2vnP3x/kqV1NuMQaQ/6qFdWsrA0zLxIg6HPjc7vwuAWv24XHJXjcrsFtHpcgojcFKec4\ncwP3pv+CJ+6F5ACEaqG3CbylcPOXYN2ds98eNWNH2vr48aZGntrVxO5T2c+jUlbiZcW8ECtqw6yq\nDbGyNsyK2jDVIZ8GdFWQ8rHm5A3A1wA38G1jzOcmOn7WA3c6DU/fB3/8Giy/Bm7+F6hcamXdP/8g\nHP49vPY+eOVHZq9Ng21LwUs/hpcfg4FOmH8+XPROzLxzSKQTpEyKgDugwSQLPdEEB1r6aO+L0W/f\n0ZlMGRJp+zGVJpk2JJJpTnZH2dfUw96m3mFzr1SW+qivClJZ6qeq1EfA60JEcLusL4/LyuC9diYf\nDnipCvmoDvmpCfmpDvsI+nSmCJV7OQ3cIuIG9gLXAo3AX4A7jDEvj/eeWQ3ciSj87G9h5yPQ8C64\n8Yvgtv5hGWMYiHXT+Yv307XvV3Re8m7M6psp9ZUS8oaI+CJE/BH87hHD0JIxEtsfovHIcxw2MY6W\nL+RoSYijPcdoi7aRSCVwi5vKkkqqAlXML51vfQWtx9rSWioDlaS6jtP503dyuGU7hyK1HCwJcTjW\nxlGPmxNeLyms773X5aWqpIr6SD3Ly5ezrGwZiyOLWVi6kPml8/G5fVZfkgO0R9tpj7bTOtBK60Ar\nbQNttAy00Jvoxevy4nP7qApUMS84j3nBedSU1FBVUkVloBKf2ze8m+kk3fFuumPd9CX6GEgOMJAc\nIJaK4XF5CHgCBNwBgt4gIW+IUm8pQW8Qr2v4HZSFzBhDS0+MvU297GnqYV9TD0fb+2nvi9PRHyeW\nTJNKG4yBpP0LIJme+N9E0Oe2AnnYT7Ud1IM+Nz6PC5/bevTa5Rqvx4XXZZdv3C68dknH4xZ8bhd+\nj8t6n8cq82Se+91uvB7BJZkvawrezKMqPrkO3JcB9xljrrdf3wtgjPnn8d4z3cD9Lw9eTxRIi5AS\nIeUNkg5WkSwpJ40hZVKk0inSJk3SJEknoqRObScV6yZdsYRkqJa0SRNNRemMddIV6yKRHn+mu4yA\nO0DEHyHii0BigO6eE7RLmuSQfyBlRlhcsYLq8CL8bj/JdHIwgJ7qO0U8PXy2PLe4SJnhw9+CniBL\nQotY0t9FXfM+gpHFuNa8mS5StPa3cqjrEAe6DjCQHH7zimD94x15PmsfVPjKCPnLSKaTg31Pj3Fs\nyBvC4/JAKkEiGaWP6S2cEHAHKPWWUuoNUioeTDJKPDlAPJ0knk6SwJASSAMpA2kMIoLfHcDv9uP3\n+PG7/QTs1x4MpJNIOoWkE0gqBekEJpXApBMYcYG4weXGuNwYxN4miLitwIbBZRh8FGNwmTRuQEwa\nlzEIp3/WZfBRhjy3XhuM9cpAyhhSKUMybQX4ZNqQTBnr+eB2SBsr+BvMkDOObeTe06/N6ddj/LMc\n/T77k2T4MeNuy/x3xL7Mk6Hnl+F7B98z3vFjHsPw7+/INlnnEMbYPaqdUzHR22Ssb2wOBdwlfP5d\nv5jWe6cSuLP5m28RMHRdq0bgkjE+dAOwAWDx4sXZfPZw6TS/jJ0kZv9jc9v/8Dzt4HK5cftCuH1h\nXL4gbgR3rAd390ncJo1r3mo8pfPwixuXy0Wtu5a1/rWU+cso85dR7i+nzFNK2R+/jvvYn+m97L30\nLrmU7lg33fFuumJddEfb6Tq+EVf7IcLuAFXLXsvSpVdTH17CksYtlD35aWjaDG/9FCy6cFjTjTF0\nxDo41XfK+uo4QMvGb+If6KGs4W4WL7mK5eXLqQ3Wns6Wtj9slXHavwlv+R8462Lr22DSNPU10djb\nyPHe4zR1HiZx5A+kT7xIOJmkonIZleVnU42b6pZ9VDS+iNd1As65BS55N5x1CUmTom2gjeb+Zisz\nj7bS1nGIjsY/k27ZA4k+PC4fkcgiIqXziXhKCMX7Kek4RqDtAIF0kiQuolVLic5bTX+wgj63m95E\nP33RDvr6munt6aAv3kifgMsYvIDfGHwGvCaN2w6ibsBlrDgUdXuIe0uIut3E0imiJknMpIgLZMLq\nyH9WIi7EpBHM8IBgrH+gKYEkkEZIA2mxnht7nxmyD8DI6c8YGsqHfu7QYxhrP1j3HLvGjK9DjpNh\nr8YKG2aMKDPWOUe3J7uolt8wNTufORd9mK5wenb+Gsom434zcL0x5m/s128DLjbGfGC89+SsVNLX\nBoeehR2PwP6nIBkdvr/uYnj916D23OzOl4zBQ3fB/ifhonfCqz4GwUrY/b/w9Geh6xg03G3VwwOR\n4e9tPwj/8wbo74C7fgxLLhv7M/rb4b9vgbZ91nFLrxy/PSe3w8Nvg67jcNU9cNn7wFti7Yt2w8Zv\nwR//DWLdsOY2eM29UFE/ol2HYON/wpYfQKwL5l8A57zB+p6IG9oPwP6n4eCzYNJw9mth/V2w8kbw\njjHELt4HxzfDkefh6PPQ+BeI957eL27r+kHVCqhZBTWroWYlVC4HfwRcLkgMWPX8/jboa4beFuhr\nOf08OWBdNPYFIVAOJRXjfJWD24tVx4hCrBfiPdb/x2TUekwMgLjA47e+3JlH34hHv9U2pQqUY0sl\nE4r3WxcZW/eBywNnvQIWXjj1v6eScXjms/CnrzPsd/m886zRJ0suH/+9Xcet4N11DG79Lqy+efj+\nvjZrf+teuONBOPuaydsz0AGPvR92Pw6BMuuXUSoGxzZawWnVTXD130PteROfJ95nZfGbvwcntw3v\nW9UKOOd1cNFfjw78kzEGYj3WLw9vEPxhK5gqpXIq14Hbg3Vx8hrgONbFyTuNMTvHe48jbsBpPwh7\nf20FvEUXWZmxK4u79Ppa4YG3wIktcOUn4IoPWZnjkefh0XdbQxBvfyC7oD3U4T/A1gegaYf1i2lR\nA6y7Axaun3rfBjqg47A1oqXsLAjrVKtKFbp8DAe8Cfgq1nDA7xpj/mmi4x0RuGci3ge/+JA1zM8b\ntDLlnpMQWQRv/t5gvVoppbKV64uTGGN+CfxyRq0qJr5S+Ktvwyv+BnY+CtEuK2tfd6e1Tyml8kjv\nJJiJxZdaX0opNYv0MrtSSjmMBm6llHIYDdxKKeUwGriVUsphNHArpZTDaOBWSimH0cCtlFIOo4Fb\nKaUcJi9Ll4lIC3Bkmm+vBlpz2BwnONP6fKb1F7TPZ4qZ9HmJMaYmmwPzErhnQkQ2ZXu/frE40/p8\npvUXtM9nitnqs5ZKlFLKYTRwK6WUwxRi4L5/rhswB860Pp9p/QXt85liVvpccDVupZRSEyvEjFsp\npdQENHArpZTDFEzgFpEbRGSPiOwXkXvmuj0zISLfFZFmEdkxZFuliDwpIvvsxwp7u4jIv9r93i4i\nFw55zzvs4/eJyDvmoi/ZEpGzROS3IrJLRHaKyIfs7UXbbxEJiMhGEdlm9/kz9valIvKC3f4fiYjP\n3u63X++399cPOde99vY9InL93PQoOyLiFpEtIvK4/brY+3tYRF4Ska0issneNrc/18aYOf/CWsvy\nALAM8AHbgHPnul0z6M+VwIXAjiHbvgDcYz+/B/i8/fwm4FeAAJcCL9jbK4GD9mOF/bxirvs2QZ8X\nABfaz8NYC0yfW8z9ttsesp97gRfsvjwM3G5v/w/gb+3n7wX+w35+O/Aj+/m59s+8H1hq/1twz3X/\nJuj3R4EHgMft18Xe38NA9Yhtc/pzPeffFLtTlwG/HvL6XuDeuW7XDPtUPyJw7wEW2M8XAHvs598C\n7hh5HHAH8K0h24cdV+hfwGPAtWdKv4Eg8CJwCdadcx57++DPNvBr4DL7ucc+Tkb+vA89rtC+gDrg\naeBq4HG7/UXbX7t9YwXuOf25LpRSySLg2JDXjfa2YlJrjDkJYD/Os7eP13fHfk/sP4nXY2WgRd1v\nu2ywFWgGnsTKHjuNMUn7kKHtH+ybvb8LqMJZff4q8Akgbb+uorj7C2CA34jIZhHZYG+b05/rQlks\nWMbYdqaMUxyv7478nohICPgp8GFjTLfIWN2wDh1jm+P6bYxJAetEpBx4FDhnrMPsR0f3WUReBzQb\nYzaLyFWZzWMcWhT9HeIKY8wJEZkHPCkiuyc4dlb6XCgZdyNw1pDXdcCJOWpLvjSJyAIA+7HZ3j5e\n3x33PRERL1bQ/qEx5hF7c9H3G8AY0wk8i1XXLBeRTFI0tP2DfbP3lwHtOKfPVwC3iMhh4CGscslX\nKd7+AmCMOWE/NmP9cr6YOf65LpTA/RdghX112od1IePnc9ymXPs5kLmS/A6sGnBm+9vtq9GXAl32\nn16/Bq4TkQr7ivV19raCJFZq/R1glzHmy0N2FW2/RaTGzrQRkRLgtcAu4LfArfZhI/uc+V7cCjxj\nrILnz4Hb7VEYS4EVwMbZ6UX2jDH3GmPqjDH1WP9GnzHG3EWR9hdAREpFJJx5jvXzuIO5/rme68L/\nkGL9TVgjEQ4An5rr9sywLw8CJ4EE1m/au7Fqe08D++zHSvtYAf7d7vdLQMOQ87wL2G9/vXOu+zVJ\nn1+J9affdmCr/XVTMfcbWANssfu8A/gHe/syrEC0H/gx4Le3B+zX++39y4ac61P292IPcONc9y2L\nvl/F6VElRdtfu2/b7K+dmdg01z/Xesu7Uko5TKGUSpRSSmVJA7dSSjmMBm6llHIYDdxKKeUwGriV\nUsphNHArpZTDaOBWSimH+f+AH9DnI7ibjwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10906c2d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot them\n",
    "ellsP = Fobj.ellsP\n",
    "\n",
    "plt.plot(ellsP,crossClsP[0,0],label='TT')\n",
    "plt.plot(ellsP,crossClsP[0,1],label='TE')\n",
    "plt.plot(ellsP,crossClsP[1,1],label='EE')\n",
    "plt.legend()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# the familiar shape of Cl^TT here shows that it really is Dl=ell*(ell+1)*Cl/2pi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(2, 2, 4999)\n"
     ]
    }
   ],
   "source": [
    "# convert the Dl to Cl in crossCls, the re-compute everything that follows from it\n",
    "DlFac = ellsP*(ellsP+1)/(2*np.pi)\n",
    "crossClsMod = crossClsP/DlFac\n",
    "print crossClsMod.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEDCAYAAADOc0QpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAFlxJREFUeJzt3X2QXfV93/H3d1dPYMkGpHWNtRDJ\nQSlRYp5mg0ndaRTbxIJxoe04sTTEIQmJJmlJ23HcWow7QMg/xswkHk9wsCZhSDIEjIkbK4xs+Qlq\n94GHpTwYgWUUcMpG2FoE2MEGC4lv/7hnxWW5Z8/V7l3d/V29XzM7ex5+e87vt3P10dnv/d1zIjOR\nJA2WoX53QJLUe4a7JA0gw12SBpDhLkkDyHCXpAFkuEvSAOpruEfEjRGxLyIe6dHxvhgRz0fEHdO2\nr42IeyLi8Yj4TEQs6cX5JGmh6veV+03Axh4e7zrggx22Xwv8cWauA54DLuvhOSVpwelruGfm14Fn\n27dFxE9WV+D3R8Q3IuL0IzjeV4F/mna8AN4F3F5t+gvg38yt55K0sC3qdwc62Ab8TmY+HhHvAD5F\nK5xnayXwfGYerNYngNVz7KMkLWgLKtwjYjnwL4DPti64AVha7ft3wDUdfuwfM/O9Mx22wzbvuSBp\noC2ocKdVJno+M8+aviMzPwd8bhbHfAY4ISIWVVfvo8DeuXVTkha2fr+h+hqZ+QPgyYj4ZWjVyyPi\nzDkeM4E7gfdXmy4FPj+njkrSAhf9vCtkRNwCbABWAd8DrgK+BvwpcDKwGLg1MzuVYzod7xvA6cBy\nYD9wWWbujIi3AbcCJwEPAL+amT/u7WgkaeHoa7hLkubHgirLSJJ6o29vqK5atSrXrFnTr9NLUpHu\nv//+ZzJzpKld38J9zZo1jI+P9+v0klSkiPiHbtpZlpGkAWS4S9IAMtwlaQAttE+oSlLPvPzyy0xM\nTPDSSy/1uytHbNmyZYyOjrJ48eJZ/bzhLmlgTUxMsGLFCtasWUPb/aoWvMxk//79TExMsHbt2lkd\nw7KMpIH10ksvsXLlyqKCHSAiWLly5Zz+4jDcJQ200oJ9ylz7XVy4f/t7/8QffWk3z7zgrWEkqU5x\nNffHv/cCn/zaHt535ltZtXxpv7sjSR3t37+fd7/73QB897vfZXh4mJNOOonh4WEOHDjAs88+y/Dw\nMCMjrQ+b3nvvvSxZ0rvHOxcX7pJUgpUrV/Lggw8CcPXVV7N8+XI+/OEPH97faVsvFVeWkSQ1K/bK\n3TsVSzoSf/B3u3h07w96esz1b30jV/3rn+npMXuluCv3Qt/4lqSjqtgrd0k6Egv1Cnu+FHflPiWx\nLiNJdYoLd6syktTMsowkzbOrr766q229VNyVuySpWbHh7lRISapXXLg7FVKSmjWGe0TcGBH7IuKR\nhnY/FxGHIuL9veueJGk2urlyvwnYOFODiBgGrgV29qBPXbEsI0n1GsM9M78OPNvQ7PeAvwH29aJT\nM7MuI0lN5jwVMiJWA/8WeBfwcw1ttwBbAE499dS5nlqSFqxOt/ydur3vQw89xJlnnnm47aZNm9i6\ndWtPz9+Lee6fAD6SmYeanhySmduAbQBjY2MWViQNrJlu+bt8+fLD++ZLL8J9DLi1CvZVwIURcTAz\n/7YHx67l7Qckqd6cwz0zDz+aOyJuAu6Yz2B3KqSkWfnCVvjuN3t7zLe8HS742BH/2IsvvshZZ511\neP2KK67gAx/4QC971hzuEXELsAFYFRETwFXAYoDMvKGnvZGkY8Bxxx3X/7JMZm7u9mCZ+etz6s0R\ncCqkpCMyiyvskpX3CdV+d0CSCuBdISXpKJtec9+4cSMf+1hv/7Iw3CVpnk2/ve+hQ4fm/ZzFlWUk\nSc2KC/emD0pJkgoMd0lSs2LD3amQklSvuHC3KCNJzYoLd0lSM6dCStI8GR4e5u1vf/vh9alb+27Y\nsIGnn36a4447DoDTTjuN22+/vafnLjbcvSukpIVupnvI3HzzzYyNjc3buYsryzgTUpKaFXvlLklH\n4tp7r+Vbz36rp8c8/aTT+ci5H6ndP9OtfS+55JLDZZnzzz+f6667rqd9KzbcnQopaaHrZ1mmuHC3\nLCNpNma6wh5ExdXcJUnNirtyl6RSzHRr3/aa+6pVq/jKV77S03N385i9G4H3Afsy82c77L8EmPp7\n5wXgdzPzoZ72sgNL7pIWurpb+951113zfu5uyjI3ARtn2P8k8AuZeQbwh8C2HvSrVngDAklq1M0z\nVL8eEWtm2P+/21bvBkbn3i1J0lz0+g3Vy4Av9PiYHaVzISV1odSsmGu/exbuEfGLtMK9dr5RRGyJ\niPGIGJ+cnJzliWb3Y5KOPcuWLWP//v3FBXxmsn//fpYtWzbrY/RktkxEnAH8GXBBZu6va5eZ26hq\n8mNjY2X9tiUVZ3R0lImJCWZ9MdlHy5YtY3R09lXuOYd7RJwKfA74YGZ+e67Hk6ReWbx4MWvXru13\nN/qim6mQtwAbgFURMQFcBSwGyMwbgCuBlcCnquebHszM+ftMbcXLfkmq181smc0N+38L+K2e9aiB\nJXdJaubtByRpABUb7oW9+S1JR1Vx4R7eFlKSGhUX7pKkZoa7JA2ggsPdorsk1Sku3K24S1Kz4sJd\nktSs2HB3KqQk1Ssu3J0JKUnNigt3SVIzw12SBlCx4W7JXZLqFRfuPiBbkpoVF+6SpGbFhrtTISWp\nXnHh7lRISWrWGO4RcWNE7IuIR2r2R0R8MiL2RMTDEXFO77spSToS3Vy53wRsnGH/BcC66msL8Kdz\n75YkaS4awz0zvw48O0OTi4G/zJa7gRMi4uRedXCGfs33KSSpWL2oua8Gnmpbn6i2zQtL7pLUrBfh\n3ilvO15WR8SWiBiPiPHJyckenFqS1Ekvwn0COKVtfRTY26lhZm7LzLHMHBsZGZnTSS3KSFK9XoT7\nduDXqlkz5wHfz8yne3DczqzLSFKjRU0NIuIWYAOwKiImgKuAxQCZeQOwA7gQ2AP8CPiN+eqsJKk7\njeGemZsb9ifwH3rWI0nSnBX3CdUpzoSUpHrFhbt3hZSkZsWFuySpWbHhnk6GlKRaxYW7d4WUpGbF\nhbskqZnhLkkDqNxwt+QuSbWKC3dL7pLUrLhwlyQ1KzbcrcpIUr3iwj2cCylJjYoLd0lSM8NdkgZQ\nseHuXSElqV5x4W7JXZKaFRfukqRmXYV7RGyMiN0RsScitnbYf2pE3BkRD0TEwxFxYe+7+lreFVKS\n6jWGe0QMA9cDFwDrgc0RsX5as/8G3JaZZwObgE/1uqOH+zNfB5akAdLNlfu5wJ7MfCIzDwC3AhdP\na5PAG6vlNwF7e9dFSdKRanxANrAaeKptfQJ4x7Q2VwNfiojfA94AvKcnvZMkzUo3V+6dKiHTC96b\ngZsycxS4EPiriHjdsSNiS0SMR8T45OTkkfe2vQOW3CWpVjfhPgGc0rY+yuvLLpcBtwFk5v8BlgGr\nph8oM7dl5lhmjo2MjMyqw06FlKRm3YT7fcC6iFgbEUtovWG6fVqb/we8GyAifppWuM/t0lySNGuN\n4Z6ZB4HLgZ3AY7RmxeyKiGsi4qKq2e8Dvx0RDwG3AL+eOb+FE6syklSvmzdUycwdwI5p265sW34U\neGdvu1bHuowkNfETqpI0gAx3SRpAxYb7PJf0JaloxYW7UyElqVlx4S5JalZsuFuUkaR6xYW7VRlJ\nalZcuEuSmhnukjSAyg13i+6SVKu4cA/nQkpSo+LCXZLUrNhw9wHZklSvuHC3KCNJzYoLd0lSM8Nd\nkgZQseHuTSElqV5x4e5MSElq1lW4R8TGiNgdEXsiYmtNm1+JiEcjYldE/HVvuylJOhKNz1CNiGHg\neuB8YAK4LyK2V89NnWqzDrgCeGdmPhcRb56vDk+xLCNJ9bq5cj8X2JOZT2TmAeBW4OJpbX4buD4z\nnwPIzH297earwsmQktSom3BfDTzVtj5RbWv3U8BPRcT/ioi7I2JjpwNFxJaIGI+I8cnJydn1WJLU\nqJtw73SpPL0osghYB2wANgN/FhEnvO6HMrdl5lhmjo2MjBxpXyVJXeom3CeAU9rWR4G9Hdp8PjNf\nzswngd20wn7eWHKXpHrdhPt9wLqIWBsRS4BNwPZpbf4W+EWAiFhFq0zzRC87OsWpkJLUrDHcM/Mg\ncDmwE3gMuC0zd0XENRFxUdVsJ7A/Ih4F7gT+S2bun69OS5Jm1jgVEiAzdwA7pm27sm05gQ9VX0dF\nOhdSkmoV9wlVSVIzw12SBpDhLkkDqNhwt+IuSfWKC3enQkpSs+LCXZLUrNhwdyakJNUrLty9K6Qk\nNSsu3CVJzQx3SRpABYe7RXdJqlNcuDsVUpKaFRfukqRmxYa7UyElqV5x4W5ZRpKaFRfukqRmhrsk\nDaCuwj0iNkbE7ojYExFbZ2j3/ojIiBjrXRc7s+QuSfUawz0ihoHrgQuA9cDmiFjfod0K4D8C9/S6\nk685j7cfkKRG3Vy5nwvsycwnMvMAcCtwcYd2fwh8HHiph/2TJM1CN+G+GniqbX2i2nZYRJwNnJKZ\nd8x0oIjYEhHjETE+OTl5xJ1t51RISarXTbh3qoMcjtaIGAL+GPj9pgNl5rbMHMvMsZGRke572d4Z\nqzKS1KibcJ8ATmlbHwX2tq2vAH4WuCsivgOcB2w/Gm+qSpI66ybc7wPWRcTaiFgCbAK2T+3MzO9n\n5qrMXJOZa4C7gYsyc3xeeixJatQY7pl5ELgc2Ak8BtyWmbsi4pqIuGi+O1jbLydDSlKtRd00yswd\nwI5p266sabth7t2qZ8ldkpr5CVVJGkDFhrtTISWpXnHh7lRISWpWXLhLkpoZ7pI0gIoNd0vuklSv\nwHC36C5JTQoMd0lSk2LDPZ0LKUm1igt3p0JKUrPiwl2S1Mxwl6QBVFy4T1VlLLlLUr3iwn2oKrq/\nYrpLUq3iwn14qBXuh14x3CWpTnHhPlSFuxfuklSvvHCviu6HTHdJqtVVuEfExojYHRF7ImJrh/0f\niohHI+LhiPhqRPxE77vaMhyWZSSpSWO4R8QwcD1wAbAe2BwR66c1ewAYy8wzgNuBj/e6o1NeLcsY\n7pJUp5sr93OBPZn5RGYeAG4FLm5vkJl3ZuaPqtW7gdHedvNVQ165S1KjbsJ9NfBU2/pEta3OZcAX\nOu2IiC0RMR4R45OTk933ss3hsozZLkm1ugn3Tndz6RitEfGrwBhwXaf9mbktM8cyc2xkZKT7XrYZ\nGjp8rFn9vCQdCxZ10WYCOKVtfRTYO71RRLwH+CjwC5n549507/Usy0hSs26u3O8D1kXE2ohYAmwC\ntrc3iIizgU8DF2Xmvt5381WHP8Tklbsk1WoM98w8CFwO7AQeA27LzF0RcU1EXFQ1uw5YDnw2Ih6M\niO01h5t7h8MPMUlSk27KMmTmDmDHtG1Xti2/p8f9qnX4Q0yWZSSpVnGfUPXeMpLUrLhwjwginC0j\nSTMpLtyhVXf3DVVJqldkuA9HcOiVfvdCkhauIsN9aMiyjCTNpMxwj/ANVUmaQZHhPmzNXZJmVGS4\nDw2FH2KSpBmUGe7hPHdJmkmR4T48NMTBV5wuI0l1igz345cM8+KBQ/3uhiQtWMWG+w8Nd0mqVWS4\nv2HpIn504GC/uyFJC1aR4X78kmF++GOv3CWpTpHhvnzpIl74sVfuklSnyHB/y5uW8fTzL3oLAkmq\nUWS4n3rS8fzwwCGeeeFAv7siSQtSV+EeERsjYndE7ImIrR32L42Iz1T774mINb3uaLszRk8A4H98\ne3I+TyNJxWp8zF5EDAPXA+cDE8B9EbE9Mx9ta3YZ8FxmnhYRm4BrgQ/MR4f54X7Ofvp2Tn/zP+ea\nv9vFPz73Ij/z1jcysmIpb1g6zPFLFrF00RCLhoYYHg4WDbW+hoeCqJ6/KkmDrptnqJ4L7MnMJwAi\n4lbgYqA93C8Grq6Wbwf+JCIi56Eo/j8f2Ma1j/0Fi980zOgbgy8+CV98cvbHM+4lHW3nLTuTKz94\n87yeo5twXw081bY+Abyjrk1mHoyI7wMrgWfaG0XEFmALwKmnnjqrDq847b2cvvcb5IEXADiUyaFD\nyaFMMlv3eZ/6HyVbHXp1+dXVjvuOxOt+xvd2JXXpxGVvnvdzdBPunS5up0dZN23IzG3ANoCxsbFZ\nxeGZbzmHM39lx2x+VJKOGd28oToBnNK2PgrsrWsTEYuANwHP9qKDkqQj10243wesi4i1EbEE2ARs\nn9ZmO3Bptfx+4GvzUW+XJHWnsSxT1dAvB3YCw8CNmbkrIq4BxjNzO/DnwF9FxB5aV+yb5rPTkqSZ\ndVNzJzN3ADumbbuybfkl4Jd72zVJ0mwV+QlVSdLMDHdJGkCGuyQNIMNdkgZQ9GvGYkRMAv8wyx9f\nxbRPvx4DHPOxwTEfG+Yy5p/IzJGmRn0L97mIiPHMHOt3P44mx3xscMzHhqMxZssykjSADHdJGkCl\nhvu2fnegDxzzscExHxvmfcxF1twlSTMr9cpdkjQDw12SBlBx4d70sO6SRMSNEbEvIh5p23ZSRHw5\nIh6vvp9YbY+I+GQ17ocj4py2n7m0av94RFza6VwLQUScEhF3RsRjEbErIv5TtX2Qx7wsIu6NiIeq\nMf9BtX1t9TD5x6uHyy+pttc+bD4irqi2746I9/ZnRN2LiOGIeCAi7qjWB3rMEfGdiPhmRDwYEePV\ntv69tjOzmC9atxz+e+BtwBLgIWB9v/s1h/H8K+Ac4JG2bR8HtlbLW4Frq+ULgS/QeurVecA91faT\ngCeq7ydWyyf2e2w14z0ZOKdaXgF8G1g/4GMOYHm1vBi4pxrLbcCmavsNwO9Wy/8euKFa3gR8plpe\nX73elwJrq38Hw/0eX8PYPwT8NXBHtT7QYwa+A6yatq1vr+2+/0KO8Jf388DOtvUrgCv63a85jmnN\ntHDfDZxcLZ8M7K6WPw1snt4O2Ax8um37a9ot5C/g88D5x8qYgeOB/0vrGcTPAIuq7Ydf17Sem/Dz\n1fKiql1Mf623t1uIX7Se2PZV4F3AHdUYBn3MncK9b6/t0soynR7WvbpPfZkv/ywznwaovk89Sbdu\n7EX+Tqo/vc+mdSU70GOuyhMPAvuAL9O6An0+Mw9WTdr7/5qHzQNTD5svaszAJ4D/CrxSra9k8Mec\nwJci4v6I2FJt69tru6uHdSwgXT2Ie0DVjb2430lELAf+BvjPmfmDiE5DaDXtsK24MWfmIeCsiDgB\n+O/AT3dqVn0vfswR8T5gX2beHxEbpjZ3aDowY668MzP3RsSbgS9HxLdmaDvvYy7tyr2bh3WX7nsR\ncTJA9X1ftb1u7EX9TiJiMa1gvzkzP1dtHugxT8nM54G7aNVYT4jWw+Thtf2ve9h8SWN+J3BRRHwH\nuJVWaeYTDPaYycy91fd9tP4TP5c+vrZLC/duHtZduvaHjV9Kqy49tf3XqnfZzwO+X/2ZtxP4pYg4\nsXon/peqbQtOtC7R/xx4LDP/qG3XII95pLpiJyKOA94DPAbcSeth8vD6MXd62Px2YFM1s2QtsA64\n9+iM4shk5hWZOZqZa2j9G/1aZl7CAI85It4QESumlmm9Jh+hn6/tfr8JMYs3LS6kNcvi74GP9rs/\ncxzLLcDTwMu0/se+jFat8avA49X3k6q2AVxfjfubwFjbcX4T2FN9/Ua/xzXDeP8lrT8xHwYerL4u\nHPAxnwE8UI35EeDKavvbaAXVHuCzwNJq+7JqfU+1/21tx/po9bvYDVzQ77F1Of4NvDpbZmDHXI3t\noepr11Q29fO17e0HJGkAlVaWkSR1wXCXpAFkuEvSADLcJWkAGe6SNIAMd0kaQIa7JA2g/w+08WtZ\nxOHBcwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10a196090>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW4AAAEDCAYAAAAVyO4LAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzt3Xd4XFeZ+PHvO12jmVG3XBRbtuOS\n5pKIVAghIR0CDxtIoyxkf4altx8kyz5LYJ/9LW0puywsWWDZBZIQICGQhUAKIZQQY8cldty7XNS7\nNP38/rh3ZHWNpBlp7vj9PI8yM/feuXOOIr969d5zzxFjDEoppZzDNdcNUEopNTUauJVSymE0cCul\nlMNo4FZKKYfRwK2UUg6jgVsppRwmb4FbRL4rIs0isiNH53tCRDpF5PER25eKyAsisk9EfiQivlx8\nnlJKFap8ZtzfA27I4fm+CLxtjO2fB75ijFkBdAB35/AzlVKq4OQtcBtjngPah24TkeV25rxZRH4v\nIquncL6ngZ4R5xPgauAn9qb/Bt44s5YrpVRh88zy590PvMcYs09ELgG+gRV4p6sK6DTGJO3XjcCi\nGbZRKaUK2qwFbhEJAZcDP7YSZQD89r43AZ8d423HjTHXT3TaMbbpPfxKqaI2mxm3Cys7XjdyhzHm\nEeCRaZyzFSgXEY+dddcBJ2bWTKWUKmyzNhzQGNMNHBKRN4NVnxaRtTM8pwF+C9xqb3oH8NiMGqqU\nUgVO8jU7oIg8CFwFVANNwKeBZ4BvAgsAL/CQMWasEslY5/s9sBoIAW3A3caYX4vIMuAhoBLYArzV\nGBPLbW+UUqpw5C1wK6WUyg+9c1IppRwmLxcnq6urTX19fT5OrZRSRWnz5s2txpiabI7NS+Cur69n\n06ZN+Ti1UkoVJRE5ku2xWipRSimH0cCtlFIOo4FbKaUcZrbnKlFKqRlLJBI0NjYSjUbnuilTFggE\nqKurw+v1TvscGriVUo7T2NhIOBymvr6eIXMfFTxjDG1tbTQ2NrJ06dJpnyerUomIfEREdorIDhF5\nUEQC0/5EpZSaoWg0SlVVlaOCNoCIUFVVNeO/FCYN3CKyCPgg0GCMOR9wA7fP6FOVUmqGnBa0M3LR\n7mwvTnqAEhHxAEEKaAa+rcc6eWZ301w3QymlZs2kNW5jzHER+RJwFBgAfmOM+c3I40RkA7ABYPHi\nxblu55i6owne9I0/kjbw8/dfwZq68ln5XKXUma2trY1rrrkGgFOnTuF2u6msrMTtdhOPx2lvb8ft\ndlNTY90IuXHjRny+3C2HO2ngFpEK4A3AUqATayGEtxpjfjD0OGPM/Vgr3NDQ0DArM1f9+UAbafuT\nnthxSgO3UmpWVFVVsXXrVgDuu+8+QqEQH//4xwf3j7Utl7IplbwWOGSMaTHGJLAWPLg8L62Zos1H\nOvC5XayYF2Lrsc65bo5SSs2KbIYDHgUuFZEgVqnkGqAgJiLZdaqHlfNDnL+wjCdf1jq3Umeiz/xi\nJy+f6M7pOc9dGOHTrz8vp+fMpUkzbmPMC1irqL8IvGS/5/48tysrR9v6WFJVSn11KW19cbqjiblu\nklJK5V1WN+AYYz6NtYJNwUim0jR2DHDjBQuorwoCcKS1nwvqyua4ZUqp2VTImXG+OHaukpNdUZJp\nw5LKIEuqSgE40t43x61SSqn8c2zgPtbeD8DiyiALy0oAaOrWpSaVUsXPsXOVNPdYQbq2LECkxIPf\n46K523kTziilnO2+++7LalsuOTbjbrEDd03Yj4hQGwnQpIFbKXUGcG7g7o3h97gI+60/Gmojfi2V\nKKXOCM4N3D2xwWwbYF4kQFOPZtxKqeLn2MDd2hujOuQffD0v7KdZM26l1BnAsYE7k3FnVJX66I0l\niSVTc9gqpZTKv6IJ3BWl1sxbnf1696RSqrg5cjhgIpWmvT9OzZBSSUXQCtwd/XFqI7pAj1Iqf8aa\n1jUzheu2bdtYu3bt4LG3334799xzT04/35GBu7M/gTFQFTo9v20mcLf3xeeqWUqpM8RE07qGQqHB\nffniyFJJ14AVnMtKTq+SXGmXSjr6tFSilCpujsy4uwas4FweHJJxl1pBvL1fM26lzii/ugdOvZTb\nc86/AG783LTeOjAwwLp16wZf33vvvdx22225ahng0MCduQBZPiTjLi/JZNwauJVSc6ekpCTvpRJn\nB+7g6cDts++i7NCMW6kzyzQzYyebtMYtIqtEZOuQr24R+fBsNG48mVLJ0Bo3WEMCNeNWShW7bFZ5\n3wOsAxARN3AceDTP7ZpQ50ACEQgHRgfuNg3cSqk5NLLGfcMNN/C5z+X2r4KplkquAQ4YY47ktBVT\n1NUfJxLw4nbJsO1lJd7BbFwppWbDyClcU6n837091eGAtwMPjrVDRDaIyCYR2dTS0jLzlk2gcyAx\nqkwCEAl46NHArZQqclkHbhHxAbcAPx5rvzHmfmNMgzGmIXMHUb50DSSGXZjMiJR46Y4m8/rZSik1\n16aScd8IvGiMacpXY7LV2T92xh0OeHSld6VU0ZtK4L6Dccoks83KuH2jtkcCXuLJNNGEzhColCpe\nWQVuEQkC1wKP5Lc52ekaSFBWMvq6asTOwjXrVkoVs6wCtzGm3xhTZYzpyneDJpNOGzr744N3Sg4V\nCVjBvHtA69xKqeLluDsne+NJ0oZxL06CZtxKqfxzu91ccMEFg68z07deddVVnDx5kpKSEgDOPvts\nfvKTn+T0sx0XuLvt4X7hwBilEvuGnB4dWaKUyrOJ5iT54Q9/SENDQ94+23HTuvbGrKA88q5JYLDu\n3a1juZVSRcxxGXevnU2H/KObngnmWipR6szx+Y2fZ3f77pyec3Xlaj558ScnPGai6VvvuuuuwVLJ\ntddeyxe/+MWcts9xgbvHzrhDE5RK9OKkUirf5rJU4rjAncm4w2Nk3AGvC69bNONW6gwyWWZcjIqq\nxi0iRAJerXErpYqaYzPusUolYA0J1FElSql8m2j61qE17urqap566qmcfrbjAndPLIkIBL3uMfdH\ndL4SpdQsGG/61meffTbvn+28Ukk0ScjnwTViLu6MsJZKlFJFznmBO5YYt0wCECnx6NSuSqmi5sDA\nnRxzDHdGyO8ZrIMrpYqXMWaumzAtuWi34wJ3TzQ5YcYd8nvpi2ngVqqYBQIB2traHBe8jTG0tbUR\nCARmdB7HXZycNOMOeKyJqNJm3Dq4UsrZ6urqaGxsJN/LJOZDIBCgrq5uRudwXuCOJllQNv5vq7Df\ngzHQn0hNGOCVUs7l9XpZunTpXDdjzmS7kEK5iPxERHaLyC4RuSzfDRvPZBl3qb1P69xKqWKVbUr6\nNeAJY8yt9qLBwTy2aUK90SQh/+i7JjMy9e/eWAKYWR1JKaUK0aSBW0QiwJXAXwMYY+JAPL/NGls6\nbeiNT3xxMjOHid49qZQqVtmUSpYBLcB/icgWEfm2iJSOPEhENojIJhHZlK8LBn3xJMaMPcFURiao\n98V0wWClVHHKJnB7gAuBbxpj1gN9wD0jDzLG3G+MaTDGNNTU1OS4mZbeCaZ0zcjUv61SiVJKFZ9s\nAncj0GiMecF+/ROsQD7rJlpEISOkpRKlVJGbNHAbY04Bx0Rklb3pGuDlvLZqHBMtopBxOuPWwK2U\nKk7Zjir5APBDe0TJQeCd+WvS+DIZd2SCwK3DAZVSxS6rwG2M2Qrkbx2eLA3WuCcYDujzuPB7XPTG\nNXArpYqTo+YqmWwRhQydaEopVcwcFbgHa9yT3MoeCni0xq2UKlqOCtzZjCrJ7NeMWylVrJwVuGMJ\ngj437klm/Qv5PYPZuVJKFRuHBe6JJ5jKCAc8Oie3UqpoOSpwT7aIQkapX2vcSqni5ajA3RtLTjhP\nSYbWuJVSxcxZgTvLjDsU0Bq3Uqp4OStwZ1vj9nuIJ9PEkjpDoFKq+DgqcPdMsohCRia469SuSqli\n5KjA3RtLEs7y4iSgI0uUUkXJMYHbGDOl4YCgU7sqpYqTYwL3QCJFKm2yuzhpl1N0SKBSqhg5JnBn\ne7s7jFwwWCmliotjAndmeF82NW5dBUcpVcyymo9bRA4DPUAKSBpjZn1u7kzGPZXAraNKlFLFKNsV\ncABeY4xpzVtLJpHNIgoZWipRShUz55RKplDjDnrdiOjyZUqp4pRt4DbAb0Rks4hsGOsAEdkgIptE\nZFNLS0vuWmjrzXIRBQCXSwj59LZ3pVRxyjZwX2GMuRC4EXifiFw58gBjzP3GmAZjTENNTU1OGwnQ\nG7XKHtnUuMGeIVAzbqVUEcoqcBtjTtiPzcCjwMX5bNRYMhl3aRYZN+jyZUqp4jVp4BaRUhEJZ54D\n1wE78t2wkXpiSXweFz5Pdn8khHRObqVUkcomfa0FHhWRzPEPGGOeyGurxtAbzW4u7oywZtxKqSI1\naSQ0xhwE1s5CWybUF8tuLu6MkN/Dqa5oHluklFJzwzHDAbOdYCpDly9TShUrxwRuay7uqWXcOqpE\nKVWMHBO4s52LOyMc8NAbT5JOmzy2SimlZp+jAne2QwHBCtzGQF9cs26lVHFxTuCeYqkkHLDmNOnW\ncolSqsg4JnD3THFUScQO3D1RnWhKKVVcHBG4Y8kU8WR6yuO4QefkVkoVH0cE7sy82lMrlWQCt2bc\nSqni4ojAPbhsWWDyubgzwoOlEs24lVLFxRmBewpTumZE7IxbL04qpYqNowL31MZx68VJpVRxckjg\ntoLvVMZxB7wuPC7RUolSqug4InBPZdmyDBEhUuLVjFspVXQcEbinUyrJHK8Zt1Kq2DgjcE8j4wYr\ncHcPaMatlCouWQduEXGLyBYReTyfDRpLbyyJCAR97im9L+z3asatlCo6U8m4PwTsyldDJpKZ0tVe\nhSdrWipRShWjrAK3iNQBNwPfzm9zxtYbm9qyZRnhgF6cVEoVn2wz7q8CnwDSeWzLuHqjU5vSNUMz\nbqVUMcpmlffXAc3GmM2THLdBRDaJyKaWlpacNRDsZcumOKIErLsndTEFpVSxySbjvgK4RUQOAw8B\nV4vID0YeZIy53xjTYIxpqKmpyWkjp7reZEakxIsx0KuLKSilisikgdsYc68xps4YUw/cDjxjjHlr\n3ls2xFSXLcvQqV2VUsXIMeO4p5Nx63wlSqliNKVoaIx5Fng2Ly2ZgFUqyX5K14xMxt09oBm3Uqp4\nFHzGnU4bO3BP7eYb0IxbKVWcCj5w99jzlERKpp9xa41bKVVMCj5wZ+YamU7gnuqCwU/vauLTj+3g\nROfAlD9LKaVmy9Sv+M2ybjvoRqawbFlGmR3sO/snD9x7m3p4zw82k0gZthzr5GfvvQKXa2q32Btj\n+ObvDrDzeDefvGE1i6uCU26zUkpNxgEZd6ZUMvXfMT6Pi1Kfm44sAvd//fEwHpeLv7tpNdsbu3j+\nYNuUP++nLx7nC0/s4X9fOsl7H9isN/4opfKi8AP3DDJugPKgj87++ITHJFJpHt92ghsvmM/bL6sn\nEvDw0xcbp/Q56bThm8/u5/xFEb5w6xp2HO/mz9MI/kopNZnCD9x2jbtsGjVugIpSLx2TBO7tjV30\nxJK89pxaAl43r1k9j+f2tkwpY95yrJMDLX2847J6blm7kEjAw6Nbjk+rzUopNZHCD9z2iJBpZ9wl\nvklLJc8faAXg0mVVAFy5oobW3jgvn+zO+nOe3tWE2yVcd+58Al43r1pRw3P7WjBGyyVKqdwq/MBt\nZ9zTmWQKoDzopWuSVXD+criDVbVhKkt9AFx+thXANx5qz/pzntndzCvqKygLWr9grlxZTVN3jL1N\nvdNqt1JKjafwA3c0QdjvwT3FER4ZFUHfpKWSPad6OG9hZPD1grISFpQF2HKsM6vP6OiLs/tUD69a\ncXpyrUz2/uLRjmm0Wimlxlf4gXsgOa0x3BkVdsadGqde3dWf4FR3lJXzw8O2r19czpYsg+7Wxs7B\n92QsrgxSHvSyLcvgr5RS2Sr8wB1NTGtmwIzyoA9jGHfR4D1NPQCsGhm4z6qgsWOA5p7opJ+x9Wgn\nLoE1dacDt4iwpq6cbY1d0267UkqNpfAD90BiZhl3qfXe8colg4G7dnTGDVZQnsyWY52srA2PmsFw\nbV0Ze5t6iCZSU263UkqNp/ADdzQ57RElYGXcwLgjS/ac6iYc8LCgLDBs+3kLy3C7hO2TZMzptGHb\nsc5hZZKM1fMjpNKG/c16gVIplTuFH7gHEtO6azKjwg7c492Es/dUL6tqw6NWkC/xuVlZG2Zb48QZ\n96G2ProGEqw7a3TgXlkbAtDArZTKqcIP3NHEjDLuimCmVDI64zbGsPtU96j6dsbaujK2N3ZNOBY7\nU0pZv7hi1L766lI8LmGvXY5RSqlcyGax4ICIbBSRbSKyU0Q+MxsNg9Nzcc+kxp0Zm93aGxu1r6k7\nRnc0OX7gPqucroEER9r6xz3/lmMdhP0ezq4JjdrndbtYWl2qY7mVUjmVTcYdA642xqwF1gE3iMil\n+W2WpSeWxBhrtfbpCvk9BH1umrtHB+7dp6w7I0demMxYU1cGMGG5ZOuxTtacVTbuTIIra8Psb9aM\nWymVO9ksFmyMMZmU0Wt/zcp93DOZiztDRJgX9o85rC9Twlg5TuBeWRvG73GNe4FyIJ5i18ke1p81\nukySsaI2xJH2fh1ZopTKmaxq3CLiFpGtQDPwpDHmhTGO2SAim0RkU0tLS04al7lVfSYZN8C8cIDm\nnrEy7h7mhf1U2OWUkbxuF+cvKmP7OBn3jhNdpNJmzAuTGSvmhTEGDrRouUQplRtZBW5jTMoYsw6o\nAy4WkfPHOOZ+Y0yDMaahpqZm9EmmoWtwZsCxA2u2aiJ+WsYI3Hubesatb2esqSvjpeNdJFPpUfsy\nd1auG2MoYMbS6lIADrX2TaXJSik1rimNKjHGdGKt8n5DXlozQuammcxNNNM1L+ynuXt4qSSVNuxr\n6h23vp2xtq6caCLNvjGG9L14pJMlVUGqQ/5x319fba2CM9EFTqWUmopsRpXUiEi5/bwEeC2wO98N\ng9ND+DJjsadrXjhAXzxFX+z0osFH2vqIJdNZZdzAqHKJMYYXj3awfoIyCUDQ56E24teMWymVM9lk\n3AuA34rIduAvWDXux/PbLEtnn5VxlwdnnnEDw+rce06NPUfJSPVVpUQCnlFzjpzsitLcExtz/PZY\n5zisgVsplSOTXvUzxmwH1s9CW0bp6E8Q9Lnxe9wzOs+Ccut29hOdA4M1512nenDJ+CNKMlwua7Ko\nkRn35iN2fXuSjBusOvdTu5qm03SllBqloO+c7OyPz7hMAtYUqzC8zrzrZDdLq0sJeCf/pbCmrozd\nJ4dPFvWnA62EA55h83iPZ0lVKa29cXqiky9arJRSkynowN3RH5/xhUmwFkbwuV0caT9drth1sptz\nFkwedMHKqpP2ZFJg1bef29vKZcuq8Lgn/xYutS9QHm7VC5RKqZkr8MCdyEnG7XYJdZUlHLUz7q6B\nBI0dA1kH7suWV+Fzu3jyZavcsb+5l+OdA7xqRXVW76/PDAls0zq3UmrmCjpwd/bHB6dlnakllcHB\nUsnLJ6xb3c/NoswBEA54eeWKan614xTptOGxrSdwCVx/3vwsP9sK3HqBUimVCwUduK2Me+alEoAV\ntWH2t/QST6bZfMRaBPjCCW5VH+mN6xdxvHOAH206xgMbj/LqlTXMiwQmfyPWFLELygIauJVSOVGw\ngTuVNnRHEznLuC9YVEY8mWZvUw8b7VXdy6bwS+HmCxawen6Yex95ie6BBB+5duWUPr++qpTDWipR\nSuXAzCYByaOugQTGkLOMOzNs73d7W3jhYBt3XLx4Su93u4T/ftfFfO9Ph7l8edWw9SWzUV8d5Nc7\ndUigUmrmCjZwD97unqOMu66ihFW1Yb746z0A3LxmwZTPURsJ8MkbVk/r8+urSmnvi9M1kKBsBrMd\nKqVUwZZKMkuNzfSuyQwR4f9cuQyAK86uomFJ9vXtXMiMLNE6t1Jqpgo34+7LzTwlQ916UR2XLK2k\nNhIYtcZkvg2dJXBtFndbKqXUeAo3cOe4VJJxln0X5WxbUhXEJXBQ5+VWSs1QwZZK2vtyM6VrofB7\n3CyuDHKgRUslSqmZKdjA3dobw+9xEfIX7B8FU7a8JqQr4SilZqyAA3ec6pB/1mvR+bR8XoiDrX2k\n0rOyZKdSqkgVcOCOUR0ef2UZJ1peU0o8meZ4x8BcN0Up5WDZrIBzloj8VkR2ichOEfnQbDSspSdG\nTSi3Fybn2rKaEKALByulZiabjDsJfMwYcw5wKfA+ETk3v806XSopJss1cCulcmDSwG2MOWmMedF+\n3gPsAhbls1GptKG9L1Z0gbuy1EdF0KuBWyk1I1OqcYtIPdYyZi+MsW+DiGwSkU0tLS0zalRHf5y0\ngeoiK5VAZmSJDglUSk1f1oFbRELAT4EPG2O6R+43xtxvjGkwxjTU1NTMqFGtvdaivsV2cRKswL2v\nqQdjdGSJUmp6sgrcIuLFCto/NMY8kt8mQWuPdfNNsZVKwFq8oaM/wanu6Fw3RSnlUNmMKhHgO8Au\nY8yX89+kIRl3EQbu8xdZq+7sOD7qjxallMpKNhn3FcDbgKtFZKv9dVM+G5UJ3DVFGLjPWRBBBHYc\n75r02K7+BJsOtw9bXV4ppSa9n9wY8wdgVm9fbOmN4XO7iJQUz+3uGUGfh+U1IXaemDjjfurlJj74\n0Bb64ynqKkr4/t2XDM4wqJQ6sxXknZOtPXGqQr6iut19qPMXRnjpeOe4Fyhbe2N85OGtLKsp5au3\nraM/nuI9399MLKmZt1KqQAN3U3c064V4nejCJRU0dcc41j72re9feGI3A/EUX71tPW9cv4gvvXkN\ne5p6+P7zR2a5pUqpQlSQgftUd5QFRRy4L19eBcDzB1tH7dtytIOHNzVy9yuXcvY8607Lq1fXcvny\nKv7jdwe13q2UKtDA3RVlflnxBu7lNSFqwn7+dKBt2PZU2vAPj+1kXtjPB65ZMWzf+19zNq29MX6x\n7cRsNlUpVYAKLnD3RBP0xpJFHbhFhFetqObZPS3D6tYPbzrGS8e7+NTN54yah/yy5VUsrynlgY1H\nZ7u5SqkCU3CBu8m+MWVBEQdugDesW0TXQILf7m4GrBV/vvDEbi6ur+SWtQtHHS8i3HnJErYc7WTX\nSR0DrtSZrKAC9z/978uDF+Bqi7jGDXDF8ioWlZfwjWcPkEyl+fufvURvLMln33jeuKNp3rR+ET63\nix9vapzl1iqlCklBBe4HNx7jZ1utGm6xZ9wet4tP3LCK7Y1dXP65Z/jlS6f46LWrWD0/Mu57Kkp9\nXHtuLT/bepx4Mj2LrVVKFZKCCtwlPjddAwmg+DNusMol//jG81lWU8p9rz+X97x62aTvubWhjva+\nOM/sbpqFFiqlClFB3ZoY9LkBa97qgNc9x62ZHW+7dAlvu3RJ1sdfuaKG2oifhzc1csP5C/LYMqVU\noSqojDvos36PLK4MznFLCpfbJfzVhXU8u6eZZp1hUKkzUoEFbivLXlKlgXsit15UR9rAI1uOz3VT\nlFJzoKACt8dljaZYohn3hJbVhHhFfQUPbzqW9YIMqbQu3KBUsSiowJ0Zw720RmfBm8ybLzqLgy19\nvHi0c8LjfrPzFNd95Xcs/7tfcv1XnuOP+0ffZq+UcpaCCtwN9ZUAXLS4co5bUvhuWrOAoM/Nw385\nNu4x3/vjITZ8fzOC8N6rlpNIpXn7dzfy5Ms6IkUpJ8tmBZzvikiziOzId2M+/fpzeeS9l7NYa9yT\nCvk9vHH9Ih7Z0six9v5R+7/zh0Pc94uXuf68Wn7xgVfyiRtW89j7r+C8hRE+/NCWMd+jlHKGbDLu\n7wE35LkdAIQDXi5cXDEbH1UUPnj1Ctwu4bOPvzys1v2t3x3gHx9/mRvPn8/X77wQn8f63xwOePnm\nWy/CAP/w2A5dsFgph5o0cBtjngPaZ6EtaormlwX42LWrePLlJv7fL3dxqLWPTz+2g3/+1W5et2YB\n/3rHerzu4f+LF5WX8NFrV/LbPS08tat5jlqulJoJySbrEpF64HFjzPkTHLMB2ACwePHii44c0Un/\nZ4Mxhk/9bAcPvHB61sB3XbGUT918Dm7X2HOeJFNprv3Kc/g9Ln75wVfhGuc4pdTsEZHNxpiGrI7N\nVeAeqqGhwWzatCmbQ1UOGGN48WgHe0710lBfwcra8KTveWzrcT700Fa+fud6Xrdm9GyESqnZNZXA\nXVCjStT0iAgXLankzksWZxW0AV63ZiEra0N8+cm9JFM6YZVSTqKB+wzldgkfvXYlB1v6eGyrrqqj\nlJNkMxzwQeB5YJWINIrI3flvlpoN1583n/MXRfjq03t1mlilHCSbUSV3GGMWGGO8xpg6Y8x3ZqNh\nKv9EhI9dt4pj7QM8vGn8G3mUUoVFSyVnuKtW1nDRkgr+7Zl9uoK8Ug6hgfsMJyJ8/LpVNHXH+Pbv\nD851c5RSWdDArbhseRU3XTCff316P/ube+a6OUqpSWjgVgB85pbzCfrdvOcHLw4uHzeWtt4Y+5p6\n6OiLz2LrlFJDFdTSZWru1IT9fOOuC3nHdzdy27ee51/espbzFpYRT6bZ1tjJs3uaeWZ3C7tOdg++\n59wFEd71yqW8af0ivftSqVmU1Z2TU6V3TjrXc3tb+MCDW+gaSFAR9NIXTxFPpnG7hIYlFVy1ah6L\nKko43jHAz7edYNfJbi5aUsHXbl9HXYXO6qjUdOX8lvep0sDtbF39CR7d0sje5l5Cfg8XLi7nsuXV\nlJV4hx2XThse2XKcz/x8Jy6X8OW3rOWac2rnqNVKOZsGbjWrDrf28b4HXmTniW7e/epl/N/rVuFx\n6+UTpaZC5ypRs6q+upSf/u3l3HnJYr71u4Pc+Z8vcKJzYK6bpVTR0oxb5dSjWxr5u0d2YDC87dIl\nvGHdIlbND+N1u+iPJzneMcCJrigtPTE6+uIEfG4WlQdYW1dOVcg/181Xas5oqUTNqWPt/XzpN3t4\nfPtJUmmDS8DjchGfZBbCS5dV8uaLzuL1axcOrtqj1JlCA7cqCK29MZ7b28Kh1j4SKUM44KGuooSF\n5SXMC/spD/qIJVIcau3j+YNtPLb1BIda+1hYFmDDlcu47RWLKfG557obSs0KDdzKkYwxPLu3hX9/\nZj+bjnRQEfRy1yVLePtlS5iz7Z7mAAANxklEQVQXCcx185TKqzM7cHefhG0PQPtBqFwGa++EyIK5\naYuato2H2vnP3x/kqV1NuMQaQ/6qFdWsrA0zLxIg6HPjc7vwuAWv24XHJXjcrsFtHpcgojcFKec4\ncwP3pv+CJ+6F5ACEaqG3CbylcPOXYN2ds98eNWNH2vr48aZGntrVxO5T2c+jUlbiZcW8ECtqw6yq\nDbGyNsyK2jDVIZ8GdFWQ8rHm5A3A1wA38G1jzOcmOn7WA3c6DU/fB3/8Giy/Bm7+F6hcamXdP/8g\nHP49vPY+eOVHZq9Ng21LwUs/hpcfg4FOmH8+XPROzLxzSKQTpEyKgDugwSQLPdEEB1r6aO+L0W/f\n0ZlMGRJp+zGVJpk2JJJpTnZH2dfUw96m3mFzr1SW+qivClJZ6qeq1EfA60JEcLusL4/LyuC9diYf\nDnipCvmoDvmpCfmpDvsI+nSmCJV7OQ3cIuIG9gLXAo3AX4A7jDEvj/eeWQ3ciSj87G9h5yPQ8C64\n8Yvgtv5hGWMYiHXT+Yv307XvV3Re8m7M6psp9ZUS8oaI+CJE/BH87hHD0JIxEtsfovHIcxw2MY6W\nL+RoSYijPcdoi7aRSCVwi5vKkkqqAlXML51vfQWtx9rSWioDlaS6jtP503dyuGU7hyK1HCwJcTjW\nxlGPmxNeLyms773X5aWqpIr6SD3Ly5ezrGwZiyOLWVi6kPml8/G5fVZfkgO0R9tpj7bTOtBK60Ar\nbQNttAy00Jvoxevy4nP7qApUMS84j3nBedSU1FBVUkVloBKf2ze8m+kk3fFuumPd9CX6GEgOMJAc\nIJaK4XF5CHgCBNwBgt4gIW+IUm8pQW8Qr2v4HZSFzBhDS0+MvU297GnqYV9TD0fb+2nvi9PRHyeW\nTJNKG4yBpP0LIJme+N9E0Oe2AnnYT7Ud1IM+Nz6PC5/bevTa5Rqvx4XXZZdv3C68dknH4xZ8bhd+\nj8t6n8cq82Se+91uvB7BJZkvawrezKMqPrkO3JcB9xljrrdf3wtgjPnn8d4z3cD9Lw9eTxRIi5AS\nIeUNkg5WkSwpJ40hZVKk0inSJk3SJEknoqRObScV6yZdsYRkqJa0SRNNRemMddIV6yKRHn+mu4yA\nO0DEHyHii0BigO6eE7RLmuSQfyBlRlhcsYLq8CL8bj/JdHIwgJ7qO0U8PXy2PLe4SJnhw9+CniBL\nQotY0t9FXfM+gpHFuNa8mS5StPa3cqjrEAe6DjCQHH7zimD94x15PmsfVPjKCPnLSKaTg31Pj3Fs\nyBvC4/JAKkEiGaWP6S2cEHAHKPWWUuoNUioeTDJKPDlAPJ0knk6SwJASSAMpA2kMIoLfHcDv9uP3\n+PG7/QTs1x4MpJNIOoWkE0gqBekEJpXApBMYcYG4weXGuNwYxN4miLitwIbBZRh8FGNwmTRuQEwa\nlzEIp3/WZfBRhjy3XhuM9cpAyhhSKUMybQX4ZNqQTBnr+eB2SBsr+BvMkDOObeTe06/N6ddj/LMc\n/T77k2T4MeNuy/x3xL7Mk6Hnl+F7B98z3vFjHsPw7+/INlnnEMbYPaqdUzHR22Ssb2wOBdwlfP5d\nv5jWe6cSuLP5m28RMHRdq0bgkjE+dAOwAWDx4sXZfPZw6TS/jJ0kZv9jc9v/8Dzt4HK5cftCuH1h\nXL4gbgR3rAd390ncJo1r3mo8pfPwixuXy0Wtu5a1/rWU+cso85dR7i+nzFNK2R+/jvvYn+m97L30\nLrmU7lg33fFuumJddEfb6Tq+EVf7IcLuAFXLXsvSpVdTH17CksYtlD35aWjaDG/9FCy6cFjTjTF0\nxDo41XfK+uo4QMvGb+If6KGs4W4WL7mK5eXLqQ3Wns6Wtj9slXHavwlv+R8462Lr22DSNPU10djb\nyPHe4zR1HiZx5A+kT7xIOJmkonIZleVnU42b6pZ9VDS+iNd1As65BS55N5x1CUmTom2gjeb+Zisz\nj7bS1nGIjsY/k27ZA4k+PC4fkcgiIqXziXhKCMX7Kek4RqDtAIF0kiQuolVLic5bTX+wgj63m95E\nP33RDvr6munt6aAv3kifgMsYvIDfGHwGvCaN2w6ibsBlrDgUdXuIe0uIut3E0imiJknMpIgLZMLq\nyH9WIi7EpBHM8IBgrH+gKYEkkEZIA2mxnht7nxmyD8DI6c8YGsqHfu7QYxhrP1j3HLvGjK9DjpNh\nr8YKG2aMKDPWOUe3J7uolt8wNTufORd9mK5wenb+Gsom434zcL0x5m/s128DLjbGfGC89+SsVNLX\nBoeehR2PwP6nIBkdvr/uYnj916D23OzOl4zBQ3fB/ifhonfCqz4GwUrY/b/w9Geh6xg03G3VwwOR\n4e9tPwj/8wbo74C7fgxLLhv7M/rb4b9vgbZ91nFLrxy/PSe3w8Nvg67jcNU9cNn7wFti7Yt2w8Zv\nwR//DWLdsOY2eM29UFE/ol2HYON/wpYfQKwL5l8A57zB+p6IG9oPwP6n4eCzYNJw9mth/V2w8kbw\njjHELt4HxzfDkefh6PPQ+BeI957eL27r+kHVCqhZBTWroWYlVC4HfwRcLkgMWPX8/jboa4beFuhr\nOf08OWBdNPYFIVAOJRXjfJWD24tVx4hCrBfiPdb/x2TUekwMgLjA47e+3JlH34hHv9U2pQqUY0sl\nE4r3WxcZW/eBywNnvQIWXjj1v6eScXjms/CnrzPsd/m886zRJ0suH/+9Xcet4N11DG79Lqy+efj+\nvjZrf+teuONBOPuaydsz0AGPvR92Pw6BMuuXUSoGxzZawWnVTXD130PteROfJ95nZfGbvwcntw3v\nW9UKOOd1cNFfjw78kzEGYj3WLw9vEPxhK5gqpXIq14Hbg3Vx8hrgONbFyTuNMTvHe48jbsBpPwh7\nf20FvEUXWZmxK4u79Ppa4YG3wIktcOUn4IoPWZnjkefh0XdbQxBvfyC7oD3U4T/A1gegaYf1i2lR\nA6y7Axaun3rfBjqg47A1oqXsLAjrVKtKFbp8DAe8Cfgq1nDA7xpj/mmi4x0RuGci3ge/+JA1zM8b\ntDLlnpMQWQRv/t5gvVoppbKV64uTGGN+CfxyRq0qJr5S+Ktvwyv+BnY+CtEuK2tfd6e1Tyml8kjv\nJJiJxZdaX0opNYv0MrtSSjmMBm6llHIYDdxKKeUwGriVUsphNHArpZTDaOBWSimH0cCtlFIOo4Fb\nKaUcJi9Ll4lIC3Bkmm+vBlpz2BwnONP6fKb1F7TPZ4qZ9HmJMaYmmwPzErhnQkQ2ZXu/frE40/p8\npvUXtM9nitnqs5ZKlFLKYTRwK6WUwxRi4L5/rhswB860Pp9p/QXt85liVvpccDVupZRSEyvEjFsp\npdQENHArpZTDFEzgFpEbRGSPiOwXkXvmuj0zISLfFZFmEdkxZFuliDwpIvvsxwp7u4jIv9r93i4i\nFw55zzvs4/eJyDvmoi/ZEpGzROS3IrJLRHaKyIfs7UXbbxEJiMhGEdlm9/kz9valIvKC3f4fiYjP\n3u63X++399cPOde99vY9InL93PQoOyLiFpEtIvK4/brY+3tYRF4Ska0issneNrc/18aYOf/CWsvy\nALAM8AHbgHPnul0z6M+VwIXAjiHbvgDcYz+/B/i8/fwm4FeAAJcCL9jbK4GD9mOF/bxirvs2QZ8X\nABfaz8NYC0yfW8z9ttsesp97gRfsvjwM3G5v/w/gb+3n7wX+w35+O/Aj+/m59s+8H1hq/1twz3X/\nJuj3R4EHgMft18Xe38NA9Yhtc/pzPeffFLtTlwG/HvL6XuDeuW7XDPtUPyJw7wEW2M8XAHvs598C\n7hh5HHAH8K0h24cdV+hfwGPAtWdKv4Eg8CJwCdadcx57++DPNvBr4DL7ucc+Tkb+vA89rtC+gDrg\naeBq4HG7/UXbX7t9YwXuOf25LpRSySLg2JDXjfa2YlJrjDkJYD/Os7eP13fHfk/sP4nXY2WgRd1v\nu2ywFWgGnsTKHjuNMUn7kKHtH+ybvb8LqMJZff4q8Akgbb+uorj7C2CA34jIZhHZYG+b05/rQlks\nWMbYdqaMUxyv7478nohICPgp8GFjTLfIWN2wDh1jm+P6bYxJAetEpBx4FDhnrMPsR0f3WUReBzQb\nYzaLyFWZzWMcWhT9HeIKY8wJEZkHPCkiuyc4dlb6XCgZdyNw1pDXdcCJOWpLvjSJyAIA+7HZ3j5e\n3x33PRERL1bQ/qEx5hF7c9H3G8AY0wk8i1XXLBeRTFI0tP2DfbP3lwHtOKfPVwC3iMhh4CGscslX\nKd7+AmCMOWE/NmP9cr6YOf65LpTA/RdghX112od1IePnc9ymXPs5kLmS/A6sGnBm+9vtq9GXAl32\nn16/Bq4TkQr7ivV19raCJFZq/R1glzHmy0N2FW2/RaTGzrQRkRLgtcAu4LfArfZhI/uc+V7cCjxj\nrILnz4Hb7VEYS4EVwMbZ6UX2jDH3GmPqjDH1WP9GnzHG3EWR9hdAREpFJJx5jvXzuIO5/rme68L/\nkGL9TVgjEQ4An5rr9sywLw8CJ4EE1m/au7Fqe08D++zHSvtYAf7d7vdLQMOQ87wL2G9/vXOu+zVJ\nn1+J9affdmCr/XVTMfcbWANssfu8A/gHe/syrEC0H/gx4Le3B+zX++39y4ac61P292IPcONc9y2L\nvl/F6VElRdtfu2/b7K+dmdg01z/Xesu7Uko5TKGUSpRSSmVJA7dSSjmMBm6llHIYDdxKKeUwGriV\nUsphNHArpZTDaOBWSimH+f+AH9DnI7ibjwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10aa563d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot them\n",
    "ellsP = Fobj.ellsP\n",
    "\n",
    "plt.plot(ellsP,crossClsMod[0,0],label='TT')\n",
    "plt.plot(ellsP,crossClsMod[0,1],label='TE')\n",
    "plt.plot(ellsP,crossClsMod[1,1],label='EE')\n",
    "plt.legend()\n",
    "plt.show()\n",
    "\n",
    "plt.plot(ellsP,crossClsMod[0,0]*DlFac,label='TT')\n",
    "plt.plot(ellsP,crossClsMod[0,1]*DlFac,label='TE')\n",
    "plt.plot(ellsP,crossClsMod[1,1]*DlFac,label='EE')\n",
    "plt.legend()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "# now, are the units for Cl and Nl compatible?\n",
    "# sensitivities given in units of microK-arcmin should be converted to microK-radian for microK^2 Cl\n",
    "\n",
    "# convert unitless spectra to microK^2\n",
    "TCMB = 2.7260 #Fixsen\n",
    "crossClsMod *=(TCMB*1e6)**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD8CAYAAAB+UHOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4xLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvAOZPmwAAIABJREFUeJzt3XmcXNV14PHfqb33vSWhbu3NJvbI\nLMbGGLEZLzg2xMTYlgkJExsSZ2aSGJJJIHYc43g+tkNm4oTY2DgDBowXFALGMptXFoGEWITQrm61\nWt3qfauu7cwf71VTrd6qu6u3V+f7+dSnqu6779W9pdY7de99715RVYwxxuQf33wXwBhjzPywAGCM\nMXnKAoAxxuQpCwDGGJOnLAAYY0yesgBgjDF5ygKAMcbkKQsAxhiTpywAGGNMngrMdwEmUl1dratW\nrZrvYhhjzKLy0ksvHVPVmsnyLegAsGrVKrZu3TrfxTDGmEVFRA5mk8+6gIwxJk9ZADDGmDyVVQAQ\nkXIReVhE3hSRnSJygYhUisgWEdntPle4eUVE7hKRPSKyQ0TOyTjOJjf/bhHZNFuVMsYYM7lsxwD+\nCfipql4jIiGgEPgr4ElVvVNEbgVuBT4PvA9ocB/nAd8EzhORSuB2YAOgwEsisllVO3NaI2OMmYJ4\nPE5TUxPRaHS+izJlkUiEuro6gsHgtPafNACISClwEfBpAFWNATERuRq42M12L/AMTgC4GvieOgsN\nPOe2Hpa5ebeoaod73C3AlcD3p1VyY4zJgaamJkpKSli1ahUiMt/FyZqq0t7eTlNTE6tXr57WMbLp\nAloDtAHfEZFtIvItESkClqjqEbcgR4BaN/9yoDFj/yY3bbx0Y4yZN9FolKqqqkV18gcQEaqqqmbU\ncskmAASAc4BvqurZQD9Od8+45RojTSdIH7mzyE0islVEtra1tWVRPGOMmZnFdvJPm2m5swkATUCT\nqj7vvn8YJyAcdbt2cJ9bM/LXZ+xfBzRPkD6Cqt6tqhtUdUNNzaT3MeTU468eYeeRnjn9TGOMmS+T\njgGoaouINIrISaq6C9gIvOE+NgF3us+PuLtsBm4RkQdwBoG7VfWIiDwB/EP6aiHgcuC23FZn+p7b\n185n7nuZUMDHq3dcTjjgn+8iGWPyQHt7Oxs3bgSgpaUFv99PZWUlfr+fWCxGR0cHfr+f9A/iF154\ngVAolJPPzvYqoD8B7nOvANoH3IDTenhIRG4EDgHXunkfA64C9gADbl5UtUNEvgi86Ob7QnpAeCH4\n2etHAYglUjy/r4OLTpzb1ocxJj9VVVWxfft2AO644w6Ki4v58z//8+HtY6XlSlYBQFW341y+ebyN\nY+RV4OZxjnMPcM9UCjhXXjrYwSnLSnmzpYftjV0WAIwxnreg5wKaK6mU8mZLL588fyWDsQS7Wnrn\nu0jGmHnwd//5Om8053Yc8NQTSrn9g+tzesxcsakggNbeIYYSKVZWF7GquogD7f3zXSRjjJl11gIA\nDron/BWVhayqKuLF/R2o6qK9NMwYMz0L9Zf6bLEWAHCwYwCAlZWFrKwqpD+WpL0/Ns+lMsaY2WUB\nAGjqGMAnsLyigGVlBQAc7Vl884IYY8xUWBcQzhhAVXGYoN/HktKwk9YzxPoT5rlgxpi8cscdd2SV\nlivWAgDaeoeoKXZO/EtKI4C1AIwx3mcBAGjrG6KmxAkA6eejPUPzWSRjjJl1FgBwWwDuiT/o91Fd\nHOJor7UAjDHelvcBQFU51jdEtdsFBFBTEqHVuoCMMR6X9wGgezBOPKnDLQCAqqIQHXYZqDHG4/I+\nALT1On39mQGgoihE50B8vopkjDFzIu8vAx0OABldQBWFQToHrAVgjJl9Y00HnZ76+ZVXXuHMM88c\nznvddddx660Trcc1NXkfANJ3/FYVvz2/dkVhiO7BOIlkioA/7xtJxphZNNF00MXFxcPbZkPen926\nBp2unvKC4HBaZVEIVWd8wBhjvCrvWwA97km+NCMAVBQ5rYHOgRhVGV1DxhiPe/xWaHk1t8dcejq8\n785p7To4OMhZZ501/P62227jYx/7WK5KZgGgayBGQdBPJPj2EpAVhU4w6Oi3FoAxZv4UFBTMaheQ\nBYCBOOWFwRFpFYVvtwCMMXlkmr/UF6u8HwPoHoxTVjAyAFSmu4DsXgBjjIdZC2CCAGBrAhhj5tPx\nYwBXXnkld96Zu1ZK3geA7oE4q6oLR6RFgn5CAR89URsDMMbMneOnfk4mk7P6eXnfBdQ1GBvVAgAo\njQTpGUzMQ4mMMWZu5H0A6B6MU14YGpVeWhCg11oAxhgPy+sAEI0nicZTY7YASiJBeqLWAjDGeFdW\nAUBEDojIqyKyXUS2ummVIrJFRHa7zxVuuojIXSKyR0R2iMg5GcfZ5ObfLSKbZqdK2Uvf6Xv8ZaAA\npZHA8E1ixhjjRVNpAbxXVc9S1Q3u+1uBJ1W1AXjSfQ/wPqDBfdwEfBOcgAHcDpwHnAvcng4a8yUd\nAMYcAygI2iCwMcbTZtIFdDVwr/v6XuDDGenfU8dzQLmILAOuALaoaoeqdgJbgCtn8Pkz1jWQngdo\njDEAGwQ2xnhctpeBKvAzEVHg31T1bmCJqh4BUNUjIlLr5l0ONGbs2+SmjZc+b7rcO33H7AIqCFgL\nwBgzJ/x+P6effvrw+/S0zxdffDFHjhyhoKAAgHXr1vHwww/n7HOzDQAXqmqze5LfIiJvTpBXxkjT\nCdJH7ixyE07XEStWrMiyeNOTHuQtiYz+GkojQWKJFNF4csQ8QcYYk2sTzflz3333sWHDhjG3zVRW\nXUCq2uw+twI/xunDP+p27eA+t7rZm4D6jN3rgOYJ0o//rLtVdYOqbkgvijBb+txf+CWRsccAAHrt\nSiBjjEdN2gIQkSLAp6q97uvLgS8Am4FNwJ3u8yPuLpuBW0TkAZwB3263i+gJ4B8yBn4vB27LaW2m\nqG/IObkXhUf/wi91WwU90fiI5SKNMd71lRe+wpsdE3VwTN3JlSfz+XM/P2GeiaZ9vv7664e7gC67\n7DK++tWv5qxs2XQBLQF+LCLp/Per6k9F5EXgIRG5ETgEXOvmfwy4CtgDDAA3AKhqh4h8EXjRzfcF\nVe3IWU2moXcoQSjgIxwYKwA4LQC7FNQYM9vmqwto0gCgqvuAM8dIbwc2jpGuwM3jHOse4J6pF3N2\n9EUTlITH/gpKC9ItAOsCMiZfTPZL3Wvy+k7gvqHEmAPAYC0AY4z35fVsoH3RBMXjBQAbBDbGzJGJ\npn3OHAOorq7m5z//ec4+N68DQO9QguLxuoDSLQC7F8AYM8vGm/b5mWeemdXPze8uoGiC4vDoS0AB\nIkEfAZ9YF5AxxrPyOwBMMAYgIjYfkDHG0/I+AIzXBQRQHA7QZ2MAxniec/Hi4jPTcud3AJhgEBjc\nADA0u0uyGWPmVyQSob29fdEFAVWlvb2dSCQy7WPk7SDwUCJJLJmauAUQCdA3ZF1AxnhZXV0dTU1N\ntLW1zXdRpiwSiVBXVzft/fM2APRNMBFcWkk4wNHe6FwVyRgzD4LBIKtXr57vYsyLvO0CSs8DNFEL\noMjGAIwxHpa3ASB9g9fkXUAWAIwx3mQBYJIuILsT2BjjVXkbANK/7EvGuREMnNbBUCJFPJmaq2IZ\nY8ycyeMA4FzdM+FloO62fusGMsZ4UP4GgCzGAIrcbdYNZIzxorwNAL1D2V0GCthAsDHGk/I2APRF\nEwT9Qjgw/leQ7gKyAGCM8aL8DQDuPEDuUpdjKrYWgDHGw/I3AEwyDxBkBAAbAzDGeFDeBgBnMZjx\nLwEF6wIyxnhb3gaAiRaET7MWgDHGy/I3AAxN3gVUFHIvA7UWgDHGg/I7AEzSAvD5hOJwwG4EM8Z4\nUtYBQET8IrJNRB51368WkedFZLeIPCgiITc97L7f425flXGM29z0XSJyRa4rMxW9WQwCAxSF/dYF\nZIzxpKm0AD4H7Mx4/xXg66raAHQCN7rpNwKdqroO+LqbDxE5FbgOWA9cCfyLiPhnVvzp6xuKTzoG\nAOlVwSwAGGO8J6sAICJ1wPuBb7nvBbgEeNjNci/wYff11e573O0b3fxXAw+o6pCq7gf2AOfmohJT\nFU+miMYnXg0srTgStDEAY4wnZdsC+Abwl0B6WswqoEtV02fGJmC5+3o50Ajgbu928w+nj7HPnEr3\n6WfTBVQSDtAXtWUhjTHeM2kAEJEPAK2q+lJm8hhZdZJtE+2T+Xk3ichWEdk6W2t0ZrMYTJp1ARlj\nvCqbFsCFwIdE5ADwAE7XzzeAchFJn0HrgGb3dRNQD+BuLwM6MtPH2GeYqt6tqhtUdUNNTc2UK5SN\nviwmgksrCgfoH0rOSjmMMWY+TRoAVPU2Va1T1VU4g7hPqer1wNPANW62TcAj7uvN7nvc7U+pqrrp\n17lXCa0GGoAXclaTKXi7BTDxncDgBIle6wIyxnjQ5D+Bx/d54AER+XtgG/BtN/3bwH+IyB6cX/7X\nAajq6yLyEPAGkABuVtV5+WmdzWIwaekuIFWdcOI4Y4xZbKYUAFT1GeAZ9/U+xriKR1WjwLXj7P8l\n4EtTLWSuTWkMIBIgpTAYT1IYmkm8NMaYhSUv7wSeyhiAzQdkjPGq/AwA0WkEALsSyBjjMfkZAIYS\n+AQKgpPfiGwBwBjjVXkZAHqjk68Glja8JoB1ARljPCYvA0A2M4GmpfPZdBDGGK/JzwAQTVASmfwe\nALBBYGOMd+VnABhKUBTObiJSWxbSGONVeRkAeocSFE+1BWABwBjjMXkZAPqi2a0FABAO+Aj6xQKA\nMcZz8jIA9A8lsx4EFnGWhbQxAGOM1+RlAMhmQfhMRTYltDHGg/IuAKRSOqXLQMEZB+i1FoAxxmPy\nLgD0x7KfBiKtJBIYnkHUGGO8Iu8CQLorp2gKLYCSSNBaAMYYz8m/ADCFqaDTSiIBemxRGGOMx+Rd\nAOidwoLwaaXWAjDGeFDeBYDhqaCn2ALojTqrghljjFfkXwCYRgugJBIkmVIG47Y4vDHGO/IvAExz\nDACwbiBjjKfkXwBILwcZzm4uIMgMADYQbIzxjrwNANnOBgrOIDBAj7UAjDEekpcBIBL0EfBnX3Xr\nAjLGeFHeBQBnOcjsu38ASguC7r7WBWSM8Y68CwB9Q4kpTQMB1gIwxnjTpAFARCIi8oKIvCIir4vI\n37npq0XkeRHZLSIPikjITQ+77/e421dlHOs2N32XiFwxW5WaSF80PqUrgIDh5SN7Bq0FYIzxjmxa\nAEPAJap6JnAWcKWInA98Bfi6qjYAncCNbv4bgU5VXQd83c2HiJwKXAesB64E/kVEsh+JzZGpzgQK\nUBTy4xNrARhjvGXSAKCOPvdt0H0ocAnwsJt+L/Bh9/XV7nvc7RtFRNz0B1R1SFX3A3uAc3NSiyno\njU5tLQB4e1EYGwMwxnhJVmMAIuIXke1AK7AF2At0qWr6J3ETsNx9vRxoBHC3dwNVmelj7DNn+oYS\nU5oGIs1mBDXGeE1WAUBVk6p6FlCH86v9lLGyuc8yzrbx0kcQkZtEZKuIbG1ra8umeFPSN5SY0lTQ\nac6MoBYAjDHeMaWrgFS1C3gGOB8oF5H0mbQOaHZfNwH1AO72MqAjM32MfTI/425V3aCqG2pqaqZS\nvGzKT980uoAgPSOodQEZY7wjm6uAakSk3H1dAFwK7ASeBq5xs20CHnFfb3bf425/Sp1pNDcD17lX\nCa0GGoAXclWRbAwlUiRSOuVBYIDSAlsW0hjjLdmcCZcB97pX7PiAh1T1URF5A3hARP4e2AZ8283/\nbeA/RGQPzi//6wBU9XUReQh4A0gAN6vqnE6vOTwP0DRaACWRIL1DvbkukjHGzJtJz4SqugM4e4z0\nfYxxFY+qRoFrxznWl4AvTb2YuTGdmUDT0msCGGOMV+TVncDDawHMIADYojDGGK/IqwCQXtd3egHA\nWRRmIGaLwhhjvCGvAkC6Cyc9udtU2HxAxhivyasAkJ7Lp2waASC9JkA2l4KqKt/65T6+9rNdRG0Z\nSWPMAjX1vpBFLH0jV/pkPhXpoNGVxYRwD7zYyN//104ABmJJ/tcHTp3y5/VG43z58Tfxi/BXV51C\nQWjOp00yxnhcfgUA9+Q9nRvBKgpDAHT2xybMp6r8+y/2cVZ9OWuqi/h/zx/kc5c2DM8omq2/+883\nePilJgBCAR9/M40gYowxE8mvLqBonJJwAL9vrFkpJlZe6LYABiZuAexo6mbfsX4+ft4KPnnBSqLx\nFI+/1jKlzzrSPciPtx3mDy5czUfPqeP+5w/ZXcjGmJzLrwAwmJjWADBARZHbAhiYuAXwqz3HALjk\n5FrOqi9naWmEZ3dNbU6jH718mGRKueHCVXzi/BUMxpP87PWj0yq3McaMJ78CQDQ+rbuAwVkTIOAT\nOidpAfx2bzsnLy2hujiMiPDuhmp+ubuNZCr7+wee3HmU05eXUV9ZyJl15VQXh3j2rdxPjGeMyW/5\nFQAG49NuAYgI5YUhugfHbwEkU8pLBzs5b3XlcNo711XRE02wqyW7aSTa+4bY1tjFxlNqAfD5hHc3\n1PCrPcfsJjRjTE7lVwCIJqZ1BVBaRWGQzv7xWwCNHQMMxpOsP6FsOO2cFRUAbGvszOozXjzQiSq8\nu+HtmVDPX1NJR3+MA+0D0yy5McaMll8BYDBOacH0L3yqKAxNOAbwpvsr/8SlJcNpKyoLqSwKse1Q\nV1afsa2xk6BfWH9C6XDaGXXlALzSmN0xjDEmG/kVAKLxGbUAyguDE14F9NZRNwAsKR5OExHOri9n\n26HsWgDbD3Vx6gllRIJvX/ffUFtMQdDPK00WAIwxuZM3ASCVUvqGpn8VEEzeAth1tJcVlYUUhka2\nMs5eUc7etn66JxlATiRT7Gjq5uz68hHpAb+P05aXsqOpe9plN8aY4+VNAOgdSqAKpdO8CgigvMhp\nAYw3GLurpZeTMrp/0s50T+ivHp74BP7W0T4G40nOXlE+atvJS0t5q6XXBoKNMTmTNwEgfRfwTFsA\nsWRqzBlBhxJJ9h/r56QlowPAGcvdPvxJunDSA8Vn1Y8OACcuKaZ3KMHRnqHpFN0YY0bJnwDg3kk7\n06uAYOybwfa29pNM6ZgtgLLCIKuqCicdxN1+qIvKohArKgtHbWtwA0t6nMEYY2YqfwLAYHoq6Ol3\nAVUWhQE41jc6AKRPzGMFAHC6gSbrw9/W2MXZ9eWIjJ6qoqG2eMTnGGPMTOVPAMhBC2BJqRMAWnui\no7a92dJL0C+sri4ac98z6spp6YmOuS9A92CcPa19Y3b/AFQVh6kqCrGntW+apTfGmJHyJwDMYC2A\ntNqSCACtvaP74d862svammKC/rG/0jPrnJvDXhmnFbDDHR84271xbCwNS4qtBWCMyZm8CQDdbgCY\n7lxAANXFIUTGDgC7Wno5cYwB4LT1J5Th98nwif542w91IQJn1JeNuR2gobaE3a19diWQMSYn8ioA\niDDlefkzBfw+qopCtPWO7MbpjcY53DU4bv8/QEHIz4lLStg+zkDwtsYu1tUUT9hFtbq6iN5ogvZJ\n1iQwxphs5E0A6ByIUVYQnNZaAJlqSiK0Hncp5vAA8AQtAHC6gV493D3qF7yqsu1Q5/C8QeNJjy8c\nbO+farGNMWaUSQOAiNSLyNMislNEXheRz7nplSKyRUR2u88VbrqIyF0iskdEdojIORnH2uTm3y0i\nm2avWqN1DsSHV/WaidqS8KguoF0tzsDsRC0AcAaCuwbiHOoYOanbgfYBOgfiY94AlmmVGwD2H7NJ\n4YwxM5dNCyAB/E9VPQU4H7hZRE4FbgWeVNUG4En3PcD7gAb3cRPwTXACBnA7cB5wLnB7OmjMha6B\n2PCqXjPhBICRXUC7WnooCvlZXl4w4b5n1o89EJyeJ2iiAWCAuooC/D7hwDFrARhjZm7SAKCqR1T1\nZfd1L7ATWA5cDdzrZrsX+LD7+mrge+p4DigXkWXAFcAWVe1Q1U5gC3BlTmszgc7+3LQAlpUX0NY7\nRCyRGk7b6U4B4Zuke+nEJSWEAz52HDcO8NLBTorDAdbVFo+zpyPo91FfUcB+6wIyxuTAlMYARGQV\ncDbwPLBEVY+AEySAWjfbcqAxY7cmN2289DmRqxbAispCUgpNnU43jKqy80gPpywrnWRP5wS+/oTR\nk7r9Zm87566uzGp8YmVVkbUAjDE5kXUAEJFi4IfAn6lqz0RZx0jTCdKP/5ybRGSriGxta8vdMoid\nA3Eqc9ACWFXlTNNw0O3HP9w1SG80kVUAADirvoIdh7uIxp35hBo7Bth/rJ93ravOav/V1U4AsEtB\njTEzlVUAEJEgzsn/PlX9kZt81O3awX1uddObgPqM3euA5gnSR1DVu1V1g6puqKmpOX7ztETjSQbj\nyeGF3WdihRsADrmrc+084lwBlG0AuPikGqLxFL/a7Swe/8wu52t7d0N2AWBVVSH9sSRtfTYpnDFm\nZrK5CkiAbwM7VfVrGZs2A+kreTYBj2Skf8q9Guh8oNvtInoCuFxEKtzB38vdtFmXXsQlF11ANcVh\nCkN+DroB4LXD3YjAyZNcAZR2/poqSiMBHnvtCACPbG+mobZ40v7/tPSVQAfsSiBjzAxl0wK4EPgk\ncImIbHcfVwF3ApeJyG7gMvc9wGPAPmAP8O/AZwFUtQP4IvCi+/iCmzbr0rN35mIQWERYV1vMmy1O\nL9hLBzs5ZWkpReHs7jAOBXx84MwTePSVI2x+pZmtBzu55nfqxpwAbiyrhwOAjQMYY2Zm0rOWqv6K\nsfvvATaOkV+Bm8c51j3APVMpYC6kA0AuWgAApy8vY/P2ZmKJFC8f6uTa36mb0v6fec9afvzyYf70\n+9tYWhrhE+evzHrf5eUFBHzCAbsSyBgzQ3lxJ3C6CygXLQBwpnbuHUrwvd8eYCCW5IK12fXfp9VX\nFnL/H53Hf7toDff/0XlZtx7AmY6ivrLQAoAxZsamPzPaIpLLLiCA955Ui98n/P1/7aQo5Ofik6Y+\nWH32iopJb/waz6qqQrsb2BgzY3nVAshVF1BNSZiPnO3cwvDH71lLJOjPyXGztaq6iIPtdimoMWZm\n8qMF0B+jIOjP6Yn6Kx89g8++d93wfQFzaXV1EQOxJEd7hlhaFpnzzzfGeENetACcieBy8+s/zedz\nVv/K9uqdXFpb41wyuq/NVgczxkxfXgSAjv6hnNwEtlCkA8BeCwDGmBnIiwBwrC9GdXF4vouRM0tK\nwxSF/OxtsyuBjDHTlycBYMhTAUBEWFtbbC0AY8yMeD4AqCrtfTGqS7zTBQRON9DeVgsAxpjp83wA\n6BlMEEumqPFQCwBgTXURzd1R+ocS810UY8wi5fkAkJ4100tdQABr3cnj9tucQMaYafJ8ADjm1QBg\nVwIZY2YofwKAx8YAVlYV4hPsSiBjzLR5PwD0erMFEAn6qa8sZPfR3vkuijFmkfJ+AOiL4ZPcTQS3\nkJy6rJQ3jky0OqcxxowvDwLAEJVF4awWXF9sTltexsH2AboH4/NdFGPMIpQXAaC62Hu//gHWn+Cs\nQ/xG8+StgMaOAXY0ddkMosaYYZ6fDbStL0ZNibf6/9PWn1AGwOvN3VywtmrcfF/72S7uemoP4Cw+\n/2+f/B0KQ57/pzfGTML7LYBeb00DkammJMzS0givHu4eN88vd7dx11N7+N2zl/PXV53Cr/cc44uP\nvjGHpTTGLFSe/hmYSimtvVFqS70ZAADOWVnOi/s7UNVRU1PHEilu3/w6q6oKufOjpxMO+GnrG+Lu\nX+zj+vNWctrysnkqtTFmIfB0C6BjIEY8qSwr9e6iKResraa5O8qhjtFLRH7n1/vZ19bP7R9cTzjg\nLIZzyyXrKI0E+MbPd891UY0xC4ynA0BLdxSApWUF81yS2fNOt+//N3vbR6S3dEe568ndXHpKLe89\nuXY4vTQS5NPvXMXPdx7lULutK2xMPsuTAODdFsCa6iKWlxfws9dbRqR/+fGdxFPK33zg1FH7/P55\nK/AJfP/FQ3NVTGPMAuTpAHCkxwkAyzwcAESED511Ar/YfYw2967n3+5t55HtzfzxRWtYWVU0ap9l\nZQVsPGUJP9jaSDyZmusiG2MWiEkDgIjcIyKtIvJaRlqliGwRkd3uc4WbLiJyl4jsEZEdInJOxj6b\n3Py7RWTT7FTH0dgxwB2bX+dXu9vw+8SzVwGlffScOlKq/Nuze+kejHPrj3aworKQz1y8btx9rntH\nPcf6Yjz9ZuscltQYs5Bk0wL4LnDlcWm3Ak+qagPwpPse4H1Ag/u4CfgmOAEDuB04DzgXuD0dNGZD\n92Cc7/7mAM/t66C2xJt3AWdaV1vMde+o51u/2s+lX3uW5q5B/ve1Z1IQ8o+7z3tOrKGmJMxDW5vm\nsKTGmIVk0gCgqr8AOo5Lvhq41319L/DhjPTvqeM5oFxElgFXAFtUtUNVO4EtjA4qOZM+8XUPxj3d\n/5/p9g+u54YLV3Hy0hK+8+lzOXd15YT5A34fHzlnOU/vaqW1NzpHpTTGLCTTHQNYoqpHANzn9GUm\ny4HGjHxNbtp46bOiMOOXr5f7/zNFgn5u/+B6/uPG83hXQ3VW+1z7O/UkU8pPth2e5dIZYxaiXA8C\nj9XXohOkjz6AyE0islVEtra1tU2rEJnTHKyoHD0Iahzraos5Z0U5P9jaZHMEGZOHphsAjrpdO7jP\n6ZHEJqA+I18d0DxB+iiqereqblDVDTU1NdMqXGYLYGVV4bSOkS+u3VDP7tY+XmkafzoJY4w3TTcA\nbAbSV/JsAh7JSP+UezXQ+UC320X0BHC5iFS4g7+Xu2mzIuh/u1orKy0ATOQDZywjEvTx0NbGyTMD\nqkoyZa0FY7wgm8tAvw/8FjhJRJpE5EbgTuAyEdkNXOa+B3gM2AfsAf4d+CyAqnYAXwRedB9fcNNm\n3eoa6wKaSEkkyFWnLeM/tzczGEuOmy+RTPHNZ/Zy7j88ScNfP8ame16gcYzpJ4wxi8ekk8Gp6u+P\ns2njGHkVuHmc49wD3DOl0s1AcThA31CCpR6eByhXrt1Qz4+2Heanrx/hd8+uG7U9kUzxZw9u59Ed\nR7j4pBoaaot58MVGfvdffs2PP3sh9dbKMmZRkoU8+LdhwwbdunXrtPZt7BggkVJWV1sLYDKplHLZ\n158F4Ik/u4hARhdaPJnicw+QcuhOAAARfElEQVRs47FXW/irq07mpovWArD7aC8f/eZvqKso5JFb\nLhzR7WaMmV8i8pKqbpgsn2f/19ZXFtrJP0s+n/CXV57M3rZ+vvubA8PpQ4kkt9z/Mo+92sL/ev8p\nwyd/gIYlJfzjNWfwxpEevvPr/fNQamPMTHl6PQCTvctPXcLGk2v5yk/fpKo4xElLSvnCo6/z3L4O\n7vjgqXz6wtWj9rli/VIuPaWWr2/ZzdVnLWeJdbcZs6h4tgVgpkZE+NrvncX6E8r47w++wlV3/ZId\nTd187ffOHPPkn97nbz+wnngyxb88vWeOS2yMmSlrAZhhZYVBfvDHF/D0m610D8Z5z0k11JZM/Kt+\nRVUh126o4/svNHLTe9ayvNy7ay8Y4zXWAjAjBP0+Ll+/lGs31E968k+75ZIGAP75SVtlzJjFxAKA\nmbHl5QV8/LwV/OClJg4c65/v4hhjsmQBwOTEZ9+7lqBf+CdrBRizaFgAMDlRWxJh0ztX8ZPth3nr\naO98F8cYkwULACZn/viitRSFAnx9y1vzXRRjTBYsAJicqSgKceO7VvP4ay28eGBOpnoyxsyABQCT\nUzddtIa6igI+//AOovHxJ5czxsw/CwAmp4rCAb78kdPZd6yfz/9wx7gLzagqTZ0D7GntnXAWUmPM\n7LEbwUzOvbuhhr+44iS++sQuEinl9g+eSm1JhJ5onOf3dfD0rlaeebOV5m5nLWKfwEUn1nDze9fx\njlUTr2VsjMkdz84GauaXqvKvz+7jq0+8iQJVRSE6+mOkFIpCft7VUM27GmooCQd4s6WXH73cRGvv\nEB/bUM8dH1pPQcaqbsaYqcl2NlALAGZW7T/Wz+btzbT0RKktCXPemko2rKwkFBjZ+zgQS3DXk3v4\nt1/s5cTaEv7v9eewrrZ4nkptzOJmAcAsSs++1cZ/f3A7Q/EkX/7oGXzozBPmu0jGLDp5vx6AWZze\nc2IN//Wn7+LkZaX86fe38Tc/ec0GiY2ZJTYIbBacZWUFPHDT+fzjT9/k33+5ny1vHOUP372aK9Yv\npa7CmW20ezBOU+cgLd1R2vqG6IsmKI4EWFNdxOl1ZRSG7E/bmMlYF5BZ0F7Y38FXfvomLx3sBCDg\nE0Qgnhz/7zYccGY0/f131HPB2ipEZK6Ka8yCYGMAxlP2tPby3L4ODncNogrVxSHqKgpZWhahtiRM\ncSRAz2CcXS29PLOrjf/c0UzXQJwz68u55b3r2HhyLT6fBQKTHywAmLwWjSf54ctN/Ouze2nsGGRt\nTRE3XLiaj55TZ5eYGs+zAJALzdvh9R/BYBcsPwfO+BgEbcWrxSSRTPHojiN8+1f7efVwN0UhP+9u\nqOG8NZWsrSmmojBEQchH0O8j4PcR9AtBn4+AXwj6nXS/tRzMImMBYCaScXjsL+Cl74AvCJFSGGiH\nshXwe/c6wcAsKqrKiwc6+cn2wzy1s5WWnmjW+y4tjdCwpJgTl5RwovvcsKSE4rANNJuFacEGABG5\nEvgnwA98S1XvHC/vvASAaDc89CnY9wxccAu85y8hXAoHfgk/+Sz0H4OPPwhr3jO35QKnJbL1Htj/\nLKjCygthwx+QKqoilowhIoT94bkv1yKjqhzri7GvrY+eaIKBWIJEUkmkUsSTSiKZIpFS4kllKJHk\nYPsAbx3tZU9rH0OJ1PBxlpcXsLyigKqiEBVFIYI+wecTfCIEfDKiFRHy+ygrDFJTHKa6OEx1SYiq\novCoG+KMyYUFGQBExA+8BVwGNAEvAr+vqm+MlX/OA0B3E9x3LRx7Cz54F5x9/fCmlKbo7TxAzwPX\n0dXXQvdltxNYup7iYDFFwSLKwmWUhkoJ+I77VdjXxsC2ezl49BUOBP0cqlzJIUnQ2NtIb6yXeCpO\nQaCAinAFNYU1LClcwtKipSMeJcEShg78imOP/Df2D7Wzr2oF+33KwVg3B0Mh2vxvn0QKAgUsKVzC\n6rLVrC1fy5qyNdSV1LGsaBk1BTX4fX6nLrFeOqOdtEfbOTZ4bMQjlowR9AWJBCLUFNRQW1jLksIl\nVBVUUVVQRUW4Ar9vZD/6UHKInqEeemI9DMQHGEwMEk1GiSVjhPwhIv4IkUCEomARRcEiioPFFAYL\n8cniOQEmU0pjhxMMnEcfLT1ROvpjdA3EiCeVlCqplJJwH8nUxP+/ygqCVBeHqC4OU1MSpqooRCTo\nJxRwA4f7HPQLgeGuKed10A0w6a6rcMBHyO/sO/zw+wgHnWcAnwg+cZ5FsCukPGqhBoALgDtU9Qr3\n/W0AqvrlsfJPNwA0N2/le0/9BYlAiBRCSnwkwiWkiqtJhIpIaYqUpkikEqQ0RVKTJAc7SbW8SkJT\npGpOJBkuJplK0h/vp3uom+5YNylNTfrZRcEiykJllIRKiA0co2egjXb/yJNcrYSor15PRUEVQV+Q\nwcQg7dF22gbbaB1oHfU5AfGT0JE3Q1VGKllZUMuK9kMs6z5CuO4dpNZdSle8l5b+FvZ17eNgz0ES\nmhixn098oEqK0f/uAXxUFlQRCRQQS8UYTAzSPdQ9Kp9PfBQHi/EhkIwRTUaJjnG8bBQGCt0gWkgB\nQiLWTyw5RDyVIKZJ4ihJICWQUiWJ4vcFiPgjhANhwn7nEfFHCPtD+FIpSMWRVMJ9JNFkDFIJNJVA\nxQ8+Pyo+1OcH8aEIiA+fz4/gNE1FFZ8qPtxn97VoCp8qmadNAcRNEfdBRpqioEIylRoOCk6Lww0S\nbusjvU0VdxbViU/OY21Nf+KI9zpWnrHLn7lxuC7pTW6RJCOTHJff2X788SfOP2aZjkvP/H4zX2S+\nH/E54xxvqibaTab5N5+tlaWn8D8+9s/T2jfbADDXnZjLgcaM903AeZkZROQm4CaAFStWTOtDunqb\n2Jw4hj+u+BQCKL6+FP5j4PcH8YVKCIRL8AUL8KdS+Ac78fUdxe8P4a89jVCoBL/Pj0981JXUUR4u\npzRUSnm4nPJIOWWJBGVb7iAx1Ev/xr+hp6SGnpjz67dnqIfu3sP0Nj1PuK+NksIaljZ8kFUnvINV\nBbWs2PETCn79TzBYBJ/4OhRVjSh7IpXg2OAxWvpbnEfzi3Rtu5fiYDEVF/wpq5ecxZqyNZRHyp0d\nkgn4+e3w2/8D3X1w7XehdBkA8VScxt5Gmvuaae5rpr39LRL7n4G2tyiXIOU1J1NdWk91PEb14e2U\nHduLL9IJZ38S3nEjVK5hKDlE60ArrQOttA+2Oy2G1tfoObwVOvahqTgFwSJKS+spLaym1BeiaLCb\nSPs+CrqaCKLE/CGiNScRrV5Hf6SEfoH+WB99g8foGzjGQH8rffEBBnxCUJWQKmFVgkBQFb86/9n8\ngE8hKTDkDxENFhATiGqSIU0SJUkCUPe/ZuZ/T+dk40M0hfD2CVwyMiUE5weDu28q431KnLRk+tju\nAdT9vLdfvy0zz3DacX+rCk7E8Y/eNjKfjHinY5yadIyz1ZifN+Gxx/v8uTcbn7lwRzxHi3UPzfpn\nzHUL4FrgClX9Q/f9J4FzVfVPxsqfsy4gVeg6CG/9DF77ITQ+z6jTw5nXwRX/AIVZTkfc1QjfvQr6\nWmHj38LZn4DEELx8L/zya84xL70D3vGH4Duum+OtJ5xxhopV8KlHoGTp2J/RvA3uvdop0w2PQekE\n8+K8+jBs/hMIFcNVX4VTr377Z0/nQfjFV2H7/c5VTBfcAhfc7AxuZ35HjS/A89+ENzaDpmDdpdBw\nGVSshkQUjr4Gux6Hlh0QLITTPgJnXQ8rLhj7J1b/MTj0HBz6rfM48gqkMlokwUKoWgtVDVBzMtSc\n5DzKV759tVWsDwY7nWP1tzmPvta3X6tCqBCCRc73VFDhPsozXldAuMz5d0ilID7gHDfW79QrEXX+\n7RJRZ9A/EHYe/jAEQu5zGPyht5+t68QsYHndBTSpvjZnILW7yTlRrHkvVKycxnFa4SefgT0/H5l+\n0vvhfXdC+QQtmP2/gPuvc1oAn/gRVDeM3H74ZfiPDzsnrhseg/L6ycvTuhMe/gNofcP57CWnQW8L\nNL/snLQ23Ajv/p9QXDPxcXqa4aXvwvbvQ/ehjA0CdRvgtGucgFlQPnmZMqVSEO2C+KBz0k6flI0x\nObVQA0AAZxB4I3AYZxD446r6+lj55/0+gGykfzkf/LVzkl1zMSw9Lbt9D78M9/8epJLOr/b1HwEU\ntt8Hj38eimrg049OHEiOl0w4rZw3HnFaPQUVsPoiOOvjUFY39fp1NzlBxB90WiyRsqkfwxgzpxZk\nAAAQkauAb+D0fN6jql8aL++iCAAz1bEPfnADHNkOBZVO10u0C1a8E665Z7g/3xhjsrVQB4FR1ceA\nx+b6cxesyjXwR085v9j3PgXic/rdT3q/dY8YY2aV3cq4EPj8zoDqaR+Z75IYY/KI/cQ0xpg8ZQHA\nGGPylAUAY4zJUxYAjDEmT1kAMMaYPGUBwBhj8pQFAGOMyVMWAIwxJk8t6CUhRaQNODjN3auBYzks\nzmJgdc4PVuf8MJM6r1TVSWZ9XOABYCZEZGs2c2F4idU5P1id88Nc1Nm6gIwxJk9ZADDGmDzl5QBw\n93wXYB5YnfOD1Tk/zHqdPTsGYIwxZmJebgEYY4yZgCcDgIhcKSK7RGSPiNw63+WZCRG5R0RaReS1\njLRKEdkiIrvd5wo3XUTkLrfeO0TknIx9Nrn5d4vIpvmoSzZEpF5EnhaRnSLyuoh8zk33cp0jIvKC\niLzi1vnv3PTVIvK8W/4HRSTkpofd93vc7asyjnWbm75LRK6YnxplT0T8IrJNRB5133u6ziJyQERe\nFZHtIrLVTZu/v21V9dQDZ6nJvcAaIAS8Apw63+WaQX0uAs4BXstI+0fgVvf1rcBX3NdXAY8DApwP\nPO+mVwL73OcK93XFfNdtnPouA85xX5fgrCF9qsfrLECx+zoIPO/W5SHgOjf9X4HPuK8/C/yr+/o6\n4EH39anu33sYWO3+P/DPd/0mqfv/AO4HHnXfe7rOwAGg+ri0efvb9mIL4Fxgj6ruU9UY8ABw9TyX\nadpU9RdAx3HJVwP3uq/vBT6ckf49dTwHlIvIMuAKYIuqdqhqJ7AFuHL2Sz91qnpEVV92X/cCO4Hl\neLvOqqp97tug+1DgEuBhN/34Oqe/i4eBjSIibvoDqjqkqvuBPTj/HxYkEakD3g98y30veLzO45i3\nv20vBoDlQGPG+yY3zUuWqOoRcE6YQK2bPl7dF+V34jbzz8b5RezpOrtdIduBVpz/0HuBLlVNuFky\nyz9cN3d7N1DFIqsz8A3gL4GU+74K79dZgZ+JyEsicpObNm9/215cE1jGSMuXS53Gq/ui+05EpBj4\nIfBnqtrj/NgbO+sYaYuuzqqaBM4SkXLgx8ApY2Vznxd9nUXkA0Crqr4kIhenk8fI6pk6uy5U1WYR\nqQW2iMibE+Sd9Tp7sQXQBNRnvK8DmuepLLPlqNsUxH1uddPHq/ui+k5EJIhz8r9PVX/kJnu6zmmq\n2gU8g9PnWy4i6R9pmeUfrpu7vQynm3Ax1flC4EMicgCnm/YSnBaBl+uMqja7z604gf5c5vFv24sB\n4EWgwb2aIIQzYLR5nsuUa5uB9Mj/JuCRjPRPuVcPnA90u03KJ4DLRaTCvcLgcjdtwXH7db8N7FTV\nr2Vs8nKda9xf/ohIAXApztjH08A1brbj65z+Lq4BnlJndHAzcJ17xcxqoAF4YW5qMTWqepuq1qnq\nKpz/o0+p6vV4uM4iUiQiJenXOH+TrzGff9vzPSo+Gw+c0fO3cPpR/3q+yzPDunwfOALEcSL/jTh9\nn08Cu93nSjevAP/XrferwIaM4/wBzgDZHuCG+a7XBPV9F05zdgew3X1c5fE6nwFsc+v8GvC3bvoa\nnJPZHuAHQNhNj7jv97jb12Qc66/d72IX8L75rluW9b+Yt68C8myd3bq94j5eT5+b5vNv2+4ENsaY\nPOXFLiBjjDFZsABgjDF5ygKAMcbkKQsAxhiTpywAGGNMnrIAYIwxecoCgDHG5CkLAMYYk6f+P1iL\nkZnf+9WFAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x109098f10>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(ellsP,crossClsMod[0,0]*DlFac,label='TT')\n",
    "plt.plot(ellsP,crossClsMod[0,1]*DlFac,label='TE')\n",
    "plt.plot(ellsP,crossClsMod[1,1]*DlFac,label='EE')\n",
    "plt.legend()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "getting (primary CMB) detector noise...\n"
     ]
    }
   ],
   "source": [
    "# grab some FisherCl\n",
    "import noiseCl as ncl\n",
    "\n",
    "#arcminPerRad = 60*180/np.pi\n",
    "\n",
    "print 'getting (primary CMB) detector noise...'\n",
    "fwhm = 1; ST = 1; SP = ST*1.414\n",
    "#fwhm = 2; ST = 1; SP = ST*1.414\n",
    "#ST /= arcminPerRad; SP /= arcminPerRad\n",
    "noiseCMBS4_TT1 = ncl.noisePower(ST,ST,fwhm,ellsP)\n",
    "noiseCMBS4_PP1 = ncl.noisePower(SP,SP,fwhm,ellsP)\n",
    "\n",
    "# shape like crossCls\n",
    "noiseClsP = np.zeros(crossClsP.shape)\n",
    "noiseClsP[0,0] = noiseCMBS4_TT1\n",
    "noiseClsP[1,1] = noiseCMBS4_PP1\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "# also scale crossClsPPlus, crossClsPMinus, dClVecsP\n",
    "crossClsPPlus = Fobj.crossClsPPlus*(TCMB*1e6)**2/DlFac\n",
    "crossClsPMinus = Fobj.crossClsPMinus*(TCMB*1e6)**2/DlFac\n",
    "dClVecsP = Fobj.dClVecsP*(TCMB*1e6)**2/DlFac"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "starting covariance set  1  of  2 ... \n",
      "starting covariance set  2  of  2 ... \n"
     ]
    }
   ],
   "source": [
    "# insert these into Fobj and make covar\n",
    "Fobj.crossClsP = crossClsMod\n",
    "Fobj.noiseClsP = noiseClsP\n",
    "Fobj.crossClsPPlus = crossClsPPlus\n",
    "Fobj.crossClsPMinus = crossClsPMinus\n",
    "Fobj.dClVecsP = dClVecsP\n",
    "\n",
    "Fobj.covarP,Fobj.invCovP,Fobj.ellsP,Fobj.obsListP = \\\n",
    "            Fobj.makeCovar(Fobj.crossClsP,Fobj.noiseClsP,Fobj.lminP,Fobj.lmaxP)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n"
     ]
    }
   ],
   "source": [
    "# make Fisher matrices\n",
    "Fobj.FijTE = Fobj.makeFisher(Fobj.lminP,TE=True)\n",
    "#FijTE = Fobj.makeFisher(30,TE=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "creating noise arrays... \n",
      "getting (EB) lensing reconstruction noise... \n",
      "[00:00:00] calculating full-sky noise level for estimator of type <class 'quicklens.qest.lens.phi_EB'>\n",
      "getting galaxy shot noise... \n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/ipykernel_launcher.py:52: ComplexWarning: Casting complex values to real discards the imaginary part\n"
     ]
    }
   ],
   "source": [
    "# make a similar adjustment for Nl^kk\n",
    "# grab some FisherCl\n",
    "import crosspower as cp\n",
    "\n",
    "print 'creating noise arrays... '\n",
    "\n",
    "print 'getting (EB) lensing reconstruction noise... '\n",
    "# noise levels from a possible CMB-S4 design:\n",
    "nlev_t     = 1.   # temperature noise level, in uK.arcmin.\n",
    "nlev_p     = 1.414   # polarization noise level, in uK.arcmin.\n",
    "beam_fwhm  = 1.   # Gaussian beam full-width-at-half-maximum.\n",
    "\n",
    "#beam_fwhm  = 2.   # Gaussian beam full-width-at-half-maximum.\n",
    "#beam_fwhm  = 4.   # Gaussian beam full-width-at-half-maximum. \n",
    "#From Hu & Okamoto's \"near-perfect\" experiment\n",
    "\n",
    "ells,EB_noise = ncl.getRecNoise(Fobj.lmax,nlev_t,nlev_p,beam_fwhm)\n",
    "\n",
    "# convert Nl^dd to Nl^kk\n",
    "# this formula is a reasonable guess based on Fourier analog.\n",
    "Nlkk = EB_noise * ells*(ells+1)/4\n",
    "\n",
    "print 'getting galaxy shot noise... '\n",
    "# From Schaan et. al.: LSST n_source = 26/arcmin^2 for full survey\n",
    "#nbar = 26 # arcmin^-2\n",
    "nbar = 66 # 66 arcmin^-2 to match Bye's value\n",
    "\n",
    "# the selection of beesBins must be consistent with that which is selected in cp.tophat\n",
    "#beesBins = True\n",
    "beesBins = False\n",
    "if beesBins:\n",
    "    binEdges = [0.0,0.5,1.0,2.0,3.0,4.0,7.0]\n",
    "    nBins = 6\n",
    "else:\n",
    "    binEdges = np.linspace(Fobj.zmin,Fobj.zmax,Fobj.nBins+1)\n",
    "    nBins = Fobj.nBins\n",
    "\n",
    "# the selection of dndz function must be consistent with that which is selected in cp.getDNDZinterp\n",
    "# myDNDZ must be a function only of z\n",
    "#myDNDZ = lambda z: cp.modelDNDZ(z,z0)\n",
    "myDNDZ = lambda z: cp.modelDNDZ3(z,Fobj.z0)\n",
    "\n",
    "N_gg = ncl.shotNoise(nbar,binEdges,myDNDZ=myDNDZ)\n",
    "\n",
    "\n",
    "\n",
    "# create noiseCls array to accompany crossCls: \n",
    "#  Nl^kk will be reconstruction noise\n",
    "#  Nl^gigi will be shot noise\n",
    "#  Nl^kg, Nl^gigj are zero because noise is uncorrelated\n",
    "noiseCls = np.zeros(Fobj.crossCls.shape)\n",
    "noiseCls[0,0] = Nlkk[Fobj.lmin:Fobj.lmax+1]\n",
    "for binNum in range(nBins):\n",
    "    noiseCls[binNum+1,binNum+1] = N_gg[binNum]*np.ones(Fobj.lmax-Fobj.lmin+1)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "starting covariance set  1  of  2 ... \n",
      "starting covariance set  2  of  2 ... \n"
     ]
    }
   ],
   "source": [
    "# put into Fobj and make covar\n",
    "Fobj.noiseCls = noiseCls\n",
    "\n",
    "Fobj.covar,Fobj.invCov,Fobj.ells,Fobj.obsList = \\\n",
    "        Fobj.makeCovar(Fobj.crossCls,Fobj.noiseCls,Fobj.lmin,Fobj.lmax)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Back to our regularly scheduled programming"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "adjusting lmax for TT... \n",
      "(6, 6) (6, 6) (6, 6)\n"
     ]
    }
   ],
   "source": [
    "# get TE fisher\n",
    "myLmin = 30\n",
    "myFij = Fobj.makeFisher(myLmin,TE=True)\n",
    "\n",
    "# fix w, mnu \n",
    "# nCosParams should be 6 to fix w, mnu\n",
    "Fij_S4_TE = myFij[:nCosParams,:nCosParams] \n",
    "Fij_plikHM_lowTEB = Fij_plikHM_lowTEB[:nCosParams,:nCosParams]\n",
    "Fij_plikHM_TTTEEE_lowTEB = Fij_plikHM_TTTEEE_lowTEB[:nCosParams,:nCosParams]\n",
    "print Fij_S4_TE.shape, Fij_plikHM_lowTEB.shape, Fij_plikHM_TTTEEE_lowTEB.shape\n",
    "\n",
    "# convert As to 10^9 As\n",
    "Fij_S4_TE[:,AsIndex] *= 1e-9\n",
    "Fij_S4_TE[AsIndex,:] *= 1e-9\n",
    "\n",
    "# combine using appropriate fsky weighting\n",
    "fsky_S4 = 0.4\n",
    "fsky_Pl = 0.2 # this is in addition to area covered by S4\n",
    "\n",
    "# fsky is in denominator of covariance matrix, putting it in numerator of Fij matrix\n",
    "Fij_Sum = (Fij_S4_TE+Fij_plikHM_lowTEB)*fsky_S4 + Fij_plikHM_TTTEEE_lowTEB*fsky_Pl\n",
    "Fij_Sum_Pl = Fij_plikHM_TTTEEE_lowTEB*(fsky_S4+fsky_Pl)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "# add sigma = 0.01 prior on tau\n",
    "tauIndex = 5\n",
    "tauPrior = 1/0.01**2\n",
    "Fij_Sum[tauIndex,tauIndex] += tauPrior\n",
    "Fij_Sum_Pl[tauIndex,tauIndex] += tauPrior"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get sigmas\n",
    "sigmas_Sum = getSigmas(Fij_Sum)\n",
    "sigmas_Sum_Pl = getSigmas(Fij_Sum_Pl)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "# one more conversion before printout and comparison\n",
    "ombh2Fac = 100\n",
    "sigmas_Sum[0] *= ombh2Fac\n",
    "sigmas_Sum_Pl[0] *= ombh2Fac\n",
    "paramList[0] = '100 ombh2'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Planck only and Planck+CMB-S4\n",
      "+-----------+---------+-----------+\n",
      "| Parameter |  Planck | S4+Planck |\n",
      "+-----------+---------+-----------+\n",
      "| 100 ombh2 |  0.0174 |   0.0025  |\n",
      "|   omch2   | 0.00171 |  0.000712 |\n",
      "|     H0    |  0.759  |   0.286   |\n",
      "|  10^9 As  |  0.0393 |   0.0407  |\n",
      "|     ns    |  0.0048 |  0.00196  |\n",
      "|    tau    |  0.009  |  0.00922  |\n",
      "+-----------+---------+-----------+\n"
     ]
    }
   ],
   "source": [
    "# show sigmas\n",
    "from prettytable import PrettyTable\n",
    "\n",
    "print 'Planck only and Planck+CMB-S4'\n",
    "t = PrettyTable(['Parameter', 'Planck','S4+Planck'])\n",
    "for i,p in enumerate(paramList[:nCosParams]):\n",
    "    Zi = sigmas_Sum_Pl[i]\n",
    "    Bi = sigmas_Sum[i]\n",
    "    #diff = (Bi-Zi)/Zi *100\n",
    "    t.add_row([p,\"%5.3g\" % Zi,\"%5.3g\" % Bi] )\n",
    "print t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+-----------+-------------+----------+--------+\n",
      "| Parameter | Planck_book | Planck_Z | % diff |\n",
      "+-----------+-------------+----------+--------+\n",
      "| 100 ombh2 |    0.017    |  0.0174  |  2.61  |\n",
      "|   omch2   |    0.0014   | 0.00171  |  21.8  |\n",
      "|     H0    |      0.7    |  0.759   |  8.38  |\n",
      "|  10^9 As  |    0.039    |  0.0393  |  0.69  |\n",
      "|     ns    |    0.004    |  0.0048  |  20.1  |\n",
      "|    tau    |     0.01    |  0.009   | -9.98  |\n",
      "+-----------+-------------+----------+--------+\n",
      "+-----------+----------------+-------------+--------+\n",
      "| Parameter | S4+Planck_book | S4+Planck_Z | % diff |\n",
      "+-----------+----------------+-------------+--------+\n",
      "| 100 ombh2 |     0.003      |    0.0025   | -16.6  |\n",
      "|   omch2   |     0.0006     |   0.000712  |  18.7  |\n",
      "|     H0    |      0.24      |    0.286    |  19.1  |\n",
      "|  10^9 As  |     0.021      |    0.0407   |  93.6  |\n",
      "|     ns    |     0.002      |   0.00196   | -1.87  |\n",
      "|    tau    |     0.006      |   0.00922   |  53.7  |\n",
      "+-----------+----------------+-------------+--------+\n"
     ]
    }
   ],
   "source": [
    "# print out with CMB-S4 table 8.1 values\n",
    "sigmas_S4book_Planck = [0.017,0.0014,0.7,0.039,0.004,0.01]\n",
    "sigmas_S4book_PlanckS4 = [0.003,0.0006,0.24,0.021,0.002,0.006]\n",
    "\n",
    "t = PrettyTable(['Parameter','Planck_book', 'Planck_Z','% diff'])\n",
    "for i,p in enumerate(paramList[:nCosParams]):\n",
    "    Bi = sigmas_Sum_Pl[i]\n",
    "    Zi = sigmas_S4book_Planck[i]\n",
    "    diff = (Bi-Zi)/Zi *100\n",
    "    t.add_row([p,\"%5.3g\" % Zi,\"%5.3g\" % Bi,\"%5.3g\" % diff] )\n",
    "print t\n",
    "\n",
    "t = PrettyTable(['Parameter','S4+Planck_book', 'S4+Planck_Z','% diff'])\n",
    "for i,p in enumerate(paramList[:nCosParams]):\n",
    "    Bi = sigmas_Sum[i]\n",
    "    Zi = sigmas_S4book_PlanckS4[i]\n",
    "    diff = (Bi-Zi)/Zi *100\n",
    "    t.add_row([p,\"%5.3g\" % Zi,\"%5.3g\" % Bi,\"%5.3g\" % diff] )\n",
    "print t"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# that was without kappa.  Now get kappa into the S4 part"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(3, 3, 2999)\n",
      "(3, 9, 2999)\n",
      "2 3000\n"
     ]
    }
   ],
   "source": [
    "print Fobj.covar.shape\n",
    "print Fobj.dClVecs.shape\n",
    "print Fobj.lmin,Fobj.lmax"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "parameters:  ['ombh2', 'omch2', 'H0', 'As', 'ns', 'tau', 'mnu', 'w']\n",
      "kappa-only, noise-included Fisher matrix:\n",
      "[[ 1.42206889e+09 -1.58420001e+09  8.98451527e+05 -4.55824080e+16\n",
      "  -1.16209723e+08  2.51574140e+05  2.62964188e+07  3.69206028e+07]\n",
      " [-1.58420001e+09  1.76520321e+09 -1.00131934e+06  5.08187823e+16\n",
      "   1.29261301e+08 -2.80174719e+05 -2.93128138e+07 -4.11520360e+07]\n",
      " [ 8.98451527e+05 -1.00131934e+06  5.68266873e+02 -2.88532688e+13\n",
      "  -7.30323351e+04  1.58439597e+02  1.66376565e+04  2.33534966e+04]\n",
      " [-4.55824080e+16  5.08187823e+16 -2.88532688e+13  1.46614355e+24\n",
      "   3.69042144e+15 -8.02590267e+12 -8.45136302e+14 -1.18592061e+15]\n",
      " [-1.16209723e+08  1.29261301e+08 -7.30323351e+04  3.69042144e+15\n",
      "   9.80918158e+06 -2.10753653e+04 -2.13467173e+06 -3.00251394e+06]\n",
      " [ 2.51574140e+05 -2.80174719e+05  1.58439597e+02 -8.02590267e+12\n",
      "  -2.10753653e+04  4.56985488e+01  4.63891840e+03  6.51963222e+03]\n",
      " [ 2.62964188e+07 -2.93128138e+07  1.66376565e+04 -8.45136302e+14\n",
      "  -2.13467173e+06  4.63891840e+03  4.87275091e+05  6.83848789e+05]\n",
      " [ 3.69206028e+07 -4.11520360e+07  2.33534966e+04 -1.18592061e+15\n",
      "  -3.00251394e+06  6.51963222e+03  6.83848789e+05  9.59871214e+05]]\n"
     ]
    }
   ],
   "source": [
    "# make kk fisher with noise\n",
    "nCosParams = Fobj.nCosParams\n",
    "dClkk = Fobj.dClVecs[0] # the dCl/dkk part\n",
    "cov_kk2 = Fobj.covar[0,0] # the kk part\n",
    "Fij_kk2 = np.zeros((nCosParams,nCosParams))\n",
    "minEllInd = myLmin-Fobj.lmin  # eg. 30-2\n",
    "for pNum1 in range(nCosParams):\n",
    "    for pNum2 in range(pNum1,nCosParams):\n",
    "        # the next line does the summation over ell\n",
    "        Fij_kk2[pNum1,pNum2] = np.sum(dClkk[pNum1,minEllInd:]*dClkk[pNum2,minEllInd:]/cov_kk2[minEllInd:])\n",
    "        Fij_kk2[pNum2,pNum1] = Fij_kk2[pNum1,pNum2]\n",
    "print 'parameters: ',Fobj.paramList[:nCosParams]\n",
    "print 'kappa-only, noise-included Fisher matrix:'\n",
    "print Fij_kk2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "# convert As to 10^9 As\n",
    "Fij_kk2[:,AsIndex] *= 1e-9\n",
    "Fij_kk2[AsIndex,:] *= 1e-9\n",
    "\n",
    "# fix w, mnu\n",
    "Fij_kk2 = Fij_kk2[:6,:6]\n",
    "\n",
    "# get sigmas\n",
    "sigmas_Sum2 = getSigmas(Fij_Sum+Fij_kk2*fsky_S4)\n",
    "\n",
    "# one more conversion before printout and comparison\n",
    "sigmas_Sum2[0] *= ombh2Fac"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Planck only and Planck+CMB-S4\n",
      "+-----------+---------+-----------+\n",
      "| Parameter |  Planck | S4+Planck |\n",
      "+-----------+---------+-----------+\n",
      "| 100 ombh2 |  0.0174 |  0.00229  |\n",
      "|   omch2   | 0.00171 |  0.00055  |\n",
      "|     H0    |  0.759  |   0.219   |\n",
      "|  10^9 As  |  0.0393 |   0.0205  |\n",
      "|     ns    |  0.0048 |  0.00172  |\n",
      "|    tau    |  0.009  |  0.00543  |\n",
      "+-----------+---------+-----------+\n"
     ]
    }
   ],
   "source": [
    "# show sigmas\n",
    "from prettytable import PrettyTable\n",
    "\n",
    "print 'Planck only and Planck+CMB-S4'\n",
    "t = PrettyTable(['Parameter', 'Planck','S4+Planck'])\n",
    "for i,p in enumerate(paramList[:nCosParams]):\n",
    "    Zi = sigmas_Sum_Pl[i]\n",
    "    Bi = sigmas_Sum2[i]\n",
    "    #diff = (Bi-Zi)/Zi *100\n",
    "    t.add_row([p,\"%5.3g\" % Zi,\"%5.3g\" % Bi] )\n",
    "print t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+-----------+-------------+----------+--------+\n",
      "| Parameter | Planck_book | Planck_Z | % diff |\n",
      "+-----------+-------------+----------+--------+\n",
      "| 100 ombh2 |    0.017    |  0.0174  |  2.61  |\n",
      "|   omch2   |    0.0014   | 0.00171  |  21.8  |\n",
      "|     H0    |      0.7    |  0.759   |  8.38  |\n",
      "|  10^9 As  |    0.039    |  0.0393  |  0.69  |\n",
      "|     ns    |    0.004    |  0.0048  |  20.1  |\n",
      "|    tau    |     0.01    |  0.009   | -9.98  |\n",
      "+-----------+-------------+----------+--------+\n",
      "+-----------+----------------+-------------+--------+\n",
      "| Parameter | S4+Planck_book | S4+Planck_Z | % diff |\n",
      "+-----------+----------------+-------------+--------+\n",
      "| 100 ombh2 |     0.003      |   0.00229   | -23.5  |\n",
      "|   omch2   |     0.0006     |   0.00055   | -8.35  |\n",
      "|     H0    |      0.24      |    0.219    | -8.75  |\n",
      "|  10^9 As  |     0.021      |    0.0205   | -2.42  |\n",
      "|     ns    |     0.002      |   0.00172   | -13.8  |\n",
      "|    tau    |     0.006      |   0.00543   | -9.47  |\n",
      "+-----------+----------------+-------------+--------+\n"
     ]
    }
   ],
   "source": [
    "# print out with CMB-S4 table 8.1 values\n",
    "sigmas_S4book_Planck = [0.017,0.0014,0.7,0.039,0.004,0.01]\n",
    "sigmas_S4book_PlanckS4 = [0.003,0.0006,0.24,0.021,0.002,0.006]\n",
    "\n",
    "t = PrettyTable(['Parameter','Planck_book', 'Planck_Z','% diff'])\n",
    "for i,p in enumerate(paramList[:nCosParams]):\n",
    "    Bi = sigmas_Sum_Pl[i]\n",
    "    Zi = sigmas_S4book_Planck[i]\n",
    "    diff = (Bi-Zi)/Zi *100\n",
    "    t.add_row([p,\"%5.3g\" % Zi,\"%5.3g\" % Bi,\"%5.3g\" % diff] )\n",
    "print t\n",
    "\n",
    "t = PrettyTable(['Parameter','S4+Planck_book', 'S4+Planck_Z','% diff'])\n",
    "for i,p in enumerate(paramList[:nCosParams]):\n",
    "    Bi = sigmas_Sum2[i]\n",
    "    Zi = sigmas_S4book_PlanckS4[i]\n",
    "    diff = (Bi-Zi)/Zi *100\n",
    "    t.add_row([p,\"%5.3g\" % Zi,\"%5.3g\" % Bi,\"%5.3g\" % diff] )\n",
    "print t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+-----------+----------------+---------------+--------+\n",
      "| Parameter | S4+Planck_book | S4+Planck_Bye | % diff |\n",
      "+-----------+----------------+---------------+--------+\n",
      "| 100 ombh2 |     0.003      |    0.00307    |  2.33  |\n",
      "|   omch2   |     0.0006     |    0.000664   |  10.7  |\n",
      "|     H0    |      0.24      |     0.253     |  5.42  |\n",
      "|  10^9 As  |     0.021      |     0.0204    | -2.86  |\n",
      "|     ns    |     0.002      |    0.00203    |   1.5  |\n",
      "|    tau    |     0.006      |    0.00557    | -7.17  |\n",
      "+-----------+----------------+---------------+--------+\n"
     ]
    }
   ],
   "source": [
    "# compare to Byeonghee's result\n",
    "#from prettytable import PrettyTable\n",
    "#sigmas_S4book_Planck = [0.017,0.0014,0.7,0.039,0.004,0.01]\n",
    "#sigmas_S4book_PlanckS4 = [0.003,0.0006,0.24,0.021,0.002,0.006]\n",
    "#paramNamesB = ['100 ombh2','omch2','H0','10^9 As','ns','tau']\n",
    "\n",
    "# from Bye:\n",
    "\"\"\"\n",
    "{'As': 2.0374125419013829e-11,\n",
    " 'H0': 0.25302588229880729,\n",
    " 'ns': 0.0020335606322583613,\n",
    " 'ombh2': 3.0755699027795784e-05,\n",
    " 'omch2': 0.00066427378612649382,\n",
    " 'tau': 0.0055703760229570097}\n",
    " \"\"\"\n",
    "sigmas_Byeonghee = [0.00307,0.000664,0.253,0.0204,0.00203,0.00557]\n",
    "\n",
    "t = PrettyTable(['Parameter','S4+Planck_book', 'S4+Planck_Bye','% diff'])\n",
    "for i,p in enumerate(paramList[:nCosParams]):\n",
    "    Bi = sigmas_Byeonghee[i]\n",
    "    Zi = sigmas_S4book_PlanckS4[i]\n",
    "    diff = (Bi-Zi)/Zi *100\n",
    "    t.add_row([p,\"%5.3g\" % Zi,\"%5.3g\" % Bi,\"%5.3g\" % diff] )\n",
    "print t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# now compare my result directly with Bye's"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "compare my result directly to Bye's\n",
      "+-----------+-------------+---------------+--------+\n",
      "| Parameter | S4+Planck_Z | S4+Planck_Bye | % diff |\n",
      "+-----------+-------------+---------------+--------+\n",
      "| 100 ombh2 |   0.00229   |    0.00307    |  33.8  |\n",
      "|   omch2   |   0.00055   |    0.000664   |  20.7  |\n",
      "|     H0    |    0.219    |     0.253     |  15.5  |\n",
      "|  10^9 As  |    0.0205   |     0.0204    | -0.443 |\n",
      "|     ns    |   0.00172   |    0.00203    |  17.8  |\n",
      "|    tau    |   0.00543   |    0.00557    |  2.55  |\n",
      "+-----------+-------------+---------------+--------+\n"
     ]
    }
   ],
   "source": [
    "print 'compare my result directly to Bye\\'s'\n",
    "t = PrettyTable(['Parameter','S4+Planck_Z', 'S4+Planck_Bye','% diff'])\n",
    "for i,p in enumerate(paramList[:nCosParams]):\n",
    "    Bi = sigmas_Byeonghee[i]\n",
    "    Zi = sigmas_Sum2[i]\n",
    "    diff = (Bi-Zi)/Zi *100\n",
    "    t.add_row([p,\"%5.3g\" % Zi,\"%5.3g\" % Bi,\"%5.3g\" % diff] )\n",
    "print t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Next try with CosmicFish results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 4.30400590e+09 -6.59605228e+08 -1.95219475e+08  7.37773638e+07\n",
      "   9.94588306e+07 -1.42641941e+08]\n",
      " [-6.59605228e+08  2.06770781e+08  5.40309866e+07 -1.12167620e+07\n",
      "  -1.24332342e+07  2.26744356e+07]\n",
      " [-1.95219475e+08  5.40309866e+07  1.54081445e+07 -3.36576138e+06\n",
      "  -3.86681499e+06  6.28546909e+06]\n",
      " [ 7.37773638e+07 -1.12167620e+07 -3.36576138e+06  3.44653161e+06\n",
      "   3.14441968e+06 -6.42453437e+06]\n",
      " [ 9.94588306e+07 -1.24332342e+07 -3.86681499e+06  3.14441968e+06\n",
      "   3.66513392e+06 -5.99246391e+06]\n",
      " [-1.42641941e+08  2.26744356e+07  6.28546909e+06 -6.42453437e+06\n",
      "  -5.99246391e+06  1.22193648e+07]]\n"
     ]
    }
   ],
   "source": [
    "cfDir = 'cosmicfish_results/'\n",
    "FijFilename = '3_CMBS4_TEB_fisher_matrix_cls.dat'  # TEB only\n",
    "\n",
    "# The parameters of this Fisher matrix are:\n",
    "#\n",
    "#    1    omegabh2    \\Omega_b h^2       0.2220000000000000E-01\n",
    "#    2    omegach2    \\Omega_c h^2       0.1200000000000000E+00\n",
    "#    3    h    h       0.6900000000000001E+00\n",
    "#    4    logA    {\\rm{ln}}(10^{10} A_s)       0.3091042453358316E+01\n",
    "#    5    ns    n_s       0.9660000000000000E+00\n",
    "#    6    tau    \\tau       0.6000000000000000E-01\n",
    "#\n",
    "\n",
    "Fij_CF = np.loadtxt(cfDir+FijFilename)\n",
    "print Fij_CF"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Transform to $10^9 As$ via Jacobian"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [],
   "source": [
    "# only one partial derivative needed\n",
    "# d( ln(10^10*As)) / d(10^9*As) = 1/(10^9*As)*d(10*10^9*As)/d(10^9*As) = 1/(10^9*As)\n",
    "TenToTheNineAs_fid = 2.2\n",
    "AsIndex = 3\n",
    "Fij_CF[:,AsIndex] *= 1/TenToTheNineAs_fid\n",
    "Fij_CF[AsIndex,:] *= 1/TenToTheNineAs_fid"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Transform to $H_0$ via Jacobian"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [],
   "source": [
    "# d(h)/d(100h) = 1/100\n",
    "HIndex = 2\n",
    "Fij_CF[:,HIndex] /= 100\n",
    "Fij_CF[HIndex,:] /= 100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[3.83938345e-05 8.56005710e-04 3.23719005e-01 2.75374827e-02\n",
      " 2.48307919e-03 7.44515332e-03]\n"
     ]
    }
   ],
   "source": [
    "# get sigmas\n",
    "getSigmas = lambda Fij: np.sqrt(np.diag(np.linalg.inv(Fij)))\n",
    "sigmas_CosmicFish = getSigmas(Fij_CF)\n",
    "print sigmas_CosmicFish"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [],
   "source": [
    "# combine with others\n",
    "# fsky is in denominator of covariance matrix, putting it in numerator of Fij matrix\n",
    "#Fij_Sum = (Fij_S4_TE+Fij_plikHM_lowTEB)*fsky_S4 + Fij_plikHM_TTTEEE_lowTEB*fsky_Pl\n",
    "#Fij_Sum_Pl = Fij_plikHM_TTTEEE_lowTEB*(fsky_S4+fsky_Pl)\n",
    "\n",
    "# fsky already applied to Fij_CF\n",
    "Fij_Sum_CF = Fij_CF + Fij_plikHM_TTTEEE_lowTEB*fsky_Pl\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [],
   "source": [
    "# add sigma = 0.01 prior on tau\n",
    "tauIndex = 5\n",
    "tauPrior = 1/0.01**2\n",
    "#Fij_Sum[tauIndex,tauIndex] += tauPrior\n",
    "#Fij_Sum_Pl[tauIndex,tauIndex] += tauPrior\n",
    "\n",
    "Fij_Sum_CF[tauIndex,tauIndex] += tauPrior\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Transform to $100 \\Omega_b h^2$ via Jacobian"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [],
   "source": [
    "# d(ombh2)/d(100obmh2) = 1/100\n",
    "ObIndex = 0\n",
    "Fij_CF[:,ObIndex] /= 100\n",
    "Fij_CF[ObIndex,:] /= 100\n",
    "\n",
    "Fij_Sum_CF[:,ObIndex] /= 100\n",
    "Fij_Sum_CF[ObIndex,:] /= 100\n",
    "\n",
    "Fij_kk2[:,ObIndex] /= 100\n",
    "Fij_kk2[ObIndex,:] /= 100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+-----------+-------------+----------+--------+\n",
      "| Parameter | Planck_book | Planck_Z | % diff |\n",
      "+-----------+-------------+----------+--------+\n",
      "| 100 ombh2 |    0.017    |  0.0174  |  2.61  |\n",
      "|   omch2   |    0.0014   | 0.00171  |  21.8  |\n",
      "|     H0    |      0.7    |  0.759   |  8.38  |\n",
      "|  10^9 As  |    0.039    |  0.0393  |  0.69  |\n",
      "|     ns    |    0.004    |  0.0048  |  20.1  |\n",
      "|    tau    |     0.01    |  0.009   | -9.98  |\n",
      "+-----------+-------------+----------+--------+\n",
      "+-----------+----------------+----------------+--------+\n",
      "| Parameter | S4+Planck_book | S4_CF+Planck_Z | % diff |\n",
      "+-----------+----------------+----------------+--------+\n",
      "| 100 ombh2 |     0.003      |    0.00372     |  23.9  |\n",
      "|   omch2   |     0.0006     |    0.000675    |  12.5  |\n",
      "|     H0    |      0.24      |     0.256      |  6.84  |\n",
      "|  10^9 As  |     0.021      |     0.0214     |   1.8  |\n",
      "|     ns    |     0.002      |    0.00213     |   6.3  |\n",
      "|    tau    |     0.006      |    0.00575     | -4.15  |\n",
      "+-----------+----------------+----------------+--------+\n",
      "This from CosmicFish with no kappa\n"
     ]
    }
   ],
   "source": [
    "# Show the sigmas\n",
    "# get sigmas\n",
    "#sigmas_Sum = getSigmas(Fij_Sum)\n",
    "#sigmas_Sum_Pl = getSigmas(Fij_Sum_Pl)\n",
    "sigmas_Sum_CF = getSigmas(Fij_Sum_CF)\n",
    "\n",
    "t = PrettyTable(['Parameter','Planck_book', 'Planck_Z','% diff'])\n",
    "for i,p in enumerate(paramList[:nCosParams]):\n",
    "    Bi = sigmas_Sum_Pl[i]\n",
    "    Zi = sigmas_S4book_Planck[i]\n",
    "    diff = (Bi-Zi)/Zi *100\n",
    "    t.add_row([p,\"%5.3g\" % Zi,\"%5.3g\" % Bi,\"%5.3g\" % diff] )\n",
    "print t\n",
    "\n",
    "t = PrettyTable(['Parameter','S4+Planck_book', 'S4_CF+Planck_Z','% diff'])\n",
    "for i,p in enumerate(paramList[:nCosParams]):\n",
    "    Bi = sigmas_Sum_CF[i]\n",
    "    Zi = sigmas_S4book_PlanckS4[i]\n",
    "    diff = (Bi-Zi)/Zi *100\n",
    "    t.add_row([p,\"%5.3g\" % Zi,\"%5.3g\" % Bi,\"%5.3g\" % diff] )\n",
    "print t\n",
    "print \"This from CosmicFish with no kappa\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {},
   "outputs": [],
   "source": [
    "# now add in kappa\n",
    "Fij_Sum_CF_k = Fij_Sum_CF + Fij_kk2*fsky_S4\n",
    "sigmas_Sum_CF_k = getSigmas(Fij_Sum_CF_k)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+-----------+-------------+----------+--------+\n",
      "| Parameter | Planck_book | Planck_Z | % diff |\n",
      "+-----------+-------------+----------+--------+\n",
      "| 100 ombh2 |    0.017    |  0.0174  |  2.61  |\n",
      "|   omch2   |    0.0014   | 0.00171  |  21.8  |\n",
      "|     H0    |      0.7    |  0.759   |  8.38  |\n",
      "|  10^9 As  |    0.039    |  0.0393  |  0.69  |\n",
      "|     ns    |    0.004    |  0.0048  |  20.1  |\n",
      "|    tau    |     0.01    |  0.009   | -9.98  |\n",
      "+-----------+-------------+----------+--------+\n",
      "+-----------+----------------+---------------------+--------+\n",
      "| Parameter | S4+Planck_book | S4_CF+Planck_Z+kk_Z | % diff |\n",
      "+-----------+----------------+---------------------+--------+\n",
      "| 100 ombh2 |     0.003      |       0.00323       |  7.53  |\n",
      "|   omch2   |     0.0006     |       0.000567      | -5.52  |\n",
      "|     H0    |      0.24      |        0.212        | -11.5  |\n",
      "|  10^9 As  |     0.021      |        0.0201       | -4.08  |\n",
      "|     ns    |     0.002      |       0.00196       | -2.22  |\n",
      "|    tau    |     0.006      |       0.00546       | -8.92  |\n",
      "+-----------+----------------+---------------------+--------+\n"
     ]
    }
   ],
   "source": [
    "t = PrettyTable(['Parameter','Planck_book', 'Planck_Z','% diff'])\n",
    "for i,p in enumerate(paramList[:nCosParams]):\n",
    "    Bi = sigmas_Sum_Pl[i]\n",
    "    Zi = sigmas_S4book_Planck[i]\n",
    "    diff = (Bi-Zi)/Zi *100\n",
    "    t.add_row([p,\"%5.3g\" % Zi,\"%5.3g\" % Bi,\"%5.3g\" % diff] )\n",
    "print t\n",
    "\n",
    "t = PrettyTable(['Parameter','S4+Planck_book', 'S4_CF+Planck_Z+kk_Z','% diff'])\n",
    "for i,p in enumerate(paramList[:nCosParams]):\n",
    "    Bi = sigmas_Sum_CF_k[i]\n",
    "    Zi = sigmas_S4book_PlanckS4[i]\n",
    "    diff = (Bi-Zi)/Zi *100\n",
    "    t.add_row([p,\"%5.3g\" % Zi,\"%5.3g\" % Bi,\"%5.3g\" % diff] )\n",
    "print t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}