Commit fd8c2e98 authored by Riccardo Barbieri's avatar Riccardo Barbieri
Browse files

Merge branch 'schech_params_commandline' into 'pending_review'

Schech params commandline

See merge request cbc-cosmo/gwcosmo!21
parents 0c204878 87b79b6e
......@@ -54,6 +54,14 @@ parser = OptionParser(
help="Hard redshift cut to apply to the galaxy catalogue (default=None)"),
Option("--mth", default=None,
help="Override the apparent magnitude threshold of the catalogue, if provided (default=None)"),
Option("--schech_alpha", default=None,
help="Override the default value for slope of schechter function for given band, if provided (default=None)"),
Option("--schech_Mstar", default=None,
help="Override the default value for Mstar of schechter function for given band, if provided (default=None)"),
Option("--schech_Mmin", default=None,
help="Override the default value for Mmin of schechter function for given band, if provided (default=None)"),
Option("--schech_Mmax", default=None,
help="Override the default value for Mmax of schechter function for given band, if provided (default=None)"),
Option("--nside", default='32', type=int,
help="skymap nside choice for reading in galaxies from the overlap of catalogue and skymap (default=32)"),
Option("--sky_area", default='0.999', type=float,
......@@ -89,6 +97,26 @@ if opts.mth is None:
else:
mth_str = f"--mth {opts.mth}"
if opts.schech_alpha is None:
schech_alpha_str = ""
else:
schech_alpha_str = f"--schech_alpha {opts.schech_alpha}"
if opts.schech_Mstar is None:
schech_Mstar_str = ""
else:
schech_Mstar_str = f"--schech_Mstar {opts.schech_Mstar}"
if opts.schech_Mmin is None:
schech_Mmin_str = ""
else:
schech_Mmin_str = f"--schech_Mmin {opts.schech_Mmin}"
if opts.schech_Mmax is None:
schech_Mmax_str = ""
else:
schech_Mmax_str = f"--schech_Mmax {opts.schech_Mmax}"
if opts.zcut is None:
zcut_str = ""
else:
......@@ -149,6 +177,10 @@ common_args = f"\
--zmax {opts.zmax} \
{zcut_str} \
{mth_str} \
{schech_alpha_str} \
{schech_Mstar_str} \
{schech_Mmin_str} \
{schech_Mmax_str} \
--nside {opts.nside} \
--sky_area {opts.sky_area} \
--min_pixels {opts.min_pixels} \
......
#!/usr/bin/env python3
"""
This script computes posterior on H0 using a single gravitational wave event
and an electromagnetic counterpart.
......@@ -106,6 +107,14 @@ parser = OptionParser(
help="Hard redshift cut to apply to the galaxy catalogue (default=%default)"),
Option("--mth", default=None,
help="Override the apparent magnitude threshold of the catalogue, if provided (default=None)"),
Option("--schech_alpha", default=None,
help="Override the default value for slope of schechter function for given band, if provided (default=None)"),
Option("--schech_Mstar", default=None,
help="Override the default value for Mstar of schechter function for given band, if provided (default=None)"),
Option("--schech_Mmin", default=None,
help="Override the default value for Mmin of schechter function for given band, if provided (default=None)"),
Option("--schech_Mmax", default=None,
help="Override the default value for Mmax of schechter function for given band, if provided (default=None)"),
Option("--nside", default=32, type=int,
help="skymap nside choice for reading in galaxies from the overlap of catalogue and skymap (default=32)"),
Option("--sky_area", default=0.999, type=float,
......@@ -304,8 +313,25 @@ else:
zcut = float(opts.zcut)
else:
zcut = None
if opts.schech_alpha is not None:
schech_alpha=float(opts.schech_alpha)
else:
schech_alpha = None
if opts.schech_Mstar is not None:
schech_Mstar=float(opts.schech_Mstar)
else:
schech_Mstar = None
if opts.schech_Mmin is not None:
schech_Mmin=float(opts.schech_Mmin)
else:
schech_Mmin = None
if opts.schech_Mmax is not None:
schech_Mmax=float(opts.schech_Mmax)
else:
schech_Mmax = None
sp = SchechterParams(band)
sp = SchechterParams(band, schech_alpha, schech_Mstar, schech_Mmin, schech_Mmax)
print("Schechter function with parameters: alpha={}, Mstar={}, Mmin={}, Mmax={}, ".format(sp.alpha, sp.Mstar, sp.Mmin, sp.Mmax))
p_M = SchechterMagFunction(Mstar_obs=sp.Mstar,alpha=sp.alpha)
if galaxy_weighting:
ps_M = gwcosmo.gwcosmo.LuminosityWeighting()
......@@ -314,7 +340,10 @@ else:
if opts.method == 'statistical':
me = gwcosmo.gwcosmo.WholeSkyGalaxyCatalogLikelihood(catalog, skymap, band, cosmo, px_zH0, pdet.pD_zH0_eval, zprior, ps_z, p_M, ps_M, Kcorr=Kcorr, mth=mth, zcut=zcut, zmax=zmax,zuncert=zuncert, complete_catalog=complete_catalog, sky_thresh = opts.sky_area, nside=opts.nside)
me = gwcosmo.gwcosmo.WholeSkyGalaxyCatalogLikelihood(catalog, skymap, band, sp, cosmo, px_zH0, pdet.pD_zH0_eval, zprior, ps_z, p_M, ps_M, Kcorr=Kcorr, mth=mth, zcut=zcut, zmax=zmax,zuncert=zuncert, complete_catalog=complete_catalog, sky_thresh = opts.sky_area, nside=opts.nside)
if opts.method == 'pixel':
......@@ -345,7 +374,7 @@ else:
print("Setting up p(x|z,H0)")
px_zH0 = pixelated_samples.make_los_px_function(samp_ind,H0,hyper_params_dict,name=mass_distribution,reweight_samples=reweight_samples)
me = gwcosmo.gwcosmo.SinglePixelGalaxyCatalogLikelihood(pixel_index, catalog, skymap, band, cosmo,\
me = gwcosmo.gwcosmo.SinglePixelGalaxyCatalogLikelihood(pixel_index, catalog, skymap, band, sp, cosmo,\
px_zH0, pdet.pD_zH0_eval, zprior, ps_z, p_M, \
ps_M, outputfile, Kcorr=Kcorr, mth=mth, zcut=zcut, \
zmax=zmax,zuncert=zuncert, complete_catalog=complete_catalog, \
......
......@@ -133,6 +133,8 @@ class GalaxyCatalogLikelihood(gwcosmoLikelihood):
----------
skymap : gwcosmo.likelihood.skymap.skymap object
provides p(x|Omega) and skymap properties
sp : gwcosmo.utilities.schechter_params.SchechterParams class
Class that stores the schechter function parameters alpha, Mstar, Mmin, Mmax
fast_cosmology : gwcosmo.utilities.standard_cosmology.fast_cosmology object
Cosmological model
Kcorr : bool, optional
......@@ -150,7 +152,7 @@ class GalaxyCatalogLikelihood(gwcosmoLikelihood):
"""
def __init__(self, skymap, observation_band, fast_cosmology, px_zH0, pD_zH0, zprior, zrates, luminosity_prior, luminosity_weights, Kcorr=False, zmax=10.):
def __init__(self, skymap, observation_band, sp, fast_cosmology, px_zH0, pD_zH0, zprior, zrates, luminosity_prior, luminosity_weights, Kcorr=False, zmax=10.):
"""
Parameters
----------
......@@ -187,7 +189,6 @@ class GalaxyCatalogLikelihood(gwcosmoLikelihood):
self.Kcorr = Kcorr
self.band = observation_band
sp = SchechterParams(self.band)
self.Mmin_obs = sp.Mmin
self.Mmax_obs = sp.Mmax
......@@ -478,6 +479,8 @@ class SinglePixelGalaxyCatalogLikelihood(GalaxyCatalogLikelihood):
The galaxy catalogue
skymap : gwcosmo.likelihood.skymap.skymap object
provides p(x|Omega) and skymap properties
sp : gwcosmo.utilities.schechter_params.SchechterParams class
Class that stores the schechter function parameters alpha, Mstar, Mmin, Mmax
fast_cosmology : gwcosmo.utilities.standard_cosmology.fast_cosmology object
Cosmological model
Kcorr : bool, optional
......@@ -497,7 +500,7 @@ class SinglePixelGalaxyCatalogLikelihood(GalaxyCatalogLikelihood):
"""
def __init__(self, pixel_index, galaxy_catalog, skymap, observation_band, fast_cosmology, px_zH0, pD_zH0, zprior, zrates, luminosity_prior, luminosity_weights, outputfile, Kcorr=False, mth=None, zcut=None, zmax=10.,zuncert=True, complete_catalog=False, nside=32, nside_low_res = None):
def __init__(self, pixel_index, galaxy_catalog, skymap, observation_band, sp, fast_cosmology, px_zH0, pD_zH0, zprior, zrates, luminosity_prior, luminosity_weights, outputfile, Kcorr=False, mth=None, zcut=None, zmax=10.,zuncert=True, complete_catalog=False, nside=32, nside_low_res = None):
"""
Parameters
----------
......@@ -540,7 +543,7 @@ class SinglePixelGalaxyCatalogLikelihood(GalaxyCatalogLikelihood):
The high-resolution value of nside to subdivide the current pixel into
"""
super().__init__(skymap, observation_band, fast_cosmology, px_zH0, pD_zH0, zprior, zrates, luminosity_prior, luminosity_weights, Kcorr=Kcorr, zmax=zmax)
super().__init__(skymap, observation_band, sp, fast_cosmology, px_zH0, pD_zH0, zprior, zrates, luminosity_prior, luminosity_weights, Kcorr=Kcorr, zmax=zmax)
self.nside_low_res = nside_low_res
self.zcut = zcut
self.complete_catalog = complete_catalog
......@@ -786,7 +789,7 @@ class WholeSkyGalaxyCatalogLikelihood(GalaxyCatalogLikelihood):
catalogue method.
"""
def __init__(self, galaxy_catalog, skymap, observation_band, fast_cosmology, px_zH0, pD_zH0, zprior, zrates, luminosity_prior, luminosity_weights, Kcorr=False, mth=None, zcut=None, zmax=10.,zuncert=True, complete_catalog=False, sky_thresh = 0.999, nside=32):
def __init__(self, galaxy_catalog, skymap, observation_band, sp, fast_cosmology, px_zH0, pD_zH0, zprior, zrates, luminosity_prior, luminosity_weights, Kcorr=False, mth=None, zcut=None, zmax=10.,zuncert=True, complete_catalog=False, sky_thresh = 0.999, nside=32):
"""
Parameters
----------
......@@ -796,6 +799,8 @@ class WholeSkyGalaxyCatalogLikelihood(GalaxyCatalogLikelihood):
The GW skymap
observation_band : str
Observation band (eg. 'B', 'K', 'u', 'g')
sp : gwcosmo.utilities.schechter_params.SchechterParams class
Class that stores the schechter function parameters alpha, Mstar, Mmin, Mmax
fast_cosmology : object
Fast cosmology
px_zH0 : object
......@@ -828,7 +833,7 @@ class WholeSkyGalaxyCatalogLikelihood(GalaxyCatalogLikelihood):
TODO: UPDATE this for new catalog classes
"""
super().__init__(skymap, observation_band, fast_cosmology, px_zH0, pD_zH0, zprior, zrates, luminosity_prior, luminosity_weights, Kcorr=Kcorr, zmax=zmax)
super().__init__(skymap, observation_band, sp, fast_cosmology, px_zH0, pD_zH0, zprior, zrates, luminosity_prior, luminosity_weights, Kcorr=Kcorr, zmax=zmax)
self.mth = mth
self.zcut = zcut
......
......@@ -16,7 +16,7 @@ class SchechterParams():
(note paper quotes M*=-23.55 but means M*=-23.55 + 5log10(h))
"""
def __init__(self, band):
def __init__(self, band, schech_alpha=None, schech_Mstar=None, schech_Mmin=None, schech_Mmax=None):
"""
Parameters
----------
......@@ -28,23 +28,32 @@ class SchechterParams():
self.Mmin = None
self.Mmax = None
self.alpha, self.Mstar, self.Mmin, self.Mmax = self.values(band)
self.alpha, self.Mstar, self.Mmin, self.Mmax = self.default_values(band) #initialize to default values
def values(self, band):
if schech_alpha is not None: #change to input value if provided
self.alpha=schech_alpha
if schech_Mstar is not None:
self.Mstar=schech_Mstar
if schech_Mmin is not None:
self.Mmin=schech_Mmin
if schech_Mmax is not None:
self.Mmax=schech_Mmax
def default_values(self, band):
if band == 'B':
return -1.21, -19.70, -22.96, -12.96
elif band == 'K':
return -1.02, -23.55, -27.0, -12.96
elif band == 'u':
return -0.92, -17.93, -21.93, -15.54 #TODO check Mmin and Mmax
elif band == 'u': #These values are actually u', g', r', i', zi, a.k.a. redshifted to the median redshift (0.1) of SDSS.
return -0.92, -17.93, -21.93, -15.54
elif band == 'g':
return -0.89, -19.39, -23.38, -16.10 #TODO check Mmin and Mmax
return -0.89, -19.39, -23.38, -16.10
elif band == 'r':
return -1.05, -20.44, -24.26, -16.11 #TODO check Mmin and Mmax
return -1.05, -20.44, -24.26, -16.11
elif band == 'i':
return -1.00, -20.82, -23.84, -17.07 #TODO check Mmin and Mmax
return -1.00, -20.82, -23.84, -17.07
elif band == 'z':
return -1.08, -21.18, -24.08, -17.34 #TODO check Mmin and Mmax
return -1.08, -21.18, -24.08, -17.34
elif band == 'W1':
return -1.12, -24.09, -28, -16.6 # https://iopscience.iop.org/article/10.1088/0004-637X/697/1/506/pdf Tab 3 (All 3.6)
# https://arxiv.org/pdf/1702.07829.pdf (Mmin Mmax Fig2)
......
Supports Markdown
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