#!/usr/bin/env python
"""
Create residuals maps.

1. make mapcube with gtbin of the source region
2. run gtsrcmaps,then gtmodel to make a irf-convolved map using fit source model xml file
3. use farith to subtract the maps

@author R. Dubois <richard@slac.stanford.edu>
"""
#

##from setPaths import *
from gt_apps import evtbin
from optparse import OptionParser
from GtApp import GtApp
#import pyfits
import os

srcmaps = GtApp('gtsrcmaps','Likelihood')
model = GtApp('gtmodel','Likelihood')
farith = GtApp('farith')

def run(clean=False):

    exposure = GtApp('gtexposure','pyExposure')

##
## Define command line arguments
##
    parser = OptionParser(version="%prog $Revision 0.1",
                          description="run gtbin")
    parser.add_option("-s","--scfile",action="store",type="string",dest="scfile",
                      default="",help="Filespec of spacecraft file")
    parser.add_option("-e","--evfile",action="store",type="string",dest="evfile",
                      default="",help="Filespec of events file")
    parser.add_option("--alg",action="store",type="string",dest="algorithm",
                      default="CMAP",help="choice of algorithm")
    parser.add_option("--obfile",action="store",type="string",dest="obfile",                  default="outBin.fits",help="gtbin output file name")
    parser.add_option("--ebinalg",action="store",type="string",dest="ebinalg",
                      default="LOG",help="energy binning algorithm")
    parser.add_option("--enumbins",action="store",type="int",dest="enumbins",                         default=10,help="number of energy bins")
    parser.add_option("--emin",action="store",type="float",dest="emin",                      default=100.,help="minimum energy")
    parser.add_option("--emax",action="store",type="float",dest="emax",                      default=300000.,help="minimum energy")
    parser.add_option("--binsz",action="store",type="float",dest="binsz",                      default=0.2,help="angular bin size in image center (deg)")
    parser.add_option("--nxpix",action="store",type="int",dest="nxpix",                         default=50,help="number x pixels")
    parser.add_option("--nypix",action="store",type="int",dest="nypix",                         default=50,help="number y pixels")
    parser.add_option("-c","--coord",action="store",type="string",dest="coord",                      default="CEL",help="coordinate system")
    parser.add_option("--proj",action="store",type="string",dest="proj",                      default="CAR",help="coordinate projection")
    parser.add_option("--ra",action="store",type="float",dest="ra",                      default=0.,help="right ascension")
    parser.add_option("--dec",action="store",type="float",dest="dec",                      default=0.,help="declination")


    parser.add_option("--expCube",action="store",type="string",dest="expCube",                      default="",help="exposure cube")
    parser.add_option("--sourceModel",action="store",type="string",dest="sourceModel",                      default="",help="source model xml file")
    parser.add_option("--bexpmap",action="store",type="string",dest="bexpmap",                      default="bexpmap.fits",help="binned exposure map file")
    parser.add_option("--osfile",action="store",type="string",dest="osfile",                  default="outSrc.fits",help="gtsrcmaps output file name")
    parser.add_option("--irfs",action="store",type="string",dest="irfs",                      default="P6_V1_DIFFUSE",help="irfs")

    parser.add_option("--omfile",action="store",type="string",dest="omfile",                  default="outMod.fits",help="gtmodel output file name")
    parser.add_option("--modelOnly",action="store",type="string",dest="modOnly",                  default="no",help="gtmodel output file name")

    parser.add_option("--orfile",action="store",type="string",dest="orfile",                  default="outRes.fits",help="residuals output file name")

    (options, args) = parser.parse_args()

#    x = pyfits.open(options.evfile);


    evtbin['algorithm'] = 'CCUBE'
    evtbin['evfile'] = options.evfile
    evtbin['outfile'] = options.obfile
    evtbin['scfile'] = options.scfile
    evtbin['ebinalg'] = options.ebinalg
    evtbin['nxpix'] = options.nxpix
    evtbin['nypix'] = options.nypix
    evtbin['binsz'] = options.binsz
    evtbin['coordsys'] = options.coord
    evtbin['xref'] = options.ra
    evtbin['yref'] = options.dec
    evtbin['axisrot'] = 0.0
    evtbin['emin'] = options.emin
    evtbin['emax'] = options.emax
    evtbin['enumbins'] = options.enumbins
    evtbin['denergy'] = 0.0
    evtbin['proj'] =  options.proj
    if (options.modOnly != "yes"):
        evtbin.run()
        evtbin['algorithm'] = 'CMAP'
        evtbin['outfile'] = 'gtbin-map.fits'
        evtbin.run()
    
    srcmaps['scfile'] = options.scfile
    srcmaps['expcube'] = options.expCube
    srcmaps['cmap'] = options.obfile
    srcmaps['srcmdl'] = options.sourceModel
    srcmaps['irfs'] = options.irfs
    srcmaps['bexpmap'] = options.bexpmap
    srcmaps['outfile'] = options.osfile
    if (options.modOnly != "yes"):
        srcmaps.run()

    model['srcmaps'] = srcmaps['outfile']
    model['srcmdl'] = srcmaps['srcmdl']
    model['irfs'] = srcmaps['irfs']
    model['outfile'] = options.omfile
    model['bexpmap'] = srcmaps['bexpmap']
    model['expcube'] = srcmaps['expcube']
    model.run()

    farith['infil1'] = 'gtbin-map.fits'
    farith['infil2'] = model['outfile']
    farith['ops'] = 'sub'
    farith['clobber'] = 'yes'
    farith['datatype'] = '-'
    farith['outfil'] = options.orfile

    cmdFarith = "farith infil1=gtbin-map.fits infil2=" + model['outfile'] + ' outfil=' + options.orfile + ' ops="SUB" datatype="-" overflow=no blank=-999.0 null=no copyprime=yes evenvec=yes clobber=yes'
    print cmdFarith
    os.system(cmdFarith)
#    farith.run()

if __name__ == "__main__":
    run()