quickcat calibration

This notebook is the quickcat calibration script.

Its input is a redshift catalog merged with a target list and a truth table from simulations.
Its output is a set of coefficients saved in a yaml file to be copied to desisim/py/desisim/data/quickcat.yaml .

This notebook does sequencially :

open the merged redshift catalog
run quickcat on it
for each target class
- fit a model for redshift efficiency, and display input, quickcat, best fit model
- fit a model for redshift uncertainty, and display input, (quickcat,) best fit model
save the best fit parameters in a yaml file (to be copied in desisim/data/quickcat.yaml)

Please first edit first the following path:



In [1]:

    
# input merged catalog (from simulations for now)
simulation_catalog_filename="/home/guy/Projets/DESI/analysis/quickcat/20180926/zcatalog-redwood-target-truth.fits"
# output quickcat parameter file that this code will write
quickcat_param_filename="/home/guy/Projets/DESI/analysis/quickcat/20180926/quickcat.yaml"
# output quickcat catalog (same input target and truth)
quickcat_catalog_filename="/home/guy/Projets/DESI/analysis/quickcat/20180926/zcatalog-redwood-target-truth-quickcat.fits"



In [2]:

    
import os.path
import numpy as np
import matplotlib.pyplot as plt
import astropy.io.fits as pyfits
import scipy.optimize
from pkg_resources import resource_filename
import yaml
from desisim.quickcat import eff_model,get_observed_redshifts
from desisim.simexp import reference_conditions

def eff_err(k,n) :
    # given k and n
    # the most probable efficiency is k/n but the uncertainty is complicated
    # I choose to define error bar as FWHM/2.35 , converging to sigma for large k,n-k,and n
    # this is the Bayesian probability
    # P(eff|k,n) = gamma(n+2)/(gamma(k+1)*gamma(n-k+1)) * eff**k*(1-eff)**(n-k)
    
    if k>10 and n-k>10 and n>10 :
        return np.sqrt(k*(1-k/n))/n
    
    ns=300
    e=np.arange(ns)/ns
    p=e**k*(1-e)**(n-k)
    xc=float(k)/n
    i=int(ns*xc+0.5)
    if i>ns-1 : i=ns-1
    p/=p[i]
    if k==0 :
        xl=0
    else :
        xl = np.interp(0.5*p[i],p[:i],e[:i])
    if k==n :
        xh=1
    else :
        xh = np.interp(0.5*p[i],p[i:][::-1],e[i:][::-1])
    sigma =  (xh-xl)/2.35
    return sigma

def efficiency(x,selection,bins=40) :
    h0,bins=np.histogram(x,bins=bins)
    hx,bins=np.histogram(x,bins=bins,weights=x)
    h1,bins=np.histogram(x[selection],bins=bins)
    ii=(h0>1)
    n=h0[ii]
    k=h1[ii]
    meanx=hx[ii]/n
    eff=k/n
    err=np.zeros(eff.shape)
    for i in range(err.size) :
        err[i] = eff_err(k[i],n[i])
    return meanx,eff,err

def prof(x,y,bins=40) :
    h0,bins=np.histogram(x,bins=bins)
    hx,bins=np.histogram(x,bins=bins,weights=x)
    hy,bins=np.histogram(x,bins=bins,weights=y)
    hy2,bins=np.histogram(x,bins=bins,weights=y**2)
    ii=(h0>1)
    n=h0[ii]
    x=hx[ii]/n
    y=hy[ii]/n
    y2=hy2[ii]/n
    var=y2-y**2
    err=np.zeros(x.size)
    err[var>0]=np.sqrt(var[var>0])
    return x,y,err,n

def efficiency2d(x,y,selection,bins=20) :
    h0,xx,yy=np.histogram2d(x,y,bins=bins)
    h1,xx,yy=np.histogram2d(x[selection],y[selection],bins=(xx,yy))
    shape=h0.shape
    n=h0.ravel()
    k=h1.ravel()
    eff=np.zeros(n.size)
    err=np.zeros(n.size)
    for i in range(n.size) :
        if n[i]==0 :
            err[i]=1000.
        else :
            eff[i]=k[i]/n[i]
            err[i]=eff_err(k[i],n[i])
    return xx,yy,eff.reshape(shape),err.reshape(shape),n.reshape(shape)

def prof2d(x,y,z,bins=20) :
    h0,xx,yy=np.histogram2d(x,y,bins=bins)
    hz,xx,yy=np.histogram2d(x,y,bins=(xx,yy),weights=z)
    hz2,xx,yy=np.histogram2d(x,y,bins=(xx,yy),weights=z**2)
    n=(h0+(h0==0)).astype(float)
    z=hz/n
    z2=hz2/n
    var=z2-z**2
    err=np.sqrt(var*(var>0))
    x=xx[:-1]+(xx[1]-xx[0])/2.
    y=yy[:-1]+(yy[1]-yy[0])/2.
    return x,y,z,err



In [3]:

    
## open input file
hdulist = pyfits.open(simulation_catalog_filename)
table = hdulist["ZCATALOG"].data
print(table.dtype.names)









    



('TARGETID', 'CHI2', 'COEFF', 'Z', 'ZERR', 'ZWARN', 'NPIXELS', 'SPECTYPE', 'SUBTYPE', 'NCOEFF', 'DELTACHI2', 'BRICKNAME', 'NUMEXP', 'NUMTILE', 'RA', 'DEC', 'MAG', 'BRICKID', 'BRICK_OBJID', 'FLUX_G', 'FLUX_R', 'FLUX_Z', 'FLUX_W1', 'FLUX_W2', 'MW_TRANSMISSION_G', 'MW_TRANSMISSION_R', 'MW_TRANSMISSION_Z', 'MW_TRANSMISSION_W1', 'MW_TRANSMISSION_W2', 'PSFDEPTH_G', 'PSFDEPTH_R', 'PSFDEPTH_Z', 'GALDEPTH_G', 'GALDEPTH_R', 'GALDEPTH_Z', 'PSFDEPTH_W1', 'PSFDEPTH_W2', 'SHAPEDEV_R', 'SHAPEDEV_E1', 'SHAPEDEV_E2', 'SHAPEEXP_R', 'SHAPEEXP_E1', 'SHAPEEXP_E2', 'SUBPRIORITY', 'DESI_TARGET', 'BGS_TARGET', 'MWS_TARGET', 'HPXPIXEL', 'MOCKID', 'CONTAM_TARGET', 'TRUEZ', 'TRUESPECTYPE', 'TEMPLATETYPE', 'TEMPLATESUBTYPE', 'TEMPLATEID', 'SEED', 'VDISP', 'OIIFLUX', 'HBETAFLUX', 'TEFF', 'LOGG', 'FEH')



In [4]:

    
# quickcat parameters
quickcat_params=dict()
# quickcat output table (for display purpose only)
qtable=None



In [5]:

    
if True :
    # run the quickcat simulation in this cell (don't necessarily have to do this to 
    # follow the rest of the notebook)
    # use default parameters or the ones in the file specified above 
    # (and probably obtained with a previous run of this script) if exist
    
    input_quickcat_param_filename = None
    if os.path.isfile(quickcat_param_filename) :
        input_quickcat_param_filename = quickcat_param_filename
    
    # dummy tiles
    targets_in_tile=dict()
    targets_in_tile[0]=table["TARGETID"]

    # dummy obs. conditions
    tmp = reference_conditions['DARK']
    tmp['TILEID']=0
    obsconditions=dict()
    for k in tmp :
        obsconditions[k]=np.array([tmp[k],])
    
    #qtable = table.copy()
    hdulist = pyfits.open(simulation_catalog_filename)
    qtable = hdulist["ZCATALOG"].data

    # run quickcat
    # ignore_obsconditions because it only add extra noise
    z, zerr, zwarn = get_observed_redshifts(qtable,qtable,targets_in_tile,obsconditions,
                                                parameter_filename=quickcat_param_filename,
                                                ignore_obscondition=True)
    # replace z,zwarn and write quickcat
    qtable["Z"]=z
    qtable["ZWARN"]=zwarn
    hdulist["ZCATALOG"].data = qtable
    hdulist.writeto(quickcat_catalog_filename,overwrite=True)
    print("done")









    



INFO:quickcat.py:444:get_observed_redshifts: SKY use constant sigmav = 9999 km/s
WARNING:quickcat.py:248:get_redshift_efficiency: using default redshift efficiency of 0.98 for SKY
INFO:quickcat.py:491:get_observed_redshifts: SKY n_failed/n_tot=0/7187=0.000
INFO:quickcat.py:412:get_observed_redshifts: ELG sigma=0.00064 index=0.94 median zerr=0.00016
INFO:quickcat.py:491:get_observed_redshifts: ELG n_failed/n_tot=575/29469=0.020
INFO:quickcat.py:432:get_observed_redshifts: LOWZ QSO sigma=0.00111 index=0.54 median zerr=0.00131
INFO:quickcat.py:439:get_observed_redshifts: LYA QSO sigma=0.00046 index=0.17 median zerr=0.00140
INFO:quickcat.py:230:get_redshift_efficiency: LOWZ QSO eff = sigmoid with cutoff = 23.158 fudge = 0.311
INFO:quickcat.py:238:get_redshift_efficiency: LYA QSO eff = sigmoid with cutoff = 23.206 fudge = 0.086
INFO:quickcat.py:491:get_observed_redshifts: QSO n_failed/n_tot=192/14536=0.013
INFO:quickcat.py:444:get_observed_redshifts: BGS use constant sigmav = 37.7 km/s
INFO:quickcat.py:491:get_observed_redshifts: BGS n_failed/n_tot=37/12042=0.003
INFO:quickcat.py:422:get_observed_redshifts: LRG sigma=0.00036 index=0.60 median zerr=0.00037
INFO:quickcat.py:217:get_redshift_efficiency: LRG eff = sigmoid with cutoff = 26.819 fudge = 1.047
INFO:quickcat.py:491:get_observed_redshifts: LRG n_failed/n_tot=35/21177=0.002
INFO:quickcat.py:444:get_observed_redshifts: STAR use constant sigmav = 51.51 km/s
WARNING:quickcat.py:248:get_redshift_efficiency: using default redshift efficiency of 0.98 for STAR
INFO:quickcat.py:491:get_observed_redshifts: STAR n_failed/n_tot=262/5399=0.049
done



In [6]:

    
# open quickcat catalog
if qtable is None and os.path.isfile(quickcat_catalog_filename) :
    qcat_hdulist = pyfits.open(quickcat_catalog_filename)
    qtable = qcat_hdulist["ZCATALOG"].data

ELG redshift efficiency

We assume the ELG redshift efficiency is a function of

the S/N in the emission lines, approximately proportional to OII flux.
the S/N in the continuum, approximately proportional to the r-band flux.
the redshift

We know that for a given ELG, the S/N in the lines varies with redshift according to the flux limit defined in the FDR. So, we will scale the OII flux with this flux limit to account for some of the redshift dependency. We ignore the evolution of the continuum S/N with redshift for fixed r-band magnitude.

We model the efficiency with an error function,

$ Eff(SNR) = \frac{1}{2} \left( 1+Erf \left( \frac{SNR-3}{b \sqrt{2}} \right) \right) $

with

$SNR = \sqrt{ \left( 7 \frac{OII flux}{fluxlimit} \right)^2 + \left( a \times rflux \right)^2 }$

$a$ is the continuum $SNR$ normalization, which is proportionnal to the r-band flux.

