Spectro-Photometric Fitting

Aaron Robotham

2025-09-30

It is becoming more common to have access to both flux calibrated spectra and good quality broad-band photometry covering a very wide range of rest-frame wavelength (certainly optical to near-infrared). Modern programmes with JWST (NIRcam, MIRI and NIRspec) and surveys like WAVES (4MOST spectra, plus UV-FIR photometry) ideally wish to combine information from both spectra and photometry when fitting SED models to galaxies.

Early on ProSpect was focussed on broad-band photometry modelling only, and in about 2023 it had capacity to fit spectra in a separate mode. Now (finally) you can fit both at once! Happy days. In this vignette we will look at a simple mock example of how this works in practice.

Making Mock Data

Load the libraries we will need:

library(ProSpect)
#> Loading required package: ProSpectData
#> Loading required package: Rcpp
library(LaplacesDemon)
library(foreach)
library(celestial)
#> Loading required package: RANN
#> Loading required package: NISTunits
#> Loading required package: pracma
#> 
#> Attaching package: 'pracma'
#> The following objects are masked from 'package:LaplacesDemon':
#> 
#>     logit, loglog, Mode
library(magicaxis)

Load some data we want to use:

data("BC03lr") #BC03 spectral library
data("Dale_NormTot") #Normalised Dale templates
data("pivwave") # Pivot/effective wavelengths
set.seed(1)
redshift = 0.1

filters=c('FUV_GALEX', 'NUV_GALEX', 'u_SDSS', 'g_SDSS', 'r_SDSS', 'i_SDSS', 'Z_VISTA',
          'Y_VISTA', 'J_VISTA', 'H_VISTA', 'K_VISTA', 'W1_WISE' , 'W2_WISE', 'W3_WISE',
          'W4_WISE', 'P100_Herschel', 'P160_Herschel', 'S250_Herschel' , 'S350_Herschel',
          'S500_Herschel')

filtout=foreach(i = filters)%do%{approxfun(getfilt(i))}

temppiv=pivwave[pivwave$filter %in% filters,]

agemax = 13.3e9 - cosdistTravelTime(z=redshift, H0 = 67.8, OmegaM = 0.308)*1e9

Create our model galaxy We will be using a snorm_trunc SFH function and a massmap_lin metallicity with the following parameters:

inpar=c(mSFR = 0, #log-space
        mpeak = 0.7, #log-space
        mperiod = 0.3, #log-space
        mskew = 0.3,
        tau_birth = 0, #log-space
        tau_screen = -0.5, #log-space
        alpha_SF_birth = 1,
        alpha_SF_screen = 3
        )

Simple function to view our target SFH:

plotSFH=function(par, magemax=13.3, add=FALSE,col='black',ylim=NULL,...){
 magcurve(massfunc_snorm_trunc(age=x,mSFR=10^par[1],mpeak=10^par[2],mperiod=10^par[3],
                               mskew=par[4], magemax=magemax),0,13.8e9,add=add,col=col,
                               ylim=ylim,xlab='Age (Yr)', ylab='SFR (Msol / Yr)',...)
}

Let’s take a look:

plotSFH(inpar)

Now we can generate our galaxies SED:

genSED=ProSpectSED(massfunc = massfunc_snorm_trunc,
                   mSFR = 10^inpar[1],
                   mpeak = 10^inpar[2],
                   mperiod = 10^inpar[3],
                   mskew = inpar[4],
                   tau_birth = 10^inpar[5], 
                   tau_screen = 10^inpar[6], 
                   alpha_SF_birth = inpar[7], 
                   alpha_SF_screen = inpar[8],
                   z = redshift,
                   Z = Zfunc_massmap_lin,
                   filtout = filtout,
                   Dale = Dale_NormTot,
                   speclib = BC03hr,
                   agemax = agemax,
                   magemax = agemax/1e9,
                   sparse = 1,
                   waveout = NULL
)

We also want a small bit of the spectrum:

spec_input = genSED$FinalFlux[genSED$FinalFlux$wave > 3800 & genSED$FinalFlux$wave < 1e4,]
spec_input[,'fluxerr'] = spec_input[,'flux']*0.05

At this point we create our mock photometry catalogue entry with a fractional error in flux of 0.1 assumed (roughly 0.092 mag error). The photom output of ProSpectSED is always Jansky, which is what we want for fitting (you shouldn’t fit in magnitude space!).

phot_input = data.frame(filter=temppiv$filter, pivwave=temppiv$pivwave, flux=genSED$Photom, fluxerr=genSED$Photom*0.1)
print(phot_input)
#>           filter     pivwave         flux      fluxerr
#> 1      FUV_GALEX    1535.080 3.310725e-06 3.310725e-07
#> 2      NUV_GALEX    2300.785 5.123891e-06 5.123891e-07
#> 3         u_SDSS    3567.001 1.105813e-05 1.105813e-06
#> 4         g_SDSS    4734.851 3.227777e-05 3.227777e-06
#> 5         r_SDSS    6194.748 5.530228e-05 5.530228e-06
#> 6         i_SDSS    7509.465 7.145967e-05 7.145967e-06
#> 7        Z_VISTA    8833.013 8.735418e-05 8.735418e-06
#> 8        Y_VISTA   10223.626 9.981373e-05 9.981373e-06
#> 9        J_VISTA   12545.988 1.117815e-04 1.117815e-05
#> 10       H_VISTA   16476.995 1.282882e-04 1.282882e-05
#> 11       K_VISTA   21549.195 1.228454e-04 1.228454e-05
#> 12       W1_WISE   33897.093 7.050606e-05 7.050606e-06
#> 13       W2_WISE   46406.422 4.584160e-05 4.584160e-06
#> 14       W3_WISE  125700.270 2.425128e-04 2.425128e-05
#> 15       W4_WISE  223141.602 5.831915e-04 5.831915e-05
#> 16 P100_Herschel  983273.314 6.693945e-03 6.693945e-04
#> 17 P160_Herschel 1546806.337 8.262902e-03 8.262902e-04
#> 18 S250_Herschel 2482806.324 5.258834e-03 5.258834e-04
#> 19 S350_Herschel 3484355.144 2.752377e-03 2.752377e-04
#> 20 S500_Herschel 5004155.150 9.370676e-04 9.370676e-05

Fitting a Galaxy

In principle if you have a real galaxy to fit and it is already in the required format of flux_input, you can start from this point.

To speed things up we will also pre-compute the Luminosity distance to our source (since this will be the same for every iteration, and will save some computation time):

LumDist_Mpc = cosdistLumDist(z=redshift, H0=67.8, OmegaM=0.308)
Data=list(flux = phot_input,
          spec = spec_input,
          mode = 'both',
          arglist=list(z=redshift, massfunc=massfunc_snorm_trunc, agemax=agemax, magemax=agemax/1e9,
                       Z=Zfunc_massmap_lin, LumDist_Mpc=LumDist_Mpc, sparse=1, waveout=NULL),
          speclib=BC03hr, 
          Dale=Dale_NormTot, 
          filtout=filtout, 
          SFH=SFHfunc, # the preferred functional form of the SFH (eg either SFHfunc, SFHburst)
          parm.names=c('mSFR','mpeak','mperiod','mskew','tau_birth','tau_screen',
                       'alpha_SF_birth','alpha_SF_screen'), # which parameters to fit for
          logged=c(T,T,T,F,T,T,F,F), # fit parameters in logged or linear space
          intervals=list(lo=c(-4,-2,-1,-0.5,-2.5,-2.5,0.0625,0.0625), hi=c(3,1,1,1,1,1,4,4)), # fitting range for parameters
          fit = 'LD', # specifies the way in which the SED should be fitted ('LD', 'optim', 'CMA', or 'check')
          mon.names=c('LP','masstot','SFRburst'),
          N=length(filters), # number of observed filters
          like='norm',
          verbose=FALSE
  )
ProSpectSEDlike(inpar, Data)
#> $LP
#> [1] -5197.516
#> 
#> $Dev
#> [1] 10395.03
#> 
#> $Monitor
#>            LP       masstot      SFRburst 
#> -5.197516e+03  4.984556e+09  3.152976e-01 
#> 
#> $yhat
#> [1] 1
#> 
#> $parm
#>            mSFR           mpeak         mperiod           mskew       tau_birth 
#>             0.0             0.7             0.3             0.3             0.0 
#>      tau_screen  alpha_SF_birth alpha_SF_screen 
#>            -0.5             1.0             3.0
Data$fit = 'check' #we have to set the fit type to 'check' to get the required outputs
inputs = ProSpectSEDlike(parm=inpar, Data=Data)
plot(inputs)
plot(inputs$SEDout, xlim=c(1e3,2e6), ylim=c(1e3,5e6))

So now we can fit both the spectra and the broadband photometry at the same time!

Data$fit = 'optim'
library(Highlander)
badpar = (Data$intervals$lo + Data$intervals$hi) / 2
Highout = Highlander(badpar, Data=Data, likefunc=ProSpectSEDlike, Niters=c(1e2,1e2), optim_iters=5)
Data$fit = 'check' #we have to set the fit type to 'check' to get the required outputs
outputs = ProSpectSEDlike(parm=Highout$parm, Data=Data)
plot(outputs)
plot(outputs$SEDout, xlim=c(1e3,2e6), ylim=c(1e3,5e6))