dH0_vs_N_v2.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 287
    },
    "colab_type": "code",
    "id": "EZt-5M22LXdd",
    "outputId": "eebea36e-625b-40d0-f603-a0e846de7135"
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "%matplotlib notebook\n",
    "import matplotlib.pyplot as plt\n",
    "golden_ratio=(1.+np.sqrt(5))/2.\n",
    "ysize = 3.5 #inches\n",
    "plt.rcParams['figure.figsize']=(ysize * golden_ratio, ysize)\n",
    "try:\n",
    "    from matplotlib import rc\n",
    "    rc('font',**{'family':'sans','serif':['Times New Roman']})\n",
    "    rc('text', usetex=True)\n",
    "except:\n",
    "    print(\"Not using Latex\")\n",
    "    \n",
    "plt.rcParams[\"font.family\"] = \"Times New Roman\"\n",
    "\n",
    "# from Belgacem+2019 with 15 short GRB / yr (45 in 3 yrs)\n",
    "\n",
    "def dH0(Nb,P,N):\n",
    "    # Nb = number of short GRB from Belgacem et al. 2019 that provide precision P on H0\n",
    "    # P = precision on H0 given the number of short GRB Nb\n",
    "    # N = number of short GRB that we want to test\n",
    "    return P*(np.sqrt(Nb)/np.sqrt(N))\n",
    "\n",
    "\n",
    "# === INPUT DATA ===\n",
    "\n",
    "# Belgacem values\n",
    "xb_ET=150\n",
    "P_ET=0.41\n",
    "xb_ET_CE_CE=255\n",
    "P_ET_CE_CE=0.235\n",
    "\n",
    "#Zhao and van den broek et al\n",
    "xz=1000\n",
    "Pz=0.33\n",
    "\n",
    "# SGRB with Theseus (1/3, 1/2, 2/3 and tot) \n",
    "#pz=0.63   #with z (63%)\n",
    "pz=0.73   #with z (25%+38%) + 0.75*(13%)\n",
    "#pz=1.0\n",
    "\n",
    "Non=12   # onaxis short GRB with THESEUS\n",
    "\n",
    "Nonoff=Non*1.768   # onaxis + offaxis short GRB with THESEUS\n",
    "\n",
    "#N_ET_1yr=round(2.3*pz,1)\n",
    "#N_ET_4yr=round(8*pz,1)\n",
    "N_ET_1yr=round(0.63*Nonoff*pz,1)   # percentage detected with ET for on+offaxis view\n",
    "N_ET_4yr=round(N_ET_1yr*3.45,1)\n",
    "\n",
    "#N_ET_CE_1yr=round(6*pz,1)\n",
    "#N_ET_CE_4yr=round(21*pz,1)\n",
    "N_ET_CE_1yr=round(0.76*Nonoff*pz,1)\n",
    "N_ET_CE_4yr=round(N_ET_CE_1yr*3.45,1)\n",
    "N_ET_2CE_1yr=round(0.83*Nonoff*pz,1)\n",
    "N_ET_2CE_4yr=round(N_ET_2CE_1yr*3.45,1)\n",
    "\n",
    "\n",
    "Ntot=20\n",
    "Nnew=40\n",
    "\n",
    "\n",
    "\n",
    "# === PLOT ===\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x=np.arange(300)+1\n",
    "\n",
    "# compute Belgacem accuracy \n",
    "y_ET=dH0(xb_ET,P_ET,x)\n",
    "y_ET_CE_CE=dH0(xb_ET_CE_CE,P_ET_CE_CE,x)\n",
    "yz=dH0(xz,Pz,x)\n",
    "\n",
    "# compute true accuracy\n",
    "y_ET_4yr=dH0(xb_ET,P_ET,x*0+N_ET_4yr)\n",
    "y_ET_1yr=dH0(xb_ET,P_ET,x*0+N_ET_1yr)\n",
    "\n",
    "y_ET_CE_CE_1yr=dH0(xb_ET_CE_CE,P_ET_CE_CE,x*0+N_ET_2CE_1yr)\n",
    "y_ET_CE_CE_4yr=dH0(xb_ET_CE_CE,P_ET_CE_CE,x*0+N_ET_2CE_4yr)\n",
    "\n",
    "\n",
    "# theseus short GRB detection with IRT and from XGIS only\n",
    "old=[12,14,30,36]\n",
    "new=[14,17,36,43]\n",
    "\n",
    "for i in old:\n",
    "    print('With '+str(i)+' short GRB:')\n",
    "    print('Old H0 accuracy with THESEUS+ET and THESEUS+ET+CE+CE')\n",
    "    print(dH0(xb_ET,P_ET,x*0+i)[0])\n",
    "    print(dH0(xb_ET_CE_CE,P_ET_CE_CE,x*0+i)[0])\n",
    "for i in new:\n",
    "    print('With '+str(i)+' short GRB:')\n",
    "    print('New H0 accuracy with THESEUS+ET')\n",
    "    print(dH0(xb_ET,P_ET,x*0+i)[0])    \n",
    "    print(dH0(xb_ET_CE_CE,P_ET_CE_CE,x*0+i)[0])\n",
    "\n",
    "\n",
    "plt.xlabel('Number of short GRB with measured redshift')\n",
    "\n",
    "if plt.rcParams[\"text.usetex\"]:\n",
    "    plt.ylabel('$H_0$ 1$\\sigma$ precision (\\%)')\n",
    "else:\n",
    "    plt.ylabel('$H_0$ 1$\\sigma$ precision (%)')\n",
    "\n",
    "\n",
    "# red dots - Belgacem H0 accuracy\n",
    "yb_ET=dH0(xb_ET,P_ET,xb_ET)\n",
    "yb_ET_CE_CE=dH0(xb_ET_CE_CE,P_ET_CE_CE,xb_ET_CE_CE)\n",
    "#plt.loglog(xb_ET,yb_ET,'ro',linewidth=2, markersize=8, label='Belgacem+2019 (10yr)')\n",
    "#plt.loglog(x,y_ET,'--r',label='1/SQRT(N)')\n",
    "plt.semilogx(xb_ET,yb_ET,'ro',linewidth=2, markersize=8, label='Belgacem+2019 (ET,10yr)')\n",
    "plt.semilogx(xb_ET_CE_CE,yb_ET_CE_CE,'bD',linewidth=2, markersize=8, label='Belgacem+2019 (ET+2CE,10yr)')\n",
    "\n",
    "# red crosses - ET Real accuracy\n",
    "#plt.semilogx(Ntrue,y_ET_true[0],'kx',linewidth=2, markersize=8,label='THESEUS '+str(Npesss)+','+str(Npess)+','+str(Ntrue)+','+str(Ntot)+','+str(Nnew)+','+' SGRB')\n",
    "#plt.semilogx(Ntrue,y_ET_true[0],'kx',linewidth=2, markersize=8,label='THESEUS '+str(Npesss)+','+str(Ntrue)+','+' SGRB+GW')\n",
    "#plt.semilogx(Ntrue,y_ET_true[0],'kx',linewidth=2, markersize=8,label='THESEUS '+str(Npesss)+','+str(Ntrue)+','+' SGRB+GW')\n",
    "plt.semilogx(N_ET_1yr,y_ET_1yr[0],'rx',linewidth=8, markersize=10)\n",
    "#plt.arrow(N_ET_1yr,y_ET_1yr[0],0,-0.1,width=0.1)\n",
    "#plt.semilogx(Ntot,y_ET_tot[0],'kx',linewidth=2, markersize=8)\n",
    "#plt.semilogx(Npess,y_ET_pess[0],'kx',linewidth=2, markersize=8)\n",
    "#plt.semilogx(N_ET_4yr,y_ET_4yr[0],'rx',linewidth=8, markersize=10, label='THESEUS+ET (1-3.45yr): '+str(N_ET_1yr)+'-'+str(N_ET_4yr)+' SGRBs+GW with z')\n",
    "plt.semilogx(N_ET_4yr,y_ET_4yr[0],'rx',linewidth=8, markersize=10, label='THESEUS+ET (1-3.45yr): '+str(N_ET_1yr)+'-'+str(N_ET_4yr)+' SGRBs+GW with z')\n",
    "#plt.arrow(N_ET_4yr,y_ET_4yr[0],1.5,-0.3,color='r',width=0.04,head_width=0.5)\n",
    "#plt.annotate('text',(N_ET_4yr,y_ET_4yr[0]))\n",
    "#plt.semilogx(Nnew,y_ET_new[0],'kx',linewidth=2, markersize=8)\n",
    "\n",
    "plt.semilogx(N_ET_2CE_1yr,y_ET_CE_CE_1yr[0],'b+',linewidth=8, markersize=10, \n",
    "             label='THESEUS+ET+2CE (1-3.45yr): '+str(N_ET_2CE_1yr)+'-'+str(N_ET_CE_4yr)+' SGRBs+GW with z')\n",
    "#plt.semilogx(Ngoal,y_ET_goal[0],'kx',linewidth=2, markersize=8)\n",
    "#plt.semilogx(Ntot,y_ET_CE_CE_tot[0],'kx',linewidth=2, markersize=8)\n",
    "#plt.semilogx(Npess,y_ET_CE_CE_pess[0],'kx',linewidth=2, markersize=8)\n",
    "#plt.semilogx(Npesss,y_ET_CE_CE_pesss[0],'kx',linewidth=2, markersize=8)\n",
    "plt.semilogx(N_ET_2CE_4yr,y_ET_CE_CE_4yr[0],'b+',linewidth=8, markersize=10)\n",
    "\n",
    "# 1/SQRT(N) line\n",
    "plt.semilogx(x,y_ET,':k',label='$N^{-1/2}$')\n",
    "plt.semilogx(x,y_ET_CE_CE,':k')\n",
    "#plt.loglog(xb_ET_CE_CE,yb_ET_CE_CE,'ro',linewidth=2, markersize=8)\n",
    "\n",
    "# realistic accuracy\n",
    "#plt.loglog(x,y_ET, label='THESEUS+ET')\n",
    "#plt.loglog(x,y_ET_true,'k', label='THESEUS+ET with 12 SGRB in 3 yr')\n",
    "#plt.loglog(x,y_ET_goal,'--k', label='THESEUS+ET with 30 SGRB in 3 yr')\n",
    "#plt.semilogx(x,y_ET_true,color='silver', label='THESEUS+ET with '+str(Npess)+','+str(Ntrue)+','+str(Ntot)+','+str(Ngoal)+' SGRB')\n",
    "#plt.semilogx(x,y_ET_true,color='silver', label='THESEUS+ET with '+str(Npess)+','+str(Ntrue)+','+str(Ntot)+','+' SGRB')\n",
    "#plt.semilogx(x,y_ET_CE_CE_true,color='silver', label='THESEUS+ET+CE+CE with '+str(Npess)+','+str(Ntrue)+','+str(Ntot)+','+' SGRB')\n",
    "\n",
    "# plot horizontal lines\n",
    "#plt.semilogx(x,y_ET_tot,color='silver')\n",
    "#plt.semilogx(x,y_ET_CE_CE_tot,color='silver')\n",
    "#plt.semilogx(x,y_ET_pess,color='silver')\n",
    "#plt.semilogx(x,y_ET_CE_CE_pess,color='silver')\n",
    "\n",
    "\n",
    "# plot vertical lines\n",
    "#plt.axvline(x=Ntrue,color='silver')\n",
    "#plt.axvline(x=Ntot,color='silver')\n",
    "#plt.axvline(x=Npess,color='silver')\n",
    "\n",
    "yplanck=x*0+0.7\n",
    "#plt.loglog(x,yplanck,'m', label='Planck')\n",
    "plt.semilogx(x,yplanck,'m', label='Planck')\n",
    "ySNIa=x*0+2.4\n",
    "#plt.loglog(x,ySNIa,'c', label='SNIa')\n",
    "plt.semilogx(x,ySNIa,'c', label='SNIa')\n",
    "y1=x*0+1.\n",
    "#plt.loglog(x,y1,'--r', label='~1% with 15-25 sGRB')\n",
    "\n",
    "\n",
    "#plt.legend(loc=1,fontsize=(8))\n",
    "plt.legend(loc='upper right', bbox_to_anchor=(1.0, 1.0),fontsize=(7))\n",
    "plt.xlim(1.0,380,10)\n",
    "plt.ylim(0,4.9,0.5)\n",
    "#plt.yticks(np.arange(0,2,0.2))\n",
    "plt.grid()\n",
    "\n",
    "#plt.savefig('H0_plot.tiff', dpi=600, bbox_inches = \"tight\")\n",
    "plt.savefig('H0_plot.png', dpi=600, bbox_inches = \"tight\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.font_manager\n",
    "from IPython.core.display import HTML\n",
    "\n",
    "def make_html(fontname):\n",
    "    return \"<p>{font}: <span style='font-family:{font}; font-size: 24px;'>{font}</p>\".format(font=fontname)\n",
    "\n",
    "code = \"\\n\".join([make_html(font) for font in sorted(set([f.name for f in matplotlib.font_manager.fontManager.ttflist]))])\n",
    "\n",
    "HTML(\"<div style='column-count: 2;'>{}</div>\".format(code))"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "name": "dH0_vs_N.ipynb",
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}