$b$ is a fudge factor. One would have $b = 1$ if $SNR$ was the variable that determines the redshift efficiency.

However $SNR$ is only a proxy that is not 100% correlated with the efficiency, so we expect $b>1$.



In [7]:

    
# OII flux limit (FDR), the has-build version should be recomputed but is probably not very different
filename = resource_filename('desisim', 'data/elg_oii_flux_threshold_fdr.txt')
fdr_z, fdr_flux_limit = np.loadtxt(filename, unpack=True)

plt.figure()
plt.plot(fdr_z, fdr_flux_limit)
plt.ylim([0,1.5e-16])
plt.xlabel("Redshift")
plt.ylabel("OII flux limit (ergs/s/cm2)")
plt.grid()

Measured ELG efficiency as a function of rmag and oii flux



In [8]:

    
######################
elgs=(table["TEMPLATETYPE"]=="ELG")&(table["TRUEZ"]>0.6)&(table["TRUEZ"]<1.6)
z=table["Z"][elgs]
tz=table["TRUEZ"][elgs]
dz=z-tz
good=(table["ZWARN"][elgs]==0)
rflux=table["FLUX_R"][elgs]
print("Number of ELGs={}".format(rflux.size))
rflux=rflux*(rflux>0)+0.00001*(rflux<=0)
oiiflux=table["OIIFLUX"][elgs]
oiiflux=oiiflux*(oiiflux>0)+1e-20*(oiiflux<=0)

qgood=None
if qtable is not None : #quickcat output
    qgood=(qtable["ZWARN"][elgs]==0)

######################

#good=oiiflux>8e-17 #debug to verify indexation

bins2d=20
rmag=-2.5*np.log10(rflux)+22.5
xx,yy,eff2d,err2d,nn2d = efficiency2d(np.log10(oiiflux),rmag,good,bins=bins2d)


plt.figure()
plt.imshow(eff2d.T,origin=0,extent=(xx[0],xx[-1],yy[0],yy[-1]),vmin=0.2,aspect="auto")
plt.xlabel("log10([OII] flux)")
plt.ylabel("rmag")
plt.colorbar()









    



Number of ELGs=29110






    Out[8]:





<matplotlib.colorbar.Colorbar at 0x7f4bd7105ef0>

Model



In [9]:

    
# model ELG efficiency vs rflux and oiiflux
oiiflux = table["OIIFLUX"][elgs]
oiiflux = oiiflux*(oiiflux>=0)+0.00001*(oiiflux<=0)
fluxlimit=np.interp(z,fdr_z,fdr_flux_limit)
fluxlimit[fluxlimit<=0]=1e-20
snr_lines=7*oiiflux/fluxlimit

def elg_efficiency_model_2d(params,log10_snr_lines,rmag) :
    p = params
    snr_tot  = np.sqrt( (p[0]*10**log10_snr_lines)**2 + (p[1]*10**(-0.4*(rmag-22.5))) **2 )
    return 0.5*(1.+np.erf((snr_tot-3)/(np.sqrt(2.)*p[2])))
     
    
def elg_efficiency_2d_residuals(params,log10_snr_lines,mean_rmag,eff2d,err2d) :
    
    model = elg_efficiency_model_2d(params,log10_snr_lines,mean_rmag)
    #res = (eff2d-model)
    res = (eff2d-model)/err2d #np.sqrt(err2d**2+(0.1*(eff2d>0.9))**2)
    res = res[(err2d<2)&(mean_rmag>22)]
    #chi2 = np.sum(res**2)
    #print("params={} chi2/ndata={}/{}={}".format(params,chi2,res.size,chi2/res.size))
    return res



# 2d fit
#good=snr_lines>4. # debug
#good=rmag<22. # debug
xx,yy,eff2d_bis,err2d_bis,nn = efficiency2d(np.log10(snr_lines),rmag,good)
x1d = xx[:-1]+(xx[1]-xx[0])
y1d = yy[:-1]+(yy[1]-yy[0]) 
x2d=np.tile(x1d,(y1d.size,1)).T
y2d=np.tile(y1d,(x1d.size,1))

#elg_efficiency_params=[1,3,2]
elg_efficiency_params=[1,2.,1,0,0]#,1,2,1,0,0]
if 0 :
    meff2d=elg_efficiency_model_2d(elg_efficiency_params,x2d,y2d)
    i=(y2d.ravel()>22.)&(y2d.ravel()<22.4)&(err2d.ravel()<1)
    plt.plot(x2d.ravel()[i],eff2d.ravel()[i],"o")
    plt.plot(x2d.ravel()[i],meff2d.ravel()[i],"o")
    
    
result=scipy.optimize.least_squares(elg_efficiency_2d_residuals,elg_efficiency_params,args=(x2d,y2d,eff2d_bis,err2d_bis))
elg_efficiency_params=result.x

quickcat_params["ELG"]=dict()
quickcat_params["ELG"]["EFFICIENCY"]=dict()
quickcat_params["ELG"]["EFFICIENCY"]["SNR_LINES_SCALE"]=float(elg_efficiency_params[0])
quickcat_params["ELG"]["EFFICIENCY"]["SNR_CONTINUUM_SCALE"]=float(elg_efficiency_params[1])
quickcat_params["ELG"]["EFFICIENCY"]["SIGMA_FUDGE"]=float(elg_efficiency_params[2])

print("Best fit parameters for ELG efficiency model:") 
print(elg_efficiency_params)
print("SNR_lines = {:4.3f} * 7 * OIIFLUX/limit".format(elg_efficiency_params[0]))
print("SNR_cont  = {:4.3f} * R_FLUX".format(elg_efficiency_params[1]))
print("sigma fudge = {:4.3f}".format(elg_efficiency_params[2]))









    



Best fit parameters for ELG efficiency model:
[0.52175753 6.03820505 2.11466464 0.         0.        ]
SNR_lines = 0.522 * 7 * OIIFLUX/limit
SNR_cont  = 6.038 * R_FLUX
sigma fudge = 2.115



In [10]:

    
#params[0]=0.001 # no dependence on rmag
meff=elg_efficiency_model_2d(elg_efficiency_params,np.log10(snr_lines),rmag)
xx,yy,meff2d,merr=prof2d(np.log10(oiiflux),rmag,meff,bins=bins2d)
#plt.imshow(meff2d.T,aspect="auto")

plt.imshow(meff2d.T,origin=0,extent=(xx[0],xx[-1],yy[0],yy[-1]),aspect="auto")
plt.colorbar()


    

if 1 :
    plt.figure()
    print("meff2d.shape=",meff2d.shape)
    ii=np.arange(meff2d.shape[0])
    y1=eff2d[ii,-ii]
    e1=err2d[ii,-ii]
    y2=meff2d[ii,-ii]
    ok=(e1<1)
    plt.errorbar(ii[ok],y1[ok],e1[ok],fmt="o",label="input")
    plt.plot(ii[ok],y2[ok],"-",label="model")
    plt.legend(loc="lower right")
    plt.xlabel("linear combination of log10([OII] flux) and rmag")
    plt.ylabel("efficiency")

plt.figure()
bins1d=20
x,eff1d,err1d   = efficiency(rmag,good,bins=bins1d)
x,meff1d,merr,nn = prof(rmag,meff,bins=bins1d)
plt.errorbar(x,eff1d,err1d,fmt="o",label="input")
plt.plot(x,meff1d,"-",label="model")

if qgood is not None : #quickcat output
    x,eff1d,err1d   = efficiency(rmag,qgood,bins=bins1d)
    plt.errorbar(x,eff1d,err1d,fmt="x",label="qcat run")

plt.legend(loc="lower left")
plt.xlabel("rmag")
plt.ylabel("efficiency")


plt.figure()
bins1d=20
x,eff1d,err1d   = efficiency(np.log10(oiiflux),good,bins=bins1d)
x,meff1d,merr,nn = prof(np.log10(oiiflux),meff,bins=bins1d)
plt.errorbar(x,eff1d,err1d,fmt="o",label="input")
plt.plot(x,meff1d,"-",label="model")

if qgood is not None : #quickcat output
    x,eff1d,err1d   = efficiency(np.log10(oiiflux),qgood,bins=bins1d)
    plt.errorbar(x,eff1d,err1d,fmt="x",label="qcat run")

plt.legend(loc="lower right")
plt.xlabel("log10(oiiflux)")
plt.ylabel("efficiency")


plt.figure()
fcut=8e-17
mcut=22.5
s=(oiiflux<fcut)&(rmag>mcut) # select faint ones to increase contrast in z
bins=100
x,eff1d,err1d   = efficiency(tz[s],good[s],bins=bins)
x,meff1d,merr,nn = prof(tz[s],meff[s],bins=bins)
plt.errorbar(x,eff1d,err1d,fmt="o",label="input")
plt.plot(x,meff1d,"-",label="model")

if qgood is not None : #quickcat output
    x,eff1d,err1d   = efficiency(tz[s],qgood[s],bins=bins1d)
    plt.errorbar(x,eff1d,err1d,fmt="x",label="qcat run")



plt.legend(loc="upper left",title="Faint ELGs with [OII] flux<{} and rmag>{}".format(fcut,mcut))
plt.xlabel("redshift")
plt.ylabel("efficiency")
plt.ylim([0.,1.4])









    



meff2d.shape= (20, 20)






    Out[10]:





(0.0, 1.4)

ELG redshift uncertainty

Power law of [OII] flux (proxy for all lines)



In [11]:

    
#ELG redshift uncertainty
######################
elgs=(table["TEMPLATETYPE"]=="ELG")&(table["TRUEZ"]>0.6)&(table["TRUEZ"]<1.6)
z=table["Z"][elgs]
dz=z-table["TRUEZ"][elgs]
good=(table["ZWARN"][elgs]==0)&(np.abs(dz/(1+z))<0.003)
rflux=table["FLUX_R"][elgs]
print("Number of ELGs={}".format(rflux.size))
rflux=rflux*(rflux>0)+0.00001*(rflux<=0)
oiiflux=table["OIIFLUX"][elgs]
oiiflux=oiiflux*(oiiflux>0)+1e-20*(oiiflux<=0)
lflux=np.log10(oiiflux)

qz=None
qdz=None
if qtable is not None : # quickcat output
    qz=qtable["Z"][elgs]
    qdz=qz-qtable["TRUEZ"][elgs]
    qgood=(qtable["ZWARN"][elgs]==0)&(np.abs(qdz/(1+qz))<0.003)

######################

bins=20

binlflux,var,err,nn=prof(lflux[good],((dz/(1+z))**2)[good],bins=bins)
binflux=10**(binlflux)
var_err = np.sqrt(2/nn)*var
rms=np.sqrt(var)
rmserr=0.5*var_err/rms

def redshift_error(params,flux) :
    return params[0]/(1e-9+flux)**params[1]

def redshift_error_residuals(params,flux,rms,rmserror) :
    model = redshift_error(params,flux)
    res = (rms-model)/np.sqrt(rmserror**2+1e-6**2)
    return res
 
#plt.plot(binlflux,rms,"o",label="meas")
plt.errorbar(binlflux,rms,rmserr,fmt="o",label="sim")
params=[0.0006,1.]
binoiiflux=np.array(10**binlflux)
result=scipy.optimize.least_squares(redshift_error_residuals,params,args=(binoiiflux*1e17,rms,rmserr))
params=result.x
elg_uncertainty_params = params
print("ELG redshift uncertainty parameters = ",params)
quickcat_params["ELG"]["UNCERTAINTY"]=dict()
quickcat_params["ELG"]["UNCERTAINTY"]["SIGMA_17"]=float(elg_uncertainty_params[0])
quickcat_params["ELG"]["UNCERTAINTY"]["POWER_LAW_INDEX"]=float(elg_uncertainty_params[1])

