diff --git a/ligo/em_bright/computeDiskMass.py b/ligo/em_bright/computeDiskMass.py index 26e14669779501d62505bf69d51e6b695824f2ca..4e0545d5ae70ac3e8a6e3ee7e7e8bf663ec56fd2 100644 --- a/ligo/em_bright/computeDiskMass.py +++ b/ligo/em_bright/computeDiskMass.py @@ -64,8 +64,17 @@ def _compactness_baryon_mass(m_ns, r_ns): return [C_ns, m2_b] -def computeCompactness(M_ns, eosname='2H', R_ns=None, - max_mass=2.834648092299807): +def max_mass_from_eosname(eosname): + if eosname == "2H": + max_mass = 2.834648092299807 + else: + eos = lalsim.SimNeutronStarEOSByName(eosname) + fam = lalsim.CreateSimNeutronStarFamily(eos) + max_mass = lalsim.SimNeutronStarMaximumMass(fam)/c.M_sun.value + return max_mass + + +def computeCompactness(M_ns, eosname='2H', max_mass=None): ''' Return the neutron star compactness as a function of mass and equation of state or radius @@ -76,8 +85,6 @@ def computeCompactness(M_ns, eosname='2H', R_ns=None, Neutron star mass in solar masses eosname : str or interp1d Neutron star equation of state to be used - R_ns : array_like - Neutron star radius in m. max_mass : float Maximum mass of neutron star. @@ -91,7 +98,7 @@ def computeCompactness(M_ns, eosname='2H', R_ns=None, ----- The radius and maximum mass of the neutron star is inferred based on the equation of state supplied. - Supply explicitly if `eosname=None`. + Max mass only needs to be supplied for EoS marginalization. Examples -------- @@ -102,12 +109,8 @@ def computeCompactness(M_ns, eosname='2H', R_ns=None, >>> m_ns = np.array([1.1, 1.2, 1.3]) >>> computeDiskMass.computeCompactness(m_ns, eosname='AP4') [array([0.141, 0.154, 0.167]), array([1.199, 1.318, 1.439]), 2.212] - >>> computeCompactness(2.1, R_ns=16000, eosname=None, max_mass=2.5) - [0.193, 2.362, 2.5] ''' - if R_ns: - C_ns, m2_b = _compactness_baryon_mass(M_ns, R_ns) - elif isinstance(eosname, interp1d): + if isinstance(eosname, interp1d): # find R as a function of M R_ns = eosname(M_ns) C_ns, m2_b = _compactness_baryon_mass(M_ns, R_ns) @@ -118,7 +121,7 @@ def computeCompactness(M_ns, eosname='2H', R_ns=None, else: with open(PACKAGE_FILENAMES['equil_2H.dat'], 'rb') as f: M_g, M_b, Compactness = np.loadtxt(f, unpack=True) - max_mass = 2.834648092299807 + max_mass = max_mass_from_eosname("2H") s = UnivariateSpline(M_g, Compactness, k=5) s_b = UnivariateSpline(M_g, M_b, k=5) C_ns = s(M_ns) @@ -244,7 +247,7 @@ def computeDiskMass(m1, m2, chi1, chi2, eosname='2H', kerr=False, chi_secondary = chi1 [C_ns, m2_b, max_mass] = computeCompactness(m_secondary, eosname=eosname, - R_ns=R_ns, max_mass=max_mass) + max_mass=max_mass) eta = m_primary*m_secondary/(m_primary + m_secondary)**2 BBH = m_secondary > max_mass diff --git a/ligo/em_bright/em_bright.py b/ligo/em_bright/em_bright.py index 835acd4031fcc538e7b3d8a476d9e8ed350f86ea..475f1079505d012ed4b05fe177d9eebf504894e5 100644 --- a/ligo/em_bright/em_bright.py +++ b/ligo/em_bright/em_bright.py @@ -157,9 +157,6 @@ def source_classification_pe(posterior_samples_file, **kwargs): posterior_samples_file : str Posterior samples file - threshold : float, optional - Maximum neutron star mass for `HasNS` computation - num_eos_draws : int providing an int here runs eos marginalization with the value determining how many eos's to draw @@ -198,9 +195,6 @@ def source_classification_pe_from_table(table, **kwargs): table : numpy.recarray, dict table containing the posterior samples - threshold : float, optional - Maximum neutron star mass for `HasNS` computation - num_eos_draws : int providing an int here runs eos marginalization with the value determining how many eos's to draw @@ -245,9 +239,8 @@ def source_classification_pe_from_table(table, **kwargs): def source_classification_pe_from_samples(mass_1_source, mass_2_source, - spin_1z, spin_2z, threshold=3.0, - num_eos_draws=10000, eos_seed=None, - eosname=None): + spin_1z, spin_2z, eosname=None, + num_eos_draws=10000, eos_seed=None): """ Compute ``HasNS``, ``HasRemnant``, and ``HasMassGap`` probabilities from samples. @@ -268,17 +261,14 @@ def source_classification_pe_from_samples(mass_1_source, mass_2_source, Samples for the spin component aligned with the orbital angular momentum for the secondary object - threshold : float, optional - Maximum neutron star mass for `HasNS` computation - - num_eos_draws : int + num_eos_draws : int, optional providing an int here runs eos marginalization with the value determining how many eos's to draw - eos_seed : int + eos_seed : int, optional seed for random eos draws - eosname : str + eosname : str, optional Equation of state name, inferred from ``lalsimulation``. Supersedes eos marginalization method when provided. @@ -288,10 +278,12 @@ def source_classification_pe_from_samples(mass_1_source, mass_2_source, (HasNS, HasRemnant, HasMassGap) predicted values. """ if eosname: - M_rem = computeDiskMass.computeDiskMass(mass_1_source, mass_2_source, + M_rem = computeDiskMass.computeDiskMass(mass_1_source, + mass_2_source, spin_1z, spin_2z, eosname=eosname) - prediction_ns = np.sum(mass_2_source <= threshold)/len(mass_2_source) + max_mass = computeDiskMass.max_mass_from_eosname(eosname) + prediction_ns = np.sum(mass_2_source <= max_mass)/len(mass_2_source) prediction_em = np.sum(M_rem > 0)/len(M_rem) else: diff --git a/ligo/em_bright/tests/test_em_bright.py b/ligo/em_bright/tests/test_em_bright.py index bdbba3212f383727350f3dae7ee86233d7f77603..cc375d4b53e83fba29eeec34fb446e2c28904a96 100644 --- a/ligo/em_bright/tests/test_em_bright.py +++ b/ligo/em_bright/tests/test_em_bright.py @@ -19,62 +19,71 @@ def test_version(): @pytest.mark.parametrize( - 'posteriors, dtype, result, result_eos', + 'posteriors, dtype, result_2H, result_SLy, result_eos', [[[(1.2, 1.0, 0.0, 0.0, 0.0, 0.0, 100.0), (2.0, 0.5, 0.99, 0.99, 0.0, 0.0, 150.0)], [('chirp_mass', '<f8'), ('mass_ratio', '<f8'), ('a_1', '<f8'), ('a_2', '<f8'), ('tilt_1', '<f8'), ('tilt_2', '<f8'), ('luminosity_distance', '<f8')], - (1.0, 1.0, 0.5), (1.0, 1.0, 0.5)], + (1.0, 1.0, 0.5), (1.0, 1.0, 0.5), (1.0, 1.0, 0.5)], [[(1.2, 1.0, 0.0, 0.0, 100.0), (2.0, 0.5, 0.99, 0.99, 150.0)], [('chirp_mass', '<f8'), ('mass_ratio', '<f8'), ('a_1', '<f8'), ('a_2', '<f8'), ('luminosity_distance', '<f8')], - (1.0, 1.0, 0.5), (1.0, 0.5, 0.5)], + (1.0, 1.0, 0.5), (1.0, 0.5, 0.5), (1.0, 0.529, 0.5)], [[(1.4, 1.4, 0.0, 0.0, 100.0), (2.0, 0.5, 0.99, 0.99, 150.0)], [('mass_1', '<f8'), ('mass_2', '<f8'), ('a_1', '<f8'), ('a_2', '<f8'), ('luminosity_distance', '<f8')], - (1.0, 1.0, 0.0), (1.0, 1.0, 0.0)], + (1.0, 1.0, 0.0), (1.0, 1.0, 0.0), (1.0, 1.0, 0.0)], [[(1.4, 1.4, 100.0), (2.0, 0.5, 150.0)], [('mass_1', '<f8'), ('mass_2', '<f8'), ('luminosity_distance', '<f8')], - (1.0, 1.0, 0.0), (1.0, 1.0, 0.0)], + (1.0, 1.0, 0.0), (1.0, 1.0, 0.0), (1.0, 1.0, 0.0)], [[(1.4, 1.4, 1.4, 1.4, 0.0, 0.0, 100.0), (2.0, 0.5, 2.0, 0.5, 0.99, 0.99, 150.0)], [('mass_1_source', '<f8'), ('mass_2_source', '<f8'), ('mass_1', '<f8'), ('mass_2', '<f8'), ('a_1', '<f8'), ('a_2', '<f8'), ('luminosity_distance', '<f8')], - (1.0, 1.0, 0.0), (1.0, 1.0, 0.0)], + (1.0, 1.0, 0.0), (1.0, 1.0, 0.0), (1.0, 0.99985, 0.0)], [[(4.5, -0.1, 200.0, 100000, 1.4, 1.4), (1.6, 0.3, 201.0, 100000, 1.5, 1.3)], [('ra', '<f8'), ('dec', '<f8'), ('luminosity_distance', '<f8'), ('time', '<f8'), ('mass_1', '<f8'), ('mass_2', '<f8')], - (1.0, 1.0, 0.0), (1.0, 1.0, 0.0)], + (1.0, 1.0, 0.0), (1.0, 1.0, 0.0), (1.0, 1.0, 0.0)], [[(4.5, -0.1, 200.0, 100000, 4.5, 4.4), (1.6, 0.3, 201.0, 100000, 4.3, 4.2)], [('ra', '<f8'), ('dec', '<f8'), ('luminosity_distance', '<f8'), ('time', '<f8'), ('mass_1', '<f8'), ('mass_2', '<f8')], - (0.0, 0.0, 1.0), (0.0, 0.0, 1.0)], + (0.0, 0.0, 1.0), (0.0, 0.0, 1.0), (0.0, 0.0, 1.0)], [[(4.5, -0.1, 200.0, 100000, 40.5, 4.4), (1.6, 0.3, 201.0, 100000, 40.3, 4.2)], [('ra', '<f8'), ('dec', '<f8'), ('luminosity_distance', '<f8'), ('time', '<f8'), ('mass_1', '<f8'), ('mass_2', '<f8')], - (0.0, 0.0, 1.0), (0.0, 0.0, 1.0)], + (0.0, 0.0, 1.0), (0.0, 0.0, 1.0), (0.0, 0.0, 1.0)], [[(4.5, -0.1, 200.0, 1.0, 2.0, 2.0, 0.0, 0.0), (1.6, 0.3, 201.0, 1.0, 5.0, 2.0, 0.0, 0.0)], [('ra', '<f8'), ('dec', '<f8'), ('luminosity_distance', '<f8'), ('time', '<f8'), ('mass_1', '<f8'), ('mass_2', '<f8'), ('spin_1z', '<f8'), ('spin_2z', '<f8')], - (1.0, 0.5, 0.5), (1.0, 0.5, 0.5)], + (1.0, 0.5, 0.5), (1.0, 0.5, 0.5), (1.0, 0.50045, 0.5)], [[(4.5, -0.1, 200.0, 1.0, 2.0, 2.0, 0.99, 0.0), (1.6, 0.3, 201.0, 1.0, 5.0, 2.0, 0.99, 0.0)], [('ra', '<f8'), ('dec', '<f8'), ('luminosity_distance', '<f8'), ('time', '<f8'), ('mass_1', '<f8'), ('mass_2', '<f8'), ('spin_1z', '<f8'), ('spin_2z', '<f8')], - (1.0, 1.0, 0.5), (1.0, 1.0, 0.5)]] + (1.0, 1.0, 0.5), (1.0, 1.0, 0.5), (1.0, 1.0, 0.5)], + [[(3.0, 2.8, 0.0, 0.0, 100.0)], + [('mass_1_source', '<f8'), ('mass_2_source', '<f8'), + ('a_1', '<f8'), ('a_2', '<f8'), ('luminosity_distance', '<f8')], + (1.0, 0.0, 1.0), (0.0, 0.0, 1.0), (0.0034, 0.0, 1.0)], + [[(2.5, 1.8, 0.0, 0.0, 100.0)], + [('mass_1_source', '<f8'), ('mass_2_source', '<f8'), + ('a_1', '<f8'), ('a_2', '<f8'), ('luminosity_distance', '<f8')], + (1.0, 1.0, 0.0), (1.0, 0.0, 0.0), (1.0, 0.0426, 0.0)]] ) -def test_source_classification_pe(posteriors, dtype, result, result_eos): +def test_source_classification_pe(posteriors, dtype, result_2H, + result_SLy, result_eos): """Test em_bright classification from posterior samples - both aligned and precessing cases. """ @@ -89,10 +98,13 @@ def test_source_classification_pe(posteriors, dtype, result, result_eos): 'posterior_samples', data=data ) - r = em_bright.source_classification_pe(filename, eosname='2H') - r_eos = em_bright.source_classification_pe(filename, num_eos_draws=5, + r_2H = em_bright.source_classification_pe(filename, eosname='2H') + r_SLy = em_bright.source_classification_pe(filename, eosname='SLy') + r_eos = em_bright.source_classification_pe(filename, + num_eos_draws=10000, eos_seed=0) - assert r == result + assert r_2H == result_2H + assert r_SLy == result_SLy assert r_eos == result_eos