From ba0d18b6bcfd98aa262c699c1893f638de7bbb73 Mon Sep 17 00:00:00 2001 From: Anonymous Date: Tue, 20 Mar 2018 18:57:49 +0000 Subject: [PATCH] Added more work towards ringdowns, and moving xml reading to glue rather than weird bits of lalsuite. --- minke/mdctools.py | 64 ++++++++++++++++++++++++++--------------- minke/sources.py | 29 +++++++++++++++++-- tests/test_minke.py | 69 ++++++++++++++++++++++++++++++++++++++++++++- 3 files changed, 137 insertions(+), 25 deletions(-) diff --git a/minke/mdctools.py b/minke/mdctools.py index d51c6a8..320dd3c 100644 --- a/minke/mdctools.py +++ b/minke/mdctools.py @@ -31,6 +31,9 @@ from pylal.xlal.datatypes.ligotimegps import LIGOTimeGPS from glue.ligolw.utils import process import glue +import glue.ligolw +import gzip + import lal, lalframe from pylal import Fr @@ -42,6 +45,7 @@ import matplotlib matplotlib.use('Agg') import matplotlib.pyplot as plt import re +import random table_types = { @@ -57,7 +61,8 @@ table_types = { # Long Duration "adi" : lsctables.SimBurstTable, # Ringdown - "rng" : lsctables.SimRingdownTable + "rng" : lsctables.SimRingdownTable, + "gng" : lsctables.SimRingdownTable, } tables = { @@ -99,7 +104,8 @@ class MDCSet(): # Long-duration 'adi' : 'ADI', # Ringdown - 'rng' : "Ringdown", + 'rng' : "BBHRingdown", + 'gng' : "Ringdown", } inj_families_abb = dict((v,k) for k,v in inj_families_names.iteritems()) @@ -182,17 +188,19 @@ class MDCSet(): waveform_row = waveform._row(sim) waveform_row.process_id = procrow.process_id except: - row = sim.RowType() - - for a in lsctables.SimBurstTable.validcolumns.keys(): - setattr(row, a, getattr(waveform, a)) + for a in self.table_type.validcolumns.keys(): + if not hasattr(waveform, a): + setattr(row, a, 0) + else: + setattr(row, a, getattr(waveform, a)) row.waveform = waveform.waveform - # Fill in the time - row.set_time_geocent(GPS(float(waveform.time))) - # Get the sky locations - row.ra, row.dec, row.psi = waveform.ra, waveform.dec, waveform.psi + if self.table_type == lsctables.SimBurstTable: + # Fill in the time + row.set_time_geocent(GPS(float(waveform.time))) + # Get the sky locations + row.ra, row.dec, row.psi = waveform.ra, waveform.dec, waveform.psi row.simulation_id = waveform.simulation_id row.waveform_number = random.randint(0,int(2**32)-1) ### !! This needs to be updated. @@ -235,24 +243,30 @@ class MDCSet(): file. This should be fixed so that the object works symmetrically. """ i = 0 - sim_burst_table = lalburst.SimBurstTableFromLIGOLw(filename, start, stop) - while True: + #sim_burst_table = lalburst.SimBurstTableFromLIGOLw(filename, start, stop) + + xml = glue.ligolw.utils.load_filename("test_simbursttable.xml.gz", + contenthandler = glue.ligolw.ligolw.LIGOLWContentHandler, + verbose = True) + sim_burst_table = glue.ligolw.table.get_table(xml, self.table_type.tableName) + + for i,simrow in enumerate(sim_burst_table): # This is an ugly kludge to get around the poor choice of wavform name in the xmls, and - if sim_burst_table.waveform[:3]=="s15": + if simrow.waveform[:3]=="s15": self.numrel_file = str(sim_burst_table.waveform) sim_burst_table.waveform = "Dimmelmeier+08" - self.waveforms.append(sim_burst_table) + self.waveforms.append(simrow)#_burst_table) if full: self._measure_hrss(i) self._measure_egw_rsq(i) - self.times = np.append(self.times, float(sim_burst_table.time_geocent_gps)) + self.times = np.append(self.times, float(simrow.time_geocent_gps)) #np.insert(self.times, len(self.times), sim_burst_table.time_geocent_gps) - if sim_burst_table.next is None: break - sim_burst_table = sim_burst_table.next - i += 1 - del(sim_burst_table) + #if sim_burst_table.next is None: break + #sim_burst_table = sim_burst_table.next + #i += 1 + #del(sim_burst_table) def _generate_burst(self,row,rate=16384.0): """ @@ -280,7 +294,7 @@ class MDCSet(): # This is a temporary kludge to allow LALSimulation to # be bypassed for pre-calculated waveforms. # A more robust solution should be considered. - exceptions = ["Ott+13", "Mueller+12", "Scheidegger+10"] + exceptions = ["Ott+13", "Mueller+12", "Scheidegger+10", "ADI"] row = self.waveforms[row] swig_row = lalburst.CreateSimBurst() for a in lsctables.SimBurstTable.validcolumns.keys(): @@ -290,8 +304,14 @@ class MDCSet(): except TypeError: print a, getattr(row,a) continue # the structure is different than the TableRow - theta, phi = np.cos(swig_row.incl), swig_row.phi - swig_row.numrel_data = row.numrel_data + theta, phi = np.cos(swig_row.incl), swig_row.phi + try: + swig_row.numrel_data = row.numrel_data + except AttributeError: + # Don't fail badly if there's a problem with numerical relativity data unless this is a supernova injection + swig_row.numrel_data = "" + if row.waveform in exceptions: + raise AttributeError hp, hx = lalburst.GenerateSimBurst(swig_row, 1.0/rate) hp0, hx0 = lalburst.GenerateSimBurst(swig_row, 1.0/rate) return hp, hx, hp0, hx0 diff --git a/minke/sources.py b/minke/sources.py index c305bc5..8db705a 100644 --- a/minke/sources.py +++ b/minke/sources.py @@ -44,7 +44,8 @@ class Waveform(object): """ sim = lsctables.New(lsctables.SimBurstTable) - + table_type = lsctables.SimBurstTable + numrel_data = [] waveform = "Generic" expnum = 1 @@ -205,7 +206,7 @@ class Waveform(object): if not sim: sim = self.sim row = sim.RowType() - for a in lsctables.SimBurstTable.validcolumns.keys(): + for a in self.table_type.validcolumns.keys(): setattr(row, a, self.params[a]) row.waveform = self.waveform @@ -276,6 +277,7 @@ class SineGaussian(Waveform): self.params['pol_ellipse_e'], self.params['pol_ellipse_angle'] = self.parse_polarisation(self.polarisation) + class Gaussian(Waveform): """ A class to represent a Gaussian injection. @@ -960,6 +962,7 @@ class Ringdown(Waveform): """ A class to handle Rindown waveforms. """ + table_type = lsctables.SimRingdownTable waveform = "Ringdown" pass @@ -1091,3 +1094,25 @@ class ADI(LongDuration): output[:,2] = strainc_new return output + +class BBHRingdown(Ringdown): + """ + A class to represent BBH ringdowns. + """ + #lalsimfunction = SimBlackHoleRingdown + + def __init__(self, time, hrss, phi0, deltaT, mass, spin, massloss, distance, inclination, l, m, sky_dist=uniform_sky): + self._clear_params() + self.time = self.v_start_time_ns = time + self.sky_dist = sky_dist + self.params['simulation_id'] = self.simulation_id = self.sim.get_next_id() + self.params['phi0'] = phi0 + self.params['deltaT'] = deltaT + self.params['mass'] = mass + self.params['spin'] = spin + self.params['massloss'] = massloss + self.params['eff_dist_l'] = self.eff_dist_l = distance + self.params['hrss'] = self.hrss = hrss + self.params['inclination'] = inclination + self.params['l'] = l + self.params['m'] = m diff --git a/tests/test_minke.py b/tests/test_minke.py index 3eec3eb..03ca052 100755 --- a/tests/test_minke.py +++ b/tests/test_minke.py @@ -11,7 +11,7 @@ Tests for `minke` module. import unittest import minke - +from minke import mdctools, sources class TestMinke(unittest.TestCase): @@ -21,6 +21,73 @@ class TestMinke(unittest.TestCase): def tearDown(self): pass + # Check that failure occurs if the wrong waveforms are inserted into the wrong tables + + def test_insert_wrong_waveform_to_table(self): + """Test whether adding a ringdown waveform to a SimBurstTable throws + an error. + """ + mdcset = mdctools.MDCSet(["L1"], table_type = "burst") + ring = sources.Ringdown() + + with self.assertRaises(mdctools.TableTypeError): + mdcset + ring + + def test_insert_burst_waveform_to_burst_table(self): + """ + Test whether inserting the correct type of waveform into a table works + """ + mdcset = mdctools.MDCSet(["L1"], table_type = "burst") + waveform = sources.Gaussian(0,0,0) + + mdcset + waveform + + self.assertEqual(len(mdcset.waveforms), 1) + + # Check that the correct XML files are produced for different table types + + def test_write_simbursttable(self): + """ + Write out a simburst table xml file + """ + mdcset = mdctools.MDCSet(["L1"], table_type = "burst") + waveform = sources.Gaussian(0.1,1e-23,1000) + + mdcset + waveform + + mdcset.save_xml("test_simbursttable.xml.gz") + + def test_verify_simbursttable(self): + """ + Read-in the xml simburst table. + """ + mdcset = mdctools.MDCSet(["L1"]) + mdcset.load_xml("test_simbursttable.xml.gz", full = False) + + self.assertEqual(len(mdcset.waveforms), 1) + + def test_write_simringdowntable(self): + """ + Write out a simburst table xml file + """ + mdcset = mdctools.MDCSet(["L1"], table_type = "ringdown") + waveform = sources.BBHRingdown(1000, 1e-22, 0 ,0.1, 10, 0.1, 0.01, 10, 0, 2, 2,) + + mdcset + waveform + + mdcset.save_xml("test_simringdowntable.xml.gz") + + def test_verify_simringdowntable(self): + """ + Read-in the xml simburst table. + """ + mdcset = mdctools.MDCSet(["L1"]) + mdcset.load_xml("test_simringdowntable.xml.gz", full = False) + + self.assertEqual(len(mdcset.waveforms), 1) + + + if __name__ == '__main__': import sys sys.exit(unittest.main()) -- GitLab