The purpose of this notebook is to demonstrate how to generate spectra and apply target selection cuts for various mock catalogs and target types. Here we generate spectra for targets in a single healpixel with no constraints on the target density (relative to the expected target density) or contaminants.
For code to generate large numbers of spectra over significant patches of sky and to create a representative DESI dataset (with parallelism), see desitarget/bin/select_mock_targets (as well as its MPI-parallelized cousin, desitarget/bin/select_mock_targets) and desitarget.mock.build.targets_truth.
Finally, note that the various python Classes instantiated here (documented in desitarget.mock.mockmaker) are easily extensible to other mock catalogs and galaxy/QSO/stellar physics. Please contact @desi-data if you have specific suggestions, requirements, or desired features.
John Moustakas
Siena College
2018 September
In [1]:
import os
import sys
import numpy as np
import matplotlib.pyplot as plt
In [2]:
from desiutil.log import get_logger, DEBUG
log = get_logger()
In [3]:
import seaborn as sns
sns.set(style='white', font_scale=1.1, palette='Set2')
In [4]:
%matplotlib inline
In [5]:
healpixel = 26030
nside = 64
In [6]:
seed = 555
rand = np.random.RandomState(seed)
In [7]:
def plot_subset(wave, flux, truth, objtruth, nplot=16, ncol=4, these=None,
xlim=None, loc='right', targname='', objtype=''):
"""Plot a random sampling of spectra."""
nspec, npix = flux.shape
if nspec < nplot:
nplot = nspec
nrow = np.ceil(nplot / ncol).astype('int')
if loc == 'left':
xtxt, ytxt, ha = 0.05, 0.93, 'left'
else:
xtxt, ytxt, ha = 0.93, 0.93, 'right'
if these is None:
these = rand.choice(nspec, nplot, replace=False)
these = np.sort(these)
ww = (wave > 5500) * (wave < 5550)
fig, ax = plt.subplots(nrow, ncol, figsize=(2.5*ncol, 2*nrow), sharey=False, sharex=True)
for thisax, indx in zip(ax.flat, these):
thisax.plot(wave, flux[indx, :] / np.median(flux[indx, ww]))
if objtype == 'STAR' or objtype == 'WD':
thisax.text(xtxt, ytxt, r'$T_{{eff}}$={:.0f} K'.format(objtruth['TEFF'][indx]),
ha=ha, va='top', transform=thisax.transAxes, fontsize=13)
else:
thisax.text(xtxt, ytxt, 'z={:.3f}'.format(truth['TRUEZ'][indx]),
ha=ha, va='top', transform=thisax.transAxes, fontsize=13)
thisax.xaxis.set_major_locator(plt.MaxNLocator(3))
if xlim:
thisax.set_xlim(xlim)
for thisax in ax.flat:
thisax.yaxis.set_ticks([])
thisax.margins(0.2)
fig.suptitle(targname)
fig.subplots_adjust(wspace=0.05, hspace=0.05, top=0.93)
In [8]:
from desitarget.mock.mockmaker import QSOMaker
In [9]:
QSO = QSOMaker(seed=seed)
read methods return a dictionary with (hopefully self-explanatory) target- and mock-specific quantities.Because most mock catalogs only come with (cosmologically accurate) 3D positions (RA, Dec, redshift), we use Gaussian mixture models trained on real data to assign other quantities like shapes, magnitudes, and colors, depending on the target class. For more details see the gmm-dr7.pynb Python notebook.
In [10]:
dir(QSOMaker)
Out[10]:
In [11]:
data = QSO.read(healpixels=healpixel, nside=nside)
In [12]:
for key in sorted(list(data.keys())):
print('{:>20}'.format(key))
In [13]:
%time flux, wave, targets, truth, objtruth = QSO.make_spectra(data)
In [14]:
print(flux.shape, wave.shape)
The truth catalog contains the target-type-agnostic, known properties of each object (including the noiseless photometry), while the objtruth catalog contains different information depending on the type of target.
In [15]:
truth
Out[15]:
In [16]:
objtruth
Out[16]:
In [17]:
QSO.select_targets(targets, truth)
In [18]:
targets
Out[18]:
In [19]:
from desitarget.targetmask import desi_mask
isqso = (targets['DESI_TARGET'] & desi_mask.QSO) != 0
print('Identified {} / {} QSO targets.'.format(np.count_nonzero(isqso), len(targets)))
In [20]:
plot_subset(wave, flux, truth, objtruth, targname='QSO')
In [21]:
from desitarget.mock.mockmaker import LYAMaker
In [22]:
mockfile='/project/projectdirs/desi/mocks/lya_forest/london/v9.0/v9.0.0/master.fits'
In [23]:
LYA = LYAMaker(seed=seed, balprob=0.25)
In [24]:
lyadata = LYA.read(mockfile=mockfile,healpixels=healpixel, nside=nside)
In [25]:
%time lyaflux, lyawave, lyatargets, lyatruth, lyaobjtruth = LYA.make_spectra(lyadata)
In [26]:
lyaobjtruth
Out[26]:
In [27]:
plot_subset(lyawave, lyaflux, lyatruth, lyaobjtruth, xlim=(3500, 5500), targname='LYA')
In [28]:
#Now lets generate the same spectra but including the different features and the new continum model.
#For this we need to reload the desitarget module, for some reason it seems not be enough with defining a diferen variable for the LYAMaker
del sys.modules['desitarget.mock.mockmaker']
In [29]:
from desitarget.mock.mockmaker import LYAMaker
In [30]:
LYA = LYAMaker(seed=seed,sqmodel='lya_simqso_model_develop',balprob=0.25)
In [31]:
lyadata_continum = LYA.read(mockfile=mockfile,healpixels=healpixel, nside=nside)
In [32]:
%time lyaflux_cont, lyawave_cont, lyatargets_cont, lyatruth_cont, lyaobjtruth_cont = LYA.make_spectra(lyadata_continum)
In [33]:
plt.figure(figsize=(20, 10))
indx=rand.choice(len(lyaflux),9)
for i in range(9):
plt.subplot(3, 3, i+1)
plt.plot(lyawave,lyaflux[indx[i]],label="Old Continum")
plt.plot(lyawave_cont,lyaflux_cont[indx[i]],label="New Continum")
plt.legend()
Out[33]:
In [34]:
plt.plot(lyatruth["FLUX_W1"],lyatruth_cont["FLUX_W1"]/lyatruth["FLUX_W1"]-1,'.')
plt.xlabel("FLUX_W1")
plt.ylabel(r"FLUX_W1$^{new}$/FLUX_W1-1")
Out[34]:
In [35]:
plt.plot(lyatruth["FLUX_W2"],lyatruth_cont["FLUX_W2"]/lyatruth["FLUX_W2"]-1,'.')
plt.xlabel("FLUX_W2")
plt.ylabel(r"(FLUX_W2$^{new}$/FLUX_W2)-1")
Out[35]:
In [46]:
plt.hist(lyatruth["FLUX_W1"],bins=100,label="Old Continum",alpha=0.7)
plt.hist(lyatruth_cont["FLUX_W1"],bins=100,label="New Continum",histtype='step',linestyle='--')
plt.xlim(0,100) #Limiting to 100 to see it better.
plt.xlabel("FLUX_W1")
plt.legend()
Out[46]:
In [47]:
plt.hist(lyatruth["FLUX_W2"],bins=100,label="Old Continum",alpha=0.7)
plt.hist(lyatruth_cont["FLUX_W2"],bins=100,label="New Continum",histtype='step',linestyle='--')
plt.xlim(0,100) #Limiting to 100 to see it better.
plt.xlabel("FLUX_W2")
plt.legend()
Out[47]:
In [48]:
del sys.modules['desitarget.mock.mockmaker']
from desitarget.mock.mockmaker import LYAMaker ##Done in order to reload the desitarget, it doesn't seem to be enough with initiating a diferent variable for the LYAMaker class.
In [49]:
LYA = LYAMaker(seed=seed,sqmodel='lya_simqso_model',balprob=0.25,add_dla=True,add_metals="all",add_lyb=True)
lyadata_all= LYA.read(mockfile=mockfile,healpixels=healpixel, nside=nside)
%time lyaflux_all, lyawave_all, lyatargets_all, lyatruth_all, lyaobjtruth_all = LYA.make_spectra(lyadata_all)
In [50]:
plot_subset(lyawave_all, lyaflux_all, lyatruth_all, lyaobjtruth_all, xlim=(3500, 5500), targname='LYA')
In [51]:
def demo_mockmaker(Maker, seed=None, nrand=16, loc='right'):
TARGET = Maker(seed=seed)
log.info('Reading the mock catalog for {}s'.format(TARGET.objtype))
tdata = TARGET.read(healpixels=healpixel, nside=nside)
log.info('Generating {} random spectra.'.format(nrand))
indx = rand.choice(len(tdata['RA']), np.min( (nrand, len(tdata['RA'])) ) )
tflux, twave, ttargets, ttruth, tobjtruth = TARGET.make_spectra(tdata, indx=indx)
log.info('Selecting targets')
TARGET.select_targets(ttargets, ttruth)
plot_subset(twave, tflux, ttruth, tobjtruth, loc=loc,
targname=tdata['TARGET_NAME'], objtype=TARGET.objtype)
In [52]:
from desitarget.mock.mockmaker import LRGMaker
In [53]:
%time demo_mockmaker(LRGMaker, seed=seed, loc='left')
In [54]:
from desitarget.mock.mockmaker import ELGMaker
In [55]:
%time demo_mockmaker(ELGMaker, seed=seed, loc='left')
In [56]:
from desitarget.mock.mockmaker import BGSMaker
In [57]:
%time demo_mockmaker(BGSMaker, seed=seed)
In [58]:
from desitarget.mock.mockmaker import MWS_MAINMaker
In [59]:
%time demo_mockmaker(MWS_MAINMaker, seed=seed, loc='left')
In [60]:
from desitarget.mock.mockmaker import MWS_NEARBYMaker
In [61]:
%time demo_mockmaker(MWS_NEARBYMaker, seed=seed, loc='left')
In [62]:
from desitarget.mock.mockmaker import WDMaker
In [63]:
%time demo_mockmaker(WDMaker, seed=seed, loc='right')
In [64]:
from desitarget.mock.mockmaker import SKYMaker
In [65]:
SKY = SKYMaker(seed=seed)
In [66]:
skydata = SKY.read(healpixels=healpixel, nside=nside)
In [67]:
skyflux, skywave, skytargets, skytruth, objtruth = SKY.make_spectra(skydata)
In [68]:
SKY.select_targets(skytargets, skytruth)
In [ ]: