forked from tdelubac/eBOSSLens
-
Notifications
You must be signed in to change notification settings - Fork 0
/
SDSSObject.py
219 lines (204 loc) · 8.12 KB
/
SDSSObject.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
import os
import numpy as np
import pyfits as pf
from scipy.integrate import quad
from scipy.optimize import minimize
from utils import chi2g, chi2D
class SDSSObject():
'''
SDSSObject
==========
A class for SDSS BOSS or eBOSS objects. Note that each fiberid will also
be separated for better serialization
'''
def __init__(self, plate, mjd, fid, dataVersion, baseDir="/SCRATCH"):
# Data version check
if dataVersion == "v5_7_0" or dataVersion == "v5_7_2":
flag = "BOSS"
elif dataVersion == "v5_10_0":
flag = "eBOSS"
else:
raise ValueError("Undefined data version")
# Construct dirs for loading data
dir1 = os.path.join(baseDir, flag, "data", dataVersion, str(plate))
suffix = str(plate) + "-" + str(mjd) + ".fits"
spFile = os.path.join(dir1, "spPlate-" + suffix)
zbFile = os.path.join(dir1, dataVersion, "spZbest-" + suffix)
zlFile = os.path.join(dir1, dataVersion, "spZline-" + suffix)
# Match fiberid with the index in the plate
hdulist = pf.open(zbFile)
fiberidList = hdulist[1].data.field('FIBERID')
hdulist.close()
findex = np.argwhere(fiberidList == fid)[0][0]
self.fiberid = fid
self.mjd = mjd
self.plate = plate
self.dataVersion = dataVersion
self.baseDir = baseDir
# Load data from spPlate
hdulist = pf.open(spFile)
self.c0 = hdulist[0].header['coeff0']
self.c1 = hdulist[0].header['coeff1']
self.npix = hdulist[0].header['naxis1']
self.wave = 10. ** (self.c0 + self.c1 * np.arange(self.npix))
self.flux = hdulist[0].data[findex]
self.ivar = hdulist[1].data[findex]
# Also load for previous and next fiber
fluxPre = []
fluxNxt = []
if findex != 0:
fluxPre = hdulist[0].data[findex - 1]
self.ivarPre = hdulist[1].data[findex - 1]
if findex != 999:
fluxNxt = hdulist[0].data[findex + 1]
self.ivarNxt = hdulist[1].data[findex + 1]
hdulist.close()
# Load data from spZbest
hdulist = pf.open(zbFile)
self.vdisp = hdulist[1].data.field('VDISP')[findex]
self.synflux = hdulist[2].data[findex]
self.RA = hdulist[1].data.field('PLUG_RA')[findex]
self.DEC = hdulist[1].data.field('PLUG_DEC')[findex]
self.obj_id = hdulist[1].data.field('OBJID')[findex]
self.obj_class = hdulist[1].data.field('CLASS')[findex]
self.obj_type = hdulist[1].data.field('OBJTYPE')[findex]
self.z = hdulist[1].data.field('Z')[findex]
self.zwarning = hdulist[1].data.field('ZWARNING')[findex]
self.z_err = hdulist[1].data.field('Z_ERR')[findex]
self.spectroflux = hdulist[1].data.field('SPECTROFLUX')[findex]
self.rchi2 = hdulist[1].data.field('RCHI2')[findex]
# Also load reduced flux / ivar for previous fiber and next fiber
if findex != 0:
self.reduced_flux_pre = fluxPre - hdulist[2].data[findex - 1]
if findex != 999:
self.reduced_flux_nxt = fluxNxt - hdulist[2].data[findex + 1]
hdulist.close()
# Load data from spZline
hdulist = pf.open(zlFile)
zlineList = hdulist[1].data
hdulist.close()
self.zline = zlineList[np.where(zlineList['fiberid'] == self.fiberid)]
# Additional processing
self.reduced_flux = self.flux - self.synflux
# Additional holder for future QSO continuum removal if needed
self.reduced_flux_QSO = self.reduced_flux
self.nMax = len(self.flux)
self.snSpectra = np.nanmean(self.flux * self.ivar)
# Holder for later of SN array
self.SN = np.zeros(len(self.wave))
def wave2bin(self, waveLength):
'''
SDSSObject.wave2bin(waveLength)
===============================
Return the accordinf bin of given waveLength
Paramters:
waveLength: given wavelength to convert (float)
Returns:
b: the bin of the given wavelength (same type as waveLength)
'''
b = int((np.log10(waveLength) / self.c1) - self.c0 / self.c1)
if b < 0:
b = 0
elif b > self.nMax:
b = self.nMax
return b
def mask(self, lineList):
'''
SDSSObject.mask(lineList)
=========================
Mask part of the spectra by setting ivar = 0
Parameters:
lineList: 2-d list with [:, 0] as start and [:, 1] as end
Returns:
None
'''
for each in lineList:
self.ivar[self.wave2bin(each[0]): self.wave2bin(each[1])] = 0
def nearLine(self, x0, width=5.0):
'''
nearLine(mjd, plate, i, x0, obj, width=10.0)
============================================
Whether given wavelength is near a masked or not
Parameters:
x0: given wavelength
width: default 10.0, a parameter for nearness
Returns:
a boolean of whether near or not
'''
crit = abs(self.zline['linewave'] * (1 + self.z) - x0) < width
if np.any(crit):
return True
else:
return False
def singletFit(self, bounds, initParam, paramLim):
'''
SDSSObject.singleFit(bounds, initParam, paramLim)
====================================================
Fit singlet
Parameters:
bounds: the range that will be fitted
initParam: initial parameter
paramLim: limit of parameter to search for
Returns:
res.x: parameter value after fitting
res.fun: chisquare of the fitting
'''
res = minimize(chi2g, initParam, args=(self.wave[bounds],
self.reduced_flux[bounds],
self.ivar[bounds]),
method='SLSQP', bounds=paramLim)
return res.x, res.fun
def doubletFit(self, bounds, initParam, paramLim, NxtPre=None):
'''
SDSSObject.singleFit(bounds, initParam, paramLim)
====================================================
Fit doublet
Parameters:
bounds: the range that will be fitted
initParam: initial parameter
paramLim: limit of parameter to search for
NxtPre: use next, present or previous fiber data
Returns:
res.x: parameter value after fitting
res.fun: chisquare of the fitting
'''
if NxtPre is None:
rf = self.reduced_flux[bounds]
iv = self.ivar[bounds]
elif NxtPre == "nxt" and self.fiberid != 1000:
rf = self.reduced_flux_nxt[bounds]
iv = self.ivarNxt[bounds]
elif NxtPre == "pre" and self.fiberid != 1:
rf = self.reduced_flux_pre[bounds]
iv = self.ivarPre[bounds]
else:
raise ValueError("Invalid NxtPre value")
res = minimize(chi2D, initParam, args=(self.wave[bounds], rf, iv),
method='SLSQP', bounds=paramLim)
return res.x, res.fun
def radEinstein(self, zs):
'''
SDSSObject.radEinstein(z1, z2)
==============================
Estimated Einstein Radius from Single Isothermal Sphere (SIS) model
Parameters:
zs: source redshift
Returns:
None
'''
# Necessary parameter
H0 = 72e3
c = 299792458.0
OmegaM = 0.258
OmegaL = 0.742
vdisp = self.vdisp * 1000.0
# Function needed for numerical computation of angular diameter distance
def x12(z):
return 1.0 / np.sqrt((1.0 - OmegaM - OmegaL) * (1.0 + z) *
(1.0 + z) + OmegaM * (1.0 + z) ** 3.0 + OmegaL)
# compute ang. diam. distances
Ds = ((c / H0) * quad(x12, 0.0, zs)[0]) / (1.0 + zs)
Dls = ((c / H0) * quad(x12, self.z, zs)[0]) / (1.0 + zs)
# return value in arcsec
coeff = 3600.0 * (180.0 / np.pi)
return coeff * 4.0 * np.pi * vdisp * vdisp * Dls / (c * c * Ds)