nodding.py 44.1 KB
Newer Older
Monica Rainer's avatar
Monica Rainer committed
1
2
3
4
5
6
7
8
9
10
11
12
"""
Reduction of the single AB nodding:
- check the image quality (signal in well defined region)
- check that the exposure times of A and B are the same
- remove bad pixel using the bad pixel mask
- create nodding A-B
- straigthen the nodding
- divide by the masterflat (if required)
- perform optimal extraction of A-B and B-A
- use the optimal profiles with UNe lamp and do wavelength calibration of A-B and B-A
- save A, B and the average of A+B
Reduction of the nodding group:
Monica Rainer's avatar
   
Monica Rainer committed
13
- average the A-B noddings
Monica Rainer's avatar
Monica Rainer committed
14
15
16
17
18
19
20
21
22
23
- straigthen the average nodding
- divide by the masterflat (if required)
- perform optimal extraction of A-B and B-A
- use the optimal profiles with UNe lamp and do wavelength calibration of A-B and B-A
- save A, B and the average of A+B
"""


from drslib.config import CONFIG
from drslib import db, varie
Andrea Bignamini's avatar
Andrea Bignamini committed
24
from drslib import metadata
Monica Rainer's avatar
Monica Rainer committed
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40

from astropy import units as u
from astropy.io import ascii, fits

import warnings
from astropy.utils.exceptions import AstropyWarning
import ccdproc

import numpy as np
import math, os, subprocess, time, shutil

from collections import OrderedDict
#import matplotlib.pyplot as plt


class GBNodding():
Monica Rainer's avatar
Monica Rainer committed
41
    def __init__(self, nodding, group, dbconn, dbnight):
Monica Rainer's avatar
Monica Rainer committed
42
43
44
        self.nodding = nodding
        self.group = group
        self.dbconn = dbconn
Monica Rainer's avatar
Monica Rainer committed
45
        self.dbnight = dbnight
Monica Rainer's avatar
Monica Rainer committed
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
        self.quality = []
        self.messages = []
        self.nodlist = []
        self.nodcorr = {}

    def qualitycheck(self):
        """
        Check image's quality: check the signal in a well-defined region.
        If one of the nodding images fails the check, the other will be discarded too.
        Check that the exposure time of the nodding is the same for the two images.
        """

        expt = []
        for frame in self.nodding:
            #print frame

            nod = ccdproc.CCDData.read(frame, unit=u.adu)
            expt.append(nod.header[CONFIG['KEYS']['EXPTIME']])

            # check the signal in a well-defined zone

            zone = nod.data[CONFIG['SCIENCECHECK'][0]:CONFIG['SCIENCECHECK'][1],CONFIG['SCIENCECHECK'][2]:CONFIG['SCIENCECHECK'][3]]
            mean = np.mean(zone)

            if mean < CONFIG['NODSIGNAL']:
                self.messages.append('Science frame %s failed quality check: signal too low (%s). Neither of the nodding images (%s and %s) will be reduced.' % (str(os.path.basename(frame)),str(mean),str(os.path.basename(self.nodding[0])),str(os.path.basename(self.nodding[1]))))
                nod = None

                return False

            else:
                self.nodlist.append(nod)

        if expt[0] != expt[1]:
            self.messages.append('The two nodding images %s and %s have different exposure times and they can not be reduced.' % (str(os.path.basename(self.nodding[0])),str(os.path.basename(self.nodding[1]))))
            return False

        return True


86
    def create_nodcorr(self):
Monica Rainer's avatar
Monica Rainer committed
87
        """
88
        Create nodcorr
Monica Rainer's avatar
Monica Rainer committed
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
        """
        badpix = ccdproc.CCDData.read(CONFIG['BADPIX_MASK'], unit=u.adu)
        bad_mask=badpix.data
        inverse_mask=np.logical_not(bad_mask)

        for nod in self.nodlist:

            # mask the image using the badpix mask
            nod_corrected = varie.badpix(nod.data,bad_mask,inverse_mask)
            self.messages.append('Bad pixel correction done.')

            nod_corr = ccdproc.CCDData(nod_corrected,unit=u.adu)
            nod_corr.header = nod.header

            self.nodcorr[nod.header[CONFIG['KEYS']['SLIT']]] = nod_corr
            nod = None

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
        return


    def check_nodcorr(self):
        """
        Check if nodcorr is not corrupted
        """

        try:
            nodA = self.nodcorr[CONFIG['A']]
        except Exception as e:
            nodA = None

        try:
            nodB = self.nodcorr[CONFIG['B']]
        except Exception as e:
            nodB = None

        if not nodA and not nodB:
            self.messages.append('The nodding pair is corrupted. Both slits A and B are missing. Skipping pair processing.')
            return False
        else:
            if not nodA:
                self.messages.append('The nodding pair is corrupted. Slit A is missing. Skipping pair processing.')
                return False
            if not nodB:
                self.messages.append('The nodding pair is corrupted. Slit B is missing. Skipping pair processing.')
                return False

        return True


    def createAB(self):
        """
        Create nodding images A-B and save them in temporary directory if
        the keyword SPEXTMODE is set to GRPAVG_EXT.
        """

        hea = self.nodlist[0].header

Monica Rainer's avatar
Monica Rainer committed
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
        # Create A-B image

        nodA = self.nodcorr[CONFIG['A']]
        nodB = self.nodcorr[CONFIG['B']]
        AB = nodA.data - nodB.data
        nodAB = ccdproc.CCDData(AB, unit=u.adu)
        nodAB.data = np.asarray(nodAB.data, dtype='float32')
        nodAB.header = hea
        heaA = nodA.header
        heaB = nodB.header
        #qui = str(os.path.basename(self.nodding[0])).rindex('.')
        #ABnome = '_'.join((str(os.path.basename(self.nodding[0]))[0:qui], 'AB.fits'))
        nomebase = os.path.splitext(os.path.basename(self.nodding[0]))[0]
        ABnome = '_'.join((nomebase, 'AB.fits'))

        keyA = ''.join((CONFIG['SPEC_USED'][0],'1','A'))
        mjdA = ''.join((CONFIG['SPEC_MJD'][0],'1','A'))
        keyB = ''.join((CONFIG['SPEC_USED'][0],'1','B'))
        mjdB = ''.join((CONFIG['SPEC_MJD'][0],'1','B'))

        for h in (heaA,heaB,nodAB.header):
            h[CONFIG['GAIN_EFF'][0]] = (CONFIG['GAIN'],CONFIG['GAIN_EFF'][1])
            h[CONFIG['RON_EFF'][0]] = (CONFIG['RON']*math.sqrt(2),CONFIG['RON_EFF'][1])
            h[CONFIG['KEYS']['NCOMBINE']] = 2
            h[keyA] = (nodA.header[CONFIG['KEYS']['IMANAME']],CONFIG['SPEC_USED'][1])
            h[mjdA] = (nodA.header[CONFIG['KEYS']['MJD']],CONFIG['SPEC_MJD'][1])
            h[keyB] = (nodB.header[CONFIG['KEYS']['IMANAME']],CONFIG['SPEC_USED'][1])
            h[mjdB] = (nodB.header[CONFIG['KEYS']['MJD']],CONFIG['SPEC_MJD'][1])

        nodAB.header[CONFIG['KEYS']['FILENAME']] = ABnome
176
        nodAB.header[CONFIG['KEYS']['IMANAME']] = ABnome
Monica Rainer's avatar
Monica Rainer committed
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
        ABtmpnome = os.path.join(CONFIG['TMP_DIR'],ABnome)


        try:
            nodAB.header[CONFIG['KEYS']['EXTMODE']]
        except:
            nodAB.header[CONFIG['KEYS']['EXTMODE']] = CONFIG['EXTDEFAULT']
            heaA[CONFIG['KEYS']['EXTMODE']] = CONFIG['EXTDEFAULT']
            heaB[CONFIG['KEYS']['EXTMODE']] = CONFIG['EXTDEFAULT']

        if nodAB.header[CONFIG['KEYS']['EXTMODE']] == CONFIG['EXTAVG']:
            self.group['noddings'].append(ABtmpnome)


        hdu = fits.PrimaryHDU(data=nodAB.data,header=nodAB.header)
        nodABfits = fits.HDUList([hdu])
        nodABfits.writeto(ABtmpnome,clobber=True)

        self.messages.append('A-B nodding created.')
        nodA = None
        nodB = None

        return ABtmpnome, heaA, heaB

    def reduce(self,fitsfile,slit_pos,hea):
        """
        Straighten, divide by the masterflat, optimal extraction
        """
Monica Rainer's avatar
Monica Rainer committed
205
        dbreduced = {}
Monica Rainer's avatar
Monica Rainer committed
206
207
208
209
210
211
212
213
214

        if slit_pos == CONFIG['A_POS']:
            slit = 'A'
        else:
            slit = 'B'

        # straighten

        straight = fitsfile.replace('.fits','_str.fits')
215

216
217
218
219
220
        args = [CONFIG['STRAIGHT'],fitsfile,straight]
        args.extend(CONFIG['STRAIGHT_OPT'])
        # search for shift defined in the straighten options in config.py
        dy = True
        for opt in CONFIG['STRAIGHT_OPT']:
221
            try:
222
223
224
225
                dy = opt.rindex('DY=')
                ypos = int(opt[dy-2:])
                shift = CONFIG['SHIFT_Y'] + ypos
                dy = False
226
            except:
227
228
229
                pass

        if dy:
Unknown's avatar
Unknown committed
230
231
232
233
234
            try:
                shift = db.extract_dbfile(self.dbconn,'shiftY')
            except:
                try:
                    cal_flat = db.extract_dbfile(self.dbconn,'flat')
235
236
237
                    mflat = ccdproc.CCDData.read(cal_flat, unit=u.adu)
                    shift = varie.shiftY(mflat.data)
                    db.insert_dbfile(self.dbconn,'shiftY',shift)
Unknown's avatar
Unknown committed
238
                except:
Monica Rainer's avatar
   
Monica Rainer committed
239
240
241
242
243
244
245
                    try:
                        db.copy_dbfile(self.dbconn,'shiftY')
                        shift = db.extract_dbfile(self.dbconn,'shiftY')
                    except:
                        shift = CONFIG['SHIFT_Y']
                        self.messages.append('No flat-field or shift value present in the calibration database, no shift will be applied.')

246
247
            if not shift:
                shift = CONFIG['SHIFT_Y']
248

249
250
251
            shiftY = [''.join(('DY=',str(shift - CONFIG['SHIFT_Y'])))]
            #print shiftY
            args.extend(shiftY)
252

Monica Rainer's avatar
Monica Rainer committed
253
        subprocess.call(args)
254

Andrea Bignamini's avatar
Andrea Bignamini committed
255
        # Read straight file
Monica Rainer's avatar
   
Monica Rainer committed
256
        #imstr = ccdproc.CCDData.read(straight, unit=u.adu)
Andrea Bignamini's avatar
Andrea Bignamini committed
257
258
259
260

        # Update FILENAME in header then
        # add metadata to header and save straight file

Monica Rainer's avatar
   
Monica Rainer committed
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277

        with fits.open(straight, 'update') as s:
            s[0].header[CONFIG['KEYS']['FILENAME']] = os.path.basename(straight)
            s[0].header = metadata.add_metadata(s[0].header)
            s.flush()




        #imstr.header[CONFIG['KEYS']['FILENAME']] = os.path.basename(straight)
        #imstr.header = metadata.add_metadata(imstr.header)

        #hdu = fits.PrimaryHDU(data=imstr.data, header=imstr.header)
        #str_fits = fits.HDUList([hdu])
        #str_fits.writeto(straight, overwrite=True)

        imstr = ccdproc.CCDData.read(straight, unit=u.adu)
Andrea Bignamini's avatar
Andrea Bignamini committed
278

Monica Rainer's avatar
Monica Rainer committed
279
280
281
282
283
284
285
286
287
        str_file = os.path.join(CONFIG['RED_STR'],os.path.basename(straight))
        try: shutil.copyfile(straight,str_file)
        except: pass

        self.messages.append('%s: orders straightened (nodding %s).' % (str(os.path.basename(fitsfile)),slit,))

        hea_ima = hea
        for key in CONFIG['STRAIGHT_PAR']:
            hea_ima[CONFIG['STRAIGHT_PAR'][key]] = imstr.header[CONFIG['STRAIGHT_PAR'][key]]
Monica Rainer's avatar
Monica Rainer committed
288
        imflat = imstr.data.copy()
Monica Rainer's avatar
Monica Rainer committed
289
290
291
292
293
294
295
296
297
298
299
300
301

        # read the number of images averaged to obtain the current image
        # in order to compute the SNR
        try: nspec = hea_ima[CONFIG['KEYS']['NCOMBINE']]
        except: nspec = 1

        # use only the regions of the orders
        try:
            goodmask = ccdproc.CCDData.read(CONFIG['MASK_C'], unit=u.adu)
        except:
            try:
                masterflat = db.extract_dbfile(self.dbconn,'flatstr')
            except:
Monica Rainer's avatar
Monica Rainer committed
302
303
                masterflat = False
            if not masterflat:
Monica Rainer's avatar
   
Monica Rainer committed
304
305
306
307
308
309
310
#                db.copy_dbfile(self.dbconn,'flatstr')
#                masterflat = db.extract_dbfile(self.dbconn,'flatstr')
#                self.messages.append('No masterflat found for this night, it will be taken from the calibration database: %s' % (os.path.basename(masterflat)))
                try:
                    db.copy_dbfile(self.dbconn,'flatstr')
                    masterflat = db.extract_dbfile(self.dbconn,'flatstr')
                    self.messages.append('No masterflat found for this night, it will be taken from the calibration database: %s' % (os.path.basename(masterflat)))
Monica Rainer's avatar
Monica Rainer committed
311

Monica Rainer's avatar
   
Monica Rainer committed
312
313
314
315
                except:
                    self.messages.append('No masterflat found in the calibration database, it is not possible to identify the orders. The spectra will not be reduced.')
                    return

Monica Rainer's avatar
Monica Rainer committed
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
            mflat = ccdproc.CCDData.read(masterflat, unit=u.adu)
            varie.buildMaskC(mflat.data)
            self.messages.append('The extraction mask was created.')
            goodmask = ccdproc.CCDData.read(CONFIG['MASK_C'], unit=u.adu)

        gmask = goodmask.data


        roneff = hea_ima[CONFIG['RON_EFF'][0]]
        gaineff = hea_ima[CONFIG['GAIN_EFF'][0]]


        if CONFIG['USE_FLAT']['global']:
            try:
                masterflat = db.extract_dbfile(self.dbconn,'flatstr')
            except:
                masterflat = False

            if not masterflat:
                db.copy_dbfile(self.dbconn,'flatstr')
                masterflat = db.extract_dbfile(self.dbconn,'flatstr')
                self.messages.append('No masterflat found for this night, it will be taken from the calibration database: %s' % (os.path.basename(masterflat)))

            flat = ccdproc.CCDData.read(masterflat, unit=u.adu)
            froneff = flat.header[CONFIG['RON_EFF'][0]]
            fgaineff = flat.header[CONFIG['GAIN_EFF'][0]]
            meanflat = np.mean(flat.data)
            norflat = np.true_divide(flat.data,meanflat)
            with np.errstate(divide='ignore', invalid='ignore'):
                imflat = np.true_divide(imflat,norflat)

        elif CONFIG['USE_FLAT']['order']:
            try:
                masterflat = db.extract_dbfile(self.dbconn,'flatstr')
            except:
                masterflat = False

            if not masterflat:
                db.copy_dbfile(self.dbconn,'flatstr')
                masterflat = db.extract_dbfile(self.dbconn,'flatstr')
                self.messages.append('No masterflat found for this night, it will be taken from the calibration database: %s' % (os.path.basename(masterflat)))
            flat = ccdproc.CCDData.read(masterflat, unit=u.adu)
            froneff = flat.header[CONFIG['RON_EFF'][0]]
            fgaineff = flat.header[CONFIG['GAIN_EFF'][0]]


        elif CONFIG['USE_FLAT']['nor']:
            try:
                masterflat = db.extract_dbfile(self.dbconn,'flatnor')
            except:
                masterflat = False

            if not masterflat:
                db.copy_dbfile(self.dbconn,'flatnor')
                masterflat = db.extract_dbfile(self.dbconn,'flatnor')
                self.messages.append('No masterflat found for this night, it will be taken from the calibration database: %s' % (os.path.basename(masterflat)))

            flat = ccdproc.CCDData.read(masterflat, unit=u.adu)
            froneff = flat.header[CONFIG['RON_EFF'][0]]
            fgaineff = flat.header[CONFIG['GAIN_EFF'][0]]
            meanflat = 1.0
            norflat = flat.data
            with np.errstate(divide='ignore', invalid='ignore'):
                imflat = np.true_divide(imflat,norflat)


        try:
            masterlamp = db.extract_dbfile(self.dbconn,'une_str')
        except:
385
            masterlamp = False
Monica Rainer's avatar
Monica Rainer committed
386

387
        if not masterlamp:
Monica Rainer's avatar
   
Monica Rainer committed
388
389
390
391
392
393
            try:
                db.copy_dbfile(self.dbconn,'une_str')
                masterlamp = db.extract_dbfile(self.dbconn,'une_str')
                db.copy_dbfile(self.dbconn,'une_calib')
                self.messages.append('No calibration lamp found for this night, it will be taken from the calibration database: %s' % (os.path.basename(masterlamp)))
            except:
394
395
                masterlamp = False
            if not masterlamp:
Monica Rainer's avatar
   
Monica Rainer committed
396
397
                self.messages.append('No calibration lamp found in the calibration database, the spectra will not be reduced.')
                return
Monica Rainer's avatar
Monica Rainer committed
398
399
400
401
402
403
404
405
406
407
408

        mlamp = ccdproc.CCDData.read(masterlamp, unit=u.adu)

        lroneff = mlamp.header[CONFIG['RON_EFF'][0]]
        lgaineff = mlamp.header[CONFIG['GAIN_EFF'][0]]

        # read the lines to use in the wavelength calibration
        select_lines, all_lines = varie.UNe_linelist()

        # prepare the structure for the calibrated results
        heacal = OrderedDict()
Monica Rainer's avatar
Monica Rainer committed
409
410
        #stdSpectrum = np.zeros((CONFIG['N_ORD'],CONFIG['YCCD']))
        optSpectrum = np.zeros((CONFIG['N_ORD'],CONFIG['YCCD']))
Monica Rainer's avatar
Monica Rainer committed
411
        snrSpectrum = np.zeros((CONFIG['N_ORD'],CONFIG['YCCD']))
Monica Rainer's avatar
Monica Rainer committed
412
413
        fsnr = np.zeros((CONFIG['N_ORD'],CONFIG['YCCD']))
        snr = np.zeros((CONFIG['N_ORD'],CONFIG['YCCD']))
Monica Rainer's avatar
Monica Rainer committed
414

Monica Rainer's avatar
Monica Rainer committed
415
416
417
418
        x0 = np.zeros(CONFIG['N_ORD'])
        x1 = np.zeros(CONFIG['N_ORD'],dtype=int)
        x2 = np.zeros(CONFIG['N_ORD'],dtype=int)

Monica Rainer's avatar
Monica Rainer committed
419
420
421
422
        all_cosmics = 0
        for x in xrange(CONFIG['N_ORD']):
            start = x*CONFIG['W_ORD']
            end = min(start+CONFIG['W_ORD'],CONFIG['YCCD'])
Monica Rainer's avatar
Monica Rainer committed
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442

            # select only the rows wit the signal using the appropriate mask

            omask = gmask[start:end]

            order = imflat[start:end]
            if slit_pos == CONFIG['B_POS']:
                order = -order

            if CONFIG['USE_FLAT']['order']:
                ordflat = flat.data[start:end]
                # divide by masterflat normalized by its average value
                meanflat = np.mean(ordflat)
                norflat = np.true_divide(ordflat,meanflat)
                with np.errstate(divide='ignore', invalid='ignore'):
                    order = np.true_divide(order,norflat)

            ordermasked = np.ma.MaskedArray(order,mask=omask)
            goodorder = np.ma.compress_rows(ordermasked)

Monica Rainer's avatar
Monica Rainer committed
443
            optSpectrum[x],varOptFlux,profile,x0[x],x1[x],x2[x],cosmics = varie.optExtract(goodorder,gaineff,roneff,slit_pos,x)
Monica Rainer's avatar
Monica Rainer committed
444
            all_cosmics = all_cosmics + cosmics
Monica Rainer's avatar
Monica Rainer committed
445
            #print x0
Monica Rainer's avatar
Monica Rainer committed
446
447
448
449
450
451
452
453
454
455
            #snr.append(round(max( (np.mean(optSpectrum[x][1000:1050])/np.std(optSpectrum[x][1000:1050])) ,0),2))
            #optSpectrum[x] = optSpectrum[x][::-1]
            #varOptFlux = varOptFlux[::-1]

            #t3 = time.time()

            olamp = mlamp.data[start:end]
            orderlamp = np.ma.MaskedArray(olamp,mask=omask)
            goodlamp = np.ma.compress_rows(orderlamp)

Monica Rainer's avatar
Monica Rainer committed
456
457
458
            #print x1[x]
            extlamp = varie.extract(goodlamp, optSpectrum[x], x1[x], x2[x], profile, lgaineff, lroneff)
            #print 2
Monica Rainer's avatar
Monica Rainer committed
459

Monica Rainer's avatar
Monica Rainer committed
460
461
462
463
464
            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)
Monica Rainer's avatar
Monica Rainer committed
465
            snrSpectrum[x] = varie.extract(goodsnr, optSpectrum[x], x1[x], x2[x], profile, gaineff, roneff)
Monica Rainer's avatar
Monica Rainer committed
466
467
468



Monica Rainer's avatar
Monica Rainer committed
469
            if any(CONFIG['USE_FLAT'].values()) is True:
Monica Rainer's avatar
Monica Rainer committed
470
                extflat = varie.extract(norflat, optSpectrum[x], x1[x], x2[x], profile, fgaineff, froneff)
Monica Rainer's avatar
Monica Rainer committed
471
472
                #print extflat
                with np.errstate(divide='ignore', invalid='ignore'):
Monica Rainer's avatar
Monica Rainer committed
473
474
                    #fsnr[x] = (np.sqrt(froneff**2 + (fgaineff*(extflat*meanflat))))/(fgaineff*(extflat*meanflat))
                    fsnr[x] = np.true_divide(np.sqrt(froneff**2 + (fgaineff*(extflat*meanflat))),fgaineff*(extflat*meanflat))
Monica Rainer's avatar
Monica Rainer committed
475
476
477
478
                    fsnr[fsnr==np.inf] = 0
                    fsnr[fsnr==-np.inf] = 0
                    fsnr = np.nan_to_num(fsnr)
                #print fsnr
Monica Rainer's avatar
Monica Rainer committed
479
480
481
482
483
484
485


                #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)
Monica Rainer's avatar
Monica Rainer committed
486
                #snrSpectrum[x] = varie.extract(goodsnr, optSpectrum[x], x1[x], x2[x], profile, gaineff, roneff)
Monica Rainer's avatar
Monica Rainer committed
487
488
489
490
491
492
493
494
495



                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))

Monica Rainer's avatar
Monica Rainer committed
496
497
498
            else:
                fsnr[x] = np.zeros(len(optSpectrum[x]))

Monica Rainer's avatar
Monica Rainer committed
499
500
501
502
                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))
Monica Rainer's avatar
Monica Rainer committed
503
504
505
506
507
508

            calib_failed, coeffs, comments = varie.UNe_calibrate(extlamp,x+32,select_lines[x+32],all_lines[x+32])

            for comment in comments:
                self.messages.append(comment)

Monica Rainer's avatar
Monica Rainer committed
509
            keyfail = ''.join((CONFIG['CAL_FAILED'][0],str(x+32)))
Monica Rainer's avatar
Monica Rainer committed
510
511

            if calib_failed:
Monica Rainer's avatar
Monica Rainer committed
512
                heacal[keyfail] = (False,CONFIG['CAL_FAILED'][1])
Monica Rainer's avatar
Monica Rainer committed
513
514
515
516
517
518
519
520
521
                self.messages.append(' *** WARNING ***')
                self.messages.append('The default wavelength calibration for the order %s will be taken from the database and as such it will not be optimal for the night.' % (str(x+32),))
                wcalib = db.extract_dbfile(self.dbconn,'une_calib')
                wlc = ccdproc.CCDData.read(wcalib, unit=u.adu)
                for key in coeffs:
                    keyword = ''.join((CONFIG['WLCOEFFS'][key][0],str(x+32)))
                    heacal[keyword] = wlc.header[keyword]

            else:
Monica Rainer's avatar
Monica Rainer committed
522
                heacal[keyfail] = (True,CONFIG['CAL_FAILED'][1])
Monica Rainer's avatar
Monica Rainer committed
523
524
525
526
527
528
529
530
                for key in coeffs:
                    keyword = ''.join((CONFIG['WLCOEFFS'][key][0],str(x+32)))
                    heacal[keyword] = (coeffs[key],CONFIG['WLCOEFFS'][key][1])


        optSpectrum = np.asarray(optSpectrum, dtype='float32')

        self.messages.append('The nodding %s was extracted.' % str(slit),)
Monica Rainer's avatar
Monica Rainer committed
531
        if all_cosmics == (CONFIG['NCOSMIC']+1)*CONFIG['N_ORD']:
Monica Rainer's avatar
Monica Rainer committed
532
533
534
535
            self.messages.append('%s cosmics were removed (maximum iteration reached).' % str(all_cosmics),)
        else:
            self.messages.append('%s cosmics were removed.' % str(all_cosmics),)

Monica Rainer's avatar
Monica Rainer committed
536
537
538
539
540
541
542

        if slit_pos == CONFIG['A_POS']:
            keyA = ''.join((CONFIG['SPEC_USED'][0],'1','A'))
            aname = hea_ima[keyA]
            nomebase = os.path.splitext(aname)[0]
            #qui = aname.rindex('.')
            if 'grp' in fitsfile:
543
544
545
546
                msfx = ''.join(('Agrp_',CONFIG['UNMERGED'],'.fits'))
                sfx = ''.join(('Agrp_',CONFIG['MERGED'],'.fits'))
                #calname = '_'.join((nomebase,msfx))
                #calname1d = '_'.join((nomebase,sfx))
Monica Rainer's avatar
Monica Rainer committed
547
            else:
548
549
550
551
552
553
554
                msfx = ''.join(('A_',CONFIG['UNMERGED'],'.fits'))
                sfx = ''.join(('A_',CONFIG['MERGED'],'.fits'))
                #calname = '_'.join((nomebase,'A_e2ds.fits'))
                #calname1d = '_'.join((nomebase,'A_s1d.fits'))

            calname = '_'.join((nomebase,msfx))
            calname1d = '_'.join((nomebase,sfx))
Monica Rainer's avatar
Monica Rainer committed
555
556
557
558
559
560
561
562

            #calname = str(os.path.basename(fitsfile)).replace('_AB','_A')
        elif slit_pos == CONFIG['B_POS']:
            keyB = ''.join((CONFIG['SPEC_USED'][0],'1','B'))
            bname = hea_ima[keyB]
            nomebase = os.path.splitext(bname)[0]
            #qui = bname.rindex('.')
            if 'grp' in fitsfile:
563
564
565
566
                msfx = ''.join(('Bgrp_',CONFIG['UNMERGED'],'.fits'))
                sfx = ''.join(('Bgrp_',CONFIG['MERGED'],'.fits'))
                #calname = '_'.join((nomebase,msfx))
                #calname1d = '_'.join((nomebase,sfx))
Monica Rainer's avatar
Monica Rainer committed
567
            else:
568
569
570
571
572
573
574
                msfx = ''.join(('B_',CONFIG['UNMERGED'],'.fits'))
                sfx = ''.join(('B_',CONFIG['MERGED'],'.fits'))
                #calname = '_'.join((nomebase,'B_e2ds.fits'))
                #calname1d = '_'.join((nomebase,'B_s1d.fits'))

            calname = '_'.join((nomebase,msfx))
            calname1d = '_'.join((nomebase,sfx))
Monica Rainer's avatar
Monica Rainer committed
575
576
            #calname = str(os.path.basename(fitsfile)).replace('_AB','_B')
        else: 
Monica Rainer's avatar
Monica Rainer committed
577
578
579
            #print 'Wrong slit position!'
            self.messages.append('Wrong slit position!')
            return
Monica Rainer's avatar
Monica Rainer committed
580
581
582
583


        heaspe = fits.Header(hea_ima)
        heaspe[CONFIG['KEYS']['FILENAME']] = calname
584
        heaspe[CONFIG['KEYS']['IMANAME']] = calname
Monica Rainer's avatar
Monica Rainer committed
585
        drs_mjd = float(heaspe[CONFIG['KEYS']['MJD']]) + (float(heaspe[CONFIG['KEYS']['EXPTIME']])/(2.0*86400.0))
Monica Rainer's avatar
Monica Rainer committed
586
587
        heaspe[CONFIG['DRS_MJD'][0]] = (drs_mjd,CONFIG['DRS_MJD'][1])

Monica Rainer's avatar
Monica Rainer committed
588
        try:
589
590
591
            heaspe[CONFIG['MASTERFLAT'][0]] = (os.path.basename(masterflat),CONFIG['MASTERFLAT'][1])
        except:
            heaspe[CONFIG['MASTERFLAT'][0]] = ('None',CONFIG['MASTERFLAT'][1])
Monica Rainer's avatar
Monica Rainer committed
592
        try:
593
594
595
            heaspe[CONFIG['MASTERLAMP'][0]] = (os.path.basename(masterlamp),CONFIG['MASTERLAMP'][1])
        except:
            heaspe[CONFIG['MASTERLAMP'][0]] = ('None',CONFIG['MASTERLAMP'][1])
596

Monica Rainer's avatar
Monica Rainer committed
597
598
599
600
        snr[snr==np.inf] = 0
        snr[snr==-np.inf] = 0
        snr = np.nan_to_num(snr)

Monica Rainer's avatar
Monica Rainer committed
601
602
603
604
605
606
607
608
        snry = round(np.mean(snr[41]),2)
        snrj = round(np.mean(snr[29]),2)
        snrh = round(np.mean(snr[14]),2)
        snrk = round(np.mean(snr[3]),2)
        self.messages.append('Nodding %s: SNR[Y band, order=73, wl=1050 nm] = %s' % (str(slit),str(snry),))
        self.messages.append('Nodding %s: SNR[J band, order=61, wl=1250 nm] = %s' % (str(slit),str(snrj),))
        self.messages.append('Nodding %s: SNR[H band, order=46, wl=1650 nm] = %s' % (str(slit),str(snrh),))
        self.messages.append('Nodding %s: SNR[K band, order=35, wl=2200 nm] = %s' % (str(slit),str(snrk),))
Monica Rainer's avatar
Monica Rainer committed
609

Monica Rainer's avatar
Monica Rainer committed
610

Monica Rainer's avatar
Monica Rainer committed
611
        for o in xrange(CONFIG['N_ORD']):
Monica Rainer's avatar
Monica Rainer committed
612
613
614
615
616
617
618
            key_snr = ''.join((CONFIG['SNR'][0],str(o+32)))
            heaspe[key_snr] = (round(np.mean(snr[o]),2),CONFIG['SNR'][1])

        heaspe[CONFIG['WLFIT'][0]] = (CONFIG['WLFIT_FUNC'],CONFIG['WLFIT'][1])
        for hea in heacal:
            heaspe[hea] = heacal[hea]

Monica Rainer's avatar
Monica Rainer committed
619
        barycorr, hjd, bjd = varie.berv_corr(heaspe)
Monica Rainer's avatar
Monica Rainer committed
620
621
        heaspe[CONFIG['BERV'][0]] = (barycorr,CONFIG['BERV'][1])
        heaspe[CONFIG['HJD'][0]] = (hjd,CONFIG['HJD'][1])
Monica Rainer's avatar
Monica Rainer committed
622
        heaspe[CONFIG['BJD'][0]] = (bjd,CONFIG['BJD'][1])
Monica Rainer's avatar
Monica Rainer committed
623

Monica Rainer's avatar
Monica Rainer committed
624
625
        waves = np.zeros((CONFIG['N_ORD'],CONFIG['YCCD']))
        for o in xrange(CONFIG['N_ORD']):
Monica Rainer's avatar
Monica Rainer committed
626
627
            waves[o] = varie.wcalib(heaspe,o)

628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
        #spefits = fits.PrimaryHDU(optSpectrum,header=heaspe)
        #wavefits = fits.ImageHDU(waves,name='WAVE')
        #snrfits = fits.ImageHDU(snr,name='SNR')

        #results = fits.HDUList([spefits,wavefits,snrfits])

        #calname = os.path.join(CONFIG['RED_DIR'],calname)
        #results.writeto(calname,clobber=True)

        orders=np.arange(CONFIG['N_ORD'])+32
        c1 = fits.Column(name='ORDER', format='I', array=orders)
        c2 = fits.Column(name='WAVE', format=''.join((str(CONFIG['YCCD']),'D')), unit='nm', array=waves)
        c3 = fits.Column(name='FLUX', format=''.join((str(CONFIG['YCCD']),'D')), array=optSpectrum)
        c4 = fits.Column(name='SNR', format=''.join((str(CONFIG['YCCD']),'D')), array=snr)

Monica Rainer's avatar
Monica Rainer committed
643
        heaspe[CONFIG['DRS_VERSION'][0]] = (CONFIG['VERSION'], CONFIG['DRS_VERSION'][1])
Monica Rainer's avatar
Monica Rainer committed
644
645
646
647
648
649
650
651
652
        if slit_pos == CONFIG['A_POS']:
            heaspe[CONFIG['AYPOS'][0]] = (round((np.mean(x0)),2), CONFIG['AYPOS'][1])
            heaspe[CONFIG['AYPOS_LOW'][0]] = (round((np.mean(x1)),2), CONFIG['AYPOS_LOW'][1])
            heaspe[CONFIG['AYPOS_UP'][0]] = (round((np.mean(x2)),2), CONFIG['AYPOS_UP'][1])
        elif slit_pos == CONFIG['B_POS']:
            heaspe[CONFIG['BYPOS'][0]] = (round((np.mean(x0)),2), CONFIG['BYPOS'][1])
            heaspe[CONFIG['BYPOS_LOW'][0]] = (round((np.mean(x1)),2), CONFIG['BYPOS_LOW'][1])
            heaspe[CONFIG['BYPOS_UP'][0]] = (round((np.mean(x2)),2), CONFIG['BYPOS_UP'][1])

Monica Rainer's avatar
Monica Rainer committed
653

Andrea Bignamini's avatar
Andrea Bignamini committed
654
655
656
        # Add metadata to header
        heaspe = metadata.add_metadata(heaspe)

657
658
659
660
        #tbhdu = fits.BinTableHDU.from_columns([c1, c2, c3, c4],header=heaspe)
        tbhdu = fits.BinTableHDU.from_columns([c1, c2, c3, c4])
        prihdu = fits.PrimaryHDU(data=None, header=heaspe)
        hdulist = fits.HDUList([prihdu, tbhdu])
Monica Rainer's avatar
Monica Rainer committed
661
662

        calname = os.path.join(CONFIG['RED_DIR'],calname)
663
664
        #tbhdu.writeto(calname,clobber=True)
        hdulist.writeto(calname,clobber=True)
Monica Rainer's avatar
Monica Rainer committed
665

Monica Rainer's avatar
Monica Rainer committed
666
667
668
669
        #t1 = time.time()

        #print 's1d'
        #print calname1d
Monica Rainer's avatar
Monica Rainer committed
670
671
        try: obj_name = heaspe[CONFIG['KEYS']['OBJECT']]
        except: obj_name = 'NONE'
Monica Rainer's avatar
Monica Rainer committed
672
673
674

        if CONFIG['S1D']:
            #s1d = varie.create_s1d(optSpectrum,snr,heaspe)
675
            s1d, startval = varie.create_s1d(optSpectrum,heaspe)
Monica Rainer's avatar
Monica Rainer committed
676
            heaspe[CONFIG['KEYS']['FILENAME']] = calname1d
677
            heaspe[CONFIG['KEYS']['IMANAME']] = calname1d
Monica Rainer's avatar
Monica Rainer committed
678
            heaspe['CRPIX1'] = (1.,'Reference pixel')
679
            heaspe['CRVAL1'] = (startval,'Coordinates at reference pixel')
Monica Rainer's avatar
Monica Rainer committed
680
681
682
683
684
            heaspe['CDELT1'] = (CONFIG['S1D_STEP'],'Coordinates increment per pixel')

            heaspe['CTYPE1'] = ('Nanometers','Units of coordinates')
            heaspe['BUNIT'] = ('Relative Flux','Units of data values')

Andrea Bignamini's avatar
Andrea Bignamini committed
685
686
687
            # Add metadata to header
            heaspe = metadata.add_metadata(heaspe)

688
            s1dfits = fits.PrimaryHDU(s1d,header=heaspe)
Monica Rainer's avatar
Monica Rainer committed
689
690
691
            calname1d = os.path.join(CONFIG['RED_DIR'],calname1d)
            s1dfits.writeto(calname1d,clobber=True)

Monica Rainer's avatar
Monica Rainer committed
692
693
694
695
696
            rid = varie.random_id(12)
            dbreduced['s1d'] = {'slit':slit, 'path':calname1d, 'snry':snry, 'snrj':snrj, 'snrh':snrh, 'snrk':snrk, 'type':'s1d', 'name':obj_name, 'id':rid}

        rid = varie.random_id(12)
        dbreduced['ms1d'] = {'slit':slit, 'path':calname, 'snry':snry, 'snrj':snrj, 'snrh':snrh, 'snrk':snrk, 'type':'ms1d', 'name':obj_name, 'id':rid}
Monica Rainer's avatar
Monica Rainer committed
697
698
        #t2 = time.time()
        #print 's1d spectrum: %s s' %  str(t2-t1)
Monica Rainer's avatar
Monica Rainer committed
699
700

        if hea_ima[CONFIG['KEYS']['EXTMODE']] == CONFIG['EXTPAIR']:
Monica Rainer's avatar
Monica Rainer committed
701
            return calname, fsnr, straight, dbreduced, snrSpectrum
Monica Rainer's avatar
Monica Rainer committed
702
703

        elif 'grp' in fitsfile:
Monica Rainer's avatar
Monica Rainer committed
704
            return calname, fsnr, straight, dbreduced, snrSpectrum
Monica Rainer's avatar
Monica Rainer committed
705

Monica Rainer's avatar
Monica Rainer committed
706
        return calname, fsnr, False, dbreduced, snrSpectrum
Monica Rainer's avatar
Monica Rainer committed
707
708


Monica Rainer's avatar
Monica Rainer committed
709
    def combine(self,acalib,bcalib,fsnr, asnrSpec, bsnrSpec):
Monica Rainer's avatar
Monica Rainer committed
710

Monica Rainer's avatar
Monica Rainer committed
711
712
        dbreduced = {}

713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
        #abnome = acalib.replace('_A_e2ds.fits','_AB_e2ds.fits')
        #abnome1d = acalib.replace('_A_e2ds.fits','_AB_s1d.fits')
        #abnome = abnome.replace('_Agrp_e2ds.fits','_ABgrp_e2ds.fits')
        #abnome1d = abnome1d.replace('_Agrp_e2ds.fits','_ABgrp_s1d.fits')

        old = ''.join(('_A_',CONFIG['UNMERGED'],'.fits'))
        oldgrp = ''.join(('_Agrp_',CONFIG['UNMERGED'],'.fits'))
        sfx = ''.join(('_AB_',CONFIG['MERGED'],'.fits'))
        msfx = ''.join(('_AB_',CONFIG['UNMERGED'],'.fits'))
        grp = ''.join(('_ABgrp_',CONFIG['MERGED'],'.fits'))
        mgrp = ''.join(('_ABgrp_',CONFIG['UNMERGED'],'.fits'))

        abnome = acalib.replace(old,msfx)
        abnome1d = acalib.replace(old,sfx)
        abnome = abnome.replace(oldgrp,mgrp)
        abnome1d = abnome1d.replace(oldgrp,grp)

Monica Rainer's avatar
Monica Rainer committed
730
731
        #print abnome
        #print abnome1d
Monica Rainer's avatar
Monica Rainer committed
732
733

        acal = fits.open(acalib)
734
735
736

        adata = acal[1].data
        afluxes = adata.field(2)
737
738
739
740
741
742
743

        #awaves = acal[1].header
        abhea = acal[0].header
        #roneff = math.sqrt(2)*acal[1].header[CONFIG['RON_EFF'][0]]
        #gaineff = 2*acal[1].header[CONFIG['GAIN_EFF'][0]]
        roneff = math.sqrt(2)*abhea[CONFIG['RON_EFF'][0]]
        gaineff = 2*abhea[CONFIG['GAIN_EFF'][0]]
744

Monica Rainer's avatar
Monica Rainer committed
745
        bcal = fits.open(bcalib)
746
747
748
749
750

        #bfluxes = bcal[0].data
        #bwaves = bcal[0].header
        bdata = bcal[1].data
        bfluxes = bdata.field(2)
751
752
        #bwaves = bcal[1].header
        bwaves = bcal[0].header
Monica Rainer's avatar
Monica Rainer committed
753

Monica Rainer's avatar
Monica Rainer committed
754
        abcalib = np.zeros((CONFIG['N_ORD'],CONFIG['YCCD']))
Monica Rainer's avatar
Monica Rainer committed
755

Monica Rainer's avatar
Monica Rainer committed
756
        snr = np.zeros((CONFIG['N_ORD'],CONFIG['YCCD']))
Monica Rainer's avatar
Monica Rainer committed
757

Monica Rainer's avatar
Monica Rainer committed
758
        for o in xrange(CONFIG['N_ORD']):
759
760
            #bshift = varie.rebin(awaves,bfluxes[o],bwaves,o)
            bshift = varie.rebin(abhea,bfluxes[o],bwaves,o)
Monica Rainer's avatar
Monica Rainer committed
761
            abcalib[o] = (afluxes[o]+bshift)/2.0
Monica Rainer's avatar
Monica Rainer committed
762
763
764

            bsnr = varie.rebin(abhea,bsnrSpec[o],bwaves,o)
            absnr = (asnrSpec[o]+bsnr)/2.0
Monica Rainer's avatar
Monica Rainer committed
765
766
            #snr.append(max(np.mean(abcalib[o][1000:1050])/np.std(abcalib[o][1000:1050]),0))
            with np.errstate(divide='ignore', invalid='ignore'):
Monica Rainer's avatar
Monica Rainer committed
767
                #ssnr = (np.sqrt(roneff**2 + (gaineff*2*abcalib[o])))/(gaineff*2*abcalib[o])
Monica Rainer's avatar
Monica Rainer committed
768
769
770
                #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)
Monica Rainer's avatar
Monica Rainer committed
771
772
773
774
                snr[o] = 1.0/(np.sqrt(ssnr**2 + fsnr[o]**2))

        abcalib = np.asarray(abcalib, dtype='float32')

775
        #abhea = acal[0].header
776
        #abhea = acal[1].header
Monica Rainer's avatar
Monica Rainer committed
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
        abhea[CONFIG['RON_EFF'][0]] = (roneff,CONFIG['RON_EFF'][1])
        abhea[CONFIG['GAIN_EFF'][0]] = (gaineff,CONFIG['GAIN_EFF'][1])
        #abhea[CONFIG['TEXP_EFF'][0]] = (,CONFIG['TEXP_EFF'][1])
        abhea[CONFIG['KEYS']['SLIT']] = 'AB'
        n = abhea[CONFIG['KEYS']['NCOMBINE']]
        key_mjdA = ''.join((CONFIG['SPEC_MJD'][0],'1','A'))
        mjdA = float(abhea[key_mjdA]) + (float(abhea[CONFIG['KEYS']['EXPTIME']])/(2*86400))
        key_mjdB = ''.join((CONFIG['SPEC_MJD'][0],str(int(n/2)),'B'))
        mjdB = float(abhea[key_mjdB]) + (float(abhea[CONFIG['KEYS']['EXPTIME']])/(2*86400))
        mjd = (mjdA+mjdB)/2.0
        abhea[CONFIG['DRS_MJD'][0]] = (mjd,CONFIG['DRS_MJD'][1])

        snr[snr==np.inf] = 0
        snr[snr==-np.inf] = 0
        snr = np.nan_to_num(snr)

Monica Rainer's avatar
Monica Rainer committed
793
        for o in xrange(CONFIG['N_ORD']):
Monica Rainer's avatar
Monica Rainer committed
794
795
796
797
798
799
            key_snr = ''.join((CONFIG['SNR'][0],str(o+32)))
            abhea[key_snr] = (round(np.mean(snr[o]),2),CONFIG['SNR'][1])


        heaspe = fits.Header(abhea)

800
801
802
        heaspe[CONFIG['KEYS']['FILENAME']] = os.path.basename(abnome)
        heaspe[CONFIG['KEYS']['IMANAME']] = os.path.basename(abnome)

Monica Rainer's avatar
Monica Rainer committed
803
        barycorr, hjd, bjd = varie.berv_corr(heaspe)
Monica Rainer's avatar
Monica Rainer committed
804
805
        heaspe[CONFIG['BERV'][0]] = (barycorr,CONFIG['BERV'][1])
        heaspe[CONFIG['HJD'][0]] = (hjd,CONFIG['HJD'][1])
Monica Rainer's avatar
Monica Rainer committed
806
        heaspe[CONFIG['BJD'][0]] = (bjd,CONFIG['BJD'][1])
Monica Rainer's avatar
Monica Rainer committed
807
808
809


        try:
810
811
            #am_a = acal[0].header[CONFIG['AIRMASS'][0]]
            #am_b = bcal[0].header[CONFIG['AIRMASS'][0]]
812
813
814
815
            #am_a = acal[1].header[CONFIG['AIRMASS'][0]]
            #am_b = bcal[1].header[CONFIG['AIRMASS'][0]]
            am_a = abhea[CONFIG['AIRMASS'][0]]
            am_b = bwaves[CONFIG['AIRMASS'][0]]
Monica Rainer's avatar
Monica Rainer committed
816
        except:
817
818
            #am_a = acal[0].header[CONFIG['KEYS']['AM']]
            #am_b = bcal[0].header[CONFIG['KEYS']['AM']]
819
820
            am_a = abhea[CONFIG['KEYS']['AM']]
            am_b = bwaves[CONFIG['KEYS']['AM']]
Monica Rainer's avatar
Monica Rainer committed
821
822
823
824
825
826
827
828
829
830

        am = (am_a+am_b)/2.0

        heaspe[CONFIG['AIRMASS'][0]] = (am,CONFIG['AIRMASS'][1])

        try:
            heaspe[CONFIG['AIRMASS'][0]]
        except:
            heaspe[CONFIG['AIRMASS'][0]] = (heaspe[CONFIG['KEYS']['AM']],CONFIG['AIRMASS'][1])

Monica Rainer's avatar
Monica Rainer committed
831

Monica Rainer's avatar
Monica Rainer committed
832
833
        waves = np.zeros((CONFIG['N_ORD'],CONFIG['YCCD']))
        for o in xrange(CONFIG['N_ORD']):
Monica Rainer's avatar
Monica Rainer committed
834
835
            waves[o] = varie.wcalib(heaspe,o)

836
837
838
839
840
841
842
843
        #spefits = fits.PrimaryHDU(abcalib,header=heaspe)
        #wavefits = fits.ImageHDU(waves,name='WAVE')
        #snrfits = fits.ImageHDU(snr,name='SNR')

        #results = fits.HDUList([spefits,wavefits,snrfits])


        #results.writeto(abnome,clobber=True)
Monica Rainer's avatar
Monica Rainer committed
844

Monica Rainer's avatar
Monica Rainer committed
845
846
847
848
849
850
851
852
853
        snry = round(np.mean(snr[41]),2)
        snrj = round(np.mean(snr[29]),2)
        snrh = round(np.mean(snr[14]),2)
        snrk = round(np.mean(snr[3]),2)
        self.messages.append('Nodding %s: SNR[Y band, order=73, wl=1050 nm] = %s' % ('AB',str(snry)),)
        self.messages.append('Nodding %s: SNR[J band, order=61, wl=1250 nm] = %s' % ('AB',str(snrj)),)
        self.messages.append('Nodding %s: SNR[H band, order=46, wl=1650 nm] = %s' % ('AB',str(snrh)),)
        self.messages.append('Nodding %s: SNR[K band, order=35, wl=2200 nm] = %s' % ('AB',str(snrk)),)

Monica Rainer's avatar
Monica Rainer committed
854

855
856
857
858
859
860
        orders=np.arange(CONFIG['N_ORD'])+32
        c1 = fits.Column(name='ORDER', format='I', array=orders)
        c2 = fits.Column(name='WAVE', format=''.join((str(CONFIG['YCCD']),'D')), unit='nm', array=waves)
        c3 = fits.Column(name='FLUX', format=''.join((str(CONFIG['YCCD']),'D')), array=abcalib)
        c4 = fits.Column(name='SNR', format=''.join((str(CONFIG['YCCD']),'D')), array=snr)

Monica Rainer's avatar
Monica Rainer committed
861
862
        heaspe[CONFIG['DRS_VERSION'][0]] = (CONFIG['VERSION'], CONFIG['DRS_VERSION'][1])

Monica Rainer's avatar
Monica Rainer committed
863
864
865
866
        heaspe[CONFIG['BYPOS'][0]] = (bwaves[CONFIG['BYPOS'][0]], CONFIG['BYPOS'][1])
        heaspe[CONFIG['BYPOS_LOW'][0]] = (bwaves[CONFIG['BYPOS_LOW'][0]], CONFIG['BYPOS_LOW'][1])
        heaspe[CONFIG['BYPOS_UP'][0]] = (bwaves[CONFIG['BYPOS_UP'][0]], CONFIG['BYPOS_UP'][1])

Andrea Bignamini's avatar
Andrea Bignamini committed
867
868
869
        # Add metadata to header
        heaspe = metadata.add_metadata(heaspe)

870
871
872
873
        #tbhdu = fits.BinTableHDU.from_columns([c1, c2, c3, c4],header=heaspe)
        tbhdu = fits.BinTableHDU.from_columns([c1, c2, c3, c4])
        prihdu = fits.PrimaryHDU(data=None, header=heaspe)
        hdulist = fits.HDUList([prihdu, tbhdu])
874

875
876
        #tbhdu.writeto(abnome,clobber=True)
        hdulist.writeto(abnome,clobber=True)
Monica Rainer's avatar
Monica Rainer committed
877
878


Monica Rainer's avatar
Monica Rainer committed
879
880
881
        try: obj_name = heaspe[CONFIG['KEYS']['OBJECT']]
        except: obj_name = 'NONE'

Monica Rainer's avatar
Monica Rainer committed
882
        if CONFIG['S1D']:
883
884
            heaspe[CONFIG['KEYS']['FILENAME']] = os.path.basename(abnome1d)
            heaspe[CONFIG['KEYS']['IMANAME']] = os.path.basename(abnome1d)
Monica Rainer's avatar
Monica Rainer committed
885
            #s1d = varie.create_s1d(abcalib,snr,heaspe)
886
            s1d, startval = varie.create_s1d(abcalib,heaspe)
Monica Rainer's avatar
Monica Rainer committed
887
            heaspe['CRPIX1'] = (1.,'Reference pixel')
888
            heaspe['CRVAL1'] = (startval,'Coordinates at reference pixel')
Monica Rainer's avatar
Monica Rainer committed
889
890
891
892
893
            heaspe['CDELT1'] = (CONFIG['S1D_STEP'],'Coordinates increment per pixel')

            heaspe['CTYPE1'] = ('Nanometers','Units of coordinates')
            heaspe['BUNIT'] = ('Relative Flux','Units of data values')

Andrea Bignamini's avatar
Andrea Bignamini committed
894
895
896
            # Add metadata to header
            heaspe = metadata.add_metadata(heaspe)

897
            s1dfits = fits.PrimaryHDU(s1d,header=heaspe)
Monica Rainer's avatar
Monica Rainer committed
898
899
900

            s1dfits.writeto(abnome1d,clobber=True)

Monica Rainer's avatar
Monica Rainer committed
901
902
903
904
905
            rid = varie.random_id(12)
            dbreduced['s1d'] = {'slit':'AB', 'path':abnome1d, 'snry':snry, 'snrj':snrj, 'snrh':snrh, 'snrk':snrk, 'type':'s1d', 'name':obj_name, 'id':rid}

        rid = varie.random_id(12)
        dbreduced['ms1d'] = {'slit':'AB', 'path':abnome, 'snry':snry, 'snrj':snrj, 'snrh':snrh, 'snrk':snrk, 'type':'ms1d', 'name':obj_name, 'id':rid}
Monica Rainer's avatar
Monica Rainer committed
906

Monica Rainer's avatar
Monica Rainer committed
907

Monica Rainer's avatar
Monica Rainer committed
908

Monica Rainer's avatar
Monica Rainer committed
909
910
911
912
913
914
        #calibrated = np.vstack((np.concatenate(np.flipud(awaves)),np.concatenate(np.flipud(abcalib))))
        #print calibrated

        #extract = abnome.replace('.fits','.txt')
        #ascii.write(np.transpose(calibrated),extract)

Monica Rainer's avatar
Monica Rainer committed
915
        return dbreduced
Monica Rainer's avatar
Monica Rainer committed
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937


    def group_avg(self):
        noddings = []
        headers = []
        am = []
        #print self.group['noddings']
        for n in self.group['noddings']:
            nod = ccdproc.CCDData.read(n, unit=u.adu)
            noddings.append(nod)
            headers.append(nod.header)
            am.append(nod.header[CONFIG['KEYS']['AM']])

        combine_nod = ccdproc.Combiner(noddings)
        nodA = combine_nod.average_combine()
        nodA.data = np.asarray(nodA.data, dtype='float32')

        nodA.header = headers[0]

        ABnome = str(os.path.basename(self.group['noddings'][0])).replace('_AB.fits','_ABgrp.fits')

        nodA.header[CONFIG['KEYS']['FILENAME']] = ABnome
938
        nodA.header[CONFIG['KEYS']['IMANAME']] = ABnome
