Skip to content
Snippets Groups Projects
Commit 926ea2fd authored by ChiWai Chan's avatar ChiWai Chan
Browse files

gstlal_inspiral_calc_snr & svd_bank_snr.py:

	--move template making step to svd_bank_snr.FIR_SNR.make_template()
	--calculte the template duration instead of inputting template duration
	--correct snr's epoch by adding time offset due to conditioning of templates
parent 6366807e
No related branches found
No related tags found
No related merge requests found
......@@ -43,7 +43,6 @@ Typical Usages:
2. Template options:
--table (require)
--approximant (require)
--template-duration (require)
--sample-rate (default = 2048Hz)
--f-low (default = 10)
--f-high (optional)
......@@ -67,7 +66,6 @@ import numpy
from optparse import OptionParser, OptionGroup, IndentedHelpFormatter
import os
from gstlal import cbc_template_fir
from gstlal import datasource
from gstlal import pipeparts
from gstlal import reference_psd
......@@ -209,27 +207,16 @@ def parse_command_line():
missing_required_options.append("--outdir")
if options.instrument is None:
missing_required_options.append("--instrument")
if options.template_duration is None:
missing_required_options.append("--template-duration")
if options.template_psd is None:
missing_required_options.append("--template-psd")
# Raise VauleError is required option(s) is/are missing
if missing_required_options:
raise ValueError("Missing required option(s) %s" % ", ".join(sorted(missing_required_options)))
xmldoc = ligolw_utils.load_url(options.table, contenthandler = ContentHandler, verbose = options.verbose)
template_table = table.get_table(xmldoc, lsctables.SnglInspiralTable.tableName)
row = [row for row in template_table if row.ifo == options.instrument]
if len(row) != 1 :
raise ValueError("Expecting only one template for --instrument=%s or cannot find template for --instrument=%s" %(options.instrument, options.instrument))
template_psd = lal.series.read_psd_xmldoc(ligolw_utils.load_url(options.template_psd, contenthandler = lal.series.PSDContentHandler))
# work around for building a single whitened template
time_slice = numpy.array([(options.sample_rate, 0, options.template_duration)], dtype = [("rate", "int"),("begin", "float"), ("end", "float")])
workspace = cbc_template_fir.templates_workspace(template_table, options.approximant, template_psd[options.instrument], options.f_low, time_slice, autocorrelation_length = None, fhigh = options.f_high)
template, autocorrelation, sigma = workspace.make_whitened_template(row[0])
return options, gw_data_source_info, template, psd
return options, gw_data_source_info, template_table, psd
# Unknown mode
else:
......@@ -257,9 +244,16 @@ if options.mode == 0:
verbose = options.verbose
)
for index, snr in enumerate(lloid_snr(COMPLEX = options.complex, row_number = options.row_number, start = options.start, end = options.end)):
snrdict = {options.instrument: [snr]}
svd_bank_snr.write_url(svd_bank_snr.make_xmldoc(snrdict), os.path.join(options.outdir, "%s-SNR_%d-%d-%d.xml.gz" % (options.instrument, index, int(snr.epoch), int(snr.data.length * snr.deltaT))), verbose = options.verbose)
if options.row_number is None:
for index, snr in enumerate(lloid_snr(COMPLEX = options.complex, row_number = options.row_number, start = options.start, end = options.end)):
snr.epoch += bank.sngl_inspiral_table[index].end
snrdict = {options.instrument: [snr]}
svd_bank_snr.write_url(svd_bank_snr.make_xmldoc(snrdict), os.path.join(options.outdir, "%s-SNR_%d-%d-%d.xml.gz" % (options.instrument, index, int(snr.epoch), int(snr.data.length * snr.deltaT))), verbose = options.verbose)
else:
lloidsnr = lloid_snr(COMPLEX = options.complex, row_number = options.row_number, start = options.start, end = options.end)
lloidsnr[0].epoch += bank.sngl_inspiral_table[options.row_number].end
snrdict = {options.instrument: lloidsnr}
svd_bank_snr.write_url(svd_bank_snr.make_xmldoc(snrdict), os.path.join(options.outdir, "%s-SNR_%d-%d-%d.xml.gz" % (options.instrument, options.row_number, int(lloidsnr[0].epoch), int(lloidsnr[0].data.length * lloidsnr[0].deltaT))), verbose = options.verbose)
#
# uncomment to save all snrs in one single XML
#
......@@ -267,6 +261,9 @@ if options.mode == 0:
#svd_bank_snr.write_url(svd_bank_snr.make_xmldoc(snrdict), os.path.join(options.outdir, "%s-SNR-%d-%d.xml.gz" % (options.instrument, int(snrdict.[options.instrument][0].epoch), int(snrdict[options.instrument][0].data.length * snrdict[options.instrument][0].deltaT))), verbose = options.verbose)
elif options.mode == 1:
template_table = template
template, time_offset = svd_bank_snr.FIR_SNR.make_template(template_table, options.template_psd, options.sample_rate, options.approximant, options.instrument, options.f_low, f_high = options.f_high, verbose = options.verbose)
#FIXME: proper handle for latency
fir_snr = svd_bank_snr.FIR_SNR(
gw_data_source_info,
......@@ -281,5 +278,7 @@ elif options.mode == 1:
verbose = options.verbose
)
snrdict = {options.instrument : fir_snr(COMPLEX = options.complex, start = options.start, end = options.end)}
firsnr = fir_snr(COMPLEX = options.complex, start = options.start, end = options.end)
firsnr[0].epoch += time_offset
snrdict = {options.instrument : firsnr}
svd_bank_snr.write_url(svd_bank_snr.make_xmldoc(snrdict),os.path.join(options.outdir, "%s-SNR-%d-%d.xml.gz" % (options.instrument, int(snrdict[options.instrument][0].epoch), int(snrdict[options.instrument][0].data.length * snrdict[options.instrument][0].deltaT))), verbose = options.verbose)
......@@ -14,6 +14,7 @@ from gi.repository import GObject, Gst, GstAudio
GObject.threads_init()
Gst.init(None)
from gstlal import cbc_template_fir
from gstlal import datasource
from gstlal import lloidparts
from gstlal import multirate_datasource
......@@ -25,11 +26,13 @@ from gstlal import simplehandler
import lal
from lal import LIGOTimeGPS
import lal.series
import lalsimulation as lalsim
from ligo.lw import array as ligolw_array
from ligo.lw import ligolw
from ligo.lw import param as ligolw_param
from ligo.lw import utils as ligolw_utils
from ligo.lw import lsctables
@ligolw_array.use_in
@ligolw_param.use_in
......@@ -296,6 +299,24 @@ class FIR_SNR(SNR_Pipeline):
self.snr_info["data"] = numpy.vectorize(complex)(self.snr_info["data"][:,0], self.snr_info["data"][:,1])
self.snr_info["data"].shape = len(self.snr_info["data"]), 1
@staticmethod
def make_template(template_table, template_psd, sample_rate, approximant, instrument, f_low, f_high = None, autocorrelation_length = None, verbose = False):
row = [row for row in template_table if row.ifo == instrument]
if len(row) != 1 :
raise ValueError("Expecting only one template for instrument=%s or cannot find template for instrument=%s" % (instrument, instrument))
template_psd = lal.series.read_psd_xmldoc(ligolw_utils.load_url(template_psd, contenthandler = lal.series.PSDContentHandler))
if instrument not in set(template_psd):
raise ValueError("No such instrument: %s in template psd: (%s)"% (instrument, ", ".join(set(template_psd))))
# work around for building a single whitened template
template_duration = lalsim.SimInspiralChirpTimeBound(f_low, row[0].mass1 * lal.MSUN_SI, row[0].mass2 * lal.MSUN_SI, 0., 0.)
time_slice = numpy.array([(sample_rate, 0, template_duration)], dtype = [("rate", "int"),("begin", "float"), ("end", "float")])
workspace = cbc_template_fir.templates_workspace(template_table, approximant, template_psd[instrument], f_low, time_slice, autocorrelation_length = None, fhigh = f_high)
template, autocorrelation, sigma = workspace.make_whitened_template(row[0])
return template, row[0].end
def __call__(self, COMPLEX = False, row_number = 0 , start = None, end = None):
return self.get_snr_series(COMPLEX, row_number, start, end)
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment