update and merge README

ComputePosteriors.sh

@@ -0,0 +1,21 @@
+#!/bin/bash
+
+
+echo "    Recomputing Gravitational Wave posteriors on f_PBH (may take ~ 1 hour...)"
+echo " "
+
+for NOBS in 1 80
+do
+    for PRIOR in J LF
+	do
+	    python src/tabulate_posteriors_O3.py -N_obs $NOBS -prior $PRIOR -M_PBH 0.5
+	done
+done
+
+for NOBS in 1 24000
+do
+    for PRIOR in J LF
+	do
+	    python src/tabulate_posteriors_ET.py -N_obs $NOBS -prior $PRIOR
+	done
+done

README.md

@@ -18,11 +18,10 @@ The figures can be remade with the following commands:
 
     `python plot_frequentist_limits.py -plot_ps_diff`
 
-Computing the point-source limits requires making tables containing `p_gamma`, the probability that an individual PBH is detectable by Fermi-LAT. The tables used for the analysis in the paper are contained in the directory `data/p_gammas/`. These may take a few hours to recompute. They can be regenerated with
+Computing the point-source limits requires making tables containing `p_gamma`, the probability that an individual PBH is detectable by Fermi-LAT. The tables used for the analysis in the paper are contained in the directory [`data/p_gammas/`](data/p_gammas/). These may take a few hours to recompute. They can be regenerated with `python generate_p_gamma_tables.py -test`. With the `-test` flag, the `p_gamma` tables will be written to [`data/p_gammas/test/`](data/p_gammas/test/) instead of overwriting the precomputed tables. The script can also be used to generate `p_gamma` tables over different `(m_dm, sv)` grids.
 
-    `python generate_p_gamma_tables.py -test`
+All limit calculations require tabulated posteriors of `f_pbh`. These are provided in [`data/posteriors_f`](data/posteriors_f/) but they can be recomputed by running the script [`ComputePosteriors.sh`](ComputePosteriors.sh). This script calls a number of scripts in [`src/`](src/) which calculate the posteriors from gravitational wave observations. Run `python src/tabulate_posteriors_O3.py -h` and `python src/tabulate_posteriors_ET.py -h` for more detailed usage of these scripts. Scripts for calculating posteriors from SKA will be added shortly.
 
-With the `-test` flag, the `p_gamma` tables will be written to `data/p_gammas/test/` instead of overwriting the precomputed tables. The script can also be used to generate `p_gamma` tables over different `(m_dm, sv)` grids. 
 
 ## `bayesian` branch
 

data/ET_redshift.txt

@@ -0,0 +1,5 @@
+# Total source frame mass 20 M_sun
+# Columns: Redshift, Fraction detectable
+7.678380919212519	0.5
+28.52205602981411	0.1
+95.87373195587614	0

data/SubSolar_Horizon.txt

@@ -0,0 +1,11 @@
+# Figure 1 (dashed green line) of https://arxiv.org/pdf/1808.04771.pdf
+# Columns: Component Mass [M_sun],  Horizon [Mpc]
+0.1981205951448708 27.286681430014227
+0.3008613938919342 38.37156747012553
+0.4004698512137823 48.75339548973085
+0.5000783085356305 58.784038259116684
+0.6003132341425216 68.46363328229641
+0.7005481597494128 77.70424674270178
+0.800156617071261 86.8569263865389
+0.9003915426781521 95.83415090927964
+1.0025058731401724 104.8117879440606

data/posteriors_f/Posterior_f_ET_Prior_J_M=10.0_N=1.txt

@@ -1,4 +1,4 @@
-# Posterior distribution for f, given N = 1 merger events in ET. PBH Mass: 10.0 M_sunPrior: jeffreys
+# Posterior distribution for f, given N = 1 merger events in ET. PBH Mass: 10.0 M_sunPrior: J
 # Columns: f, P(f|N)
 9.999999999999999547e-07 7.246589174428889635e+01
 1.071891319205128528e-06 8.324571817517015404e+01

data/posteriors_f/Posterior_f_O3_Prior_J_M=0.5_N=1.txt

@@ -1,4 +1,4 @@
-# Posterior distribution for f, given N = 1 merger events at LIGO O3. PBH Mass: 0.5 M_sunPrior: jeffreys
+# Posterior distribution for f, given N = 1 merger events at LIGO O3. PBH Mass: 0.5 M_sunPrior: J
 # Columns: f, P(f|N)
 9.999999999999999547e-07 8.416139272344972731e-06
 1.071891319205128528e-06 9.669730117647979165e-06

requirements.txt

@@ -1,3 +1,5 @@
 matplotlib==3.0.3
 scipy==1.2.1
 numpy==1.16.2
+tqdm==4.19.9
+

src/Cosmo.py

@@ -0,0 +1,75 @@
+""" Cosmo.py
+
+Functions/constants for calculating cosmological quantities
+such as the Hubble rate and various distances.
+
+"""
+
+import numpy as np
+
+from scipy.integrate import cumtrapz, quad, dblquad
+from scipy.interpolate import interp1d
+
+#--- Some constants ------------
+#-------------------------------
+
+G_N = 4.302e-3 #(pc/solar mass) (km/s)^2
+G_N_Mpc = 1e-6*4.302e-3 #(Mpc/solar mass) (km/s)^2
+
+h = 0.678
+Omega_DM = 0.1186/(h**2)
+H0 = 100.0*h #(km/s) Mpc^-1
+H0_peryr = 67.8*(3.24e-20)*(60*60*24*365)
+ageUniverse = 13.799e9 #y
+Omega_L = 0.692
+Omega_m = 0.308
+Omega_r = 9.3e-5
+
+D_H = (3000/h) #Mpc
+
+c = 3e5 #km/s
+
+
+#-------------------------------
+
+#--- Useful cosmological functions ---
+#-------------------------------------
+
+
+rho_critical = 3.0*H0**2/(8.0*np.pi*G_N_Mpc) #Solar masses per Mpc^3
+
+def Hubble(z):
+    return H0_peryr*np.sqrt(Omega_L + Omega_m*(1+z)**3 + Omega_r*(1+z)**4)
+
+def Hubble2(z):
+    return H0*np.sqrt(Omega_L + Omega_m*(1+z)**3 + Omega_r*(1+z)**4)
+
+def HubbleLaw(age):
+    return H0_peryr*age
+
+def rho_z(z):
+    return 3.0*Hubble2(z)**2/(8*np.pi*G_N)
+
+def t_univ(z):
+    integ = lambda x: 1.0/((1+x)*Hubble(x))
+    return quad(integ, z, np.inf)[0]
+
+def Omega_PBH(f):  
+    return f*Omega_DM
+
+#Luminosity distance (Mpc)
+def calcdL(z):
+    c = 3.06594845e-7
+    return c*(1+z)*quad(lambda x: Hubble(x)**-1, 0, z)[0]
+
+#https://arxiv.org/pdf/astro-ph/9905116.pdf
+#Comoving distance (in a flat universe, in Mpc)
+def D_C(z):
+    return c*quad(lambda z: 1/Hubble2(z), 0, z)[0]
+
+#Comoving volume out to redshift z (in Mpc^3)
+def V_C(z):
+    return (4*np.pi/3)*(D_C(z))**3
+
+def dVdz(z):
+    return 4*np.pi*D_C(z)**2*(c/Hubble2(z))