<jupyter_start><jupyter_text>Set where to write Lya ensemble dflux.<jupyter_code>outdir = '/global/homes/m/mjwilson/sandbox/lya-signal/desimodel/0.14.0/data/tsnr'<jupyter_output><empty_output><jupyter_text>Load a set of Vi'd Lya QSOs.<jupyter_code>dat = fits.open('/project/projectdirs/desi/spectro/redux/cascades/tiles/80609/deep/coadd-0-80609-deep.fits')
dat.info()
vi = pandas.read_csv('/project/projectdirs/desi/sv/vi/TruthTables/Blanc/QSO/desi-vi_QSO_tile80609_nightdeep_merged_all_210210_ADDING_object_info.csv')
vi
isin = (vi['best_spectype'] == 'QSO') & (vi['best_quality'] >= 2.5) & (vi['best_z'] >= 2.1)
vi = vi[isin]
vi
tids = vi['TARGETID']
gauss_kernel = Gaussian1DKernel(15)
isin = np.isin(dat['FIBERMAP'].data['TARGETID'], tids)
fmap_ids = dat['FIBERMAP'].data['TARGETID'][isin]
nin = np.count_nonzero(fmap_ids)
gmags = 22.5 - 2.5*np.log10(dat['FIBERMAP'].data['FLUX_G'][isin] / mwdust_transmission(dat['FIBERMAP'].data['EBV'][isin], 'G', dat['FIBERMAP'].data['PHOTSYS'][isin]))
gmags<jupyter_output><empty_output><jupyter_text>Our QSOs<jupyter_code>fig, axes = plt.subplots(nin, 1, figsize=(5, 5 * nin))

for band in ['B','R','Z']:
    for i, x in enumerate(dat['{}_FLUX'.format(band)].data[isin]):
        axes[i].plot(dat['{}_WAVELENGTH'.format(band)].data,  convolve(x, gauss_kernel), lw=0.5)
        axes[i].set_ylim(bottom=-0.5)<jupyter_output><empty_output><jupyter_text>Take closest to g=22 to be our reference.<jupyter_code>idx = np.where(np.abs(gmags - 22.) == np.abs(gmags - 22.).min())[0][0]
# Force 7 
idx = 7 
# Closest to 22.
master_fluxes = {'gmag': gmags[idx], 'tid': fmap_ids[idx]}

for band in ['B', 'R', 'Z']:
    master_fluxes[band] = {'wave':                dat['{}_WAVELENGTH'.format(band)].data,
                           'smoothflux': convolve(dat['{}_FLUX'.format(band)].data[isin][idx], gauss_kernel),
                           'ivar':                dat['{}_IVAR'.format(band)].data[isin][idx]}
master_fluxes['tid']
master_fluxes['gmag']
vi[vi['TARGETID'] == master_fluxes['tid']]
master_fluxes['z'] = vi[vi['TARGETID'] == master_fluxes['tid']]['best_z']
master_fluxes['continuum'] = 0.43
pl.plot(master_fluxes['B']['wave'], master_fluxes['B']['smoothflux'])
pl.axhline(master_fluxes['continuum'], c='k', lw=0.5)

pl.xlabel('Wavelength [A]')
pl.ylabel('1.e-17 ergs/s/cm2/A')<jupyter_output><empty_output><jupyter_text>Later we use this (by eye) 'continuum' as our asymptotic 'signal' normalization at the blue end.  Get a QSO n(z)<jupyter_code># https://desi.lbl.gov/svn/code/desimodel/tags/0.14.0/data/targets/nz_qso.dat; 
# Number per sq. deg. per dz=0.1
# Note: Cascades
zlo, zhi, Nz = np.loadtxt('/global/common/software/desi/cori/desiconda/20200801-1.4.0-spec/code/desimodel/0.14.0/data/targets/nz_qso.dat', unpack=True)
zmid = 0.5 * (zlo + zhi)
Nz /= Nz.max()
pl.plot(zmid, Nz, c='k', lw=0.5)
pl.xlabel('z')
zs = np.random.uniform(0.0, 5.0, 500000)
zs = np.sort(zs)
# pl.hist(zs, bins=np.arange(0.0, 5.0, 0.1))
draws = np.random.uniform(0.0, 1.0, 500000)
idx = np.digitize(zs, bins=np.arange(0.0, 5.1, 0.1))
probs = np.zeros_like(idx, dtype=np.float)
for i, uid in enumerate(np.unique(idx)[:-1]):
    probs[idx == uid] = Nz[i]
draws
probs
isin = draws <= probs
qso_zs = zs[isin]<jupyter_output><empty_output><jupyter_text>Here we've drawn an ensemble of zs from this distribution.<jupyter_code>pl.plot(zmid, 5000. * Nz, c='k', lw=0.5)
pl.hist(qso_zs, bins=np.arange(0.0, 5.0, 0.05), alpha=0.5)
pl.xlabel('z')
lya_zs = qso_zs[qso_zs > 2.1]
lya_zs 
# lya_zs = lya_zs[:2]
# 1216. * (1. + lya_zs)
nlya = len(lya_zs)<jupyter_output><empty_output><jupyter_text>Our 'signal' will be unity bluer than Lya for a given redshift (zero otherwise).  We then stack across the ensemble.<jupyter_code>tracer        = 'LYA'

hdr           = fits.Header()
hdr['NMODEL'] = nlya
hdr['TRACER'] = tracer
hdr['FILTER'] = 'decam2014-g'
hdr['ZLO']    = 2.1
        
hdu_list      = [fits.PrimaryHDU(header=hdr)]

for band in ['b', 'r', 'z']:
    wave   = dat['{}_WAVELENGTH'.format(band)].data
    nwave  = wave[:,None] * np.ones(nlya, dtype=float)[None,:]
    
    weight = np.zeros(shape=(len(wave), nlya), dtype=float)
    
    for i, z in zip(range(nlya), lya_zs):
        weight[nwave[:,i] < (1. + z) * 1216., i] = 1.0

    mweight = np.mean(weight, axis=1)    
    
    zpivot  = 2.4 
    zfactor = (wave / (1. + zpivot) / 1216.)**0.95    
    zweight = zfactor * mweight
 
    mweight = np.expand_dims(master_fluxes['continuum'] * mweight, axis=0)
    zweight = np.expand_dims(master_fluxes['continuum'] * zweight, axis=0)

    if band =='b':
        pl.plot(wave, mweight[0], c='k', linestyle='--', label='No z weight')
        pl.plot(wave, zweight[0], c='k', label='z weight')
        
    else:
        pl.plot(wave, mweight[0], c='k', linestyle='--', label='')
        pl.plot(wave, zweight[0], c='k', label='')
        
    hdu_list.append(fits.ImageHDU(wave, name='WAVE_{}'.format(band.upper())))
    hdu_list.append(fits.ImageHDU(zweight, name='DFLUX_{}'.format(band.upper())))

hdu_list = fits.HDUList(hdu_list)
hdu_list.writeto('{}/tsnr-ensemble-{}.fits'.format(outdir, tracer.lower()), overwrite=True)

pl.xlabel('Wavelength [A]')
pl.ylabel('1.e-17 ergs/s/cm2/A')
pl.legend(frameon=False, loc=1)

print('Written to {}/tsnr-ensemble-{}.fits'.format(outdir, tracer.lower()))<jupyter_output>Written to /global/homes/m/mjwilson/sandbox/lya-signal/desimodel/0.14.0/data/tsnr/tsnr-ensemble-lya.fits<jupyter_text>Finally, here we've used our reference continuum from above as the blue end normalization and write to disk at outdir.  Check against QSO tsnr.<jupyter_code>ens = fits.open('/global/common/software/desi/cori/desiconda/20200801-1.4.0-spec/code/desimodel/0.14.0/data/tsnr/tsnr-ensemble-qso.fits')
ens.info()
ens['DFLUX_B'].shape<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.19720978298312092, 'fraction_numerical': 0.03496383052015157, 'mean_word_length': 4.452582159624413, 'pattern_counts': {'":': 0, '<': 34, '<?xml version=': 0, '>': 35, 'https://': 1, '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': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '649563', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 2631, 'n_tokens_neox': 2461, 'n_words': 387, 'path': 'doc/nb/Lya-tsnr-signal.ipynb', 'repo_name': 'echaussidon/desispec', 'size': '5806'}