diff --git a/gstlal-inspiral/bin/gstlal_inspiral_mass_model b/gstlal-inspiral/bin/gstlal_inspiral_mass_model index 1ecb18a9fbda24ed5ce9ba60db8d5a45ba3ae830..661099a1effb9669806c7c939f3b2eacb6f5df6f 100755 --- a/gstlal-inspiral/bin/gstlal_inspiral_mass_model +++ b/gstlal-inspiral/bin/gstlal_inspiral_mass_model @@ -32,28 +32,34 @@ from lal import rate class LIGOLWContentHandler(ligolw.LIGOLWContentHandler): pass +def chirpmass(m1, m2): + return (m1 * m2)**.6 / (m1 + m2)**.2 + parser = argparse.ArgumentParser(description = "Create analytic mass models for prior weighting of templates") parser.add_argument("--template-bank", metavar='name', type=str, help='The input template bank file name.', required = True) parser.add_argument("--output", metavar='name', type=str, help='The output file name', default = "inspiral_mass_model.h5") -parser.add_argument("--model", metavar='name', type=str, help='Mass model. Options are: salpeter. If you want another one, submit a patch.') +parser.add_argument("--model", metavar='name', type=str, help='Mass model. Options are: salpeter, ligo. If you want another one, submit a patch.') parser.add_argument("--verbose", help='Be verbose', action="store_true") options = parser.parse_args() # Read the template bank file xmldoc = ligolw_utils.load_filename(options.template_bank, verbose = options.verbose, contenthandler = LIGOLWContentHandler) sngl_inspiral_table = lsctables.SnglInspiralTable.get_table(xmldoc) -mass1 = sngl_inspiral_table.get_column("mass1") -mass2 = sngl_inspiral_table.get_column("mass2") -num_templates = len(mass1) -num_bins = max(2, int((num_templates / 100.)**.5)) -min_mass = min(min(mass1), min(mass2)) - 1.e-6 -max_mass = max(max(mass1), max(mass2)) + 1.e-6 -massBA = rate.BinnedDensity(rate.NDBins((rate.LogarithmicBins(min_mass, max_mass, num_bins), rate.LogarithmicBins(min_mass, max_mass, num_bins)))) -print min_mass, max_mass -for m1, m2 in zip(mass1, mass2): - massBA.count[(m1, m2)] += 1 - massBA.count[(m2, m1)] += 1 -rate.filter_array(massBA.array, rate.gaussian_window(1.5, 1.5, sigma = 5)) + +# +# Someday if noise is actually a pdf in mass this might matter +# +# mass1 = sngl_inspiral_table.get_column("mass1") +# mass2 = sngl_inspiral_table.get_column("mass2") +# num_templates = len(mass1) +# num_bins = max(2, int((num_templates / 100.)**.5)) +# min_mass = min(min(mass1), min(mass2)) - 1.e-6 +# max_mass = max(max(mass1), max(mass2)) + 1.e-6 +# massBA = rate.BinnedDensity(rate.NDBins((rate.LogarithmicBins(min_mass, max_mass, num_bins), rate.LogarithmicBins(min_mass, max_mass, num_bins)))) +# for m1, m2 in zip(mass1, mass2): +# massBA.count[(m1, m2)] += 1 +# massBA.count[(m2, m1)] += 1 +# rate.filter_array(massBA.array, rate.gaussian_window(1.5, 1.5, sigma = 5)) # Assign the proper mass probabilities ids = {} @@ -61,14 +67,26 @@ tmplt_ids = [] for row in sngl_inspiral_table: assert row.template_id not in ids tmplt_ids.append(int(row.template_id)) + primary = max(row.mass1, row.mass2) if options.model == "salpeter": - ids[row.template_id] = numpy.log(row.mass1**-2.35 / massBA[row.mass1, row.mass2]) + ids[row.template_id] = primary**-2.35# / massBA[row.mass1, row.mass2] + elif options.model == "ligo": + # assume a 0.15 solar mass std deviation, this should capture both population distribution and snr effects + sigma = 0.15 + mean = 1.2 + bnsprob = 1. / (2 * numpy.pi * sigma**2)**.5 * numpy.exp(-(chirpmass(row.mass1, row.mass2) - mean)**2 / 2. / sigma**2) + # normalised over 5 -- 45 Msun + bbhprob = 0.46 * primary**-1.6 + # From: https://www.lsc-group.phys.uwm.edu/ligovirgo/cbcnote/RatesAndSignificance/O1O2CatalogRates + bns_to_bbh_rate = 916. / 56. + ids[row.template_id] = (bns_to_bbh_rate * bnsprob + bbhprob)# / massBA[row.mass1, row.mass2] else: raise ValueError("Invalid mass model") +norm = sum(ids.values()) coefficients = numpy.zeros((1, 1, max(ids)+1), dtype=float) for tid in ids: - coefficients[0,0,tid] = ids[tid] + coefficients[0,0,tid] = numpy.log(ids[tid] / norm * len(sngl_inspiral_table)) # Write it out f = h5py.File(options.output, "w")