make-hisograms_v2.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "%matplotlib notebook\n",
    "import matplotlib.pyplot as plt\n",
    "\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(\"Latex not available\")\n",
    "plt.rcParams[\"font.family\"] = \"Times New Roman\"\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "binstep=0.1\n",
    "\n",
    "# read output from read_zlim_vs_theta.py\n",
    "#dbtheta=pd.read_csv('thetabin_zbin01.log')\n",
    "dbtheta=pd.read_csv('thetabin_zbin01.log')\n",
    "#dbtheta=pd.read_csv('thetabin_zbin002.log')\n",
    "\n",
    "thetaview=dbtheta['theta']\n",
    "zview=dbtheta['z']\n",
    "\n",
    "\n",
    "# read output from BNS simulations by Jan Delta_z=0.3\n",
    "\n",
    "dbeff5deg=pd.read_csv('Efficiency_BNS_SNR8.0_deg5_z.txt',header=None,sep=' ')\n",
    "Pbns_ET2CE5deg_03=dbeff5deg[1]\n",
    "Pbns_ETCE5deg_03=dbeff5deg[2]\n",
    "Pbns_ET5deg_03=dbeff5deg[3]\n",
    "\n",
    "dbeff10deg=pd.read_csv('Efficiency_BNS_SNR8.0_deg10_z.txt',header=None,sep=' ')\n",
    "zeff=dbeff10deg[0].values\n",
    "Pbns_ET2CE10deg_03=dbeff10deg[1]\n",
    "Pbns_ETCE10deg_03=dbeff10deg[2]\n",
    "Pbns_ET10deg_03=dbeff10deg[3]\n",
    "\n",
    "dbeff=pd.read_csv('Efficiency_BNS_SNR8.0_z.txt',header=None,sep=' ')\n",
    "Pbns_ET2CE_all_03=dbeff[1]\n",
    "Pbns_ETCE_all_03=dbeff[2]\n",
    "Pbns_ET_all_03=dbeff[3]\n",
    "\n",
    "Pbns_ET2CE_03=[]\n",
    "Pbns_ETCE_03=[]\n",
    "Pbns_ET_03=[]\n",
    "for i in range(0,len(zeff)):\n",
    "    if zeff[i]<0.3:\n",
    "        Pbns_ET2CE_03.append(Pbns_ET2CE_all_03[i])\n",
    "        Pbns_ETCE_03.append(Pbns_ETCE_all_03[i])\n",
    "        Pbns_ET_03.append(Pbns_ET_all_03[i])\n",
    "    if zeff[i]>=0.3:\n",
    "        Pbns_ET2CE_03.append(Pbns_ET2CE10deg_03[i])\n",
    "        Pbns_ETCE_03.append(Pbns_ETCE10deg_03[i])\n",
    "        Pbns_ET_03.append(Pbns_ET10deg_03[i])\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# read output from MOS (i.e. the redshift of simulated short GRBs)\n",
    "dfx=pd.read_csv('shortGRBsDetectedBySXIXGIS.txt',header=None)\n",
    "zx=dfx[0]\n",
    "\n",
    "\n",
    "# define z grid\n",
    "zmax=20\n",
    "\n",
    "l=int((zmax-0)/binstep)\n",
    "\n",
    "bins=[]\n",
    "for i in range(0,l):\n",
    "    bins.append(round(binstep*i,1))\n",
    "    #bins.append(round(binstep*i,2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#fig = plt.figure()\n",
    "\n",
    "hist_x_density=np.histogram(zx,bins=bins,density=True)\n",
    "\n",
    "#hist_x_density=plt.hist(zx,bins=bins,density=True,alpha=0)\n",
    "hist_x_number=np.histogram(zx,bins=bins)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# distribuzione delle simulazioni (400)\n",
    "#Fx=hist_x_density[0]\n",
    "Nx0=hist_x_number[0]\n",
    "zbinx=hist_x_density[1]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# deg up to which a grb at zbin can be detected offaxis\n",
    "thetajet=4.0\n",
    "theta=np.asarray(Nx0)*0+thetajet\n",
    "\n",
    "#for i in range(0,32):\n",
    "for i in range(0,17):\n",
    "#for i in range(0,90):\n",
    "    theta[i]=thetaview[i]\n",
    "\n",
    "thetaoffaxisvsz=theta*np.pi/180\n",
    "thetonaxis=thetajet*np.pi/180\n",
    "increasefact=(1-np.cos(thetaoffaxisvsz))/(1-np.cos(thetonaxis))\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Fx0=Nx0/sum(Nx0)\n",
    "Nx=Nx0*increasefact\n",
    "Fx=Nx/sum(Nx)\n",
    "\n",
    "\n",
    "# efficienza di detection in bin di ampiezza 0.02\n",
    "#Pbns_ETCE=np.repeat(Pbns_ETCE_03,15)\n",
    "#Pbns_ET=np.repeat(Pbns_ET_03,15)\n",
    "\n",
    "# efficienza di detection in bin di ampiezza 0.1\n",
    "Pbns_ET2CE=np.repeat(Pbns_ET2CE_03,3)\n",
    "Pbns_ETCE=np.repeat(Pbns_ETCE_03,3)\n",
    "Pbns_ET=np.repeat(Pbns_ET_03,3)\n",
    "\n",
    "# efficienza di detection in bin di ampiezza 0.05 (ogni valore della prob. è ripetuto 3 volte)\n",
    "#Pbns_ETCE=(np.asarray([0.115,0.115,0.115,0.115,0.115,0.115,0.125,0.125,0.125,0.125,0.125,0.125,0.118,0.118,0.118,0.118,0.118,0.118,0.08,0.08,0.08,0.08,0.08,0.08,0.055,0.055,0.055,0.055,0.055,0.055,0.04,0.04,0.04,0.04,0.04,0.04,0.03,0.03,0.03,0.03,0.03,0.03,0.02,0.02,0.02,0.02,0.02,0.02,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01,0.01])/0.125)\n",
    "#Pbns_ET=(np.asarray([0.115,0.115,0.115,0.115,0.115,0.115,0.12,0.12,0.12,0.12,0.12,0.12,0.07,0.07,0.07,0.07,0.07,0.07,0.025,0.025,0.025,0.025,0.025,0.025,0.021,0.021,0.021,0.021,0.021,0.021,0.005,0.005,0.005,0.005,0.005,0.005,0.01,0.01,0.01,0.01,0.01,0.01,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0])/0.125)\n",
    "\n",
    "\n",
    "\n",
    "# percentuale vista anche con GW\n",
    "#P_TH_ETCE=PGRBinz_x*Pbns_ETCE\n",
    "#P_TH_ET=PGRBinz_x*Pbns_ET\n",
    "leff=len(Pbns_ET)\n",
    "P_TH_ET2CE=Fx[0:leff]*Pbns_ET2CE\n",
    "P_TH_ETCE=Fx[0:leff]*Pbns_ETCE\n",
    "P_TH_ET=Fx[0:leff]*Pbns_ET\n",
    "#P_TH_ETCE=Fx[0:150]*Pbns_ETCE\n",
    "#P_TH_ET=Fx[0:150]*Pbns_ET"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.bar(zbinx[0:-1],100*Fx,width=binstep,align='edge',log=False,linewidth=2,color=\"None\",edgecolor='cornflowerblue',alpha=0.9,label='THESEUS short GRB')\n",
    "\n",
    "plt.bar(zbinx[0:leff],100*P_TH_ET2CE,width=binstep,log=False,align='edge',hatch='-',color='violet',alpha=0.5,label='THESEUS+ET+2CE short GRB')\n",
    "plt.bar(zbinx[0:leff],100*P_TH_ETCE,width=binstep,log=False,align='edge',hatch='\\\\',color='magenta',alpha=0.5,label='THESEUS+ET+CE short GRB')\n",
    "plt.bar(zbinx[0:leff],100*P_TH_ET,width=binstep,log=False,align='edge',hatch='/',color='lime',alpha=0.5,label='THESEUS+ET short GRB')\n",
    "#plt.bar(zbinx[0:150],100*P_TH_ETCE,width=binstep,log=False,align='edge',hatch='\\\\',color='magenta',alpha=0.5,label='THESEUS+ET+CE on+off axis short GRB')\n",
    "#plt.bar(zbinx[0:150],100*P_TH_ET,width=binstep,log=False,align='edge',hatch='/',color='lime',alpha=0.5,label='THESEUS+ET on+offaxis short GRB')\n",
    "#plt.bar((x1+x2)/2.,P_TH_offaxis_ET,width=(x2-x1),align='center',linewidth=5,alpha=0.2,edgecolor='green',label='THESEUS+ET short GRB including off-axis ')\n",
    "#plt.bar((x1+x2)/2.,P_TH_offaxis_ETCE,width=(x2-x1),align='center',linewidth=5,alpha=0.3,edgecolor='magenta',label='THESEUS+ET+CE short GRB including off-axis ')\n",
    "#plt.xlim(0.,13)\n",
    "plt.xlim(0.,6)\n",
    "plt.ylim(0,7.5)\n",
    "plt.title('Aligned+misaligned short GRBs')\n",
    "plt.xlabel('redshift')\n",
    "plt.ylabel(' short GRB [$\\%$ in each redshift bin]')\n",
    "\n",
    "#plt.legend(fontsize='small')\n",
    "plt.legend()\n",
    "\n",
    "\n",
    "#plt.savefig('SGRB_hist_ET_CE_offaxis_bin01.pdf')\n",
    "plt.savefig('SGRB_hist_ET_CE_offaxis_bin01.png',dpi=600, bbox_inches = \"tight\")\n",
    "#plt.savefig('SGRB_hist_ET_2CE_offaxis_bin01.tif',dpi=600, bbox_inches = \"tight\")\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "k=sum(Nx)/sum(Nx0)\n",
    "\n",
    "print('Independently on GW det. efficiency, the total number of on- and off-axis GBR is ',k,' times larger than on-axis only GRBs (i.e. 12/yr on-axis and ',12*k,'/yr on- and off-axis)')\n",
    "print('')\n",
    "\n",
    "print('The fraction of on- + off-axis SGRB detected with GW interferometers are:')\n",
    "print('THESEUS+ET:',100*sum(P_TH_ET),'%')\n",
    "print('THESEUS+ET+CE:',100*sum(P_TH_ETCE),'%')\n",
    "print('THESEUS+ET+2CE:',100*sum(P_TH_ET2CE),'%')\n",
    "\n",
    "print('')\n",
    "print('Total GRB THESEUS (on+off axis) in 1(3.45)yr:', 12*k,'(',12*k*3.45,')')\n",
    "print('Total GRB THESEUS+ET (on+off axis) in 1(3.45)yr:', 12*k*sum(P_TH_ET),'(',12*k*3.45*sum(P_TH_ET),')')\n",
    "print('Total GRB THESEUS+ET+CE (on+off axis) in 1(3.45)yr:', 12*k*sum(P_TH_ETCE),'(',12*k*3.45*sum(P_TH_ETCE),')')\n",
    "print('Total GRB THESEUS+ET+2CE (on+off axis) in 1(3.45)yr:', 12*k*sum(P_TH_ET2CE),'(',12*k*3.45*sum(P_TH_ET2CE),')')\n",
    "\n",
    "#up to z=3.0\n",
    "#k03=sum(Nx[0:10])/sum(Nx0[0:10])\n",
    "# GRB with ET too\n",
    "\n",
    "\n",
    "# only if bin <=0.1\n",
    "print('')\n",
    "\n",
    "#print('Total GRB THESEUS+2G (at z<0.1) in 3.45 yr: ',40*k*Fx[0])\n",
    "#print('Total GRB THESEUS+2G (at z<0.08) in 3.45 yr: ',sum(40*k*Fx[0:4]))\n",
    "#print('Total GRB THESEUS+2G (at z<0.06) in 3.45 yr: ',sum(40*k*Fx[0:3]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# ---- PLOT ON-AXIS SGRB DISTRIBUTION ------\n",
    "\n",
    "dfr=pd.read_csv('shortGRBsDetectedByIRT.txt',header=None)\n",
    "#dfx=pd.read_csv('shortGRBsDetectedBySXIXGIS.txt',header=None)\n",
    "zr=dfr[0]\n",
    "#zx=dfx[0]\n",
    "\n",
    "#zmax=20\n",
    "\n",
    "#l=int((zmax-0)/binstep)\n",
    "\n",
    "#bins=[]\n",
    "#for i in range(0,l):\n",
    "#    bins.append(round(binstep*i,1))\n",
    "    #bins.append(round(binstep*i,2))\n",
    "\n",
    "#hist_x_density=plt.hist(zx,bins=bins,density=True,color='blue',linewidth=5,histtype='step',alpha=0)\n",
    "#hist_x_density=np.histogram(zx,bins=bins,density=True)\n",
    "\n",
    "#hist_x_density=plt.hist(zx,bins=bins,density=True,alpha=0)\n",
    "#hist_x_number=np.histogram(zx,bins=bins)\n",
    "# distribuzione delle simulazioni (400)\n",
    "#Fx=hist_x_density[0]\n",
    "#zbinx=hist_x_density[1]\n",
    "#Nx=hist_x_number[0]\n",
    "#Fx=Nx/sum(Nx)     # Fx con stepbin =1 riproduce le percentuali del MOS quindi ok\n",
    "\n",
    "# con plt.hist credo si deve moltiplicare per stebin ogni bin (con 1 viene automatico ma con 0.3 no)\n",
    "\n",
    "# percentuale di casi per ciascun zbin fino a z=3\n",
    "\n",
    "#PGRBinz_x=Fx[0:10]\n",
    "\n",
    "hist_r_density=np.histogram(zr,bins=bins,density=True)\n",
    "hist_r_number=np.histogram(zr,bins=bins)\n",
    "# distribuzione delle simulazioni (400)\n",
    "# percentuale di casi per ciascun zbin\n",
    "#Fr=hist_r_density[0]\n",
    "zbin_r=hist_r_density[1]\n",
    "Nr=hist_r_number[0]\n",
    "Fr=Nr/sum(Nr)\n",
    "\n",
    "frac=len(zr)/len(zx)\n",
    "\n",
    "\n",
    "# efficienza di detection in bin di ampiezza 0.02\n",
    "#Pbns_ET2CE=np.repeat(Pbns_ET2CE_03,15)\n",
    "#Pbns_ETCE=np.repeat(Pbns_ETCE_03,15)\n",
    "#Pbns_ET=np.repeat(Pbns_ET_03,15)\n",
    "\n",
    "# efficienza di detection in bin di ampiezza 0.1\n",
    "Pbns_ET2CE5deg=np.repeat(Pbns_ET2CE5deg_03,3)\n",
    "Pbns_ETCE5deg=np.repeat(Pbns_ETCE5deg_03,3)\n",
    "Pbns_ET5deg=np.repeat(Pbns_ET5deg_03,3)\n",
    "\n",
    "\n",
    "#x1=x1[0:6]\n",
    "#x2=x2[0:6]\n",
    "\n",
    "# percentuale vista anche con GW\n",
    "#P_TH_ETCE=PGRBinz_x*Pbns_ETCE\n",
    "#P_TH_ET=PGRBinz_x*Pbns_ET\n",
    "P_TH_ET2CE5deg=Fx0[0:leff]*Pbns_ET2CE5deg\n",
    "P_TH_ETCE5deg=Fx0[0:leff]*Pbns_ETCE5deg\n",
    "P_TH_ET5deg=Fx0[0:leff]*Pbns_ET5deg\n",
    "#P_TH_ETCE=Fx[0:150]*Pbns_ETCE\n",
    "#P_TH_ET=Fx[0:150]*Pbns_ET\n",
    "#P_TH_ETCE=Pbns_ETCE\n",
    "#P_TH_ET=Pbns_ET\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "# plot histogram\n",
    "#plt.bar((x1+x2)/2,Pbns_ETCE,align='center',width=(x2-x1),label='THESEUS+ET+CE') # A bar chart\n",
    "#plt.hist(zr,bins=bins,density=True,color='red',linewidth=8,histtype='step',alpha=0.7,label='THESEUS short GRB with IRT detection')\n",
    "\n",
    "#plt.bar(zbinx[0:10],100*PGRBinz_x,width=binstep,align='edge',log=False,linewidth=4,color=\"None\",edgecolor='blue',alpha=0.9,label='THESEUS short GRB')\n",
    "#plt.bar(zbin_r[0:10],100*PGRBinz_xr,width=binstep,align='edge',log=False,linewidth=4,color=\"None\",edgecolor='red',alpha=0.9,label='THESEUS short GRB with IRT')\n",
    "plt.bar(zbinx[0:-1],100*Fx0,width=binstep,align='edge',log=False,linewidth=2,color=\"None\",edgecolor='cornflowerblue',alpha=0.9,label='THESEUS short GRB')\n",
    "\n",
    "\n",
    "#plt.bar((x1+x2)/2.,100*P_TH_ETCE,width=(x2-x1),log=False,align='center',hatch='\\\\',color='magenta',alpha=0.5,label='THESEUS+ET+CE short GRB')\n",
    "#plt.bar((x1+x2)/2.,100*P_TH_ET,width=(x2-x1),log=False,align='center',hatch='/',color='lime',alpha=0.5,label='THESEUS+ET short GRB')\n",
    "#plt.bar(zbinx[0:10],100*P_TH_ETCE,width=(x2-x1),log=False,align='edge',hatch='\\\\',color='magenta',alpha=0.5,label='THESEUS+ET+CE short GRB')\n",
    "#plt.bar(zbinx[0:10],100*P_TH_ET,width=(x2-x1),log=False,align='edge',hatch='/',color='lime',alpha=0.5,label='THESEUS+ET short GRB')\n",
    "plt.bar(zbinx[0:leff],100*P_TH_ET2CE5deg,width=binstep,log=False,align='edge',hatch='-',color='violet',alpha=0.5,label='THESEUS+ET+2CE short GRB')\n",
    "plt.bar(zbinx[0:leff],100*P_TH_ETCE5deg,width=binstep,log=False,align='edge',hatch='\\\\',color='magenta',alpha=0.5,label='THESEUS+ET+CE short GRB')\n",
    "plt.bar(zbinx[0:leff],100*P_TH_ET5deg,width=binstep,log=False,align='edge',hatch='/',color='lime',alpha=0.5,label='THESEUS+ET short GRB')\n",
    "#plt.bar(zbinx[0:150],100*P_TH_ETCE,width=binstep,log=False,align='edge',hatch='\\\\',color='magenta',alpha=0.5,label='THESEUS+ET+CE short GRB')\n",
    "#plt.bar(zbinx[0:150],100*P_TH_ET,width=binstep,log=False,align='edge',hatch='/',color='lime',alpha=0.5,label='THESEUS+ET short GRB')\n",
    "plt.bar(zbin_r[0:-1],100*Fr*frac,width=binstep,align='edge',log=False,linewidth=2,color=\"None\",edgecolor='indigo',alpha=0.9,label='THESEUS short GRB with IRT')\n",
    "\n",
    "#plt.xlim(0.,13)\n",
    "plt.xlim(0.,6)\n",
    "plt.ylim(0,7.5)\n",
    "plt.title('Aligned short GRBs')\n",
    "plt.xlabel('redshift')\n",
    "plt.ylabel(' short GRB [$\\%$ in each redshift bin]')\n",
    "\n",
    "#plt.legend(fontsize='small')\n",
    "plt.legend()\n",
    "\n",
    "#plt.savefig('SGRB_hist_ET_CE_onaxis_bin01.pdf')\n",
    "plt.savefig('SGRB_hist_ET_CE_onaxis_bin01.png',dpi=600, bbox_inches = \"tight\")\n",
    "#plt.savefig('SGRB_hist_ET_2CE_onaxis_bin01.tif',dpi=600, bbox_inches = \"tight\")\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print('--- ON-AXIS SHORT GRBs ---')\n",
    "\n",
    "onaxisGRBperyear=Fx0*12\n",
    "\n",
    "print('')\n",
    "print('THESEUS+ET:',100*sum(P_TH_ET5deg),'%')\n",
    "print('THESEUS+ET+CE:',100*sum(P_TH_ETCE5deg),'%')\n",
    "print('THESEUS+ET+2CE:',100*sum(P_TH_ET2CE5deg),'%')\n",
    "print('')\n",
    "\n",
    "print('Total GRB THESEUS (on axis) in 1(3.45)yr:', 12,'(',12*3.45,')')\n",
    "print('Total GRB THESEUS+ET (on axis) in 1(3.45)yr:', 12*sum(P_TH_ET5deg),'(',12*3.45*sum(P_TH_ET5deg),')')\n",
    "print('Total GRB THESEUS+ET+CE (on axis) in 1(3.45)yr:', 12*sum(P_TH_ETCE5deg),'(',12*3.45*sum(P_TH_ETCE5deg),')')\n",
    "print('Total GRB THESEUS+ET+2CE (on axis) in 1(3.45)yr:', 12*sum(P_TH_ET2CE5deg),'(',12*3.45*sum(P_TH_ET2CE5deg),')')\n",
    "\n",
    "# only if bin <=0.1\n",
    "print('')\n",
    "#print('Total GRB THESEUS+2G (at z<0.1) in 3.45 yr: ',sum(40*Fx[0]))\n",
    "#print('Total GRB THESEUS+2G (at z<0.08) in 3.45 yr: ',sum(40*Fx[0]))\n",
    "#print('Total GRB THESEUS+2G (at z<0.06) in 3.45 yr: ',sum(40*Fx[0]))\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}