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