added rates and pop codes to gstlal-ugly

gstlal-ugly/bin/gstlal_lvc_rates_injection_dag

+#!/usr/bin/env python3
+#
+#  Copyright (C) 2018 Sarah Caudill
+#
+#  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 2 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 with program; see the file COPYING. If not, write to the
+#  Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston,
+#  MA  02111-1307  USA
+#
+
+## @file lvc_rates_injection_dag
+# The rates injection dag generator.
+#
+# ### Usage examples
+# - Please add some!
+#
+# ### Command line interface
+
+__author__ = 'Sarah Caudill <sarah.caudill@ligo.org>'
+
+import os
+
+import random
+
+from optparse import OptionParser
+from configparser import ConfigParser
+import tempfile
+from glue import pipeline
+
+import injection_generation_dag
+
+def parse_command_line():
+        parser = OptionParser(description = __doc__)
+        parser.add_option("--tmp-dir", metavar = "directory", help = "Set tmpdir")
+        parser.add_option("--tag", metavar = "str", help = "output dag name")
+        parser.add_option("--bns-params-ini", metavar = "filename", help = "Set path to ini file with inj params.")
+        parser.add_option("--bns-gaussian-params-ini", metavar = "filename", help = "Set path to ini file with inj params.")
+        parser.add_option("--nsbh-params-ini", metavar = "filename", help = "Set path to ini file with inj params.")
+        parser.add_option("--bbh-params-ini", metavar = "filename", help = "Set path to ini file with inj params.")
+        parser.add_option("--nosplit-params-ini", metavar = "filename", help = "Set path to ini file with inj params.")
+        parser.add_option("--analysis-ini", metavar = "filename", help = "Set path to ini file with start and end times of each chunk and reference psd files.")
+        parser.add_option("--condor-command", action = "append", default = [], metavar = "command=value", help = "set condor commands of the form command=value; can be given multiple times")
+
+        options, filenames = parser.parse_args()
+        return options, filenames
+
+try:
+        os.mkdir("logs")
+except:
+        pass
+
+options, filenames = parse_command_line()
+
+ccommands = {}
+for o in options.condor_command:
+	osplit = o.split("=")
+	k = osplit[0]
+	v = "=".join(osplit[1:])
+	ccommands.update([(k, v)])
+
+
+dag = injection_generation_dag.DAG("rates_injections_%s" %(options.tag), logpath = options.tmp_dir)
+analysis_parser = ConfigParser()
+analysis_parser.read(options.analysis_ini)
+
+### BNS Injections
+if options.bns_params_ini is not None:	
+	bnsinjparser = ConfigParser()
+	bnsinjparser.read(options.bns_params_ini)
+	bnsInjJob = injection_generation_dag.InjectionJob("lvc_rates_injections", tag_base = "bns_rates_injections", condor_commands = ccommands)
+	for analysis_section in analysis_parser.sections():
+		for inj_section in bnsinjparser.sections():
+
+			node_opts = {}
+
+			node_opts['random-seed'] = random.randint(1, 100000000)
+
+			for name, value in bnsinjparser.items(inj_section):
+				node_opts[name] = value
+			for name, value in analysis_parser.items(analysis_section):
+				node_opts[name] = value
+
+
+			bnsinjnode = injection_generation_dag.InjectionNode(bnsInjJob, dag, p_node = [], opts = node_opts)
+			bnsinjnode.set_priority(99)
+
+### BNS Gaussian Injections
+if options.bns_gaussian_params_ini is not None:
+	bnsgaussianinjparser = ConfigParser()
+	bnsgaussianinjparser.read(options.bns_gaussian_params_ini)
+	bnsGaussianInjJob = injection_generation_dag.InjectionJob("lvc_rates_injections", tag_base = "bns_gaussian_rates_injections", condor_commands = ccommands)
+	for analysis_section in analysis_parser.sections():
+		for inj_section in bnsgaussianinjparser.sections():
+
+			node_opts = {}
+
+			node_opts['random-seed'] = random.randint(1, 100000000)
+
+			for name, value in bnsgaussianinjparser.items(inj_section):
+				node_opts[name] = value
+			for name, value in analysis_parser.items(analysis_section):
+				node_opts[name] = value
+
+
+			bnsgaussianinjnode = injection_generation_dag.InjectionNode(bnsGaussianInjJob, dag, p_node = [], opts = node_opts)
+			bnsgaussianinjnode.set_priority(99)
+
+
+### NSBH Injections
+if options.nsbh_params_ini is not None:
+	nsbhinjparser = ConfigParser()
+	nsbhinjparser.read(options.nsbh_params_ini)
+	nsbhInjJob = injection_generation_dag.InjectionJob("lvc_rates_injections", tag_base = "nsbh_rates_injections", condor_commands = ccommands)
+	for analysis_section in analysis_parser.sections():
+		for inj_section in nsbhinjparser.sections():
+
+			node_opts = {}
+
+			node_opts['random-seed'] = random.randint(1, 100000000)
+
+			for name, value in nsbhinjparser.items(inj_section):
+				node_opts[name] = value
+			for name, value in analysis_parser.items(analysis_section):
+				node_opts[name] = value
+
+
+			nsbhinjnode = injection_generation_dag.InjectionNode(nsbhInjJob, dag, p_node = [], opts = node_opts)
+			nsbhinjnode.set_priority(99)
+
+
+### BBH Injections
+if options.bbh_params_ini is not None:
+	bbhinjparser = ConfigParser()
+	bbhinjparser.read(options.bbh_params_ini)
+	bbhInjJob = injection_generation_dag.InjectionJob("lvc_rates_injections", tag_base = "bbh_rates_injections", condor_commands = ccommands)
+	for analysis_section in analysis_parser.sections():
+		for inj_section in bbhinjparser.sections():
+
+			node_opts = {}
+
+			node_opts['random-seed'] = random.randint(1, 100000000)
+
+			for name, value in bbhinjparser.items(inj_section):
+				node_opts[name] = value
+			for name, value in analysis_parser.items(analysis_section):
+				node_opts[name] = value
+
+
+			bbhinjnode = injection_generation_dag.InjectionNode(bbhInjJob, dag, p_node = [], opts = node_opts)
+			bbhinjnode.set_priority(99)
+
+### NOSPLIT Injections
+if options.nosplit_params_ini is not None:
+	nosplitinjparser = ConfigParser()
+	nosplitinjparser.read(options.nosplit_params_ini)
+	nosplitInjJob = injection_generation_dag.InjectionJob("lvc_rates_injections", tag_base = "nosplit_rates_injections", condor_commands = ccommands)
+	for analysis_section in analysis_parser.sections():
+		for inj_section in nosplitinjparser.sections():
+
+			node_opts = {}
+
+			node_opts['random-seed'] = random.randint(1, 100000000)
+
+			for name, value in nosplitinjparser.items(inj_section):
+				node_opts[name] = value
+			for name, value in analysis_parser.items(analysis_section):
+				node_opts[name] = value
+
+			nosplitinjnode = injection_generation_dag.InjectionNode(nosplitInjJob, dag, p_node = [], opts = node_opts)
+			nosplitinjnode.set_priority(99)
+
+			#all_inj_nodes[(analysis_section, inj_section)].append(inj_nodes)
+
+dag.write_sub_files()
+dag.write_dag()
+dag.write_script()

gstlal-ugly/bin/gstlal_lvc_rates_injections

+#!/usr/bin/env python3
+#
+#  Copyright (C) 2017 Jolien Creighton
+#
+#  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 2 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 with program; see the file COPYING. If not, write to the
+#  Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston,
+#  MA  02111-1307  USA
+#
+
+## @file lvc_rates_injections
+# The rates injection generator.
+#
+# ### Usage examples
+# - Please add some!
+#
+# ### Command line interface
+
+__author__ = 'Jolien Creighton <jolien.creighton@ligo.org>'
+
+import sys
+from optparse import OptionParser
+import numpy
+import numpy.random
+import scipy.integrate
+from ligo.lw import utils
+from ligo.lw import ligolw
+from ligo.lw import table
+from ligo.lw import lsctables
+from ligo.lw.utils import process as ligolw_process
+import lal
+import lalsimulation
+
+import cosmology_utils as cutils
+import injection_utils as injutils
+
+def parse_command_line():
+    parser = OptionParser(description = __doc__)
+    parser.add_option("--gps-start-time", type = "int", help = "Start time of the injection data to be created.")
+    parser.add_option("--gps-end-time", type = "int", help = "End time of the injection data to be created.")
+    parser.add_option("--use-segments", default = False, action = "store_true", help = "Not currently implemented. Cut on segments.")
+    parser.add_option("--output-tag", metavar = "tag", default = "injections", help = "Set the string to write the injection xml.")
+
+    # Quantities for constructing injection parameter distributions
+    parser.add_option("--random-seed", default = 7, type = "int", help = "Set the seed of the random number generator to get reproducible results, default = 7.")
+    parser.add_option("--max-redshift", type = "float", help = "Set the maximum redshift at which injections will be placed.")
+    parser.add_option("--redshift-power", type = "float", help = "Power of (1+z) to multiply the constant-comoving-rate distribution by.")
+    parser.add_option("--time-step", type = "float", help = "Set the time step interval between injections.")
+    parser.add_option("--time-interval", type = "float", help = "Set the size of the interval used to randomly place injections around the --time-step. For example, 10 will result in +-10s interval on either side of time-step.")
+    parser.add_option("--randomize-start-time", default = False, action = "store_true", help = "Add a small random jitter to the gps-start-time. In the case of many injection sets, this results in jittering the time-intervals a bit.")
+    parser.add_option("--mass-distribution", default = None, type = "str", help = "Set the desired component mass distribution. Select one of the following: " +
+                      "1. IMF_PAIR - First component from Salpeter IMF distribution and second drawn uniformly between min_mass and mass of first component. " +
+                      "2. UNIFORM_PAIR - Both components drawn from uniform distribution so that min_mass <= m1,m2 < max_mass. " +
+                      "3. UNIFORMLNM_PAIR - Both components drawn from uniform-in-log distribution. " +
+                      "4. NORMAL_PAIR - Both components drawn from normal distribution of mean mean_mass and standard deviation sigma_mass." +
+                      "5. UNIFORM_DISTINCT - Components will be drawn from two distinct uniform mass distributions. Must specify min_mass1, max_mass1, " +
+                      "min_mass2, max_mass2." +
+                      "6. UNIFORMLNM_DISTINCT - Components will be drawn from two distinct uniform-in-log mass distributions. Must specify min_mass1, max_mass1, " +
+                      "min_mass2, max_mass2." +
+                      "7. UNIFORM_UNIFORMLNM_DISTINCT - Component 1 will be drawn from a uniform mass distribution and component 2 will be drawn " +
+                      "from a uniform-in-log mass distribution. Must specify min_mass1, max_mass1, min_mass2, max_mass2." +
+                      "8. NORMAL_IMF_DISTINCT - Component 1 will be drawn from a normal mass distribution and component 2 will be drawn " +
+                      "from a Salpeter IMF distribution." + 
+                      "9. UNIFORM_IMF_DISTINCT - Primary drawn from a uniform mass distribution, secondary drawn from a Salpeter IMF distribution. " + 
+                      "10. POWER_PAIR - Primary from power law, secondary from another power law conditional on m2 < m1. " +
+                      "Default power laws m1**2.35, m2**2 reproduce the IMF_M2SQ distribution.")
+    parser.add_option("--min-mass", default = None, type = "float", help = "Set the minimum component mass of the two compact objects.")
+    parser.add_option("--max-mass", default = None, type = "float", help = "Set the maximum component mass of the two compact objects.")
+    parser.add_option("--max-mtotal", default = None, type = "float", help = "Set a cutoff on the allowed total mass. Currently only implemented for " +
+                      "IMF_PAIR, UNIFORMLNM_PAIR, UNIFORM_DISTINCT, UNIFORMLNM_DISTINCT, UNIFORM_UNIFORMLNM_DISTINCT.")
+    parser.add_option("--mean-mass", default = None, type = "float", help = "Mean of normal mass distribution. Required if using NORMAL_PAIR.")
+    parser.add_option("--sigma-mass", default = None, type = "float", help = "Stdev of normal mass distribution. Required if using NORMAL_PAIR.")
+    parser.add_option("--min-mass1", default = None, type = "float", help = "Set the minimum component mass of the lighter of two compact objects.")
+    parser.add_option("--max-mass1", default = None, type = "float", help = "Set the maximum component mass of the lighter of two compact objects.")
+    parser.add_option("--min-mass2", default = None, type = "float", help = "Set the minimum component mass of the heavier of two compact objects.")
+    parser.add_option("--max-mass2", default = None, type = "float", help = "Set the maximum component mass of the heavier of two compact objects.")
+    parser.add_option("--mass1-power", default = 2.35, type = "float", help = "Power law index for primary mass. Default 2.35")
+    parser.add_option("--mass2-power", default = 2.0, type = "float", help = "Power law index for secondary mass. Default 2")
+    parser.add_option("--spin-distribution", type = "str", help = "Set the desired component spin distribution. Select one of the following: " +
+                      "1. ALIGNED - Spinz uniformly random in (-max_spin, +max_spin) for each component. " +
+                      "2. ISOTROPIC - (s_x, s_y, s_z) isotropically distributed with uniform magnitude distribution for each component." +
+                      "3. ALIGNED_ALIGNED - Components uniformly random in (-max_spin1, +max_spin1) and (-max_spin2, +max_spin2)" +
+                      "4. ISOTROPIC_ALIGNED - Component 2 spinz uniformly random in (-max_spin1, +max_spin1), component 1 isotropically distributed." +
+                      "5. ISOTROPIC_ISOTROPIC - Both components isotropically distributed with uniform magnitude distribution." + 
+                      "6. ALIGNED_EQUAL - Spinz uniform in (-max_spin, +max_spin) for component 1, spin2z = spin1z." +
+                      "7. SALPETER_PRIMARY_SPIN - Both components (1/2 chi_max) ln(chi_max/s1z), needs --max-chi")
+    parser.add_option("--max-spin", type = "float", help = "Set the maximum component spin magnitude. Required if using ALIGNED, ALIGNED_EQUAL or ISOTROPIC.")
+    parser.add_option("--max-chi", default = 0.99, type = "float", help = "Maximum chi spin magnitude. Takes 0.99. Required if using SALPETER_PRIMARY_SPIN.")
+    parser.add_option("--max-spin1", type = "float", help = "Set the maximum component spin magnitude of mass_1. Required if using ALIGNED_ISOTROPIC.")
+    parser.add_option("--max-spin2", type = "float", help = "Set the maximum component spin magnitude of mass_2. Required if using ALIGNED_ISOTROPIC.")
+    parser.add_option("--waveform", type = "str", help = "Set the waveform family to populate the sim_inspiral table.")
+
+    # Quantities for computing expected SNRs
+    parser.add_option("--snr-calculation", type = "str", help = "Set the desired method for computing the approximate injection snr. Select one of the following: " +
+                      "1. OPTIMALLY_ORIENTED_1MPC - Assumes inclination = 0, distance = 1Mpc. Computes only hp." +
+                      "2. INJ_PARAMS - Computes hp and hc using injections' parameters." +
+                      "3. GENERIC - Computes hp and hc using injections' parameters and can take FD and TD waveforms.")
+    parser.add_option("--snr-threshold", default = 3.0, type = "float", help = "Set the single detector snr threshold to determine " +
+                      "whether an injection is expected to be found or not, default = 3.0")
+    parser.add_option("--min-frequency", default = 15.0, type = "float", help = "Lower frequency cutoff for SNR calculation, default = 15 Hz.")
+    parser.add_option("--max-frequency", default = 1500.0, type = "float", help = "Upper frequency cutoff for SNR calculation, default = 1500 Hz.")
+    parser.add_option("--delta-frequency", default = 1.0, type = "float", help = "Step in frequency for SNR calculation, default = 1 Hz.")
+    parser.add_option("--approximant", default = 'SEOBNRv2_ROM_DoubleSpin', type = "str", help = "Set the approximant waveform to use in the SNR" +
+                      "calculation, default = SEOBNRv2_ROM_DoubleSpin")
+#    parser.add_option("--reference-spectrum-file", metavar = "filename", help = "Set the full path to the file containing two columns for the frequency " +
+#                      "and the spectral data. It can be ASD or PSD data. The code will try to figure out which data you have provided.")
+    parser.add_option("--h1-reference-spectrum-file", metavar = "filename", help = "Full path to the H1 PSD as txt")
+    parser.add_option("--l1-reference-spectrum-file", metavar = "filename", help = "Full path to the L1 PSD as txt")
+    parser.add_option("--v1-reference-spectrum-file", metavar = "filename", help = "Full path to the V1 PSD as txt")
+    #parser.add_option("--k1-reference-spectrum-file", metavar = "filename", help = "Full path to the K1 PSD as txt")
+
+    parser.add_option("-v", "--verbose", action = "store_true", help = "Be verbose.")
+
+    options, filenames = parser.parse_args()
+
+    process_params = options.__dict__.copy()
+
+    required_options = ("gps_start_time", "gps_end_time", "max_redshift", "time_interval", "time_step", "mass_distribution",
+			"spin_distribution", "waveform", "approximant", 
+			"h1_reference_spectrum_file", "l1_reference_spectrum_file", "v1_reference_spectrum_file")
+			#"h1_reference_spectrum_file", "l1_reference_spectrum_file", "v1_reference_spectrum_file", "k1_reference_spectrum_file")
+    missing_options = [option for option in required_options if getattr(options, option) is None]
+    if missing_options:
+        raise ValueError("missing required option(s) %s" % ", ".join("--%s" % option.replace("_", "-") for option in missing_options))
+
+    return options, process_params, filenames
+
+
+options, process_params, filenames = parse_command_line()
+
+# seed the random number generator to get reproducible results
+numpy.random.seed(options.random_seed)
+
+# psd and related quantities for computing SNRs
+approximant = lalsimulation.SimInspiralGetApproximantFromString(options.approximant)
+h1_psd = lal.CreateREAL8FrequencySeries(
+    'PSD',
+    lal.LIGOTimeGPS(0),
+    0.0,
+    options.delta_frequency,
+    lal.SecondUnit,
+    int(round(options.max_frequency / options.delta_frequency)) + 1)
+lalsimulation.SimNoisePSDFromFile(h1_psd, options.min_frequency, options.h1_reference_spectrum_file)
+# FIXME: Write code to determine whether we have ASD or PSD
+# SimNoisePSDFromFile expects ASD in the file, but this one
+# contains the PSD, so take the square root
+#h1_psd.data.data = h1_psd.data.data ** 0.5
+
+l1_psd = lal.CreateREAL8FrequencySeries(
+    'PSD',
+    lal.LIGOTimeGPS(0),
+    0.0,
+    options.delta_frequency,
+    lal.SecondUnit,
+    int(round(options.max_frequency / options.delta_frequency)) + 1)
+lalsimulation.SimNoisePSDFromFile(l1_psd, options.min_frequency, options.l1_reference_spectrum_file)
+#l1_psd.data.data = l1_psd.data.data ** 0.5
+
+v1_psd = lal.CreateREAL8FrequencySeries(
+    'PSD',
+    lal.LIGOTimeGPS(0),
+    0.0,
+    options.delta_frequency,
+    lal.SecondUnit,
+    int(round(options.max_frequency / options.delta_frequency)) + 1)
+lalsimulation.SimNoisePSDFromFile(v1_psd, options.min_frequency, options.v1_reference_spectrum_file)
+#v1_psd.data.data = v1_psd.data.data ** 0.5
+
+#k1_psd = lal.CreateREAL8FrequencySeries(
+#    'PSD',
+#    lal.LIGOTimeGPS(0),
+#    0.0,
+#    options.delta_frequency,
+#    lal.SecondUnit,
+#    int(round(options.max_frequency / options.delta_frequency)) + 1)
+#lalsimulation.SimNoisePSDFromFile(k1_psd, options.min_frequency, options.k1_reference_spectrum_file)
+# k1_psd.data.data = k1_psd.data.data ** 0.5
+
+omega = cutils.get_cosmo_params()
+
+VT = cutils.surveyed_spacetime_volume(options.gps_start_time, options.gps_end_time, options.max_redshift, omega)
+print ('surveyed spacetime volume: %g Gpc^3 yr' % VT)
+
+xmldoc = ligolw.Document()
+xmldoc.appendChild(ligolw.LIGO_LW())
+sim_inspiral_table = lsctables.New(lsctables.SimInspiralTable)
+
+process = ligolw_process.register_to_xmldoc(
+    xmldoc,
+    program="lvc_rates_injections",
+    paramdict=process_params,
+    comment="Injection parameter generator for rates calculations")
+
+random_sim_inspiral_rows = iter(injutils.draw_sim_inspiral_row(process_params, h1_psd, l1_psd, v1_psd, omega))
+
+for accept, reject, row in random_sim_inspiral_rows:
+
+    row.process_id = process.process_id
+    sim_inspiral_table.append(row)
+    if (accept % 100) == 0:
+        print ('(accept,reject)=(%d,%d)\r' % (accept, reject))
+        sys.stdout.flush()
+
+ligolw_process.set_process_end_time(process)
+
+print ('(accept,reject)=(%d,%d)' % (accept, reject))
+sys.stdout.flush()
+
+acceptance_rate = float(accept)/float(accept + reject)
+print ('%s acceptance rate: %g' % (options.output_tag, acceptance_rate))
+if process_params['redshift_power'] is None:  # constant rate in comoving volume-time
+    print ('%s actual surveyed spacetime volume: %g Gpc^3 yr' % (options.output_tag, VT * acceptance_rate))
+sys.stdout.flush()
+
+VT_params = {
+    'accept' : accept,
+    'reject' : reject,
+    'acceptance_rate' : acceptance_rate,
+    'VT' : VT
+    }
+
+ligolw_process.register_to_xmldoc(
+    xmldoc,
+    program="lvc_rates_injections_vtparams",
+    paramdict=VT_params,
+    comment="Acceptance and rejection related statistics for VT calculations.")
+
+xmldoc.childNodes[-1].appendChild(sim_inspiral_table)
+output_file = options.output_tag + '-' + str(options.gps_start_time) + '-' + str(options.gps_end_time) + '.xml.gz'
+utils.write_filename(xmldoc, output_file,  verbose = options.verbose)
+xmldoc.unlink()

