Commit 0431b383 authored by Monica Rainer's avatar Monica Rainer
Browse files

Fix bug in SNR computation

SNR was overestimated by a factor sqrt(2). Other minor issues were fixed.
parent 6ef82986
......@@ -241,6 +241,7 @@ CONFIG['HWTM'] = 6. # hwtm limit for contamination between A and B
CONFIG['S1D'] = True # flag for creation of s1d output (True/False)
CONFIG['S1D_NORM'] = False # flag for normalization of s1d output (True/False)
CONFIG['S1D_STEP'] = 0.001 # constant step to use in s1d spectra
CONFIG['NCOSMIC'] = 20 # NCOSMIC+1 is maximum number of cosmic/outliers to be removed for each orders
#
# Wavelength calibration
......
......@@ -9,7 +9,7 @@ from drslib import db, rawfiles, darkframes, flatframes, wlframes, nodding, star
def gofio_init(docopt_args):
if docopt_args['--cfg']:
if docopt_args['--cfg'] or docopt_args['-g']:
try:
USRCFG = read_usr_config(docopt_args['-g'])
except:
......
......@@ -225,7 +225,7 @@ class GBNodding():
hea_ima = hea
for key in CONFIG['STRAIGHT_PAR']:
hea_ima[CONFIG['STRAIGHT_PAR'][key]] = imstr.header[CONFIG['STRAIGHT_PAR'][key]]
imflat = imstr.data
imflat = imstr.data.copy()
# read the number of images averaged to obtain the current image
# in order to compute the SNR
......@@ -350,6 +350,7 @@ class GBNodding():
heacal = OrderedDict()
#stdSpectrum = np.zeros((CONFIG['N_ORD'],CONFIG['YCCD']))
optSpectrum = np.zeros((CONFIG['N_ORD'],CONFIG['YCCD']))
snrSpectrum = np.zeros((CONFIG['N_ORD'],CONFIG['YCCD']))
fsnr = np.zeros((CONFIG['N_ORD'],CONFIG['YCCD']))
snr = np.zeros((CONFIG['N_ORD'],CONFIG['YCCD']))
......@@ -392,6 +393,15 @@ class GBNodding():
extlamp = varie.extract(goodlamp, optSpectrum[x], x1, x2, profile, lgaineff, lroneff)
osnr = imstr.data[start:end]
if slit_pos == CONFIG['B_POS']:
osnr = -osnr
ordersnr = np.ma.MaskedArray(osnr,mask=omask)
goodsnr = np.ma.compress_rows(ordersnr)
snrSpectrum[x] = varie.extract(goodsnr, optSpectrum[x], x1, x2, profile, gaineff, roneff)
if any(CONFIG['USE_FLAT'].values()) is True:
extflat = varie.extract(norflat, optSpectrum[x], x1, x2, profile, fgaineff, froneff)
#print extflat
......@@ -402,13 +412,30 @@ class GBNodding():
fsnr[fsnr==-np.inf] = 0
fsnr = np.nan_to_num(fsnr)
#print fsnr
#osnr = imstr.data[start:end]
#if slit_pos == CONFIG['B_POS']:
# osnr = -osnr
#ordersnr = np.ma.MaskedArray(osnr,mask=omask)
#goodsnr = np.ma.compress_rows(ordersnr)
#snrSpectrum[x] = varie.extract(goodsnr, optSpectrum[x], x1, x2, profile, gaineff, roneff)
with np.errstate(divide='ignore', invalid='ignore'):
#ssnr = (np.sqrt(roneff**2 + (gaineff*nspec*optSpectrum[x])))/(gaineff*nspec*optSpectrum[x])
ssnr = np.true_divide(np.sqrt(roneff**2 + (gaineff*snrSpectrum[x])),gaineff*snrSpectrum[x])
#ssnr = np.true_divide(np.sqrt(np.abs(x2-x1)*roneff**2 + (gaineff*snrSpectrum[x])),gaineff*snrSpectrum[x])
snr[x] = 1.0/(np.sqrt(ssnr**2 + fsnr[x]**2))
else:
fsnr[x] = np.zeros(len(optSpectrum[x]))
with np.errstate(divide='ignore', invalid='ignore'):
#ssnr = (np.sqrt(roneff**2 + (gaineff*nspec*optSpectrum[x])))/(gaineff*nspec*optSpectrum[x])
ssnr = np.true_divide(np.sqrt(roneff**2 + (gaineff*nspec*optSpectrum[x])),gaineff*nspec*optSpectrum[x])
snr[x] = 1.0/(np.sqrt(ssnr**2 + fsnr[x]**2))
with np.errstate(divide='ignore', invalid='ignore'):
#ssnr = np.true_divide(np.sqrt(roneff**2 + (gaineff*nspec*optSpectrum[x])),gaineff*nspec*optSpectrum[x])
ssnr = np.true_divide(np.sqrt(roneff**2 + (gaineff*snrSpectrum[x])) , gaineff*snrSpectrum[x])
snr[x] = 1.0/(np.sqrt(ssnr**2 + fsnr[x]**2))
calib_failed, coeffs, comments = varie.UNe_calibrate(extlamp,x+32,select_lines[x+32],all_lines[x+32])
......@@ -590,15 +617,15 @@ class GBNodding():
#print 's1d spectrum: %s s' % str(t2-t1)
if hea_ima[CONFIG['KEYS']['EXTMODE']] == CONFIG['EXTPAIR']:
return calname, fsnr, straight, dbreduced
return calname, fsnr, straight, dbreduced, snrSpectrum
elif 'grp' in fitsfile:
return calname, fsnr, straight, dbreduced
return calname, fsnr, straight, dbreduced, snrSpectrum
return calname, fsnr, False, dbreduced
return calname, fsnr, False, dbreduced, snrSpectrum
def combine(self,acalib,bcalib,fsnr):
def combine(self,acalib,bcalib,fsnr, asnrSpec, bsnrSpec):
dbreduced = {}
......@@ -651,10 +678,15 @@ class GBNodding():
#bshift = varie.rebin(awaves,bfluxes[o],bwaves,o)
bshift = varie.rebin(abhea,bfluxes[o],bwaves,o)
abcalib[o] = (afluxes[o]+bshift)/2.0
bsnr = varie.rebin(abhea,bsnrSpec[o],bwaves,o)
absnr = (asnrSpec[o]+bsnr)/2.0
#snr.append(max(np.mean(abcalib[o][1000:1050])/np.std(abcalib[o][1000:1050]),0))
with np.errstate(divide='ignore', invalid='ignore'):
#ssnr = (np.sqrt(roneff**2 + (gaineff*2*abcalib[o])))/(gaineff*2*abcalib[o])
ssnr = np.true_divide(np.sqrt(roneff**2 + (gaineff*2*abcalib[o])),gaineff*2*abcalib[o])
#ssnr = np.true_divide(np.sqrt(roneff**2 + (gaineff*2*abcalib[o])),gaineff*2*abcalib[o])
#ssnr = np.true_divide(np.sqrt(roneff**2 + (gaineff*abcalib[o])),gaineff*abcalib[o])
ssnr = np.true_divide(np.sqrt(roneff**2 + (gaineff*absnr)),gaineff*absnr)
snr[o] = 1.0/(np.sqrt(ssnr**2 + fsnr[o]**2))
abcalib = np.asarray(abcalib, dtype='float32')
......@@ -883,7 +915,7 @@ class GBNodding():
#print 'Bad pixels, create nodding: %s s' % str(t2-t1)
try:
acalib, fsnr, straight, dbreduced = self.reduce(ab,CONFIG['A_POS'], heaA)
acalib, fsnr, straight, dbreduced, asnrSpec = self.reduce(ab,CONFIG['A_POS'], heaA)
except:
straight = False
......@@ -898,7 +930,7 @@ class GBNodding():
pass
try:
bcalib, fsnr, straight, dbreduced = self.reduce(ab,CONFIG['B_POS'], heaB)
bcalib, fsnr, straight, dbreduced, bsnrSpec = self.reduce(ab,CONFIG['B_POS'], heaB)
except:
straight = False
......@@ -912,7 +944,7 @@ class GBNodding():
pass
try:
dbreduced = self.combine(acalib,bcalib,fsnr)
dbreduced = self.combine(acalib,bcalib,fsnr, asnrSpec, bsnrSpec)
except:
pass
......@@ -942,9 +974,9 @@ class GBNodding():
try:
if self.group['noddings'] and len(self.group['noddings'])>1 :
ab, heaA, heaB = self.group_avg()
acalib, fsnr, straight, dbreduced = self.reduce(ab,CONFIG['A_POS'], heaA)
bcalib, fsnr, straight, dbreduced = self.reduce(ab,CONFIG['B_POS'], heaB)
dbreduced = self.combine(acalib,bcalib,fsnr)
acalib, fsnr, straight, dbreduced, asnrSpec = self.reduce(ab,CONFIG['A_POS'], heaA)
bcalib, fsnr, straight, dbreduced, bsnrSpec = self.reduce(ab,CONFIG['B_POS'], heaB)
dbreduced = self.combine(acalib,bcalib,fsnr, asnrSpec, bsnrSpec)
acalib = None
bcalib = None
os.remove(ab)
......@@ -962,9 +994,9 @@ class GBNodding():
try:
if self.group['noddings'] and len(self.group['noddings'])>1:
ab, heaA, heaB = self.group_avg()
acalib, fsnr, straight, dbreduced = self.reduce(ab,CONFIG['A_POS'], heaA)
bcalib, fsnr, straight, dbreduced = self.reduce(ab,CONFIG['B_POS'], heaB)
dbreduced = self.combine(acalib,bcalib,fsnr)
acalib, fsnr, straight, dbreduced, asnrSpec = self.reduce(ab,CONFIG['A_POS'], heaA)
bcalib, fsnr, straight, dbreduced, bsnrSpec = self.reduce(ab,CONFIG['B_POS'], heaB)
dbreduced = self.combine(acalib,bcalib,fsnr, asnrSpec, bsnrSpec)
acalib = None
bcalib = None
os.remove(ab)
......
......@@ -393,6 +393,7 @@ class GBStare():
self.messages.append('%s: orders straightened.' % str(os.path.basename(fitsfile)),)
imstr = ccdproc.CCDData.read(straight, unit=u.adu)
imflat = imstr.data.copy()
try: nspec = imstr.header[CONFIG['KEYS']['NCOMBINE']]
except: nspec = 1
......@@ -537,7 +538,8 @@ class GBStare():
omask = gmask[start:end]
order = imstr.data[start:end]
#order = imstr.data[start:end]
order = imflat[start:end]
if CONFIG['USE_FLAT']['order']:
ordflat = flat.data[start:end]
......@@ -569,14 +571,32 @@ class GBStare():
fsnr[fsnr==np.inf] = 0
fsnr[fsnr==-np.inf] = 0
fsnr = np.nan_to_num(fsnr)
#print fsnr
osnr = imstr.data[start:end]
if slit_pos == CONFIG['B_POS']:
osnr = -osnr
ordersnr = np.ma.MaskedArray(osnr,mask=omask)
goodsnr = np.ma.compress_rows(ordersnr)
snrSpectrum = varie.extract(goodsnr, optSpectrum[x], x1, x2, profile, gaineff, roneff)
with np.errstate(divide='ignore', invalid='ignore'):
#ssnr = (np.sqrt(roneff**2 + (gaineff*nspec*optSpectrum[x])))/(gaineff*nspec*optSpectrum[x])
#ssnr = np.true_divide(np.sqrt(roneff**2 + (gaineff*nspec*snrSpectrum)),gaineff*nspec*snrSpectrum)
ssnr = np.true_divide(np.sqrt(roneff**2 + (gaineff*snrSpectrum)),gaineff*snrSpectrum)
snr[x] = 1.0/(np.sqrt(ssnr**2 + fsnr[x]**2))
else:
fsnr[x] = np.zeros(len(optSpectrum[x]))
#fsnr[x] = np.zeros(len(optSpectrum[x]))
with np.errstate(divide='ignore', invalid='ignore'):
#ssnr = np.true_divide(np.sqrt(roneff**2 + (gaineff*nspec*optSpectrum[x])),gaineff*nspec*optSpectrum[x])
ssnr = np.true_divide(np.sqrt(roneff**2 + (gaineff*optSpectrum[x])),gaineff*optSpectrum[x])
snr[x] = 1.0/(np.sqrt(ssnr**2 + fsnr[x]**2))
with np.errstate(divide='ignore', invalid='ignore'):
#ssnr = (np.sqrt(roneff**2 + (gaineff*nspec*optSpectrum[x])))/(gaineff*nspec*optSpectrum[x])
ssnr = np.true_divide(np.sqrt(roneff**2 + (gaineff*nspec*optSpectrum[x])),gaineff*nspec*optSpectrum[x])
snr[x] = 1.0/(np.sqrt(ssnr**2 + fsnr[x]**2))
#with np.errstate(divide='ignore', invalid='ignore'):
# ssnr = np.true_divide(np.sqrt(roneff**2 + (gaineff*nspec*optSpectrum[x])),gaineff*nspec*optSpectrum[x])
# snr[x] = 1.0/(np.sqrt(ssnr**2 + fsnr[x]**2))
calib_failed, coeffs, comments = varie.UNe_calibrate(extlamp,x+32,select_lines[x+32],all_lines[x+32])
......
......@@ -303,6 +303,7 @@ def optExtract(data,gain,ron,slit_pos,ordine):
#plt.show()
#plt.plot(profile)
#plt.plot(profile[x1:x2])
#plt.show()
# enforce normalization
......@@ -345,22 +346,32 @@ def optExtract(data,gain,ron,slit_pos,ordine):
model = OptFlux*profile
variance = var(model)
#variance = var(model[x1:x2])
#plt.plot(np.sum(data[x1:x2], axis=0),'k:')
#plt.plot(OptFlux,'b-.')
#plt.plot(np.sum(model[x1:x2], axis=0),'r--')
#plt.show()
stop = 0
cosmic = True
while cosmic:
with np.errstate(divide='ignore', invalid='ignore', over='ignore'):
outliers = np.true_divide((data-model)**2,np.abs(variance))
#outliers[outliers==np.inf] = 50
#outliers = np.true_divide((data-model)**2,np.abs(variance))
outliers = np.true_divide((data[x1:x2]-model[x1:x2])**2,np.abs(variance[x1:x2]))
outliers[np.isnan(outliers)] = 50
#if outliers[np.isnan(outliers)]:
# print len(outliers[np.isnan(outliers)])
#outliers = np.nan_to_num(outliers)
if np.amax(outliers[x1:x2]) > 25:
rworst = np.unravel_index(np.argmax(outliers[x1:x2]),outliers[x1:x2].shape)[0] + x1
cworst = np.unravel_index(np.argmax(outliers[x1:x2]),outliers[x1:x2].shape)[1]
#if np.amax(outliers[x1:x2]) > 25:
# rworst = np.unravel_index(np.argmax(outliers[x1:x2]),outliers[x1:x2].shape)[0] + x1
# cworst = np.unravel_index(np.argmax(outliers[x1:x2]),outliers[x1:x2].shape)[1]
if np.amax(outliers) > 25:
rworst = np.unravel_index(np.argmax(outliers),outliers.shape)[0] + x1
cworst = np.unravel_index(np.argmax(outliers),outliers.shape)[1]
data[rworst,cworst] = np.median([data[rworst,max(cworst-4,0):min(cworst+5,CONFIG['XCCD'])]])
profile[rworst,cworst] = np.median([profile[rworst,max(cworst-4,0):min(cworst+5,CONFIG['XCCD'])]])
......@@ -383,7 +394,8 @@ def optExtract(data,gain,ron,slit_pos,ordine):
variance = var(model)
stop += 1
if stop > 20:
# if stop > 20:
if stop > CONFIG['NCOSMIC']:
cosmic = False
else:
......
......@@ -3,10 +3,10 @@ GOFIO DRS
Written by Monica Rainer
Usage:
gofioDRS.py -h
gofioDRS.py [--cfg=<config>]
gofioDRS.py <date> [--cfg=<config>]
gofioDRS.py <date> <calib_date> [--cfg=<config>]
gofioDRS.py -h | --help
gofioDRS.py [(-g <config> | --cfg=<config>)]
gofioDRS.py <date> [(-g <config> | --cfg=<config>)]
gofioDRS.py <date> <calib_date> [(-g <config> | --cfg=<config>)]
gofioDRS.py <date> [--dark --flat --une --fp --all_calib --only_calib --use_flat=<flag> --s1d=<s1d> --group]
gofioDRS.py <date> <calib_date> [--dark --flat --une --fp --all_calib --only_calib --use_flat=<flag> --s1d=<s1d> --group]
......@@ -14,6 +14,8 @@ Options:
-h,--help : show this screen
--cfg=<config> : path of the configuration file (optional)
'filepath/filename'
-g <config> : path of the configuration file (optional)
'filepath/filename'
date : date to be reduced (always first input)
calib_date : calibration date to be used (always second input)
--dark : darks are reduced
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment