def __init__(self, name, filename):  # name means datid!!
        self.name = name
        sau.load_pha(name, filename)
        self.data = sau.get_data(name)
        self.arf = sau.get_data(name)
        self.rmf = sau.get_data(name)

        try:  # Read keywords from pha header
            self.threshold = self.data.header['ETH']
        except KeyError:
            print(" ! WARNING: no threshold found, using 200 GeV")
            self.threshold = 2e8  # default value 200 GeV
        self.emax = 1e11  # default value 100 TeV

        try:
            self.zenith = self.data.header['ZENITH']
        except KeyError:
            print("WARNING: no mean zenith angle found, using 45 deg")
            self.zenith = 45.0  # default value 200 GeV

        try:
            self.offset = self.data.header['OFFSET']
        except KeyError:
            print("WARNING: no offset angle found, using 1.0 deg")
            self.offset = 1.0  # default value 200 GeV

        try:
            self.telcode = self.data.header['TELCODE']
        except KeyError:
            print("WARNING: no telcode found, using 0")
            self.telcode = 0  # default value 200 GeV

    def __init__(self, name, filename=None):
        self.name = name
        if filename is not None:
            sau.load_pha(name, filename)
            self.data = sau.get_data(name)
            self.arf = self.data.get_arf()
            self.rmf = self.data.get_rmf()

            # Read keywords from pha header
            try:
                self.threshold = self.data.header['ETH']
            except KeyError:
                print("WARNING: no threshold found using 200 GeV")
                self.threshold = 2e8  # default value 200 GeV
                self.emax = 1e11  # default value 100 TeV

            try:
                self.zenith = self.data.header['ZENITH']
            except KeyError:
                print("WARNING: no mean zenith angle found using 45 deg")
                self.zenith = 45.0  # default value 45 deg

            try:
                self.offset = self.data.header['OFFSET']
            except KeyError:
                print("WARNING: no offset angle found using 1.0 deg")
                self.offset = 1.0  # default value 1 deg

            try:
                self.n_tels = self.data.header['N_TELS']
            except KeyError:
                print("WARNING: no number of telescopes found using 0")
                self.n_tels = 0  # default value

            try:
                self.eff = self.data.header['EFFICIEN']
            except KeyError:
                print("WARNING: no efficiency found using 1.0")
                self.eff = 1.00  # default value

            try:
                self.tstart = self.data.header['TSTART']
            except KeyError:
                print("WARNING: no tstart found using 0")
                self.tstart = 0.  # default value

            try:
                self.tstop = self.data.header['TSTOP']
            except KeyError:
                print("WARNING: no tstop found using tsart+1800")
                self.tstop = self.tstart + 1800  # default value

        else:
            self.data = sau.get_data(name)
            self.arf = self.data.get_arf()
            self.rmf = self.data.get_rmf()

 def query(self, func):
     output = []
     for dataset in self.filter_datasets():
         id = dataset['id']
         if func(ui.get_data(id)):
             output.append(id)
     return output

    def setUp(self):

        self._old_logger_level = logger.getEffectiveLevel()
        logger.setLevel(logging.ERROR)

        ui.set_stat('wstat')

        infile = self.make_path('3c273.pi')
        ui.load_pha(1, infile)

        # Change the backscale value slightly so that the
        # results are different to other runs with this file.
        #
        nbins = ui.get_data(1).get_dep(False).size
        bscal = 0.9 * np.ones(nbins) * ui.get_backscal(1)
        ui.set_backscal(1, backscale=bscal)

        ui.set_source(1, ui.powlaw1d.pl)

        # The powerlaw slope and normalization are
        # intended to be "a reasonable approximation"
        # to the data, just to make sure that any statistic
        # calculation doesn't blow-up too much.
        #
        ui.set_par("pl.gamma", 1.7)
        ui.set_par("pl.ampl", 1.7e-4)

    def validate_pha(self, idval):
        """Check that the PHA dataset in id=idval is
        as expected.
        """

        self.assertEqual(ui.list_data_ids(), [idval])

        pha = ui.get_data(idval)
        self.assertIsInstance(pha, DataPHA)

        arf = ui.get_arf(idval)
        self.assertIsInstance(arf, ARF1D)

        rmf = ui.get_rmf(idval)
        self.assertIsInstance(rmf, RMF1D)

        bpha = ui.get_bkg(idval, bkg_id=1)
        self.assertIsInstance(bpha, DataPHA)

        barf = ui.get_arf(idval, bkg_id=1)
        self.assertIsInstance(barf, ARF1D)

        brmf = ui.get_rmf(idval, bkg_id=1)
        self.assertIsInstance(brmf, RMF1D)

        # normally the background data set would have a different name,
        # but this is a  PHA Type 3 file.
        # self.assertEqual(pha.name, bpha.name)
        self.assertEqual(arf.name, barf.name)
        self.assertEqual(rmf.name, brmf.name)

 def test_load_table_fits(self):
     # QUS: why is this not in the sherpa-test-data repository?
     this_dir = os.path.dirname(os.path.abspath(__file__))
     ui.load_table(1, os.path.join(this_dir, 'data', 'two_column_x_y.fits.gz'))
     data = ui.get_data(1)
     self.assertEqualWithinTol(data.x, [1, 2, 3])
     self.assertEqualWithinTol(data.y, [4, 5, 6])

    def setUp(self):
        # hide warning messages from file I/O
        self._old_logger_level = logger.level
        logger.setLevel(logging.ERROR)

        self._id = 1
        fname = self.make_path('3c273.pi')
        ui.load_pha(self._id, fname)
        self._pha = ui.get_data(self._id)

    def _check_stat(self, nbins, expected):

        # check the filter sizes (mainly so that these tests
        # get flagged as in need of a look if anything changes
        # in other parts of the code, such as filtering and binning
        #
        self.assertEqual(nbins, ui.get_data().get_dep(True).size)

        stat = ui.calc_stat()
        self.assertAlmostEqual(expected, stat, places=7)

    def testReadImplicit(self):
        """Exclude .gz from the file name"""

        idval = "13"
        fname = self.head + '_pha3.fits'
        ui.load_pha(idval, fname)

        self.validate_pha(idval)

        pha = ui.get_data(idval)
        bpha = ui.get_bkg(idval, bkg_id=1)
        self.assertEqual(pha.name, bpha.name)

 def center_psf(self):
     """Set xpos and ypos of the PSF to the dataspace center"""
     import sherpa.astro.ui as sau
     try:
         ny, nx = sau.get_data().shape
         for par in sau.get_psf().kernel.pars:
             if par.name is 'xpos':
                 par.val = (nx + 1) / 2.
             elif par.name is 'ypos':
                 par.val = (ny + 1) / 2.
     except:
         logging.warning('PSF is not centered.')

 def center_psf(self):
     """Set ``xpos`` and ``ypos`` of the PSF to the dataspace center."""
     import sherpa.astro.ui as sau
     try:
         ny, nx = sau.get_data().shape
         for par in sau.get_psf().kernel.pars:
             if par.name is 'xpos':
                 par.val = (nx + 1) / 2.
             elif par.name is 'ypos':
                 par.val = (ny + 1) / 2.
     except:
         raise Exception('PSF is not centered.')

    def center_psf(self):
        """Set ``xpos`` and ``ypos`` of the PSF to the dataspace center."""
        import sherpa.astro.ui as sau
        try:
            ny, nx = sau.get_data().shape
            for _ in ['psf1', 'psf2', 'psf3']:
                par = sau.get_par(_ + '.xpos')
                par.val = nx / 2.

                par = sau.get_par(_ + '.ypos')
                par.val = ny / 2.
        except:
            raise Exception('PSF is not centered.')

 def containment_fraction(self, theta, npix=1000):
     """Compute fraction of PSF contained inside theta."""
     import sherpa.astro.ui as sau
     sau.dataspace2d((npix, npix))
     self.set()
     # x_center = get_psf().kernel.pars.xpos
     # y_center = get_psf().kernel.pars.ypos
     x_center, y_center = sau.get_psf().model.center
     x_center, y_center = x_center + 0.5, y_center + 0.5  # shift seen on image.
     x, y = sau.get_data().x0, sau.get_data().x1
     # @note Here we have to use the source image, before I used
     # get_model_image(), which returns the PSF-convolved PSF image,
     # which is a factor of sqrt(2) ~ 1.4 too wide!!!
     p = sau.get_source_image().y.flatten()
     p /= np.nansum(p)
     mask = (x - x_center) ** 2 + (y - y_center) ** 2 < theta ** 2
     fraction = np.nansum(p[mask])
     if 0:  # debug
         sau.get_data().y = p
         sau.save_data('psf_sherpa.fits', clobber=True)
         sau.get_data().y = mask.astype('int')
         sau.save_data('mask_sherpa.fits', clobber=True)
     return fraction

    def testReadExplicit(self):
        """Include .gz in the file name"""

        idval = 12
        fname = self.head + '_pha3.fits.gz'
        ui.load_pha(idval, fname)

        self.validate_pha(idval)

        # TODO: does this indicate that the file name, as read in,
        #       should have the .gz added to it to match the data
        #       read in, or left as is?
        pha = ui.get_data(idval)
        bpha = ui.get_bkg(idval, bkg_id=1)
        self.assertEqual(pha.name, bpha.name + '.gz')

 def test_xmm2(self):
     self.run_thread('xmm2')
     self.assertEqualWithinTol(ui.get_data().channel[0], 1.0, 1e-4)
     self.assertEqual(ui.get_rmf().detchans, 800)
     self.assertEqual(len(ui.get_rmf().energ_lo), 2400)
     self.assertEqual(len(ui.get_rmf().energ_hi), 2400)
     self.assertEqual(len(ui.get_rmf().n_grp), 2400)
     self.assertEqual(len(ui.get_rmf().f_chan), 2394)
     self.assertEqual(len(ui.get_rmf().n_chan), 2394)
     self.assertEqual(len(ui.get_rmf().matrix), 1281216)
     self.assertEqual(ui.get_rmf().offset, 0)
     self.assertEqual(len(ui.get_rmf().e_min), 800)
     self.assertEqual(len(ui.get_rmf().e_max), 800)
     self.assertEqual(len(ui.get_arf().energ_lo), 2400)
     self.assertEqual(len(ui.get_arf().energ_hi), 2400)
     self.assertEqual(len(ui.get_arf().specresp), 2400)

    def __init__(self, dataids):
        self.datasets = []
        self.tot_excess = None
        self.tot_expo = None

        for dataid in dataids:
            spec = spectral_data(sau.get_data(dataid), sau.get_bkg(dataid))
            self.datasets.append(spec)
        self.bkg = np.concatenate([a.fit_bkg for a in self.datasets])
        self.alpha = np.concatenate([a.fit_alpha for a in self.datasets])
        self.Ttot = np.concatenate([a.fit_Ttot for a in self.datasets])
        self.ONexpo = np.concatenate([a.fit_ONexpo for a in self.datasets])

        for a in self.datasets:
            # Carefull we are assuming that all the spectra have the same binning
            if self.tot_excess is None:
                self.tot_excess = np.zeros_like(a.excess)
            if self.tot_expo is None:
                self.tot_expo = np.zeros_like(a.excess)

            self.tot_excess += a.excess
            self.tot_expo += a.full_expo

    def query(self, func):
        """Return the data sets identified by a function.

        Parameters
        ----------
        func
           A function which accepts a Sherpa Data object and
           returns ``True`` if the data set matches.

        Returns
        -------
        matches : list of int or str
           A list of the data set identifiers that match. This
           list may be empty.

        See Also
        --------
        query_by_header_keyword, query_by_obsid

        Examples
        --------

        This query function selects those data sets which have
        been grouped, or had their background subtracted, or both.

        >>> myquery = lambda d: d.subtracted or d.grouped
        >>> query(myquery)
        []

        """

        output = []
        for dataset in self.filter_datasets():
            id = dataset['id']
            if func(ui.get_data(id)):
                output.append(id)
        return output

        def _load(self, *args, **kwargs):
            """Load a dataset and add to the datasets for stacked analysis.
            """
            if self.getitem_ids:
                dataid = self.getitem_ids[0]
                self.getitem_ids = None
            else:
                dataid = 1
                while dataid in _all_dataset_ids:
                    dataid += 1
                
            if dataid in self.dataset_ids:
                raise ValueError('Data ID = {0} is already in the DataStack'.format(dataset['id']))

            logger.info('Loading dataset id %s' % dataid)
            out = load_func(dataid, *args, **kwargs)

            dataset = dict(id=dataid, args=args, model_comps={}, data=ui.get_data(dataid))
            dataset.update(kwargs)  # no sherpa load func 'args' keyword so no conflict
            _all_dataset_ids[dataid] = dataset
            self.dataset_ids[dataid] = dataset
            self.datasets.append(dataset)

            return out

    def load_pha(self, specfile, annulus):
        """Load a pha file and add to the datasets for stacked analysis.

        It is required that datasets for all annuli are loaded before
        the source model is created (to ensure that components are
        created for each annulus).

        Parameters
        ----------
        specfile : str or sherpa.astro.data.DataPHA object
            If a string, the name of the file containing the source spectrum,
            which must be in PHA format (the data is expected to be extracted
            on the PI column). If a DataPHA object, then this is used (and
            is assumed to contain any needed background data).
        annulus : int
            The annulus number for the data.

        Returns
        -------
        dataid : int
            The Sherpa dataset identifier used for this spectrum.

        Examples
        --------

        Load the data for four annuli from the files 'ann1.pi' to 'ann4.pi'.

        >>> dep.load_pha('ann1.pi', 0)
        >>> dep.load_pha('ann2.pi', 1)
        >>> dep.load_pha('ann3.pi', 2)
        >>> dep.load_pha('ann4.pi', 3)

        Load in the PHA files into Sherpa DataPHA objects, and then use
        these objects:

        >>> s1 = ui.unpack_pha('src1.pi')
        >>> s2 = ui.unpack_pha('src2.pi')
        >>> s3 = ui.unpack_pha('src3.pi')
        >>> dep.load_pha(s1, 0)
        >>> dep.load_pha(s2, 1)
        >>> dep.load_pha(s3, 2)

        """

        dataid = len(self.datasets)

        # If the input has a counts attribute then assume a DataPHA
        # style object.
        #
        if hasattr(specfile, 'counts'):
            print('Using spectrum {} '.format(specfile.name) +
                  ' as dataset id {}'.format(dataid))
            ui.set_data(dataid, specfile)

        else:
            print('Loading spectrum file {} '.format(specfile) +
                  ' as dataset id {}'.format(dataid))
            ui.load_pha(dataid, specfile)

        data = ui.get_data(dataid)
        try:
            obsid = int(data.header['OBS_ID'])
        except (KeyError, TypeError, ValueError):
            obsid = 0

        dataset = dict(file=specfile,
                       obsid=obsid,
                       id=dataid,
                       annulus=annulus
                       )
        self.datasets.append(dataset)
        self.obsids.add(obsid)
        return dataid

def mainloop(mymodel, fwhm, id = None, maxiter = 5, mindist = 0., do_plots = 0):
    
    if id is None:
        id = ui.get_default_id()
    data = ui.get_data(id)
    wave = data.get_indep()[0]
    error = data.get_error()[0]
    
    # model could habe been initalized with arbitrary values
    ui.fit(id) 

    for i in range(maxiter):
        oldmodel = smh.get_model_parts(id)
        res_flux = ui.get_resid_plot(id).y
        if smoothwindow is not None:
            fwhminpix = int(fwhm / np.diff(wave).mean())
            y = smooth(res_flux/error, window_len = 3*fwhminpix, window = smoothwindow)
        else:
            y = res_flux/error
        peaks = findlines(wave, y, fwhm, smoothwindow = None, sigma_threshold = sigma_threshold)
        if has_mpl and (do_plots > 2):
            plt.figure()
            plt.plot(wave, res_flux/error, 's')
            for pos in mymodel.line_value_list('pos'):
                plt.plot([pos, pos], plt.ylim(),'k:')
            for peak in peaks:
                plt.plot([wave[peak], wave[peak]], plt.ylim())
            plt.plot(wave, y)
            plt.draw()
            
        for peak in peaks:
            if (len(mymodel.line_value_list('pos')) == 0) or (min(np.abs(mymodel.line_value_list('pos') - wave[peak])) >= mindist):
                mymodel.add_line(**mymodel.guess(wave, smooth(res_flux, window_len = 3*fwhminpix, window = smoothwindow), peak, fwhm = fwhm))
        newmodel = smh.get_model_parts(id)
        print 'Iteration {0:3n}: {1:3n} lines added'.format(i, len(newmodel) - len(oldmodel))
        
        if set(newmodel) == set(oldmodel):
            print 'No new lines added this step - fitting finished'
            break
        # Now do the fitting in Sherpa
        #ui.set_method('simplex')
        ui.fit(id)
        #ui.set_method('moncar')
        #ui.fit(id)
        
        if has_mpl and (do_plots > 0):
            if do_plots > 1:
                plt.figure()
            else:
                plt.clf()
            ui.plot_fit(id)
            for pos in mymodel.line_value_list('pos'):
                plt.plot([pos, pos], plt.ylim(),'k:')
            for peak in peaks:
                plt.plot([wave[peak], wave[peak]], plt.ylim())
            plt.plot(wave, res_flux)
            plt.draw()
        

    else:
        print 'Max number of iterations reached'
    #model.cleanup() #remove lines running to 0 etc.
    return mymodel