m=redshift_error(params,10**binlflux*1e17)
plt.plot(binlflux,m,"-",label="model")

if qz is not None :
    qbinlflux,qvar,qerr,nn=prof(lflux[qgood],((qdz/(1+qz))**2)[qgood],bins=bins)
    qbinflux=10**(qbinlflux)
    qvar_err = np.sqrt(2/nn)*qvar
    qrms=np.sqrt(qvar)
    qrmserr=0.5*qvar_err/qrms
    plt.errorbar(qbinlflux,qrms,qrmserr,fmt="x",label="quickcat")

plt.legend(loc="upper right",title="ELG")
plt.xlabel("log10([oII] flux)")
plt.ylabel("rms dz/(1+z)")









    



Number of ELGs=29110
ELG redshift uncertainty parameters =  [6.43012258e-04 9.35233498e-01]






    Out[11]:





Text(0,0.5,'rms dz/(1+z)')

ELG catastrophic failure rate

Fraction of targets with ZWARN=0 and $|\Delta z/(1+z)|>0.003$



In [12]:

    
nbad = np.sum((table["ZWARN"][elgs]==0)&(np.abs(dz/(1+z))>0.003))
ntot = np.sum(table["ZWARN"][elgs]==0)
frac = float(nbad/float(ntot))
print("ELG catastrophic failure rate={}/{}={:4.3f}".format(nbad,ntot,frac))
quickcat_params["ELG"]["FAILURE_RATE"]=frac

qnbad = np.sum((qtable["ZWARN"][elgs]==0)&(np.abs(qdz/(1+qz))>0.003))
qntot = np.sum(qtable["ZWARN"][elgs]==0)
qfrac = float(qnbad/float(qntot))
print("quickcat run ELG catastrophic failure rate={}/{}={:4.3f}".format(qnbad,qntot,qfrac))









    



ELG catastrophic failure rate=564/26720=0.021
quickcat run ELG catastrophic failure rate=563/26430=0.021

LRG redshift efficiency

Sigmoid function of the r-band magnitude

$Eff = \frac{1}{1+exp (( rmag - a ) / b))}$



In [13]:

    
# simply use RFLUX for snr
######################
lrgs=(table["TEMPLATETYPE"]=="LRG")
z=table["Z"][lrgs]
tz=table["TRUEZ"][lrgs]
dz=z-tz
good=(table["ZWARN"][lrgs]==0)
rflux=table["FLUX_R"][lrgs]
print("Number of LRGs={}".format(rflux.size))
rflux=rflux*(rflux>0)+0.00001*(rflux<=0)
rmag=-2.5*np.log10(rflux)+22.5

qgood=None
if qtable is not None : #quickcat output
    qgood=(qtable["ZWARN"][lrgs]==0)


######################

bins=15
bin_rmag,eff,err=efficiency(rmag,good,bins=bins)
print("eff=",eff)
print("err=",err)
plt.errorbar(bin_rmag,eff,err,fmt="o",label="sim")


def sigmoid(params,x) :
    return 1/(1+np.exp((x-params[0])/params[1]))

def sigmoid_residuals(params,x,y,err) :
    m = sigmoid(params,x)
    res = (m-y)/err
    return res

lrg_efficiency_params=[26.,1.]
result=scipy.optimize.least_squares(sigmoid_residuals,lrg_efficiency_params,args=(bin_rmag,eff,err)) 
lrg_efficiency_params=result.x
plt.plot(bin_rmag,sigmoid(lrg_efficiency_params,bin_rmag),"-",label="model")

if qgood is not None:
    bin_rmag,eff,err=efficiency(rmag,qgood,bins=bins)
    plt.errorbar(bin_rmag,eff,err,fmt="x",label="quickcat run")

plt.xlabel("rmag")
plt.ylabel("efficiency")
plt.legend(loc="lower left")

print("LRG redshift efficiency parameters = ",lrg_efficiency_params)

quickcat_params["LRG"]=dict()
quickcat_params["LRG"]["EFFICIENCY"]=dict()
quickcat_params["LRG"]["EFFICIENCY"]["SIGMOID_CUTOFF"]=float(lrg_efficiency_params[0])
quickcat_params["LRG"]["EFFICIENCY"]["SIGMOID_FUDGE"]=float(lrg_efficiency_params[1])









    



Number of LRGs=21177
eff= [1.         1.         1.         0.9976247  1.         0.99685535
 0.99852071 0.99737945 0.9954181  0.99274942 0.99169694 0.98456416
 0.98784722 0.99161426 0.97619048]
err= [0.00752226 0.00344551 0.00283688 0.00141844 0.00283688 0.00141844
 0.00141844 0.00141844 0.00131965 0.00144485 0.00146168 0.0021447
 0.00263579 0.00443925 0.0241891 ]
LRG redshift efficiency parameters =  [26.81863486  1.04715984]



In [14]:

    
meff=sigmoid(lrg_efficiency_params,rmag)
plt.figure()
mcut=22.
s=(rmag>mcut) # select faint ones to increase contrast in z
bins=50
x,eff1d,err1d   = efficiency(tz[s],good[s],bins=bins)
x,meff1d,merr,nn = prof(tz[s],meff[s],bins=bins)
plt.errorbar(x,eff1d,err1d,fmt="o",label="sim")
plt.plot(x,meff1d,"-",label="model")
plt.legend(loc="upper left",title="Faint LRGs with rmag>{}".format(mcut))
plt.xlabel("redshift")
plt.ylabel("efficiency")









    Out[14]:





