Commit 3b1f8e98 authored by Ignacio Magana's avatar Ignacio Magana 💬

Merge branch 'pdets' into 'master'

Pdets

See merge request lscsoft/gwcosmo!20
parents b9a6e030 60c4baee
......@@ -58,16 +58,18 @@ parser = OptionParser(
help="Set number of H0 Posterior bins"),
Option("--posterior_samples", default=None,
help="LALinference posterior samples file in format (.dat or hdf5) or use GW170817, GW170814, GW170818"),
Option("--pdet", default=None,
help="Provide a precomputed Pdet.p selection function."),
Option("--mass_distribution", default=None,
help="Choose between BNS or BBH mass distributions for Pdet calculation."),
help="Choose between BNS-gaussian or BBH-powerlaw mass distributions for default Pdet calculations."),
Option("--psd", default=None,
help="Select between 'O1' and 'O2' PSDs, for default Pdet calculations. By default we use aLIGO at design sensitivity."),
Option("--powerlaw_slope", default='1.6', type=float,
help="Set powerlaw slope for BBH powelaw mass distribution for Pdet calculation."),
Option("--skymap", default=None,
help="LALinference 3D skymap file in format (.fits)"),
Option("--galaxy_catalog", default=None,
help="Load galaxy catalog in pickle format"),
Option("--psd", default=None,
help="Select between 'O1' and 'O2' PSDs, by default we use aLIGO at design sensitivity (default=None)."),
Option("--galaxy_weighting", default='False',
help="Set galaxy catalog weighting"),
Option("--completion", default='False',
......@@ -123,26 +125,26 @@ if opts.plot is not None:
outputfile = filename[:-4]
else:
if (opts.posterior_samples is None and
opts.skymap is None):
parser.error('Provide either posterior samples or skymap.')
if (opts.posterior_samples is None and opts.skymap is None):
parser.error('Provide either posterior samples or skymap.')
if opts.mass_distribution is None:
parser.error('Provide a mass distribution to use for Pdet calculation.')
if (opts.mass_distribution is not None and opts.psd is not None and opts.pdet is None):
mass_distribution = str(opts.mass_distribution)
psd = str(opts.psd)
elif (opts.mass_distribution is None and opts.psd is None and opts.pdet is not None):
mass_distribution = None
psd = None
else:
parser.error('Provide a mass distribution or precomputed Pdet to use for calculation.')
if (opts.galaxy_catalog is None and opts.method == 'statistical'):
parser.error('The statistical method requires a galaxy catalog. Please provide one.')
if opts.psd is None:
parser.error('Please provide a PSD.')
if opts.posterior_samples is not None:
posterior_samples = str(opts.posterior_samples)
if opts.skymap is not None:
skymap = str(opts.skymap)
posterior_samples = str(opts.posterior_samples)
if opts.mass_distribution is not None:
mass_distribution = str(opts.mass_distribution)
if opts.skymap is not None:
skymap = str(opts.skymap)
if opts.galaxy_catalog is not None:
galaxy_catalog = str(opts.galaxy_catalog)
......@@ -175,7 +177,6 @@ else:
linear = str2bool(opts.linear_cosmology)
basic = str2bool(opts.basic_pdet)
Lambda = float(opts.Lambda)
psd = str(opts.psd)
alpha = float(opts.powerlaw_slope)
options_string = opts.method
......@@ -239,20 +240,24 @@ else:
catalog = None
allsky = True
radeclims = None
if mass_distribution == 'BBH-powerlaw':
pdet_path = data_path + '{}PSD_{}_alpha_{}_5000Nsamps_z_H0_pD_array.p'.format(psd, mass_distribution, alpha)
if os.path.isfile(pdet_path) is True:
pdet = pickle.load(open(pdet_path, 'rb'))
print('Loading precomputed pdet with a {} mass distribution (alpha = {}) at {} sensitivity'.format(mass_distribution, alpha, psd))
else:
pdet = gwcosmo.detection_probability.DetectionProbability(mass_distribution=mass_distribution, psd=psd, basic=basic)
if (mass_distribution == 'BNS-gaussian' or mass_distribution == 'BNS-uniform'):
if mass_distribution == 'BNS-gaussian':
pdet_path = data_path + '{}PSD_{}_5000Nsamps_z_H0_pD_array.p'.format(psd, mass_distribution)
if os.path.isfile(pdet_path) is True:
pdet = pickle.load(open(pdet_path, 'rb'))
print('Loading precomputed pdet with a {} mass distribution at {} sensitivity'.format(mass_distribution, psd))
if mass_distribution is None:
pdet_path = str(opts.pdet)
print(pdet_path)
if os.path.isfile(pdet_path) is True:
pdet = pickle.load(open(pdet_path, 'rb'))
print('Loading provided pdet.')
me = gwcosmo.gwcosmo.gwcosmoLikelihood(samples,skymap,catalog,pdet,EM_counterpart=counterpart,
linear=linear,weighted=galaxy_weighting,whole_cat=allsky,
radec_lim=radeclims,uncertainty=uncertainty,basic=basic,
......
......@@ -93,7 +93,7 @@ class gwcosmoLikelihood(object):
basic=False, uncertainty=False, Lambda=0, area=0.999):
self.pdet = pdet
self.event_type = pdet.mass_distribution
self.psd = pdet.psd
self.psd = pdet.asd
self.Omega_m = Omega_m
self.linear = linear
self.weighted = weighted
......
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