-
Notifications
You must be signed in to change notification settings - Fork 3
/
systematics_QSO.py
268 lines (244 loc) · 9.06 KB
/
systematics_QSO.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
import fitsio #needed to read data
import numpy as np
from math import *
from optimize import fmin
import os
import pyfits
#try: mkesample_ver = os.environ['MKESAMPLE_VER']
#except: mkesample_ver = None
#mkesample_ver='1.0'
#dirsys = os.environ['MKESAMPLE_DIR']+'/inputFiles/' #change to local directory where maps are
#dir = os.environ['EBOSS_LSS_CATALOGS']+'/'+mkesample_ver+'/' #change to where your catalog files are
#dir = os.environ['MKESAMPLE_SCRATCH']+'/lsscats/'+mkesample_ver+'/' #change to where your catalog files are
dir = '/Users/ashleyross/fitsfiles/'
outdir = '/Users/ashleyross/eboss/'
#outdir=os.environ['EBOSS_LSS_CATALOGS']+'/'+mkesample_ver+'/' #to where files will be written
def ngvsys_ran(sampl,NS,ver,sys,sysmin,sysmax,zmin,zmax,band=-1,wm=''):
#sample is the sample being used, e.g. 'lrg'
#uses files with imaging properties filled
#NS is either 'N' or 'S'
#ver is the version, e.g., 'v1.0_IRt'
#sys is a string containing the name of the systematic to be tested
#sysmin is the minimum value of the systematic to be tested, ~25 is a good value to use for stars
#sysmax is the maximum value of the systematic to be tested, ~200 is a good value to use for stars
#res is Nside for the healpix map, default should be 512 for sky,seeing,air and 256 for extinction and stars
#zmin is the minimum redshift to use, 0.6 is minimum used for lrgs, 0.9 for qsos
#zmax is the maximum redshift to used, 1.0 for lrgs and 2.2 for qsos
stl = []
wstl = []
errl = []
if wm == 'wgdepthext': #this is used to set the weights for extinction
b,m = findlinmb(sampl,NS,ver,'IMAGE_DEPTH_EXT1',.8,2.2,wm='nosys')
if wm == 'wgdepthextext': #this can be used to test that the weights work for everything
b,m = findlinmb(sampl,NS,ver,'IMAGE_DEPTH_EXT1',.8,2.2)
be,me = findlinmb(sampl,NS,ver,'EB_MINUS_V-1',.8,2.2,wm='wgdepthext')
print 'fit parameters for g-band depth'
print b,m
print 'fit parameters for EB_MINUS_V (after weighting for depth)'
print be,me
binng = []
binnr = []
nsysbin = 10 #10 bins are being used
for i in range(0,nsysbin):
binng.append(0)
binnr.append(0)
f = fitsio.read(dir+'eboss_'+ver+'-'+sampl+'-'+NS+'-eboss_'+ver+'.ran.fits') #read random file
print 'reading randoms from ' + dir+'eboss_'+ver+'-'+sampl+'-'+NS+'-eboss_'+ver+'.ran.fits'
nr = 0
sysm = float(nsysbin)/(sysmax-sysmin)
bs = 0
bsr = 0
extc = [4.239,3.303,2.285,1.698,1.263]
for i in range (0,len(f)):
if band == -1:
sysv = f[i][sys]
else:
if sys == 'IMAGE_DEPTH_EXT':
sysv = f[i]['IMAGE_DEPTH'][band]
sysv = luptm(sysv,band)-extc[band]*f[i]['EB_MINUS_V']
else:
sysv = f[i][sys][band]
if sys == 'IMAGE_DEPTH':
sysv = luptm(sysv,band)
bins = int((sysv-sysmin)*sysm)
if bins >= 0 and bins < nsysbin:
binnr[bins] += 1.
else:
bsr += 1.
nr += 1.
print nr
f = fitsio.read(dir+'eboss_'+ver+'-'+sampl+'-'+NS+'-eboss_'+ver+'.dat.fits') #read quasar file
print 'reading quasars from: ' + dir+'eboss_'+ver+'-'+sampl+'-'+NS+'-eboss_'+ver+'.dat.fits'
no = 0
zm = 0
nt = 0
for i in range (0,len(f)):
z = f[i]['Z']
gc = True
um = f[i]['MODELMAG'][0]-f[i]['EXTINCTION'][0]
gm = f[i]['MODELMAG'][1]-f[i]['EXTINCTION'][1]
rm = f[i]['MODELMAG'][2]-f[i]['EXTINCTION'][2]
im = f[i]['MODELMAG'][3]-f[i]['EXTINCTION'][3]
if z > zmin and z < zmax:
no += 1
#w = 1.
#if wm == '':
w = (f[i]['WEIGHT_NOZ']+f[i]['WEIGHT_CP']-1.)*f[i]['WEIGHT_FKP']*f[i]['WEIGHT_SYSTOT'] #weight in current catalog
if wm == 'nosys':
w = (f[i]['WEIGHT_NOZ']+f[i]['WEIGHT_CP']-1.)*f[i]['WEIGHT_FKP'] #compare to result with no systematic weights
if sys == 'SKYFLUX':
sys = 'SKY_FLUX'
if band == -1:
sysv = f[i][sys]
else:
if sys == 'IMAGE_DEPTH_EXT':
sysv = f[i]['IMAGE_DEPTH'][band]
sysv = luptm(sysv,band)-extc[band]*f[i]['EB_MINUS_V']
else:
sysv = f[i][sys][band]
if sys == 'IMAGE_DEPTH':
sysv = luptm(sysv,band)
bins = int((sysv-sysmin)*sysm)
if wm == 'wgdepthext':
sysw = f[i]['IMAGE_DEPTH'][1]
sysw = luptm(sysw,1)-extc[1]*f[i]['EB_MINUS_V']
w = (f[i]['WEIGHT_NOZ']+f[i]['WEIGHT_CP']-1.)*f[i]['WEIGHT_FKP']*(1./(b+m*sysw))
if wm == 'wgdepthextext':
sysw = f[i]['IMAGE_DEPTH'][1]
sysw = luptm(sysw,1)-extc[1]*f[i]['EB_MINUS_V']
w = (f[i]['WEIGHT_NOZ']+f[i]['WEIGHT_CP']-1.)*f[i]['WEIGHT_FKP']*(1./(b+m*sysw))
ext = f[i]['EB_MINUS_V']
w = w*(1./(be+me*ext))
if bins >= 0 and bins < nsysbin:
binng[bins] += 1.*w
else:
bs += w #count numbers outside of sysmin/sysmax
zm += w*z
nt += w
print 'total number, weighted number'
print no,nt
print 'mean redshift'
print zm/nt
print 'total number of randoms/objects '+str(nr)+'/'+str(nt)
print 'number of randoms/objects outside tested range '+str(bsr)+'/'+str(bs)
ave = nt/nr
print 'average number of objects per random is '+ str(ave)
fs = open(dir+'n'+'geboss'+sampl+'_'+NS+ver+'_mz'+str(zmin)+'xz'+str(zmax)+wm+'v'+sys+str(band)+'.dat','w')
xl = []
yl = []
el = []
for i in range(0,nsysbin):
sysv = sysmin + 1./(2.*sysm) + i/sysm
if binnr[i] > 0:
ns = binng[i]/binnr[i]/ave
nse = sqrt(binng[i]/(binnr[i])**2./(ave)**2.+(binng[i]/ave)**2./(binnr[i])**3.) #calculate poisson error
else:
ns = 1. #write out 1.0 1.0 if no pixels at given value of sys
nse = 1.
fs.write(str(sysv)+' '+str(ns)+' '+str(nse)+'\n')
xl.append(sysv)
yl.append(ns)
el.append(nse)
fs.close()
chin = sum((np.array(yl)-1.)**2./np.array(el)**2.)
print 'chi2 for null test'
print chin
xl = np.array(xl)
yl = np.array(yl)
el = np.array(el)
#plotvssys_simp(xl,yl,el,sys)
return xl,yl,el
def findlinmb(sampl,NS,ver,sys,zmin,zmax,wm='',res=''):
#finds linear fit parameters (depth or stellar density relationships hard-coded to expect given resolutions)
if sys == 'star':
res = '256'
if sys == 'depth':
res = '512'
d = np.loadtxt(dir+'n'+'geboss'+sampl+'_'+NS+ver+'_mz'+str(zmin)+'xz'+str(zmax)+wm+str(res)+'v'+sys+'.dat').transpose()
lf = linfit(d[0],d[1],d[2])
inl = np.array([1.,0])
b0,m0 = fmin(lf.chilin,inl)
return b0,m0
class linfit:
def __init__(self,xl,yl,el):
self.xl = xl
self.yl = yl
self.el = el
def chilin(self,bml):
chi = 0
b = bml[0]
m = bml[1]
for i in range(0,len(self.xl)):
y = b+m*self.xl[i]
chi += (self.yl[i]-y)**2./self.el[i]**2.
return chi
def luptm(nmag,bnd):
#calculates SDSS magnitudes from fluxes
b = []
b.append(1.4*10**-10.)
b.append(.9*10**-10.)
b.append(1.2*10**-10.)
b.append(1.8*10**-10.)
b.append(7.4*10**-10.)
return -2.5/log(10.)*(asinh((nmag/10.**9.)/(2.*b[bnd]))+log(b[bnd]))
def calcweights(sample,NS,version,zmin=.8,zmax=2.2,app='.fits'):
#note, zmin/zmax assume LRG sample these need to be change for QSO files
outfile = outdir+'eboss_'+version+'-'+sample+'-'+NS+'-eboss_'+version+'.dat'+app
b,m = findlinmb(sample,NS,version,'IMAGE_DEPTH_EXT1',.8,2.2,wm='nosys')
be,me = findlinmb(sample,NS,version,'EB_MINUS_V-1',.8,2.2,wm='wgdepthext')
print 'fit parameters for g-band depth'
print b,m
print 'fit parameters for EB_MINUS_V (after weighting for depth)'
print be,me
f = fitsio.read(dir+'eboss_'+version+'-'+sample+'-'+NS+'-eboss_'+version+'.dat'+app) #read galaxy/quasar file
#no = 0
#fo = open(os.environ['HOME']+'/mkEsample_work/depthextweights'+sample+'_'+NS+version+'_mz'+str(zmin)+'xz'+str(zmax)+'.dat','w')
# for i in range(0,len(f)):
# sysw = f[i]['IMAGE_DEPTH'][1]
# sysw = luptm(sysw,1)-3.303*f[i]['EB_MINUS_V']
# wd = (1./(b+m*sysw))
# ext = f[i]['EB_MINUS_V']
# we = (1./(be+me*ext))
# wtot = wd*we
# #fo.write(str(ws)+'\n')
# f[i]['WEIGHT_SYSTOT'] = wtot
# #fo.close()
# pyfits.writeto(outfile,f,clobber=True)
# print 'Wrote file with systematic weights to: ' + outfile
#print no
return True
def calcweights_ran(sample,NS,version,zmin=.8,zmax=2.2,app='.fits'):
#note, zmin/zmax assume LRG sample these need to be change for QSO files
#outfile = outdir+'eboss_'+version+'-'+sample+'-'+NS+'-eboss_'+version+'.dat'+app
b,m = findlinmb(sample,NS,version,'IMAGE_DEPTH_EXT1',.8,2.2,wm='nosys')
be,me = findlinmb(sample,NS,version,'EB_MINUS_V-1',.8,2.2,wm='wgdepthext')
print 'fit parameters for g-band depth'
print b,m
print 'fit parameters for EB_MINUS_V (after weighting for depth)'
print be,me
f = fitsio.read(dir+'eboss_'+version+'-'+sample+'-'+NS+'-eboss_'+version+'.ran'+app) #read galaxy/quasar file
#no = 0
#fo = open(os.environ['HOME']+'/mkEsample_work/depthextweights'+sample+'_'+NS+version+'_mz'+str(zmin)+'xz'+str(zmax)+'.dat','w')
fo = open(outdir+'ran'+sample+NS+version+'sysweights.dat','w')
for i in range(0,len(f)):
sysv = f[i]['IMAGE_DEPTH'][1]
sysw = luptm(sysv,1)-3.303*f[i]['EB_MINUS_V']
wd = (1./(b+m*sysw))
ext = f[i]['EB_MINUS_V']
we = (1./(be+me*ext))
wtot = wd*we
fo.write(str(f[i]['RA'])+' '+str(f[i]['DEC'])+' '+str(wtot)+'\n')
# f[i]['WEIGHT_SYSTOT'] = wtot
fo.close()
# pyfits.writeto(outfile,f,clobber=True)
# print 'Wrote file with systematic weights to: ' + outfile
#print no
return True
# if __name__ == '__main__':
# zmin = .6
# zmax = 1.
# sample = 'lrg'
# NS = 'N'
# version = 'v0.8_IRc'
# ngvsys(sample,NS,version,zmin,zmax)
# calcstarweight(sample,NS,version,zmin=.6,zmax=1.)