<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Connecting Spectra to Mocks
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
Step1:  To keep the calculations below manageable we specify a single nside=64 healpixel in an arbitrary location of the DESI footprint.
Step2:  Specifying the random seed makes our calculations reproducible.
Step4:  Define a couple wrapper routines we will use below several times.
Step5:  Tracer QSOs
Both tracer and Lya QSO spectra contain an underlying QSO spectrum, but the Lya QSOs (which we demonstrate below) also include the Lya forest (here, based on the v2.0 of the "London" mocks).
Every target class has its own dedicated "Maker" class.
Step6:  The various 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.
Step7:  Now we can generate the spectra as well as the targeting catalogs (targets) and corresponding truth table.
Step8:  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.
Step9:  Next, let's run target selection, after which point the targets catalog should look just like an imaging targeting catalog (here, using the DR7 data model).
Step10:  And indeed, we can see that only a subset of the QSOs were identified as targets (the rest scattered out of the QSO color selection boxes).
Step11:  Finally, let's plot some example spectra.
Step12:  Generating QSO spectra with cosmological Lya skewers proceeds along similar lines.
Here, we also include BALs with 25% probability.
Step13:  Lets plot together some of the spectra with the old and new continum model
Step14:  And finally we compare the colors, for the two runs with the new and old continum
Step15:  Conclusion
Step16:  Demonstrate the other extragalactic target classes
Step17:  LRGs
Step18:  ELGs
Step19:  BGS
Step20:  Next, demonstrate how to generate spectra of stars...
MWS_MAIN
Step21:  MWS_NEARBY
Step22:  White dwarfs (WDs)
Step23:  Finally demonstrate how to generate (empyt) SKY spectra.
<ASSISTANT_TASK:>
Python Code:
import os
import sys
import numpy as np
import matplotlib.pyplot as plt
from desiutil.log import get_logger, DEBUG
log = get_logger()
import seaborn as sns
sns.set(style='white', font_scale=1.1, palette='Set2')
%matplotlib inline
Explanation: Connecting Spectra to Mocks
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
End of explanation
healpixel = 26030
nside = 64
Explanation: To keep the calculations below manageable we specify a single nside=64 healpixel in an arbitrary location of the DESI footprint.
End of explanation
seed = 555
rand = np.random.RandomState(seed)
Explanation: Specifying the random seed makes our calculations reproducible.
End of explanation
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)
Explanation: Define a couple wrapper routines we will use below several times.
End of explanation
from desitarget.mock.mockmaker import QSOMaker
QSO = QSOMaker(seed=seed)
Explanation: Tracer QSOs
Both tracer and Lya QSO spectra contain an underlying QSO spectrum, but the Lya QSOs (which we demonstrate below) also include the Lya forest (here, based on the v2.0 of the "London" mocks).
Every target class has its own dedicated "Maker" class.
End of explanation
dir(QSOMaker)
data = QSO.read(healpixels=healpixel, nside=nside)
for key in sorted(list(data.keys())):
    print('{:>20}'.format(key))
Explanation: The various 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.
End of explanation
%time flux, wave, targets, truth, objtruth = QSO.make_spectra(data)
print(flux.shape, wave.shape)
Explanation: Now we can generate the spectra as well as the targeting catalogs (targets) and corresponding truth table.
End of explanation
truth
objtruth
Explanation: 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.
End of explanation
QSO.select_targets(targets, truth)
targets
Explanation: Next, let's run target selection, after which point the targets catalog should look just like an imaging targeting catalog (here, using the DR7 data model).
End of explanation
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)))
Explanation: And indeed, we can see that only a subset of the QSOs were identified as targets (the rest scattered out of the QSO color selection boxes).
End of explanation
plot_subset(wave, flux, truth, objtruth, targname='QSO')
Explanation: Finally, let's plot some example spectra.
End of explanation
from desitarget.mock.mockmaker import LYAMaker
mockfile='/project/projectdirs/desi/mocks/lya_forest/london/v9.0/v9.0.0/master.fits'
LYA = LYAMaker(seed=seed, balprob=0.25)
lyadata = LYA.read(mockfile=mockfile,healpixels=healpixel, nside=nside)
%time lyaflux, lyawave, lyatargets, lyatruth, lyaobjtruth = LYA.make_spectra(lyadata)
lyaobjtruth
plot_subset(lyawave, lyaflux, lyatruth, lyaobjtruth, xlim=(3500, 5500), targname='LYA')
#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']
from desitarget.mock.mockmaker import LYAMaker
LYA = LYAMaker(seed=seed,sqmodel='lya_simqso_model_develop',balprob=0.25)
lyadata_continum = LYA.read(mockfile=mockfile,healpixels=healpixel, nside=nside)
%time lyaflux_cont, lyawave_cont, lyatargets_cont, lyatruth_cont, lyaobjtruth_cont = LYA.make_spectra(lyadata_continum)
Explanation: Generating QSO spectra with cosmological Lya skewers proceeds along similar lines.
Here, we also include BALs with 25% probability.
End of explanation
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()
Explanation: Lets plot together some of the spectra with the old and new continum model
End of explanation
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")
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")
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()
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()
Explanation: And finally we compare the colors, for the two runs with the new and old continum
End of explanation
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. 
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)
plot_subset(lyawave_all, lyaflux_all, lyatruth_all, lyaobjtruth_all, xlim=(3500, 5500), targname='LYA')
Explanation: Conclusion: Colors are slightly affected by changing the continum model.
To Finalize the LYA section, lets generate another set of spectra now including DLAs, metals, LYB, etc.
End of explanation
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)
Explanation: Demonstrate the other extragalactic target classes: LRG, ELG, and BGS.
For simplicity let's write a little wrapper script that does all the key steps.
End of explanation
from desitarget.mock.mockmaker import LRGMaker
%time demo_mockmaker(LRGMaker, seed=seed, loc='left')
Explanation: LRGs
End of explanation
from desitarget.mock.mockmaker import ELGMaker
%time demo_mockmaker(ELGMaker, seed=seed, loc='left')
Explanation: ELGs
End of explanation
from desitarget.mock.mockmaker import BGSMaker
%time demo_mockmaker(BGSMaker, seed=seed)
Explanation: BGS
End of explanation
from desitarget.mock.mockmaker import MWS_MAINMaker
%time demo_mockmaker(MWS_MAINMaker, seed=seed, loc='left')
Explanation: Next, demonstrate how to generate spectra of stars...
MWS_MAIN
End of explanation
from desitarget.mock.mockmaker import MWS_NEARBYMaker
%time demo_mockmaker(MWS_NEARBYMaker, seed=seed, loc='left')
Explanation: MWS_NEARBY
End of explanation
from desitarget.mock.mockmaker import WDMaker
%time demo_mockmaker(WDMaker, seed=seed, loc='right')
Explanation: White dwarfs (WDs)
End of explanation
from desitarget.mock.mockmaker import SKYMaker
SKY = SKYMaker(seed=seed)
skydata = SKY.read(healpixels=healpixel, nside=nside)
skyflux, skywave, skytargets, skytruth, objtruth = SKY.make_spectra(skydata)
SKY.select_targets(skytargets, skytruth)
Explanation: Finally demonstrate how to generate (empyt) SKY spectra.
End of explanation
<END_TASK>