gstlal-ugly/bin/gstlal_lvc_rates_plot_injections

+#!/usr/bin/env python3
+
+import os
+import sys
+import math
+import matplotlib
+matplotlib.rcParams.update({
+	"font.size": 10.0,
+	"axes.titlesize": 10.0,
+	"axes.labelsize": 10.0,
+	"xtick.labelsize": 8.0,
+	"ytick.labelsize": 8.0,
+	"legend.fontsize": 8.0,
+	"figure.dpi": 600,
+	"savefig.dpi": 600,
+	"text.usetex": True
+})
+from matplotlib import figure
+from matplotlib.backends.backend_agg import FigureCanvasAgg as FigureCanvas
+import numpy
+
+from optparse import OptionParser
+
+from ligo.lw import utils
+from ligo.lw import ligolw
+from ligo.lw import table
+from ligo.lw import lsctables
+from ligo.lw.utils import process as ligolw_process
+
+import astropy.cosmology as cosmo
+import astropy.stats as ast
+import astropy.units as u
+
+__author__ = "Sarah Caudill <sarah.caudill@ligo.org>"
+__version__ = "git id %s" % ""	# FIXME
+__date__ = ""	# FIXME
+
+class ContentHandler(ligolw.LIGOLWContentHandler):
+    pass
+
+lsctables.use_in(ContentHandler)
+
+
+def create_plot(x_label = None, y_label = None, width = 165.0, aspect = None):
+    if aspect is None:
+        aspect = (1 + math.sqrt(5)) / 2
+    fig = figure.Figure()
+    FigureCanvas(fig)
+    fig.set_size_inches(width / 25.4, width / 25.4 / aspect)
+    axes = fig.gca()
+    axes.grid(True)
+    if x_label is not None:
+        axes.set_xlabel(x_label)
+    if y_label is not None:
+        axes.set_ylabel(y_label)
+    return fig, axes
+
+def parse_command_line():
+    parser = OptionParser(
+        description = __doc__
+    )
+
+    parser.add_option("-v", "--verbose", action = "store_true", help = "Be verbose.")
+    options, urls = parser.parse_args()
+
+    return options, urls
+
+options, urls = parse_command_line()
+
+for file_path in urls:
+
+    tag = os.path.basename(file_path).split(".")[0]
+#    fp = open(file_path, 'r')
+#    xmldoc, digest = utils.load_fileobj(fp, gz=True, contenthandler=ContentHandler)
+#    sim_table = table.get_table(xmldoc, 'sim_inspiral')
+    xmldoc = utils.load_filename(sys.argv[1], False, contenthandler=ContentHandler)
+    sim_table = table.get_table(xmldoc, lsctables.SimInspiralTable.tableName)
+
+    mass1 = []
+    mass2 = []
+    spin1x = []
+    spin1y = []
+    spin1z = []
+    spin2x = []
+    spin2y = []
+    spin2z = []
+    distance = []
+    time = []
+    alpha3 = []
+    for sim in sim_table:
+        mass1.append(sim.mass1)
+        mass2.append(sim.mass2)
+        spin1x.append(sim.spin1x)
+        spin1y.append(sim.spin1y)
+        spin1z.append(sim.spin1z)
+        spin2x.append(sim.spin2x)
+        spin2y.append(sim.spin2y)
+        spin2z.append(sim.spin2z)
+        distance.append(sim.distance)
+        time.append(sim.geocent_end_time + 10**-9*sim.geocent_end_time_ns)
+        alpha3.append(sim.alpha3)
+
+    fig, axes = create_plot('det mass1','Count')
+    axes.set_title(r"Plot")
+    axes.hist(mass1, 100)
+    fig.savefig('det_mass1_hist_%s.png' % str(tag))
+
+    fig, axes = create_plot('det mass2','Count')
+    axes.set_title(r"Plot")
+    axes.hist(mass2, 100)
+    fig.savefig('det_mass2_hist_%s.png' % str(tag))
+
+    fig, axes = create_plot('det mtotal','Count')
+    axes.set_title(r"Plot")
+    axes.hist(numpy.array(mass1) + numpy.array(mass2), 100)
+    fig.savefig('det_mtotal_hist_%s.png' % str(tag))
+
+    fig, axes = create_plot('src mass1','Count')
+    axes.set_title(r"Plot")
+    axes.hist(numpy.array(mass1)/(1.0 + numpy.array(alpha3)), 100, log=True)
+    fig.savefig('src_mass1_loghist_%s.png' % str(tag))
+
+    fig, axes = create_plot('src mass2','Count')
+    axes.set_title(r"Plot")
+    axes.hist(numpy.array(mass2)/(1.0 + numpy.array(alpha3)), 100, log=True)
+    fig.savefig('src_mass2_loghist_%s.png' % str(tag))
+
+    fig, axes = create_plot('src mtotal','Count')
+    axes.set_title(r"Plot")
+    axes.hist((numpy.array(mass1) + numpy.array(mass2))/(1.0 + numpy.array(alpha3)), 100, log=True)
+    fig.savefig('src_mtotal_loghist_%s.png' % str(tag))
+
+    fig, axes = create_plot('det mass1','Count')
+    axes.set_title(r"Plot")
+    axes.hist(mass1, 100, log=True)
+    fig.savefig('det_mass1_loghist_%s.png' % str(tag))
+
+    fig, axes = create_plot('det mass2','Count')
+    axes.set_title(r"Plot")
+    axes.hist(mass2, 100, log=True)
+    fig.savefig('det_mass2_loghist_%s.png' % str(tag))
+
+    fig, axes = create_plot('det mtotal','Count')
+    axes.set_title(r"Plot")
+    axes.hist(numpy.array(mass1) + numpy.array(mass2), 100, log=True)
+    fig.savefig('det_mtotal_loghist_%s.png' % str(tag))
+
+    fig, axes = create_plot('src mass1','Count')
+    axes.set_title(r"Plot")
+    axes.hist(numpy.array(mass1)/(1.0 + numpy.array(alpha3)), 100)
+    fig.savefig('src_mass1_hist_%s.png' % str(tag))
+
+    fig, axes = create_plot('src mass2','Count')
+    axes.set_title(r"Plot")
+    axes.hist(numpy.array(mass2)/(1.0 + numpy.array(alpha3)), 100)
+    fig.savefig('src_mass2_hist_%s.png' % str(tag))
+
+    fig, axes = create_plot('src mtotal','Count')
+    axes.set_title(r"Plot")
+    axes.hist((numpy.array(mass1) + numpy.array(mass2))/(1.0 + numpy.array(alpha3)), 100)
+    fig.savefig('src_mtotal_hist_%s.png' % str(tag))
+
+    fig, axes = create_plot('spin1x','Count')
+    axes.set_title(r"Plot")
+    axes.hist(spin1x, 100)
+    fig.savefig('spin1x_hist_%s.png' % str(tag))
+
+    fig, axes = create_plot('spin1y','Count')
+    axes.set_title(r"Plot")
+    axes.hist(spin1y, 100)
+    fig.savefig('spin1y_hist_%s.png' % str(tag))
+
+    fig, axes = create_plot('spin1z','Count')
+    axes.set_title(r"Plot")
+    axes.hist(spin1z, 100)
+    fig.savefig('spin1z_hist_%s.png' % str(tag))
+
+    fig, axes = create_plot('spin2x','Count')
+    axes.set_title(r"Plot")
+    axes.hist(spin2x, 100)
+    fig.savefig('spin2x_hist_%s.png' % str(tag))
+
+    fig, axes = create_plot('spin2y','Count')
+    axes.set_title(r"Plot")
+    axes.hist(spin2y, 100)
+    fig.savefig('spin2y_hist_%s.png' % str(tag))
+
+    fig, axes = create_plot('spin2z','Count')
+    axes.set_title(r"Plot")
+    axes.hist(spin2z, 100)
+    fig.savefig('spin2z_hist_%s.png' % str(tag))
+
+    fig, axes = create_plot('distance','Count')
+    axes.set_title(r"Plot")
+    axes.hist(distance, 100)
+    fig.savefig('distance_hist_%s.png' % str(tag))
+
+    fig, axes = create_plot('Time','Count')
+    axes.set_title(r"Plot")
+    axes.hist(time, 100)
+    fig.savefig('time_hist_%s.png' % str(tag))
+
+    fig, axes = create_plot('Redshift','Count')
+    axes.hist(alpha3, 100, normed=True)
+    zs = numpy.linspace(0, max(alpha3), 1000)
+    dvcdzo1pzs = cosmo.Planck13.differential_comoving_volume(zs).to(u.Gpc**3/u.sr).value/(1+zs)
+    propernorm = numpy.trapz(dvcdzo1pzs, zs)
+    axes.plot(zs, dvcdzo1pzs/propernorm, label='Cosmological Distribution')
+    fig.savefig('cosmological_distribution_%s.png' % str(tag))