Text(0,0.5,'efficiency')

LRG redshift uncertainty

Power law of broad band flux



In [15]:

    
# LRGs redshift uncertainties

######################
lrgs=(table["TEMPLATETYPE"]=="LRG")
z=table["Z"][lrgs]
tz=table["TRUEZ"][lrgs]
dz=z-tz
good=(table["ZWARN"][lrgs]==0)&(np.abs(dz/(1+z))<0.003)
rflux=table["FLUX_R"][lrgs]
print("Number of LRGs={}".format(rflux.size))
rflux=rflux*(rflux>0)+0.00001*(rflux<=0)
rmag=-2.5*np.log10(rflux)+22.5

qz=None
qdz=None
if qtable is not None : # quickcat output
    qz=qtable["Z"][lrgs]
    qdz=qz-qtable["TRUEZ"][lrgs]
    qgood=(qtable["ZWARN"][lrgs]==0)&(np.abs(qdz/(1+qz))<0.003)


######################

bins=20
binmag,var,err,nn=prof(rmag[good],((dz/(1+z))**2)[good],bins=bins)
binflux=10**(-0.4*(binmag-22.5))
var_err = np.sqrt(2/nn)*var
rms=np.sqrt(var)
rmserr=0.5*var_err/rms

params=[1.,1.2]
result=scipy.optimize.least_squares(redshift_error_residuals,params,args=(binflux,rms,rmserr))
params=result.x
print("LRG redshift error parameters = ",params)
quickcat_params["LRG"]["UNCERTAINTY"]=dict()
quickcat_params["LRG"]["UNCERTAINTY"]["SIGMA_17"]=float(params[0])
quickcat_params["LRG"]["UNCERTAINTY"]["POWER_LAW_INDEX"]=float(params[1])

model = redshift_error(params,binflux)

plt.errorbar(binmag,rms,rmserr,fmt="o",label="sim")
plt.plot(binmag,model,"-",label="model")


if qz is not None :
    qbinmag,qvar,qerr,nn=prof(rmag[qgood],((qdz/(1+qz))**2)[qgood],bins=bins)
    qvar_err = np.sqrt(2/nn)*qvar
    qrms=np.sqrt(qvar)
    qrmserr=0.5*qvar_err/qrms
    plt.errorbar(qbinmag,qrms,qrmserr,fmt="x",label="quickcat")




plt.legend(loc="upper left",title="LRG")
plt.xlabel("rmag")
plt.ylabel("rms dz/(1+z)")









    



Number of LRGs=21177
LRG redshift error parameters =  [3.57946732e-04 6.02482440e-01]






    Out[15]:





Text(0,0.5,'rms dz/(1+z)')

LRG catastrophic failure rate

Fraction of targets with ZWARN=0 and $|\Delta z/(1+z)|>0.003$



In [16]:

    
nbad = np.sum((table["ZWARN"][lrgs]==0)&(np.abs(dz/(1+z))>0.003))
ntot = np.sum(table["ZWARN"][lrgs]==0)
frac = float(nbad/float(ntot))
print("LRG catastrophic failure rate={}/{}={:4.3f}".format(nbad,ntot,frac))
quickcat_params["LRG"]["FAILURE_RATE"]=frac

qnbad = np.sum((qtable["ZWARN"][lrgs]==0)&(np.abs(qdz/(1+qz))>0.003))
qntot = np.sum(qtable["ZWARN"][lrgs]==0)
qfrac = float(qnbad/float(qntot))
print("quickcat run LRG catastrophic failure rate={}/{}={:4.3f}".format(qnbad,qntot,qfrac))









    



LRG catastrophic failure rate=39/21020=0.002
quickcat run LRG catastrophic failure rate=35/21033=0.002



In [17]:

    
# choice of redshift for splitting between "lower z / tracer" QSOs and Lya QSOs
zsplit = 2.0

QSO tracers (z<~2) redshift efficiency

Sigmoid function of the r-band magnitude

$Eff = \frac{1}{1+exp (( rmag - a ) / b))}$



In [18]:

    
# simply use RFLUX for snr
######################
qsos=(table["TEMPLATETYPE"]=="QSO")&(table["TRUEZ"]<zsplit)
z=table["Z"][qsos]
tz=table["TRUEZ"][qsos]
dz=z-tz
good=(table["ZWARN"][qsos]==0)
rflux=table["FLUX_R"][qsos]
print("Number of QSOs={}".format(rflux.size))
rflux=rflux*(rflux>0)+0.00001*(rflux<=0)
rmag=-2.5*np.log10(rflux)+22.5

qgood=None
if qtable is not None : # quickcat output
    qgood=(qtable["ZWARN"][qsos]==0)


######################

bins=30
bin_rmag,eff,err=efficiency(rmag,good,bins=bins)
plt.errorbar(bin_rmag,eff,err,fmt="o",label="sim")
qso_efficiency_params=[23.,0.3]
result=scipy.optimize.least_squares(sigmoid_residuals,qso_efficiency_params,args=(bin_rmag,eff,err)) 
qso_efficiency_params=result.x
plt.plot(bin_rmag,sigmoid(qso_efficiency_params,bin_rmag),"-",label="model")

if qgood is not None :
    bin_rmag,eff,err=efficiency(rmag,qgood,bins=bins)
    plt.errorbar(bin_rmag,eff,err,fmt="x",label="quickcat run")

plt.xlabel("rmag")
plt.ylabel("efficiency")
plt.legend(loc="lower left")

