ctools

import gammalib as gl

import obsutils as obsu

import os

import ctools as ct

if __name__=="__main__":

    #Create one observation (run)

    #-----------------------------------------------------

    obstime = 360

    ra_degree=83.633212         #Crab position from NED

    dec_degree=22.014460 

    radius=8

    offset_degree=-1            #Observing offset in dec from Crab in deg

    pntdir = gl.GSkyDir()

    pntdir.radec_deg(ra_degree, dec_degree+offset_degree) 

    ob = obsu.set_obs(pntdir,tstart=0.0,duration=obstime,deadc=0.95, \

                  emin=0.1, emax=100.0, rad=radius, irf="irf_file.fits",\

                  caldb="$PWD/caldb/data/cta/aar/bcf/DESY20140105_50h/")

    ob.id("run_1")

    #----------------------------------------------------

    # Create background model

    #----------------------------------------------------

    spec = gl.GModelSpectralPlaw(1,0,gl.GEnergy(1.0,"TeV"))

    filename="$PWD/caldb/data/cta/aar/bcf/DESY20140105_50h/irf_file.fits"

    bck = gl.GCTAModelIrfBackground(spec)

    bck.name("Background_run_1")

    bck.instruments("CTA")

    #----------------------------------------------------

    # Create crab model

    #----------------------------------------------------

    #~ crabdir = gl.GSkyDir()

    #~ crabdir.radec_deg(ra_degree, dec_degree)

    #~ crab_spatial = gl.GModelSpatialPointSource(crabdir)

    #~ crab_spectrum = gl.GModelSpectralPlaw(3.45e-17, -2.63,gl.GEnergy(1,"TeV") ) #from HESS 2006

#~ 

    #~ crab = gl.GModelSky(crab_spatial, crab_spectrum)

    #~ crab.name("CrabNebula")

    #~ crab.ids("")

    #----------------------------------------------------

    # Create GObservations and add models

    #----------------------------------------------------

    obs = gl.GObservations()

    obs.append(ob)

    # obs.models("/home/sanchez/work/CTA/ctools/CrabNebula.xml")

    obs.models("/home/sanchez/work/CTA/ctools/crab.xml")

    #~ obs.models().append(bck)

    #~ obs.models().append(crab)

    #----------------------------------------------------

    # Simulation

    #----------------------------------------------------

    obssim = obsu.sim(obs, log=False, debug=False, edisp=False, seed=0)

    # obssim[0

    import analysisutils as ana

    analyser = ana.Analyser()

    analyser.set_obs(obssim)

    # analyser.ctbin()

    #----------------------------------------------------

    # Fit

    #----------------------------------------------------

    analyser.fit()

    like = analyser.like

    n = "CrabNebula"

    analyser.PrintResults(n)

    # like = obsu.fit(obssim,log=True, debug=False, edisp=False)

    print "TS = ",like.obs().models()[n].ts()

    print 

#    print like.obs().models()

    emin = ob.events().ebounds().emin().TeV()

    emax = ob.events().ebounds().emax().TeV()

    # res = analyser.GetSrcParams(n)

    res_ebin = ana.MakeEbin(analyser,6,n)

    res_ebin.Print()    

    from plotutils import SpectrumPlotter

    spec = SpectrumPlotter(n,like,res_ebin,emin,emax,SED=False)

    spec.draw()