This notebook shows a quick demo of how to calculate the A_lambda/E(B-V) coefficients needed to remove Galactic dust from the DC2 catalogs.
In DC2, the CCM model (https://ui.adsabs.harvard.edu/abs/1989ApJ...345..245C/abstract) was assumed when calculating the dust model, we will use the lsst.sims.photUtils package to compute the effective wavelengths for the filters, and then set up the CCM model and calculate Alambda/E(B-V) for each of the LSST filters.
import os
import GCRCatalogs
import pandas as pd
import numpy as np
import scipy.interpolate
from lsst.sims import photUtils
from lsst.sims.photUtils import BandpassSet
This function grabs the "LSST_THROUGHPUTS_BASELINE" filters as currently specified in the photUtils package, then computes both lambda_eff and A_lambda/E(B-V) for each of the six LSST filters:
def compute_alambda_over_ebv(filterset='ugrizy',basepath=None):
"""
compute the effective wavelengths and alambda/E(B-V) values for a set of filters
We will grab the flat SED from the SIMS library to calculate the CCM dust model
that was assumed for DC2, and then grab the baseline ugrizy filters and calculate
their effective wavelengths, and evaluate the CCM alam_over_ebv value at those
wavelengths
inputs: filterset:
vector of filters (limited to ugrizy present for the baseline LSST filterset)
returns:
lam_eff_list:
np 1d array of filter effective wavelengths for the filters
alam_over_ebv_list:
np 1d array of alam_over_ebv values for the filters
"""
lam_eff_list = []
alam_over_ebv_list = []
sed_file = os.path.join(os.environ['SIMS_SED_LIBRARY_DIR'],'flatSED','sed_flat.txt.gz')
sed = photUtils.Sed()
sed.readSED_flambda(sed_file)
ax,bx = sed.setupCCM_ab()
ccm_model = 3.1*ax+bx
wl = sed.wavelen
ccm_spline = scipy.interpolate.interp1d(wl,ccm_model,bounds_error=True)
alam_over_ebv = 3.1*ax+bx
for filt in filterset:
if basepath is None:
bp_file = os.path.join(os.environ['LSST_THROUGHPUTS_BASELINE'],'',f'total_{filt}.dat')
else:
bp_file = os.path.join(basepath,f'total_{filt}.dat')
bandpass = photUtils.Bandpass()
bandpass.readThroughput(bp_file)
_,leff = bandpass.calcEffWavelen()
lam_eff_list.append(leff)
alam = ccm_spline(leff)
alam_over_ebv_list.append(alam)
return np.array(lam_eff_list),np.array(alam_over_ebv_list)
Let's do a quick check that we are getting the results that we expect. For DC2 we should get the following for the effective wavelengths and A_lam/E(B-V) values:
u 367.07 nm A_lambda/EBV = 4.812
g 482.69 nm A_lambda/EBV = 3.643
r 622.32 nm A_lambda/EBV = 2.699
i 754.60 nm A_lambda/EBV = 2.063
z 869.01 nm A_lambda/EBV = 1.578
y 971.03 nm A_lambda/EBV = 1.313
If these values do not match those in the next cell, check that the Baseline filter definitions have not changed!
filterlist = ['u','g','r','i','z','y']
leff_list,alamebv_list = compute_alambda_over_ebv(filterlist)
for filt,leff, alamebv in zip(filterlist,leff_list,alamebv_list):
print(f"filter {filt} lam_eff: {leff:.2f}nm alam/E(B-V): {alamebv:.3f}")
filter u lam_eff: 368.48nm alam/E(B-V): 4.802 filter g lam_eff: 480.20nm alam/E(B-V): 3.668 filter r lam_eff: 623.12nm alam/E(B-V): 2.695 filter i lam_eff: 754.17nm alam/E(B-V): 2.065 filter z lam_eff: 869.05nm alam/E(B-V): 1.578 filter y lam_eff: 973.64nm alam/E(B-V): 1.307
These numbers do not match, there are very slight differences, as the filters were updated in mid 2019. The filter curves used for DC2 are released as v1.4 cosmoDC2: https://github.com/lsst/throughputs/releases/tag/1.4. For convenience, we have saved a copy of these filters in the data/dc2_throughputs subdirectory. So, we can calculate the lambda_eff and A_lambda/E(B-V) for these filters:
leff_list,alamebv_list = compute_alambda_over_ebv(filterlist,basepath="data/dc2_throughputs")
for filt,leff, alamebv in zip(filterlist,leff_list,alamebv_list):
print(f"filter {filt} lam_eff: {leff:.2f}nm alam/E(B-V): {alamebv:.3f}")
filter u lam_eff: 367.07nm alam/E(B-V): 4.812 filter g lam_eff: 482.69nm alam/E(B-V): 3.643 filter r lam_eff: 622.32nm alam/E(B-V): 2.699 filter i lam_eff: 754.60nm alam/E(B-V): 2.063 filter z lam_eff: 869.09nm alam/E(B-V): 1.578 filter y lam_eff: 971.03nm alam/E(B-V): 1.313
Yes, these numbers do agree. So, if need A/E(B-V) coefficients for a different filter set, you can use this procedure to derive them.