print("QSO redshift efficiency parameters = ",qso_efficiency_params)
quickcat_params["QSO_ZSPLIT"]=zsplit
quickcat_params["LOWZ_QSO"]=dict()
quickcat_params["LOWZ_QSO"]["EFFICIENCY"]=dict()
quickcat_params["LOWZ_QSO"]["EFFICIENCY"]["SIGMOID_CUTOFF"]=float(qso_efficiency_params[0])
quickcat_params["LOWZ_QSO"]["EFFICIENCY"]["SIGMOID_FUDGE"]=float(qso_efficiency_params[1])









    



Number of QSOs=10249
QSO redshift efficiency parameters =  [23.15765723  0.3110151 ]



In [19]:

    
meff=sigmoid(qso_efficiency_params,rmag)
plt.figure()
mcut=22.
s=(rmag>mcut) # select faint ones to increase contrast in z
bins=50
x,eff1d,err1d   = efficiency(tz[s],good[s],bins=bins)
x,meff1d,merr,nn = prof(tz[s],meff[s],bins=bins)
plt.errorbar(x,eff1d,err1d,fmt="o",label="input")
plt.plot(x,meff1d,"-",label="model")

if qgood is not None :
    x,qeff1d,qerr1d   = efficiency(tz[s],qgood[s],bins=bins)
    plt.errorbar(x,qeff1d,qerr1d,fmt="x",label="quickcat")
    
plt.legend(loc="upper left",title="Faint tracer QSOs with rmag>{}".format(mcut))
plt.xlabel("redshift")
plt.ylabel("efficiency")
plt.ylim([0.5,1.2])









    Out[19]:





(0.5, 1.2)

QSO (z<2) redshift uncertainty

Power law of broad band flux



In [20]:

    
# QSO redshift uncertainties
qsos=(table["TEMPLATETYPE"]=="QSO")&(table["TRUEZ"]<zsplit)
z=table["Z"][qsos]
dz=z-table["TRUEZ"][qsos]
good=(table["ZWARN"][qsos]==0)&(np.abs(dz/(1+z))<0.01)
rflux=table["FLUX_R"][qsos]
print("Number of QSOs={}".format(rflux.size))
rflux=rflux*(rflux>0)+0.00001*(rflux<=0)
rmag=-2.5*np.log10(rflux)+22.5

qgood=None
qz=None
qdz=None
if qtable is not None : # quickcat output
    qz=qtable["Z"][qsos]
    qdz=qz-qtable["TRUEZ"][qsos]
    qgood=(qtable["ZWARN"][qsos]==0)&(np.abs(qdz/(1+qz))<0.01)



bins=20
binmag,var,err,nn=prof(rmag[good],((dz/(1+z))**2)[good],bins=bins)
binflux=10**(-0.4*(binmag-22.5))
var_err = np.sqrt(2/nn)*var
rms=np.sqrt(var)
rmserr=0.5*var_err/rms

params=[1.,1.2]
result=scipy.optimize.least_squares(redshift_error_residuals,params,args=(binflux,rms,rmserr))
params=result.x
print("QSO redshift error parameters = ",params)
quickcat_params["LOWZ_QSO"]["UNCERTAINTY"]=dict()
quickcat_params["LOWZ_QSO"]["UNCERTAINTY"]["SIGMA_17"]=float(params[0])
quickcat_params["LOWZ_QSO"]["UNCERTAINTY"]["POWER_LAW_INDEX"]=float(params[1])

model = redshift_error(params,binflux)

plt.errorbar(binmag,rms,rmserr,fmt="o",label="sim")
plt.plot(binmag,model,"-",label="model")

if qz is not None :
    qbinmag,qvar,qerr,nn=prof(rmag[qgood],((qdz/(1+qz))**2)[qgood],bins=bins)
    qvar_err = np.sqrt(2/nn)*qvar
    qrms=np.sqrt(qvar)
    qrmserr=0.5*qvar_err/qrms
    plt.errorbar(qbinmag,qrms,qrmserr,fmt="x",label="quickcat")

plt.legend(loc="upper left",title="Tracer QSO")
plt.xlabel("rmag")
plt.ylabel("rms dz/(1+z)")









    



Number of QSOs=10249
QSO redshift error parameters =  [0.00111248 0.53624373]






    Out[20]:





Text(0,0.5,'rms dz/(1+z)')

Tracer QSO (z<~2) catastrophic failure rate

Fraction of targets with ZWARN=0 and $|\Delta z/(1+z)|>0.003$



In [21]:

    
nbad = np.sum((table["ZWARN"][qsos]==0)&(np.abs(dz/(1+z))>0.003))
ntot = np.sum(table["ZWARN"][qsos]==0)
frac = float(nbad/float(ntot))
print("Tracer QSO catastrophic failure rate={}/{}={:4.3f}".format(nbad,ntot,frac))
quickcat_params["LOWZ_QSO"]["FAILURE_RATE"]=frac

qnbad = np.sum((qtable["ZWARN"][qsos]==0)&(np.abs(qdz/(1+qz))>0.003))
qntot = np.sum(qtable["ZWARN"][qsos]==0)
qfrac = float(qnbad/float(qntot))
print("quickcat run tracer QSO catastrophic failure rate={}/{}={:4.3f}".format(qnbad,qntot,qfrac))









    



Tracer QSO catastrophic failure rate=198/9954=0.020
quickcat run tracer QSO catastrophic failure rate=204/9959=0.020

Lya QSO (z>~2) redshift efficiency

Sigmoid function of the r-band magnitude

$Eff = \frac{1}{1+exp (( rmag - a ) / b))}$



In [22]:

    
# simply use RFLUX for snr
######################
qsos=(table["TEMPLATETYPE"]=="QSO")&(table["TRUEZ"]>zsplit)
z=table["Z"][qsos]
tz=table["TRUEZ"][qsos]
dz=z-tz
good=(table["ZWARN"][qsos]==0)
rflux=table["FLUX_R"][qsos]
print("Number of QSOs={}".format(rflux.size))
rflux=rflux*(rflux>0)+0.00001*(rflux<=0)
rmag=-2.5*np.log10(rflux)+22.5

qgood=None
if qtable is not None : # quickcat output
    qgood=(qtable["ZWARN"][qsos]==0)


######################

bins=30
bin_rmag,eff,err=efficiency(rmag,good,bins=bins)
plt.errorbar(bin_rmag,eff,err,fmt="o",label="sim")
qso_efficiency_params=[23.,0.3]
result=scipy.optimize.least_squares(sigmoid_residuals,qso_efficiency_params,args=(bin_rmag,eff,err)) 
qso_efficiency_params=result.x
plt.plot(bin_rmag,sigmoid(qso_efficiency_params,bin_rmag),"-",label="model")

if qgood is not None :
    bin_rmag,eff,err=efficiency(rmag,qgood,bins=bins)
    plt.errorbar(bin_rmag,eff,err,fmt="x",label="quickcat run")

plt.xlabel("rmag")
plt.ylabel("efficiency")
plt.legend(loc="lower left")

print("QSO redshift efficiency parameters = ",qso_efficiency_params)

quickcat_params["LYA_QSO"]=dict()
quickcat_params["LYA_QSO"]["EFFICIENCY"]=dict()
quickcat_params["LYA_QSO"]["EFFICIENCY"]["SIGMOID_CUTOFF"]=float(qso_efficiency_params[0])
quickcat_params["LYA_QSO"]["EFFICIENCY"]["SIGMOID_FUDGE"]=float(qso_efficiency_params[1])









    



Number of QSOs=4287
QSO redshift efficiency parameters =  [23.20604852  0.08628842]



In [23]:

    
meff=sigmoid(qso_efficiency_params,rmag)
plt.figure()
mcut=22.5
s=(rmag>mcut) # select faint ones to increase contrast in z
bins=50
x,eff1d,err1d   = efficiency(tz[s],good[s],bins=bins)
x,meff1d,merr,nn = prof(tz[s],meff[s],bins=bins)
plt.errorbar(x,eff1d,err1d,fmt="o",label="sim")
plt.plot(x,meff1d,"-",label="model")
plt.legend(loc="upper left",title="Faint Lya QSOs with rmag>{}".format(mcut))
plt.xlabel("redshift")
plt.ylabel("efficiency")
plt.ylim([0.,1.4])









    Out[23]:





(0.0, 1.4)

Lya QSO (z>2) redshift uncertainty

Power law of broad band flux



In [24]:

    
# QSO redshift uncertainties
qsos=(table["TEMPLATETYPE"]=="QSO")&(table["TRUEZ"]>zsplit)
z=table["Z"][qsos]
dz=z-table["TRUEZ"][qsos]
good=(table["ZWARN"][qsos]==0)&(np.abs(dz/(1+z))<0.01)
rflux=table["FLUX_R"][qsos]
print("Number of QSOs={}".format(rflux.size))
rflux=rflux*(rflux>0)+0.00001*(rflux<=0)
rmag=-2.5*np.log10(rflux)+22.5

qgood=None
qz=None
qdz=None
if qtable is not None : # quickcat output
    qz=qtable["Z"][qsos]
    qdz=qz-qtable["TRUEZ"][qsos]
    qgood=(qtable["ZWARN"][qsos]==0)&(np.abs(qdz/(1+qz))<0.01)



bins=20
binmag,var,err,nn=prof(rmag[good],((dz/(1+z))**2)[good],bins=bins)
binflux=10**(-0.4*(binmag-22.5))
var_err = np.sqrt(2/nn)*var
rms=np.sqrt(var)
rmserr=0.5*var_err/rms

params=[1.,1.2]
result=scipy.optimize.least_squares(redshift_error_residuals,params,args=(binflux,rms,rmserr))
params=result.x
print("LYA_QSO redshift error parameters = ",params)
quickcat_params["LYA_QSO"]["UNCERTAINTY"]=dict()
quickcat_params["LYA_QSO"]["UNCERTAINTY"]["SIGMA_17"]=float(params[0])
quickcat_params["LYA_QSO"]["UNCERTAINTY"]["POWER_LAW_INDEX"]=float(params[1])

model = redshift_error(params,binflux)

plt.errorbar(binmag,rms,rmserr,fmt="o",label="sim")
plt.plot(binmag,model,"-",label="model")

if qz is not None :
    qbinmag,qvar,qerr,nn=prof(rmag[qgood],((qdz/(1+qz))**2)[qgood],bins=bins)
    qvar_err = np.sqrt(2/nn)*qvar
    qrms=np.sqrt(qvar)
    qrmserr=0.5*qvar_err/qrms
    plt.errorbar(qbinmag,qrms,qrmserr,fmt="x",label="quickcat")

plt.legend(loc="upper left",title="Lya QSO")
plt.xlabel("rmag")
plt.ylabel("rms dz/(1+z)")









    



Number of QSOs=4287
LYA_QSO redshift error parameters =  [0.00045565 0.17072328]






    Out[24]:





Text(0,0.5,'rms dz/(1+z)')

Lya QSO (z>~2) catastrophic failure rate

Fraction of targets with ZWARN=0 and $|\Delta z/(1+z)|>0.003$



In [25]:

    
nbad = np.sum((table["ZWARN"][qsos]==0)&(np.abs(dz/(1+z))>0.003))
ntot = np.sum(table["ZWARN"][qsos]==0)
frac = float(nbad/float(ntot))
print("Lya QSO catastrophic failure rate={}".format(frac))
quickcat_params["LYA_QSO"]["FAILURE_RATE"]=frac









    



Lya QSO catastrophic failure rate=0.0



In [26]:

    
# write results to a yaml file
with open(quickcat_param_filename, 'w') as outfile:
    yaml.dump(quickcat_params, outfile, default_flow_style=False)



In [ ]: