Source code for pypeit.spectrographs.gtc_osiris

"""
Module for GTC OSIRIS specific methods.

.. include:: ../include/links.rst
"""
import numpy as np

from pypeit import msgs
from pypeit import telescopes
from pypeit.core import parse
from pypeit.core import framematch
from pypeit.spectrographs import spectrograph
from pypeit.images import detector_container
from astropy import wcs, units
import astropy.io.fits as fits
from astropy.time import Time
from astropy.coordinates import SkyCoord, EarthLocation
from IPython import embed




class GTCOSIRISPlusSpectrograph(spectrograph.Spectrograph):
    """
    Child of Spectrograph to handle GTC/OSIRIS specific code
    """
    ndet = 1
    name = 'gtc_osiris_plus'
    telescope = telescopes.GTCTelescopePar()
    camera = 'OSIRIS'
    url = 'http://www.gtc.iac.es/instruments/osiris/'
    header_name = 'OSIRIS'
    supported = True
    comment = 'See :doc:`gtc_osiris`'

    def __init__(self):
        super().__init__()
        self.location = EarthLocation.of_site('lapalma')



    def get_detector_par(self, det, hdu=None):
        """
        Return metadata for the selected detector.

        Detector data from `here
        <http://www.gtc.iac.es/instruments/osiris/>`__.

        Args:
            det (:obj:`int`):
                1-indexed detector number.
            hdu (`astropy.io.fits.HDUList`_, optional):
                The open fits file with the raw image of interest.  If not
                provided, frame-dependent parameters are set to a default.

        Returns:
            :class:`~pypeit.images.detector_container.DetectorContainer`:
            Object with the detector metadata.
        """
        binning = '1,1' if hdu is None else self.get_meta_value(self.get_headarr(hdu), 'binning')
        gain = 1.90 if hdu is None else self.get_headarr(hdu)[0]['GAIN']
        ronoise = 4.3 if hdu is None else self.get_headarr(hdu)[0]['RDNOISE']

        # Detector 1
        detector_dict1 = dict(
            binning         = binning,
            det             = 1,
            dataext         = 0,
            specaxis        = 1,
            specflip        = True,
            spatflip        = False,
            platescale      = 0.125,  # arcsec per pixel
            darkcurr        = 5.0,  #e-/hr/pixel
            saturation      = 65535., # ADU
            nonlinear       = 0.95,
            mincounts       = 0,
            numamplifiers   = 1,
            gain            = np.atleast_1d([gain]),
            ronoise         = np.atleast_1d([ronoise]),
            datasec         = np.atleast_1d('[180:4112,50:4096]'),
            oscansec        = np.atleast_1d('[180:4112,8:46]')  # Trim down the oscansec - looks like some bad pixels
            )

        detectors = [detector_dict1]
        # Return
        return detector_container.DetectorContainer(**detectors[det-1])





    @classmethod
    def default_pypeit_par(cls):
        """
        Return the default parameters to use for this instrument.

        Returns:
            :class:`~pypeit.par.pypeitpar.PypeItPar`: Parameters required by
            all of PypeIt methods.
        """
        par = super().default_pypeit_par()

        # Ignore PCA
        par['calibrations']['slitedges']['sync_predict'] = 'nearest'
        par['calibrations']['slitedges']['bound_detector'] = True

        # Set pixel flat combination method
        par['calibrations']['pixelflatframe']['process']['combine'] = 'median'
        # Wavelength calibration methods
        par['calibrations']['wavelengths']['method'] = 'full_template'
        par['calibrations']['wavelengths']['lamps'] = ['XeI','HgI','NeI','ArI']
        par['calibrations']['wavelengths']['reid_cont_sub'] = False

        # Set the default exposure time ranges for the frame typing
        par['scienceframe']['exprng'] = [90, None]
        par['calibrations']['biasframe']['exprng'] = [None, 0.001]
        par['calibrations']['darkframe']['exprng'] = [999999, None]     # No dark frames
        par['calibrations']['pinholeframe']['exprng'] = [999999, None]  # No pinhole frames
        par['calibrations']['arcframe']['exprng'] = [None, None]  # Long arc exposures
        par['calibrations']['arcframe']['process']['clip'] = False
        par['calibrations']['standardframe']['exprng'] = [None, 300]
        # Multiple arcs with different lamps, so can't median combine nor clip, also need to remove continuum
        par['calibrations']['arcframe']['process']['combine'] = 'mean'
        par['calibrations']['arcframe']['process']['subtract_continuum'] = True
        par['calibrations']['tiltframe']['process']['clip'] = False
        par['calibrations']['tiltframe']['process']['combine'] = 'mean'
        par['calibrations']['tiltframe']['process']['subtract_continuum'] = True

        # Increase the wave tilts order, given the longish slit
        par['calibrations']['tilts']['spat_order'] = 5
        par['calibrations']['tilts']['spec_order'] = 5

        # Only extract one object per standard frame
        par['reduce']['findobj']['maxnumber_std'] = 1

        # Turn off the 2D fit - this seems to be giving bad results for OSIRIS
        par['reduce']['skysub']['no_poly'] = True

        # Don't extrapolate the sensitivity function for the low resolution gratings
        # Sensitivity function parameters
        par['sensfunc']['extrap_blu'] = 0.0
        par['sensfunc']['extrap_red'] = 0.0
        par['fluxcalib']['extrap_sens'] = True
        par['sensfunc']['algorithm'] = 'IR'
        par['sensfunc']['polyorder'] = 13
        par['sensfunc']['IR']['maxiter'] = 2
        par['sensfunc']['IR']['telgridfile'] = 'TellPCA_3000_26000_R10000.fits'
        return par




    def init_meta(self):
        """
        Define how metadata are derived from the spectrograph files.

        That is, this associates the PypeIt-specific metadata keywords
        with the instrument-specific header cards using :attr:`meta`.
        """
        self.meta = {}
        # Required (core)
        self.meta['ra'] = dict(ext=0, card='RA')
        self.meta['dec'] = dict(ext=0, card='DEC')
        self.meta['target'] = dict(ext=0, card='object')
        self.meta['idname'] = dict(ext=0, card='obsmode')
        self.meta['decker'] = dict(ext=0, card='MASKNAME')
        self.meta['binning'] = dict(card=None, compound=True)
        self.meta['detector'] = dict(ext=0, card='detector')
        self.meta['mjd'] = dict(ext=0, card='MJD-OBS')
        self.meta['exptime'] = dict(ext=0, card='EXPTIME')
        self.meta['airmass'] = dict(ext=0, card='AIRMASS')
        self.meta['dispname'] = dict(ext=0, card='GRISM')
        self.meta['datasec'] = dict(ext=0, card='DETSIZE')
        self.meta['dichroic'] = dict(ext=0, card='FILTER1')
        self.meta['instrument'] = dict(ext=0, card='INSTRUME')
        self.meta['slitwid'] = dict(card=None, compound=True)




    def compound_meta(self, headarr, meta_key):
        """
        Methods to generate metadata requiring interpretation of the header
        data, instead of simply reading the value of a header card.

        Args:
            headarr (:obj:`list`):
                List of `astropy.io.fits.Header`_ objects.
            meta_key (:obj:`str`):
                Metadata keyword to construct.

        Returns:
            object: Metadata value read from the header(s).
        """
        if meta_key == 'binning':
            binspatial, binspec = parse.parse_binning(headarr[0]['CCDSUM'])
            binning = parse.binning2string(binspec, binspatial)[::-1]
            return binning
        elif meta_key == 'pressure':
            try:
                return headarr[0]['PRESSURE']  # Must be in astropy.units.mbar
            except KeyError:
                msgs.warn("Pressure is not in header")
                msgs.info("The default pressure will be assumed: 611 mbar")
                return 611.0
        elif meta_key == 'temperature':
            try:
                return headarr[0]['TAMBIENT']  # Must be in astropy.units.deg_C
            except KeyError:
                msgs.warn("Temperature is not in header")
                msgs.info("The default temperature will be assumed: 1.5 deg C")
                return 1.5
        elif meta_key == 'humidity':
            try:
                return headarr[0]['HUMIDITY']
            except KeyError:
                msgs.warn("Humidity is not in header")
                msgs.info("The default relative humidity will be assumed: 20 %")
                return 20.0
        elif meta_key == 'parangle':
            try:
                msgs.work("Parallactic angle is not available for MAAT - DAR correction may be incorrect")
                return headarr[0]['PARANG']  # Must be expressed in radians
            except KeyError:
                msgs.error("Parallactic angle is not in header")
        elif meta_key == 'obstime':
            return Time(headarr[0]['DATE-END'])
        elif meta_key == 'gain':
            return headarr[0]['GAIN']
        elif meta_key == 'slitwid':
            if self.name == "gtc_maat":
                msgs.warn("HACK FOR MAAT SIMS --- NEED TO GET SLICER SCALE FROM HEADER, IDEALLY")
                return 0.305 / 3600.0
            elif self.name == "gtc_osiris_plus":
                return headarr[0]['SLITW']/3600.0   # Convert slit width from arcseconds to degrees
            else:
                msgs.error("Could not determine slit width from header information")




    def configuration_keys(self):
        """
        Return the metadata keys that define a unique instrument
        configuration.

        This list is used by :class:`~pypeit.metadata.PypeItMetaData` to
        identify the unique configurations among the list of frames read
        for a given reduction.

        Returns:
            :obj:`list`: List of keywords of data pulled from file headers
            and used to constuct the :class:`~pypeit.metadata.PypeItMetaData`
            object.
        """
        return ['dispname', 'decker', 'binning']




    def raw_header_cards(self):
        """
        Return additional raw header cards to be propagated in
        downstream output files for configuration identification.

        The list of raw data FITS keywords should be those used to populate
        the :meth:`~pypeit.spectrographs.spectrograph.Spectrograph.configuration_keys`
        or are used in :meth:`~pypeit.spectrographs.spectrograph.Spectrograph.config_specific_par`
        for a particular spectrograph, if different from the name of the
        PypeIt metadata keyword.

        This list is used by :meth:`~pypeit.spectrographs.spectrograph.Spectrograph.subheader_for_spec`
        to include additional FITS keywords in downstream output files.

        Returns:
            :obj:`list`: List of keywords from the raw data files that should
            be propagated in output files.
        """
        return ['GRISM', 'MASKNAME', 'CCDSUM', 'obsmode']




    def check_frame_type(self, ftype, fitstbl, exprng=None):
        """
        Check for frames of the provided type.

        Args:
            ftype (:obj:`str`):
                Type of frame to check. Must be a valid frame type; see
                frame-type :ref:`frame_type_defs`.
            fitstbl (`astropy.table.Table`_):
                The table with the metadata for one or more frames to check.
            exprng (:obj:`list`, optional):
                Range in the allowed exposure time for a frame of type
                ``ftype``. See
                :func:`pypeit.core.framematch.check_frame_exptime`.

        Returns:
            `numpy.ndarray`_: Boolean array with the flags selecting the
            exposures in ``fitstbl`` that are ``ftype`` type frames.
        """
        good_exp = framematch.check_frame_exptime(fitstbl['exptime'], exprng)
        if ftype in ['science', 'standard']:
            return good_exp & (np.logical_not(np.char.startswith(np.char.lower(fitstbl['target']), 'arclamp'))) & \
                   (np.char.lower(fitstbl['target']) != 'spectralflat') & \
                   (np.char.lower(fitstbl['target']) != 'bias')
        if ftype in ['arc', 'tilt']:
            return good_exp & (np.char.startswith(np.char.lower(fitstbl['target']), 'arclamp'))
        if ftype in ['pixelflat', 'trace', 'illumflat']:
            return good_exp & (np.char.lower(fitstbl['target']) == 'spectralflat')
        if ftype == 'bias':
            return good_exp & (np.char.lower(fitstbl['target']) == 'bias')

        msgs.warn('Cannot determine if frames are of type {0}.'.format(ftype))
        return np.zeros(len(fitstbl), dtype=bool)




    def config_independent_frames(self):
        """
        Define frame types that are independent of the fully defined
        instrument configuration.

        Bias and dark frames are considered independent of a configuration.
        Standards are assigned to the correct configuration frame group by
        grism (i.e. ignoring that they are taken with a wider slit).
        See :func:`~pypeit.metadata.PypeItMetaData.set_configurations`.

        Returns:
            :obj:`dict`: Dictionary where the keys are the frame types that
            are configuration independent and the values are the metadata
            keywords that can be used to assign the frames to a configuration
            group.
        """
        return {'standard': 'dispname', 'bias': 'binning', 'dark': 'binning'}




    def config_specific_par(self, scifile, inp_par=None):
        """
        Modify the PypeIt parameters to hard-wired values used for
        specific instrument configurations.

        Args:
            scifile (:obj:`str`):
                File to use when determining the configuration and how
                to adjust the input parameters.
            inp_par (:class:`~pypeit.par.parset.ParSet`, optional):
                Parameter set used for the full run of PypeIt.  If None,
                use :func:`default_pypeit_par`.

        Returns:
            :class:`~pypeit.par.parset.ParSet`: The PypeIt parameter set
            adjusted for configuration specific parameter values.
        """
        # Start with instrument wide
        par = super().config_specific_par(scifile, inp_par=inp_par)

        if self.get_meta_value(scifile, 'idname') == 'OsirisMOS':
            par['reduce']['findobj']['find_trim_edge'] = [1,1]
            par['calibrations']['slitedges']['sync_predict'] = 'pca'
            par['calibrations']['slitedges']['det_buffer'] = 1
        elif self.get_meta_value(scifile, 'idname') == 'OsirisLongSlitSpectroscopy':
            # Do not tweak the slit edges for longslit
            par['calibrations']['flatfield']['tweak_slits'] = False

        # Wavelength calibration and setup-dependent parameters
        if self.get_meta_value(scifile, 'dispname') == 'R300B':
            par['calibrations']['wavelengths']['lamps'] = ['XeI','HgI','NeI']
            par['calibrations']['wavelengths']['reid_arxiv'] = 'gtc_osiris_R300B.fits'
            par['reduce']['findobj']['find_min_max'] = [750, 2051]
            par['calibrations']['slitedges']['det_min_spec_length'] = 0.25
            par['calibrations']['slitedges']['fit_min_spec_length'] = 0.25
            par['calibrations']['slitedges']['smash_range'] = [0.38, 0.62]
            par['calibrations']['flatfield']['slit_illum_finecorr'] = False
            par['reduce']['cube']['wave_min'] = 3600.0
            par['reduce']['cube']['wave_max'] = 7200.0
        elif self.get_meta_value(scifile, 'dispname') == 'R300R':
            par['calibrations']['wavelengths']['lamps'] = ['XeI','HgI','NeI']
            par['calibrations']['wavelengths']['reid_arxiv'] = 'gtc_osiris_R300R.fits'
            par['reduce']['findobj']['find_min_max'] = [750, 2051]
            par['calibrations']['slitedges']['det_min_spec_length'] = 0.25
            par['calibrations']['slitedges']['fit_min_spec_length'] = 0.25
            par['calibrations']['slitedges']['smash_range'] = [0.38, 0.62]
            par['calibrations']['flatfield']['slit_illum_finecorr'] = False
            par['reduce']['cube']['wave_min'] = 4800.0
            par['reduce']['cube']['wave_max'] = 10000.0
        elif self.get_meta_value(scifile, 'dispname') == 'R500B':
            par['calibrations']['wavelengths']['lamps'] = ['HgI','NeI']
            par['calibrations']['wavelengths']['reid_arxiv'] = 'gtc_osiris_R500B.fits'
            par['reduce']['findobj']['find_min_max'] = [500, 2051]
            par['reduce']['cube']['wave_min'] = 3600.0
            par['reduce']['cube']['wave_max'] = 7200.0
        elif self.get_meta_value(scifile, 'dispname') == 'R500R':
            par['calibrations']['wavelengths']['lamps'] = ['XeI','HgI','NeI']
            par['calibrations']['wavelengths']['reid_arxiv'] = 'gtc_osiris_R500R.fits'
            par['reduce']['findobj']['find_min_max'] = [450, 2051]
            par['reduce']['cube']['wave_min'] = 4800.0
            par['reduce']['cube']['wave_max'] = 10000.0
        elif self.get_meta_value(scifile, 'dispname') == 'R1000B':
            par['calibrations']['wavelengths']['lamps'] = ['ArI','HgI','NeI']
            par['calibrations']['wavelengths']['reid_arxiv'] = 'gtc_osiris_R1000B.fits'
        elif self.get_meta_value(scifile, 'dispname') == 'R1000R':
            par['calibrations']['wavelengths']['lamps'] = ['XeI','HgI','NeI']
            par['calibrations']['wavelengths']['reid_arxiv'] = 'gtc_osiris_R1000R.fits'
        elif self.get_meta_value(scifile, 'dispname') == 'R2000B':
            par['calibrations']['wavelengths']['fwhm'] = 15.0
            par['calibrations']['wavelengths']['lamps'] = ['XeI','HgI','NeI']
            par['calibrations']['wavelengths']['reid_arxiv'] = 'gtc_osiris_R2000B.fits'
        elif self.get_meta_value(scifile, 'dispname') == 'R2500U':
            par['calibrations']['wavelengths']['lamps'] = ['XeI','HgI']
            par['calibrations']['wavelengths']['reid_arxiv'] = 'gtc_osiris_R2500U.fits'
            par['calibrations']['wavelengths']['nsippet'] = 1
        elif self.get_meta_value(scifile, 'dispname') == 'R2500V':
            par['calibrations']['wavelengths']['lamps'] = ['HgI','NeI','XeI']
            par['calibrations']['wavelengths']['reid_arxiv'] = 'gtc_osiris_R2500V.fits'
        elif self.get_meta_value(scifile, 'dispname') == 'R2500R':
            par['calibrations']['wavelengths']['lamps'] = ['ArI','HgI','NeI']
            par['calibrations']['wavelengths']['reid_arxiv'] = 'gtc_osiris_R2500R.fits'
        elif self.get_meta_value(scifile, 'dispname') == 'R2500I':
            par['calibrations']['wavelengths']['lamps'] = ['ArI','XeI','NeI']
            par['calibrations']['wavelengths']['reid_arxiv'] = 'gtc_osiris_R2500I.fits'
            par['sensfunc']['algorithm'] = 'IR'
            par['sensfunc']['IR']['telgridfile'] = "TellPCA_3000_26000_R10000.fits"
        else:
            msgs.warn('gtc_osiris.py: template arc missing for this grism! Trying holy-grail...')
            par['calibrations']['wavelengths']['method'] = 'holy-grail'

        # Return
        return par




    def bpm(self, filename, det, shape=None, msbias=None):
        """
        Generate a default bad-pixel mask.

        Even though they are both optional, either the precise shape for
        the image (``shape``) or an example file that can be read to get
        the shape (``filename`` using :func:`get_image_shape`) *must* be
        provided.

        Args:
            filename (:obj:`str` or None):
                An example file to use to get the image shape.
            det (:obj:`int`):
                1-indexed detector number to use when getting the image
                shape from the example file.
            shape (tuple, optional):
                Processed image shape
                Required if filename is None
                Ignored if filename is not None
            msbias (`numpy.ndarray`_, optional):
                Processed bias frame used to identify bad pixels. **This is
                ignored for KCWI.**

        Returns:
            `numpy.ndarray`_: An integer array with a masked value set
            to 1 and an unmasked value set to 0.  All values are set to
            0.
        """
        # Call the base-class method to generate the empty bpm; msbias is always set to None.
        bpm_img = super().bpm(filename, det, shape=shape, msbias=None)

        # Extract some header info
        head0 = fits.getheader(filename, ext=0)
        binning = self.get_meta_value([head0], 'binning')

        msgs.warn("Bad pixel mask is not available for det={0:d} binning={1:s}".format(det, binning))
        # Construct a list of the bad columns
        bc = []
        # TODO :: Add BPM

        # Apply these bad columns to the mask
        for bb in range(len(bc)):
            bpm_img[bc[bb][2]:bc[bb][3] + 1, bc[bb][0]:bc[bb][1] + 1] = 1

        return bpm_img




    def tweak_standard(self, wave_in, counts_in, counts_ivar_in, gpm_in, meta_table, trim_std_pixs=None,
                       log10_blaze_function=None, debug=False):
        """

        This routine is for performing instrument/disperser specific tweaks to standard stars so that sensitivity
        function fits will be well behaved. For example, masking second order light. For instruments that don't
        require such tweaks it will just return the inputs, but for instruments that do this function is overloaded
        with a method that performs the tweaks.

        Parameters
        ----------
        wave_in: `numpy.ndarray`_
            Input standard star wavelengths (:obj:`float`, ``shape = (nspec,)``)
        counts_in: `numpy.ndarray`_
            Input standard star counts (:obj:`float`, ``shape = (nspec,)``)
        counts_ivar_in: `numpy.ndarray`_
            Input inverse variance of standard star counts (:obj:`float`, ``shape = (nspec,)``)
        gpm_in: `numpy.ndarray`_
            Input good pixel mask for standard (:obj:`bool`, ``shape = (nspec,)``)
        meta_table: :obj:`dict`
            Table containing meta data that is slupred from the :class:`~pypeit.specobjs.SpecObjs`
            object.  See :meth:`~pypeit.specobjs.SpecObjs.unpack_object` for the
            contents of this table.
        trim_std_pixs: :obj:`list` or :obj:`tuple`, optional
            List or tuple of two integers specifying the number of pixels to
            trim from the start and end of the standard star spectrum. If None,
            no trimming is applied. Default=None.
        log10_blaze_function: `numpy.ndarray`_ or None
            Input blaze function to be tweaked, optional. Default=None.


        Returns
        -------
        wave_out: `numpy.ndarray`_
            Output standard star wavelengths (:obj:`float`, ``shape = (nspec,)``)
        counts_out: `numpy.ndarray`_
            Output standard star counts (:obj:`float`, ``shape = (nspec,)``)
        counts_ivar_out: `numpy.ndarray`_
            Output inverse variance of standard star counts (:obj:`float`, ``shape = (nspec,)``)
        gpm_out: `numpy.ndarray`_
            Output good pixel mask for standard (:obj:`bool`, ``shape = (nspec,)``)
        log10_blaze_function_out: `numpy.ndarray`_ or None
            Output blaze function after being tweaked.
        """
        # If trim_std_pixs is provided, use the base class method to trim the standard star
        if trim_std_pixs is not None:
            return super().tweak_standard(wave_in, counts_in, counts_ivar_in, gpm_in, meta_table,
                                          trim_std_pixs=trim_std_pixs, log10_blaze_function=log10_blaze_function)

        # Could check the wavelenghts here to do something more robust to header/meta data issues
        if 'R300R' in meta_table['DISPNAME']:
            wave_blue = 4800.0  # blue wavelength below which there is contamination
            wave_red = 9300.0  # red wavelength above which the spectrum is contaminated
        elif 'R500R' in meta_table['DISPNAME']:
            wave_blue = 4800.0  # blue wavelength below which there is contamination
            wave_red = 9300.0  # red wavelength above which the spectrum is contaminated
        elif 'R300B' in meta_table['DISPNAME']:
            wave_blue = 3400.0  # blue wavelength below which there is contamination
            wave_red = 6500.0  # red wavelength above which the spectrum is contaminated
        elif 'R500B' in meta_table['DISPNAME']:
            wave_blue = 3400.0  # blue wavelength below which there is contamination
            wave_red = 6500.0  # red wavelength above which the spectrum is contaminated
        elif 'R1000B' in meta_table['DISPNAME']:
            wave_blue = 3500.0  # blue wavelength below which there is contamination
            wave_red = 6300.0  # red wavelength above which the spectrum is contaminated
        else:
            # keep everything the same
            wave_blue = 0.0
            wave_red = np.inf

        second_order_region = (wave_in < wave_blue) | (wave_in > wave_red)
        wave = wave_in.copy()
        counts = counts_in.copy()
        counts_ivar = counts_ivar_in.copy()
        gpm = gpm_in.copy()

        wave[second_order_region] = 0.0
        counts[second_order_region] = 0.0
        counts_ivar[second_order_region] = 0.0
        # By setting the wavelengths to zero, we guarantee that the sensitvity function will only be computed
        # over the valid wavelength region. While we could mask, this would still produce a wave_min and wave_max
        # for the zeropoint that includes the bad regions, and the polynomial fits will extrapolate crazily there
        gpm[second_order_region] = False

        if log10_blaze_function is not None:
            log10_blaze_function_out = log10_blaze_function.copy()
            log10_blaze_function_out[second_order_region] = 0.0
        else:
            log10_blaze_function_out = None

        if debug:
           from matplotlib import pyplot as plt
           from pypeit import utils
           counts_sigma = np.sqrt(utils.inverse(counts_ivar_in))
           plt.plot(wave_in, counts, color='red', alpha=0.7, label='apodized flux')
           plt.plot(wave_in, counts_in, color='black', alpha=0.7, label='flux')
           plt.plot(wave_in, counts_sigma, color='blue', alpha=0.7, label='flux')
           plt.axvline(wave_blue, color='blue')
           plt.axvline(wave_red, color='red')
           plt.legend()
           plt.show()

        return wave, counts, counts_ivar, gpm, log10_blaze_function_out






class GTCMAATSpectrograph(GTCOSIRISPlusSpectrograph):
    pypeline = 'SlicerIFU'
    name = 'gtc_maat'



    def init_meta(self):
        super().init_meta()
        self.meta['obstime'] = dict(card=None, compound=True, required=False)
        self.meta['pressure'] = dict(card=None, compound=True, required=False)
        self.meta['temperature'] = dict(card=None, compound=True, required=False)
        self.meta['humidity'] = dict(card=None, compound=True, required=False)
        self.meta['parangle'] = dict(card=None, compound=True, required=False)




    @classmethod
    def default_pypeit_par(cls):
        par = super().default_pypeit_par()

        # LACosmics parameters
        par['scienceframe']['process']['sigclip'] = 4.0
        par['scienceframe']['process']['objlim'] = 1.5
        par['scienceframe']['process']['use_illumflat'] = False  # illumflat is applied when building the relative scale image in reduce.py, so should be applied to scienceframe too.
        par['scienceframe']['process']['use_specillum'] = False  # apply relative spectral illumination
        par['scienceframe']['process']['spat_flexure_correct'] = False  # don't correct for spatial flexure - varying spatial illumination profile could throw this correction off. Also, there's no way to do astrometric correction if we can't correct for spatial flexure of the contbars frames
        par['scienceframe']['process']['use_biasimage'] = False
        par['scienceframe']['process']['use_darkimage'] = False
        par['calibrations']['flatfield']['slit_illum_finecorr'] = False
        # Don't do 1D extraction for 3D data - it's meaningless because the DAR correction must be performed on the 3D data.
        par['reduce']['extraction']['skip_extraction'] = True  # Because extraction occurs before the DAR correction, don't extract

        # Decrease the wave tilts order, given the shorter slits of the IFU
        par['calibrations']['tilts']['spat_order'] = 1
        par['calibrations']['tilts']['spec_order'] = 1

        # Tweak the slit edges using the gradient method for SlicerIFU
        par['calibrations']['slitedges']['pad'] = 0  # Do not pad the slits - this ensures that the tweak_edges method=gradient guarantees that the edges are defined at the maximum gradient.
        par['calibrations']['flatfield']['tweak_slits'] = True  # Tweak the slit edges
        par['calibrations']['flatfield']['tweak_method'] = 'gradient'  # The gradient method is better for SlicerIFU.
        par['calibrations']['flatfield']['tweak_slits_thresh'] = 0.0  # Make sure the full slit is used (i.e. when the illumination fraction is > 0.5)
        par['calibrations']['flatfield']['tweak_slits_maxfrac'] = 0.0  # Make sure the full slit is used (i.e. no padding)

        # Make sure that this is reduced as a slit (as opposed to fiber) spectrograph
        par['reduce']['cube']['slit_spec'] = True
        par['reduce']['cube']['combine'] = False  # Make separate spec3d files from the input spec2d files

        # Sky subtraction parameters
        par['reduce']['skysub']['no_poly'] = True
        par['reduce']['skysub']['bspline_spacing'] = 0.6
        par['reduce']['skysub']['joint_fit'] = False
        par['reduce']['findobj']['skip_skysub'] = True
        par['reduce']['findobj']['skip_final_global'] = True

        # Don't correct flexure by default, but you should use slitcen,
        # because this is a slit-based IFU where no objects are extracted.
        par['flexure']['spec_method'] = 'skip'
        par['flexure']['spec_maxshift'] = 3  # Just in case someone switches on spectral flexure, this needs to be minimal

        # Flux calibration parameters
        par['sensfunc']['UVIS']['extinct_correct'] = False  # This must be False - the extinction correction is performed when making the datacube

        return par




    def get_wcs(self, hdr, slits, platescale, wave0, dwv, spatial_scale=None):
        """
        Construct/Read a World-Coordinate System for a frame.

        Args:
            hdr (`astropy.io.fits.Header`_):
                The header of the raw frame. The information in this
                header will be extracted and returned as a WCS.
            slits (:class:`~pypeit.slittrace.SlitTraceSet`):
                Slit traces.
            platescale (:obj:`float`):
                The platescale of an unbinned pixel in arcsec/pixel (e.g.
                detector.platescale).
            wave0 (:obj:`float`):
                The wavelength zeropoint.
            dwv (:obj:`float`):
                Change in wavelength per spectral pixel.

        Returns:
            `astropy.wcs.WCS`_: The world-coordinate system.
        """
        msgs.info("Calculating the WCS")
        # Get the x and y binning factors, and the typical slit length
        binspec, binspat = parse.parse_binning(self.get_meta_value([hdr], 'binning'))

        # Get the pixel and slice scales
        pxscl = platescale * binspat / 3600.0  # Need to convert arcsec to degrees
        slscl = self.get_meta_value([hdr], 'slitwid')
        if spatial_scale is not None:
            if pxscl > spatial_scale / 3600.0:
                msgs.warn("Spatial scale requested ({0:f}'') is less than the pixel scale ({1:f}'')".format(spatial_scale, pxscl*3600.0))
            # Update the pixel scale
            pxscl = spatial_scale / 3600.0  # 3600 is to convert arcsec to degrees

        # Get the typical slit length (this changes by ~0.3% over all slits, so a constant is fine for now)
        slitlength = int(np.round(np.median(slits.get_slitlengths(median=True))))

        # Get RA/DEC
        raval = self.get_meta_value([hdr], 'ra')
        decval = self.get_meta_value([hdr], 'dec')

        # Create a coordinate
        coord = SkyCoord(raval, decval, unit=(units.deg, units.deg))

        # Get rotator position
        msgs.warn("HACK FOR MAAT SIMS --- NEED TO FIGURE OUT RPOS and RREF FOR MAAT FROM HEADER INFO")
        if 'ROTPOSN' in hdr:
            rpos = hdr['ROTPOSN']
        else:
            rpos = 0.
        if 'ROTREFAN' in hdr:
            rref = hdr['ROTREFAN']
        else:
            rref = 0.
        # Get the offset and PA
        rotoff = 0.0  # IFU-SKYPA offset (degrees)
        skypa = rpos + rref  # IFU position angle (degrees)
        crota = np.radians(-(skypa + rotoff))

        # Calculate the fits coordinates
        cdelt1 = -slscl
        cdelt2 = pxscl
        if coord is None:
            ra = 0.
            dec = 0.
            crota = 1
        else:
            ra = coord.ra.degree
            dec = coord.dec.degree
        # Calculate the CD Matrix
        cd11 = cdelt1 * np.cos(crota)                          # RA degrees per column
        cd12 = abs(cdelt2) * np.sign(cdelt1) * np.sin(crota)   # RA degrees per row
        cd21 = -abs(cdelt1) * np.sign(cdelt2) * np.sin(crota)  # DEC degress per column
        cd22 = cdelt2 * np.cos(crota)                          # DEC degrees per row
        # Get reference pixels (set these to the middle of the FOV)
        crpix1 = 11   # i.e. see get_datacube_bins (11 is used as the reference point - somewhere in the middle of the FOV)
        crpix2 = slitlength / 2.
        crpix3 = 1.
        # Get the offset
        msgs.warn("HACK FOR MAAT SIMS --- Need to obtain offset from header?")
        off1 = 0.
        off2 = 0.
        off1 /= binspec
        off2 /= binspat
        crpix1 += off1
        crpix2 += off2

        # Create a new WCS object.
        msgs.info("Generating MAAT WCS")
        w = wcs.WCS(naxis=3)
        w.wcs.equinox = hdr['EQUINOX']
        w.wcs.name = 'MAAT'
        w.wcs.radesys = 'FK5'
        # Insert the coordinate frame
        w.wcs.cname = ['MAAT RA', 'MAAT DEC', 'MAAT Wavelength']
        w.wcs.cunit = [units.degree, units.degree, units.Angstrom]
        w.wcs.ctype = ["RA---TAN", "DEC--TAN", "WAVE"]
        w.wcs.crval = [ra, dec, wave0]  # RA, DEC, and wavelength zeropoints
        w.wcs.crpix = [crpix1, crpix2, crpix3]  # RA, DEC, and wavelength reference pixels
        w.wcs.cd = np.array([[cd11, cd12, 0.0], [cd21, cd22, 0.0], [0.0, 0.0, dwv]])
        w.wcs.lonpole = 180.0  # Native longitude of the Celestial pole
        w.wcs.latpole = 0.0  # Native latitude of the Celestial pole

        return w




    def get_datacube_bins(self, slitlength, minmax, num_wave):
        r"""
        Calculate the bin edges to be used when making a datacube.

        Args:
            slitlength (:obj:`int`):
                Length of the slit in pixels
            minmax (`numpy.ndarray`_):
                An array with the minimum and maximum pixel locations on each
                slit relative to the reference location (usually the centre
                of the slit). Shape must be :math:`(N_{\rm slits},2)`, and is
                typically the array returned by
                :func:`~pypeit.slittrace.SlitTraceSet.get_radec_image`.
            num_wave (:obj:`int`):
                Number of wavelength steps.  Given by::
                    int(round((wavemax-wavemin)/delta_wave))

        Args:
            :obj:`tuple`: Three 1D `numpy.ndarray`_ providing the bins to use
            when constructing a histogram of the spec2d files. The elements
            are :math:`(x,y,\lambda)`.
        """
        xbins = np.arange(1 + 23) - 11.0 - 0.5  # 23 is for 23 slices, and 11 is the reference slit
        ybins = np.linspace(np.min(minmax[:, 0]), np.max(minmax[:, 1]), 1+slitlength) - 0.5
        spec_bins = np.arange(1+num_wave) - 0.5
        return xbins, ybins, spec_bins






class GTCOSIRISSpectrograph(spectrograph.Spectrograph):
    """
    Child of Spectrograph to handle GTC/OSIRIS specific code (old detector: MAT-44-82)
    """
    ndet = 2
    name = 'gtc_osiris'
    telescope = telescopes.GTCTelescopePar()
    camera = 'OSIRIS'
    url = 'http://www.gtc.iac.es/instruments/osiris/'
    header_name = 'OSIRIS'
    supported = True
    comment = 'See :doc:`gtc_osiris`'



    def get_detector_par(self, det, hdu=None):
        """
        Return metadata for the selected detector.

        Detector data from `here
        <http://www.gtc.iac.es/instruments/osiris/>`__.

        Args:
            det (:obj:`int`):
                1-indexed detector number.
            hdu (`astropy.io.fits.HDUList`_, optional):
                The open fits file with the raw image of interest.  If not
                provided, frame-dependent parameters are set to a default.

        Returns:
            :class:`~pypeit.images.detector_container.DetectorContainer`:
            Object with the detector metadata.
        """
        binning = '1,1' if hdu is None else self.get_meta_value(self.get_headarr(hdu), 'binning')
        #detwin2 = '[1:4102,300:600]' if self.par['rdx']['quicklook'] else '[1:4102,52:1920]'

        # Detector 1
        detector_dict1 = dict(
            binning         = binning,
            det             = 1,
            dataext         = 1,  # Not sure this is used
            specaxis        = 0,
            specflip        = False,
            spatflip        = False,
            platescale      = 0.127,  # arcsec per pixel
            darkcurr        = 0.0,  # e-/pixel/hour
            saturation      = 65535., # ADU
            nonlinear       = 0.95,
            mincounts       = 0,
            numamplifiers   = 1,
            gain            = np.atleast_1d([0.95]),
            ronoise         = np.atleast_1d([4.5]),
            datasec         = np.atleast_1d('[1:4102,280:2048]'),
            oscansec        = np.atleast_1d('[1:4102,6:44]')
            )
        # Detector 2
        detector_dict2 = dict(
            binning         = binning,
            det             = 2,
            dataext         = 2,  # Not sure this is used
            specaxis        = 0,
            specflip        = False,
            spatflip        = False,
            platescale      = 0.127,
            darkcurr        = 0.0,  # e-/pixel/hour
            saturation      = 65535., # ADU
            nonlinear       = 0.95,
            mincounts       = 0,
            numamplifiers   = 1,
            gain            = np.atleast_1d([0.95]),
            ronoise         = np.atleast_1d([4.5]),
            datasec         = np.atleast_1d('[1:4102,52:1920]'),
            oscansec        = np.atleast_1d('[1:4102,6:40]')
            )

        detectors = [detector_dict1, detector_dict2]
        # Return
        return detector_container.DetectorContainer(**detectors[det-1])





    @classmethod
    def default_pypeit_par(cls):
        """
        Return the default parameters to use for this instrument.

        Returns:
            :class:`~pypeit.par.pypeitpar.PypeItPar`: Parameters required by
            all of PypeIt methods.
        """
        par = super().default_pypeit_par()

        # Ignore PCA
        par['calibrations']['slitedges']['sync_predict'] = 'nearest'
        par['calibrations']['slitedges']['bound_detector'] = True

        # Set pixel flat combination method
        par['calibrations']['pixelflatframe']['process']['combine'] = 'median'
        # Wavelength calibration methods
        par['calibrations']['wavelengths']['method'] = 'full_template'
        par['calibrations']['wavelengths']['lamps'] = ['XeI','HgI','NeI','ArI']
        par['calibrations']['wavelengths']['reid_cont_sub'] = False

        # Set the default exposure time ranges for the frame typing
        par['scienceframe']['exprng'] = [90, None]
        par['calibrations']['biasframe']['exprng'] = [None, 0.001]
        par['calibrations']['darkframe']['exprng'] = [999999, None]     # No dark frames
        par['calibrations']['pinholeframe']['exprng'] = [999999, None]  # No pinhole frames
        par['calibrations']['arcframe']['exprng'] = [None, None]  # Long arc exposures
        par['calibrations']['arcframe']['process']['clip'] = False
        par['calibrations']['standardframe']['exprng'] = [None, 180]
        # Multiple arcs with different lamps, so can't median combine nor clip, also need to remove continuum
        par['calibrations']['arcframe']['process']['combine'] = 'mean'
        par['calibrations']['arcframe']['process']['subtract_continuum'] = True
        par['calibrations']['tiltframe']['process']['clip'] = False
        par['calibrations']['tiltframe']['process']['combine'] = 'mean'
        par['calibrations']['tiltframe']['process']['subtract_continuum'] = True

        # Increase the wave tilts order, given the longish slit
        par['calibrations']['tilts']['spat_order'] = 5
        par['calibrations']['tilts']['spec_order'] = 5

        #Only extract one object per standard frame
        par['reduce']['findobj']['maxnumber_std']=1

        # Turn off the 2D fit - this seems to be giving bad results for OSIRIS
        par['reduce']['skysub']['no_poly'] = True

        # Don't extrapolate the sensitivity function for the low resolution gratings
        # Sensitivity function parameters
        par['sensfunc']['extrap_blu'] = 0.0
        par['sensfunc']['extrap_red'] = 0.0
        par['fluxcalib']['extrap_sens'] = True
        par['sensfunc']['algorithm'] = 'IR'
        par['sensfunc']['polyorder'] = 13
        par['sensfunc']['IR']['maxiter'] = 2
        par['sensfunc']['IR']['telgridfile'] = 'TellPCA_3000_26000_R10000.fits'
        return par




    def init_meta(self):
        """
        Define how metadata are derived from the spectrograph files.

        That is, this associates the PypeIt-specific metadata keywords
        with the instrument-specific header cards using :attr:`meta`.
        """
        self.meta = {}
        # Required (core)
        self.meta['ra'] = dict(ext=0, card='RA')
        self.meta['dec'] = dict(ext=0, card='DEC')
        self.meta['target'] = dict(ext=0, card='object')
        self.meta['idname'] = dict(ext=0, card='obsmode')
        self.meta['decker'] = dict(ext=0, card='MASKNAME')
        self.meta['binning'] = dict(card=None, compound=True)  # Uses CCDSUM
        self.meta['detector']=dict(ext=0,card='detector')
        self.meta['mjd'] = dict(ext=0, card='MJD-OBS')
        self.meta['exptime'] = dict(ext=0, card='EXPTIME')
        self.meta['airmass'] = dict(ext=0, card='AIRMASS')
        self.meta['dispname'] = dict(ext=0, card='GRISM')
        self.meta['datasec'] = dict(ext=1, card='DATASEC')
        self.meta['dichroic'] = dict(ext=0, card='FILTER1')
        self.meta['instrument'] = dict(ext=0, card='INSTRUME')




    def compound_meta(self, headarr, meta_key):
        """
        Methods to generate metadata requiring interpretation of the header
        data, instead of simply reading the value of a header card.

        Args:
            headarr (:obj:`list`):
                List of `astropy.io.fits.Header`_ objects.
            meta_key (:obj:`str`):
                Metadata keyword to construct.

        Returns:
            object: Metadata value read from the header(s).
        """
        if meta_key == 'binning':
            binspatial, binspec = parse.parse_binning(headarr[0]['CCDSUM'])
            binning = parse.binning2string(binspec, binspatial)
            return binning
        elif meta_key == 'obstime':
            return Time(headarr[0]['DATE-END'])




    def configuration_keys(self):
        """
        Return the metadata keys that define a unique instrument
        configuration.

        This list is used by :class:`~pypeit.metadata.PypeItMetaData` to
        identify the unique configurations among the list of frames read
        for a given reduction.

        Returns:
            :obj:`list`: List of keywords of data pulled from file headers
            and used to constuct the :class:`~pypeit.metadata.PypeItMetaData`
            object.
        """
        return ['dispname', 'decker', 'binning']




    def raw_header_cards(self):
        """
        Return additional raw header cards to be propagated in
        downstream output files for configuration identification.

        The list of raw data FITS keywords should be those used to populate
        the :meth:`~pypeit.spectrographs.spectrograph.Spectrograph.configuration_keys`
        or are used in :meth:`~pypeit.spectrographs.spectrograph.Spectrograph.config_specific_par`
        for a particular spectrograph, if different from the name of the
        PypeIt metadata keyword.

        This list is used by :meth:`~pypeit.spectrographs.spectrograph.Spectrograph.subheader_for_spec`
        to include additional FITS keywords in downstream output files.

        Returns:
            :obj:`list`: List of keywords from the raw data files that should
            be propagated in output files.
        """
        return ['GRISM', 'MASKNAME', 'CCDSUM', 'obsmode']




    def check_frame_type(self, ftype, fitstbl, exprng=None):
        """
        Check for frames of the provided type.

        Args:
            ftype (:obj:`str`):
                Type of frame to check. Must be a valid frame type; see
                frame-type :ref:`frame_type_defs`.
            fitstbl (`astropy.table.Table`_):
                The table with the metadata for one or more frames to check.
            exprng (:obj:`list`, optional):
                Range in the allowed exposure time for a frame of type
                ``ftype``. See
                :func:`pypeit.core.framematch.check_frame_exptime`.

        Returns:
            `numpy.ndarray`_: Boolean array with the flags selecting the
            exposures in ``fitstbl`` that are ``ftype`` type frames.
        """
        good_exp = framematch.check_frame_exptime(fitstbl['exptime'], exprng)
        if ftype in ['science','standard']:
            return good_exp & (fitstbl['target'] != 'ArcLamp_Xe') \
            & (fitstbl['target'] != 'ArcLamp_HgAr') \
            & (fitstbl['target'] != 'ArcLamp_Ne') \
            & (fitstbl['target'] != 'SpectralFlat') \
                    & (fitstbl['target'] != 'BIAS')
        if ftype in ['arc', 'tilt']:
            return good_exp & ((fitstbl['target'] == 'ArcLamp_Xe') \
            | (fitstbl['target'] == 'ArcLamp_HgAr') \
            | (fitstbl['target'] == 'ArcLamp_Ne'))
        if ftype in ['pixelflat', 'trace', 'illumflat']:
            return good_exp & (fitstbl['target'] == 'SpectralFlat')
        if ftype == 'bias':
            return good_exp & (fitstbl['target'] == 'BIAS')

        msgs.warn('Cannot determine if frames are of type {0}.'.format(ftype))
        return np.zeros(len(fitstbl), dtype=bool)




    def config_independent_frames(self):
        """
        Define frame types that are independent of the fully defined
        instrument configuration.

        Bias and dark frames are considered independent of a configuration.
        Standards are assigned to the correct configuration frame group by
        grism (i.e. ignoring that they are taken with a wider slit).
        See :func:`~pypeit.metadata.PypeItMetaData.set_configurations`.

        Returns:
            :obj:`dict`: Dictionary where the keys are the frame types that
            are configuration independent and the values are the metadata
            keywords that can be used to assign the frames to a configuration
            group.
        """
        return {'standard': 'dispname','bias': 'binning', 'dark': 'binning'}




    def config_specific_par(self, scifile, inp_par=None):
        """
        Modify the PypeIt parameters to hard-wired values used for
        specific instrument configurations.

        Args:
            scifile (:obj:`str`):
                File to use when determining the configuration and how
                to adjust the input parameters.
            inp_par (:class:`~pypeit.par.parset.ParSet`, optional):
                Parameter set used for the full run of PypeIt.  If None,
                use :func:`default_pypeit_par`.

        Returns:
            :class:`~pypeit.par.parset.ParSet`: The PypeIt parameter set
            adjusted for configuration specific parameter values.
        """
        # Start with instrument wide
        par = super().config_specific_par(scifile, inp_par=inp_par)

        if self.get_meta_value(scifile, 'idname') == 'OsirisMOS':
            par['reduce']['findobj']['find_trim_edge'] = [1,1]
            par['calibrations']['slitedges']['sync_predict'] = 'pca'
            par['calibrations']['slitedges']['det_buffer'] = 1
        elif self.get_meta_value(scifile, 'idname') == 'OsirisLongSlitSpectroscopy':
            # Do not tweak the slit edges for longslit
            par['calibrations']['flatfield']['tweak_slits'] = False

        # Wavelength calibrations
        if self.get_meta_value(scifile, 'dispname') == 'R300B':
            par['calibrations']['wavelengths']['lamps'] = ['XeI','HgI','NeI']
            par['calibrations']['wavelengths']['reid_arxiv'] = 'gtc_osiris_R300B.fits'
            par['reduce']['findobj']['find_min_max']=[750,2051]
        elif self.get_meta_value(scifile, 'dispname') == 'R300R':
            par['calibrations']['wavelengths']['lamps'] = ['XeI','HgI','NeI']
            par['calibrations']['wavelengths']['reid_arxiv'] = 'gtc_osiris_R300R.fits'
            par['reduce']['findobj']['find_min_max']=[750,2051]
        elif self.get_meta_value(scifile, 'dispname') == 'R500B':
            par['calibrations']['wavelengths']['lamps'] = ['HgI','NeI']
            par['calibrations']['wavelengths']['reid_arxiv'] = 'gtc_osiris_R500B.fits'
            par['reduce']['findobj']['find_min_max']=[500,2051]
        elif self.get_meta_value(scifile, 'dispname') == 'R500R':
            par['calibrations']['wavelengths']['lamps'] = ['XeI','HgI','NeI']
            par['calibrations']['wavelengths']['reid_arxiv'] = 'gtc_osiris_R500R.fits'
            par['reduce']['findobj']['find_min_max']=[450,2051]
        elif self.get_meta_value(scifile, 'dispname') == 'R1000B':
            par['calibrations']['wavelengths']['lamps'] = ['ArI','HgI','NeI']
            par['calibrations']['wavelengths']['reid_arxiv'] = 'gtc_osiris_R1000B.fits'
        elif self.get_meta_value(scifile, 'dispname') == 'R1000R':
            par['calibrations']['wavelengths']['lamps'] = ['XeI','HgI','NeI']
            par['calibrations']['wavelengths']['reid_arxiv'] = 'gtc_osiris_R1000R.fits'
        elif self.get_meta_value(scifile, 'dispname') == 'R2000B':
            par['calibrations']['wavelengths']['lamps'] = ['XeI','HgI','NeI']
            par['calibrations']['wavelengths']['reid_arxiv'] = 'gtc_osiris_R2000B.fits'
        elif self.get_meta_value(scifile, 'dispname') == 'R2500U':
            par['calibrations']['wavelengths']['lamps'] = ['XeI','HgI']
            par['calibrations']['wavelengths']['reid_arxiv'] = 'gtc_osiris_R2500U.fits'
            par['calibrations']['wavelengths']['nsippet'] = 1
        elif self.get_meta_value(scifile, 'dispname') == 'R2500V':
            par['calibrations']['wavelengths']['lamps'] = ['HgI','NeI','XeI']
            par['calibrations']['wavelengths']['reid_arxiv'] = 'gtc_osiris_R2500V.fits'
        elif self.get_meta_value(scifile, 'dispname') == 'R2500R':
            par['calibrations']['wavelengths']['lamps'] = ['ArI','HgI','NeI']
            par['calibrations']['wavelengths']['reid_arxiv'] = 'gtc_osiris_R2500R.fits'
        elif self.get_meta_value(scifile, 'dispname') == 'R2500I':
            par['calibrations']['wavelengths']['lamps'] = ['ArI','XeI','NeI']
            par['calibrations']['wavelengths']['reid_arxiv'] = 'gtc_osiris_R2500I.fits'
            par['sensfunc']['algorithm'] = 'IR'
            par['sensfunc']['IR']['telgridfile'] = "TellPCA_3000_26000_R10000.fits"
        else:
            msgs.warn('gtc_osiris.py: template arc missing for this grism! Trying holy-grail...')
            par['calibrations']['wavelengths']['method'] = 'holy-grail'

        # Return
        return par




    def bpm(self, filename, det, shape=None, msbias=None):
        """
        Generate a default bad-pixel mask.

        Even though they are both optional, either the precise shape for
        the image (``shape``) or an example file that can be read to get
        the shape (``filename`` using :func:`get_image_shape`) *must* be
        provided.

        Args:
            filename (:obj:`str` or None):
                An example file to use to get the image shape.
            det (:obj:`int`):
                1-indexed detector number to use when getting the image
                shape from the example file.
            shape (tuple, optional):
                Processed image shape
                Required if filename is None
                Ignored if filename is not None
            msbias (`numpy.ndarray`_, optional):
                Processed bias frame used to identify bad pixels. **This is
                ignored for KCWI.**

        Returns:
            `numpy.ndarray`_: An integer array with a masked value set
            to 1 and an unmasked value set to 0.  All values are set to
            0.
        """
        # Call the base-class method to generate the empty bpm; msbias is always set to None.
        bpm_img = super().bpm(filename, det, shape=shape, msbias=None)

        # Extract some header info
        head0 = fits.getheader(filename, ext=0)
        binning = head0['CCDSUM']

        # Construct a list of the bad columns
        bc = []
        if det == 1:
            # No bad pixel columns on detector 1
            pass
        elif det == 2:
            if binning == '1 1':
                # The BPM is based on 2x2 binning data, so the 2x2 numbers are just multiplied by two
                msgs.warn("BPM is likely over-estimated for 1x1 binning")
                bc = [[220, 222, 3892, 4100],
                      [952, 954, 2304, 4100]]
            elif binning == '2 2':
                bc = [[110, 111, 1946, 2050],
                      [476, 477, 1154, 2050]]
        else:
            msgs.warn("Bad pixel mask is not available for det={0:d} binning={1:s}".format(det, binning))
            bc = []

        # Apply these bad columns to the mask
        for bb in range(len(bc)):
            bpm_img[bc[bb][2]:bc[bb][3] + 1, bc[bb][0]:bc[bb][1] + 1] = 1

        return bpm_img




    def tweak_standard(self, wave_in, counts_in, counts_ivar_in, gpm_in, meta_table, trim_std_pixs=None,
                       log10_blaze_function=None, debug=False):
        """
        This routine is for performing instrument/disperser specific tweaks to standard stars so that sensitivity
        function fits will be well behaved. For example, masking second order light. For instruments that don't
        require such tweaks it will just return the inputs, but for instruments that do this function is overloaded
        with a method that performs the tweaks.

        Parameters
        ----------
        wave_in: `numpy.ndarray`_
            Input standard star wavelengths (:obj:`float`, ``shape = (nspec,)``)
        counts_in: `numpy.ndarray`_
            Input standard star counts (:obj:`float`, ``shape = (nspec,)``)
        counts_ivar_in: `numpy.ndarray`_
            Input inverse variance of standard star counts (:obj:`float`, ``shape = (nspec,)``)
        gpm_in: `numpy.ndarray`_
            Input good pixel mask for standard (:obj:`bool`, ``shape = (nspec,)``)
        meta_table: :obj:`dict`
            Table containing meta data that is slupred from the :class:`~pypeit.specobjs.SpecObjs`
            object.  See :meth:`~pypeit.specobjs.SpecObjs.unpack_object` for the
            contents of this table.
        trim_std_pixs: :obj:`list` or :obj:`tuple`, optional
            List or tuple of two integers specifying the number of pixels to
            trim from the start and end of the standard star spectrum. If None,
            no trimming is applied. Default=None.
        log10_blaze_function: `numpy.ndarray`_ or None
            Input blaze function to be tweaked, optional. Default=None.

        Returns
        -------
        wave_out: `numpy.ndarray`_
            Output standard star wavelengths (:obj:`float`, ``shape = (nspec,)``)
        counts_out: `numpy.ndarray`_
            Output standard star counts (:obj:`float`, ``shape = (nspec,)``)
        counts_ivar_out: `numpy.ndarray`_
            Output inverse variance of standard star counts (:obj:`float`, ``shape = (nspec,)``)
        gpm_out: `numpy.ndarray`_
            Output good pixel mask for standard (:obj:`bool`, ``shape = (nspec,)``)
        log10_blaze_function_out: `numpy.ndarray`_ or None
            Output blaze function after being tweaked.
        """
        # If trim_std_pixs is provided, use the base class method to trim the standard star
        if trim_std_pixs is not None:
            return super().tweak_standard(wave_in, counts_in, counts_ivar_in, gpm_in, meta_table,
                                          trim_std_pixs=trim_std_pixs, log10_blaze_function=log10_blaze_function)

        # Could check the wavelenghts here to do something more robust to header/meta data issues
        if 'R300R' in meta_table['DISPNAME']:
            wave_blue = 4800.0  # blue wavelength below which there is contamination
            wave_red = 9300.0  # red wavelength above which the spectrum is contaminated
        elif 'R500R' in meta_table['DISPNAME']:
            wave_blue = 4800.0  # blue wavelength below which there is contamination
            wave_red = 9300.0  # red wavelength above which the spectrum is contaminated
        elif 'R300B' in meta_table['DISPNAME']:
            wave_blue = 3400.0  # blue wavelength below which there is contamination
            wave_red = 6500.0  # red wavelength above which the spectrum is contaminated
        elif 'R500B' in meta_table['DISPNAME']:
            wave_blue = 3400.0  # blue wavelength below which there is contamination
            wave_red = 6500.0  # red wavelength above which the spectrum is contaminated
        elif 'R1000B' in meta_table['DISPNAME']:
            wave_blue = 3500.0  # blue wavelength below which there is contamination
            wave_red = 6300.0  # red wavelength above which the spectrum is contaminated
        else:
            # keep everything the same
            wave_blue = 0.0
            wave_red = np.inf

        second_order_region = (wave_in < wave_blue) | (wave_in > wave_red)
        wave = wave_in.copy()
        counts = counts_in.copy()
        counts_ivar = counts_ivar_in.copy()
        gpm = gpm_in.copy()

        wave[second_order_region] = 0.0
        counts[second_order_region] = 0.0
        counts_ivar[second_order_region] = 0.0
        # By setting the wavelengths to zero, we guarantee that the sensitvity function will only be computed
        # over the valid wavelength region. While we could mask, this would still produce a wave_min and wave_max
        # for the zeropoint that includes the bad regions, and the polynomial fits will extrapolate crazily there
        gpm[second_order_region] = False

        if log10_blaze_function is not None:
            log10_blaze_function_out = log10_blaze_function.copy()
            log10_blaze_function_out[second_order_region] = 0.0
        else:
            log10_blaze_function_out = None

        if debug:
           from matplotlib import pyplot as plt
           from pypeit import utils
           counts_sigma = np.sqrt(utils.inverse(counts_ivar_in))
           plt.plot(wave_in, counts, color='red', alpha=0.7, label='apodized flux')
           plt.plot(wave_in, counts_in, color='black', alpha=0.7, label='flux')
           plt.plot(wave_in, counts_sigma, color='blue', alpha=0.7, label='flux')
           plt.axvline(wave_blue, color='blue')
           plt.axvline(wave_red, color='red')
           plt.legend()
           plt.show()

        return wave, counts, counts_ivar, gpm, log10_blaze_function_out