gstlal-ugly/python/cosmology_utils.py

+#  Copyright (C) 2017 Jolien Creighton
+#
+#  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 2 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 with program; see the file COPYING. If not, write to the
+#  Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston,
+#  MA  02111-1307  USA
+
+## @file
+# The python module for utilities needed for cosmological calculations.
+
+#
+# =============================================================================
+#
+#                                   Preamble
+#
+# =============================================================================
+#
+
+import numpy
+import scipy.integrate
+
+import lal
+
+
+#
+# =============================================================================
+#
+#                      Cosmological and Other Constants
+#
+# =============================================================================
+#
+
+# FIXME: go from astropy instead
+def get_cosmo_params():
+    # From Planck2015, Table IV
+    omega = lal.CreateCosmologicalParametersAndRate().omega
+    lal.SetCosmologicalParametersDefaultValue(omega)
+    omega.h = 0.679
+    omega.om = 0.3065
+    omega.ol = 0.6935
+    omega.ok = 1.0 - omega.om - omega.ol
+    omega.w0 = -1.0
+    omega.w1 = 0.0
+    omega.w2 = 0.0
+
+    return omega
+
+#
+# =============================================================================
+#
+#                               Cosmology Utilities
+#
+# =============================================================================
+#
+
+
+def surveyed_spacetime_volume(gps_start_time, gps_end_time, max_redshift, omega):
+    '''
+    Returns the total spacetime volume surveyed:
+
+        <VT> = T \int dz \frac{dV_c}{dz} \frac{1}{1+z}
+
+    Results are given in units Gpc^3 yr(Julian).
+    '''
+
+    # Note: LAL's cosmology routines returns distances in Mpc
+    def integrand(z, omega):
+        '''
+        Returns the integrand
+
+            (1 + z) D_A^2(z) / E(z)
+
+        in units of Mpc^2.  Multiply the integral by 4 * pi * D_H
+        to get the desired integral.
+        '''
+
+        return (1.0 + z) * lal.AngularDistance(omega, z)**2 * lal.HubbleParameter(z, omega)
+
+    I, _ = scipy.integrate.quad(integrand, 0.0, max_redshift, args=omega)
+
+    # multiply by remaining factors and scale to Gpc^3
+    V = 4.0 * numpy.pi * lal.HubbleDistance(omega) * I / (1e3)**3
+    
+    # surveyed time in Julian years
+    T = (gps_end_time - gps_start_time) / lal.YRJUL_SI
+
+    return V * T
+
+

