Commit 9052ca74 authored by Monica Rainer's avatar Monica Rainer
Browse files

Modified MASKC and added another wavelength calibration method

parent 4193320d
backup
calibrations
gofio_latest
gofio_latest.tar
reduced
varie
vradlib
......
......@@ -138,9 +138,6 @@ CONFIG['KEY_DRS'] = 'HIERARCH TNG DRS'
CONFIG['RON_EFF'] = (' '.join((CONFIG['KEY_DRS'],'RON')),'Effective RON after image reduction')
CONFIG['GAIN_EFF'] = (' '.join((CONFIG['KEY_DRS'],'GAIN')),'Effective gain after image reduction')
CONFIG['TEXP_EFF'] = (' '.join((CONFIG['KEY_DRS'],'EXPTIME')),'Effective exposure time after image reduction')
CONFIG['WLFIT'] = (' '.join((CONFIG['KEY_DRS'],'CAL WLFUNC')),'Function used for wavelength calibration')
CONFIG['WLFIT_FUNC'] = 'l0 + k1*(x-xc) + k2*(x-xc)**2 + k3*(x-xc)**3'
CONFIG['WLCOEFFS'] = {'k1':(' '.join((CONFIG['KEY_DRS'],'CAL K1')),'k1 coefficient of the calibration function'),'k2':(' '.join((CONFIG['KEY_DRS'],'CAL K2')),'k2 coefficient of the calibration function'),'k3':(' '.join((CONFIG['KEY_DRS'],'CAL K3')),'k3 coefficient of the calibration function'),'l0':(' '.join((CONFIG['KEY_DRS'],'CAL L0')),'l0 coefficient of the calibration function'),'xc':(' '.join((CONFIG['KEY_DRS'],'CAL XC')),'xc coefficient of the calibration function'),'rms':(' '.join((CONFIG['KEY_DRS'],'CAL RMSE')),'[m/s] RMS error of the calibration')}
CONFIG['SPEC_USED'] = (' '.join((CONFIG['KEY_DRS'],'SPEC')),'Images combined')
CONFIG['SPEC_MJD'] = (' '.join((CONFIG['KEY_DRS'],'MJD')),'MJD of the relative image')
CONFIG['SNR'] = (' '.join((CONFIG['KEY_DRS'],'SNR')),'SNR in the middle of the order',)
......@@ -152,6 +149,18 @@ CONFIG['AIRMASS'] = (' '.join((CONFIG['KEY_DRS'],'AIRMASS')),'Airmass of the com
CONFIG['MASTERFLAT'] = (' '.join((CONFIG['KEY_DRS'],'FLAT')),'Flat-field used for the reduction')
CONFIG['MASTERLAMP'] = (' '.join((CONFIG['KEY_DRS'],'LAMP')),'Calibration lamp used for the reduction')
# Wavelength calibration
CONFIG['CAL_FUNC'] = {'Oliva':True, 'Poly3':False} # function used for wavelength calibration (Tino Oliva Polynomial or 3-degree polynomial)
CONFIG['WLFIT'] = (' '.join((CONFIG['KEY_DRS'],'CAL WLFUNC')),'Function used for wavelength calibration')
if CONFIG['CAL_FUNC']['Oliva']:
CONFIG['WLFIT_FUNC'] = 'l0 + k1*(x-xc) + k2*(x-xc)**2 + k3*(x-xc)**3'
#CONFIG['WLCOEFFS'] = {'k1':(' '.join((CONFIG['KEY_DRS'],'CAL K1')),'k1 coefficient of the calibration function'),'k2':(' '.join((CONFIG['KEY_DRS'],'CAL K2')),'k2 coefficient of the calibration function'),'k3':(' '.join((CONFIG['KEY_DRS'],'CAL K3')),'k3 coefficient of the calibration function'),'l0':(' '.join((CONFIG['KEY_DRS'],'CAL L0')),'l0 coefficient of the calibration function'),'xc':(' '.join((CONFIG['KEY_DRS'],'CAL XC')),'xc coefficient of the calibration function'),'rms':(' '.join((CONFIG['KEY_DRS'],'CAL RMSE')),'[m/s] RMS error of the calibration')}
else:
CONFIG['WLFIT_FUNC'] = 'c0 + c1*x + c2*x**2 + c3*x**3'
CONFIG['WLCOEFFS'] = {'k1':(' '.join((CONFIG['KEY_DRS'],'CAL K1')),'k1 coefficient of the calibration function'),'k2':(' '.join((CONFIG['KEY_DRS'],'CAL K2')),'k2 coefficient of the calibration function'),'k3':(' '.join((CONFIG['KEY_DRS'],'CAL K3')),'k3 coefficient of the calibration function'),'l0':(' '.join((CONFIG['KEY_DRS'],'CAL L0')),'l0 coefficient of the calibration function'),'xc':(' '.join((CONFIG['KEY_DRS'],'CAL XC')),'xc coefficient of the calibration function'),'rms':(' '.join((CONFIG['KEY_DRS'],'CAL RMSE')),'[m/s] RMS error of the calibration'), 'c0':(' '.join((CONFIG['KEY_DRS'],'CAL C0')),'c0 coefficient of the calibration function'),'c1':(' '.join((CONFIG['KEY_DRS'],'CAL C1')),'c1 coefficient of the calibration function'),'c2':(' '.join((CONFIG['KEY_DRS'],'CAL C2')),'c2 coefficient of the calibration function'),'c3':(' '.join((CONFIG['KEY_DRS'],'CAL C3')),'c3 coefficient of the calibration function'),'rms_poly':(' '.join((CONFIG['KEY_DRS'],'CAL RMSE')),'[m/s] RMS error of the calibration')}
#
# DB configuration
#
......@@ -192,6 +201,7 @@ CONFIG['B_POS'] = 23. # guess pixel position for B
CONFIG['C_POS'] = 15. # guess pixel position for C
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
#
......
......@@ -25,6 +25,7 @@ from astropy import coordinates as coord
from astropy.io import fits
import numpy as np
import numpy.polynomial.polynomial as poly
import math
import warnings
......@@ -238,7 +239,8 @@ def optExtract(data,gain,ron,slit_pos,ordine):
stop = 0
outlier = True
while outlier:
fitprofile = np.polyval(np.polyfit(np.arange(columns),profile[row],deg=2,w=1./np.sqrt(variance[row])),np.arange(columns))
#fitprofile = np.polyval(np.polyfit(np.arange(columns),profile[row],deg=2,w=1./np.sqrt(variance[row])),np.arange(columns))
fitprofile = poly.polyval(np.arange(columns),poly.polyfit(np.arange(columns),profile[row],deg=2,w=1./np.sqrt(variance[row])))
sigma = np.mean((profile[row]-fitprofile)**2)
......@@ -462,7 +464,7 @@ def UNe_linelist():
#--------------------- UNe calibration -------------------
def UNe_calibrate(lamp,order,select_lines,all_lines):
def UNe_calibrate(lamp,order,select_lines,all_lines,use_oliva=CONFIG['CAL_FUNC']['Oliva'],use_poly=CONFIG['CAL_FUNC']['Poly3']):
warnings.simplefilter('error',RuntimeWarning)
warnings.simplefilter("ignore", optimize.OptimizeWarning)
......@@ -650,7 +652,10 @@ def UNe_calibrate(lamp,order,select_lines,all_lines):
# if fitting, check that the values are reasonable
# then apply the calibration to the whole pixel range
coeffs = OrderedDict()
if used > 1:
xc0 = CONFIG['XC_GUESS'][order]
p0 = (lambda0,xc0)
pars, pcov = optimize.curve_fit(lambdafit,all_pixels,all_waves,p0)
......@@ -659,12 +664,15 @@ def UNe_calibrate(lamp,order,select_lines,all_lines):
#print order
#print ' **** FIT order %s : xc0 fitted - xc0 tabulated: %s' % (str(order),str(xc0-CONFIG['XC_GUESS'][order]))
#fitpoly3 = np.polyfit(all_pixels,all_waves,deg=3)
if abs(xc0-CONFIG['XC_GUESS'][order]) > CONFIG['WAVE_FIT']['xc_range']:
#print ' **** FIT def. order %s : xc0 fitted - xc0 tabulated: %s' % (str(order),str(xc0-CONFIG['XC_GUESS'][order]))
#calib = np.zeros(len(pixrange))
messages.append('Calibration failed for order %s.' % (str(order),))
coeffs = OrderedDict({'k1':None,'k2':None,'k3':None,'l0':None,'xc':None,'rms':None})
if use_oliva:
#coeffs = OrderedDict({'k1':None,'k2':None,'k3':None,'l0':None,'xc':None,'rms':None})
coeffs.update({'k1':None,'k2':None,'k3':None,'l0':None,'xc':None,'rms':None})
calib_failed = True
#coeffs = OrderedDict()
......@@ -682,15 +690,25 @@ def UNe_calibrate(lamp,order,select_lines,all_lines):
rvrmse=np.sqrt(rvchisq/dof)*const.c
#messages.append('RMS of the calibration for order %s: %s' % (str(order),str(rvrmse),))
#print 'RMS of the calibration for order %s: %s' % (str(order),str(rvrmse),)
coeffk1 = float('%.5e' % k1)
coeffk2 = float('%.5e' % k2)
coeffk3 = float('%.5e' % k3)
coeffl0 = round(lambda0,5)
coeffxc0 = round(xc0,5)
#coeffk1 = float('%.5e' % k1)
#coeffk2 = float('%.5e' % k2)
#coeffk3 = float('%.5e' % k3)
#coeffl0 = round(lambda0,5)
#coeffxc0 = round(xc0,5)
#coeffrms = round(rvrmse.value,2)
coeffk1 = float(k1)
coeffk2 = float( k2)
coeffk3 = float(k3)
coeffl0 = lambda0
coeffxc0 = xc0
coeffrms = round(rvrmse.value,2)
coeffs = OrderedDict({'k1':coeffk1,'k2':coeffk2,'k3':coeffk3,'l0':coeffl0,'xc':coeffxc0,'rms':coeffrms})
if use_oliva:
#coeffs = OrderedDict({'k1':coeffk1,'k2':coeffk2,'k3':coeffk3,'l0':coeffl0,'xc':coeffxc0,'rms':coeffrms})
coeffs.update({'k1':coeffk1,'k2':coeffk2,'k3':coeffk3,'l0':coeffl0,'xc':coeffxc0,'rms':coeffrms})
calib_failed = False
#calib = lambdafit(pixrange,lambda0,xc0)
#coeffs = OrderedDict({'k1':k1,'k2':k2,'k3':k3,'l0':lambda0,'xc':xc0,'rms':rvrmse.value})
......@@ -698,20 +716,52 @@ def UNe_calibrate(lamp,order,select_lines,all_lines):
#plt.plot(np.arange(len(pixrange)),lambdafit(np.arange(len(pixrange))),'r-')
#plt.show()
#calib = np.polyval(fitpoly3,pixrange)
if use_poly:
fitpoly3 = poly.polyfit(all_pixels,all_waves,deg=3)
#fitpoly3 = np.polyfit(all_pixels,all_waves,deg=3)
#polycalib = np.polyval(fitpoly3,pixrange)
#check_calib = np.mean(np.absolute(polycalib-calib))
#if check_calib < CONFIG['CHECK_CALIB']:
chisq=((poly.polyval(all_pixels,fitpoly3)-all_waves)**2).sum()
rvchisq=(((poly.polyval(all_pixels,fitpoly3)-all_waves)/all_waves)**2).sum()
#chisq=((np.polyval(fitpoly3,all_pixels)-all_waves)**2).sum()
#rvchisq=(((np.polyval(fitpoly3,all_pixels)-all_waves)/all_waves)**2).sum()
#dof=len(all_pixels)-4
#rmse=np.sqrt(chisq/dof)
#rvrmse=np.sqrt(rvchisq/dof)*const.c
dof=max(len(all_pixels)-4,1)
rmse=np.sqrt(chisq/dof)
rvrmse=np.sqrt(rvchisq/dof)*const.c
#print 'RMS of the calibration for order %s: %s' % (str(order),str(rvrmse),)
#coeffc0 = round(fitpoly3[0],5)
#coeffc1 = round(fitpoly3[1],5)
#coeffc2 = round(fitpoly3[2],5)
#coeffc3 = round(fitpoly3[3],5)
#coeffrms = round(rvrmse.value,2)
coeffc0 = fitpoly3[0]
coeffc1 = fitpoly3[1]
coeffc2 = fitpoly3[2]
coeffc3 = fitpoly3[3]
coeffrms = round(rvrmse.value,2)
calib_failed = False
#coeffs = OrderedDict({'c0':coeffc0,'c1':coeffc1,'c2':coeffc2,'c3':coeffc3,'rms':coeffrms})
coeffs.update({'c0':coeffc0,'c1':coeffc1,'c2':coeffc2,'c3':coeffc3,'rms_poly':coeffrms})
else:
print ' **** WARNING **** Order %s: not enough lines for the calibration!' % (str(order))
#calib = np.zeros(len(lamp))
messages.append('Calibration failed for order %s, not enough lines.' % (str(order),))
coeffs = OrderedDict({'k1':None,'k2':None,'k3':None,'l0':None,'xc':None,'rms':None})
if use_oliva:
#coeffs = OrderedDict({'k1':None,'k2':None,'k3':None,'l0':None,'xc':None,'rms':None})
coeffs.append({'k1':None,'k2':None,'k3':None,'l0':None,'xc':None,'rms':None})
if use_poly:
#coeffs = OrderedDict({'c0':None,'c1':None,'c2':None,'c3':None,'rms':None})
coeffs.append({'c0':None,'c1':None,'c2':None,'c3':None,'rms_poly':None})
calib_failed = True
......@@ -741,6 +791,11 @@ def UNe_calibrate(lamp,order,select_lines,all_lines):
#--------------------- Apply wavelength calibration -------------------
def wcalib(heawave,o):
#pixrange = -np.arange(CONFIG['YCCD'])+CONFIG['YCCD']
pixrange = np.arange(CONFIG['YCCD'])+1
if CONFIG['CAL_FUNC']['Oliva']:
def lambdafit(x,lambda0,xc0):
return lambda0 + k1*(x-xc0) + k2*(x-xc0)**2 + k3*(x-xc0)**3
......@@ -755,10 +810,28 @@ def wcalib(heawave,o):
keyxc = ''.join((CONFIG['WLCOEFFS']['xc'][0],str(o+32)))
xc = float(heawave[keyxc])
#pixrange = -np.arange(CONFIG['YCCD'])+CONFIG['YCCD']
pixrange = np.arange(CONFIG['YCCD'])+1
wave = lambdafit(pixrange,l0,xc)
else:
#def lambdafit(x):
# return c0 + c1*x + c2*(x**2) + c3*(x**3)
keyc0 = ''.join((CONFIG['WLCOEFFS']['c0'][0],str(o+32)))
c0 = float(heawave[keyc0])
keyc1 = ''.join((CONFIG['WLCOEFFS']['c1'][0],str(o+32)))
c1 = float(heawave[keyc1])
keyc2 = ''.join((CONFIG['WLCOEFFS']['c2'][0],str(o+32)))
c2 = float(heawave[keyc2])
keyc3 = ''.join((CONFIG['WLCOEFFS']['c3'][0],str(o+32)))
c3 = float(heawave[keyc3])
fitpoly3 = np.array([c0,c1,c2,c3])
wave = poly.polyval(pixrange,fitpoly3)
#wave = np.polyval(fitpoly3,pixrange)
#wave = lambdafit(pixrange)
return wave
......@@ -1054,6 +1127,7 @@ def create_s1d(spectrum, header):
#wfit = -99999
flux_old = spectrum[o][::-1]
if CONFIG['S1D_NORM']:
# use the central part of the order to normalize it
chunk = int(len(flux_old)/4)
norm = flux_old[chunk:(chunk*3)]
......
......@@ -51,7 +51,7 @@ class GBWls():
def masterlamp(self,unefp):
if self.wllist:
exptime = self.wllist[0].header[CONFIG['KEYS']['EXPTIME']]
exptime = float(self.wllist[0].header[CONFIG['KEYS']['EXPTIME']])
darkname = 'dark' + str(int(exptime))
use_dark = True
......@@ -185,7 +185,7 @@ class GBWls():
ordini[x] = np.sum(wldata[row-5:row+5],axis=0)
if unefp == 'une':
calib_failed, coeffs, comments = varie.UNe_calibrate(ordini[x],x+32,select_lines[x+32],all_lines[x+32])
calib_failed, coeffs, comments = varie.UNe_calibrate(ordini[x],x+32,select_lines[x+32],all_lines[x+32],True,True)
for comment in comments:
self.messages.append(comment)
......
......@@ -29,6 +29,10 @@ try:
CONFIG['LOG_FILE'] = os.path.join(CONFIG['RED_DIR'], 'drs.log')
#CONFIG['DO_CALIB'] = {'dark':True,'flat':True,'une':True,'fp':True,'only_calib':True}
CONFIG['DO_CALIB'] = {'dark':False,'flat':False,'une':False,'fp':False,'only_calib':False}
maskc = ''.join(('GIANOB_MASKC_',CONFIG['DATE'],'.fits'))
CONFIG['MASK_C'] = os.path.join(CONFIG['RES_DIR'], maskc)
except:
pass
......@@ -370,6 +374,8 @@ if __name__ == "__main__":
if db.check_raw(dbramps,raws[n]):
GianoBreduce(raws[n], rawlists, dbcalib, working, group, do_calib, stop)
break
try: os.remove(CONFIG['MASK_C'])
except: pass
drslogger.log('Pipeline stopped.')
sys.exit(0)
dbcalib.close()
......
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