From 8c3ce90b6cc70cfd17b19a05050980659588ecda Mon Sep 17 00:00:00 2001
From: Leo Tsukada <leo.tsukada@ligo.org>
Date: Mon, 21 Mar 2022 13:45:22 -0700
Subject: [PATCH 1/2] gstlal_inspiral_bank_splitter : implement 2D-sorting by
 orthogonalized 0/1PN order phase terms

---
 .../bin/gstlal_inspiral_bank_splitter         | 46 +++++++++++++++++--
 1 file changed, 42 insertions(+), 4 deletions(-)

diff --git a/gstlal-inspiral/bin/gstlal_inspiral_bank_splitter b/gstlal-inspiral/bin/gstlal_inspiral_bank_splitter
index b32e53b74f..023ac238f5 100755
--- a/gstlal-inspiral/bin/gstlal_inspiral_bank_splitter
+++ b/gstlal-inspiral/bin/gstlal_inspiral_bank_splitter
@@ -31,6 +31,7 @@ from ligo.lw import utils as ligolw_utils
 from ligo.lw.utils import process as ligolw_process
 import lal
 import numpy, scipy.interpolate
+from astropy import constants as const
 
 from gstlal import chirptime
 from gstlal import datafind
@@ -71,6 +72,33 @@ class LIGOLWContentHandler(ligolw.LIGOLWContentHandler):
 	pass
 lsctables.use_in(LIGOLWContentHandler)
 
+def calc_mu(mass1, mass2, spin1z, spin2z, fref=100, mu="mu1"):
+	"""
+	Calculate the first orthogonal PN phase coefficient
+	see https://arxiv.org/abs/2007.09108
+	"""
+
+	M = mass1 + mass2
+	mchirp = (mass1 * mass2)**0.6/M**0.2
+	eta = mass1 * mass2 / M**2
+	beta = ((113. * (mass1 / M)**2 + 75. * eta) * spin1z + (113. * (mass2 / M)**2 + 75. * eta) * spin2z) / 12.
+
+	norm = const.G.value * const.M_sun.value / const.c.value**3
+	v = numpy.pi * mchirp * fref * norm
+	psi0 = 3. / 4 / (8 * v)**(5./3)
+	psi2 = 20. / 9 * (743./ 336 + 11. / 4 * eta) / eta**(0.4) * v**(2./3) * psi0
+	psi3 = (4 * beta - 16 * numpy.pi) / eta**0.6 * v * psi0
+
+	# FIXME : the following linear combinations are taken from the ones in the
+	# paper above, but this will need to be re-computed with o4 representitive
+	# psd.
+	if mu=="mu1":
+		return 0.974 * psi0 + 0.209 * psi2 + 0.0840 * psi3
+	elif mu=="mu2":
+		return -0.221 * psi0 + 0.823 * psi2 + 0.524 * psi3
+	else:
+		raise ValueError("%s is not implemented, so cannnot be computed." % mu)
+
 
 def group_templates(templates, n, overlap = 0):
 	"""
@@ -191,12 +219,19 @@ if options.output_full_bank_file is not None:
 
 outputrows = []
 
-# Bin by Chi
-sngl_inspiral_table.sort(key = lambda row: spawaveform.computechi(row.mass1, row.mass2, row.spin1z, row.spin2z))
+if options.sort_by == "mu":
+	# First partitioning by mu1
+	sngl_inspiral_table.sort(key = lambda row: calc_mu(row.mass1, row.mass2, row.spin1z, row.spin2z, mu="mu1"))
+else:
+	# First partitioning by chi
+	sngl_inspiral_table.sort(key = lambda row: spawaveform.computechi(row.mass1, row.mass2, row.spin1z, row.spin2z))
 for chirows in group_templates(sngl_inspiral_table, len(sngl_inspiral_table) / options.group_by_chi, overlap = 0):
 
 	def sort_func(row, column = options.sort_by):
-		return getattr(row, column)
+		if column=="mu":
+			return calc_mu(row.mass1, row.mass2, row.spin1z, row.spin2z, mu="mu2")
+		else:
+			return getattr(row, column)
 
 	chirows.sort(key=sort_func, reverse = True if options.sort_by == "template_duration" else False)
 
@@ -211,7 +246,10 @@ for chirows in group_templates(sngl_inspiral_table, len(sngl_inspiral_table) / o
 
 # A sort of the groups of templates so that the sub banks are ordered.
 def sort_func(row_rows_tup, column = options.sort_by):
-	return getattr(row_rows_tup[0], column)
+	if column=="mu":
+		return calc_mu(row_rows_tup[0].mass1, row_rows_tup[0].mass2, row_rows_tup[0].spin1z, row_rows_tup[0].spin2z, mu="mu2")
+	else:
+		return getattr(row_rows_tup[0], column)
 
 outputrows.sort(key=sort_func, reverse = True if options.sort_by == "template_duration" else False)
 svd_dirs = []
-- 
GitLab


From 3e15b66a8ad25e9cc3d75f6db694a1c5bbbce055 Mon Sep 17 00:00:00 2001
From: Leo Tsukada <leo.tsukada@ligo.org>
Date: Mon, 21 Mar 2022 15:16:29 -0700
Subject: [PATCH 2/2] gstlal_inspiral_bank_splitter : make lal constants and
 add --group-by-mu for consistentcy

---
 gstlal-inspiral/bin/gstlal_inspiral_bank_splitter | 8 +++++---
 1 file changed, 5 insertions(+), 3 deletions(-)

diff --git a/gstlal-inspiral/bin/gstlal_inspiral_bank_splitter b/gstlal-inspiral/bin/gstlal_inspiral_bank_splitter
index 023ac238f5..0f4b4f0637 100755
--- a/gstlal-inspiral/bin/gstlal_inspiral_bank_splitter
+++ b/gstlal-inspiral/bin/gstlal_inspiral_bank_splitter
@@ -31,7 +31,6 @@ from ligo.lw import utils as ligolw_utils
 from ligo.lw.utils import process as ligolw_process
 import lal
 import numpy, scipy.interpolate
-from astropy import constants as const
 
 from gstlal import chirptime
 from gstlal import datafind
@@ -83,7 +82,7 @@ def calc_mu(mass1, mass2, spin1z, spin2z, fref=100, mu="mu1"):
 	eta = mass1 * mass2 / M**2
 	beta = ((113. * (mass1 / M)**2 + 75. * eta) * spin1z + (113. * (mass2 / M)**2 + 75. * eta) * spin2z) / 12.
 
-	norm = const.G.value * const.M_sun.value / const.c.value**3
+	norm = lal.G_SI * lal.MSUN_SI / lal.C_SI**3
 	v = numpy.pi * mchirp * fref * norm
 	psi0 = 3. / 4 / (8 * v)**(5./3)
 	psi2 = 20. / 9 * (743./ 336 + 11. / 4 * eta) / eta**(0.4) * v**(2./3) * psi0
@@ -134,6 +133,7 @@ def parse_command_line():
 	parser.add_option("--approximant", type = "string", action = "append", help = "Must specify an approximant given as mchirp_min:mchirp_max:string")
 	parser.add_option("--f-low", type = "float", metavar = "frequency", help = "Lower frequency cutoff. Required")
 	parser.add_option("--group-by-chi", type = "int", metavar = "N", default = 1, help = "group templates into N groups of chi - helps with SVD. Default 1")
+	parser.add_option("--group-by-mu", type = "int", metavar = "N", default = 20, help = "group templates into N groups of mu1, the 1PN phase term, to help with SVD. Default 20")
 	parser.add_option("--num-banks", metavar = "str", help = "The number of parallel subbanks per SVD bank. can be given as a list like 1,2,3,4 then it will split up the bank into N groups with M banks each. (required)") # FIXME The second half of this help message is incomprehensible
 	options, filenames = parser.parse_args()
 
@@ -222,10 +222,12 @@ outputrows = []
 if options.sort_by == "mu":
 	# First partitioning by mu1
 	sngl_inspiral_table.sort(key = lambda row: calc_mu(row.mass1, row.mass2, row.spin1z, row.spin2z, mu="mu1"))
+	num_group = options.group_by_mu
 else:
 	# First partitioning by chi
 	sngl_inspiral_table.sort(key = lambda row: spawaveform.computechi(row.mass1, row.mass2, row.spin1z, row.spin2z))
-for chirows in group_templates(sngl_inspiral_table, len(sngl_inspiral_table) / options.group_by_chi, overlap = 0):
+	num_group = options.group_by_chi
+for chirows in group_templates(sngl_inspiral_table, len(sngl_inspiral_table) / num_group, overlap = 0):
 
 	def sort_func(row, column = options.sort_by):
 		if column=="mu":
-- 
GitLab