Commit 0d24d483 authored by Ignacio Magana's avatar Ignacio Magana 💬

Merge branch 'master' into 'master'

Minor Pdet changes

See merge request lscsoft/gwcosmo!33
parents a975a92f 236d9de7
......@@ -18,17 +18,13 @@ import time
import progressbar
import pkg_resources
import os
"""
We want to create a function for $p(D|z,H_{0},I)$, so that when it is
passed a value of $z$ and $H_0$, a probability of detection is returned.
This involves marginalising over masses, inclination,
polarisation, and sky location.
"""
class DetectionProbability(object):
"""
Class to compute p(det | d_L, detectors, m1, m2, ...)
Class to compute the detection probability p(D |z, H0) as a function of z and H0
by marginalising over masses, inclination, polarisation, and sky location for a
set of detectors at a given sensitivity.
Parameters
----------
......@@ -62,7 +58,8 @@ class DetectionProbability(object):
"""
def __init__(self, mass_distribution, asd, detectors=['H1', 'L1'],
Nsamps=5000, network_snr_threshold=12, Omega_m=0.308,
linear=False, basic=False, alpha=1.6, M1=50., M2=50., constant_H=False, inspiral=True):
linear=False, basic=False, alpha=1.6, M1=50., M2=50.,
constant_H=False, inspiral=True):
self.data_path = pkg_resources.resource_filename('gwcosmo', 'data/')
self.mass_distribution = mass_distribution
self.asd = asd
......@@ -95,7 +92,7 @@ class DetectionProbability(object):
# TODO: For higher values of z (z=10) this goes
# outside the range of the psds and gives an error
self.z_array = np.logspace(-4.0, 0.5, 500)
# set up the samples for monte carlo integral
N = self.Nsamps
self.RAs = np.random.rand(N)*2.0*np.pi
......@@ -147,13 +144,12 @@ class DetectionProbability(object):
#to solve this for interpolation purposes, set logit(prob=1)=100, so then expit(100)=logit^-1(100)=1
#instead of 100 anything from 35 to sys.float_info.max can be set as in this range expit is effectively 1
#yet: higher values make interpolation less effective
logit_prob=logit(self.prob)
for i in range (len(logit_prob)):
logit_prob[i]=np.where(logit_prob[i]==float('+inf'), 100, logit_prob[i])
self.interp_average = interp2d(self.z_array, self.H0vec, logit_prob, kind='linear')
def mchirp(self, m1, m2):
"""
Calculates the source chirp mass
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment