ProTools

First this to do is get all the data you need for this workshop, it is at ProTools Workshop.

Next install and load the CRAN packages we need:

install.packages('remotes')
install.packages('data.table')
install.packages('Rcpp')
install.packages('foreach')
install.packages('magicaxis')
install.packages('celestial')

library(remotes)
library(data.table)
library(Rcpp)
library(foreach)
library(magicaxis)
library(celestial)

## Loading required package: RANN

## Loading required package: NISTunits

## Loading required package: pracma

All of the ProTools packages here are on my GitHub at asgr. To install a package you do not have run:

install_github('asgr/PUT_PACKAGE_NAME_HERE')
library(PUT_PACKAGE_NAME_HERE)

Where obviously PUT_PACKAGE_NAME_HERE should be replaced with the target packages name, like Rfits and ProFound below:

install_github('asgr/Rfits')
install_github('asgr/ProFound')

Assuming we have them name, we can load them.

library(Rfits)
library(ProFound)

Below we will generally load the packages at the point we need them. Either go ahead and get all of the ones you do not have now, or just wait until we get there. The other ones we will use are ProPane (including Rwcs), ProSpect (including ProSpectData), and Highlander (including LaplacesDemon and cmaeshpc). Some of these also have their own extra dependencies though (and these will be printed out when you are trying to install them).

Get Images Into Shape

Now we will set-up a few things for our particular example data:

setwd('~/Dropbox (Personal)/ProTools_Workshop/') #set as locally appropriate

GAMA_filters = c('FUV', 'NUV', 'u', 'g', 'r', 'i', 'Z', 'Y', 'J', 'H', 'Ks', 'W1', 'W2',
                 'W3', 'W4', 'P100', 'P160', 'S250', 'S350', 'S500')

GAMA_magzero = c(18.82, 20.08, 0, 0, 0, 0, 30, 30, 30, 30, 30, 23.16, 22.82, 23.24, 19.6,
                 8.9, 8.9, 11.68, 11.67, 11.62)

GAMA_piv = c(1535.08, 2300.785, 3567.001, 4734.851, 6194.748, 7509.465, 8833.013, 10223.63,
             12545.99, 16477, 21549.2, 33897.09, 46406.42, 125700.3, 223141.6, 983273.3,
             1546806, 2482806, 3484355, 5004155)

target_stub = 'gama_dr2_202627'
target_redshift = 0.05141
filepaths = paste0('Data_Image/', target_stub, '_', GAMA_filters, '.fits')

We can extract useful info from the various images with Rfits_key_scan:

ImageInfo = Rfits_key_scan(filelist = filepaths, get_all = TRUE)
ImageInfo

##                                                                                       full
##                                                                                     <char>
##  1:  /Users/aaron/Dropbox (Personal)/ProTools_Workshop/Data_Image/gama_dr2_202627_FUV.fits
##  2:  /Users/aaron/Dropbox (Personal)/ProTools_Workshop/Data_Image/gama_dr2_202627_NUV.fits
##  3:    /Users/aaron/Dropbox (Personal)/ProTools_Workshop/Data_Image/gama_dr2_202627_u.fits
##  4:    /Users/aaron/Dropbox (Personal)/ProTools_Workshop/Data_Image/gama_dr2_202627_g.fits
##  5:    /Users/aaron/Dropbox (Personal)/ProTools_Workshop/Data_Image/gama_dr2_202627_r.fits
##  6:    /Users/aaron/Dropbox (Personal)/ProTools_Workshop/Data_Image/gama_dr2_202627_i.fits
##  7:    /Users/aaron/Dropbox (Personal)/ProTools_Workshop/Data_Image/gama_dr2_202627_Z.fits
##  8:    /Users/aaron/Dropbox (Personal)/ProTools_Workshop/Data_Image/gama_dr2_202627_Y.fits
##  9:    /Users/aaron/Dropbox (Personal)/ProTools_Workshop/Data_Image/gama_dr2_202627_J.fits
## 10:    /Users/aaron/Dropbox (Personal)/ProTools_Workshop/Data_Image/gama_dr2_202627_H.fits
## 11:   /Users/aaron/Dropbox (Personal)/ProTools_Workshop/Data_Image/gama_dr2_202627_Ks.fits
## 12:   /Users/aaron/Dropbox (Personal)/ProTools_Workshop/Data_Image/gama_dr2_202627_W1.fits
## 13:   /Users/aaron/Dropbox (Personal)/ProTools_Workshop/Data_Image/gama_dr2_202627_W2.fits
## 14:   /Users/aaron/Dropbox (Personal)/ProTools_Workshop/Data_Image/gama_dr2_202627_W3.fits
## 15:   /Users/aaron/Dropbox (Personal)/ProTools_Workshop/Data_Image/gama_dr2_202627_W4.fits
## 16: /Users/aaron/Dropbox (Personal)/ProTools_Workshop/Data_Image/gama_dr2_202627_P100.fits
## 17: /Users/aaron/Dropbox (Personal)/ProTools_Workshop/Data_Image/gama_dr2_202627_P160.fits
## 18: /Users/aaron/Dropbox (Personal)/ProTools_Workshop/Data_Image/gama_dr2_202627_S250.fits
## 19: /Users/aaron/Dropbox (Personal)/ProTools_Workshop/Data_Image/gama_dr2_202627_S350.fits
## 20: /Users/aaron/Dropbox (Personal)/ProTools_Workshop/Data_Image/gama_dr2_202627_S500.fits
##                                                                                       full
##                          file                 stub
##                        <char>               <char>
##  1:  gama_dr2_202627_FUV.fits  gama_dr2_202627_FUV
##  2:  gama_dr2_202627_NUV.fits  gama_dr2_202627_NUV
##  3:    gama_dr2_202627_u.fits    gama_dr2_202627_u
##  4:    gama_dr2_202627_g.fits    gama_dr2_202627_g
##  5:    gama_dr2_202627_r.fits    gama_dr2_202627_r
##  6:    gama_dr2_202627_i.fits    gama_dr2_202627_i
##  7:    gama_dr2_202627_Z.fits    gama_dr2_202627_Z
##  8:    gama_dr2_202627_Y.fits    gama_dr2_202627_Y
##  9:    gama_dr2_202627_J.fits    gama_dr2_202627_J
## 10:    gama_dr2_202627_H.fits    gama_dr2_202627_H
## 11:   gama_dr2_202627_Ks.fits   gama_dr2_202627_Ks
## 12:   gama_dr2_202627_W1.fits   gama_dr2_202627_W1
## 13:   gama_dr2_202627_W2.fits   gama_dr2_202627_W2
## 14:   gama_dr2_202627_W3.fits   gama_dr2_202627_W3
## 15:   gama_dr2_202627_W4.fits   gama_dr2_202627_W4
## 16: gama_dr2_202627_P100.fits gama_dr2_202627_P100
## 17: gama_dr2_202627_P160.fits gama_dr2_202627_P160
## 18: gama_dr2_202627_S250.fits gama_dr2_202627_S250
## 19: gama_dr2_202627_S350.fits gama_dr2_202627_S350
## 20: gama_dr2_202627_S500.fits gama_dr2_202627_S500
##                          file                 stub
##                                                             path   ext length
##                                                           <char> <num>  <int>
##  1: /Users/aaron/Dropbox (Personal)/ProTools_Workshop/Data_Image     1     28
##  2: /Users/aaron/Dropbox (Personal)/ProTools_Workshop/Data_Image     1     28
##  3: /Users/aaron/Dropbox (Personal)/ProTools_Workshop/Data_Image     1     28
##  4: /Users/aaron/Dropbox (Personal)/ProTools_Workshop/Data_Image     1     28
##  5: /Users/aaron/Dropbox (Personal)/ProTools_Workshop/Data_Image     1     28
##  6: /Users/aaron/Dropbox (Personal)/ProTools_Workshop/Data_Image     1     28
##  7: /Users/aaron/Dropbox (Personal)/ProTools_Workshop/Data_Image     1     30
##  8: /Users/aaron/Dropbox (Personal)/ProTools_Workshop/Data_Image     1     30
##  9: /Users/aaron/Dropbox (Personal)/ProTools_Workshop/Data_Image     1     30
## 10: /Users/aaron/Dropbox (Personal)/ProTools_Workshop/Data_Image     1     30
## 11: /Users/aaron/Dropbox (Personal)/ProTools_Workshop/Data_Image     1     30
## 12: /Users/aaron/Dropbox (Personal)/ProTools_Workshop/Data_Image     1     28
## 13: /Users/aaron/Dropbox (Personal)/ProTools_Workshop/Data_Image     1     28
## 14: /Users/aaron/Dropbox (Personal)/ProTools_Workshop/Data_Image     1     28
## 15: /Users/aaron/Dropbox (Personal)/ProTools_Workshop/Data_Image     1     28
## 16: /Users/aaron/Dropbox (Personal)/ProTools_Workshop/Data_Image     1     28
## 17: /Users/aaron/Dropbox (Personal)/ProTools_Workshop/Data_Image     1     28
## 18: /Users/aaron/Dropbox (Personal)/ProTools_Workshop/Data_Image     1     28
## 19: /Users/aaron/Dropbox (Personal)/ProTools_Workshop/Data_Image     1     28
## 20: /Users/aaron/Dropbox (Personal)/ProTools_Workshop/Data_Image     1     28
##                                                             path   ext length
##     dim_1 dim_2 dim_3 dim_4 centre_RA centre_Dec rotation_North rotation_East
##     <int> <int> <int> <int>     <num>      <num>          <num>         <num>
##  1:    80    80    NA    NA  130.6180 -0.2716429              0           180
##  2:    80    80    NA    NA  130.6180 -0.2716429              0           180
##  3:   600   600    NA    NA  130.6178 -0.2715875              0           180
##  4:   600   600    NA    NA  130.6178 -0.2715875              0           180
##  5:   600   600    NA    NA  130.6178 -0.2715875              0           180
##  6:   600   600    NA    NA  130.6178 -0.2715875              0           180
##  7:   354   354    NA    NA  130.6178 -0.2715797              0           180
##  8:   354   354    NA    NA  130.6178 -0.2715797              0           180
##  9:   354   354    NA    NA  130.6178 -0.2715797              0           180
## 10:   354   354    NA    NA  130.6178 -0.2715797              0           180
## 11:   354   354    NA    NA  130.6178 -0.2715797              0           180
## 12:   120   120    NA    NA  130.6178 -0.2716429              0           180
## 13:   120   120    NA    NA  130.6178 -0.2716429              0           180
## 14:   120   120    NA    NA  130.6178 -0.2716429              0           180
## 15:   120   120    NA    NA  130.6178 -0.2716429              0           180
## 16:    40    40    NA    NA  130.6178 -0.2716429              0           180
## 17:    30    30    NA    NA  130.6178 -0.2713660              0           180
## 18:    20    20    NA    NA  130.6178 -0.2716429              0           180
## 19:    15    15    NA    NA  130.6189 -0.2713664              0           180
## 20:    10    10    NA    NA  130.6194 -0.2724742              0           180
##     dim_1 dim_2 dim_3 dim_4 centre_RA centre_Dec rotation_North rotation_East
##     corner_BL_RA corner_BL_Dec corner_TL_RA corner_TL_Dec corner_TR_RA
##            <num>         <num>        <num>         <num>        <num>
##  1:     130.6346    -0.2882632     130.6346    -0.2550349     130.6014
##  2:     130.6346    -0.2882632     130.6346    -0.2550349     130.6014
##  3:     130.6344    -0.2882077     130.6343    -0.2549795     130.6012
##  4:     130.6344    -0.2882077     130.6343    -0.2549795     130.6012
##  5:     130.6344    -0.2882077     130.6343    -0.2549795     130.6012
##  6:     130.6344    -0.2882077     130.6343    -0.2549795     130.6012
##  7:     130.6344    -0.2882008     130.6344    -0.2549709     130.6013
##  8:     130.6344    -0.2882008     130.6344    -0.2549709     130.6013
##  9:     130.6344    -0.2882008     130.6344    -0.2549709     130.6013
## 10:     130.6344    -0.2882008     130.6344    -0.2549709     130.6013
## 11:     130.6344    -0.2882008     130.6344    -0.2549709     130.6013
## 12:     130.6344    -0.2882631     130.6343    -0.2550349     130.6012
## 13:     130.6344    -0.2882631     130.6343    -0.2550349     130.6012
## 14:     130.6344    -0.2882631     130.6343    -0.2550349     130.6012
## 15:     130.6344    -0.2882631     130.6343    -0.2550349     130.6012
## 16:     130.6344    -0.2882631     130.6343    -0.2550349     130.6012
## 17:     130.6344    -0.2879862     130.6343    -0.2547580     130.6012
## 18:     130.6344    -0.2882631     130.6343    -0.2550349     130.6012
## 19:     130.6355    -0.2879866     130.6355    -0.2547583     130.6023
## 20:     130.6360    -0.2890945     130.6360    -0.2558662     130.6029
##     corner_BL_RA corner_BL_Dec corner_TL_RA corner_TL_Dec corner_TR_RA
##     corner_TR_Dec corner_BR_RA corner_BR_Dec   min_RA    min_Dec   max_RA
##             <num>        <num>         <num>    <num>      <num>    <num>
##  1:    -0.2550233     130.6014    -0.2882501 130.6014 -0.2882632 130.6346
##  2:    -0.2550233     130.6014    -0.2882501 130.6014 -0.2882632 130.6346
##  3:    -0.2549679     130.6012    -0.2881946 130.6012 -0.2882077 130.6344
##  4:    -0.2549679     130.6012    -0.2881946 130.6012 -0.2882077 130.6344
##  5:    -0.2549679     130.6012    -0.2881946 130.6012 -0.2882077 130.6344
##  6:    -0.2549679     130.6012    -0.2881946 130.6012 -0.2882077 130.6344
##  7:    -0.2549593     130.6013    -0.2881877 130.6013 -0.2882008 130.6344
##  8:    -0.2549593     130.6013    -0.2881877 130.6013 -0.2882008 130.6344
##  9:    -0.2549593     130.6013    -0.2881877 130.6013 -0.2882008 130.6344
## 10:    -0.2549593     130.6013    -0.2881877 130.6013 -0.2882008 130.6344
## 11:    -0.2549593     130.6013    -0.2881877 130.6013 -0.2882008 130.6344
## 12:    -0.2550232     130.6012    -0.2882500 130.6012 -0.2882631 130.6344
## 13:    -0.2550232     130.6012    -0.2882500 130.6012 -0.2882631 130.6344
## 14:    -0.2550232     130.6012    -0.2882500 130.6012 -0.2882631 130.6344
## 15:    -0.2550232     130.6012    -0.2882500 130.6012 -0.2882631 130.6344
## 16:    -0.2550232     130.6012    -0.2882500 130.6012 -0.2882631 130.6344
## 17:    -0.2547464     130.6012    -0.2879731 130.6012 -0.2879862 130.6344
## 18:    -0.2550232     130.6012    -0.2882500 130.6012 -0.2882631 130.6344
## 19:    -0.2547467     130.6023    -0.2879735 130.6023 -0.2879866 130.6355
## 20:    -0.2558545     130.6029    -0.2890813 130.6029 -0.2890945 130.6360
##     corner_TR_Dec corner_BR_RA corner_BR_Dec   min_RA    min_Dec   max_RA
##        max_Dec range_RA range_Dec   pixscale      pixarea
##          <num>    <num>     <num>      <num>        <num>
##  1: -0.2550233 1.988248  1.994393  1.4929570   2.22891471
##  2: -0.2550233 1.988248  1.994393  1.4929570   2.22891471
##  3: -0.2549679 1.988247  1.994392  0.1990608   0.03962511
##  4: -0.2549679 1.988247  1.994392  0.1990608   0.03962511
##  5: -0.2549679 1.988247  1.994392  0.1990608   0.03962511
##  6: -0.2549679 1.988247  1.994392  0.1990608   0.03962511
##  7: -0.2549593 1.988347  1.994492  0.3374081   0.11384395
##  8: -0.2549593 1.988347  1.994492  0.3374081   0.11384395
##  9: -0.2549593 1.988347  1.994492  0.3374081   0.11384395
## 10: -0.2549593 1.988347  1.994492  0.3374081   0.11384395
## 11: -0.2549593 1.988347  1.994492  0.3374081   0.11384395
## 12: -0.2550232 1.988247  1.994392  0.9953042   0.99062783
## 13: -0.2550232 1.988247  1.994392  0.9953042   0.99062783
## 14: -0.2550232 1.988247  1.994392  0.9953042   0.99062783
## 15: -0.2550232 1.988247  1.994392  0.9953042   0.99062783
## 16: -0.2550232 1.988247  1.994392  2.9859130   8.91565357
## 17: -0.2547464 1.988247  1.994392  3.9812180  15.85005565
## 18: -0.2550232 1.988247  1.994392  5.9718277  35.66263304
## 19: -0.2547467 1.988253  1.994395  7.9624569  63.40055600
## 20: -0.2558545 1.988256  1.994398 11.9437000 142.65160072
##        max_Dec range_RA range_Dec   pixscale      pixarea

Notice how the pixel resolution varies quite a lot.

We can load in all the FITS files easily:

ImageList = Rfits_make_list(filelist = filepaths, pointer = FALSE)
names(ImageList) = GAMA_filters

And look at some:

for(i in seq_along(ImageList)){
  par(mar=c(3.1,3.1,1.1,1.1))
  plot(ImageList[[i]])
  legend('topleft', legend=names(ImageList)[i])
}

We can run a detection on an image:

pro_Z = profoundProFound(ImageList$Z, plot=TRUE, sky=0, magzero = 30)

It’s not perfect (some of the floculent structure to the North is missed) but in terms of flux we have captured most the galaxy light. Note we set the initial sky to 0 since in theory the image has already been background subtracted, you can remove the sky argument if you want profoundProFound to calculate the sky itself.

You can try experimenting with some of the options that affect the source finding (and see if you can do better!) The main ones to explore are skycut, tolerance, reltol, sigma and cliptol (see the help at ?profoundProFound for details on what these do).

profoundProFound(ImageList$Z, plot=TRUE, sky=0, magzero = 30, skycut=1, tolerance = 1)

profoundProFound(ImageList$Z, plot=TRUE, sky=0, magzero = 30, skycut=2, tolerance = 10)

Notice how we extract more sources with lower skycut, and we tend to blend things together more with higher tolerance.

We are going to warp the FUV-W2 images to the same resolution as our Z band:

library(ProPane)

MatchRes = foreach(i = 1:13)%do%{
  GAMA_filter = GAMA_filters[i]
  if(i %in% c('Z', 'Y', 'J', 'H','Ks')){
    return(ImageList[[i]])
  }else{
    return(propaneWarp(ImageList[[i]], keyvalues_out = ImageList$Z$keyvalues, magzero_in = GAMA_magzero[i], magzero_out = 8.9))
  }
}
names(MatchRes) = GAMA_filters[1:13]

Now we can run in multiband mode, where we can also detect on a stack of the data (Z, Y, J, H, Ks in our case)

pro_MatchRes = profoundMultiBand(
  inputlist = MatchRes,
  detectbands = c('Z', 'Y', 'J', 'H', 'Ks'),
  multibands = names(MatchRes),
  magzero = 8.9,
  sky = 0
  )

## *** Will use Z+Y+J+H+Ks for source detection ***

## *** Will use detection image for colour photometry segim_orig ***

## *** Will use FUV NUV u g r i Z Y J H Ks W1 W2 for multi band photometry ***

## *** Magzero: FUV=8.9 NUV=8.9 u=8.9 g=8.9 r=8.9 i=8.9 Z=8.9 Y=8.9 J=8.9 H=8.9 Ks=8.9 W1=8.9 W2=8.9  ***

## *** Gain: not specified for any bands, so shot-noise will be ignored ***

## *** Will compute total multi band photometry ***

## *** Will compute isophotal colour multi band photometry ***

## *** Will compute grouped segment multi band photometry ***

## *** Currently processing detection band Z ***

## *** Currently processing detection band Y ***

## *** Currently processing detection band J ***

## *** Currently processing detection band H ***

## *** Currently processing detection band Ks ***

## *** Currently processing stacked detection image ***

## *** Currently processing multi band FUV ***

## * Processing total photometry *

## * Processing colour photometry *

## * Processing group photometry *

## *** Currently processing multi band NUV ***

## * Processing total photometry *

## * Processing colour photometry *

## * Processing group photometry *

## *** Currently processing multi band u ***

## * Processing total photometry *

## * Processing colour photometry *

## * Processing group photometry *

## *** Currently processing multi band g ***

## * Processing total photometry *

## * Processing colour photometry *

## * Processing group photometry *

## *** Currently processing multi band r ***

## * Processing total photometry *

## * Processing colour photometry *

## * Processing group photometry *

## *** Currently processing multi band i ***

## * Processing total photometry *

## * Processing colour photometry *

## * Processing group photometry *

## *** Currently processing multi band Z ***

## * Processing total photometry *

## * Processing colour photometry *

## * Processing group photometry *

## *** Currently processing multi band Y ***

## * Processing total photometry *

## * Processing colour photometry *

## * Processing group photometry *

## *** Currently processing multi band J ***

## * Processing total photometry *

## * Processing colour photometry *

## * Processing group photometry *

## *** Currently processing multi band H ***

## * Processing total photometry *

## * Processing colour photometry *

## * Processing group photometry *

## *** Currently processing multi band Ks ***

## * Processing total photometry *

## * Processing colour photometry *

## * Processing group photometry *

## *** Currently processing multi band W1 ***

## * Processing total photometry *

## * Processing colour photometry *

## * Processing group photometry *

## *** Currently processing multi band W2 ***

## * Processing total photometry *

## * Processing colour photometry *

## * Processing group photometry *

TotDT = data.table(pro_MatchRes$cat_tot) #just because data.table is nicer

We can get the ID of the central source easily:

CentralID = pro_MatchRes$pro_detect$segim[177,177]

And now get the photometry for that objects:

CentralPhotom = TotDT[segID == CentralID, ]
flux_names = paste0('flux_', names(MatchRes), 't')
flux_err_names = paste0('flux_err_', names(MatchRes), 't')

Flux = as.numeric(CentralPhotom[,..flux_names])
Flux_err = as.numeric(CentralPhotom[,..flux_err_names])
Flux_err = sqrt(Flux_err^2 + (0.05*Flux)^2) #gives a 5% error floor

Output_res = data.frame(filter = names(MatchRes), pivwave=GAMA_piv[1:13], flux=Flux, fluxerr=Flux_err)

magplot(Output_res$pivwave, Output_res$flux, log='xy', xlab='Wavelength / Ang', ylab='Flux Density / Jansky')
magerr(Output_res$pivwave, Output_res$flux, ylo=Output_res$fluxerr)

For the other bands the Flux is the sum of what is in our cutout, and the error will come from the skyRMS.

In reality for these bands we would normally use profoundFitMagPSF to extract the flux correctly (with errors), but this will do for now.

Flux_unres = foreach(i = 14:20, .combine=c)%do%{
  sum(ImageList[[i]]$imDat, na.rm=TRUE) * 10^(0.4*(8.9 - GAMA_magzero[i]))
}

Flux_err_unres = foreach(i = 14:20, .combine=c)%do%{
  sqrt(diff(quantile(ImageList[[i]]$imDat, c(0.16, 0.5)))^2 * prod(dim(ImageList[[i]]$imDat))) * 10^(0.4*(8.9 - GAMA_magzero[i]))
}

Flux_err_unres = sqrt(Flux_err_unres^2 + (0.1*Flux_unres)^2)

Output_unres = data.frame(filter = names(ImageList)[14:20], pivwave=GAMA_piv[14:20], flux=Flux_unres, fluxerr=Flux_err_unres)

Output_comb = rbind(Output_res, Output_unres)

Let us plot it!

magplot(Output_comb$pivwave, Output_comb$flux, log='xy', xlab='Wavelength / Ang', ylab='Flux Density / Jansky')
magerr(Output_comb$pivwave, Output_comb$flux, ylo=Output_comb$fluxerr)

Fit All The Things!

Now for some ProSpect

library(ProSpect)

## Loading required package: ProSpectData

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

Load the filters we need:

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))}
GAMA_piv = foreach(i = filters, .combine='c')%do%{pivwavefunc(getfilt(i))}

We actually typed in the pivot wavelengths manually earlier for the plot, but the above is a more generic and pipeline friendly way to extract them from the filters directly (and should ensure we make fewer mistake).

We will set up our data structure using BC03:

Data_BC03 = list(flux = Output_comb,
          arglist = list(z=target_redshift, massfunc=massfunc_snorm_trunc, agemax=agemax,
                         magemax=agemax/1e9, Z=Zfunc_massmap_lin, Zfinal=0.02,
                         LumDist_Mpc=LumDist_Mpc),
          speclib = BC03lr, 
          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',paste('flux.',Output_comb$filter,sep='')),
          N = length(filters), # number of observed filters
          like = 'norm',
          verbose = FALSE
  )

Note we are using the linear mapping form for our metallicity evolution, but setting Zfinal to be solar (Z=0.02). We could also fit for the value of Zfinal by adding it to our parameter list (this is what we do for GAMA and DEVILS), but for this example we will keep things simpler.

We will fit this with Highlander:

library(Highlander) #Note LaplacesDemon and cmaeshpc are both on GitHub/asgr
badpar = (Data_BC03$intervals$lo + Data_BC03$intervals$hi) / 2
Highout_BC03 = Highlander(badpar, Data=Data_BC03, likefunc=ProSpectSEDlike, Niters=c(1e3,1e3))
saveRDS(Highout_BC03, 'Data_SSP/BC03_fit.rds')

Note I have saved my fitting run into Data_SSP/BC03_fit.rds, to load it you can run:

Highout_BC03 = readRDS('Data_SSP/BC03_fit.rds')

Plot the output!

Data_BC03$fit = 'check' #we have to set the fit type to 'check' to get the required outputs
fit_outputs_BC03 = ProSpectSEDlike(parm=Highout_BC03$parm, Data=Data_BC03)
plot(fit_outputs_BC03)

plot(fit_outputs_BC03$SEDout)

Our preferred solution shows the SFH steadily growing over time.

Let’s have a go at loading a ProGeny SSP and using that instead:

PG_Ch_Mi_C3K = speclib_FITSload('Data_SSP/PG_Ch_Mi_C3K.fits')
PG_Ch_Mi_C3K = speclibReBin(PG_Ch_Mi_C3K, BC03lr$Wave) #lower res

Insert this into our data structure:

Data_PG = Data_BC03
Data_PG$speclib = PG_Ch_Mi_C3K

And away we fit!

Highout_PG = Highlander(badpar, Data=Data_PG, likefunc=ProSpectSEDlike, Niters=c(1e3,1e3))
saveRDS(Highout_PG, 'Data_SSP/PG_fit.rds')

Note I have saved my fitting run into Data_SSP/PG_fit.rds, to load it you can run:

Highout_PG = readRDS('Data_SSP/PG_fit.rds')

Data_PG$fit = 'check' #we have to set the fit type to 'check' to get the required outputs
fit_outputs_PG = ProSpectSEDlike(parm=Highout_PG$parm, Data=Data_PG)
plot(fit_outputs_PG)

plot(fit_outputs_PG$SEDout)

Note the BC03 fit is actually a tiny bit better (log-likelihood [LL] -34.312 versus -34.659, where a larger absolute LL is more likely).

If we want to extract star formation rates and stellar masses we need to get our parameters into the correct (linear) format and then use a helper:

mSFR = 10^Highout_PG$parm[1]
mpeak = 10^Highout_PG$parm[2]
mperiod = 10^Highout_PG$parm[3]
mskew = Highout_PG$parm[4]

We can now run our helper function with the correct inputs:

PG_SMout = SMstarfunc(massfunc_snorm_trunc, mSFR=mSFR, mpeak=mpeak, mperiod=mperiod,
                      mskew=mskew, magemax=agemax/1e9, speclib=PG_Ch_Mi_C3K,
                      Z = Zfunc_massmap_lin, Zfinal = 0.02)

Now we can compare the mass formed with the mass remaining:

log10(PG_SMout['TotSMform'])

## TotSMform 
##  10.81222

log10(PG_SMout['TotSMstar'])

## TotSMstar 
##  10.60297

Some people new to the world of SED fitting might wonder why the total mass formed and the mass remaining in stellar material can differ this much. The answer is stellar recycling, things moving off the main sequence and blowing out mass. Different SSPs will predict different functions for this behaviour and the conversion depends on the specific SFH (hence it is not trivial to convert mass formed):

magplot(PG_Ch_Mi_C3K$Age, PG_Ch_Mi_C3K$Zevo[[13]][,"SMstar"], log='x', type='l',
        xlab='Age / Yrs', ylab='Fraction of Mass in Stars', ylim=c(0,1))

The older GAMA value for the stellar mass formed was around 10.46 +/- 0.13. Compared to that work ProFound tends to extract more flux (so implying a more massive galaxy) and ProSpect itself fits an older SFH with higher mass to light when using the specific SFH form used here.

Perhaps more conveniently the stellar mass extraction can all be done in one go with ParmOff:

library(ParmOff)
SFH_best = Highout_PG$parm
names(SFH_best) = Data_PG$parm.names
PG_SMout2 = ParmOff(.func = SMstarfunc, #the function we want to run
                    .args = SFH_best, #the superset of potential matching parameters
                    .logged = c('mSFR', 'mpeak', 'mperiod'), #parameters we want to log
                    massfunc = massfunc_snorm_trunc, #the SFH function to use
                    magemax = agemax/1e9, #specify agemax (defined above)
                    speclib = PG_Ch_Mi_C3K, #specify the SSP to use
                    Z = Zfunc_massmap_lin, #the form of the metallicity history assumed
                    Zfinal = 0.02 #the setting for Zfinal (solar)
                    )

sum(PG_SMout - PG_SMout2)

## [1] 0

The sum of the difference is 0, so things are agreeing.

A Different Kind of Fit

Now we are going to fit using a different kind of star formation history: an exponentially declining SFH with a recent burst (massfunc_exp_burst):

Data_PG2 = list(flux = Output_comb,
          arglist = list(z=target_redshift, massfunc=massfunc_exp_burst, agemax=agemax,
                         Z=Zfunc_massmap_lin, LumDist_Mpc=LumDist_Mpc),
          speclib = PG_Ch_Mi_C3K, 
          Dale = Dale_NormTot, 
          filtout = filtout, 
          SFH = SFHfunc, # the preferred functional form of the SFH (eg either SFHfunc, SFHburst)
          parm.names = c('mburst', 'mSFR','mtau','magemax','tau_birth','tau_screen',
                       'alpha_SF_birth','alpha_SF_screen'), # which parameters to fit for
          logged = c(T,T,F,F,T,T,F,F), # fit parameters in logged or linear space
          intervals = list(lo=c(-4,-4,0,1,-2.5,-2.5,0.0625,0.0625), hi=c(2,3,10,10,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',paste('flux.',Output_comb$filter,sep='')),
          N = length(filters), # number of observed filters
          like = 'norm',
          verbose = FALSE
  )

And away we fit!

badpar2 = (Data_PG2$intervals$lo + Data_PG2$intervals$hi) / 2
Highout_PG2 = Highlander(badpar2, Data=Data_PG2, likefunc=ProSpectSEDlike, Niters=c(1e3,1e3))
saveRDS(Highout_PG2, 'Data_SSP/PG_fit2.rds')

Note I have saved my fitting run into Data_SSP/PG_fit2.rds, to load it you can run:

Highout_PG2 = readRDS('Data_SSP/PG_fit2.rds')

Data_PG2$fit = 'check' #we have to set the fit type to 'check' to get the required outputs
fit_outputs_PG2 = ProSpectSEDlike(parm=Highout_PG2$parm, Data=Data_PG2)
plot(fit_outputs_PG2)

plot(fit_outputs_PG2$SEDout)

mburst = 10^Highout_PG2$parm[1]
mSFR = 10^Highout_PG2$parm[2]
mtau = Highout_PG2$parm[3]
magemax = Highout_PG2$parm[4]

We can now run our helper function with the correct inputs:

PG_SMout3 = SMstarfunc(massfunc_exp_burst, mburst=mburst, mSFR=mSFR, mtau=mtau,
                      magemax=magemax, speclib=PG_Ch_Mi_C3K, Z = Zfunc_massmap_lin,
                      Zfinal = 0.02)

And now we can compare the masses:

log10(PG_SMout3["TotSMform"])

## TotSMform 
##   10.9569

log10(PG_SMout3["TotSMstar"])

## TotSMstar 
##  10.73313

This is much larger stellar mass. This will largely be because of the presence of a strong recent burst. In general it is dangerous to include a burst for this reason - the data will give quite degenerate solutions (unless very very high S/N), and the harm it does is notable. Almost certainly the snorm solution is the more physical one, even though this is technically a better fit.

As a last example let’s fit a stepwise SFH:

Data_PG3 = list(flux = Output_comb,
          arglist = list(z=target_redshift, massfunc=massfunc_b5, agemax=agemax,
                         magemax=agemax/1e9, Z=Zfunc_massmap_lin, LumDist_Mpc=LumDist_Mpc),
          speclib = PG_Ch_Mi_C3K, 
          Dale = Dale_NormTot, 
          filtout = filtout, 
          SFH = SFHfunc, # the preferred functional form of the SFH (eg either SFHfunc, SFHburst)
          parm.names = c('m1', 'm2','m3','m4', 'm5','tau_birth','tau_screen',
                       'alpha_SF_birth','alpha_SF_screen'), # which parameters to fit for
          logged = c(T,T,T,T,T,T,T,F,F), # fit parameters in logged or linear space
          intervals = list(lo=c(-4,-4,-4,-4,-4,-2.5,-2.5,0.0625,0.0625), hi=c(3,3,3,3,3,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',paste('flux.',Output_comb$filter,sep='')),
          N = length(filters), # number of observed filters
          like = 'norm',
          verbose = FALSE
  )

And away we fit!

badpar3 = (Data_PG3$intervals$lo + Data_PG3$intervals$hi) / 2
Highout_PG3 = Highlander(badpar3, Data=Data_PG3, likefunc=ProSpectSEDlike, Niters=c(1e3,1e3))
saveRDS(Highout_PG3, 'Data_SSP/PG_fit3.rds')

Note I have saved my fitting run into Data_SSP/PG_fit3.rds, to load it you can run:

Highout_PG3 = readRDS('Data_SSP/PG_fit3.rds')

Data_PG3$fit = 'check' #we have to set the fit type to 'check' to get the required outputs
fit_outputs_PG3 = ProSpectSEDlike(parm=Highout_PG3$parm, Data=Data_PG3)
plot(fit_outputs_PG3)

plot(fit_outputs_PG3$SEDout)

m1 = 10^Highout_PG3$parm[1]
m2 = 10^Highout_PG3$parm[2]
m3 = 10^Highout_PG3$parm[3]
m4 = 10^Highout_PG3$parm[4]
m5 = 10^Highout_PG3$parm[5]

We can now run our helper function with the correct inputs:

PG_SMout3 = SMstarfunc(massfunc_b5, m1=m1, m2=m2, m3=m3, m4=m4, m5=m5,
                      magemax=agemax/1e9, speclib=PG_Ch_Mi_C3K, Z = Zfunc_massmap_lin,
                      Zfinal = 0.02)

And now we can compare the masses:

log10(PG_SMout3["TotSMform"])

## TotSMform 
##  10.77023

log10(PG_SMout3["TotSMstar"])

## TotSMstar 
##  10.55775

This is pretty consistent with the earlier snorm derived stellar mass.