nodding.py

"""
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:
- average the A noddings
- average the B noddings
- create average nodding A-B
- 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

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():
    def __init__(self, nodding, group, dbconn, dbnight):
        self.nodding = nodding
        self.group = group
        self.dbconn = dbconn
        self.dbnight = dbnight
        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


    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

        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

        # 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
        #nodAB.header[CONFIG['KEYS']['IMANAME']] = ABnome
        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
        """
        dbreduced = {}

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

        # straighten

        straight = fitsfile.replace('.fits','_str.fits')
        #args = [CONFIG['STRAIGHT'],fitsfile,straight,CONFIG['STRAIGHT_OPT']]
        args = [CONFIG['STRAIGHT'],fitsfile,straight]
        args.extend(CONFIG['STRAIGHT_OPT'])
        subprocess.call(args)
        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,))

        imstr = ccdproc.CCDData.read(straight, unit=u.adu)
        hea_ima = hea
        for key in CONFIG['STRAIGHT_PAR']:
            hea_ima[CONFIG['STRAIGHT_PAR'][key]] = imstr.header[CONFIG['STRAIGHT_PAR'][key]]
        imflat = imstr.data

        # 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:
                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)))
            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:
            masterlamp = False

        if not masterlamp:
            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)))

        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()
        #stdSpectrum = np.zeros((CONFIG['N_ORD'],CONFIG['YCCD']))
        optSpectrum = np.zeros((CONFIG['N_ORD'],CONFIG['YCCD']))
        fsnr = np.zeros((CONFIG['N_ORD'],CONFIG['YCCD']))
        snr = np.zeros((CONFIG['N_ORD'],CONFIG['YCCD']))

        all_cosmics = 0
        for x in xrange(CONFIG['N_ORD']):
            start = x*CONFIG['W_ORD']
            end = min(start+CONFIG['W_ORD'],CONFIG['YCCD'])

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

            optSpectrum[x],varOptFlux,profile,x1,x2,cosmics = varie.optExtract(goodorder,gaineff,roneff,slit_pos,x)
            all_cosmics = all_cosmics + cosmics

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

            extlamp = varie.extract(goodlamp, optSpectrum[x], x1, x2, profile, lgaineff, lroneff)

            if any(CONFIG['USE_FLAT'].values()) is True:
                extflat = varie.extract(norflat, optSpectrum[x], x1, x2, profile, fgaineff, froneff)
                #print extflat
                with np.errstate(divide='ignore', invalid='ignore'):
                    fsnr[x] = (np.sqrt(froneff**2 + (fgaineff*(extflat*meanflat))))/(fgaineff*(extflat*meanflat))
                    fsnr[fsnr==np.inf] = 0
                    fsnr[fsnr==-np.inf] = 0
                    fsnr = np.nan_to_num(fsnr)
                #print fsnr
            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])
                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])

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


            if calib_failed:
                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:
                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),)
        if all_cosmics == 1050:
            self.messages.append('%s cosmics were removed (maximum iteration reached).' % str(all_cosmics),)
        else:
            self.messages.append('%s cosmics were removed.' % str(all_cosmics),)


        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:
                msfx = ''.join(('Agrp_',CONFIG['UNMERGED'],'.fits'))
                sfx = ''.join(('Agrp_',CONFIG['MERGED'],'.fits'))
                #calname = '_'.join((nomebase,msfx))
                #calname1d = '_'.join((nomebase,sfx))
            else:
                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))

            #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:
                msfx = ''.join(('Bgrp_',CONFIG['UNMERGED'],'.fits'))
                sfx = ''.join(('Bgrp_',CONFIG['MERGED'],'.fits'))
                #calname = '_'.join((nomebase,msfx))
                #calname1d = '_'.join((nomebase,sfx))
            else:
                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))
            #calname = str(os.path.basename(fitsfile)).replace('_AB','_B')
        else: 
            print 'Wrong slit position!'


        heaspe = fits.Header(hea_ima)
        heaspe[CONFIG['KEYS']['FILENAME']] = calname
        drs_mjd = float(heaspe[CONFIG['KEYS']['MJD']]) + (float(heaspe[CONFIG['KEYS']['EXPTIME']])/(2*86400))
        heaspe[CONFIG['DRS_MJD'][0]] = (drs_mjd,CONFIG['DRS_MJD'][1])

        try:
            heaspe[CONFIG['MASTERFLAT'][0]] = (os.path.basename(masterflat),CONFIG['MASTERFLAT'][1])
        except:
            heaspe[CONFIG['MASTERFLAT'][0]] = ('None',CONFIG['MASTERFLAT'][1])
        try:
            heaspe[CONFIG['MASTERLAMP'][0]] = (os.path.basename(masterlamp),CONFIG['MASTERLAMP'][1])
        except:
            heaspe[CONFIG['MASTERLAMP'][0]] = ('None',CONFIG['MASTERLAMP'][1])

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

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


        for o in xrange(CONFIG['N_ORD']):
            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]

        barycorr, hjd = varie.berv_corr(heaspe)
        heaspe[CONFIG['BERV'][0]] = (barycorr,CONFIG['BERV'][1])
        heaspe[CONFIG['HJD'][0]] = (hjd,CONFIG['HJD'][1])

        waves = np.zeros((CONFIG['N_ORD'],CONFIG['YCCD']))
        for o in xrange(CONFIG['N_ORD']):
            waves[o] = varie.wcalib(heaspe,o)

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

        tbhdu = fits.BinTableHDU.from_columns([c1, c2, c3, c4],header=heaspe)

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

        #t1 = time.time()

        #print 's1d'
        #print calname1d
        try: obj_name = heaspe[CONFIG['KEYS']['OBJECT']]
        except: obj_name = 'NONE'

        if CONFIG['S1D']:
            #s1d = varie.create_s1d(optSpectrum,snr,heaspe)
            s1d, startval = varie.create_s1d(optSpectrum,heaspe)
            heaspe[CONFIG['KEYS']['FILENAME']] = calname1d
            heaspe['CRPIX1'] = (1.,'Reference pixel')
            heaspe['CRVAL1'] = (startval,'Coordinates at reference pixel')
            heaspe['CDELT1'] = (CONFIG['S1D_STEP'],'Coordinates increment per pixel')

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

            s1dfits = fits.PrimaryHDU(s1d,header=heaspe)
            calname1d = os.path.join(CONFIG['RED_DIR'],calname1d)
            s1dfits.writeto(calname1d,clobber=True)

            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}
        #t2 = time.time()
        #print 's1d spectrum: %s s' %  str(t2-t1)

        if hea_ima[CONFIG['KEYS']['EXTMODE']] == CONFIG['EXTPAIR']:
            return calname, fsnr, straight, dbreduced

        elif 'grp' in fitsfile:
            return calname, fsnr, straight, dbreduced

        return calname, fsnr, False, dbreduced


    def combine(self,acalib,bcalib,fsnr):

        dbreduced = {}

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

        #print abnome
        #print abnome1d

        acal = fits.open(acalib)

        #afluxes = acal[0].data
        #awaves = acal[0].header
        #roneff = math.sqrt(2)*acal[0].header[CONFIG['RON_EFF'][0]]
        #gaineff = 2*acal[0].header[CONFIG['GAIN_EFF'][0]]
        adata = acal[1].data
        #print adata.shape
        afluxes = adata.field(2)
        #print afluxes
        #print afluxes.shape
        awaves = acal[1].header
        roneff = math.sqrt(2)*acal[1].header[CONFIG['RON_EFF'][0]]
        gaineff = 2*acal[1].header[CONFIG['GAIN_EFF'][0]]

        bcal = fits.open(bcalib)

        #bfluxes = bcal[0].data
        #bwaves = bcal[0].header
        bdata = bcal[1].data
        bfluxes = bdata.field(2)
        bwaves = bcal[1].header

        abcalib = np.zeros((CONFIG['N_ORD'],CONFIG['YCCD']))

        snr = np.zeros((CONFIG['N_ORD'],CONFIG['YCCD']))

        for o in xrange(CONFIG['N_ORD']):
            bshift = varie.rebin(awaves,bfluxes[o],bwaves,o)
            abcalib[o] = (afluxes[o]+bshift)/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])
                snr[o] = 1.0/(np.sqrt(ssnr**2 + fsnr[o]**2))

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

        #abhea = acal[0].header
        abhea = acal[1].header
        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)

        for o in xrange(CONFIG['N_ORD']):
            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)

        barycorr, hjd = varie.berv_corr(heaspe)
        heaspe[CONFIG['BERV'][0]] = (barycorr,CONFIG['BERV'][1])
        heaspe[CONFIG['HJD'][0]] = (hjd,CONFIG['HJD'][1])


        try:
            #am_a = acal[0].header[CONFIG['AIRMASS'][0]]
            #am_b = bcal[0].header[CONFIG['AIRMASS'][0]]
            am_a = acal[1].header[CONFIG['AIRMASS'][0]]
            am_b = bcal[1].header[CONFIG['AIRMASS'][0]]
        except:
            #am_a = acal[0].header[CONFIG['KEYS']['AM']]
            #am_b = bcal[0].header[CONFIG['KEYS']['AM']]
            am_a = acal[1].header[CONFIG['KEYS']['AM']]
            am_b = bcal[1].header[CONFIG['KEYS']['AM']]

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


        waves = np.zeros((CONFIG['N_ORD'],CONFIG['YCCD']))
        for o in xrange(CONFIG['N_ORD']):
            waves[o] = varie.wcalib(heaspe,o)

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

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


        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)

        tbhdu = fits.BinTableHDU.from_columns([c1, c2, c3, c4],header=heaspe)

        tbhdu.writeto(abnome,clobber=True)


        try: obj_name = heaspe[CONFIG['KEYS']['OBJECT']]
        except: obj_name = 'NONE'

        if CONFIG['S1D']:

            #s1d = varie.create_s1d(abcalib,snr,heaspe)
            s1d, startval = varie.create_s1d(abcalib,heaspe)
            heaspe['CRPIX1'] = (1.,'Reference pixel')
            heaspe['CRVAL1'] = (startval,'Coordinates at reference pixel')
            heaspe['CDELT1'] = (CONFIG['S1D_STEP'],'Coordinates increment per pixel')

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

            s1dfits = fits.PrimaryHDU(s1d,header=heaspe)

            s1dfits.writeto(abnome1d,clobber=True)

            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}



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

        return dbreduced


    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
        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'))]
        heaB[CONFIG['KEYS']['SLIT']] = CONFIG['B']
        heaB[CONFIG['KEYS']['FILENAME']] = heaA[''.join((CONFIG['SPEC_USED'][0],'1','B'))]


        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)


        for n in self.group['noddings']:
            os.remove(n)

        return Atmpnome, heaA, heaB


    def pair_process(self):
        #reduced = ','.join(map(os.path.basename,self.nodding))
        stamp = time.time()
        #db.insert_dbnight(self.dbnight, reduced, stamp)
        db.insert_dbnight(self.dbnight, self.nodding, stamp)

        warnings.simplefilter('ignore', category=AstropyWarning)
        if self.qualitycheck():
            #t1 = time.time()
            ab, heaA, heaB = self.createAB()
            #t2 = time.time()
            #print 'Bad pixels, create nodding: %s s' %  str(t2-t1)
            acalib, fsnr, straight, dbreduced = self.reduce(ab,CONFIG['A_POS'], heaA)

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

            bcalib, fsnr, straight, dbreduced = self.reduce(ab,CONFIG['B_POS'], heaB)
            try:
                db.insert_dbreduced(self.dbnight, dbreduced['s1d'], stamp)
            except:
                pass
            db.insert_dbreduced(self.dbnight, dbreduced['ms1d'], stamp)

            dbreduced = self.combine(acalib,bcalib,fsnr)
            try:
                db.insert_dbreduced(self.dbnight, dbreduced['s1d'], stamp)
            except:
                pass
            db.insert_dbreduced(self.dbnight, dbreduced['ms1d'], stamp)

            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()
        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 = None
                bcalib = None
                os.remove(ab)
                os.remove(straight)
            else:
                self.messages.append('There are no available spectra in this nodding group.')
        except:
            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[:] = []
        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 = None
                bcalib = None
                os.remove(ab)
                os.remove(straight)
            else:
                self.messages.append('There are no available spectra in this incomplete nodding group.')
        except:
            self.messages.append('There are no available spectra in this incomplete nodding group.')
        self.group.clear()
        return