Monica Rainer's avatar
Monica Rainer committed
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
        nodA.header[CONFIG['GAIN_EFF'][0]] = (len(self.group['noddings'])*headers[0][CONFIG['GAIN_EFF'][0]],CONFIG['GAIN_EFF'][1])
        nodA.header[CONFIG['RON_EFF'][0]] = (headers[0][CONFIG['RON_EFF'][0]]*math.sqrt(len(self.group['noddings'])),CONFIG['RON_EFF'][1])

        nodA.header[CONFIG['KEYS']['NCOMBINE']] = len(self.group['noddings'])*2
        for n in xrange(len(self.group['noddings'])):
            keyA = ''.join((CONFIG['SPEC_USED'][0],str(n+1),'A'))
            keyB = ''.join((CONFIG['SPEC_USED'][0],str(n+1),'B'))


            mjdA = ''.join((CONFIG['SPEC_MJD'][0],str(n+1),'A'))
            mjdB = ''.join((CONFIG['SPEC_MJD'][0],str(n+1),'B'))

            readA = ''.join((CONFIG['SPEC_USED'][0],'1','A'))
            readB = ''.join((CONFIG['SPEC_USED'][0],'1','B'))

            read_mjdA = ''.join((CONFIG['SPEC_MJD'][0],'1','A'))
            read_mjdB = ''.join((CONFIG['SPEC_MJD'][0],'1','B'))

            value_keyA = headers[n][readA]
            value_keyB = headers[n][readB]

            value_mjdA = headers[n][read_mjdA]
            value_mjdB = headers[n][read_mjdB]

            nodA.header[keyA] = (value_keyA,CONFIG['SPEC_USED'][1])
            nodA.header[mjdA] = (value_mjdA,CONFIG['SPEC_MJD'][1])

            nodA.header[keyB] = (value_keyB,CONFIG['SPEC_USED'][1])
            nodA.header[mjdB] = (value_mjdB,CONFIG['SPEC_MJD'][1])

        am = np.average(np.asarray(am))
        nodA.header[CONFIG['AIRMASS'][0]] = (am,CONFIG['AIRMASS'][1])

        heaA = nodA.header
        heaB = nodA.header
        heaA[CONFIG['KEYS']['SLIT']] = CONFIG['A']
        heaA[CONFIG['KEYS']['FILENAME']] = heaA[''.join((CONFIG['SPEC_USED'][0],'1','A'))]
976
        heaA[CONFIG['KEYS']['IMANAME']] = heaA[''.join((CONFIG['SPEC_USED'][0],'1','A'))]
Monica Rainer's avatar
Monica Rainer committed
977
978
        heaB[CONFIG['KEYS']['SLIT']] = CONFIG['B']
        heaB[CONFIG['KEYS']['FILENAME']] = heaA[''.join((CONFIG['SPEC_USED'][0],'1','B'))]
979
        heaB[CONFIG['KEYS']['IMANAME']] = heaA[''.join((CONFIG['SPEC_USED'][0],'1','B'))]
Monica Rainer's avatar
Monica Rainer committed
980
981
982
983
984
985
986
987


        Atmpnome = os.path.join(CONFIG['TMP_DIR'],ABnome)
        #nodAfits = nodA.to_hdu()
        hdu = fits.PrimaryHDU(data=nodA.data,header=nodA.header)
        nodAfits = fits.HDUList([hdu])
        nodAfits.writeto(Atmpnome,clobber=True)

Monica Rainer's avatar
Monica Rainer committed
988

Monica Rainer's avatar
Monica Rainer committed
989
990
991
992
993
994
995
        for n in self.group['noddings']:
            os.remove(n)

        return Atmpnome, heaA, heaB


    def pair_process(self):
996
        #reduced = ','.join(map(os.path.basename,self.nodding))
Monica Rainer's avatar
Monica Rainer committed
997
        stamp = time.time()
998
999
        #db.insert_dbnight(self.dbnight, reduced, stamp)
        db.insert_dbnight(self.dbnight, self.nodding, stamp)
Monica Rainer's avatar
Monica Rainer committed
1000

Monica Rainer's avatar
Monica Rainer committed
1001
1002
1003
        warnings.simplefilter('ignore', category=AstropyWarning)
        if self.qualitycheck():
            #t1 = time.time()
1004
1005
1006
1007
            # Create and check nodcorr
            self.create_nodcorr()
            if not self.check_nodcorr():
                return
Monica Rainer's avatar
Monica Rainer committed
1008
1009
1010
            ab, heaA, heaB = self.createAB()
            #t2 = time.time()
            #print 'Bad pixels, create nodding: %s s' %  str(t2-t1)
Monica Rainer's avatar
   
Monica Rainer committed
1011
1012

            try:
Monica Rainer's avatar
Monica Rainer committed
1013
                acalib, fsnr, straight, dbreduced, asnrSpec = self.reduce(ab,CONFIG['A_POS'], heaA)
Monica Rainer's avatar
   
Monica Rainer committed
1014
1015
            except:
                straight = False
Monica Rainer's avatar
Monica Rainer committed
1016
1017
1018
1019
1020
1021

            try:
                db.insert_dbreduced(self.dbnight, dbreduced['s1d'], stamp)
            except:
                pass

Monica Rainer's avatar
   
Monica Rainer committed
1022
1023
1024
1025
1026
1027
            try:
                db.insert_dbreduced(self.dbnight, dbreduced['ms1d'], stamp)
            except:
                pass

            try:
Monica Rainer's avatar
Monica Rainer committed
1028
                bcalib, fsnr, straight, dbreduced, bsnrSpec = self.reduce(ab,CONFIG['B_POS'], heaB)
Monica Rainer's avatar
   
Monica Rainer committed
1029
1030
1031
            except:
                straight = False

Monica Rainer's avatar
Monica Rainer committed
1032
1033
1034
1035
            try:
                db.insert_dbreduced(self.dbnight, dbreduced['s1d'], stamp)
            except:
                pass
Monica Rainer's avatar
   
Monica Rainer committed
1036
1037
1038
1039
1040
1041
            try:
                db.insert_dbreduced(self.dbnight, dbreduced['ms1d'], stamp)
            except:
                pass

            try:
Monica Rainer's avatar
Monica Rainer committed
1042
                dbreduced = self.combine(acalib,bcalib,fsnr, asnrSpec, bsnrSpec)
Monica Rainer's avatar
   
Monica Rainer committed
1043
1044
            except:
                pass
Monica Rainer's avatar
Monica Rainer committed
1045
1046
1047
1048
1049

            try:
                db.insert_dbreduced(self.dbnight, dbreduced['s1d'], stamp)
            except:
                pass
Monica Rainer's avatar
   
Monica Rainer committed
1050
1051
1052
1053
            try:
                db.insert_dbreduced(self.dbnight, dbreduced['ms1d'], stamp)
            except:
                pass
Monica Rainer's avatar
Monica Rainer committed
1054

Monica Rainer's avatar
Monica Rainer committed
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
            if straight:
                os.remove(ab)
                os.remove(straight)
            acalib = None
            bcalib = None
        self.nodding[:] = []
        self.nodlist[:] = []
        if not self.group['noddings']:
            self.group.clear()
        return

    def group_process(self):
        warnings.simplefilter('ignore', category=AstropyWarning)
        self.pair_process()
1069
1070
1071
        try:
            if self.group['noddings'] and len(self.group['noddings'])>1 :
                ab, heaA, heaB = self.group_avg()
Monica Rainer's avatar
Monica Rainer committed
1072
1073
1074
                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)
1075
1076
1077
1078
1079
1080
1081
                acalib = None
                bcalib = None
                os.remove(ab)
                os.remove(straight)
            else:
                self.messages.append('There are no available spectra in this nodding group.')
        except:
Monica Rainer's avatar
Monica Rainer committed
1082
1083
1084
1085
1086
1087
1088
            self.messages.append('There are no available spectra in this nodding group.')
        self.group.clear()
        return

    def ingroup_process(self):
        warnings.simplefilter('ignore', category=AstropyWarning)
        self.nodding[:] = []
1089
1090
1091
        try:
            if self.group['noddings'] and len(self.group['noddings'])>1:
                ab, heaA, heaB = self.group_avg()
Monica Rainer's avatar
Monica Rainer committed
1092
1093
1094
                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)
1095
1096
1097
1098
1099
1100
1101
                acalib = None
                bcalib = None
                os.remove(ab)
                os.remove(straight)
            else:
                self.messages.append('There are no available spectra in this incomplete nodding group.')
        except:
Monica Rainer's avatar
Monica Rainer committed
1102
1103
1104
1105
1106
1107
            self.messages.append('There are no available spectra in this incomplete nodding group.')
        self.group.clear()
        return