<jupyter_start><jupyter_text>In this notebook we will simulate Lyman alpha quasar spectra using the DESI code quickqusars  First let us define some useful functions for visualization<jupyter_code>#A simple filter was added to retain some of the noise and check this is reproduced,
#but still be able to spot differences due to features by eye. 
#Plotspec routine was taken https://github.com/desihub/tutorials/blob/master/Intro_to_DESI_spectra.ipynb. 
#It allows to add the true spectra although we will not use this option here 
#These plots will show the three arms of the spectrograph in colors blue ['b'], red ['r'] and yellow ['z'].


def plotspec(spectra, i, truewave=None, trueflux=None, color=['b','r','y'],nfilter=1,label=None):
    plt.plot(spectra.wave['b'], medfilt(spectra.flux['b'][i], nfilter), color=color[0], alpha=0.5,label=label)
    plt.plot(spectra.wave['r'], medfilt(spectra.flux['r'][i], nfilter), color=color[1], alpha=0.5)
    plt.plot(spectra.wave['z'], medfilt(spectra.flux['z'][i], nfilter), color=color[2], alpha=0.5)
    if truewave is not None and trueflux is not None:
        plt.plot(truewave, trueflux[i], 'k-')

    plt.axhline(0, color='k', alpha=0.2)
        
    ymax = np.max(spectra.flux['b'][i])
    ymin=0
        
    plt.ylim(ymin, ymax+3)
    plt.ylabel('Flux')
    plt.xlabel('wavelength [A]')
    plt.legend(loc='best')
    return
    
#This will compare two spectra by dividing the first over the second
def plotspec_ratio(spectra1,spectra2,i,label=None):
    ratio_b=spectra1.flux['b'][i]/spectra2.flux['b'][i]
    ratio_r=spectra1.flux['r'][i]/spectra2.flux['r'][i]
    ratio_z=spectra1.flux['z'][i]/spectra2.flux['z'][i]
    plt.plot(spectra1.wave['b'],ratio_b, color='b', alpha=0.5,label=label)
    plt.plot(spectra1.wave['r'], ratio_r, color='r', alpha=0.5)
    plt.plot(spectra1.wave['z'], ratio_z, color='y', alpha=0.5)
    plt.xlabel('wavelength [A]')
    plt.legend(loc='best')
    
#This will compare two spectra by substracting the second from the first one
def plotspec_diff(spectra1,spectra2,i,label=None):
    ratio_b=spectra1.flux['b'][i]-spectra2.flux['b'][i]
    ratio_r=spectra1.flux['r'][i]-spectra2.flux['r'][i]
    ratio_z=spectra1.flux['z'][i]-spectra2.flux['z'][i]
    plt.plot(spectra1.wave['b'],ratio_b, color='b', alpha=0.5)
    plt.plot(spectra1.wave['r'], ratio_r, color='r', alpha=0.5)
    plt.plot(spectra1.wave['z'], ratio_z, color='y', alpha=0.5)
    plt.xlabel('wavelength [A]')
    if label!=None:
        plt.legend(loc='best')

#This is only to compute the observed wavelenght of given line that corresponds to the quasar redshift. e.g 

def draw_obsline(lbda,z):
    return (1.+z)*lbda<jupyter_output><empty_output><jupyter_text>Let's see what are the parameters needed to run a quickquasars simulation<jupyter_code>!quickquasars --help<jupyter_output><empty_output><jupyter_text>We will start by running quickquasars in a single pixel, for this we define some of the required parameters. We will run it using different exposures (1,2,3,4) and have one with a very large exposure time for comparison. Also this first runs will not have any contaminant, they only have continum LYA and LYB absorption.<jupyter_code>zmin = 1.8 # We use to add it, although it is already set as minimum in the transmission files.
nmax = 10 # maximun number of spectra to simulate in a given pixel, chosen as a small number for speedness in this tutorial.
seed=123 #this ensure reproducibility
path_mocks = '/global/cfs/projectdirs/desi/mocks/lya_forest/london/v9.0' # path to london mocks v9.0
ifile=path_mocks + '/v9.0.0/0/0/transmission-16-0.fits.gz'   #Choose one pixel to simulate the quasars. 
outdir=os.path.join(os.environ['SCRATCH'], 'mocks','quickquasar','1pix_test/')  #Where do you want to save the files.
print(outdir) #Adjuts the path if this is not where you want to store the files produced by this tutorial 
os.makedirs(outdir, exist_ok=True)
#Run quickquasars for the different exposures, store each run a different directory to avoid overwrting.
exposures = ['1','2','3','4','1e5']
contaminant='0.0-' # Following somehow the convention in https://desi.lbl.gov/trac/wiki/LymanAlphaWG/LondonMocks/v9_quickquasars

