def peak_luminosity_of_merger(*args, NRfit='average', return_fits_used=False):
    """Return the peak luminosity of an aligned-spin BBH using NR fits

    Parameters
    ----------
    mass_1: float/np.ndarray
        float/array of masses for the primary object
    mass_2: float/np.ndarray
        float/array of masses for the secondary object
    spin_1z: float/np.ndarray
        float/array of primary spin aligned with the orbital angular momentum
    spin_2z: float/np.ndarray
        float/array of secondary spin aligned with the orbital angular momentum
    NRFit: str
        Name of the NR fit you wish to use if chosen method is NR
    return_fits_used: Bool, optional
        if True, return the NR fits that were used. Only used when
        NRFit='average'
    """
    from pesummary.gw.file import nrutils
    if NRfit.lower() == 'average':
        func = getattr(nrutils, 'bbh_peak_luminosity_average')
    else:
        func = getattr(nrutils, 'bbh_peak_luminosity_non_precessing_{}'.
            format(NRfit))
    if NRfit.lower() == 'average':
        return func(*args, return_fits_used=return_fits_used)
    return func(*args)