gstlal-ugly/python/injection_generation_dag.py

+from optparse import OptionParser
+import tempfile
+from glue import pipeline
+from gstlal.dags.util import pipeline_dot_py_append_opts_hack
+
+class DAG(pipeline.CondorDAG):
+	def __init__(self, name, logpath = None):
+		self.basename = name.replace(".dag","")
+		tempfile.tempdir = logpath
+		tempfile.template = self.basename + '.dag.log.'
+		logfile = tempfile.mktemp()
+		fh = open( logfile, "w" )
+		fh.close()
+		pipeline.CondorDAG.__init__(self,logfile)
+		self.set_dag_file(self.basename)
+		self.jobsDict = {}
+		self.output_cache = []
+
+	def add_node(self, node):
+		node.set_retry(1)
+		node.add_macro("macronodename", node.get_name())
+		pipeline.CondorDAG.add_node(self, node)
+
+class InjectionJob(pipeline.CondorDAGJob):
+	def __init__(self, executable, tag_base = None, universe="vanilla", condor_commands = {}, **kwargs):
+		self.__prog__ = tag_base
+		self.__executable = executable
+		self.__universe = universe
+		pipeline.CondorDAGJob.__init__(self, self.__universe, self.__executable)
+		self.add_condor_cmd('getenv','True')
+		self.tag_base = tag_base
+		self.set_sub_file(tag_base+'.sub')
+		self.set_stdout_file('logs/$(macronodename)-$(cluster)-$(process).out')
+		self.set_stderr_file('logs/$(macronodename)-$(cluster)-$(process).err')
+		self.number = 1
+		# make an output directory for files
+		self.output_path = tag_base
+		try:
+			os.mkdir(self.output_path)
+		except:
+			pass
+		for cmd,val in condor_commands.items():
+			self.add_condor_cmd(cmd, val)
+
+class InjectionNode(pipeline.CondorDAGNode):
+	def __init__(self, job, dag, p_node=[], opts={}):
+		pipeline.CondorDAGNode.__init__(self, job)
+		for p in p_node:
+			self.add_parent(p)
+
+		for opt, val in opts.items():
+			if val is None:
+				continue # not the same as val = '' which is allowed
+			if not hasattr(val, "__iter__") or isinstance(val, str): # catches list like things and strings
+				if opt == "":
+					self.add_var_arg(val)
+				else:
+					self.add_var_opt(opt, val)
+			# Must be an iterable, but not str
+			else:
+				if opt == "":
+					[self.add_var_arg(a) for a in val]
+				else:
+					self.add_var_opt(opt, pipeline_dot_py_append_opts_hack(opt, val))
+
+		self.set_name("%s_%04X" % (job.tag_base, job.number))
+		job.number += 1
+		dag.add_node(self)