Prospector not recognizing and fitting nebular lines

[gracekrahm]

I am trying to fit SEDs with nebular line emission and only a few of the lines are being with prospector. Even with the ones that are fit, the continuum is much higher than the fsps generated SED that we are fitting. This is especially evident at high redshifts. .

All files used to generate the fsps mock SED and fit it using Prospector can be found here.

[jrleja]

Hi Grace,

I took a look at your repo. It looks like you're fitting photometry and not spectroscopy. Have you checked that you're able to get a good fit to the photometry that you're using as input? It may be that you're achieving a good fit but the underlying model spectra are very different.

It looks like you're performing the filter integration in a for loop in your build_obs function - it may be that there are differences between this (spec ---> phot) procedure and the one used in Prospector (which lives in the sedpy package), so I might also suggest using the filter integration procedures in sedpy for convenience and reliability. I would also suggest using the best-fit stored spectrum in the Prospector results for post-fit evaluation purposes - it is guaranteed to be exactly what was used during the fit (likely not your problem but a good procedure for reliability).

Best, Joel

[gracekrahm]

Hi Joel, Thanks for the quick response. Could you elaborate a little more on what you mean by using the sedpy filter integration procedures? Thanks, Grace

[jrleja]

Hi Grace,

Using sedpy you can do the filter projection in the same way Prospector does, e.g. from the sedpy documentation:

from sedpy import observate
# get magnitude from a spectrum:
filt = observate.Filter("sdss_r0")
mag = filt.ab_mag(angstroms, f_lambda_cgs)
# or get several magnitudes at once
filterlist = observate.load_filters(["galex_NUV", "sdss_r0"])
mags = observate.getSED(angstroms, f_lambda_cgs, filterlist=filters)

In some cases there can be substantive effects in filter projection methods that can affect the inferred magnitudes (e.g. the classic energy-counting vs photon-counting, treatment of interpolation, etc). Best to use the same approach as Prospector in this case.

My shot-in-the-dark guess as to what's going on is that Prospector is matching the photometry to high precision but using a very different model spectrum than the input spectrum. This can either be for trivial reasons (e.g., unit bug, filter projection bug, etc, in which case sedpy might help get everything onto the same scale) or for very interesting science reasons (broadband filter degeneracies).

Best, Joel

[gracekrahm]

Thank you I'll implement this filter integration and see if that helps the fitting process.

[gracekrahm]

Hi Joel, I seem to have made everything worse in my attempt to implement the sedpy integration.

I was wondering if there was a specific point in my build_obs function that is causing this. Thanks, Grace

def build_obs(pd_dir,**kwargs):
    print('loading obs')
    import sedpy
    from astropy import units as u
    from astropy import constants
    from astropy.cosmology import FlatLambdaCDM
    cosmo = FlatLambdaCDM(H0=68, Om0=0.3, Tcmb0=2.725)

    df = pd.read_csv(pd_dir)
    print(df.head())
    wav = df['wav']
    flux = df['spec']
    #wav  = np.asarray(wav)*u.micron #wav is in micron
    #wav = wav.to(u.AA)

    wav = np.asarray(wav)*u.AA
    lum = np.asarray(flux)*u.erg/u.s

    #converting the SED to maggies to be saved in obs
    dl = 1*u.cm #for mck seds from fsps (already redshifted)
    flux = lum/(4.*3.14*dl**2.)
    nu = constants.c.cgs/(wav.to(u.cm))
    nu = nu.to(u.Hz)
    flux /= nu
    flux = flux.to(u.Jy)
    maggies = flux / 3631.

    ###implementing sedpy filter integration
    ###all code changes occur from here on

    filters_unsorted = load_filters(filternames)
    waves_unsorted = [x.wave_mean for x in filters_unsorted]
    filters = [x for _,x in sorted(zip(waves_unsorted,filters_unsorted))]
    from sedpy import observate
    mags = observate.getSED(np.asarray(df['wav']), np.asarray(df['spec']), filters)
    #dividing spectra to give photometry using filters as truth
   #https://github.com/bd-j/prospector/blob/9031e1b019e6aa87a0a7534a9a2ee8540da70033/prospect/utils/obsutils.py#L128 line 128
    phot = np.atleast_1d(mags)
    wave_eff = np.asarray([f.wave_effective for f in filters])*u.AA
    model_phot = 10. ** (phot / (-2.5))  # converts from magnitudes to flux

    #convert from cgs to maggies
    model_phot = np.asarray(model_phot)*u.erg/u.s
    model_phot = model_phot/(4.*3.14*dl**2.)
    nu = constants.c.cgs/(wave_eff.to(u.cm))
    nu = nu.to(u.Hz)
    model_phot/=nu
    model_phot = model_phot.to(u.Jy)
    model_maggies = model_phot/ 3631.
    model_maggies_unc = model_maggies*.03

    obs = {}
    obs['filters'] = filters
    obs['maggies'] = model_maggies.value
    obs['maggies_unc'] = model_maggies_unc.value
    obs['phot_mask'] = np.isfinite(model_maggies.value)
    obs['wavelength'] = None
    obs['spectrum'] = None
    obs['pd_sed'] = maggies
    obs['pd_wav'] = np.asarray(df['spec'])

    return obs

[jrleja]

Hi Grace,

I'm not sure I can help very much unfortunately! This looks like a tough issue and my suspicion is that it's related to turning simulation data into physical SEDs, which is largely outside of the Prospector machinery. I'd check that the photometry you're synthesizing is consistent with the original spectrum, and I'd also check that putting the underlying parameters into a very simple Prospector SED model (i.e., matching mass, SFH, metallicity, hopefully no dust!) produces a spectrum which is (broadly) consistent with the output from the simulation. This can help sort out some of these issues.

All Prospector can do here is try to generate a model spectrum which is in the best agreement with the input data as possible, and unless there's a suggestion that that's not happening (i.e. there is a solution that Prospector should find but isn't), there's not much we can do! Happy hunting...

Best, Joel

[gracekrahm]

While some of the underlying parameters are probably part of the larger issue, I'm more immediately concerned with the unit conversion issues right now where it seems to not really be matching anything.

    #mags to cgs
    mags = observate.getSED(np.asarray(df['wav']), np.asarray(df['spec']), filters) #wav and spec from simple fsps SED
    phot = np.atleast_1d(mags)
    wave_eff = np.asarray([f.wave_effective for f in filters])*u.AA
    model_phot = 10. ** (phot / (-2.5))  # converts from magnitudes to flux

    #cgs to maggies
    model_phot = np.asarray(model_phot)*u.erg/u.s
    model_phot = model_phot/(4.*3.14*dl**2.)
    nu = constants.c.cgs/(wave_eff.to(u.cm))
    nu = nu.to(u.Hz)
    model_phot/=nu
    model_phot = model_phot.to(u.Jy)
    model_maggies = model_phot/ 3631.

[bd-j]

Hi Grace, what are the units of df['wav'] and df['spec']? The getSED function assumes units of observed frame AA and erg/s/cm^2/AA (i.e. f_lambda) respectively in order to produce AB magnitudes. I suggest doing the unit conversions on the input spectrum before calling getSED.

[gracekrahm]

Those are the units of the input spectrum already.On Jun 18, 2024, at 8:28 AM, Ben Johnson @.***> wrote: Hi Grace, what are the units of df['wav'] and df['spec']? The getSED function assumes units of observed frame AA and erg/s/cm^2/AA (i.e. f_lambda) respectively in order to produce AB magnitudes. I suggest doing the unit conversions on the input spectrum before calling getSED.

—Reply to this email directly, view it on GitHub, or unsubscribe.You are receiving this because you authored the thread.Message ID: @.***>

[bd-j]

Ok, then you shouldn't need any of the # cgs to maggies block of code. phot, the output of getSED, should then be in AB magnitudes, you can convert to maggies just by

obs["maggies"] = 10**(-phot/2.5)
obs["maggies_unc"] = 0.03 * 10**(-phot/2.5)