Commit d8c66747 authored by Charlie Hoy's avatar Charlie Hoy
Browse files

Merge branch 'add_sky_sensitivity' into 'master'

parents fdebc6f4 83e66a50
# Copyright (C) 2018 Charlie Hoy <charlie.hoy@ligo.org>
# This program is free software; you can redistribute it and/or modify it
# under the terms of the GNU General Public License as published by the
# Free Software Foundation; either version 3 of the License, or (at your
# option) any later version.
#
# This program is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
# Public License for more details.
#
# You should have received a copy of the GNU General Public License along
# with this program; if not, write to the Free Software Foundation, Inc.,
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
__version__ = "beta"
......@@ -65,6 +65,9 @@ def command_line():
parser.add_argument("-c", "--config", dest="config",
help="configuration file associcated with each samples file.", nargs='+',
default=None)
parser.add_argument("--sensitivity", action="store_true",
help="generate sky sensitivities for HL, HLV",
default=False)
return parser
def run_checks(opts):
......@@ -218,6 +221,16 @@ def make_plots(opts, colors=None):
fig = plot._1d_histogram_plot(j, param_samples, latex_labels[j])
plt.savefig("%s/plots/1d_posterior_%s_%s.png" %(opts.webdir, i, j))
plt.close()
if opts.sensitivity:
fig = plot._sky_sensitivity(["H1", "L1"], 0.2,
combined_maxL[num])
plt.savefig("%s/plots/%s_sky_sensitivity_HL" %(opts.webdir, i))
plt.close()
fig = plot._sky_sensitivity(["H1", "L1", "V1"], 0.2,
combined_maxL[num])
plt.savefig("%s/plots/%s_sky_sensitivity_HLV" %(opts.webdir, i))
plt.close()
# if len(approximants) > 1, then we need to do comparison plots
if len(opts.approximant) > 1:
for ind, j in enumerate(parameters):
......@@ -376,6 +389,20 @@ def make_comparison_pages(opts, approximants, samples, colors, parameters):
"{}/plots/compare_waveforms.png".format(opts.baseurl)]]
html_file.make_table_of_images(headings=["sky_map", "waveform"],
contents=contents)
if opts.sensitivity:
html_file.add_content("<div class='row justify-content-center' "
"style='margin=top: 2.0em;'>"
"<p>To see the sky sensitivity for the following "
"networks, click the button</p></div>")
html_file.add_content("<div class='row justify-content-center' "
"style='margin-top: 0.2em;'>"
"<button type='button' class='btn btn-light' "
"onclick='%s.src=\"%s/plots/combined_skymap.png\"'"
"style='margin-left:0.25em; margin-right:0.25em'>Sky Map</button>"
"<button type='button' class='btn btn-light' "
"onclick='%s.src=\"%s/plots/IMRPhenomPv2_sky_sensitivity_HL.png\"'"
"style='margin-left:0.25em; margin-right:0.25em'>HL</button>"
%("combined_skymap", opts.baseurl, "combined_skymap", opts.baseurl))
html_file.make_footer(user=os.environ["USER"], rundir=opts.webdir)
# edit all of the comparison pages
pages = ["Comparison_{}".format(i) for i in parameters]
......
......@@ -153,8 +153,8 @@ def _waveform_plot(maxL_params, **kwargs):
"""
logging.info("Generating the maximum likelihood waveform plot for H1")
if not LALSIMULATION:
raise Exception("LALSimulation could not be imported. Please install "
"LALSuite to be able to use all features")
raise exception("lalsimulation could not be imported. please install "
"lalsuite to be able to use all features")
delta_frequency = kwargs.get("delta_f", 1./256)
minimum_frequency = kwargs.get("f_min", 10.)
maximum_frequency = kwargs.get("f_max", 1000.)
......@@ -264,6 +264,9 @@ def _sky_map_plot(ra, dec, **kwargs):
ax = plt.subplot(111, projection="hammer")
ax.cla()
ax.grid()
ax.set_xticklabels([r"$22^{h}$", r"$20^{h}$", r"$18^{h}$", r"$16^{h}$",
r"$14^{h}$", r"$12^{h}$", r"$10^{h}$", r"$8^{h}$",
r"$6^{h}$", r"$4^{h}$", r"$2^{h}$"])
levels = [1.0 - np.exp(-0.5), 1 - np.exp(-2), 1-np.exp(-9./2.)]
H,X,Y = np.histogram2d(ra, dec, bins=50)
......@@ -328,6 +331,9 @@ def _sky_map_comparison_plot(ra_list, dec_list, approximants, colors, **kwargs):
ax = plt.subplot(111, projection="hammer")
ax.cla()
ax.grid()
ax.set_xticklabels([r"$22^{h}$", r"$20^{h}$", r"$18^{h}$", r"$16^{h}$",
r"$14^{h}$", r"$12^{h}$", r"$10^{h}$", r"$8^{h}$",
r"$6^{h}$", r"$4^{h}$", r"$2^{h}$"])
levels = [1.0 - np.exp(-0.5), 1 - np.exp(-2), 1-np.exp(-9./2.)]
for num, i in enumerate(ra_list):
......@@ -415,3 +421,117 @@ def _make_corner_plot(opts, samples, params, approximant, latex_labels,
height *= figure.dpi
np.savetxt("{}/plots/corner/{}_axes.txt".format(opts.webdir, approximant), extent)
return figure
def __get_cutoff_indices(flow, fhigh, df, N):
"""
Gets the indices of a frequency series at which to stop an overlap
calculation.
Parameters
----------
flow: float
The frequency (in Hz) of the lower index.
fhigh: float
The frequency (in Hz) of the upper index.
df: float
The frequency step (in Hz) of the frequency series.
N: int
The number of points in the **time** series. Can be odd
or even.
Returns
-------
kmin: int
kmax: int
"""
if flow:
kmin = int(flow / df)
else:
kmin = 1
if fhigh:
kmax = int(fhigh / df )
else:
kmax = int((N + 1)/2.)
return kmin,kmax
def _sky_sensitivity(network, resolution, maxL_params, **kwargs):
"""Generate the sky sensitivity for a given network
Parameters
----------
network: list
list of detectors you want included in your sky sensitivity plot
resolution: float
resolution of the skymap
maxL_params: dict
dictionary of waveform parameters for the maximum likelihood waveform
"""
logging.info("Generating the sky sensitivity for %s" %(network))
if not LALSIMULATION:
raise Exception("LALSimulation could not be imported. Please install "
"LALSuite to be able to use all features")
delta_frequency = kwargs.get("delta_f", 1./256)
minimum_frequency = kwargs.get("f_min", 20.)
maximum_frequency = kwargs.get("f_max", 1000.)
frequency_array = np.arange(minimum_frequency, maximum_frequency,
delta_frequency)
approx = lalsim.GetApproximantFromString(maxL_params["approximant"])
mass_1 = maxL_params["mass_1"]*MSUN_SI
mass_2 = maxL_params["mass_2"]*MSUN_SI
luminosity_distance = maxL_params["luminosity_distance"]*PC_SI*10**6
iota, S1x, S1y, S1z, S2x, S2y, S2z = \
lalsim.SimInspiralTransformPrecessingNewInitialConditions(
maxL_params["iota"], maxL_params["phi_jl"], maxL_params["tilt_1"],
maxL_params["tilt_2"], maxL_params["phi_12"], maxL_params["a_1"],
maxL_params["a_2"], mass_1, mass_2, kwargs.get("f_ref", 10.),
maxL_params["phase"])
h_plus, h_cross = lalsim.SimInspiralChooseFDWaveform(mass_1, mass_2, S1x,
S1y, S1z, S2x, S2y, S2z,
luminosity_distance, iota,
maxL_params["phase"], 0.0, 0.0, 0.0, delta_frequency,
minimum_frequency, maximum_frequency,
kwargs.get("f_ref", 10.), None, approx)
h_plus = h_plus.data.data
h_cross = h_cross.data.data
h_plus = h_plus[:len(frequency_array)]
h_cross = h_cross[:len(frequency_array)]
psd = {}
psd["H1"] = psd["L1"] = np.array([lalsim.SimNoisePSDaLIGOZeroDetHighPower(i) for i in frequency_array])
psd["V1"] = np.array([lalsim.SimNoisePSDVirgo(i) for i in frequency_array])
kmin, kmax = __get_cutoff_indices(minimum_frequency, maximum_frequency,
delta_frequency, (len(h_plus)-1)*2)
ra = np.arange(-np.pi, np.pi, resolution)
dec = np.arange(-np.pi, np.pi, resolution)
X,Y = np.meshgrid(ra, dec)
N = np.zeros([len(dec), len(ra)])
indices = np.ndindex(len(ra), len(dec))
for ind in indices:
ar = {}
SNR = {}
for i in network:
ard = __antenna_response(i, ra[ind[0]], dec[ind[1]],
maxL_params["psi"], maxL_params["geocent_time"])
ar[i] = [ard[0], ard[1]]
strain = np.array(h_plus*ar[i][0]+h_cross*ar[i][1])
integrand = np.conj(strain[kmin:kmax])*strain[kmin:kmax]/psd[i][kmin:kmax]
integrand = integrand[:-1]
SNR[i] = np.sqrt(4*delta_frequency*np.sum(integrand).real)
ar[i][0] *= SNR[i]
ar[i][1] *= SNR[i]
numerator = 0.0
denominator = 0.0
for i in network:
numerator += sum(i**2 for i in ar[i])
denominator += SNR[i]**2
N[ind[1]][ind[0]] = (((numerator/denominator)**0.5))
fig = plt.figure()
ax = plt.subplot(111, projection="hammer")
ax.cla()
ax.grid()
plt.pcolormesh(X,Y,N)
ax.set_xticklabels([r"$22^{h}$", r"$20^{h}$", r"$18^{h}$", r"$16^{h}$",
r"$14^{h}$", r"$12^{h}$", r"$10^{h}$", r"$8^{h}$",
r"$6^{h}$", r"$4^{h}$", r"$2^{h}$"])
return fig
......@@ -14,6 +14,7 @@
# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
import utils
from _version import __version__
import sys
from pygments import highlight
from pygments.lexers import get_lexer_by_name
......@@ -105,6 +106,8 @@ class page():
"""
self.add_content("<div class='jumbotron text-center' style='background-color: {}; margin-bottom:0'>\n".format(colour))
self.add_content(" <h1 id={}>{}</h1>\n".format(approximant, title))
self.add_content("<h4><span class='badge badge-info'>Code Version: %s"
"</span></h4>\n" %(__version__), indent=2)
self.add_content("</div>\n")
def _footer(self, user, rundir):
......@@ -294,9 +297,12 @@ class page():
self.add_content("<tbody>\n", indent=6)
for i in contents:
self.add_content("<tr>\n", indent=8)
for j in i:
self.add_content("<th><img src='{}' ".format(self.base_url+"/plots/"+j.split("/")[-1]) +
"alt='No image available' style='width:{}px;'></td>\n".format(1050./len(i)), indent=8)
self.add_content("<td><img src='{}' ".format(self.base_url+"/plots/"+j.split("/")[-1]) +
"alt='No image available' style='width:{}px;' "
"id='{}'></td>\n".format(1050./len(i),j.split("/")[-1][:-4]), indent=10)
self.add_content("</tr>\n", indent=8)
self.add_content("</tbody>\n", indent=6)
self.add_content("</table>\n", indent=4)
......
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