for exp in exposures:
    if not os.path.exists(outdir +contaminant+exp):
        os.mkdir(outdir +contaminant+exp )
    specfile = outdir +contaminant+exp+'/spectra.fits'
    cmd ='quickquasars --exptime {exptime} -i {ifile} -o {outfile}  --zmin {zmin} \
    --zbest --bbflux --add-LYB --desi-footprint \
    --nmax {nmax} --seed {seed}'.format(exptime=float(exp)*1000,ifile=ifile,outfile=specfile,zmin=zmin,nmax=nmax,seed=seed)
    !time $cmd<jupyter_output><empty_output><jupyter_text>Now lets comapare the spectra we got from each simulation<jupyter_code>dump = []
for exp in exposures:
    dump.append(outdir +'0.0-'+exp)
path_highexp = dump[-1]
del(dump[-1])
dump.sort()
print(path_highexp)
spec = fitsio.FITS(path_highexp+'/spectra.fits')
spec
truth = fitsio.FITS(path_highexp + '/truth-16-0.fits')
spectra_highexp = desispec.io.read_spectra(path_highexp+'/spectra.fits')
index = 4
plt.figure(figsize=(10,5),dpi=100)
for i in range(len(dump)):
    plt.subplot(221+i)
    exposure = dump[i].split('/0.0-')[-1]
    spectra= desispec.io.read_spectra(dump[i]+'/spectra.fits')
    plotspec(spectra,index,label='Exposure = '+exposure)
    plotspec(spectra_highexp,index,color=['k','k','k'],label='Noisless')
    z_obj = truth[1]['Z'][index]
    lya = draw_obsline(1215.67,z_obj)
    lyb = draw_obsline(1025.722,z_obj)
    civ = draw_obsline(1549.48,z_obj)
    #mgii = draw_obsline(2799.117,z_obj)
    plt.axvline(lya,color='g')
    plt.axvline(lyb,color='g')
    plt.axvline(civ,color='g')
    #plt.axvline(mgii,color='g')
    #plt.xlim(3500,5000)
plt.subplots_adjust(hspace=0.3,wspace=.2)<jupyter_output><empty_output><jupyter_text>We can compare the differences as well<jupyter_code>index = 4
plt.figure(figsize=(10,5),dpi=100)
for i in range(len(dump)):
    plt.subplot(221+i)
    exposure = dump[i].split('/')[-1]
    spectra= desispec.io.read_spectra(dump[i]+'/spectra.fits')
    plotspec_ratio(spectra,spectra_highexp,index,label='Exposure= '+exposure)
    #plt.xlim(3500,4300)
    plt.ylim(-5,5)
plt.subplots_adjust(hspace=0.3,wspace=.2)<jupyter_output><empty_output><jupyter_text>Lets look at the the other files that were produced: zbest.fits and truth.fits  The "zbest" files are fake redshift catalogs for the simulated spectra. By default it is conformed with the true redshift but we can make it more realistic by adding some errors with the following flags:--sigma_kms_fog--gamma_kms_zfit--shift_kms_los<jupyter_code>data_zbest = fitsio.FITS(path_highexp+'/zbest-16-0.fits')
data_zbest[1]<jupyter_output><empty_output><jupyter_text>The "Truth" files stores all the information regarding the templates used to construct the spectra, as well as the magnitude of the objects, true redshift, if the spectra has contaminats, etc.<jupyter_code>data_truth = fitsio.FITS(path_highexp+'/truth-16-0.fits')
data_truth<jupyter_output><empty_output><jupyter_text>So far the spectra has no contaminants, so lets add some. We will start with DLA's   There is two ways to add DLA's, one is at random redshifts and column density, and the other is by what is already stored in the transmission file. In this tutorial we will use the DLAs specified in the transmission file<jupyter_code>exp='4'
contaminant='0.1DLA-'
outdir_dla=outdir +contaminant+exp
os.makedirs(outdir +contaminant+exp,exist_ok=True)
        
specfile = outdir +contaminant + exp +'/spectra.fits'
cmd ='quickquasars  -i {ifile} -o {outfile}  --zmin {zmin} --exptime {exptime}\
    --zbest --bbflux --add-LYB --desi-footprint \
    --nmax {nmax} --seed {seed} --dla file'.format(exptime=float(exp)*1000,ifile=ifile,outfile=specfile,zmin=zmin,nmax=nmax,seed=seed)

