Skip to content
Snippets Groups Projects
Commit d1b410b5 authored by Shio Sakon's avatar Shio Sakon
Browse files

updated gstlal-inspiral/bin/gstlal_inspiral_mass_model with Tom's mass model

parent 92fdb089
No related branches found
No related tags found
2 merge requests!159updated gstlal-inspiral/bin/gstlal_inspiral_mass_model with "salpeter_primary",!157updated gstlal-inspiral/bin/gstlal_inspiral_mass_model with Tom's mass model
Pipeline #356763 passed with warnings
......@@ -47,7 +47,7 @@ parser.add_argument("--template-bank", metavar='name', type=str, help='The input
parser.add_argument("--reference-psd", metavar='filename', help = "Load the spectrum from this LIGO light-weight XML file")
parser.add_argument("--output", metavar='name', type=str, help='The output file name', default = "inspiral_mass_model.h5")
parser.add_argument("--min-instrument", metavar='name', type=int, default=2, help='Minimum instruments operating. Specified for calculating horizon distance in uniform in template model.')
parser.add_argument("--model", metavar='name', type=str, help='Mass model. Options are: ligo, narrow-bns, ligo-bns, bbh, uniform-template. If you want another one, submit a patch.')
parser.add_argument("--model", metavar='name', type=str, help='Mass model. Options are: ligo, narrow-bns, ligo-bns, bbh, uniform-template, tom. If you want another one, submit a patch.')
parser.add_argument("--verbose", help='Be verbose', action="store_true")
options = parser.parse_args()
......@@ -84,6 +84,12 @@ bbh_peak = 100
bbh_alpha = -1.5
bbhnorm = schechter_norm(bbh_min, bbh_max, bbh_peak, bbh_alpha)
tom_min = 1
tom_max = 600
tom_peak = 600
tom_alpha = -2.35
tom_norm = schechter_norm(tom_min, tom_max, tom_peak, tom_alpha)
for row in sngl_inspiral_table:
assert row.template_id not in prob
......@@ -147,6 +153,16 @@ for row in sngl_inspiral_table:
prob[row.template_id] = (bns_rate * bnsprob + bbh_rate * bbhprob) / (bns_rate + bbh_rate)
# Divide out the template density
prob[row.template_id] /= TEMPDENS(mchirp)
elif options.model == "tom":
#
# For Tom's mass moodel found here: https://wiki.ligo.org/CBC/RatesPop/O4PastroPlan
#
prob[row.template_id] = schechter(row.mass1, tom_peak, tom_alpha) / tom_norm
# Divide out the template density
prob[row.template_id] /= TEMPDENS(mchirp)
else:
raise ValueError("Invalid mass model")
......
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