tev-fit-source.py

#!/usr/bin/env python
"""Fit the morphology of a number of sources using
initial parameters from a JSON file (for now only Gaussians)"""
#@PydevCodeAnalysisIgnore
import logging
from optparse import OptionParser
from sherpa.astro.ui import *
import morphology.utils
import morphology.psf

# ---------------------------------------------------------
# Parse command line arguments
# ---------------------------------------------------------
parser = OptionParser()
parser.add_option("--counts", type=str, default='counts.fits',
                  help="Counts FITS file name "
                  "[default=%default]")
parser.add_option("--exposure", type=str, default='exposure.fits',
                  help="Exposure FITS file name "
                  "[default=%default]")
parser.add_option("--background", type=str, default='background.fits',
                  help="Background FITS file name "
                  "[default=%default]")
parser.add_option("--psf", type=str, default='psf.json',
                  help="PSF JSON file name "
                  "[default=%default]")
parser.add_option("--sources", type=str, default='sources.json',
                  help="Sources JSON file name (contains start "
                  "values for fit of Gaussians)"
                  "[default=%default]")
parser.add_option("--roi", type=str, default='roi.reg',
                  help="Region of interest (ROI) file name (ds9 reg format)"
                  "[default=%default]")
parser.add_option("--fit_result", type=str, default=None,
                  help="Output JSON file with fit results"
                  "[default=%default]")
(options, args) = parser.parse_args()

if options.fit_result is None:
    outfile = options.sources[:-5] + '_fit_results.json'
else:
    outfile = options.fit_result

# ---------------------------------------------------------
# Load images, PSF and sources
# ---------------------------------------------------------
logging.info('Clearing the sherpa session')
clean()

logging.info('Reading counts: {0}'.format(options.counts))
load_data(options.counts)

logging.info('Reading exposure: {0}'.format(options.exposure))
load_table_model("exposure", options.exposure)

logging.info('Reading background: {0}'.format(options.background))
load_table_model("background", options.background)

logging.info('Reading PSF: {0}'.format(options.psf))
morphology.psf.Sherpa(options.psf).set()

if options.roi:
    logging.info('Reading ROI: {0}'.format(options.roi))
    notice2d(options.roi)
else:
    logging.info('No ROI selected.')

logging.info('Reading sources: {0}'.format(options.sources))
morphology.utils.read_json(options.sources, set_source)

# ---------------------------------------------------------
# Set up the full model and freeze PSF, exposure, background
# ---------------------------------------------------------
# Scale exposure by 1e-10 to get ampl or order unity and avoid some fitting problems
set_full_model('background + 1e-10 * exposure * psf (' + get_source().name + ')')
freeze(background, exposure, psf)

# ---------------------------------------------------------
# Set up the fit
# ---------------------------------------------------------
set_coord('physical')
set_stat('cash')
set_method('levmar') # levmar, neldermead, moncar
set_method_opt('maxfev', int(1e3))
set_method_opt("verbose", 10)

# ---------------------------------------------------------
# Fit and save information we care about
# ---------------------------------------------------------

#show_all() # Prints info about data and model
fit() # Does the fit
covar() # Computes symmetric errors (fast)
#conf() # Computes asymmetric errors (slow)
#image_fit() # Shows data, model, residuals in ds9
logging.info('Writing %s' % outfile)
morphology.utils.write_all(outfile)