!$cmd
spectra_dla = desispec.io.read_spectra(outdir_dla +'/spectra.fits')
spectra_4exp = desispec.io.read_spectra(outdir + '0.0-4/spectra.fits')
truth_dla = fitsio.FITS(outdir_dla+'/truth-16-0.fits')
truth_dla[3]['TARGETID'][:]
mask = np.isin(truth_dla[1]['TARGETID'][:],truth_dla[3]['TARGETID'][:])
index = np.arange(nmax)
index = index[mask]
plt.figure(figsize=(10,5),dpi=100)
for i in range(3):
    plt.subplot(221+i)
    plotspec(spectra_4exp,index[i],color=['b','b','b'],label='no DLA')
    plotspec(spectra_dla,index[i],color=['k','k','k'],label=' with DLA')
    plt.xlim(3500,6000) #Make a zoom to see the differences
plt.subplots_adjust(hspace=0.3,wspace=.2)
plt.figure(figsize=(10,5),dpi=100)
for i in range(3):
    plt.subplot(221+i)
    plotspec_diff(spectra_dla,spectra_4exp,index[i])
    plt.xlim(3500,4000)
plt.subplots_adjust(hspace=0.3,wspace=.2)<jupyter_output><empty_output><jupyter_text>Next, lets add BALs<jupyter_code>exp='4'
contaminant='0.1BAL-'
outdir_bal=outdir +contaminant+exp
os.makedirs(outdir_bal,exist_ok=True)

specfile = outdir +contaminant + exp +'/spectra.fits'
cmd ='quickquasars  -i {ifile} -o {outfile}  --zmin {zmin} --exptime {exptime}\
    --zbest --bbflux --add-LYB --desi-footprint \
    --nmax {nmax} --seed {seed} --balprob 0.5'.format(exptime=float(exp)*1000,ifile=ifile,outfile=specfile,zmin=zmin,nmax=nmax,seed=seed)
!$cmd
truth_bal = fitsio.FITS(outdir_bal+'/truth-16-0.fits')
mask = np.isin(truth_bal[1]['TARGETID'][:],truth_bal[3]['TARGETID'][:])
index = np.arange(nmax)
index = index[mask]
spectra_bal = desispec.io.read_spectra(outdir_bal +'/spectra.fits')
plt.figure(figsize=(10,5),dpi=100)
for i in range(4):
    plt.subplot(221+i)
    plotspec(spectra_4exp,index[i],color=['b','b','b'],label='no BALs')
    plotspec(spectra_bal,index[i],color=['k','k','k'],label=' with BALs')
    plt.xlim(3500,6000)
plt.subplots_adjust(hspace=0.3,wspace=.2)
plt.figure(figsize=(10,5),dpi=100)
for i in range(4):
    plt.subplot(221+i)
    plotspec_diff(spectra_bal,spectra_4exp,index[i])
    plt.xlim(3500,6000)
plt.subplots_adjust(hspace=0.3,wspace=.2)<jupyter_output><empty_output><jupyter_text>Finally, we can also add metals<jupyter_code>exp='4'
contaminant='0.1MET-'
outdir_metal=outdir +contaminant+exp
os.makedirs(outdir_metal,exist_ok=True)
specfile = outdir_metal+'/spectra.fits'
cmd ='quickquasars --exptime {exptime} -i {ifile} -o {outfile}  --zmin {zmin} \
    --zbest --bbflux --add-LYB --desi-footprint \
     --nmax {nmax} --seed {seed} --metals-from-file'.format(exptime=float(exp)*1000,ifile=ifile,outfile=specfile,zmin=zmin,nmax=nmax,seed=seed)
!time $cmd
spectra_metals = desispec.io.read_spectra(outdir_metal +'/spectra.fits')
plt.figure(figsize=(10,5),dpi=100)
for i in range(4):
    plt.subplot(221+i)
    plotspec(spectra_4exp,i,color=['b','b','b'],label='no metals')
    plotspec(spectra_metals,i,color=['k','k','k'],label=' with metals')
    plt.xlim(3500,4500)
plt.subplots_adjust(hspace=0.3,wspace=.2)
plt.figure(figsize=(10,5),dpi=100)
for i in range(4):
    plt.subplot(221+i)
    plotspec_ratio(spectra_metals,spectra_4exp,i)
    plt.ylabel('Spectra ratio')
    plt.xlim(3500,4500)
    #plt.ylim(0.8,1.01)
plt.subplots_adjust(hspace=0.3,wspace=.2)<jupyter_output><empty_output><jupyter_text>As stated at the begining, all this has been for a sinlge pixel but we can also run a quickquasar simulation for the full DESI footprint. For this, lets take a look and run to the script quickquasars.sh   Now that we have a bigger quickquasar run, we can stack the outputs for each pixel into a sinlge catalog<jupyter_code>mock_path = os.path.join(os.environ['SCRATCH'], 'mocks','quickquasar')
out_zcat = mock_path+'/tutorial-0.3-4/'
cmd='desi_zcatalog -i {ipath} -o {outfile} --fibermap'.format(ipath=out_zcat+'/spectra-16/',outfile=out_zcat+'zcat.fits')
if not os.path.isfile(out_zcat+'zcat.fits'):
    !time $cmd
else:
    print('zcat.fits already exist!')
zcat = fitsio.FITS(out_zcat+'/zcat.fits')
zcat
zcat[1]<jupyter_output><empty_output><jupyter_text>To finish this part of the tutorial create some QA plots  We can look at the redshift distribution of the zcatalog and compare to the expected one.<jupyter_code>## Retrive the expected distribution as function of redshift and magnitude, and define some functions useful for plotting them. Some code was borrowed/adapted from https://github.com/dkirkby/ArgonneLymanAlpha/blob/master/notebooks/SampleProperties.ipynb

from astropy.table import Table,Column

def bin_index(bin_centers, low_edge):
    """Find the index of the bin with the specified low edge, where bins is an array of equally-spaced bin centers.
    """
    delta = bin_centers[1] - bin_centers[0]
    min_value = bin_centers[0] - 0.5 * delta
    index = int(round((low_edge - min_value) / delta))
    if abs((low_edge - min_value) / delta - index) > 1e-5:
        raise ValueError('low_edge = {} is not aligned with specified bins.'.format(low_edge))
    return index

def marginal_mag(magbin, zbin, nqso, z_min=0.):
    """Tabulate the marginal distribution of QSOs per sq.deg. per (g or r) magnitude bin with z > z_min.
    """
    cut_index = bin_index(zbin, z_min)
    return magbin, np.sum(nqso[:, cut_index:], axis=1)

def marginal_z(magbin, zbin, nqso, mag_max):
    """Tabulate the marginal distribution of QSOs per sq.deg. per redshift bin with mag < mag_max.
    """
    cut_index = bin_index(magbin, mag_max)
    return zbin, np.sum(nqso[:cut_index], axis=0)

qlf_file=os.getenv("DESIMODEL")+'/data/targets/nzr_qso.dat'
qlftab=Table.read(qlf_file,format='ascii',names=('z','r','dNdzdr'))
zbin=np.unique(qlftab['z'].data)
rbin=np.unique(qlftab['r'].data)
nz,nr=len(zbin),len(rbin)
assert len(qlftab) == nz*nr
nqso_=qlftab['dNdzdr'].data.reshape(nz,nr).transpose()

#Plot the redshift distribution from the z catalog
z=zcat[1]['Z'][:]
bins = np.arange(min(z),max(z),0.2)
nz_,_,_=plt.hist(z,bins=bins,histtype='step',linewidth=2.);
plt.xlabel(r'$Z$');
plt.ylabel('number');
zz, nz = marginal_z(rbin, zbin, nqso_, mag_max=23)
plt.plot(zz[zz>1.8],nz[zz>1.8]/nz[zz>1.8].sum(),label='DESI')
plt.hist(z,bins=bins,weights=np.ones(len(z))/nz_.sum(),histtype='step',linewidth=2.,label='mock');
plt.xlabel(r'$Z$');
plt.ylabel(r'$n(Z)$ (normalized)');
plt.legend()<jupyter_output><empty_output><jupyter_text>Also we can look at the magnitude of each of the arms [b,r,z]. As the zcatalog does not hace the magnitude information, we can calculate it from the flux as:  $$ m = 22.5 - 2.5 \log_{10}(f) $$<jupyter_code>def magnitud(flux):
    return 22.5-2.5*np.log10(flux)
plt.figure(figsize=(15,3),dpi=100)

plt.subplot(131)
mag_b = magnitud(zcat[1]['FLUX_G'][:])
bins = np.arange(min(mag_b),max(mag_b),0.5)
plt.hist(mag_b,bins=bins,histtype='step',linewidth=1,label='g');
plt.xlabel('magnitude');
plt.ylabel('#');
plt.legend(loc='upper left')

plt.subplot(132)
mag_r = magnitud(zcat[1]['FLUX_R'][:])
bins = np.arange(min(mag_r),max(mag_r),0.5)
mr,binr,_=plt.hist(mag_r,bins=bins,histtype='step',linewidth=1,label='r');
plt.xlabel('magnitude');
plt.ylabel('#');
plt.legend(loc='upper left')

plt.subplot(133)
mag_z = magnitud(zcat[1]['FLUX_Z'][:])
bins = np.arange(min(mag_z),max(mag_z),0.5)
plt.hist(mag_z,bins=bins,histtype='step',linewidth=1,label='z');
plt.xlabel('magnitude');
plt.ylabel('#');
plt.legend(loc='upper left')


plt.subplots_adjust(wspace=.3)
#Using the function above lets compare the magnitude r to the expected one. 

r, nr = marginal_mag(rbin, zbin, nqso_, z_min=1.8)
plt.plot(r,nr/nr.sum(),label='DESI')
plt.hist(mag_r,bins=bins,weights=np.ones(len(mag_r))/mr.sum(),histtype='step',linewidth=1,label='mock ');
plt.xlabel('mag-r')
plt.ylabel('n(mag_r) (normalized)')
plt.legend()<jupyter_output><empty_output><jupyter_text>We can look at the distriubution of objects in the sky<jupyter_code>plt.scatter(zcat[1]['RA'][:],zcat[1]['DEC'][:],marker='.');
plt.xlabel('RA [deg]');
plt.ylabel('DEC [deg]');<jupyter_output><empty_output><jupyter_text>Next we can look at the object distributions for each pixels. For that we define functions to go from angles to pixels and back<jupyter_code>def get_pix(ra,dec,nside=16):
    phi = ra*np.pi / 180.
    th = np.pi/2. - dec*np.pi/180.
    pix = healpy.ang2pix(nside,th,phi,nest=True)
    return pix

def pix2ang_deg(nside,pix,nest=True):
    pix_rad_dec,pix_rad_ra = healpy.pix2ang(nside,pix,nest=nest)
    pix_ra = pix_rad_ra * 180. / np.pi
    pix_dec = (np.pi/2.-pix_rad_dec)*180./np.pi
    return pix_ra, pix_dec
pix_zcat = get_pix(zcat[1]['RA'][:],zcat[1]['DEC'][:])
ra_zcat,dec_zcat = pix2ang_deg(16,np.unique(pix_zcat),nest=True)<jupyter_output><empty_output><jupyter_text>Here we are just ploting the pixels, in the previous plot we plotted all the objects<jupyter_code>plt.scatter(ra_zcat,dec_zcat,marker='.');
plt.xlabel('RA [deg]');
plt.ylabel('DEC [deg]');<jupyter_output><empty_output><jupyter_text>Now lets plot the density of objects per pixels as well at the mean redshift in each pixel<jupyter_code>def get_density(pix,nside):
    bincounts_pix = np.bincount(pix)
    unique_pix = np.unique(pix)
    area_pix = healpy.pixelfunc.nside2pixarea(nside, degrees=True)
    density = bincounts_pix[unique_pix]/area_pix
    return density
plot1 = plt.scatter(ra_zcat,dec_zcat,c=get_density(pix_zcat,16),marker='.')
cbar1 = plt.colorbar(plot1)
cbar1.set_label(r'$Density \, [\# \, deg^-2]$')
cbar1.update_ticks()
plt.xlabel('RA [deg]')
plt.ylabel('DEC [deg]')
bincounts_pix = np.bincount(pix_zcat)
bincounts_z = np.bincount(pix_zcat,weights=zcat[1]['Z'][:])

z_pix_zcat = bincounts_z[np.unique(pix_zcat)]/bincounts_pix[np.unique(pix_zcat)]

plot1 = plt.scatter(ra_zcat,dec_zcat,c=z_pix_zcat,marker='.')
cbar1 = plt.colorbar(plot1)
cbar1.set_label(r'$<z>$',size=20)
cbar1.update_ticks()
plt.xlabel('RA [deg]')
plt.ylabel('DEC [deg]')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.16176220445654022, 'fraction_numerical': 0.02653512502126212, 'mean_word_length': 5.6183864915572235, 'pattern_counts': {'":': 0, '<': 71, '<?xml version=': 0, '>': 75, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '17', 'dir': 'jupyter-structured-clean-dedup', 'id': '250850', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 7334, 'n_tokens_neox': 6888, 'n_words': 1365, 'path': 'Dec2020/quickquasars/quickquasars.ipynb', 'repo_name': 'humnaawan/tutorials', 'size': '17637'}