def get_random_mass(numPoints, massRangeParams):
    """
    This function will generate a large set of points within the chosen mass
    and spin space, and with the desired minimum remnant disk mass (this applies
    to NS-BH systems only). It will also return the corresponding PN spin
    coefficients for ease of use later (though these may be removed at some
    future point).

    Parameters
    ----------
    numPoints : int
        Number of systems to simulate
    massRangeParams : massRangeParameters instance
        Instance holding all the details of mass ranges and spin ranges.

    Returns
    --------
    mass1 : float
        Mass of heavier body.
    mass2 : float
        Mass of lighter body.
    spin1z : float
        Spin of body 1.
    spin2z : float
        Spin of body 2.
    """

    # WARNING: We expect mass1 > mass2 ALWAYS

    # Check if EM contraints are required, i.e. if the systems must produce
    # a minimum remnant disk mass.  If this is not the case, proceed treating
    # the systems as point particle binaries
    if massRangeParams.remnant_mass_threshold is None:
        mass1, mass2, spin1z, spin2z = \
            get_random_mass_point_particles(numPoints, massRangeParams)
    # otherwise, load EOS dependent data, generate the EM constraint
    # (i.e. compute the minimum symmetric mass ratio needed to
    # generate a given remnant disk mass as a function of the NS
    # mass and the BH spin along z) and then proceed by accepting
    # only systems that can yield (at least) the desired remnant
    # disk mass and that pass the mass and spin range cuts.
    else:
        _, max_ns_g_mass = load_ns_sequence(massRangeParams.ns_eos)

        # Generate EM constraint surface: minumum eta as a function of BH spin
        # and NS mass required to produce an EM counterpart
        if not os.path.isfile('constraint_em_bright.npz'):
            logging.info("""constraint_em_bright.npz not found.
                        Generating the constraint surface for EM bright binaries
                        in the physical parameter space.  One day, this will be
                        made faster, for now be patient and wait a few minutes!""")
            generate_em_constraint_data(massRangeParams.minMass2, massRangeParams.maxMass2, massRangeParams.delta_ns_mass, \
                                        -1.0, massRangeParams.maxBHSpinMag, massRangeParams.delta_bh_spin, \
                                        massRangeParams.ns_eos, massRangeParams.remnant_mass_threshold, 0.0)
        constraint_datafile = numpy.load('constraint_em_bright.npz')
        mNS_pts = constraint_datafile['mNS_pts']
        bh_spin_z_pts = constraint_datafile['sBH_pts']
        eta_mins = constraint_datafile['eta_mins']

        # Empty arrays to store points that pass all cuts
        mass1_out = []
        mass2_out = []
        spin1z_out = []
        spin2z_out = []

        # As the EM cut can remove several randomly generated
        # binaries, track the number of accepted points that pass
        # all cuts and stop only once enough of them are generated
        numPointsFound = 0
        while numPointsFound < numPoints:
            # Generate the random points within the required mass
            # and spin cuts
            mass1, mass2, spin1z, spin2z = \
                get_random_mass_point_particles(numPoints-numPointsFound,
                                                massRangeParams)
            _, eta = pnutils.mass1_mass2_to_mtotal_eta(mass1, mass2)

            # Now proceed with cutting out EM dim systems
            # Logical mask to clean up points by removing EM dim binaries
            mask = numpy.ones(len(mass1), dtype=bool)
            # Commpute the minimum eta to generate a counterpart
            min_eta_em = min_eta_for_em_bright(spin1z, mass2, mNS_pts, bh_spin_z_pts, eta_mins)
            # Remove a point if:
            # 1) eta is smaller than the eta threshold required to have a counterpart;
            # 2) the primary is a BH (mass1 >= ns_bh_boundary_mass);
            # 3) the secondary mass does not exceed the maximum NS mass
            # allowed by the EOS (if the user runs with --use-eos-max-ns-mass
            # this last condition will always be true, otherwise the user is
            # implicitly asking to keep binaries in which the secondary may be
            # a BH).
            mask[(mass1 >= massRangeParams.ns_bh_boundary_mass) & (mass2 <= max_ns_g_mass) & (eta < min_eta_em)] = False
            # Keep only binaries that can produce an EM counterpart and add them to
            # the pile of accpeted points to output
            mass1_out = numpy.concatenate((mass1_out, mass1[mask]))
            mass2_out = numpy.concatenate((mass2_out, mass2[mask]))
            spin1z_out = numpy.concatenate((spin1z_out,spin1z[mask]))
            spin2z_out = numpy.concatenate((spin2z_out,spin2z[mask]))

            # Number of points that survived all cuts
            numPointsFound = len(mass1_out)

        # Ready to go
        mass1 = mass1_out
        mass2 = mass2_out
        spin1z = spin1z_out
        spin2z = spin2z_out

    return mass1, mass2, spin1z, spin2z